微動アレイ解析ツールBIDO

６．デモデータによるプログラムの実行例

ここではデモデータを用いて以下の３つの場合についてＢＩＤＯの使い方を解説します．

・前処理のデモ（こちら）

・ハドルテストのデモ（こちら）

・アレイデータ解析のデモ（すぐ下に続く）

アレイデータ解析のデモ

アレイ観測を想定してデモ用に人工的に合成したデータを用いてプログラムの実行方法を示します．Windowsでの実行例を示しますが、Linuxでもやり方は同じです．

まず解凍ソフトでBIDO2.0.tgzを展開してください．次のファイル類が入っているのが分かると思います．

フォルダｄemo\ synth_SN100_18mGamR0.8RV0.1N6にはデモ解析用のsynthetic dataが入っています．

Infoというフォルダはsynthetic dataに関する情報を含んでいます．S0x.dというファイルが６つありますが、これらはそれぞれS0１～S06という名前の６台の地震計を下図のように配置して得られた微動波形データです．とはいってもこれらはデモストレーションのために数値的に合成したsynthetic dataです．サンプリング時間間隔0.01秒で５分間、３成分観測を実施したと想定して合成しました．微動波形の合成にあたっての詳細はこちらをご覧ください．

S01.dを適当なエディターで開くと次のようになっていることが分かると思います．

一番左側が時間、右側へ順に鉛直（ｚ）、EW（ｘ）、NS（ｙ）成分の振幅です．

ご自分の観測データを解析する場合も、３成分の場合は同様に４列で羅列するフォーマットにしてください．水平２成分のデータの場合は鉛直成分を０で埋めるなど、適当にダミーデータを入れて同様のフォーマットにして下さい．一方、鉛直成分しか観測していない場合は、時間、鉛直成分の２列を羅列するようにお願いします（3成分と同じように4列あっても良いのですが、特にダミーデータを入れる必要はないということ）．これは、私の場合これまでに上下動のみの観測が多かったということと、水平動だけの観測経験がない（地震計の仕様により、水平動を観測する場合は必ず鉛直成分も同時に観測していました）という事情によります．

ところで、本プログラムの解析法は等時間間隔で波形をサンプリングすることを前提としています．よって、波形開始時刻とサンプリング時間間隔が与えられる限り、本来、時間の列は不要のはずです（実際のアレイデータ処理では波形開始時刻も不要）．なので、微動アレイを解析したことのある方ならば、「1行目はどうして必要なのかな？」と思われたことでしょう．実は１列目の時間の並びはダミーで、実際は読み込んでいません（したがって、実は１列目にどんな数字が入っていても、解析結果に影響しないのです）．このような列を入れてあるのは、単に、私が観測データを確認するときの描画ツールのフォーマットが時間と振幅の並びを必要としたためです（描画して確認した波形をそのまま解析に用いるのが効率的だったということです）．なお、すぐ後で述べるように、サンプリング時間間隔はseism.dというファイル（seismfile）で与えます．

ともあれ、データがこのように並んでいる限り、データの種類は速度でも加速度でも（単位が何であっても）OKで、数値の書式（指数型、浮動小数点型など）も問いません．数値間のセパレータも、空白、タブ、カンマのいずれかであればＯＫです．データファイル名のつけ方にも特にきまりはありません（拡張子の*.dもなくても構いません）．このように地震計ごとに個別のデータファイルをつくったら、それらをまとめて１つのフォルダにいれて下さい．また、すべてのデータファイルの波形の開始時刻を０に統一して下さい．

同じフォルダには、seism.dというファイルがあります．これを開くと次のように書かれています．

１行目は解析に用いる成分を表す数字、２行目は波形記録のサンプリング時間間隔です（#COMP、#DTはおまじないのようなもので、省略はできません）．３行目以降はｘ、ｙ座標（km）、データファイル名、中心点かどうかのID（中心なら１、円周上なら０）です． seism.dはファイル名の変更は不可です．また観測データと同じフォルダの下に置かなければいけません．ただし、このファイルは必ずしも予め準備しなくても構いません．データフォルダにこのファイルがない場合は、プログラムが必要情報を尋ねて自動的に生成するようになっています．

なおこの例ではｘ、ｙ座標はアレイの中心点を原点としていますが、必ずしもそうする必要はありません．プログラムは円周上に配置された最初の３点の情報を用いてアレイ解析に必要な情報すなわち中心点の位置と半径を計算する仕様としています．更に、それ以外の点がきちんと中心点や円周上にあるかをチェックするようにしています．なお円周上の点は必ずしも等間隔である必要はありませんが（理論は引用文献[2]）、精度を保つためには、できるだけ等間隔のアレイを展開するようにして下さい．円周上に配置されていない点（半径方向の誤差）については、補正理論がないのでディフォルトで５％までの誤差を許容して円周上にあるとみなすことにしています．

BIDO2.0ではハドルテスト（地震計を一箇所に並べて機器特性の一致度を確認するテスト）も実施できるようになりました（ハドルテストのデモ）．ハドルテストをしたい場合、すべてのデータ点に同一のｘ、ｙ座標を入れて下さい．その場合の解析内容は以下の説明とは違っていて、パワースペクトル密度，振幅2乗コヒーレンス，位相差が計算されることになります．

seism.dの説明はまとめてこちらに書いておきますのでご確認下さい．

それではデータセットの確認がすんだところで分析を実行してみましょう． BIDO2.0のフォルダに戻ってBIDO-win.batをダブルクリックしてください．次のような画面（ターミナル）が開きます（Linuxの場合、この操作をする必要はありません）．

このターミナルに「run.sh」と打ち込んでリターンキーをおすと、次のように表示されます．

上の画面以下、ここでは基本的に「y」かリターンキーを打ち込んでみて下さい（ディフォルトパラメータを使う）．最終的に次のような画面に到達すると思います．

「＃＃＃・・・」の行に挟まれた部分に、入力したパラメータが示されていることが分かると思います．ご自分の観測データを解析する場合は、質問事項に関する詳しい説明や、パラメータ設定のコツにも目を通しておいて下さい．

なおこのデモでは３つ目の質問「Preprocessing the data?（データの前処理をしますか？）」に対してディフォルトの回答が「n（いいえ）」となっています．観測データを用いる場合はデータを前処理しておくと効果的な場合もありますが、このデモでは理論的に合成した疑似観測波形（Synthetic data）を用いていますので、ここでは煩雑を避けるため説明を省きます．詳細はこちらをご覧ください．

というわけで、上の画面では今このような値で解析を実行して良いか聞かれているところです．とにかく「y」と打ち込んでみましょう．

するとこのように、現在のアレイジオメトリで解析を実行した時の解析内容が羅列されます．ここで「ｙ」と打ちこめば、あとは自動処理で分析してくれます．

分析開始からしばらくすると、スペクトル推定のデータ処理の前に、波形記録とスペクトル解析にどの部分を使うかという情報が描画されます（下図．波形が赤、スペクトル解析に用いられる部分が緑のプラスマークで図示されている）．

ダイアログ用の端末画面に戻ってリターンキーを押すと、この描画が消えてスペクトル解析に移ります．

解析の実行ログはこのような感じです（実行ログのリンクへ）．このログはデータファイルを置いてあるフォルダに自動生成されます．

分析が終わると、

・ラブ波の位相速度、

・レーリー波の位相速度、

・レーリー波のR／Vスペクトル、

・レーリー波のR／VスペクトルとH/Vスペクトルの比較、

・水平動に占めるレーリー波のパワーの割合、

・NS比（ＳＮ比の逆数）、

・H/Vスペクトル、

・パワースペクトル密度、

が表示されます．それではそれぞれ見ていきましょう．

① ラブ波の位相速度がエラーバー（標準偏差）付きのデータ点で示されています．横軸は周波数、縦軸は位相速度です．

CCA-L、SPAC-L、SPAC+L法という方法で解析した３つの結果が同時に示されています．原点から放射状に広がる直線は、アレイ半径ｒの２倍、５倍、・・・といった波長（参考値）です．１－２Hz程度の帯域で、与えたラブ波の位相速度（synthetic data の詳細へ）が再現されているように見えます．

ところで上のグラフでは、プロットする位相速度の範囲（最大値）は自動調整されています．ご自分で調整したい（固定的に与える）場合、\script\setpar.shのymax_for_gnuplotというパラメータのコメントアウトを外し（行頭の#を消す）、適当な数値（位相速度 [km/s]）を与えて下さい．同様に、周波数の範囲（最大値）をご自分で与えたい場合は\script\setpar.shのxmax_for_gnuplotというパラメータのコメントアウトを外し（行頭の#を消す）、適当な数値（周波数[Hz]）を与えて下さい．

下図は、参考としてsynthetic dataの計算に与えた位相速度（実線）と解析結果の比較を0.４から2.2Hz程度の範囲で拡大したものです．gnuplotでは描画領域の拡大なども簡単です．詳しくはgnuplotのマニュアルページを見てください．

② ターミナルにもどってリターンキーをおしてください．今度は次のような図が描画されます．これはレーリー波の位相速度の解析結果です．

水平動の解析結果としてCCA-R、SPAC-R、SPAC+R法、鉛直動の解析結果としてCCA、nc-CCA（CCA-lwapxと表記されています）、H0、H1、V法の解析結果が同時にプロットされています．0.5－２Hz程度の帯域で、与えたレーリー波の位相速度（synthetic data の詳細へ）が再現されているように見えます．下図は、参考としてsynthetic dataの計算に与えた位相速度（実線）と解析結果の比較を、2Hz程度までの拡大図で示したものです．ｎc-CCA法は最も低周波数まで解析性能があるように見えます．

③ またターミナルに戻ってリターンキーをおすと、今度はレーリー波のR／Vスペクトル（水平動と鉛直動の振幅比）（太実線は平均値、細線は標準偏差）が表示されます．

Synthetic dataの計算時はR／Vスペクトルを周波数にかかわらず0.1としました（synthetic data の詳細へ）．解析結果は１－２Hzを中心とする周波数帯域でその程度の値をとることが分かります．

④ またリターンキーをおすと、レーリー波のR／VスペクトルとH/Vスペクトルの比較が表示されます．このグラフではR/Vスペクトルもパワースペクトルの比として表示されています（ラディアル成分と鉛直成分の振幅比ではない）のでご注意ください．③で述べた通り振幅比としてのR/Vは0. 1ですから、パワー比はその２乗で0.01となります．与えたR/VとH/Vはそれぞれ0.01、 0.0125ですから、この解析結果はその相違を良く表現していると言えると思います．

⑤ またリターンキーをおすと、今度は水平動に占めるレーリー波のパワーの割合γ_Rが表示されます（γ_Rは水平動のパワーを１とした場合のレーリー波のパワーの割合として定義されます．同様にしてラブ波のパワーの割合をγ_Lとすると、γ_R＋γ_L＝１となります）．

Synthetic dataにはγ_R＝0.8を与えました（synthetic data の詳細へ）．解析結果も１－２Hzを中心とする周波数帯域でその程度の値をとることが分かります．

⑥ またリターンキーをおすと、水平動のNS比（SN比の逆数．インコヒーレントノイズのパワーとシグナルのパワーの比）（緑線）と鉛直動のＮＳ比（赤線）が示されます．

Synthetic dataではＳＮ比を１００と設定しました（synthetic data の詳細へ）．NS比はＳＮ比の逆数ですから0.01を与えたことになりますが、実際、解析結果も１Hz以下の低周波数帯域ではその程度の値となっていることが分かります．

またリターンキーをおすと、鉛直動のNS比（赤線）のみが示されます．ただし青とピンクで示された「Upper limit」という線も同時に示されています．

この図でピンクと青の線は、それぞれCCA法、nc-CCA法で同定した位相速度を用いて計算した限界NS比です．限界ＮＳ比とは、CCA法（nc-CCA法ではありません）の解析結果の信頼性を診断するための参考値です．観測ＮＳ比が限界ＮＳ比を上まわると、ＣＣＡ法の解析結果は一般に過小評価されます．この図の場合０．７Hz程度以下で赤＞ピンク（あるいは赤＞青）となっていますから、その帯域ではCCA法の位相速度が過小評価されている可能性が危惧されます．実際、解析結果②の図を見るとそうなっています．一方、nc-ＣＣＡ法はノイズによる過小評価を補正する方法です（ncはnoise compensatedの略）（引用文献[4]）．ＮＳ比の見積もりが正確ならば、解析結果②に見られるようにさらに低周波数の解析性能を伸ばすことができます．

このように、上下動については限界ＮＳ比をプロットしましたが、水平動ではしませんでした。上下動についてはCCA法の適用限界の算出方法が理論的に明らかにされているのに対し（引用文献[３]）、残念ながら水平動では関連する解析手法の適用限界がまだ明らかにされていないからです。

⑦ またリターンキーをおすと、再度H/Vスペクトルが表示されます．今度は周波数軸を対数目盛にして、最大周波数まで表示しています．ここでいう最大周波数とは、\script\setpar.shのfreqmax_aveというパラメータで設定された値（このパラメータは計算時間短縮のために設定されているもので、平均値の算出を設定値までしか行わないようにします．ディフォルトは50[Hz]ですので、必要があればご自分で変更して下さい）とナイキスト周波数のうち小さいほうの値のことです．対数表示をやめたい場合は、\script\setpar.shのautologscale_xというパラメータの行頭に#を記入してコメントアウトして下さい．

⑧またリターンキーをおすと、各記録のパワースペクトル密度が表示されます．これも周波数軸は対数目盛で最大周波数まで表示していますのでご注意ください．図の上に示されたタイトルの右端のカッコ（[ ]）内のu、n、eは上下、南北、東西成分を表しています．凡例の平均、標準偏差の右のカッコ内の数字はseism.dに記述されたデータを上から順にNo.１、２、・・・としたものです．間違いのないように、その横に「S01.d」というふうに、データファイル名が付されています．

以上で分析結果の描画は終わりです．ここでリターンキーをおすと描画画面が消えて、ターミナルには、

・・・と表示されます．これは分析に用いた一時ファイルが消去されて解析が終了したことを示しています．最後のreal、user、sysの横の数字は計算実行にかかった時間を表しています（計算開始からこの画面に到達するまでに実時間で17分かかったということ）．ここで、コマンドプロンプトに「exit」と入力するかC-Z（コントロールキーとｚキーを同時に押す）すれば、ターミナル画面が終了します．

以上で全解析がおわりました．すべての解析結果は、データフォルダの下に作成されるRESULTというフォルダに格納されています．RESULTの中はこのようになっています（各ファイルの内容の説明はこちら）．

この中で、上の描画に用いられたデータファイルはすべてaveというフォルダ（averageの略）の中に入っています．aveの中はこうなっています（各ファイルの内容の説明はこちら）．

aveフォルダの中にあるファイルは、１、２、・・という数字名のフォルダに収納されている個別の解析結果を統計処理したものです．数字名フォルダには、aveに入れるべくして統計処理される解析結果以外に、詳細な解析をしない限り用いないようないろいろな解析結果（作業用の一時ファイルやスペクトル比など）も入っています．例えば、１というフォルダには次のような解析結果が入っています（各ファイル内容の説明はこちら）．

ところで、もとのデータフォルダに戻ると、RESULTという解析結果のフォルダに加えてrun.logという実行ログファイルとparam.shというパラメータファイルが生成されているのが分かると思います．パラメータファイル（param.sh）は、パラメータを少しだけ変更して解析しなおすときにrun.shの引数として与えると便利です．下の図のような感じです．

入力の際、パス名が長いので入力が面倒だったり誤入力が気になったりするかもしれません．しかしここは補完機能という入力支援がありますのでご安心下さい．補完機能とは、途中までアルファベットを打ち込んでTABキーを押せばその後のスペルを補ってくれる機能です．例えば、

run.sh d [TAB]

とすると、自動で

run.sh demo/

としてくれます．その後も同様に補完が効きます．というわけでパラメータファイルの引数付きでrun.shを実行した場合、プログラム開始時のダイアログで今回実行した内容のパラメータをディフォルトとして尋ねてくれます（このように明示的にパラメータが与えられない場合、scriptフォルダの下にあるdefault_param.shを読む仕様になっています）．なお、再解析の際は、RESULTの下に1、2、・・などの数字名フォルダおよびaveフォルダがあれば自動的に消去されますのでご注意ください．これらのフォルダを消去されたくない場合、再解析前に適当な場所に移動するか名前を変えておく必要があります．

最後に、このサンプルでは中心点ありの３成分アレイデータを用いたことから、