地すべりと微動

地すべりと微動に関する文献を読んでいました。

集めた中で最も古い文献が1973年です。意外と昔から研究されてきたようですね。港湾基準に取り込まれてから流行りだしたのかな？と思っていましたが、知りませんでした。

（気になった）主な内容は、以下の通りです。
・地すべりブロック内では、3~4Hz付近にて増幅率の増す報告が複数あり。
・地すべりブロック外でも同様の固有振動数を持つが、増幅はやや小さい。
・2~5Hzの particle motion では、表層のクラック直交方向の揺れが大きく見える。
・地質毎に卓越周波数が異なる。
　シラス:2Hz、花崗岩:6～10Hz、四万十層群:15～20Hz。
・表面波探査のS波速度と1/4波長則にて、地すべり深度を推定。誤差15%程度。

残念ながら報告例が少なく、微動が地すべりに十分活用されてきたとは言えないようです。盛土では動的解析を入れて微動を再現している例もありましたが、地すべりではまだまだ。平地で測る方が楽ですし、興味を持たれる方も少ないのでしょう。

以下、読んだ論文の概要です。
******************************************
島
常時微動と地すべり地への応用
地すべり　1973 年 9 巻 4 号 p. 1-8
・N-S、E-W方向とも大体同じような傾向。卓越周期ははっきりしない。→下層と上層の速度比はあまり大きくないか、あるいは表層の粘性係数が大きいか、また両者が一緒に原因しているかどちらかである。
・常時微動の観測をふつう平地でおこなうと、100mぐらいはなれたところで大きく変化しない。ところが、地すべり地をふくめて山地になると測点間の関係がかなりくずれる。このことはボーリング結果をみても同じであるから、やはり地下構造の変化を反映しているものとみられる

和田ほか
地すべり地の Crack 群の雑微動に対する影響
地震 1973 年 26 巻 4 号 p. 316-325
・亀ノ瀬地すべり
・2~5Hzのピーク：地すべりCrackに直交。crackが振動特性に強く影響。
・低周波領域（0.5~0.7Hz）：地すべり地全域で一様な振動特性。深部構造を表している。
・10Hzのピーク：断続的な振幅の強弱。人工ノイズ
・地すべり地のような地盤構造の複雑な地域では、単一の卓越周期のみに注目し、 垂直方向のみの地盤構造と対比するのは危険。
・地盤構造のディメンジョンに相当する個々の周波数領域についてそれぞれ考察することが必要
・地域下では潜在している crack 群 の発見と地すべり運動との関連性等を提供する可能性あり。

島
長野県鬼無里村地すべり地における常時微動と表層の振動特性
地すべり　1975 年 12 巻 1 号 p. 34-42
・粘性係数が大きいほど（Q値が小さいほど）減衰大
・地すべり崩落地に近いところでは粘性係数ξ=5×106～107CGS程度(Q値1～2程度）
・直接地すべりの影響をうけないところではξ=5×105CGS程度(Q値20程度）
・地すべりの発生からある程度時間が経過すると、地盤の弾性的性質は元の安定した状態に復元するものと思われる。表層の粘性係数は小さく（Q値は大きく）なり、スペクトルからみて明瞭卓越周期がみとめられる。

泉谷ほか
奈良尾地すべり地における常時微動特性
地すべり　1978 年 15 巻 3 号 p. 17-22
・仮定1：常時微動は基盤からのSH波垂直入射に対する地盤の応答である.
・仮定2:常時微動の源は時間的、 空間的にランダムに分布しており、 基盤より入射する波のy、 z成分のスペクトルは相等しい.
・仮定3:地すべり地のクラック群による土塊の異方性に着目して常時微動記録を解析することが、地すべりの状態を知る一つの手段となり得る。
・地表層上部を異方性体で近似した地盤モデル
・クラック群の到達している平均的な深さと、 クラックの混み具合とを推定できる解析手法を見出した.

泉谷ほか
常時微動測定による地すべり地盤調査の一手法
土木学会論文報告集 1981 年 309 号 p. 159-162
・「奈良尾地すべり地における常時微動特性」の改善
・常時微動によって推定された地下構造と、ボーリング調査等が対応するかを調べるにとどまっている.

島
常時微動特性から推定される味大豆地すべりの発生機構
地すべり　1981 年 18 巻 1 号 p. 15-25_1
・これまで、地すべり付近では振幅の大きくなる振動数はあらわれないで、振動数とともに振幅が小さくなるか、白色雑音のように振動数に関係せず振幅がほぼ一定になるようなスペクトルを示すことが多い。
・地すべり地下方では3～5Hz付近の振動数でスペクトル振幅が卓越する傾向を示す。
・味大豆地すべりも同傾向。
・地すべり地のように地形や地下構造が複雑なところでは、水平・鉛直成分の組み合わせによる楕円軌跡は水平軸に対して傾斜していることが多い。
・レイリー型の表面波であるなら、傾斜角は層構造の傾斜とみることができる。
・振動数ごとに傾斜角が異なっているときは、高振動数においてはごく地表の構造の傾斜に対応し、低振動数の場合は深い構造の傾斜に対応する。
・常時微動のスペクトルから地下構造を推定する場合、成層構造の中に低速度層を仮定してみると都合がよい。
・振動方向の分布にあずかる振動数と地下構造とを対比させると、低速度層より上部にある比較的硬い層が地すべり運動に寄与していることが分った。

川邉ほか
地すべりに及ぼす地震動の影響
地すべり　1983 年 36 巻 2 号 p. 5-16
・静岡県由比地すべり地において地震および常時微動の観測を行った。
・深度別二組の地震計の記録からスペクトル比を求め、 それをもとに重複反射理論を使って地表層の密度・剛性率・Q値の諸定数を推定した。
・固有周期で顕著なものはA点0.28秒(3.6Hz、増幅率6～10倍）、B点0.10秒 (10Hz、増幅率は5倍程度）
・0.1秒付近の振幅の増幅には深さ10m程度までの表層が関与。 0.3秒付近の振幅の増幅には、より深い幅が関与しているようである。
・推定された定数と重複反射理論を用いて、 基盤から地表層にある地震波が入射した場合に地表層内に発生する加速度 ・勢断応力を計算し、 土質試験より得 られたせん断強度と比較した。
・A点における最大加速度は、地表で最大で約3倍に増幅。最大勢断応力は境界面で最大約2.5kg/cm2となった。
・B点の各深さにおける最大せん断応力と、土質試験によるそれぞれの深さでの一面せん断強度を比較すると、 入力地震波の最大加速度が477ga1以上で表層内に破壊面
ができる。
・求められた加速度・勢断応力を由比の4個所の地すべり斜面の安全解析に導入した結果、 安全率が1になるのは傾斜約25度の斜面で50～110gal、約15度では100～220ga1のときであった。

秀島ほか
地すべり斜面における常時微動観測と一考察
開発土木研究所月報第457号　1991
・地すべり地域内の各成分の増幅は、地すべり域外のものと比べていくぶん大きい
・上下動成分は、地すべり域内および域外とも水乎2成分と比ぺて増幅は小さい
・地すぺり域内外でも、スペクトルの卓越振動数に相違はほとんど認められない。洪積層（第二種地盤）の基盤の振動特性そのものが強く反映された結果と考えられる。
・ハンドパスフィルクーは、既応の解析例をもとに2-5Hz。N-S方向(すべり頂部のクラックに直交する方向）が卓越した振動。

川邉
斜面表層の振動特性と不安定化
日本地すべり学会誌　2005年 42 巻 2 号 112-114
・地質毎に卓越した周波数
　シラス:2Hz
　花崗岩:6～10Hz
　四万十層群:15～20Hz
・表層部の固有周期、地質や土質、層厚などによって決定される地盤の振動特性を反映。
・下層から入射する地震動の卓越周期との関係から、表層部での安定性をある程度説明することが可能。

森ほか
微小地震観測による地すべり土塊の三次元形状と地震応答特性の評価
士木学会論文集Al(構造・地震工学）、 Vol.68、 No. 4(地震工学論文集第31-b巻）、 1_395-1_406、2012
・当地すべり地の微動は振幅が低く、電気ノイズに埋もれてしまうことが多い．微動探査は有効に用いることができない．
・微小地震観測は、地すべり地における卓越振動数の評価に有効。
・地すべりブロックの外側の安定した地山を基準にしたときのブロック内側の移動土塊部分の水平動スペクトル比（H/H比）は地震応答特性の評価には有効である．
・表面波探査等によるS波速度と合わせて地すべり深度を推定可能。誤差15%程度（H=Vs/4f、f=4H/Vs）

芝崎ほか
複数深度での地震動観測結果に基づく地すべり土塊の固有周期
日本地すべり学会誌　2016 年 53 巻 6 号 p. 227-234
・譲原地すべりおよび由比地すべりに設置された３深度（想定すべり面よりも下位の基盤岩層、想定すべり面の直上部および地表付近）の地震動観測を行い、高速フーリエ変換 によりスペクトルの増幅率を求め、地すべり土塊の固有周期について検討を行った。
・すべり面直上部付近の最大加速度は基盤岩層の最大加速度の1.0～1.9倍の値を、地表付近のそれは１．１～3.6倍の値を示した。
・いずれの地すべりにおいても、周期が概ね1～2秒を超える加速度フーリエ振幅スペクトルは基盤岩層と想定すべり面直上部、地表付近でほぼ同じ値を示した。一方、概ね1～2秒より短い周期では、基盤岩層に比べて想定すべり面より直上部および地表付近で大きな値を示す特徴が見られた。
・譲原地すべりにおける地すべり土塊の水平方向の固有周期は0.22～0.25秒で上下方向のそれは0.10～0.13秒、由比地すべりでは水平方向が0.37秒で上下方向が0.20~0.22秒と推定された。
・固有周期と加速度フーリエ振幅スペクトルの最大値を示す周期は、譲原地すべりでは異なる値を示した。
・由比地すべりでは、2009年の駿河湾沖地震において、深度40mと深度1mの水平方向の加速度フーリエ振幅スペクトルが最大値を示す周期と固有周期は同じ値を示した。
・由比地すべりで地中変位が発生しなかった原因として加速度が地すべりの変位を生じさせるよりも小さかったことが考えられた。

相関性と分離性

sklearn.preprocessing に LabelEncoder があります。

これ、categorical な変数を数値に変換・分類してくれます。sklearn を使わずとも、その前に既に整理されている場合もあるでしょう。馴染みのある分類方法です。

LabelBinarizer や OneHotEncoder もあります。
当初は馴染みのない上、スパースなのにメモリを多く使うため、その利点がわかりませんでした。が、気づくと便利。データの前処理で target との相関分析を行う際、並びを考えなくて良くなります。e>ａ>c>bの順で相関性があった場合、これらで機械的に番号を振ると良い相関を得られません。

機械学習にかけるうえでは、相関性だけでなく、分離性についても着目すべきでしょう。random forest は相関性重視ですが、LightGBMは分離性も見ています。扱う手法が何を重視しているかによって、与えるべきデータ、前処理が異なります。難しい。

2値分類の場合、KDE 等でピークの分離性を見てみます。が、数百もあると目検索では辛い。最大ピークだけ抜き出して距離を測る方法をとっていますが、まだ納得できていません。第2、第3ピークも考慮したいのですが、まだ実装できていません。プロはどうなさっているのでしょう。
https://stackoverflow.com/questions/16255622/peak-of-the-kernel-density-estimation
http://ibis.t.u-tokyo.ac.jp/suzuki/lecture/2015/dataanalysis/L9.pdf

データを理解するためには、ツールと表現方法を実装する必要があります。
もっと多くの知見を得たいものです。

相互相関関数

先輩から「これ、よかったよ」と勧められた図書です。

ディジタル画像処理編集委員会「ディジタル画像処理 [改訂新版]」

全く興味がなかったのですが、お借りして良かったですね。
フィルタ処理や領域処理等、各種画像処理の考え方が包括的に掲載されています。画像を xy 2次元＋RGB 3次元の配列と考えると、線形代数等で扱われている処理がそのまま使えます。回転テンソル、固有値解析など。ベイズ統計やフーリエ変換も出てきました。
画像処理には何らかの数値処理を施しているとは思っていましたが、求める結果を得ることが先で、その中身まで深く考えたことはありませんでした。まさか領域分割の一つにカーネル密度関数とその極大位置が利用されているとは。興味を持たないというのは怖い。

特に興味を惹かれたのが、相互相関関数を利用した画像マッチング。ちょうど、「SPAC係数は自己相関でなく相互相関でない？」などと思案中でしたので、驚き。

残念ながらこの図書には実装できるまでの実用的な記載はありません（ソフトが処理してくれますから必要ないのです）。python だと 2D の相互相関関数もあるかな？と思って調べると、ありました。まさに画像マッチング。

scipy.signal.correlate2d
https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.correlate2d.html
コチラ↓はノーマライズしてます。地球統計学でも標準偏差で割っていた似た式がありました。
https://stackoverflow.com/questions/1819124/image-comparison-algorithm

correlate2d には mode が3種 （full, valid, same）あるようです。どのような結果になるのでしょうか？

このような結果でした↓

a=np.array([[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]])
b=np.array([[1,1,1],[1,1,1],[1,1,1]])
print(a)
print(b)

[[ 1  2  3  4]
 [ 5  6  7  8]
 [ 9 10 11 12]
 [13 14 15 16]]
[[1 1 1]
 [1 1 1]
 [1 1 1]]

cor2 = signal.correlate2d(a, b, mode='full')
print(cor2)

[[ 1  3  6  9  7  4]
 [ 6 14 24 30 22 12]
 [15 33 54 63 45 24]
 [27 57 90 99 69 36]
 [22 46 72 78 54 28]
 [13 27 42 45 31 16]]

cor2 = signal.correlate2d(a, b, mode='valid')
print(cor2)

[[54 63]
 [90 99]]

cor2 = signal.correlate2d(a, b, mode='same')
print(cor2)

[[14 24 30 22]
 [33 54 63 45]
 [57 90 99 69]
 [46 72 78 54]]

左上から順に重ね合わせる作業は3つとも同じ。返すarray の範囲が異なっているだけです。same が使えそうですね。
個人的には、最大値で該当箇所を見つけられる以上に、異なる次元のベクトルを比較して類似度を判定できる性質に、相互相関関数の価値を見出せそうです。

興味がないと割り切ってしまえば終わります。なんでも手を出してみるものですね。
先輩に感謝。

事前確率

ベイズ統計について初心者用の図書を購入しました。

松原望「やさしく知りたい先端科学シリーズ1 ベイズ統計学」

本当に優しく、2日で読み終わりました。で、感想。

「私はベイジアンだったのか」

物事が起こるたびに、期待値を変更します。それが経験則となります。その過程をベイズ更新として肯定されているので、馴染みやすかったですね。

以前、機械学習による胃がんの検出例についての疑問を書き残しています。

https://phreeqc.blogspot.com/2019/05/ai.html

これと同様の内容を、ベイズ統計の考え方で説明されている個所がありました。わかりやすかったですね。事前確率の重要性を尤度と同レベルで述べられています。

この場合、事前確率が非常に小さい、ということを見落としてしまうと、80%の精度の検査で陽性、ということに驚いてしまうかもしれません。メディアなどでもAの事前確率を隠してBの尤度（条件付き確率）だけを提示して印象を操作しようとすることがあります。これからの世の中では事前確率やベイズの定理を意識しながら、冷静に提示された数字や情報、データを読む力が肝要となってくるでしょう。

雨が降って、災害が発生するかしないか。尤度を90%程度に引き上げることは比較的容易です。問題なのは上記の例と同じくらい低い事前確率。尤度を９9％近くまで引き上げないと、「命の危険があります」と繰り返してもオオカミ少年状態。この数%が難しい。
医学分野ではどの程度でスクリーニングに「有用」と判断されているのか知りたいところです。

以前からデータを読む力は重要であったに違いありません。近年では、量・種類共に膨大なデータを加工する能力と共に、その力が必要に迫られているのだと考えます。

Kuzushiji Recognition

Kuzushiji Recognition
https://www.kaggle.com/c/kuzushiji-recognition

先日始まった kaggle のコンペです。日本の人文学オープンデータ共同利用センターなどがホストをされているようです。

このような人文学の研究団体を初めて知りました。いえ、当然あるのでしょうけど、縁がありませんでした。
HPを見てみましたが、興味の引かれるデータを多く公開されています。試してみようか？のみならず、少し読んでみようかと心動かされてしまいますね。
考えてみると、こういった読めなくなった文字を機械学習で認識させる、翻訳するというのは王道です。

土木分野でも、このようなコンペを開催すれば、物によっては比較的安く、良いアルゴリズムを得ることができると思います。プロポよりも効率的でしょう（経産省さんは昨年実施されていました）。

知らない分野でも、知らない間に機械学習の利用が進んでいます。
井の中で目標を立てることなく、外にも目を向けましょう。

地下水モデル その4

10章は不確実性。

私に欠けている知見であり意気込んで読み始めたのですが、正直よくわからないところがありました。（備忘録以上のまとめ量になったので掲載を控えます。）
理解するには手を動かさないとダメなようです。

MCMC、ベイズ統計は、機械学習でも出てきます。最近、本屋でもよく見かけます。データ同化は古くからあるようですが、まだ理解できていません。次はこのあたりの知識を仕入れないといけないでしょうか？

11章は報告書。
この章以外にも、裁判を念頭に置いた地下水モデリングについて触れられています。立ち位置含め、技術者が気を付けないといけないところでしょう。
シビルアクションを見てみましょうか。

以下、備忘録です。
************************************
11.2 モデリングの報告書
11.2.4.1 水文地質環境
・流向を模式的に矢印で示した図はモデルの周囲境界と領域内における一般的な流行を表すために用いることができる。
11.2.4.2 システムの特性
・パラメータの記述（パラメータを求めた手法、推定範囲、空間分布）
11.2.4.2 概念モデル
・画像（模式図・図表）により概念モデルを表現
11.2.5 キャリブレーション

11.4 査読
・査読チェックリスト
・対立する関係者との争点にならない部分では合意を求め、重要な争点を明らかにすることができれば、訴訟は合理化。モデリングはその枠組みを提供する。
・モデリング過程では、利害関係者と意思決定者との意思疎通が重要であるとの意思疎通が高まりつつある。

地下水モデル その3

9章はモデルキャリブレーションです。

4年前の時点で PEST の使用は当たり前、手動の試行錯誤のみはダメ、と断言されています。10章の不確実性にかかわる箇所でもあり、必須なのでしょう。そうなると、使用するソフト制限されます。

ペナルティ項の追加はタンクモデル（SCE-UA）の最適化でも使われていましたね。機械学習でも L1, L2正則化として用いられていました。目的は異なりますが、数式にすると同じです。このあたりは古くからある理論が実務で利用されているということでしょう。

V&Vは数年前に聞き始めましたが、ぼーっとしているうちに、古くなっていたようです。恐ろしい。

あと、この章だけではないですが、予算の配分についても若干の記載があります。日本産の教科書にはない点です。有益です。

全体として地下水にかかわる技術者ならば、心得ておくべき内容が多く記載されていると考えます。あたりですね。

地下水モデル その2

地表流と地下水の連成に興味を惹かれ、6章を読みました。主に MODFLOW で扱われている内容です。
以下、備忘録です。

************************************
6.2 揚水井および注入井
・3つのモデル化
　・ユーザーが層間での揚水・注水量を割り振る：単純すぎて実際の問題にはほとんど適用できない。
　・鉛直透水係数を高くする：許容範囲内の解が得られるが、計算は不安定。乱流・井戸損失は表現できない。
　・MODFLOW井戸パッケージMAW、MNW

BOX6.井戸節点周辺における接点間隔の設定指針
・有効井戸半径re=0.208a→a=4.81～6.66rw

6.3 湧き出し・吸い込み
・フラックス（L/T）で与えるか？堆積流動速度（L^3/T）で与えるか？コードによる。

6.4 排水および湧水
・トンネル排水等：水頭依存境界
・流出のコンダクタンス（L^2/T）をキャリブレーションで推定するのが通常。

6.5 河川
・計算流出量から河川水位を計算
　SFR1パッケージ：マニング式
　SFR2パッケージ：河道特性をセル毎に指定可
・地下水モデルに河川水位の計算を含むと、非線形方程式が加わることで解の安定性に影響。
・解析対象から隔たった地物に対し簡便な河川パッケージ等を適用し、近傍の河川でSFRを利用するというのが通例。

Box6.2
・地表流と不飽和流コードの連成：Hydrogeosphere→時間・空間ステップを共に細かくする必要あり。長期の計算では大きな計算資源が必要。
・計算時間の短縮→GSFLOW:地表水+1次元不飽和浸透+地下水モデル
Box6.3
・地表水：少数地点の河川ハイドロ→時間的に密、空間的に粗なキャリブ
・地下水：時間的に粗、空間的に密なキャリブ
・地下水からの湧出（飽和余剰地表流）の影響大。

地下水モデル その1

先月発売になった「地下水モデル」の気になる箇所を、約1か月かけてコツコツ読んでいました。

地下水モデル―実践的シミュレーションの基礎― 第2版
https://www.kyoritsu-pub.co.jp/bookdetail/9784320047365

和訳初版が1994年ですので、25年振りの改稿です。実際は2015年版の和訳ですので、既に4年前の情報になっています。が、学ぶべきことは多くありました。

第1版との大きな違いは、「キャリブレーション」と「不確実性評価」。この2点が新たに追加されています。後者は私にとってピンポイントで不足、必要に迫られていた内容です。また、前者も十分とは言えません。
実務的には PEST の使用が前提でしょう。PEST は9年前に知り、便利なツールだと感じていましたが、それっきり。時代が進み、何とか実施しないといけないレベルまできていることは、はっきり自覚できました。

以下、1章の備忘録です。
************************************
1.2 モデルとは
・物理モデル：タンク、カラム
・数理モデル
　・データ駆動型：ブラックボックス、ニューラルネット
　・プロセス型：確率論的モデル、決定論的モデル

1.3 モデルの目的
・不確実性が必ず伴うことを強調するため、「予測」よりも「予報」という用語を好んで用いる。
・気象「予報」は確率（例えば、降水確率）によって記述される。

1.4 モデルの限界
・地下水モデルは現実の単純化→近似に制約を受ける。
・実際の複雑な自然を表現仕切ることのできるモデルは存在しない。
・不確実性を評価し、報告しなければならない。

1.5 モデルのバイアス
・依頼主から金銭的授受を受けているモデル作成者が中立を保ちバイアスから逃げられるのか？
・中立性を保ち、プロとしての信頼性を維持することが重要。

1.7 よくあるモデリングの誤り
・プロジェクトの半分の時間と予算をキャリブレーションに当てることを推奨する。
・不確実性解析の実施→予想より長い時間がかかる。モデリングの初期段階に立ち返る必要あり。
・モデル報告に十分な時間を取る。読むに耐え得る包括的な報告書を作成→再検証に必要。

モデルビルダー ArcGIS10.6

1週間ほどArcGISに取り組んでいました。

ポリゴン同士の重なった面積を集計するだけだったのですが、いくつかの条件を満たす必要があり、意外と時間がかかりました（初心者ということが最大の原因ですが）。

その中で知ったのがモデルビルダーと Python。
モデルビルダーで組んだモデルを Python スクリプトとして export できました。QGIS でも Python をサポートしていましたが、Arcでも同様（どちらが先かは知りません）。
繰り返し処理の部分は手直しが必要です、が、形にはなっています。

疑問点はニーズ。どのようなニーズがあるのでしょうか？
モデルビルダーでできることをあらためて Python で組もうとは思わないですし、その他の必要性を感じるほどのヘビーユーザーでもありません。

ヘビーユーザーの使用例を知りたいですね。まだまだ「既知の未知」があるようです。早く必要に迫られるレベルまでたどり着きたいものです。

重力勾配テンソル その2

続きです。

空間周波数fの定義は「回/m」。
対象範囲から最大・最小波長λを出して、最小・最大周波数fに変換すれば良いだけ。先の STACK OVERFLOW に python コードの例が載っていましたが、もっと簡潔に書けます。
周波数がわからなくても、波長から角波数k=2π/λを出せますし、波数だけならフーリエ変換時に０から順に出せるので、計算すら必要ないかもしれません。

曖昧な点は、実装しながら確認することに。

まずは文献に沿ったデータを用意。
次に「新・地震動のスペクトル解析入門」のソースを確認。これでも十分に短いのですが、python では numpy で1行、np.fft.fftn だけでOK。ありがたい。2D だと fft →転置→ fft →転置と同じ結果になったので、この順で内部処理しているだけかもしれません。

実装自体は容易で、2晩ほどでできました。

が、結果が合いません。似たような重力勾配や水平微分の分布になるのですが、オーダーの異なる場合があるのと、ky の符号が逆になっているように見えます。
角波数でダメならただの波数を試したり、ナイキスト周波数以上の領域を全て考慮したり、しなかったり。組み合わせによっては、それっぽい結果になるのですが、数値を記載している文献がないため確証を得られません。3日ほど考えましたが、最終的にはあきらめて寝かすことに。詳細に書かれた文献が、出てくるかもしれません。
もっと簡単に結果を出せると思っていたのですが、残念。

理解にはまだ時間が必要なようです。
ま、重力探査の現状と自分のレベルが分かっただけでも、良しとしましょう。

******************************************
20190813追記
fftn の並びに対応する fftfreq で使う符号、波数のあたりが怪しいようです。修正すると、それらしくなりました。が、検証できないのでここまで。もう少し寝かせましょう。

重力勾配テンソル

講習で知った、重力異常値から重力勾配テンソルを導く方法について、調べていました。

理論はコチラ。
Kevin L. Mickus, Juan Homero Hinojosa
The complete gravity gradient tensor derived from the vertical component of gravity: a Fourier transform technique
https://www.sciencedirect.com/science/article/pii/S0926985101000313

重力異常を重力ポテンシャルの2階微分した重力勾配の形に直しているだけなのですが、フーリエ領域で波数を使った乗算の形にしており、水平成分も算出できるように工夫されています。これを9成分の重力勾配テンソルに整理し逆変換すれば、空間領域の重力勾配テンソルを求められるという流れ。おそらく重力偏差の測定・解釈から思いつかれたのでしょう。

シンプルな流れなのですが、理解に時間を要しました（符号の誤りもありましたし）。
特にわからなかったのが波数の考え方。今回は時間領域への変換ではなく、フーリエ（空間）領域？への変換。ここで躓きました。
2次元平面から波数を求めるイメージはできても、具体的な数値にする方法を複数思い付き、定められません。
調べてみると、同じように悩まれている方がいらっしゃいました。
https://stackoverflow.com/questions/7161417/how-to-calculate-wavenumber-domain-coordinates-from-a-2d-fft

これ、案外メジャーでした。画像処理分野で。
2次元平面での周波数を「空間周波数」と呼ぶそうです。そういえば、画像処理でもハイパスフィルターとかありましたね。画像を扱う際は見た目さえ良くなればOKなので、具体的な周波数には意識を向けていませんでした。が、内部ではまったく同じ計算をしています。うーん繋がる。

続く。

物理探査

先週、物理探査学会の講習会に参加しました。

若い方を対象にしているのかな？と思いながらの参加でしたが、意外とオジ（イ）サマも多く来られていました。
内容は基礎的なことから最新の内容まで。参加するまではもったいないかな？とも考えていましたが、講義についていけないところもあり、個人的には充実した内容の講習会でした。

全体としては、取得データ数の増大、3次元解析の充実が印象に残りました。データの大量取得、大量処理が可能となっている現在、いつまでも2次元にとどまる必要はないのでしょう。3次元解析、当たり前にできると言えるようになりたいものです。

また、処理能力の向上に関しても知らないことが多くありました。多重反射をある程度除去できるようになっているとか、重力探査でソースの走向傾斜を推定できるようになっているとか（これ、ほぼ理解できるようになるまで、3日ほど論文を読み返しました）。

年を取ると「知っている」と勘違いすることが多くなります。気を付けてはいるのですが、なかなか本当に「知らない」ことを自覚するに至りません。今回は良い機会でした。
調べて身につけましょう。

以下、個人的な備忘録です。要チェック！
*****************************************************

・レーダと電磁探査（ハンドブック図9.1）
　・マクスウェル方程式、タンデルタ

・空中電磁探査
　・逆解析で深度決定。
　・表皮深度を使う簡易法では、誤差が5倍程度。

・電気探査
　・4端子法により接地抵抗をオミット。

・SAR
　・運が良ければ2時間ほどで取得可。
　・品質証明が肝要。

・微動アレイ探査
・SPAC係数
　1.分子：クロススペクトルS12
　2.分母：パワスペクトルS11
　3.ρ12=real[S12/(√S11√S22)]
　4.ρ12aveとρ13ave、ρ14aveの算術平均を出す

・位相速度の求め方
　・スパック係数と位相速度を、第一種0次のベッセル関数J0を介して関係づける。
　・J0：アレーサイズから高周波数が決まる。低周波数側は急激に落ち始めるところ。

DETERMINATION OF SOIL SHEAR MODULE AT DEPTHS BY IN-SITU VIBRATORY TECHNIQUES, ARMY ENGINEER WATERWAYS EXPERIMENT STATION VICKSBURG MISS

・重力探査
・半自動解釈手法

• 固有値・固有ベクトルを用いた解析
  ・産総研データ
  　→フーリエ変換
  　→積分→微分→Gz、Gx,Gyがそろう
  　→半自動解釈で断層傾斜角

• 高密度体(基盤岩)の方向に重力偏差テンソルの最大固有ベクトルが向くことを利用し、断層 傾斜角を推定
• 日本でも、地熱地域を中心に、重力偏差探査(重力ポテンシャルの3次元空間微分)が実施されてきている [空中探査]
• インバージョンよりも短時間で重力異常や磁気異常の異常源を推定する半自動解釈手法は、広範囲の大雑把な構造を知る解析に向いている
• JOGMECのHPから申請。データ取得。

コサイン類似度

コサイン類似度が出てきました。

以前、自然言語処理の図書で出てきたように思います。が、無関係でしたので読み飛ばしていました。

今回、あらためて数式を見ると、分子がベクトルの内積。分母がノルム。ノルムで割って単位ベクトルにした後、内積を取っているイメージ。それが cosθ になるからコサイン類似度。方向が同じなら似ている（cosθ=1.）ということでした。

これ、高校か大学1年生レベルの数学ですよね。ただの「2ベクトルのなす角」です（と言いながら、頭に入っていなかったのですが）。

まだまだ基礎力不足。頑張りましょう。