frequency_besselコマンド マニュアル

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

◆機能・用途(Purpose)

多数観測点間のクロススペクトルを用いて周波数-ベッセル変換を計算する。
Compute the frequency-Bessel transform for the cross spectra among multiple stations.

◆概要(Overview)

周波数-ベッセル変換はWang et al. (2019)によって開発されたもので、 概要は以下の通りである。
The frequency-Bessel transform was developed by Wang et al. (2019) and is overviewed below.
地表に設置された2つの観測点の座標を \(\posx_1=(x_1,y_1,0)\)、\(\posx_2=(x_2,y_2,0)\)とする。 2観測点間のクロススペクトルは一般に位置\(\posx_1\), \(\posx_2\)に依存するが、 ここではクロススペクトルが観測点の絶対位置にはよらず 相対位置\(\myvector{r}=\posx_2-\posx_1=(\Delta x,\Delta y,0)\) のみによって決まると仮定し、 \(\tilde{C}(\myvector{r},\omega)\)と表す。 ここで\(\omega\)は角周波数である。 このクロススペクトルを\(\Delta x\), \(\Delta y\)に関してフーリエ変換すると \[\begin{eqnarray} I(\omega,k_x,k_y) &=& \frac{1}{2\pi} \int_{-\infty}^{\infty} \Delta x \int_{-\infty}^{\infty} \Delta y \tilde{C}(\myvector{r},\omega) \exp(ik_x \Delta x+ik_y\Delta y) \nonumber \\ &=& \frac{1}{2\pi} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \Delta x \Delta y \tilde{C}(\myvector{r},\omega) \exp(i\myvector{k}\cdot \myvector{r}) \Delta x \Delta y \label{eq.transform.general} \end{eqnarray}\] と書ける。ここで\(\myvector{k}=(k_x,k_y,0)\)は 水平面内の波数ベクトルである。
Let us consider two stations on the ground surface located at \(\posx_1=(x_1,y_1,0)\) and \(\posx_2=(x_2,y_2,0)\). The cross spectrum between the two stations generally depends on the locations \(\posx_1\) and \(\posx_2\); however, we here assume that the cross spectrum does not depend on the absolute station locations but only depends on a relative location \(\myvector{r}=\posx_2-\posx_1=(\Delta x,\Delta y,0)\), and express this cross spectrum as \(\tilde{C}(\myvector{r},\omega)\), where \(\omega\) is an angular frequency. The Fourier transform of this cross spectrum along \(\Delta x\) and \(\Delta y\) is given as Eq. (\ref{eq.transform.general}), where \(\myvector{k}=(k_x,k_y,0)\) is a horizontal wave number vector.
Wang et al. (2019)は\(\tilde{C}(\myvector{r},\omega)\)が 観測点間の方位には依存せずに 距離\(r=|\myvector{r}|\)のみに依存すると仮定した。 この場合、(\ref{eq.transform.general})式の被積分関数は 波数ベクトルの大きさ\(k=|\myvector{k}|\) および\(\myvector{k}\)と\(\myvector{r}\)のなす角\(\theta\)を用いて \(\tilde{C}(r,\omega)\exp(ikr\cos\theta)\)と書けるので、 (\ref{eq.transform.general})式は \[\begin{eqnarray} I(\omega,k) &=& \frac{1}{2\pi} \int_0^{\infty}dr \int_0^{2\pi}r d\theta \tilde{C}(r,\omega)\exp(ikr\cos\theta) \nonumber \\ &=& \frac{1}{2\pi} \int_0^{\infty}dr r\tilde{C}(r,\omega) \int_0^{2\pi}d\theta \exp(ikr\cos\theta) \nonumber \\ &=& \int_0^{\infty}dr r\tilde{C}(r,\omega)J_0(kr) \label{eq.transform} \end{eqnarray}\] と変形できる。ここで \[\begin{equation} J_0(kr)=\frac{1}{2\pi} \int_0^{2\pi}d\theta \exp(ikr\cos\theta) \label{eq.J0} \end{equation}\] はベッセル関数である。 このプログラムでは(\ref{eq.transform})式に基づいて \(I(\omega,k)\)を計算する。
Wang et al. (2019) assumed that \(\tilde{C}(\myvector{r},\omega)\) does not depend on the inter-station azimuth but depends only on the inter-station distance \(r=|\myvector{r}|\). In this case, the integrand of Eq. (\ref{eq.transform.general}) is expressed as \(\tilde{C}(r,\omega)\exp(ikr\cos\theta)\), where \(k=|\myvector{k}|\) is the magnitude of the wave number vector and \(\theta\) is the angle between \(\myvector{k}\) and \(\myvector{r}\). Therefore, Eq. (\ref{eq.transform.general}) is arranged as (\ref{eq.transform}), where \(J_0\) is a Bessel function (Eq. \ref{eq.J0}). This program computes \(I(\omega,k)\) using Eq. (\ref{eq.transform}).
クロススペクトルは複素数であるが、 (\ref{eq.transform})式の計算では実部を使用する。 というのも、観測点\(n\)と\(m\)のクロススペクトルを \(\tilde{C}_{n,m}(\omega)\)とするとき、 観測点\(m\)と\(n\)のクロススペクトルは \(\tilde{C}_{m,n}(\omega)=\tilde{C}_{n,m}(\omega)^{∗}\) である。ここで\(^{∗}\)は複素共役を表す。 (\ref{eq.transform})式の計算では \(\tilde{C}_{n,m}(\omega)\)と\(\tilde{C}_{m,n}(\omega)\)の 両方を用いるべきであり、 この2つのクロススペクトルを平均すれば \(\tilde{C}_{n,m}(\omega)\)の実部のみが残る。 したがって計算にはクロススペクトルの実部を用いれば良く、 \(I(\omega,k)\)も実数となる。
Although cross spectra are complex numbers, the real parts of them should be used for the computation of Eq. (\ref{eq.transform}). This is because the cross spectrum of stations \(m\) and \(n\) (denoted as \(\tilde{C}_{m,n}(\omega)\)) is the complex conjugate of that of stations \(n\) and \(m\); \(\tilde{C}_{m,n}(\omega)=\tilde{C}_{n,m}(\omega)^{∗}\), where \(^{∗}\) indicates a complex conjugate. The computation of Eq. (\ref{eq.transform}) should use both \(\tilde{C}_{n,m}(\omega)\) and \(\tilde{C}_{m,n}(\omega)\). Averaging the two cross spectra, only the real part of \(\tilde{C}_{n,m}(\omega)\) remains. Therefore, the real parts of the cross spectra should be used, and thus \(I(\omega,k)\) is a real number.

◆ソースコード(Source code)

$YMAEDA_OPENTOOL_DIR/structural_survey/src/frequency_bessel.c


◆使用方法(Usage)

コマンドライン引数でパラメータを指定する。 パラメータの一覧を下表に示す。
Specify parameters by command-line arguments. The table below shows a list of parameters.

●「-」から始まらない引数 (Arguments not beginning with “-”)

引数
Argument 		与える値
Quantity to be given
第1引数
1st argument 		読み込むクロススペクトルデータのディレクトリパス。 このディレクトリは cross_spectrum_multiコマンド から出力されたもの、 またはそれと完全に同一のディレクトリ構造およびファイル書式 を持つものでなければならない。
The directory path of the input cross spectral data. This directory must have been created by cross_spectrum_multi command, or must have exactly the same directory structure and file formats as those created by the command.
第2引数
2nd argument 		計算した周波数-ベッセル変換の出力先ファイル名。
The name of the output file for the computed frequency-Bessel transform.


●1つの「-」から始まる引数 (Arguments beginning with a single “-”)

このコマンドでは1つの「-」から始まる引数は存在しない。
This command does not have arguments beginning with a single “-”.

●「--パラメータ名=パラメータ値」の形式の引数 (Arguments of a form “--Parameter name=Parameter Value”)

「--パラメータ名=パラメータ値」の形式の引数は自由な順番で指定できる。 「-」から始まらない引数の間に挿入しても良い。 相反する指定がなされた場合には後の指定が優先される。 デフォルト値を持つパラメータは省略できる。
Arguments of a form “--Parameter name=Parameter Value” can be placed in an arbitrary order. They can even be inserted between arguments not beginning with “-”. In case of conflicting options being specified, the latter option has a higher priority. Parameters that have default values can be omitted.
パラメータ名
Parameter name 		意味
Meaning 		可能なパラメータ値
Allowed parameter values 		デフォルト値
Default value
station_list_file 解析に使用する観測点のリストファイル名。
  • 第1列に観測点コード
  • 第2列にその観測点の成分コード
  • 第3列に観測点の\(x\)座標(基準点から東方向、m)
  • 第4列に観測点の\(y\)座標(基準点から北方向、m)
  • 第5列に観測点の\(z\)座標(標高、m)
を書き並べる。 列の区切りにはタブを用いる。 空行、各行の#から後の部分はコメントとして読み飛ばされるので 自由に挿入できる。
The name of a file that lists the stations used in the analysis. Write
  • a station code in the 1st column,
  • the component code for the station in the 2nd column,
  • the \(x\)-coordinate of the station (east from the reference point, m) in the 3rd column,
  • the \(y\)-coordinate of the station (north from the reference point, m) in the 4th column, and
  • the \(z\)-coordinate of the station (altitude, m) in the 5th column. Use tabs to separate the columns. Empty lines and parts after # in each line are ignored as comments.
    1つの観測点につき1成分のみとすべきであり、 また成分は全観測点で共通のものを用いるべきである。 第2列(成分コード欄)を設けているのは 同じ方向の成分であっても観測点によって U, wU, Z(いずれも上下動を表す)のように 別々のコードが用いられる場合があることを想定したものである。
    Only one component should be used for each station, and the component should be common for all stations. The 2nd column (component code field) is prepared to express the same directional component with different codes (e.g., “U”, “wU” and “Z” for the vertical component).
    観測点の\(z\)座標(標高)は現在のところ使用しないが 将来の拡張のために第5列に入力する仕様にしている。
    Currently, the \(z\)-component (altitude) of each station is not used; the field for this value (5th column) is prepared for future extension.
 ファイル名を表す文字列。 ディレクトリパスを含んでいても良い。
A string that represents a file name, possibly including the directory path. 		省略不可
Cannot be omitted
extension 読み込むクロススペクトルのファイルの書式(拡張子)。
The extension for the input cross spectral file format. 		ymaeda_opentoolsのスペクトルファイル形式のいずれかに対応する拡張子 (独自のファイル形式参照)。
One of the extensions for the spectral file formats of ymaeda_opentools; see special file formats. 		bimseq
independent_variable 波数\(k\)と位相速度\(c\)のどちらを独立変数にするかの選択。
Choice of whether the wave number \(k\) or the phase velocity \(c\) is used as the independent variable.
  • k
    波数\(k\)を独立変数として kminからkmaxまでkincの間隔で動かして計算する。
    Compute the F-J spectrum for the wave number \(k\) as the independent variable, which is varied from kmin to kmax at an increment of kinc.

  • c
    位相速度\(c\)を独立変数として cminからcmaxまでcincの間隔で動かして計算する。
    Compute the F-J spectrum for the phase velocity \(c\) as the independent variable, which is varied from cmin to cmax at an increment of cinc.

 k
kmin \(I(\omega,k)\)を計算する波数\(k\)の最小値(1/m)。
The minimum value of the wave number \(k\) (1/m) to compute \(I(\omega,k)\).
※このパラメータはindependent_variable=kの場合にのみ用いられる。
∗This parameter is used only if independent_variable=k. 		正の実数。
A positive real number. 		パラメータkincの値。
The value of parameter kinc.
kmax \(I(\omega,k)\)を計算する波数\(k\)の最大値(1/m)。
The maximum value of the wave number \(k\) (1/m) to compute \(I(\omega,k)\).
※このパラメータはindependent_variable=kの場合にのみ用いられる。
∗This parameter is used only if independent_variable=k. 		kminよりも大きな正の実数。
A positive real number greater than kmin. 		最小のノンゼロの観測点間距離の逆数の\(\pi\)倍 (Wang et al. (2019)の(14)式)。 但し、kincの整数倍となるように切り上げる。 非常に近接した観測点ペア等によってkmaxが極端に大きくなるのを防ぐため、 kincの1000倍を超える場合はその次に小さい観測点間距離の逆数とする。
The inverse of the smallest non-zero inter-station distance multiplied by \(\pi\) (Eq. 14 of Wang et al. 2019), adjusted to an integer multiple of kinc toward a greater-value side. To prevent kmax from being an extremely large value due to extremely close station pairs, the inverse of the next smallest non-zero inter-station distance is used if kmax exceeds 100 times kinc.
kinc \(I(\omega,k)\)を計算する波数\(k\)の増分(1/m)。
The increment of the wave number \(k\) (1/m) to compute \(I(\omega,k)\).
※このパラメータはindependent_variable=kの場合にのみ用いられる。
∗This parameter is used only if independent_variable=k. 		kmax-kminよりも小さな正の実数。
A positive real number less than the difference between kmax and kmin. 		最大の観測点間距離の逆数の\(2\pi\)倍 (Wang et al. (2019)の(15)式における上限)。
The inverse of the largest inter-station distance (the upper limit of Eq. 15 of Wang et al. 2019).
cmin \(I(\omega,c)\)を計算する位相速度\(c\)の最小値(m/s)。
The minimum value of the phase velocity \(c\) (m/s) to compute \(I(\omega,c)\).
※このパラメータはindependent_variable=cの場合にのみ用いられる。
∗This parameter is used only if independent_variable=c. 		正の実数。
A positive real number. 		パラメータcincの値。
The value of parameter cinc.
cmax \(I(\omega,c)\)を計算する位相速度\(c\)の最大値(m/s)。
The maximum value of the phase velocity \(c\) (m/s) to compute \(I(\omega,c)\).
※このパラメータはindependent_variable=cの場合にのみ用いられる。
∗This parameter is used only if independent_variable=c. 		cminよりも大きな正の実数。
A positive real number greater than cmin. 		10000.0
cinc \(I(\omega,c)\)を計算する位相速度\(c\)の増分(m/s)。
The increment of the phase velocity \(c\) (m/s) to compute \(I(\omega,c)\).
※このパラメータはindependent_variable=cの場合にのみ用いられる。
∗This parameter is used only if independent_variable=c. 		cmax-cminよりも小さな正の実数。
A positive real number less than the difference between cmax and cmin. 		1.0
freq_step クロススペクトルの周波数サンプルのうち、 何サンプルに1つをF-Jスペクトルの計算に使用するか。 周波数サンプルを間引くことで計算時間が短縮される。
The number ratio of total frequency samples in the cross spectrum to the samples used for the computation of F-J spectrum. Resampling the frequencies would fasten the computation. 		正の整数。
A positive integer. 		1
normalize 規格化の方法。
Method of normalization.
  • none
    クロススペクトルを規格化しない。
    Do not normalize the cross spectra.

  • pow2_int
    0 Hzからナイキスト周波数までのクロススペクトルの絶対値2乗積分が 1になるように規格化する。
    Normalize the cross spectrum by the integral of the square of the absolute value of it from 0 Hz to the Nyquist frequency.

  • Nstack
    クロススペクトルをスタック数で割る。
    Divide the cross spectra by the number of stacks.

  • ACF
    自己相関関数のスペクトルを用いてクロススペクトルを規格化する。
    Normalize the cross spectrum by the spectra of auto correlation functions.

  • Nstack_ACF
    クロススペクトルをまずスタック数で割った上で 自己相関関数のスペクトルを用いて規格化する。
    Divide the cross spectra by the number of stacks, and then normalize by the spectra of auto correlation functions.

 Nstack_ACF


◆動作(Behaviour)

第1引数で指定したディレクトリから パラメータstation_list_fileで指定した観測点リストに従って クロススペクトルを読み込み、 \(I(\omega,k)\)を計算する。 結果は第2引数で指定したファイルに以下の書式で出力する。 出力ファイルの列の区切りにはタブが用いられる。
Read the cross spectra in the directory specified by the 1st argument refering to the list of stations given by parameter station_list_file, compute \(I(\omega,k)\), and output the results into the file specified by the 2nd argument with the following format, where tabs are used to separate the columns.
列
Column 		値
Value
1 角周波数\(\omega\) (rad/s)
Angular frequency \(\omega\) (rad/s)
2 周波数\(f=\omega/(2\pi)\) (Hz)
Frequency \(f=\omega/(2\pi)\) (Hz)
3 波数\(k\) (1/m)
Wave number \(k\) (1/m)
4 位相速度\(c=\omega/k\) (m/s)
A phase velocity \(c=\omega/k\) (m/s)
5 \(I(\omega,k)\)
6 将来の拡張のためのスペースであり、現在のバージョンでは0.0が出力される
A space for future extension; 0.0 is written in the current version
7 \(I(\omega,k)\)を 現在の\(\omega\)における\(|I(\omega,k)|\)の最大値で割った値
\(I(\omega,k)\) divided by the maximum value of \(|I(\omega,k)|\) for the current \(\omega\)
8 将来の拡張のためのスペースであり、現在のバージョンでは0.0が出力される
A space for future extension; 0.0 is written in the current version
9 \(I(\omega,k)\)を 全ての\(\omega\), \(k\)における\(|I(\omega,k)|\)の最大値で割った値
\(I(\omega,k)\) divided by the maximum value of \(|I(\omega,k)|\) for all \(\omega\) and \(k\)
10 将来の拡張のためのスペースであり、現在のバージョンでは0.0が出力される
A space for future extension; 0.0 is written in the current version


◆使用例(Example)

cross_spectrum_multi --inputfiles=structural_survey/data/%YYYY/%YY%MM%DD%hh/ 7;STATION.%COMPONENT.sac --trace_list_file=structural_survey/etc/traces.dat --start=2021-01-01.00-00-00 --end=2021-12-31.23-59-59 --outputdir=structural_survey/result2021

frequency_bessel structural_survey/result2021 F-J_result.dat --station_list_file=structural_survey/etc/stations.dat --kmax=kmax=0.01 --kinc=0.0001

ここでstructural_survey/etc/stations.datは 例えば以下のようなテキストファイルである。 ここで[TAB]はタブを表す。
The file structural_survey/etc/stations.dat in this example is a text file. Below shows an example of this file, where [TAB] indicates a tab.
STN1[TAB]U[TAB]123.4[TAB]567.8[TAB]90.1
STN2[TAB]U[TAB]-2345.6[TAB]7890.1[TAB]234.5
STN3[TAB]wU[TAB]6789.0[TAB]-1234.5[TAB]67.8
STN4[TAB]Z[TAB]901.2[TAB]-3456.7[TAB]89.0
STN5[TAB]BHZ[TAB]-12.3[TAB]-456.7[TAB]890.1


◆計算方法(Computation method)

(\ref{eq.transform})式の積分を和で近似する方法については 関数TF_calculate_FJ_varying_kのマニュアル参照。 なお、同じ距離にある観測点組合せについては平均したクロススペクトルを用いる。
An approximation method for the integral of Eq. (\ref{eq.transform}) is described in the documentation of function TF_calculate_FJ_varying_k. Before the computation, the cross spectra are averaged over station pairs at equal distances.

◆検証(Validation)

このプログラムが適切に動作することは CC-FJpy (Li et al. 2021) との比較により検証した。
The output of this program was examined by comparison with that of CC-FJpy (Li et al. 2021).
詳細はこちら (Click here for detail)


◆追加の情報 (Additional information)



◆補足(Additional remarks)

このプログラムではOpenMPによる並列化をサポートしている。 並列化したい場合は環境変数USE_OPENMPの値をyesに変更した上でコンパイルする。
This program supports OpenMP parallelization. To use this option, first modify the environmental variable USE_OPENMP from no to yes, and then compile.
bashでの例(Example for bash)
export USE_OPENMP=yes
cd $YMAEDA_OPENTOOL_DIR/structural_survey
make O3


◆謝辞 (Acknowledgements)

このプログラムの開発・テストは 浅井岬氏(名古屋大学理学部; 2023年12月現在) と共同で行った。
This program was developed and tested as a cooperative work with Misaki Asai (School of Science, Nagoya University; as of December 2023).

◆引用文献 (References)