地殻構造探査学

データ解析演習（課題その3）

バイブロサイスと自己相関、相互相関

この課題ではサンプルデータ demo_imp.su を使用する。データがなければここを参考にしてコピーする。

反射法では、

地震波形＝震源波形＊反射係数列＋ノイズ

という線形モデルを採用して、時系列解析を適用する。ここで＊はコンボリューションを表す。

まず、反射係数列 demo_imp.su を表示する。

$ suxwigb < demo_imp.su &

画面の解像度や描画の具合で反射係数のインパルスが見えなかったり、太さが異なって見えたりすることがあるかもしれない。ウインドウの大きさを変えると見えるようになる。全部で5つのインパルスがあることを確認しよう。

図：反射係数列。縦軸は時間 (s)。

バイブロサイス発震で用いられるスィープ波形（時間とともに周波数が変化する波形）を作成し表示する。

$ suvibro f1=10 f2=100 tv=1 t1=0.1 t2=0.1 dt=0.001 > tmp_vib.su
$ suxwigb < tmp_vib.su &

suvibro は tv 秒間に周波数が f1 から f2 まで変化（スィープ）する波形を作成する。
dt はサンプル秒。t1, t2 は先端と末端の振幅を小さくするウィンドウ（テーパという）の時間幅。

スィープ波の自己相関関数を表示する。

$ suacor < tmp_vib.su ntout=201 | sushw key=delrt a=-100 | suxwigb &
$ suacor < tmp_vib.su ntout=201 sym=0 | suxwigb &

suacor は自己相関関数 (auto correlation) を計算するプログラム。
自己相関は時間軸で対称であることに注目。時間軸の中心が遅れ時間（ラグ）=0 である。
sushw はヘッダに値をセットするプログラム。ここでは時間軸の中心が 0 となるように時間軸のスタート（の遅れ）位置 delrt を -0.1 (s) = -100 (ms) ずらしている。
2番目の例のように 'sym=0' とすると、遅れ時間（ラグ）が正の部分のみを出力する。
この波形は Clauder wavelet と呼ばれている。

図：（左）10〜100 Hz のスィープ波形。縦軸は時間 (s)。（右）その自己相関。縦軸はラグ時間 (s)。

反射係数列に震源波形（スィープ波）をコンボリューションする。さらにノイズを加える。これが観測される地震波形であると考えることができる。

$ suconv sufile=tmp_vib.su < demo_imp.su > tmp_cnv.su
$ suxwigb < tmp_cnv.su &
$ suaddnoise sn=10 < tmp_cnv.su > tmp_cnv_noise.su
$ suxwigb < tmp_cnv_noise.su &

中間結果を出力する必要がなければ、コマンドをパイプで接続してもよい。

$ suconv sufile=tmp_vib.su < demo_imp.su | suaddnoise sn=10 > tmp_cnv_noise.su
$ suxwigb < tmp_cnv_noise.su &

suconv はコンボリューションを行うプログラム。
suaddnoise は（ランダムな）ノイズを加えるプログラム。
sn=10 は S/N 比、大きいほどノイズが小さくなる。sn=の値を変えてみよう。

図：（左）コンボリューションした結果、ノイズなし。縦軸は時間 (s)。（右）コンボリューションした結果＋ノイズ。

バイブロサイス発震の時にはこのような波形記録が観測されることになる。すなわち、

地震波形＝震源波形＊反射係数列＋ノイズ

この観測波形記録から反射係数列を推定することがここでの目的である。

波形記録と元のスィープ波（震源波形）の相互相関を計算する。

$ suxcor sufile=tmp_vib.su < tmp_cnv.su | suxwigb &

suxcor は相互相関関数 (cross correlation) を計算するプログラム。
'sufile=' で与えた時系列と、入力の時系列との相互相関を計算する。

図：（左）相互相関の結果、ノイズなし。縦軸は時間 (s)。（右）元の反射係数列。

線形の応答モデルでは、

地震波形 ＝ 震源波形 ＊ 反射係数列 ＋ ノイズ

と仮定した。相互相関の演算を cc で表すと、

地震波形 cc 震源波形
  ＝（震源波形 ＊ 反射係数列 ＋ ノイズ） cc 震源波形
  ＝（震源波形 ＊ 反射係数列） cc 震源波形 ＋ ノイズ cc 震源波形
  ≒（震源波形 cc 震源波形） ＊ 反射係数列
  ＝ 震源波形の自己相関 (Clauder Wavelet) ＊ 反射係数列

観測波形のスィープ波が圧縮されて（Clauder waveletに変換されて）反射係数列に近い結果が得られていることがわかる。また、ノイズの影響が弱められていることがわかる。

デコンボリューション（予測デコンボリューション）

この課題ではサンプルデータ demo_dcn.su を使用する。データがなければここを参考にしてコピーする。

まず、サンプルデータ demo_dcn.su の波形と自己相関を表示する。

$ suxwigb < demo_dcn.su &
$ suacor < demo_dcn.su | sushw key=delrt a=-200 | suxwigb &

図：（左）波形。縦軸は時間 (s)。（右）自己相関。縦軸はラグ時間 (s)。

波形と自己相関を見ると、振幅の大きい繰り返し（よく見ると正負反転している）が卓越していることがわかる。これは下図に示すような、浅海の記録でよく見られる水中を多重反射した波（水中残響という）の例である（多重反射をするたびに走時が 2t だけ遅れ、海面で反射するたびに正負が反転する）。

微弱な反射波はこれらに隠されてしまうので、水中残響を取り除くことが必要である。

図：水中残響の模式図。

まず、波形をパルス（スパイク）に変換するデコンボリューションフィルタを適用する。結果と自己相関を表示する。

$ supef < demo_dcn.su maxlag=???? > tmp_decon.su
$ suxwigb < tmp_decon.su &
$ suacor < tmp_decon.su | sushw key=delrt a=-200 | suxwigb &

supef は予測誤差フィルタ (predictive error filter) を計算するプログラム。デコンボリューション (deconvolution) は予測誤差フィルタの特殊な条件の場合に相当する。
maxlag はフィルタのオペレータ長(s)。適切な設定は難しいので試行錯誤して決める。自己相関が参考になる。

maxlag=???? には0.01〜0.10 (s)程度の値を入れて試せ。適切な値だと、下記のように左の波形が右のようにパルス（スパイク）状に変換される。

図：（左）元の波形（拡大）（右）デコンボリューションの結果（拡大）。

図：（左）デコンボリューションの結果の波形。縦軸は時間 (s)。（右）その自己相関。縦軸はラグ時間 (s)。波形がパルス（スパイク）状に変換されたことは自己相関を見てもわかる。

次に、繰り返しを除去する予測誤差フィルタを適用する。

$ supef minlag=???? maxlag=???? < tmp_decon.su | suxwigb &
$ supef minlag=???? maxlag=???? < tmp_decon.su | suacor | sushw key=delrt a=-200 | suxwigb &

フィルタのパラメータ minlag, maxlag はフィールドデータでは経験的に下図のように自己相関の形状から推定することが一般的だが、必ずうまくいくとは限らないので、ある程度の試行錯誤が必要である。