判別分析

他支店の方から、SPSS を買ったと伺いました。

判別分析などを実施する目的で購入されたとのこと。
描画は遅いものの計算部分は高速で、数百万行のデータでも容易に扱えるそうです。さすがに商用ソフトは最近のハードに最適化されているのでしょう。逆に、数百万行を高速に扱えないソフトは売れないでしょうね。

聞いたところ、数量化2類が入っていないとのことでした。
調べてみると、ずいぶん前にオプション販売は終了しています。ま、判別分析が入っているのなら前処理としてダミーデータを作成（ integer encoding ）しておけば良いだけなので、オプションとしての需要がなくなったのでしょう。あるいは、GBM など性能の良い判別機が出ていますので、アルゴリズム自体をあえて利用する機会が減ったのかも。

最初、判別分析は SVM など Classification Algorithm の総称を指しているのかと思って聞いていたのですが、異なっていたようです。scikit-learn で言うところの LDA, QDA アルゴリズムを指していたようです。見る限り、分布のピークを分離しているようですね（利用を迫られた段階で調べます）。
https://scikit-learn.org/stable/modules/lda_qda.html

GUI で高速というのはありがたいでしょう。備えていなくても迫られた際にすぐに答えが出せる点で実務向き。
個人的には Python で十分なのですが、どのような機能があるのか知りたいところです。

RFEM & FORM

昨日と同じ文献です。
Man-Yu Wanga et al.(2020) Probabilistic stability analyses of multi-stage soil slopes by bivariate random fields and finite element methods

RFEM はガウス分布に従ったモンテカルロシミュレーションによって不均質性を表現し、FEM や FDM モデルとした後に SSR 法により安全率を求めるという手法（空間にかかわるパラメータの扱い方がまだ理解できていません）。
知らなかったのですが、比較的多くの事例がありますね。

RFEM に至る流れは、おそらく以下の通り。
3つのレベルで確率を考慮する手法が紹介されています。FORM は通過点でしょう。
https://ieaghg.org/docs/General_Docs/5th_Risk_Assess/VaughanGriffithsieaghg_secured.pdf

FORM において標準正規分布を盲目的に当てはめるのが乱暴で、全体の分布を把握しないといけないというのが直感的にわかる図、言葉も含まれています。
The "safer" slope has a higher "probability of failure"!

FORM は国総研資料第379号でもを扱われています。正規分布でない2次元の確率分布を正規分布で近似し、性能関数に直行する1軸で評価する手順を示されています。取り扱いを簡易にしたいということなのでしょう。ちなみに、扱われた例では second-order でも精度は変わらないという結果でした。
http://www.nilim.go.jp/lab/bcg/siryou/tnn/tnn0379.htm


10年以上前の資料まで追いかけてみましたが、知らないことが多くありました。確率論を避けていたからかもしれません。
ありがたいことに10年以上前では取り扱いが難しい不規則な2次元確率密度分布でも、今ではPython のおかげで誰でも容易に取り扱うことができるようになっています。2変量を標準化しなくても、そのままカーネル密度推定で表現できますし、ある領域化の確率密度を積分することが可能です（手を動かしたことはないですが）。
試してみましょうか。

安定計算と崩壊確率

災害や崩壊にかかわる危険度として相対値なり確率なりの導入を試みる事例を時折目にします。

斜面安定計算でも確率を扱えます。

安定計算に強くかかわるキーパラメータに ｃ・φ があります。この2軸上で確率密度分布を表現するには多くの試験が必要になります。
円弧すべりを考える場合に1軸＋簡易CUを多数実施し、ｃの1次元分布を把握する。それをカーネル密度分布などで連続関数にしてFs=1.0未満の累積確率を求める。このような手順を踏めば、ある条件下での崩壊確率は得られるでしょう。
港湾分野では他分野よりも対象とする土質が比較的均質で室内試験を多く実施する傾向にあるため、上記のような確率の考え方を導入することは容易です（平成19年版設計事例集でも FORM として触れられていました。最新版では消えましたが）。

これが斜面だと難しい。基本は水も土質も均質ではない条件下で、円弧にならないすべり面形状が前提のため、多数の観測や室内試験が必要となります。この点で難色を示される方が多くいらっしゃるでしょう。仮に、試験結果が多数得られたとしても、基準にない確率の考え方を受け入れてもらうことへの抵抗が予想されます。基準に載っており、簡単で、実績のある手法以外を導入し難い気持ちはある程度理解できます。

多量のN値があれば、その確率密度推定結果をφの分布に変換することは、場合によってはありかもしれません。が、数がないので標準正規分布といった仮定を持ち込むのは、乱暴でしょう（数値実験なら良いのですが）。ちなみに、数値実験では標準正規分布ではなく、対数正規分布が多く使われているようです（∵non negative）。例えば以下の文献。
Man-Yu Wanga et al.(2020) Probabilistic stability analyses of multi-stage soil slopes by bivariate random fields and finite element methods

仮に2次元の確率密度分布が得られたところで、それを扱うにも精度の問題が生じます。安全率では計画安全率の0.1%の不足は認められません。確率密度の積分で、その精度が必要とは思えないのですが、基準などで「確率〇％以上を要求」などと決まってしまえば、安全率同様に0.1％の誤差は許されなくなります。計算時間がかかっても、精度の良い積分方法を選択せざるを得ないでしょう。


確率を論じるにはデータが必要になり、それには費用と時間を要します（リスク発現時の対応費用を考えると、微々たる金額ですが）。
こうなると純粋な技術論ではなくマネージメントの領域になります。その判断を経験に依存する分野では、確率論との併用化を進めることは難しいでしょう。

確率を導入する流れに反対ではありません。が、実務上、クリアすべき問題点は多くあるように感じます。

出張

今年度の納品は対面であったりテレビ会議であったり。

テレビ会議は御客様側で戸惑いがあるかな？などと思っていましたが、周りを見ていると案外スムーズに進行しているようです。営業さんに話を聞くと、お客様側では次年度から実施予定でテスト中だったとのこと。それを前倒しされたようです。次年度からは増えるでしょうね。

対面を望まれた場合は伺うのですが、遠方の方は長時間の移動を気にされています。新幹線は空いていますが（感覚的に1/2～1/3程度）、咳き込む方が長時間いらっしゃると私も気になります。

今回の新型コロナウイルス感染症の件では組織の有するリスクマネージメント能力、危機管理能力が露呈したように感じています。前者に成功したという優秀な例は少ないように感じます。個人も同様。私も準備していた使い捨てマスクが底をつきそうです。リスク管理は失敗ですね。

まだ出張はあります。お子様連れの長時間移動を避けられない方もいらっしゃるでしょう。通院が必要で出かけないといけない方もいらっしゃるでしょう。病で免疫力の落ちた方は特に心配です。
個人でできる対策は限られていますが、常に配慮する意識を共有できれば良いと思います。

前処理としての CNN

Zhice Fang (2020) Integration of convolutional neural network and conventional machine learning classifiers for landslide susceptibility mapping, Computers & Geosciences, Volume 13

この文献では特徴量抽出の1手法として CNN を提示。中間層の出力を特徴量としてRF 等に与えています。この実装によりどの学習モデルでも精度向上を成し遂げた点に価値を見出せます。
コードはGitHubで公開されています。読んでみましたが、非常にシンプルで難しいことはされていません。

（良し悪しは別として）こういった方法で実施しましたという報告も参考になります。

Finally, the landslide inventory map was converted into raster data with a grid size of 25 x 25 m. A simple and commonly used sample strategy was used in this study . 70% (255) landslide locations were randomly selected as a training set, and the other 30% (109) landslides were used for validation. In addition, to maintain class banlance, the same number of non-landslide locations (255 and 109) were randomly drawn from areas free of landslide to construct the training and testing sets.

素因を利用したシンプルな報告なのですが、国内での地すべり（地形、すべり面など）に対する機械学習アプローチの報告との差を感じる文献でした。

斜面崩壊と振動

斜面崩壊の振動伝播をシミュレーションで再現できないかと考えていました。

発生源については数年前からいろいろな方にアドバイスをいただいていました。基本、伝播とは別ソフト（DEMなど）で作るほうが良いようです。特に、機械分野の方々はこの意見でした。機械分野や自動車分野では振動解析が死活問題になりかねないのでお詳しいのでしょう。

伝播部分に関しては、いくつかのソフトがあります。当初は GMS を利用するつもりでしたが、プロ曰く、発生源入力部分が難しいようでした。

振動波形とシミュレーションを用いて崩壊のメカニズムを検証する試みは、これまでいくつかなされています。
残念ながら、最近まで知りませんでした。考え選択する方法は、どなたも同じようです。それだけスタンダードな方法になるのでしょう。

最近読んだのはコチラ↓

C. Yi-fei, Y. Yan, G. Jian, et al., Landslide reconstruction

using seismic signal characteristics and numerical simulations: Case study of the 2017 “6.24” Xinmo landslide, Engineering Geology (2019)

• Stationary stage:
• Before the landslide starts.
• random background noise

• Slipping stage:
• The slipping stage corresponds to the first principal signal.
• high frequency and high energy
• Numerical simulations indicate that the landslide body moves rapidly along the slip surface with a maximum velocity of 47 m/s.

• Transition stage: 
• micro-seismic signals
• low amplitude and low-frequency
• The DEM simulations indicate that part of the landslide body is deposited on the slopewash while the landslide head continues to move downward along the slope.

• Entrainment-transportation stage:
• The amplitude of the seismic signal first increases then decreases, and the frequency range first widens then narrows. 
• Numerical simulations indicate that the landslide energy at this stage gradually increases and then decreases.

• Deposition stage:
• The signal amplitude and peak frequency gradually decrease and the frequency range gradually narrows, reflected by the gradually decreased energy of the landslide body from the simulations.
• As the landslide body moves toward the valley bottom, the terrain becomes flatter, the migration path widens, and the movement rapidly slows down, resulting in deposition of the landslide mass.

発生前から停止まで、5つのステージに区分されています。
波形と数値シミュレーションを用いた内容ですが、これを読んだだけではよくわからない部分がありました（孫引きはまだ不完全）。

Feng, Z.Y., Lo, C.M., Lin, Q.F., 2017. The characteristics of the seismic signals induced by landslides using a coupling of discrete element and finite difference methods. Landslides, 14(2), 661-674.

• The rupturing, separation, and sliding of rock mass from upslope (source area, Pt. A and B) generate larger amplitude and lower-frequency signals than those of debris flowing at the downslope (Pt. C and D).
• The simulated velocity and acceleration are significantly larger than those field data of Station SGSB due to direct impact and near-field effect. Also, more high-frequency waves exist in the simulated signals than those field data because of the same reason, the near-field effect, and because that high-frequency waves attenuate very fast when traveling to the far-field Station SGSB.
• A landslide with larger rock particles generates lowerfrequency content seismic signals; i.e., larger particles led to seismic signals with lower frequencies and also higher amplitude and stronger signals. For further study, we can observe whether lower frequency in the signal corresponds to large size rock fragments to explain landslide characteristics.

こちらは FLAC ＋PFC 2D 。先の文献とは異なる結果ですが、これはこれで納得。

シミュレーションで波形を（概ね）再現できて初めて、議論のステージに立てます。波の解釈のみ、あるいは地質屋さんが語るメカニズムのみでは議論になりません。時間の無駄とまでは言わないですが、ベクトルがズレており遠回りです。

いずれにしても、きっかけを掴むことができました。ようやくスタートラインに立てた気がします。

FloPy

FloPy Ver.3.3
https://www.usgs.gov/software/flopy-python-package-creating-running-and-post-processing-modflow-based-models
https://github.com/modflowpy/flopy

MODFLOW を Python で動かすパッケージです。USGS の HP でも紹介されていました。気づきませんでしたね。
Jupyter Notebook を利用する example が豊富です。これは GUI より過程を追いやすいですね。
MODFLOW を使いたかったわけではないので今回は眺めるだけでしたが、それでも多くの機能を扱えるツールであることはわかりました。

FloPy での地質の適用について調べていて、目に留まったのがコチラ↓
https://www.hatarilabs.com/ih-en/how-to-insert-a-3d-geology-into-a-modflow-model-with-python-and-flopy-tutorial
DL したデータには、Neural Net を使用して地質を推定する Notebook が含まれていました。 前処理に標準化も組み込まれていました。プチ衝撃ですね。こんなところにまで機械学習が適用されていようとは。
動かしてみると、学習毎に分布の予測が変化します。ランダムシードを固定すべきところがあるのでしょう。ま、固定しなくても "sequential" NN Model として透水係数場を選定するのに使えないことはないのですが、結果を見る限り日本人受けしないでしょうね。

この例では最後に PyVista といった VTK のラッパー？を使われていました。FloPy と同様に多機能です。凄い時代になりました。

知らない間に技術は進んでいます。遅れ気味ですが、その変化が面白いと感じます。将来が楽しみです。

透水係数の入替え処

PHREEQC を Python で動かしてみました。

IPhreeqc と PhreeqPy を試してみましたが、残念ながら動いたのは前者のみ。

で、動かして気づいたのですが、これ、２,３の species に着目するだけではダメですね。輸送部分を FiPy に任せるつもりだったのですが、できた solution をそのまま次のステップに持っていく必要があります。例えば、calcite の溶解を扱うにしても pH、CaCo3、Ca2+、CO2 のみではダメ。HCO3-、 CaHCO3+、CO3-2はもちろん、その他の elements, species のみならず redox や温度などの条件も次のステップに持っていかなければなりません。先の文献のように、簡単ではありませんでした。

そうなるとこの種の問題は PHAST などで輸送まで計算しておいて、ステップ間にて透水係数を変化させた方が簡単そうですね。FiPy での物理量の入れ替えまで確認できていたのですが、見通しが甘かったようです。

ま、バカなりに手を動かして、見通しを得られただけでも良しとしましょう。

FiPy

FVM を利用した PDE ソルバーです。

FiPy 3.4
https://www.ctcms.nist.gov/fipy/

良いですね。項に対応する関数（fipy.terms package）を追加するだけで PDE を解くことができます。昔、先輩に紹介いただいたソフトに近いですね。MATLAB も含め有償のソフトはこのように簡単なのでしょうか？というか、これが10年以上前に既に組まれていたのですから驚きです。

https://www.ctcms.nist.gov/fipy/documentation/FAQ.htmlより引用

例えば、拡散のみであれば以下のように書きます。
eq = DiffusionTerm(D)
非定常だと
eq = TransientTerm() == DiffusionTerm(D)
移流項を加えると
eq = TransientTerm() + AdvectionTerm(U) == DiffusionTerm(D)
風上差分と同様の解方は、convection にありました。
eq = TransientTerm() + UpwindConvectionTerm(U) == DiffusionTerm(D)
advection, convection は使い分けられていますね。違いはコチラ↓
https://physics.stackexchange.com/questions/24489/what-exactly-is-the-difference-between-advection-and-convection

これでソース項を除く式の表現ができました。地下水の流速を求めるなら別途式が必要ですが、1文で移流拡散の離散化を任せられるのは手軽ですね。
さらに、∇は以下で表現されていますので、任意の項を作成・追加可能です。
().divergence

コチラの方は日本語で解説されています。
http://penguinitis.g1.xrea.com/study/FiPy/FiPy.html
短いですね。これで２次元や３次元の近似解を求めることができるのですから、手軽としか言いようがありません。


では、試行。
公式のexample と 上記の方のコードを動かしながら、３次元まで拡張です。ここまで短く記述できるのは素晴らしい。
https://www.ctcms.nist.gov/fipy/examples/README.html

from fipy import *

nx = 100
ny = 100
nz = 10
Lx = 1.
Ly = 1.
Lz = 1.
dx = Lx/nx
dy = Ly/ny
dz = Lz/nz

mesh = Grid3D(dx, dy, dz, nx, ny, nz)
phi = CellVariable(name='phi', mesh=mesh, value=.5)

D = 1.
U = (3., 2.,1.)
phi1 = 1.
phi2 = 0.5

phi.constrain(phi1, mesh.facesLeft)
phi.constrain(phi2, mesh.facesRight)

eq = TransientTerm() + UpwindConvectionTerm(coeff=U) == DiffusionTerm(coeff=D)

T=1.
dt=0.1
for step in range(1, int(T/dt)+1):
    print('time step:', dt*step)
    eq.solve(var=phi, dt=dt)
    TSVViewer(vars=phi).plot(filename=str(step)+".tsv")
    VTKViewer(vars=phi).plot(filename=str(step)+".vtk")

結果、2次元まではOKですが、3次元では異常に時間がかかりました（実務では厳しいくらいに）。ラップトップで12分。この程度の計算で12分もかかるのはダメ。ソルバーの変更か並列化が必須でしょう。

Mayavi が動かなかったのと plotly での表示が遅かったので、可視化は別ソフトに頼ることを念頭に VTK 書き出し（plotly での表示は以下）。

import pandas as pd
import plotly.express as px
import plotly.graph_objects as go

df = pd.DataFrame({'x' : mesh.x, 'y': mesh.y, 'z':mesh.z, 'value': phi})

fig = go.Figure(data=go.Volume(
    x=mesh.x,
    y=mesh.y,
    z=mesh.z,
    value=df.value,
    opacity=0.5,
    surface_count=20
    ))

fig.update_layout(margin=dict(l=0, r=0, b=0, t=0))
fig.show()

phi が FiPy のオブジェクトになっていますので、やや扱いづらい。ま、df にはなるので書き換えが可能ならPHREEQCとの連成も容易でしょう。

鉱物沈殿による透水係数の低下

Ali et al., A coupled geochemical and fluid flow model to simulate permeability decline resulting from scale formation in porous media
https://www.sciencedirect.com/science/article/abs/pii/S0883292719301544

PDEをMATLABで、沈殿による溶質の減少をPHREEQC（IPhreeqc module）で解いています。反応速度と流速は局所平衡を仮定。式だけ見ると溶解も扱えますが、フローは沈殿のみでした。
透水係数の変化はパラスタ。これまでの手法では反応部分も含め2個のパラメータ同定が必要だったようですが、PHREEQCとの連成で1個に減ったとのこと。沈殿による透水係数の変化部分と弱連成に飛びついたのですが、やはり試験や実例がないと合わせ込みは難しいようです。ま、当然ですね。

今まで見た文献で連成部分のコードを公開しているものは確かなかったと思います。国家機密として公開されなかったり、請求しても返答なかったり。
今回も同様。 PDE で濃度を求め、IPhreeqc に渡して SI と 濃度を更新。SI＞０で透水係数を更新して次のタイムステップへ。これらを時間分 For Loop させているのでしょう。1次元なのでこれで解けます（BESTかはわかりませんが）。

MATLAB の代替えとして Python を利用している報告を目にします。今回の PDE 部分も代替えできるのでしょうか？ PHREEQC は Python で利用している例を時々見ますので、今回の文献程度のレベルであればで実装できるかもしれません。調べてみましょう。

土砂の停止位置

土砂が斜面を流下し平坦面に達した後、どの程度の距離で止まるか？

これ、計算で求めるのは難しいですよね。
計算事例が少なく、これといった一般値がまだありません（ある程度の目安はありますが）。パラスタで合わせに行くにしても、いくらか時間がかかります。

そんなことを考えながら、ふと思い出したのがコチラ。図解法で停止位置のあたりを付けられます。
國生剛治「地震地盤動力学の基礎　エネルギー的視点を含めて」

が、読み直してまた思い出しました。そういえばこれも同じでした。静止摩擦と動摩擦、内部摩擦の区別がない手法でした。
https://phreeqc.blogspot.com/2015/02/blog-post_24.html

これ、崩壊後に斜面の途中で止まっていた土砂が、仮に途中で平坦面に達していたらという事例だと簡単に位置を求められそうです。が、斜面上の停止位置を順計算のみで求めるのは「これで大丈夫」が言えないのでダメですね（少なくとも実務では）。
天然ダムを形成するような”みずみずしい”、”途中で正常の変化する”土砂もダメ。感覚的には乾いた土砂の方が合うと思います。ん？だから地震時なのでしょうか？

過去の事例を集め一斉に計算すれば、適用できる条件や適したツール、パラメータを求められるでしょう。
どこかの学会でやっていただけないでしょうか？