frequency_besselコマンドの検証

(Validation of frequency_bessel command)
Last Update: 2023/12/15

frequency_besselコマンドが適切に動作することは CC-FJpy (Li et al. 2021) との比較により検証した。
The output of frequency_bessel command was examined by comparison with that of CC-FJpy (Li et al. 2021).

◆使用したデータ (Data used)

実際のデータを用いた検証はデータサイズが大きく大変である。 そこで、プログラムの検証には以下に述べる小さなデータを作成して使用した。
The examination using an actual data set is numerically expensive because of huge data quantity. Therefore, the examination of the program used a small data set that was made as follows.

●観測点 (Stations)

観測点はSTN1-STN4の4つとし、 (0,0), (1000,0), (0,1000), (1000,2000) の4地点に配置した。 以下が観測点リストファイルである。 ここで[TAB]はタブを表す。
Four stations from STN1 to STN4 were deployed at (0,0), (1000,0), (0,1000), (1000,2000). The corresponding station list file is shown below, where [TAB] indicates a tab.
(sample_station_list_for_check.dat)
STN1[TAB]U[TAB]0.0[TAB]0.0[TAB]0.0
STN2[TAB]U[TAB]1000.0[TAB]0.0[TAB]0.0
STN3[TAB]U[TAB]0.0[TAB]1000.0[TAB]0.0
STN4[TAB]U[TAB]1000.0[TAB]2000.0[TAB]0.0


●クロススペクトルのデータ (Cross spectral data)

個々のクロススペクトルは0.0 Hzから0.7 Hzまで0.1 Hz刻みの8サンプルとした。 なお、0.4 Hzをナイキスト周波数とし、 0.5-0.7 Hzは0.1-0.3 Hzのスペクトルの折り返しの複素共役とした。 同一観測点のクロススペクトルは自己相関関数のフーリエスペクトルであるので 虚部を0とした。
The cross spectra were sampled at every 0.1 Hz from 0.0 to 0.7 Hz (eight samples for each station pair). The Nyquist frequency was assumed to be 0.4 Hz; the spectra in 0.5-0.7 Hz were the complex conjugates of the spectra in 0.1-0.3 Hz with the reversed order. The imaginary parts for the cross spectra between the same station were assumed to be zero because the spectra were the Fourier spectra of auto correlation functions.
作成したクロススペクトルのデータを以下に示す。
The cross spectral data used in the examination are shown below.
周波数(Hz)
Frequency (Hz) 		クロススペクトル
Cross spectrum
STN1-STN1 STN2-STN2 STN3-STN3 STN4-STN4
0.0 1.2 1.23 9.8 9.87
0.1 3.4 4.56 7.6 6.54
0.2 5.6 7.89 5.4 3.21
0.3 7.8 0.12 3.2 0.98
0.4 9.0 3.45 1.0 7.65

周波数(Hz)
Frequency (Hz) 		クロススペクトル
Cross spectrum
STN1-STN2 STN1-STN3 STN1-STN4
0.0 1.2 1.23 1.234
0.1 3.4+5.6i 4.56+7.89i 5.678-9.012i
0.2 7.8+9.1i 0.12-3.45i 3.456+7.891i
0.3 2.3+4.5i 6.78-9.01i 2.345-6.78i
0.4 6.7 2.34 9.01

周波数(Hz)
Frequency (Hz) 		クロススペクトル
Cross spectrum
STN2-STN3 STN2-STN4 STN3-STN4
0.0 9.8 9.87 9.876
0.1 7.6+5.4i 6.54+3.21i 5.432-1.098i
0.2 3.2+1.9i 0.98-7.65i 7.654+3.219i
0.3 8.7-6.5i 4.32+1.09i 8.765+4.321i
0.4 4.3 8.76 0.987

クロススペクトルのファイルはsample_data_for_checkサブディレクトリの下に imseq1形式で作成した。 全ては示さないが、STN1とSTN2のクロススペクトルファイルを 以下に例として示す。 他の観測点ペアに対するファイルも同じ要領で作成できる。
Files for the cross spectra were created as imseq1 format under sample_data_for_check sub directory. The file for the cross spectrum between STN1 and STN2 is shown below. Files for the other station pairs can readily be created in the same manner.
(sample_data_for_check/STN1.U_STN2.U.imseq1)
size=8
t0=0.0
dt=0.1

1.2[TAB]0.0
3.4[TAB]5.6
7.8[TAB]9.1
2.3[TAB]4.5
6.7[TAB]0.0
2.3[TAB]-4.5
7.8[TAB]-9.1
3.4[TAB]-5.6


●スタック数のデータ (Data for the number of stacks)

スタック数のファイルは以下のようにした。
A file for the number of stacks was as below.
(sample_data_for_check/Nstack.dat)
STN1[TAB]U[TAB]STN1[TAB]U[TAB]123
STN1[TAB]U[TAB]STN2[TAB]U[TAB]45
STN1[TAB]U[TAB]STN3[TAB]U[TAB]67
STN1[TAB]U[TAB]STN4[TAB]U[TAB]89
STN2[TAB]U[TAB]STN2[TAB]U[TAB]234
STN2[TAB]U[TAB]STN3[TAB]U[TAB]12
STN2[TAB]U[TAB]STN4[TAB]U[TAB]34
STN3[TAB]U[TAB]STN3[TAB]U[TAB]345
STN3[TAB]U[TAB]STN4[TAB]U[TAB]56
STN4[TAB]U[TAB]STN4[TAB]U[TAB]456

◆frequency_besselコマンドの実行 (Executing frequency_bessel command)

上記のファイルを作成後、以下のコマンドを実行した。
After creating the files described above, the following command was executed.
frequency_bessel sample_data_for_check result_ymaeda_opentools.dat --station_list_file=sample_station_list_for_check.dat --extension=imseq1 --independent_variable=c --cmin=1000.0 --cmax=4000.0 --cinc=1000.0
これにより、位相速度を独立変数として1000 m/sから4000 m/sまで 1000 m/sの刻みで動かして計算した\(I(\omega,k)\)が ファイルresult_ymaeda_opentools.datに出力される。
By this, \(I(\omega,k)\) values computed for phase velocities from 1000 m/s to 4000 m/s at 1000 m/s intervals are written into a file result_ymaeda_opentools.dat.

◆CC-FJpyによる計算 (Computation by CC-FJpy)

●CC-FJpyで必要なデータ1: 観測点間距離のリスト (Data needed by CC-FJpy 1: list of inter-station distances)

CC-FJpyでは観測点の座標ではなく観測点間距離を リストで与える方式になっている。 したがって、観測点間の距離を前もって計算しておく必要がある。 今回のテストで使用する観測点について距離を計算すると下表のようになる。
The CC-FJpy requires a list of inter-station distances instead of station coordinates. Therefore, the inter-station distances must be calculated in advance to using CC-FJpy. The table below shows the inter-station distances in this test.
STN1 STN2 STN3 STN4
(0,0) (1000,0) (0,1000) (1000,2000)
STN1 (0,0) 0 1000 1000 2236.067977
STN2 (1000,0) 1000 0 1414.213562 2000
STN3 (0,1000) 1000 1414.213562 0 1414.213562
STN4 (1000,2000) 2236.067977 2000 1414.213562 0

これを観測点距離順に並べると以下のようになる。
This list is sorted based on the distances as below.
距離(m)
Distance (m) 		観測点ペア
Station pairs
0 STN1-STN1, STN2-STN2, STN3-STN3, STN4-STN4
1000 STN1-STN2, STN1-STN3
1414.213562 STN2-STN3, STN3-STN4
2000 STN2-STN4
2236.067977 STN1-STN4


●CC-FJpyで必要なデータ2: 前処理を済ませたクロススペクトル (Data needed by CC-FJpy 2: preprocessed cross spectra)

frequency_besselコマンドを用いた上記のテストでは
  1. クロススペクトルをスタック数で割る
  2. 更にその結果を自己相関関数のフーリエスペクトルを用いて規格化する
  3. 同じ距離の観測点ペアについてはクロススペクトルを平均する
という前処理を行った上で\(I(\omega,k)\)を計算している。 ここでのテストの目的はCC-FJpyにおける前処理の仕方が これと同じになっているかを調べることではなく、 frequency_besselコマンドにおける\(I(\omega,k)\)の計算自体が 正しく行えているかをチェックすることである。 そのため、上と同様の前処理を手作業で行ったデータを CC-FJpyに渡すことにする。
The test using the frequency_bessel command above preprocessed the data as follows:
  1. divide each cross spectrum by the number of stacks;
  2. normalize the result by the Fourier spectra of auto correlation functions; and
  3. average the cross spectra at equal distances.
The purpose of this test is not to compare the preprocessing procedure to that of CC-FJpy but to check the computation of \(I(\omega,k)\) in the frequency_bessel command. Therefore, the data were manually preprocessed in the same manner as above before passing them to CC-FJpy.
まずクロススペクトルをスタック数で割ると以下のようになる。 ここで表の観測点コードの欄にある()内は分母に用いたスタック数を表す。
First, dividing each cross spectrum by the number of stacks yields the following data, where the value in () in each the station code field represents the number of stacks used as the denominator.
周波数(Hz)
Frequency (Hz) 		クロススペクトル/スタック数 (\(\times 10^{-3}\))
Cross spectrum divided by the number of stacks (\(\times 10^{-3}\))
STN1-STN1 (123) STN2-STN2 (234) STN3-STN3 (345) STN4-STN4 (456)
0.0 9.756098 5.256410 28.405797 21.644737
0.1 27.642276 19.487179 22.028986 14.342105
0.2 45.528455 33.717949 15.652174 7.039474
0.3 63.414634 0.512821 9.275362 2.149123
0.4 73.170732 14.743590 2.898551 16.776316

周波数(Hz)
Frequency (Hz) 		クロススペクトル/スタック数 (\(\times 10^{-3}\))
Cross spectrum divided by the number of stacks (\(\times 10^{-3}\))
STN1-STN2 (45) STN1-STN3 (67) STN1-STN4 (89)
0.0 26.666667 +0.000000i 18.358209 +0.000000i 13.865169 +0.000000i
0.1 75.555556 +124.444444i 68.059701 +117.761194i 63.797753 -101.258427i
0.2 173.333333 +202.222222i 1.791045 -51.492537i 38.831461 +88.662921i
0.3 51.111111 +100.000000i 101.194030 -134.477612i 26.348315 -76.179775i
0.4 148.888889 +0.000000i 34.925373 +0.000000i 101.235955 +0.000000i

周波数(Hz)
Frequency (Hz) 		クロススペクトル/スタック数 (\(\times 10^{-3}\))
Cross spectrum divided by the number of stacks (\(\times 10^{-3}\))
STN2-STN3 (12) STN2-STN4 (34) STN3-STN4 (56)
0.0 816.666667 +0.000000i 290.294118 +0.000000i 176.357143 +0.000000i
0.1 633.333333 +450.000000i 192.352941 +94.411765i 97.000000 -19.607143i
0.2 266.666667 +158.333333i 28.823529 -225.000000i 136.678571 +57.482143i
0.3 725.000000 -541.666667i 127.058824 +32.058824i 156.517857 +77.160714i
0.4 358.333333 +0.000000i 257.647059 +0.000000i 17.625000 +0.000000i

次にクロススペクトルを自己相関関数のフーリエスペクトルで規格化する。 数式で表すと以下の処理になる。 観測点\(n\), \(m\)のクロススペクトルを\(C_{nm}(\omega)\)として、 これを\(\sqrt{|C_{nn}(\omega)|^2|C_{mm}(\omega)|^2}\)で割る。 周波数0.1 HzにおけるSTN1とSTN2の例でいえば上の表より \(C_{nm}(\omega)=75.555556 +124.444444i\), \(C_{nn}(\omega)=27.642276\), \(C_{mm}(\omega)=19.487179\) であるので \((75.555556 +124.444444i)/\sqrt{27.642276\times 19.487179}\) を計算することになる。 このような計算を行うと下表の結果が得られる。
Next, normalize the cross spectra by the Fourier spectra of auto correlation functions as follows. Let \(C_{nm}(\omega)\) be the cross spectrum between \(n\)th and \(m\)th stations. The normalization is realized by dividing \(C_{nm}(\omega)\) by \(\sqrt{|C_{nn}(\omega)|^2|C_{mm}(\omega)|^2}\). For example, consider the cross spectrum between STN1 and STN2 at 0.1 Hz. The table above shows \(C_{nm}(\omega)=75.555556 +124.444444i\), \(C_{nn}(\omega)=27.642276\), and \(C_{mm}(\omega)=19.487179\). Therefore, the normalization is expressed as \((75.555556 +124.444444i)/\sqrt{27.642276\times 19.487179}\). The normalized cross spectra are as below.
周波数(Hz)
Frequency (Hz) 		規格化したクロススペクトル
The normalized cross spectrum
STN1-STN1 STN2-STN2 STN3-STN3 STN4-STN4
0.0 1 1 1 1
0.1 1 1 1 1
0.2 1 1 1 1
0.3 1 1 1 1
0.4 1 1 1 1

周波数(Hz)
Frequency (Hz) 		規格化したクロススペクトル
The normalized cross spectrum
STN1-STN2 STN1-STN3 STN1-STN4
0.0 3.723797 +0.000000i 1.102780 +0.000000i 0.954137 +0.000000i
0.1 3.255405 +5.361843i 2.758074 +4.772194i 3.204144 -5.085549i
0.2 4.423948 +5.161273i 0.067093 -1.928926i 2.169063 +4.952568i
0.3 8.962680 +17.535678i 4.172484 -5.544850i 2.256979 -6.525508i
0.4 4.533065 +0.000000i 2.398181 +0.000000i 2.889467 +0.000000i

周波数(Hz)
Frequency (Hz) 		規格化したクロススペクトル
The normalized cross spectrum
STN2-STN3 STN2-STN4 STN3-STN4
0.0 66.833886 +0.000000i 27.215571 +0.000000i 7.112355 +0.000000i
0.1 30.567548 +21.719047i 11.505828 +5.647356i 5.457177 -1.103089i
0.2 11.607824 +6.892146i 1.870882 -14.604337i 13.020960 +5.476152i
0.3 332.422086 -248.361329i 121.029578 +30.537556i 35.056442 +17.282246i
0.4 54.814460 +0.000000i 16.382315 +0.000000i 2.527495 +0.000000i

最後に、同じ距離の観測点ペアについて平均を取ると 下表のデータが得られる。 結果が実数になるのは観測点を入れ替えたペアとの平均を取るためである。 例えば観測点STN2とSTN1のクロススペクトルは 観測点STN1とSTN2のクロススペクトルの複素共役になっている。 したがって、この2つのクロススペクトルの平均を取れば 観測点STN1とSTN2のクロススペクトルの実部のみが残る。
Finally, averaging the cross spectra at equal distances results in the data shown in the table below. The results are real numbers because the cross spectrum of each station pair is averaged with that of its reciprocity pair. For example, the cross spectrum between stations STN2 and STN1 is the complex conjugate of the cross spectrum between stations STN1 and STN2. Averaging the two cross spectra results in the real parts of the cross spectrum between STN1 and STN2.
周波数(Hz)
Frequency (Hz) 		平均したクロススペクトル
The averaged cross spectrum
0 m 1000 m 1414.213562 m 2000 m 2236.067977 m
0.0 1 2.413289 36.973121 27.215571 0.954137
0.1 1 3.006740 18.012362 11.505828 3.204144
0.2 1 2.245521 12.314392 1.870882 2.169063
0.3 1 6.567582 183.739264 121.029578 2.256979
0.4 1 3.465623 28.670977 16.382315 2.889467


●CC-FJpyを用いるスクリプト (A script that uses CC-FJpy)

試行錯誤の結果、上記のデータを用いたCC-FJpyによる計算は 下記のスクリプトで動かすことができた。 ここでrは観測点間距離のリスト、cは位相速度のリスト、fは周波数のリストである。 ufがクロススペクトルであり、 周波数のループを内側、観測点間距離のループを外側とする2次元配列で与える。 赤字部分はAnacondaとCC-FJpyのインストールディレクトリパスで置き換える。
After several trials, the computation by CC-FJpy using the aforementioned data succeeded by the following script, where r, c, and f are the lists of inter-station distances, phase velocities, and frequencies, respectively. The cross spectra are expressed by uf; it is a 2-D array, where the frequency is the inner loop and the inter-station distance is the outer loop. The red part should be replaced by the installation directory paths of Anaconda and CC-FJpy.
(test_CC-FJpy.py)
#! /usr/local/anaconda/bin/python3

import sys
import numpy as np
sys.path.append(’/usr/local/CC-FJpy/CC-FJpy-master’)
import ccfj

r=np.array([0,1000,1414.213562,2000,2236.067977])
c=np.array([1000,2000,3000,4000])
f=np.array([0.0,0.1,0.2,0.3,0.4])

uf=np.array([[1, 1, 1, 1, 1],
[2.413289, 3.006740, 2.245521, 6.567582, 3.465623],
[36.973121, 18.012362, 12.314392, 183.739264, 28.670977],
[27.215571, 11.505828, 1.870882, 121.029578, 16.382315],
[0.954137, 3.204144, 2.169063, 2.256979, 2.889467]])

FJ=ccfj.fj_noise(uf,r,c,f,fstride=1,itype=1,func=0,num=20)

print(FJ)

◆結果の比較 (Comparison of the results)

frequency_besselコマンドによる計算結果 (result_ymaeda_opentools.dat) の第2列(周波数)、第4列(位相速度)、第5列(\(I(\omega,k)\))、 第7列(\(I(\omega,k)\)を周波数毎の絶対値最大値で割った値) を取り出し、CC-FJpyからの出力と比較すると下表のようになる。
The computation result by frequency_bessel command (result_ymaeda_opentools.dat) is shown below, where only the 2nd column (frequency), 4th column (phase velocity), and 5th column (\(I(\omega,k)\), and 7th column (\(I(\omega,k)\) divided by its absolute maximum of the same frequency) are shown. They are compared with the output from CC-FJpy.
周波数(Hz)
Frequency (Hz) 		位相速度(m/s)
Phase velocity (m/s) 		frequency_besselコマンド(ymaeda_opentools)の出力
Output from frequency_bessel command (ymaeda_opentools) 		CC-FJpyの出力
Output from CC-FJpy
\(I(\omega,k)\) \(I(\omega,k)\)を周波数毎の絶対値最大値で割った値
\(I(\omega,k)\) divided by its absolute maximum in each frequency
0.0 1000 -nan 0.000000e+00 nan
0.0 2000 -nan 0.000000e+00 nan
0.0 3000 -nan 0.000000e+00 nan
0.0 4000 -nan 0.000000e+00 nan
0.1 1000 1.862966e+07 7.625144e-01 0.76250154
0.1 2000 2.320123e+07 9.496291e-01 0.94961363
0.1 3000 2.410934e+07 9.867982e-01 0.9867756
0.1 4000 2.443189e+07 1.000000e+00 1
0.2 1000 3.595929e+06 3.052127e-01 0.30521154
0.2 2000 9.760579e+06 8.284514e-01 0.82844895
0.2 3000 1.123685e+07 9.537534e-01 0.95374817
0.2 4000 1.178172e+07 1.000000e+00 1
0.3 1000 -5.640023e+07 -2.842333e-01 -0.28423303
0.3 2000 1.108456e+08 5.586147e-01 0.5586142
0.3 3000 1.733517e+08 8.736189e-01 0.87361753
0.3 4000 1.984294e+08 1.000000e+00 1
0.4 1000 -1.005777e+07 -3.652991e-01 -0.36529896
0.4 2000 7.702654e+06 2.797610e-01 0.27976096
0.4 3000 2.138408e+07 7.766713e-01 0.77667063
0.4 4000 2.753298e+07 1.000000e+00 1

この比較から CC-FJpyの出力は\(I(\omega,k)\)を周波数毎の絶対値最大値で割った値 であることが分かる。 frequency_besselコマンドとCC-FJpyの結果は おおよそ有効数字4桁の精度で一致している。 ややずれが大きいのはCC-FJpyの中で単精度実数が用いられていることが 影響している可能性がある。 いずれにせよ、この程度のずれであれば frequency_besselコマンドからの出力は 分散曲線推定に十分な精度があるとみなせる。
The comparison shows that the outputs of CC-FJpy are \(I(\omega,k)\) divided by its absolute maximum in the same frequency. The results from frequency_bessel command and CC-FJpy are consistent to an accuracy of ~4 effective digits. The relatively large difference may have been due to the single precision real numbers used in CC-FJpy. Anyway, the difference is small enough to use the output of frequency_bessel command for estimation of dispersion curves.