地下水モデル その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ベクトルのなす角」です（と言いながら、頭に入っていなかったのですが）。

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