５．プログラムの内容とアルゴリズム
●プログラムの内容
\BIDO-win.bat
WindowsでMSYSを起動するためのバッチファイルです。Linuxでは用いられません。
\bin
\bin\winbinにWindows用実行ファイル、\bin\linbinにLinux用実行ファイルが入っています。\bin直下にはMSYS用実行ファイルが入っています。
\demo デモデータが入っています（WindowsとLinuxで共通）。
\etc MSYS起動用のスクリプトが入っています。Linuxでは用いられません。
\run.sh プログラム起動用のＢシェルスクリプトです（WindowsとLinuxで共通）。
\script フォートランプログラムを連結するためのBシェルスクリプトが入っています（WindowsとLinuxで共通）。
\src フォートランプログラムが入っています（WindowsとLinuxで共通）。
●アルゴリズム
以下に、全体の流れと個々の手順を説明します。括弧内は関係するＢシェルスクリプトとフォートランプログラムのファイル名です。
全体の流れ(\script\circle.sh)
①
利用可能なデータ部分を抽出する。
②
円周データを方位平均する（①の結果はここでは特に用いていない）
③
①、②の結果を用いてスペクトルを推定する。
④
スペクトル比・位相速度・ＮＳ比等を推定する。
⑤
セグメントクラスター数だけ③、④の処理を繰り返す。
⑥
④の結果を用いて平均、標準偏差を計算する。
⑦
⑥の結果を描画する。
① セグメントの自動抽出(\script\mksegment.sh／\src\evalrms.F、\src\segment.F)
以下の手順でセグメントを抽出します。
１）全地震計の全成分についてそれぞれ以下の作業を実行する。
・オリジナルの波形データから線形トレンドを差し引いてRMSを計算する（この値をRMSallとする）。
・オリジナルの波形データをセグメント長で分割したデータ各部分（各部分は半分ずつ重なり合うようにする）について、線形トレンドを差し引き、RMS値を計算する。各RMS値をRMSallで規格化する（RMSsegとする）。
２）全地震計の全成分のRMSseg について、0.1刻みでRMSsegのヒストグラムを作る。このヒストグラムから最頻値を同定する。
３）全地震計の全成分のRMSsegが同時に最頻値になるデータセグメントを利用可能なセグメントとみなしてsegmentファイルに書き込む。
② 円周データの方位平均(\script\mkcrcle*.sh／\src\mkcrcl_*.F)
本プログラムで採用する手法はすべて円周上の記録を重みつき方位平均することから始まります。重みつき方位平均は、円周沿いのフーリエ級数展開におけるフーリエ係数の計算に対応します。本プログラムではディフォルトで0、1次のフーリエ係数を計算するようになっています。円周沿いの地震計配置は理論的には不等間隔にも対応しており、CCA法では実際の適用性も検討済みです（引用文献[1]）。しかしすべての方法で適用性を詳細に検討したわけではありませんので、できるだけ等間隔のアレイを用いることをお勧めします（特にクロススペクトル密度を用いる方法の適用性には不安が残ります）。
③ スペクトルの推定 (\script\estspec*.sh／\src\estspec.F)
本プログラムでは、セグメントごとにトレンドを除去してcos型のテーパー（テーパー率0.5（Carter et al, 1973））をかけたあと、FFTによる生スペクトルのセグメント平均とスペクトルの平滑化の両方を用いてスペクトルを推定します（Bendat & Piersol, 1971）。生スペクトルをParzen
windowで平滑化した後、セグメント平均を実施しています。セグメント平均の際のセグメント長、セグメント数はそれぞれプログラムの起動時のダイアログ“Duration of data segments for the evaluation of spectra”、“Number of data segments over which averages are taken. Enter 0 or a very large number if you
wish to use all segments
simultaneously”で与えられます。Parzen windowによる平滑化のバンド幅は“Band width of the Parzen spectral window”というダイアログで与えられます。
④
スペクトル比の計算 (\script\specratio*.sh)
上の手順で推定されたスペクトル密度を用いて、分母が０になる場合のみ除外して単純に比をとっています。いろいろな種類のスペクトル比が引用文献[2、5]の定式によって位相速度に結び付けられます。SPAC法の自己相関係数も引用文献[2]の定式に従ってスペクトル比により定義されていますのでご注意下さい（通常の実務では、SPAC係数は中心点と円周上の点の複素コヒーレンスを取った上で方位平均することにより計算される）。
H/Vスペクトルの計算
Version 1.2.2から，三成分データの場合seismファイルの先頭に書かれた１点のＨ/Vスペクトルを推定する機能を追加しました（H/Vは水平成分と鉛直成分のパワースペクトルの比．水平成分のパワーはNS成分とＥＷ成分のパワーの和で定義）．したがってseismファイルに記述された総観測点数がたとえ１点であっても（アレイデータでなくても），３成分あればH/Vスペクトルの計算はなされます．計算終了後，対数平均と標準偏差が描画されます．
位相速度の計算 (\script\spec2pv*.sh／\src\sctr2pv.F)
定式化にしたがってスペクトル比をベッセル関数で表した上で、２分法と割線分割法の組み合わせによる求解法（Shampine & Watts、1970）でベッセル関数の引数であるrk（半径ｘ波数）を 、[0, rkmax]の範囲で探索します。rkmaxは関数値の最初の極大あるいは極小（または最大、最小）に対応します。得られたrkから位相速度c＝2πf／k （fは周波数）を計算します。
⑤
セグメントクラスター数だけ③、④の処理を繰り返す。
⑥平均・標準偏差(\script\mkave.sh／\src\calave.F)
セグメント数（セグメント平均法でスペクトルを推定する際に平均化に用いるセグメントの数。またはプログラム開始時のダイアログ“Number of data segments over which averages are taken. Enter 0 or a very large number if you
wish to use all segments
simultaneously”で入力する整数値）が
（segmentファイルに記載された全セグメント数）＞（セグメント数）
を満たすならば、与えられた波形データから複数のスペクトルを推定できることになります。ここで、
（セグメントクラスター数）＝（segmentファイルに記載された全セグメント数）÷（セグメント数）
と定義すると、スペクトルやそれに派生して計算できるスペクトル比、位相速度などは、セグメントクラスター数だけ推定値が存在することになります（割り切れない分は切り捨てます）。よって、本プログラムではこれらの推定値から平均値や標準偏差を計算しています。
プログラムの実行時はデータフォルダの下に\RESULT\数字のフォルダが生成されますが、この数字は各セグメントクラスターを表していて、これらの数字フォルダの下には対応する解析結果が入っています。\RESULT\aveフォルダの下に、これらの数字フォルダの解析結果を統計処理したファイルが入っています。
なお\src\PARAM.hに記載されているNROBUST4AVERAGE_INCよりもセグメントクラスター数が多い場合に限り、「外れ値（outlier）」を除外する目的で統計処理から最大値と最小値が除外されます。ディフォルトではNROBUST4AVERAGE_INCは８と設定されています。また、解析結果の種類によって単純平均（算術平均）と対数平均を使い分けています。対数平均をとるのはAmpＲV_R.d、nsr.d、nsrlim_cca.d、nsrlim_cca.lwapx.d、powratio_R2L.dというような、比の平均です。
Bendat, J. S., and A. G.
Piersol, Random Ｄata: Analysis and Ｍeasurement Procedures: John Wiley & Sons, 1971.
Carter, G. C., C. H. Knapp, and A. H. Nuttall, 1973,
Estimation of the magnitude-squared coherence function via overlapped Fast
Fourier Transform processing: IEEE Transactions on Audio Electroacoustics, AU-21,
337–344.
Shampine, L. F., and H. A. Watts,
FZERO, a root-solving code, Report SC-TM-70-631, Sandia Laboratories,1970.