文献を読んでいると、「Hilbert-Huang Transform（ヒルベルトファン変換）」が出てきました。
初めて知りました。何だろうと調べてみると、少ないながら日本語の解説がありました。
大塚ほか「Hilbert-Huang Transformによる非定常・非線形信号解析」寒地土木研究所月報, 2011
FFT では時間分解能を上げると周波数分解能が下がり、周波数分解能を上げると時間分解能が落ちます。両方上げたい場合、衝撃など短時間の波形解析の場合に利点でしょうか？
どこかの分野ではよく使われているのでしょうか？
Python にも HHT ライブラリはあるので対応可能でしょう。が、少し使いづらい？異なるデータセットに対して同じ周波数の値を得られた方が、あとで使いやすいかな。
ま、頭の片隅に残しておきましょう。
地震波干渉法について。
ノイズやコーダ波から相互相関関数を計算し（といっても、周波数領域では乗算）、グリーン関数を抽出する方法として知られているそうです。先日は cross-validation of ambient-noise として、計算していました。
その際に集めていた国内の文献です。時間ができたので、あらためて整理。
山下ほか（2010）地震波干渉法による西日本の地殻速度構造（1）
・F-net
・1ヶ月のデータ
・R, T, U ３成分
・記録長1時間、最大ラグ100秒で相互相関関数を計算後スタック
林田ほか（2014）中京堆積盆地における表面波群速度の推定-Hi-net連続地震観測記録を用いた地震干渉法に基づく検証-
・F-net と Hi-net の地震干渉法結果比較
・1年間のデータ
・毎正時から1時間毎、100秒で相互相関関数を計算後スタック
・前処理：0.05～2.0Hz、Running Absolute Mean Normalization に基づく正規化
※相互相関関数に対しMFTで群速度を算出。
周波数の刻み幅0.005Hz、帯域幅=対象周波数の0.8～1.2倍の範囲
中原（2016）解説 地震干渉法 その2 応用
・1年間程度が理想
・適当な時間長の抽出
・相互相関関数を求め、スタック
・前処理：バンドパスフィルタ、1ビット化、移動平均RMS振幅エンベロープによる割り算、周波数領域で白色化
起田ほか（2019）福島県の広帯域リニアアレイで観測された常時微動の地震波干渉解析
・F-net と Hi-net
・6ヶ月間のデータ
・3成分
・15分間の相互相関関数を4個平均して1時間重合相互相関関数とし、スタック
・前処理：1秒間の時間規格化、1ビット化、白色化
※相互相関関数に対しMFTで群速度を算出
※Takagi et al.(2014)の方法で群速度から表面波の分離
Vertical-Radial XCORR （ZR）、Radial-Vertical XCORR (RZ)の和からレイリー波の寄与を相殺し、P波の寄与のみを保持。差をとることでレイリー波の寄与のみを保持。
先日実施したのは、
・1ヶ月のデータ
・1時間（50％のオーバーラップ）で相互相関関数を計算後スタック（周波数領域）
・前処理：1ビット化、白色化
１ビット化が効いたような気がします。白色化（正規化）はまだ理解できていません。
この手法、デジタル信号を扱う音響分野でも利用されているようです。いろいろな分野で使われている手法の一つなのでしょう。お手軽ですから。
次は時系列変化、実体波と表面波の分離を目標にしましょう。
*****************************************************
20210626追記
実装内容が書かれた文献に、白色化に関する記載もありました。 正規化です。
Emanuel D. Kästle et al.(2016)Two-receiver measurements of phase velocity cross-validation of ambient-noise and earthquake-based observations
Cwhitened(ω) = (U1(ω)U*2(ω))/(|U1(ω)||U2(ω)|)
There are also other approaches to achieve this (e.g. Bensen et al. 2007; Ekstrom et al. 2009); Ekstrom (2014), however, notes that effects of whitening or removal of earthquake signals have only very small influence on the phase-velocity measurement.
ノイズが表面波主体とするなら、当然ですが、目的に応じた成分での使用を。
This procedure is applied to the Z- or R-component recordings for Rayleigh waves and T-component recordings for Love waves.
SimPEG の解き方です。
Thibaut Astic et al.(2019) A framework for petrophysically and
geologically guided geophysical inversion using a dynamic Gaussian
mixture model prior
https://academic.oup.com/gji/article-abstract/219/3/1989/5556947
PGI: petrophysically and
geologically guided geophysical inversion
GMM: Gaussian
mixture model
1 Initialization: m, Θ, z
2 while Φd > Φd* and Φs > Φs* do
3 (inversion)
objective function Descent Step (Gauss-Newton)
minimaize Φ(m) = Φd(m)+β{αsΦs(m) + ΣαvΦv(m)}
misfit Φd(m) eq.21 on d
smallness Φs(m) eq.4=20 on mref
Φs(m) = (1/2)||Ws(m − mref)||22
reference model mref, local weight matrix Ws
smoothness Φv(m) along{x, y, z}
→ update m(t)
4 (Maximum A Posteriori Expectation-Maximization eq.24 with 27, 32, 34)
fit a new GMM petrophysical distribution on m(t)
GMM global variables Θ = {πj,µj,Σj}
proportion πj = nj/n
mean μj (vector of size q)
covariance Σj (matrix of size q × q)
geological unit j (j = 1..c)
GMM simply sums the Gaussian probability distribution N weighted by their proportion:
P(si|Θ) = ∑j πjN(si|µj,Σj)
generalized version:
M(m|Θ) = Πi∑j P(zi=j)N(mi|µj,wi(-2)Σj)
→ update Θ(t)
5 Classif z(t) eq.21 on m(t)
z = argmax P(m|z)P(z) （=argmax πjN(si|µj,Σj)）
→update z(t)
→update mref eq.22 and Ws eq.23 on z(t), Θ(t)
6 update αs and β
7 Loop
流れはわかりますが、部分的に理解できないところが残りました。
これで収束するのが不思議。ベイジアンの頭の中を覗いてみたいですね。
共分散を使うので、物理量の関係を定式化不要なのが利点でしょうか。地質のようなカテゴリカルな結果も同時に得られますが、その分布確率が出るなら実務でも使い道が広がりそうです。
身に着けるには、読み込み＋手を動かさないとダメかな。
物理探査へのスパースモデリングの適用を調べていた際に、Joint Inversion というキーワードが目に留まりました。
いくつかの解き方があるようです。
複数の物理量を個別に最適化すると、それらの結果から統一された解釈を生み難くなります。本来、ある関係式で物理量が結びつけられているはずなので、それを定式化し inversion 時に考慮しよう、という流れが最も理解しやすいと思います。
pyGIMLi では、岩石の空隙率と飽和度を介し、弾性波速度と比抵抗値を結合しています。
Carsten Rücker（2017）pyGIMLi: An open-source library for modelling and inversion in geophysics
The Archie equation relates the bulk electrical resistivity of a medium ρ to the fluid resistivity ρf and saturation S depending on porosity ϕ (Archie, 1942):petrophysical は貯留性の視点からみた（空隙に目を向けた）岩の物理量でしょうか。Archie equation はリザーバーとしての砂岩に対し報告されていた式のようです。結晶質岩ではありません。透水性との関係も図に表現されています。
ρ = a・ρf・ϕ^(-m)・S^(-n) (11)
To create geophysical parameter distributions, we apply common empirical petrophysical models, e.g., Archie's equation (eq. (11)) that provides the resistivity ρ as a function of saturation S. Sonic velocity v (or its reverse, the slowness s) as a function of porosity ϕ and saturation S is given by the time-average equation (Wyllie et al., 1956):
s = 1/v = (1-ϕ)/vm + ϕS/vw + ϕ(1-S)/va (12)
電気探査（比抵抗）の3次元処理（地形の考慮）をクリアしておこうと、SimPEG を1日触っていました。
といっても、良い例題があるので、入力ファイルとコードを読めばすんなり理解できます（以前に動作確認済みですので、細部を読むだけです）。
ついでに sparse 正則化を導入 。norms が効かなかったので、2.5Dの例題を移植。結果はあまり変わりませんでした。ま、これで電気探査の 3D + sparse はクリアです。
あとは可視化。これに時間を取られました。
HPの絵のように、地形＋電極位置＋断面のような絵を作りたかったのですが、SimPEG からメッシュ＋計算値を吐き出す方法が見つかりません。あーだ、こーだやっていても時間ばかりが過ぎていきましたので、slack の コミュニティー に投げました。
気づけば、30分強で返答がありました。
「この図はPyVistaを使った」でした。 活発ですし、制作側から返答をいただけるのはありがたい。
SimPEG では、メッシュ分割に discretize を利用していますが、そのメソッドの一部、VTK関連を有効にするには PyVista を入れないとダメ、という仕様でした。いや、わからないでしょう。どこかに書いてあるのでしょうけど。
PyVista をインストールして、例題のコードを参考にすれば、できました。断面の 「scalar＝表示したいarray 名」を見つけるのにやや時間がかかりましたが、比較的あっさり図化できました。何度も結果を確認しながら、ということになれば vtu 書き出し ＋ ParaVeiw より役立つかもしれません。ま、その書き出しも PyVista のインストールが必須だったわけですが。
PyVista で地味にうれしいのが、XYZ データから直接 TIN サーフェスを作成できる点です。matplot などではサーフェスを作成する前に Z データをグリッド化しないと描けないのですが、PyVista では座標データそのままでOK。 これは、使えます。
電気探査の3次元地形補正。長年放置してきましたが、ようやく正攻法でクリアできました（開発チームに感謝！）。
Liam Toney, Kate E. Allstadt(2021)lsforce: A Python Based Single‐Force Seismic Inversion Framework for Massive Landslides
インストールは bashから。自動で新たな仮想環境が作られました。
CPSはインストール済みです。タイミングが良かった。有名なのでしょう。
https://phreeqc.blogspot.com/2021/05/2.html
https://phreeqc.blogspot.com/2021/05/cps-330.html
付属の examples (notebook)は動きました。
軽い。グリーン関数が地域の状況をまとめて表しているので、FEM などを利用した計算のように大きなマトリックスを解くような手法ではないのでしょう。これはありがたい。
データは IRISWS を利用しています。米国ですから当然でしょう。
今回は、以前落とした F-net のデータを利用したいので、IRIS を使いません。その場合、Green 関数を作成するのに CPS3.30 で作成したモデルを指定します（モデル作成のスクリプトはbinにありました）。
また、F-net データの ch＝測定方向を認識しませんので、gsacで事前に書き換えます（NE-TR変換に必要）。計算途中に Python で書き換えてもOKです。
で計算。
はい、SACのヘッダーを読めません。観測局の位置を認識しませんので、位置・震央距離を stream のヘッダーに書き込みます。
******************* ココから飛ばしてOK *******************
そのまま計算すると single force が大きすぎるようでした（移動体の質量を固定した場合、距離が１桁長くなります）。
win2sac で変換した sac はnm/s
https://seisman.github.io/HinetPy/api/HinetPy.win32.html
win2sac removes sensitivity from waveform, then multiply by 1.0e9. Thus the extracted SAC files are velocity in nm/s, or acceleration in nm/s/s.obspy は SAC のヘッダーを認識するため、stream に取り込んでも nm/s で正しく認識、図化されます。一方、IRISDMC から obspy で読み込んだ SEED はm/s。lsforce は後者を念頭に置いているため、代入する際にどちらの sream も m/s とみなしているのでしょうか？
Extract PZだそうで、SEED で落としてから obspy で読み込み、remove_responseする方針へ変更。
This function works for Hi-net network only.
F-net data users are highly recommended to use FnetPy to request waveform data in SEED format and extract instrumental responses in RESP or PZ format from SEED files.
Centroid Single Force Inversion を扱うコードはいくつかあるようです。
USGSから4月末に発表された文献です。コードは1月に公開されています。
Liam Toney, Kate E. Allstadt(2021)lsforce: A Python Based Single‐Force Seismic Inversion Framework for Massive Landslides
https://pubs.geoscienceworld.org/ssa/srl/article-abstract/doi/10.1785/0220210004/596477/
グリーン関数は IRIS から DL したモデルを使うか、CPS などで用意するか選択できるようです。
Green’s functions for a user-selected 1D Earth model are obtained automatically from the Incorporated Research Institutions for Seismology Synthetics Engine webservice or can be computed for custom 1D Earth models using Computer Programs in Seismology.正則化は L0,L1,L2 をブレンドきます。先週の joint inversion のセミナーでは P, S, 表面波を扱われていましたので、重なった領域なのでしょう。如何に自分が井の中の蛙状態であったかが痛感されます。
Regularization options include a blended zeroth-, first-, and second-order semiautomated Tikhonov regularization scheme
基礎理論の理解が浅いため自分では実装まで至りませんが、公開されているコードは使えそうです。ひとまず、手を動かしてみましょう。
土砂移動と振動にかかわる文献を読んでいますと、Centroid Single Force Inversion という手法がよく出てきました。
周辺で観測された波形から崩壊の規模や方向を Inversion で求めようといった内容です。日本人が火山分野で報告された文献がよく参照されていますので、もともとは Made in Japan の技術でしょうか。そういえば、2011年の奈良県の深層崩壊に対して、当時、京大の方が報告されていました。が、それ以降まったく聞きませんので、海外の方がここ10年ほどで昇華されたのでしょう。
ベースは地震分野で震源を求めるのに開発された手法のようです。CMT 解というキーワードはよく聞きますが、内容は知りません。良い機会なので調べてみました。キーワードは以下の通り。概要は理解できますが、手を動かせるまでには理解できていません。
**********************************************************
共立出版「地震学」p114-116
地震波はすべり量と等価な力系（単力源（シングルフォース）とモーメント）で表現される。
モーメント：力の強さ×軸と作用線の距離
作用線：力の作用点から力の向きに引かれた直線
モーメントテンソル：偶力のモーメント9成分
偶力(couple)：強さが等しく向きが反対で作用線が異なる2つの単力源 (single force)
単偶力 (single couple))：1組の偶力・・・回転する
双偶力 (double couple))：2組の偶力・・・回転を打ち消す方向
丸善出版「地震 どのように起きるのか」p37-38
グリーン関数:Gin(xi,t;ξi,ζ)
連続体の運動方程式の非同次項ρfiに強制振動としてのデルタ関数（インパルス=震源）を与えた場合の解
=ある方向のインパルスが時刻ζに瞬間的に位置（ξ1,ξ2,ξ3）に働く場合に現れる、同方向の地震動（変位）
δ(xi-ξi)：xi=ξiの位置にインパルスが発生。iは地震動(変位)の方向。
δ(ｔ-ζ)：ｔ=ζの時間にインパルスが発生。
ρfi=δinδ(xi-ξi)δ(ｔ-ζ)：体積力n=iの方向にインパルスが発生
共立出版「地震学」p120-122
発震機構の推定
・P波の極性による方法
・セントロイドモーメントテンソル（CMT）インバージョン（波形を利用した方法）
震源、観測点の位置、構造が既知であれば、シングルフォースやモーメントテンソルごとに理論的に変位波形を計算することができる。
u(xi,t)=Σj,k∫Gij,k(xi,t;ξi,τ)Mjk(ξi,τ)dτ
ξiにある点震源に働くモーメントの各成分は共通の時間関数Ｔ(t)と振幅αjkで表現(Mjk(t)=αjkT(t)) されるとすると
u(xi,t)=Σj,k[∫Gij,k(xi,t;ξi,τ)T(t)dτ]αjk
となる（変位uは振幅αの線形結合）。正確なGが既知であれば、最小二乗法を使うことでモーメントテンソル６成分の振幅を求めることができる
時間関数を仮定しなくても、両辺をフーリエ変換すれば観測変位スペクトルとモーメントテンソル各成分のスペクトルを線形結合の形にできる。
得られたモーメントテンソルをλ１=-λ3、λ2=0となるように回転させることで断層面の向きがわかる。
これらの手法では、構造の影響を正確に評価できるグリーン関数を用いる必要がある。
震源領域のモーメントの重心 (centroid) と規模も同次に求めることになる。=CMT解
朝倉書店「地震と地盤振動」p1
u(xi,t)=Σj,k∫Gij,k(xi,t;ξi,τ)Mjk(ξi,τ)dτ
周波数領域では、畳み込み積分は２関数のスペクトルの積となる
U(ω)=G(ω)M(ω)
長谷川ほか「地震学」共立出版
2.3.2 世界の観測網 p31より
ISC (International Seismological Centre)
・全世界の130の地震観測網やデータセンターと協力
・地震データの蓄積や公開
・毎年、震源などが観測資料として刊行されている
USGS
・世界各地の観測所のデータをもとに、地震の震源決定、発震機構推定を行っている
IRIS (Incorporated Research Institutions for Seismology)
・100以上の米国大学の共同研究機関
・世界各地の地震データの取得、管理、配信
・IRISとUSGSが全世界の協力のもとに展開している150以上の観測点からなるGSNの広帯域地震計データは準リアルタイムで提供されている
ORFEUS (Observatories & Research Facilities for European Seismology)
・ヨーロッパ、地中海域に展開される広帯域地震観測網などのデータが EIDA から提供されている
IRIS （というか、IRISDMC）は文献やツールでよく見かけます。多くの国の時系列データが集められているようですので、利用者が多いのでしょう。地震分野はツールもデータも豊富です。
