SPH 基本1：補間式

汎用スカラー関数fのSPH補間式（連続・離散レベルの両方）　𝑓: ℝ3→ℝ

SPH法では、f(x)はx近傍の領域Ωにおける積分として記載される。
位置 x での物理量f(x)はディラックのデルタ関数を使うと、

\begin{align*}
f({x})=\iiint_{\mathrm{\Omega}}{f(\mathbf{x}\prime)\delta(\left|{x}-{x}\prime\right|)dV}\end{align*}

h→0 でW(x,h)→δ(x) となるような関数を使って近似すると、
\begin{align*}
\langle{f({x})}\rangle=\iiint_{\mathrm{\Omega}}{f(\mathbf{x}\prime)W(\left|{x}-{x}\prime\right|,h)dV}\end{align*}

これを離散化すると、
\begin{align*}
\langle{f({x_i})}\rangle=\sum_j{f(\mathbf{x_j})W(\left|{x_i}-{x_j}\right|,h)V_j}\end{align*}

これは、中心点xiの物理量を周辺の離散点xjの物理量に重みを乗じて合計した値とみなせる。

オンラインセミナー

コロナ禍により、セミナーが一斉にオンライン化してます。

今のところ、メリットしかありません。享受しています。

・場所を問わない。移動中でも参加できる。
・疑問に思ったら、即チャットで質問可能。
・施設も紙の配布資料もなくなるためか、無料のコースが多い。
・後で録画を見直せる。
・海外のセミナーも受講できる！

先週、イタリアのパルマ大学の方が SPH 関連のオンラインセミナーを開かれました。以前から受けてみたいなあと考えていたのですが、欧州は遠い。さすがに私費では難しい。しかも適当な英語しか話せない。ハードルが高かった。
そのような私がオンライン開催を見つけた時にはすぐに申し込んでいました。SPH 関連の専門用語なら聞き取れるだろうし、数式を追えばついていけるだろうと。イタリアなのでネイティブでないだろうから聞き取りやすいだろうと。まったくのポジティブシンキングです。

申し込んだ後、なぜかしばらく忘れていました。さらに案内メールも見逃していました。で、気づいたのが当日。１日目の録画のURLが送られてきて開催日初日だったことに気づきました。
次の日から出張と重なっていたのですが、道中に１日目の動画を確認。残る４日も、ホテルで受講。幸い時差の関係で開始時刻が夕方の16時～17時でした。最初の１時間くらいは受講できない日もありましたが、毎回１５分くらい遅れて始まるので、それもありがたい。後で動画を見ておくことで、十分にフォローできました（最後の方は数式すら半分以上わからなかったのですが）。
宿題も土日をかけて提出。録画が役に立ちました。
内容的に申し分なく、しかも無料で受講させていただいたことには感謝しかありません。

海外セミナー＋チャット＋録画は、私にとって最強の組み合わせでした。大きなハードルが、一気に解消されました。
欲張りすぎたのか、頭がパンパンな感じ。でも、今日も国内のセミナーに申し込んでしまいました。

オンラインセミナー、ありがたいの一言です。

****************************************
20201031追記

final examination に合格し、3ECTSをいただきました。Grazie!

equifinality problem

昨日の図書では、数値モデルに同位体をトレーサーとして融合する長所・短所が示されていました。
長所：モデルの信頼性向上
短所：校正すべき未知パラメータの増加、固有の不確定性（δ値の時空間的な不均質性など）

今まで、短所には目を向けず長所ばかりを見ていました。が、ある程度その見方も許容されるようです。一般には同位体データを含めて校正・検証されたモデルの方が信頼性が高いと言え、その融合に期待されているとありましたので。

この長所ですが、具体的には等結果性問題への対応ということになります。
等結果性（equifinality）：異なる構造・パラメータの組み合わせが同じ結果を示す
浸透流解析では、支配方程式で右辺の各項の値が異なっても、合計が左辺の値として一致することと同義です。

これはわかりやすい話です。
一般的な3次元浸透流解析では、地質構造（ジオメトリ）を１つ決定した後、水理パラメータの調整で水位や流量等の再現を試みます。これは、ジオメトリの調整・再作成に大きなコストを必要とするため、それを固定した後、相対的に安価なパラスタにて再現を試みているということです。つまり、構造の不確実性をパラメータに一手に担わせて調整していることになります（「地下水モデル」にも書かれていましたので、世界的な流儀でしょう）。
本来は文献調査や地質調査結果から複数の３次元地質構造の可能性を地質屋さんは推定します（通常は、その中の可能性の高いと考える構造を地質屋さんが１つ選択・提示しているにすぎません）。ところが、どのジオメトリでも浸透流解析では観測値を概ね再現できます。パラメータ側でその差を調整できるのです。つまり、同じ結果を再現できるモデル（ジオメトリとパラメータの組み合わせ）が複数誕生するのです。これが等結果性。
しかしながら、将来発生するイベント（トンネル掘削など）を入力した場合、どのモデルも同じ予測値を示すとは限りません。再現計算完了時点では、どれが実現象に最も近いモデルなのかわからないのです。
そこで、それらのモデルにδ18Oなどを乗せて流し、観測位置で得られたδ値を再現できるかチェックすれば、モデルの良否が判明します。時空間的に変化しない疑似トレーサーでも可能でしょう（時空間的に変化しやイオンを選択する場合には、Geochemical Modeling との融合が必要です）。これが等結果性への対応です。

「多角的に校正・検証したモデルは、より正確な予測値を出すことが期待できる」ということです。

環境同位体の利用方法

とある大学の書店に立ち寄りました。

理学部のないキャンパスなのか、地質系の図書はありませんでした。が、環境系のコーナーに以下の図書を見つけました（こういった配置、時々あります）。

山中勤「環境同位体による水循環トレーシング」

δ18O などの利用法は、文献からの知識や経験的に身に着けたように思い返されます。ですので、この図書の5章「地下水・湧水・温泉」の区分は系統的とまでは言えないものの、頭の整理に役立ちました。 

********************************************
1.涵養標高の推定：降水のδ値は標高が高くなるほど低くなる

2.流動系区分：δ値の高低で平面的な流動系区分

3.混合割合評価：
涵養源の寄与率評価：δダイアグラム上で河川水・降水・田面水などの混在割合を区分、GIS上でプロット
温泉水： δダイアグラム上で温泉水・天水（近隣の地下水や河川水）・非天水などの混在割合を区分

4.地下水流動モデルとの融合：
直接利用：δ2H, δ18O をトレーサーとして使用
間接利用：Particle Tracking から算出した河川水寄与率を、安定同位体やイオンの混合計算結果と比較

*******************************************

特に浸透流モデルのブラッシュアップにて、移流分散と粒子追跡の利用を同じくくりにできるか考えようとしたことがありませんでした（頭が固い）。
確かに、このような説明は国内で見ません。聞けば地下水のプロ？でも知らない方がいらっしゃるようです。土木分野で普及していないだけかもしれませんが。

他にも、経験的に知っていたことが明文化されており、あらためて知識として習得できたように思える個所がありました。新たな図書として、手元に置いておきましょう。

微動計3：波形のリアルタイム表示

微動計の設定に目途がつきました。

あとは測定しながら波形を確認したい。
これは簡単だろうと思って手を付けたのですが、ハマりました。

PC と微動計を無線で接続し UDP で吐き出されるデータを受信するだけなのですが、表示ソフトの画面に出てきません。ping は通っていますし、ポートもソフトがつかんでいるのですが画面に動きなし。ソフトが悪いのか？とPythonで試すも改善せず。有線ではどうか？これもダメ。
他の PC では？と、OS をクリーンインストールしたばかりの PC を用意しましたが、こちらもダメ。
念のため、もう1台の PC で試すと、
あっさりできました。原因不明。うーん。

その1台で波形を見ながら軽く衝撃を与えてみます。サンプリング周波数を変えながら、最適な測定条件をチェックします。やはり、確認しながらの作業は効率的です。

ある程度見通しが立ったので、次回、プレテストへ。
綺麗な波形が取れることを期待しています。

微動計2：オーバーサンプリング

微動計内では、これらのフィルターをいつ適用しているのか？

キーワードはアンチエイリアシング、オーバーサンプリング、デシメーション。
メーカーさんに伺うと、製品内の処理手順は以下の旭化成エレクトロニクスさんの解説と同じ手順でした。
https://www.akm.com/jp/ja/technology/technical-tutorial/basic-knowledge-adc/oversampling-advantage-1/
https://www.akm.com/jp/ja/technology/technical-tutorial/basic-knowledge-adc/oversampling-advantage-2/

オーバーサンプリングしておけば、測定対象の周波数帯に含まれるエイリアシングの影響を小さくすることができる。
その後、測定対象よりも高い周波数帯を、再度、デジタルフィルタ（LPF）でカット・デシメーションすればOK。
https://www.toyo.co.jp/mecha/casestudy/detail/id=12894

実際、複数の微動計でサンプリング周波数を変え、同じ衝撃をはかると異なる波形が得られます。低周波数に抑えて不思議な結果に悩むなら、最初から高周波にして物事を単純化しておくべきなのでしょうね。特に私のような入門者にとっては。

どちらかというと、安価で安定したセンサーを作る目的で発達した技術のようです。
これ以上は踏み込みませんが、奥が深そうな世界です。

微動計1： 信号処理

微動計を用いた衝撃の測定に呼ばれました。

今まで他の方が頑張っていました。が、得られたデータがイマイチだったということが判明。で、地盤を対象に微動計を使っていた私も御一緒することに。

どうも、微動計に備わっていたフィルタの周波数特性をうまく活用できていなかったのが一因のようでした。
その選定から見直そうとしたのですが、私も応用が利くほど詳しくありません。調べていくうちにわからない用語が出てきました。標準モード、群遅延特性、デシメーション。
使うツールの詳細を理解していないと、正しい調査ができないのは共通。いままでプロの指導でうまくいっていただけのことであり、場を変えると途端にダメになるのは同じ。
ちょうど週末で時間があり、この機会に理解することにしました。

用語自体は信号処理分野で広く使われているようです。いくつかの図書を見ましたが、基礎中の基礎といった感じでした。
わかりやすかったのが図書よりも以下のサイト。確かに、音響ととらえると理解しやすい（ホールの設計者はこのようなことまで考えていらっしゃるのですね！）。
以下、解釈した内容を( ..)φメモメモ。


**************************************
ヤマハサウンドシステム
短期集中連載　超解説FIR！第1回～第3回
https://www.yamaha-ss.co.jp/published-articles/journals-04.html

1．用語
1－1．周波数特性・位相特性・群遅延特性
入力したインパルス信号の応答を周波数領域で見る（フーリエ変換）と各周波数のエネルギー応答がわかる。これが周波数特性。
各周波数が時間的にどのくらい遅れて到達しているかを見たものが位相特性。
位相を周波数で微分したものが群遅延特性。
周波数とその位相の変化量（遅れ）が比例していれば群遅延は一定（直線位相特性）。 　
ex.
1Hz→1波長（2π、1秒）
2Hz→2波長（4π、1秒）
3Hz→3波長（6π、1秒）

1－2．インパルス応答、伝達関数
システム（ブラックボックス）にインパルスを入力した時に出力された信号がインパルス応答。
インパルスは全周波数で同じエネルギー量持っているため、インパルス応答を見ただけでブラックボックスの特性=伝達関数が分かる。
インパルス応答を周波数の観点で見ると周波数特性が分かり、周波数ごとの到達時間を調べれば位相特性が分かる。
インパルス応答、伝達関数、周波数特性、位相特性等は、データの見方を変えただけ。

2．特性を作る（リバーブエフェクトなど）
2－1．FIR フィルタ
FIR の Tap にある乗算器の係数はフィルターのインパルス応答そのもの。
係数操作でインパルス応答を人工的に加工することが可能。IIR のように最小位相として位相が変化するインパルス応答などを作ることも可能。
目的とする周波数カーブを描いてそれを逆フーリエ変換すれば、カーブのインパルス応答が得られる。
そのインパルス応答を乗算器の係数に入れることで目的の周波数カーブのフィルターを作ることができる。
FIR フィルタは入力信号（サンプルデータ）とインパルス応答波形の２つを畳み込んで新たなデータを作り出す。つまり「畳み込みをするための演算器」。

2－2．FFT と窓関数
入力信号も IR も共に FFT して周波数領域のデータに変換。
畳み込み積分は周波数領域に変換すると単純な掛け算として非常に高速の計算が可能。
その結果を逆フーリエ変換（IFFT)して元の時間領域のデータに戻す。
ただし、FFT を行なうため処理するデータを切り出す際にサイドローブ（雑音）が発生。
窓関数でサイドローブの軽減処理をしているが、軽減しようとすると反対に周波数分解能が悪くなるという二律背反の関係がある。
周波数分解能：　　　〇　　　　 △　　　　 　×　
サイドローブ軽減：　×　　　　　△　　　　　〇
窓関数：　　　　ハミング　　ハニング　　ブラックマン

**************************************

インパルスにδ関数が使われているのは当然なのに新鮮。つながりますね。 

 次は、これらの処理順序について。

文献 その2

ためがちな文献。

最近は 土石流や SPH 関連をよく読んでいたので、それ以外は放置でした。

以下、気になったモノと箇所を書きとめておきます。

**************************************

K-NET・KiK-netのPS検層記録に基づくVs・VPおよび深さの関係

物理探査/ 72 巻 (2019)

Vs-1000m/s 以下の領域では、Vp=500m/s or Vp=1500m/sに大別される。

土質によらない。不飽和、飽和による。

地下水の既知、未知に応じた回帰式を提案。

ポアソン比：K-NET:0.45以上（深度20mまで）、KiK-net:0.33（ばらつき大）

レイリー波の分散曲線の近似計算法の提案と地下構造推定への応用

土木学会論文集/1997 巻 (1997) 577 号

多層水平地盤構造における基本モードレイリ-波の位相速度の分散曲線をS波速度と層厚のみのパラメータで近似。

深層学習による崩壊・非崩壊地の自動判読手法の開発

日本地すべり学会誌/56 巻 (2019) 5 号

SonyNNC使用。深層崩壊地のタイル画像作成・利用

地すべり運動を磁気特性から探る－リングせん断試験の供試体を用いた帯磁率および残留磁化実験による検討－

日本地すべり学会誌/56 巻 (2019) 5 号

リングせん断試験により供試体の磁性粒子が円周方向に配列。磁気特性からすべり運動の方向を再現できる可能性がある。

天然ダムの決壊による洪水流下の予測と対策

砂防学会誌/45 巻 (1992-1993) 1 号

天然ダムの決壊時刻を推定する簡易推定図

単位幅ピーク流量と単位幅貯水面積との関係による単位幅ピーク流量の簡易予測図

(1)天然ダム決壊時のピーク流量の支配的な指標はダム高さと河床勾配である。

(2)天然ダムが高いほど単位幅ピーク流量が大。河床勾配が大きいほど単位幅ピーク流量小。

 

ネット上で引っかかった文献

塩化ナトリウムの地下への浸透による水質変化に関する研究

http://jgschugoku.jp/asset/00032/GE/Vol37/GE_vol37_gijyutsukenkyu01.pdf

凍結防止剤の散布

高速道路からの距離と塩化物イオンの関係：150m付近まで基準超過

安山岩と花崗岩領域の水質変化の違い

 

昨今の災害事例から見た道路斜面防災の在り方と留意点

http://jgschugoku.jp/asset/00032/GE/Vol37/GE_vol37_tenbo.pdf

治水3法から土砂災害防止法までの変遷・・・基本的にその対象が人家等の保全対象を含むもの。道路には対象から外れており、自衛手段が求めらる。

切土のり面の崩壊

①供用後5年以内

②過去に変状を来したのり面

③急速施工

④発生しやすい地質

⑤80％以上が深度5m、平均N値23未満。

⑥背後にリニアメント（断層、素因）

盛土のり面の崩壊

①盛土を通るリニアメント

②ボトルネック地形

③道路を挟む上流側にレベルバンクや集水地形

④崩壊の半数は供用後5年以内。20年以降でやや増える。

 

geochemical modeling と地下水の流路

Isotope hydrology and hydrogeochemical modeling of Troodos Fractured Aquifer, Cyprus: the development of hydrogeological descriptions of observed water types
https://www.sciencedirect.com/science/article/pii/S0883292720302729

安定同位体と geochemical modeling を組み合わせ、 地下水の涵養源と流路（水質の進化経路）について推察された報告です。

geochemical modeling は PHREEQC を利用。地下水の起源、流路、岩石との反応プロセスを示されています。他の複雑な帯水層システムにも適用可能とされていますが、ま、妥当でしょう。
浸透流計算を外しているため、PHAST よりお手軽です。この程度なら国内でも稀に受け入れてもらえるかもしれません。

残念ながら、国内にはこのような手法を扱う研究者がほとんどいらっしゃいません。研究者のお墨付きがなかったり、基準を作成してもらえないと、なかなかお客様も首を縦に振っていただけません。
そのためか技術者側もフレームワーク自体を理解しようとせず、「水質を調べて区分しました」「地下水年代を調べて区分しました」「〇％づつ混じっている可能性があります」程度で終わっているケースをよく見ます（こちらも残念です）。ま、利益に直結するルートが開拓されていないにもかかわらず、それに時間と頭を使うのはビジネスにとって省くべき無駄というのも理解できます。海外では図書レベルなのですが。

PHREEQC にハマってから10年以上が過ぎました。次の10年も同じようなものなのかもしれません。どこかに生かせる職業はないでしょうか？

SPlisHSPlasH

SPlisHSPlasH
https://github.com/InteractiveComputerGraphics/SPlisHSPlasH

先日の PySPH よりも式が豊富です。
しかも GUI で入れ替え可能。素晴らしい。

手元のソースの状態方程式や密度の時間微分の計算式を何度か入れ替えて試算していたのですが、これがあればずいぶん楽になるのでは？

などと考えていましたが、甘かった。
example はきちんと動くものの、自分でモデルを作る方法がわかりません。複雑な粒子配置はどのようにすれば実現できるのでしょうか？うーん。

「物理ベースシミュ」と説明されていますし、「泡」への対応も済んでいるようです。Maya や Blender のプラグインも用意されています。CG分野ですね。
この分野では、やはりリアルタイム計算を目指すのでしょう。速度面への対応は以下の通り。確かに、速い。
・neighborhood search on CPU or GPU
・supports vectorization using AVX

うーん。このまま寝かせるには勿体ない。
どこかに取っ掛かりがないでしょうか？

PySPH

PySPH を見てみました。
https://github.com/pypr/pysph

使用されている式は以下の通り。境界条件がイマイチはっきりしないのですが、それ以外は一通り(固体まで)そろっています。複数の式が出典とともに並べられていますので、文献を見比べたり探さなくて済みます。これだけでも参考になりますね。
https://pysph.readthedocs.io/en/latest/reference/equations.html

サンプルを動かすだけなら問題ないのですが、新たなジオメトリを作る方法がわかりません。以下のチュートを動かして少し理解できた程度でした。
https://github.com/pypr/pysph/tree/master/docs/tutorial

ブロック感覚で式を組み合わせる方式なので手軽です。生かしたいですね。