ymaeda_opentoolsの目的別ガイド

地震波干渉法による地下構造解析

(Guide of ymaeda_opentools by purposes: Structural survey by interferometry)
観測点間の地震波形の相互相関関数は 一方の観測点を発振点、他方の観測点を受振点とするグリーン関数に近似され、 これを用いて観測点間の地下構造に関する情報を引き出すことができます。 このタイプの解析は広く地震波干渉法と呼ばれます。 ymaeda_opentoolsでは地震波干渉法において役立つツールをいくつか用意しており、 それらについて紹介します。
Cross-correlation functions (CCFs) of seismograms between two stations approximate the Green functions from one station as a virtual source to the other station as a receiver, thereby producing the information on the subsurface structure between the two stations. This type of analysis is widely called seismic interferometry. Several tools in ymaeda_opentools may be useful for the seismic interferometry. These tools are introduced below.

◆データの準備 (Preparation of data)

ymaeda_opentoolsでは Seismic Analysis Code (SAC)形式の地震波形データを用いて 地震波干渉法解析を行うことを想定しています。 WIN形式のデータは wintosacコマンド (win_dataサブパッケージ) を用いてSAC形式に変換できます。
The ymaeda_opentools assumes the seismic interferometric analyses to be performed using the seismograms in Seismic Analysis Code (SAC) format. The data in the WIN format can be converted to the SAC one using wintosac command (win_data sub package).
雑微動を用いた地震波干渉法では 地震などの大振幅のイベントの影響を弱めるために 波形の1bit化(正の値を\(+1\)に、負の値を\(-1\)に変換すること)が行われます。 ymaeda_opentoolsではSAC波形データの1bit化のためのプログラムとして sacfile_1bitコマンド (sac_dataサブパッケージ) を用意しています。
In ambient noise-based seismic interferometry studies, effects of large-amplitude events (e.g., earthquakes) are suppressed by 1-bit normalization that replaces positive and negative values by \(+1\) and \(-1\), respectively. In ymaeda_opentools, sacfile_1bit command (sac_data sub package) is available for the 1-bit normalization of the waveform data in a SAC file.
地震波形の欠損部分には多くの場合、ダミーの0.0が入っています。 欠損が多いデータを1bit化した場合、 0.0も\(+1\)や\(-1\)に変換されてしまい、 それが解析結果に悪さをする恐れがあります。 そこでymaeda_opentoolsでは 波形データを\(+1\)と\(-1\)に2値化するのではなく \(+1\), \(0\), \(-1\)に3値化するコマンド ( sacfile_3valuesコマンド (sac_dataサブパッケージ)) も用意しています。
The data defects of seismograms are usually filled with dammy 0.0. The 1-bit normalization for data with abundant defects may result in undesirable results due to conversion of the dammy 0.0 to \(+1\) or \(-1\). To avoid it, sacfile_3values command (sac_data sub package) is also available in ymaeda_opentools, which normalizes the waveform into three values (\(+1\), \(0\), and \(-1\)) instead of the two values (\(+1\) and \(-1\)).
1bit化の代わりに移動時間窓を用いて地震などの大振幅のイベントの影響を弱める方法もあり、 ymaeda_opentoolsの sacfile_normalize_by_moving_aveコマンド を用いてこの処理を実装できます。
An alternative approach to reduce the effects of large-amplitude signals is a normalization by averages in a moving window. This approach is implemented using sacfile_normalize_by_moving_ave command of ymaeda_opentools.
1bit化や移動時間窓と並ぶ地震波干渉法の重要な前処理として 白色化(スペクトルをフラットにすること)があります。 ymaeda_opentoolsでは sacfile_whitenコマンド (sac_dataサブパッケージ) を用いてSAC波形データの白色化を行うことができます。
The spectral whitening, which makes the spectrum white, is another important preprocessing for seismic interferometry. In ymaeda_opentools, the spectral whitening of the waveform data in a SAC file is implemented with sacfile_whiten command (sac_data sub package).

◆相互相関関数の計算 (Computation of the CCF)

ymaeda_opentoolsでは2つのSACファイルが表す波形データの相互相関関数は sacfile_correlationコマンド (sac_dataサブパッケージ) を用いて計算できます。 例えば1つ目の観測点の波形をstation1.sac、2つ目の観測点の波形をstation2.sac とするとき
sacfile_correlation station1.sac station2.sac correlation.sac
により2つの波形の相互相関関数をcorrelation.sacのファイル名で出力できます。
The CCF between two waveform data in SAC files can be computed using sacfile_correlation command (sac_data sub package). For example, the CCF of the waveforms at the first station (station1.sac) and the second station (station2.sac) can be computed by:
sacfile_correlation station1.sac station2.sac correlation.sac
which outputs the CCF as correlation.sac.
相互相関関数の計算時間は観測点数の2乗に比例して増大するため、 観測点数が多くなると全ての観測点ペアについて相互相関係数を計算することが困難になります。 しかし計算の順序を工夫することで観測点数増大が計算時間に及ぼす影響を最小化し、 このような計算を現実的な時間内で実施することが可能になります。 ymaeda_opentoolsではそのためのツールを用意しています。 ツールの紹介の前にまず計算手順の工夫について説明します。 例として観測点数が40点で1年分の相互相関関数を1日毎に求めてスタッキングするケース を考えてみましょう。 この場合、以下が典型的な「悪い」計算手順です。
The computation time of the CCF increases as the square of the number of stations. Therefore, computation of the CCFs for all pairs is difficult when there are many stations. However, the effects of the number of stations on the computation time can be minimized by changing the order of computation, which realizes the computation of the all CCFs. The ymaeda_opentools gives tools to implement this procedure. Before introducing the tools, let me explain the changing order of computation. For example, consider to compute and stack daily CCFs for all pairs of 40 stations over a year. Below is a typical “bad” order of computation.
  1. 全ての観測点ペアについて1日毎の相互相関関数を計算する。 具体的には各観測点ペアについて以下の手順で行う。
    Compute the daily CCFs for every pair of stations, using the following procedures for each station pair.
    1. 2つの観測点の波形をそれぞれフーリエ変換する。
      総計算時間 ∝ 365日×\(_{40}C_2\)ペア=284,700。
      Fourier transform the waveforms at the two stations.
      Total computation time is proportional to 365 days × \(_{40}C_2\) pairs = 284,700.

    2. それらのクロススペクトルを計算する。
      総計算時間 ∝ 365日×\(_{40}C_2\)ペア=284,700。
      Compute the cross spectra of them.
      Total computation time is proportional to 365 days × \(_{40}C_2\) pairs = 284,700.

    3. クロススペクトルをフーリエ逆変換して相互相関関数を得る。
      総計算時間 ∝ 365日×\(_{40}C_2\)ペア=284,700。
      Apply Fourier inverse transformation to the cross spectrum to obtain a CCF.
      Total computation time is proportional to 365 days × \(_{40}C_2\) pairs = 284,700.

  2. 得られたCCFを各観測点ペア毎にスタッキングして1年分の平均のCCFを得る。
    総計算時間 ∝ 365日×\(_{40}C_2\)ペア=284,700。
    Stack the daily CCFs for each station pair to obtain the 1-year average CCF.
    Total computation time is proportional to 365 days × \(_{40}C_2\) pairs = 284,700.

同様の計算結果は以下の(より効率的な)手順で得ることができます。
The same results can be derived by the following (more efficient) procedure.
  1. 全ての観測点ペアについて1日毎のクロススペクトル(周波数領域)を計算する。 具体的には各観測点ペアについて以下の手順で行う。
    Compute the daily cross spectrum (in the frequency domain) for every pair of stations, using the following procedures for each station pair.
    1. 全ての観測点の波形をフーリエ変換する。
      総計算時間 ∝ 365日×40観測点=14,600。
      Fourier transform the waveforms at all stations.
      Total computation time is proportional to 365 days × 40 stations = 14,600.

    2. 得られたフーリエスペクトルを用いて 全ての観測点ペアについてクロススペクトルを計算する。
      総計算時間 ∝ 365日×\(_{40}C_2\)ペア=284,700。
      Compute the cross spectra using the Fourier spectra derived in the previous step.
      Total number of computations: 365 days × \(_{40}C_2\) pairs = 284,700 times.

  2. 得られたクロススペクトルを各観測点ペア毎にスタッキングして 1年分の平均のクロススペクトルを得る。
    総計算時間 ∝ 365日×\(_{40}C_2\)ペア=284,700。
    Stack the daily cross spectra for each station pair to obtain the 1-year average cross spectrum.
    Total computation time is proportional to 365 days × \(_{40}C_2\) pairs = 284,700.

  3. それらをフーリエ逆変換して1年分の平均の相互相関関数を得る。
    総計算時間 ∝ \(_{40}C_2\)ペア=780。
    Apply Fourier inverse transformation to the cross spectra to obtain the 1-year average CCF.
    Total computation time is proportional to \(_{40}C_2\) pairs = 780 times.

前者の手順では全てのステップにおいて 「日数×観測点ペア数」に比例する回数の計算を行っていました。 それに対し、後者の手順では 「日数×観測点ペア数」に比例する回数の計算は ステップ1-iiとステップ2だけになっています。 これらはフーリエスペクトルの単なる掛け算や足し算ですので個々の計算時間が短く、 回数が増えてもあまり問題にはなりません。 時間がかかるフーリエ変換や逆変換の部分に注目すると 後者の手順の方が圧倒的に計算量が少ないことが分かります。
In the former procedures, all computations were repeated with a number proportional to “the number of days × the number of station pairs”. In the latter procedures, only the steps 1-ii and 2 are repeated with this number, and these steps are simply taking the product or summation of a Fourier spectrum that can be done in a short time. The most time-consuming Fourier transformation and inverse transformation have been significantly reduced in the latter procedures.
ymaeda_opentoolsを用いるとこの後者の手順で相互相関関数を計算できます。
Using ymaeda_opentools, the CCFs can be computed this latter procedures.
  1. まず cross_spectrum_multiコマンド を用いて後者の手順のステップ1(1日毎のクロススペクトルの計算)を行います。
    First, the step 1 of the latter procedures (computation of daily cross spectra) is implemented by cross_spectrum_multi command.

  2. 次に stack_cross_spectrum_multi_resultsコマンド を用いて後者の手順のステップ2(クロススペクトルのスタッキング)を行います。
    Second, the step 2 (stacking the cross spectra) is implemented by stack_cross_spectrum_multi_results command.

  3. 最後に imsequencefile_fiftコマンド により後者の手順のステップ3(フーリエ逆変換)を行えば1年平均の相互相互関数が求まります。
    Third, the step 3 (Fourier inverse transformation) is implemented by imsequencefile_fift command to obtain the 1-year averaged CCFs.


◆分散曲線の推定 (Estimating dispersion curves)

多数の観測点ペアについて相互相関関数が求まったら次にそのデータを用いて 周波数と表面波位相速度の関係(分散曲線)を推定します。 ymaeda_opentoolsの frequency_besselコマンド を用いると、分散曲線を推定するアルゴリズムの1つである 周波数-ベッセル変換法(F-J法; Wang et al. 2019)の計算を実装できます。
The step after the CCFs for many pairs of stations is to estimate the relation between the frequency and phase velocity of a surface wave (dispersion curves). This step is implemented using frequency_bessel command of ymaeda_opentools, which is based on the frequency-Bessel (F-J) method proposed by Wang et al. (2019).
F-J法やより古典的なSPAC法(Aki 1957)など、 分散曲線推定のアルゴリズムの多くは クロススペクトルが観測点間距離の関数としてベッセル関数で近似できる という理論に基づいて構築されています。 この近似がよく成り立つデータや周波数であれば位相速度の推定結果に信頼性があると考えられる一方、 この近似が成り立たないデータや周波数では位相速度の推定結果にもあまり信頼性が無いと思われます。 このような観点からの信頼性評価のため、ymaeda_opentoolsの cross_spectrum_distance_fixFreqコマンド を用いるとクロススペクトルを観測点間距離の関数の形に整理できます。
Most algorithms to estimate dispersion curves, including the F-J and more classical SPAC (Aki 1957) methods, rely on a theory that the cross spectra can be approximated by a Bessel function of inter-station distances. The reliability of the estimated phase velocities for each dataset or frequency depends on how well this approximation holds in the actual data. For the evaluation of the reliability of the estimated phase velocities on this viewpoint, cross_spectrum_distance_fixFreq command of ymaeda_opentools reformats the data of cross spectra to a function of distance.

◆鉛直入射するイベント波形を用いた疑似反射法イメージング (Pseudo-reflection imaging by vertically incident event waveforms)

地震波干渉法の一種に鉛直入射するイベント波形の自己相関関数を用いる手法があります。 この手法では観測点の位置で地震波を発生させて同じ地点で観測を行った場合 に相当する記録を得ることができ、 観測点と地下の反射面の間を往復する反射波を取り出せます (Claerbout, 1968)。
Computation of the auto-correlation functions (ACFs) of vertically incident event waveforms is a kind of seismic interferometry. This method yields a record that should be obtained if a seismic wave was generated at a the station site and observed the resultant waveforms at the same location, which consists of the waves reflected between the station and a subsurface reflector (Claerbout, 1968).
この解析はymaeda_opentoolsの reflection_response_eventコマンド (structural_surveyサブパッケージ) を用いて行うことができます。 イベント波形のSACファイルと初動の到着時刻のリストを前もって作成しておけば、 個々のイベント波形の自己相関関数に基づく反射応答の計算からそれらのスタッキングまで コマンド一発で行うことが可能です。 Maeda and Watanabe (2022)の手法に基づき、 反射応答の誤差レベルも同時に求めることができます。
This analysis can be conducted using the reflection_response_event command (structural_survey sub package). What you need to prepare are the SAC files for the event waveforms and a list of arrival times. Then a single run of the command gives the reflection responses from the ACFs of individual waveforms and their stacks. Errors in the reflection responses are also computed based on the method of Maeda and Watanabe (2022).
反射応答は時間(地表と地下の反射面の間の往復走時)の関数になっています。 これを反射面の深さの関数に直すには twoWayTraveltime2depthコマンド (structural_surveyサブパッケージ) を使用します。
The reflection response is a function of time, which is the round-trip time between the ground surface and a subsurface reflector. The time can be converted to the depth of the reflector using the twoWayTraveltime2depth command (structural_survey sub package).
この深さ断面データは往復走時に関して一定刻みになっているため 深さの刻みは一定ではありません。 interpolate_depthSectionコマンド (structural_surveyサブパッケージ) を用いると補間により深さの刻みを一定にすることができます。 これにより、断面図のプロットや同じ標高での値を差し引く等の処理が容易になります。 特に、観測点の標高を無視できない山岳地域での解析の場合、 観測点間でデータサンプルの標高を揃えることで様々な事後処理がやりやすくなります。
This depth-section data is uniformly sampled for the round-trip time, and thus the interval of the data samples is not uniform along the depth. The interpolate_depthSection command (structural_survey sub package) interpolates the data to make the depth interval uniform. This easens some of the postprocessings, for example plotting a cross section or subtracting the average amplitude at the same altitude over stations. Especially, in mountain regions where station altitudes are not negligible, making the sample altitudes consistent among the station is essential to various postprocessings.
推定した地下構造の妥当性の検証や反射面の解釈に当たっては、 地下構造を与えて理論計算により波形を生成するアプローチが 有効な場合があります。 FDM_1D_verticalコマンド (opentwsサブパッケージ) を用いれば、1次元構造内での地震波の鉛直伝播を 反射波を含めて計算することができます。
To examine the validity of the estimated structure or to interpret the reflectors, numerical simulation of the seismic wave for a given subsurface model is sometimes effective. The FDM_1D_vertical command (opentws sub package) simulates the vertical propagation of the seismic wave through 1-D structure including reflected waves.