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 を行なうため処理するデータを切り出す際にサイドローブ（雑音）が発生。
窓関数でサイドローブの軽減処理をしているが、軽減しようとすると反対に周波数分解能が悪くなるという二律背反の関係がある。
周波数分解能：　　　〇　　　　 △　　　　 　×　
サイドローブ軽減：　×　　　　　△　　　　　〇
窓関数：　　　　ハミング　　ハニング　　ブラックマン

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

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

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