地殻構造探査学

データ解析演習(課題その6)実データの解析

課題その5へ < 課題その6 > 課題その7へ

処理の流れ

  1. ヘッダのチェック
  2. 周波数フィルタ
  3. 振幅回復
  4. 速度フィルタ(今回は必要ない)
  5. デコンボリューション(今回は必要ない)
  6. 静補正(今回は必要ない)
  7. 速度解析(次で)
  8. NMO補正とCMP重合(次で)
  9. マイグレーション(次で)

データの確認と表示

この演習で処理するデータは Nshots.su である。あるフィールド(海域)で取得された実データのごく一部で、それでも 423 MBある大きなデータである。 実際の調査データはTBオーダーになる。

まず、どのようなデータか、データのヘッダを見る。実際のデータ解析はヘッダの編集から処理が始まるのだが、このデータはすでにヘッダの編集が済んでおり、基本的な情報が記載されている。

$ surange < Nshots.su
出力は(データが大きいのでちょっと時間がかかる)
-----                    内容
19057 traces:                              : トレースの数
tracl    39069 58125 (39069 - 58125)       : トレース番号
tracr    1 19057 (1 - 19057)               : トレース番号
fldr     1687 2012 (1687 - 2012)           : 発震点番号
tracf    28 96 (96 - 28)                   : 観測点番号
cdp      900 1300 (900 - 1300)             : CDP (CMP) 番号
cdpt     1 69 (1 - 1)
trid     1 2 (1 - 1)                       : データ種別
offset   -2435 -170 (-170 - -2435)         : オフセット距離 (m)
scalel   -10000
scalco   -10000
counit   1
muts     0 11000 (0 - 0)
ns       5500                              : サンプル数
dt       2000                              : サンプル間隔 (micro second)
day      206
hour     21 22 (21 - 22)
minute   0 59 (3 - 20)
sec      1 59 (45 - 3)
-----
ヘッダの値が取る範囲が表示される。

次に、データを表示させてみる。 データが大きいので、全部表示せずに一部を切り出す。

例えば、発震点番号(ここではヘッダ fldr に設定されている)が 1800〜1805 までのトレースを取り出して、観測点番号(ここでは tracf )の順に並び替え、tmp_1800.suというデータを作る (取り出す発震点番号の値や範囲、ファイル名は任意でよい)。

$ suwind key=fldr min=1800 max=1805 < Nshots.su > tmp.su
$ susort tracf < tmp.su > tmp_1800.su
$ rm tmp.su
いくつかのヘッダ (fldr, tracf, cdp, offset) の内容を表示する。
$ surange < tmp_1800.su
$ sugethw key=fldr,tracf,cdp,offset < tmp_1800.su
データを表示する。
$ suxwigb perc=98 < tmp_1800.su &
$ suximage perc=98 < tmp_1800.su &
これは海上探査の記録である。記録から、 といったことがわかる。

CDP (CMP) の範囲は 900 - 1300である。 事前の準備課題その1 で扱った処理後のデータ Nstack.su の surange の出力では

cdp      1 2869 (1 - 2869)
である。実は、Nshot.su のデータの観測範囲は Nstack.su のデータの観測範囲の一部である。


テストデータによる前処理のパラメータ決定

前処理のパラメータを決めるために、練習用のデータを切り出す。例えば CDP (CMP) が 1100 のトレースを取り出す。取り出す CDP (CMP) の値は例であり 900〜1300 の間で任意でよい。
$ suwind key=cdp min=1100 max=1100 < Nshots.su > test_cdp1100.su
表示してみる。
$ suximage < test_cdp1100.su perc=99 &
$ suxwigb < test_cdp1100.su perc=99 &

図:波形表示(拡大)、(左)濃淡(イメージ)表示、(右)波形表示


スペクトルと周波数フィルタ

スペクトルを計算して表示する。
$ suspecfx < test_cdp1100.su | suximage perc=98 &
$ suspecfx < test_cdp1100.su | suxgraph &
$ suspecfx < test_cdp1100.su | suop op=db | suxgraph &

図:スペクトル、(左)カラー表示、(右)グラフ(dB)表示。縦軸は周波数(Hz)。

各トレースにどのような周波数成分が含まれているかわかる。

各周波数帯域にどのような成分が含まれているかを見てみる。 簡単なシェルスクリプト create_BPF_panel.sh を使用する。 自分のファイル名に合わせてスクリプトを編集すること。 10 Hz ごとに幅 10 Hz でバンドパスフィルタをかけて並べて表示する。

$ sh create_BPF_panel.sh

図:スペクトル、(左)オリジナル、(右)それぞれ 0-10 Hz, 10-20 Hz, 20-30 Hz, 30-40 Hz,.....のバンドパスフィルタ適用後。

バンドパスフィルタをかけて表示してみる。下の例では f=10,20,70,80 としている。

$ sufilter f=10,20,70,80 amps=0,1,1,0 < test_cdp1100.su | suspecfx | suximage perc=98 &
$ sufilter f=10,20,70,80 amps=0,1,1,0 < test_cdp1100.su | suspecfx | suxgraph &
$ sufilter f=10,20,70,80 amps=0,1,1,0 < test_cdp1100.su | suspecfx | suop op=db | suxgraph &
$ sufilter f=10,20,70,80 amps=0,1,1,0 < test_cdp1100.su | suximage perc=99&
まず、海水中(0~6秒くらい)に反射波はないはずなので、この時間帯に見られる低周波数および高周波数の波群はノイズであると判断し、これらを除去することにする。 フィルタのパラメータを変えて波形表示をさせていろいろ試し、適切な値が見つかったらフィルタをかけて保存する。
$ sufilter f=?,?,?,? amps=0,1,1,0 < test_cdp1100.su > test_cdp1100_flt.su

図:スペクトルのカラー表示、(左)オリジナル、(右)バンドパスフィルタ適用後。

図:スペクトルのグラフ表示(dB表示)、(左)オリジナル、(右)バンドパスフィルタ適用後。

図:波形のイメージ表示、(左)オリジナル、(右)バンドパスフィルタ適用後。主に低周波数のノイズが低減されていることがわかる。


振幅回復

データを表示しておく。
$ suximage perc=98 < test_cdp1100_flt.su &
$ suop op=nop  < test_cdp1100_flt.su | suattributes mode=amp | suop op=db | suxgraph &
海底からの反射波が含まれている走時 6 秒以降に着目する。走時が大きくなるにつれて振幅が減少しているので、深部からの反射波を強調するためにゲイン(増幅)をかける。ここでは 6 秒以降の振幅値がほぼフラットになるようにするとよい。
$ sugain tpow=2.0 < test_cdp1100_flt.su | suattributes mode=amp | suop op=db | suxgraph &
tpow の値をいろいろ変えて試し、6秒以降の振幅がおおむね等しくなるような適切な値が見つかったら振幅を補正してデータを保存する。
$ sugain tpow=? < test_cdp1100_flt.su > test_cdp1100_flt_gain.su

図:波形の振幅表示、(左)適用前(フィルタ後)、(右)振幅回復(ゲイン)適用後。6 秒以降の振幅値がほぼフラットになっていることがわかる。

図:波形のイメージ表示、(左)適用前(フィルタ後)、(右)振幅回復(ゲイン)適用後。深部からの反射波の振幅が増加している。


速度フィルタ

f-kスペクトルを表示すると以下のようになる。
$ suspecfk < test_cdp1100_flt_gain.su | suximage perc=98 &

図:f-kスペクトル。

このデータは海上探査の記録であり、大きな表面波ノイズは記録されていない。 また、低周波数・低速度のノイズはバンドパスフィルタで除去されているので、 特に速度フィルタを適用する必要はないと思われる。省略してよい

速度フィルタをショットギャザー(発震点記録)に対して適用するには、 ショット(発震点)ごとにループを回す処理が必要になるので、 スクリプトを書いて行うことになるだろう。


デコンボリューション

このデータでは特にデコンボリューションを適用する必要がない(あまり効果がないのに対して副作用が目立つ)ので省略する。

興味のある人はここを参考にして適用し、結果の違いを比べてみるとよいだろう。


前処理のための下調べはここまで。作成したファイル、例えば test_cdp1100_flt_gain.su(デコンボリューションをした人は test_cdp1100_flt_gain_decon.su に読み替える) は以後も使用するので消さないで置いておくこと。


適用

前処理のパラメータを決めたので、実際に全体のデータセットに適用する。
  1. バンドパスフィルタを適用する。
  2. 振幅回復を行う。
  3. デコンボリューションを適用する(省略可)。
  4. 再度バンドパスフィルタを適用する(省略可)。
  5. データが大きすぎるので、リサンプルしてデータ数を少なくする。 サンプリング間隔を 2 ms から 4 msに変更する。 ダウンサンプリング時にはアンチエリアスフィルタを適用(ナイキスト周波数よりも高い周波数成分を除去する) しないといけないが、 先にバンドパスフィルタを適用したのでここでは必要ない。
処理を一つ一つ行ってもよいが、 この一連の処理はパイプを使って以下のように一行で書ける。
$ sufilter f=?,?,?,? amps=0,1,1,0 < Nshots.su | sugain tpow=? | suresamp nt=2750 dt=0.004 > Ntmp.su
? の部分に各自で決めた数値を入れる。

ちなみに、デコンボリューションを入れる人は以下のようになる。

$ sufilter f=?,?,?,? amps=0,1,1,0 < Nshots.su | sugain tpow=? | supef minlag=???? | sufilter f=?,?,?,? amps=0,1,1,0 | suresamp nt=2750 dt=0.004 > Ntmp.su

最後に、cdp 順にソートする。各CDP (CMP)の中は offset 順に並べておく。

$ susort cdp offset < Ntmp.su > Ncdps.su

これで前処理ができた。

-rw-r--r-- 1 watanabe staffs 214200680 Nov  7 12:17 Ncdps.su
-rw-r--r-- 1 watanabe staffs 423827680 Nov  7 12:17 Nshots.su
リサンプリングをしたので新しいデータは前のデータのほぼ半分のサイズになっている。 Ntmp.su は不要なので消去してよい。


課題その5へ < 課題その6 > 課題その7へ
戻る
Last modified: Fri Sep 27 14:13:55 JST 2024