相互相関関数の計算

西田究のホームページ (u-tokyo.ac.jp)
解説 (u-tokyo.ac.jp)

地震動に関する様々なシミュレーション結果をお手軽に見ることができる、素晴らしいサイトです。特に、地震波干渉法で相互相関関数が出来上がっていく過程を見たことがなかったので、非常に興味深く見ていました。

　　初期

数十分後

相互相関関数を求める際は、時間領域（直接法）か周波数領域（FFT）を利用します。配列の要素数が多ければ周波数領域での計算が早いでしょう。Python の scipy.signal.correlate では手法を指定できますし、指定しなくとも早い方を自動で選択してくれます（Ver.1.13.0）。 numpy.correlate は直接法のみです（Ver.1.26）。

以下、計算例です。設定次第なのですが、以下だと答えは①≠②=③
①直接法

sig1 = np.std(data1)
sig2 = np.std(data2)
data1 = data1-data1.mean()
data2 = data2-data2.mean()
cross_correlation = np.correlate(data1, data2, mode='same'
                                ) /sig1/sig2/N

②FFT1

f1 = np.fft.fft(data1)
f2 = np.fft.fft(data2)
cross_spectrum = f1 * np.conj(f2)
cross_correlation = np.real(np.fft.ifft(cross_spectrum)
                            ) /sig1/sig2/N

③FFT2

f1 = np.fft.rfft(data1)
f2 = np.fft.rfft(data2)
cross_spectrum2 = f1 * np.conj(f2)
cross_correlation = np.fft.irfft(cross_spectrum2
                                ) /sig1/sig2/N

0をパディング + シフトで①＝②となります。

①直接法

padded_data1 = np.pad(data1, (0, N), 'constant')
padded_data2 = np.pad(data2, (0, N), 'constant')
cross_correlation = np.correlate(padded_data1,
　　　　　　　　　　　　　　　　 padded_data2,
　　　　　　　　　　　　　　　　 mode='same'
                                ) /sig1/sig2/N

②FFT1

f1 = np.fft.fft(padded_data1)
f2 = np.fft.fft(padded_data2)
cross_spectrum = f1 * np.conj(f2)
cross_correlation = np.real(np.fft.ifft(cross_spectrum))
cross_correlation = np.fft.fftshift(cross_correlation
                                   ) /sig1/sig2/N

Dispersion Analysis

SASW (Spectral analysis of surface waves)

extract the phase from cross-power spectrum S
Θ12(ω) = arg(S12(ω)) = k(ω)(x2-x1)
V(ω) = ω(x2-x1) / Θ12(ω)

MASW (Multi source Analysis of surface waves)
https://github.com/luan-th-nguyen/PyDispersion/blob/master/src/dispersion.py
calculate dispersion curves after Park et al. 1998

for fi in range(fmax_idx): # loop over frequency range
    for ci in range(len(c)): # loop over phase velocity range
        k = 2.0*np.pi*f[fi]/(c[ci])
        img[ci,fi] = np.abs(np.dot(dx * np.exp(1.0j*k*x),
                            U[:,fi]/np.abs(U[:,fi])))

Slant stack method after McMechan and Yedlin 1981

for fi in range(fmax_idx): # loop over frequency range
    for pi in range(len(p)): # loop over slowness
        k = 2.0*np.pi*f[fi]*p[pi]
        Upf[pi,fi] = np.dot(dx * np.exp(1.0j*k*x), 
                            U[:,fi]/np.abs(U[:,fi]))	
    Utaup = np.zeros((len(t)//2, len(p)), dtype = complex)
    for pi in range(len(p)):
        Utaup[:,pi] = get_ifft(Upf[pi,:])

通信方法

近年活発な現場での通信方法（LTE以外）です。

• スターリンク
  衛星を使ったインターネット接続サービスです。お値段は高めですが、施工現場では使い勝手が良いのでしょう。
• LPWA（Sigfox, LoRa, ELTRES など）
  これは実際に使うと驚きです。通信頻度にもよりますが、バッテリーで数年持つほど省電力なのに、無線で㎞オーダーの通信が可能です。富士山での利用例もあります。ビルの中でも５階差までは受信を確認しました（それ以上は未確認）。

今後期待される通信方法です。

• HAPS
  成層圏に携帯の基地局があるようなもの、だそうです。SoftBankさんが飛行機を飛ばしている紹介動画が、なかなかカッコ良い！

Survey123 レポート作成までのフロー

Survey 123 でデータを採取し、加工してレポートにするまでの大まかな流れです。

【Survey 123】

• 現地でデータ収集。

【ArcGIS Online】

• 収集されてできたフィーチャレイヤーのシンボルを変更。「詳細」「ビジュアライゼーション」タブにてレイヤーの「プロパティ」のスタイル編集からスタイルを変更します。

【ArcGIS Pro】

• 「カタログ」ウィンドウの[ポータル]タブからonlineのフィーチャレイヤーをマップに追加。属性にID等のフィールドを加えるなど加工してからフィーチャを保存すると、online のデータに反映されます。

【ArcGIS Online】

• 「詳細」「ビジュアライゼーション」タブで平面図に表示する「ラベル」を設定し、保存。
• MapViewer で地理院地図等を追加します。追加方法はコチラ。
  https://phreeqc.blogspot.com/2024/03/survey123-web.html

【Survey 123】

• 写真の入れ替え、追加の必要があればここで実施。先の MapViewer で実施すると Pro （Basic） でフィールドを変更できなくなったり、レポートで写真が追加されなかったり制限が多いのでダメ。
• 123 でレポートを作成。Pro で設定した連番を含めると平面図との対比が容易になります。
  調査結果の印刷—ArcGIS Survey123 | ドキュメント