文献：LSM uncertainty analysis

PyLandslide: A Python tool for landslide susceptibility mapping and uncertainty analysis - ScienceDirect

AI要約

背景
地すべりリスク評価において、専門家の主観的判断による重み付けが一般的だが、バイアスや不確実性が問題。
機械学習を用いた統計的アプローチが発展しているが、重みの不確実性の定量化が不十分。
既存ツール(LSM Tool Pack、LSAT PM等)は存在するが、不確実性分析機能が限定的。

手法
1.ジニ不純度に基づく特徴量重要度(重み)の計算
1)ランダムフォレスト内の各決定木において、ある特徴量で分岐する際のジニ不純度の減少量を計算
2)全ての木について、各特徴量による不純度減少量を合計
3)各特徴量の重要度 = その特徴量による不純度減少の合計 / 全特徴量による不純度減少の総和

2.不確実性分析の手法
20,000回の反復計算を実施
各反復で地すべり・非地すべり地点の80%をランダムサンプリング
accuracy 0.75以上のモデルのみの特徴量重要度（feature importance）を記録
重要度の範囲を集計し、限られたインベントリデータに起因する不確実性を定量化

3.LSI（Landslide Susceptibility Index）の算出
ある地点が地すべりに対してどれだけ脆弱であるかを数値化した指標
PyLandslideでは、複数の要因を重み付けして統合することで、空間的な地すべり発生リスクを評価

LSI = Σ(Wi × Si)

LSI: 地すべり感受性指数
Wi: 要因iの重み（weight）
Si: 要因iのスコア層（score layer）
n: 考慮する要因の総数

スコア層（Si）の算出方法
連続データ型要因（降水量、傾斜、道路距離など）：
分位数による分類: データを11個の分位数（quantile）に分類
地すべり件数の計算: 各クラス内で発生した過去の地すべり件数を集計

降水量・傾斜: そのクラス以下で発生した地すべりの累積パーセンテージ
道路距離: そのクラス以上で発生した地すべりの累積パーセンテージ
降水量が17~80mmのクラス: 80mm以下で発生した地すべりの累積割合がスコア
道路距離が143~286mのクラス: 143m以上の距離で発生した地すべりの割合がスコア

カテゴリカルデータ型要因（土地被覆、岩相など）：
各カテゴリ内で発生した地すべり件数を計算
各カテゴリのスコア = そのカテゴリ内で発生した地すべりのパーセンテージ

標準化
全ての要因について、スコアを0~100の範囲に線形正規化:
0: その要因クラスの最低寄与度
100: その要因クラスの最高寄与度

LSI分類基準
Very low: 0 ≤ LSI ≤ 20
Low: 20 < LSI ≤ 40
Moderate: 40 < LSI ≤ 60
High: 60 < LSI ≤ 80
Extremely high: 80 < LSI

感度分析の実施回数
歴史的降水条件について200回のランダムな重みの抽出を実施
9つの気候予測について各200回ずつ実施
内訳：
歴史的降水条件（1981-2023年）：200回
将来予測（2041-2050年）：9つの気候シナリオ × 各200回 = 1,800回
合計：200 + 1,800 = 2,000回

9つの気候シナリオ：
3つの気候モデル（ACCESS-ESM1-5、MRI-ESM2-0、UKESM1-0-LL）
3つの共有社会経済経路（SSP126、SSP245、SSP585）
3モデル × 3シナリオ = 9つの組み合わせ
この結果、LSIも範囲として表現される

結果
イタリアにおける6要因の重要度(重み)
道路からの距離: 0.43〜0.52(中央値0.47)
・最も影響が強い
・過去の地すべりの85%が道路から143m以内で発生

傾斜: 0.16〜0.23(中央値0.20)
・2番目に重要
・過去の地すべりの70%が傾斜9%以上の地域で発生

土地被覆: 0.10〜0.13(中央値0.12)
地質: 0.07〜0.09(中央値0.08)
降水量: 0.05〜0.07(中央値0.06)
TWI(地形湿潤指数): 0.05〜0.08(中央値0.06)

モデル精度
訓練データでの全体精度: 0.80〜0.82
テストデータでの全体精度: 0.75〜0.80
過学習は認められない

ｌSI評価
歴史的期間(1981-2023):
「極めて高い」: 7.8〜9.5%
「高い」: 23.8〜26.8%

将来予測(2041-2050):
「極めて高い」: 5.3〜7.6%(減少傾向)
「高い」: 21.5〜28.3%

イタリア北西部と南部で大幅に低下
北東部で増加

考察
人為的要因(道路)が最大の影響
道路建設に伴う斜面切土や人間活動が斜面安定性を低下

不確実性の重要性
重みの不確実性範囲を考慮することで、より堅牢な意思決定が可能
インフラ投資計画において、重みの不確実性を考慮した堅牢な地すべりリスク評価が可能

研究の限界
気候変動の影響評価
降水量変化のみを考慮、温度・山火事・植生変化は未考慮

入力データの精度
ITALICAデータベースの位置精度は最大5.6kmの誤差
より高精度なデータが望ましい

閾値の設定
accuracy0.75の閾値は本研究で設定したが、絶対的な基準ではない
ユーザーのニーズに応じて調整可能

国外では、LSM作成において機械学習を用いるアプローチは既に一般化していると言えますが、そこに不確実性評価が入っています。といっても、重要度の示し得る範囲を利用しているだけですから、新たなツールは不要です。
この文献のLSIの定義が正解かどうかはわかりませんが、少なくとも確率を導入する方針については同意です。

AIの利用法

【優勝🥇】防衛省サイバーコンテストをAIで攻略した話 #Security - Qiita

このような時代なのですね。AIで攻撃、AIで防衛。

勝敗はAIで解けない超高難度の問題で決まります。

うーん。考えさせられます。

文献：Probabilistic rainfall thresholds for landslide occurrence using a Bayesian approach

 Probabilistic rainfall thresholds for landslide occurrence using a Bayesian approach - Berti - 2012 - Journal of Geophysical Research: Earth Surface - Wiley Online Library

AI要約＋α

背景
従来の地すべり予測手法は、地すべりを引き起こした降雨イベントのみを考慮し、決定論的閾値を提供してきた。しかし、斜面安定性は降雨だけでなく、多数の要因の組み合わせで決まるため、同じ降雨条件でも異なる結果(地すべり発生/非発生)が生じる。経験的モデル:降雨強度-継続時間(ID)閾値が一般的だが、不確実性の定量化が不十分である。

研究の対象地域
イタリア・エミリア=ロマーニャ州の山岳地域(12,000 km²)
1939年以降の4,000件以上の地すべり履歴データ
176箇所の雨量計による日降水量データ

手法
1. ベイズの基本手法(1次元解析)
P(A|B) = [P(B|A) × P(A)] / P(B)
P(A|B): 事後確率 - 降雨Bが発生したときに地すべりAが起こる確率(求めたい値)
P(B|A): 尤度 - 地すべりが発生したときに降雨Bが観測される確率
P(A): 事前確率 - 降雨に関係なく地すべりが発生する確率
P(B): 周辺確率 - 地すべりに関係なく降雨Bが観測される確率

実際の計算では、以下のように相対度数で近似:
P(A) ≈ NA / NR  (地すべり発生数 / 全降雨イベント数)
P(B) ≈ NB / NR  (降雨Bの発生数 / 全降雨イベント数)
P(B|A) ≈ N(B|A) / NA  (降雨Bで発生した地すべり数 / 全地すべり数)P
(A|B) ≈ N(B|A) / NB

一変量例（強度 I>40 mm/日）
・事前 landslide 確率：P(A)=5/20=0.25
・強度大の割合：P(B)=P(I>40)=9/20=0.45
・発生時の強度大割合（尤度）：P(B|A)=4/5=0.8

→ P(A|I>40)=P(B|A)·P(A)/P(B)=0.8·0.25/0.45≃0.44

降雨B発生時の地すべり確率は44%である。これは尤度P(B|A)=80%とは異なることに注意が必要。事前確率と周辺確率を考慮する必要がある。

2. 2次元ベイズ解析
2つの変数B(降雨強度I)とC(降雨継続時間D)を考慮:
P(A|B,C) = [P(B,C|A) × P(A)] / P(B,C)
ここで、P(B,C)は2変数の同時確率(ある強度と継続時間の組み合わせが観測される確率)を表す。

各セルでの確率計算
二変量例（I>50 mm/day かつ D<=1 day）
・事前 landslide 確率：P(A) = 5/20 = 0.25
・降雨条件を満たす割合：P(B,C) = P(I,D) = 4/20 = 0.20
・発生時の当該領域割合（尤度）：全地すべり5件中、この条件で発生2件
   P(B,C|A) = P(I,D|A) = 2/5 = 0.40

→ P(A|I>50, D≤1)= P(A|B,C) = [P(B,C|A) × P(A)] / P(B,C) =0.40·0.25/0.20=0.50

I-D平面全体の確率分布ヒストグラムを作成。
強度と継続時間の相互作用効果を把握。
確率の空間分布パターンを視覚化。
高リスク領域の明確な識別。

結果
1次元ベイズ解析の結果
有意な変数:総降雨量E、継続時間D、平均強度I
P(B|A)とP(B)の分布が明確に異なる
P(A|B)が事前確率P(A)=0.005を大きく上回る
特に降雨強度Iが最も有意:I>100mm/dayでP(A|I)=0.28

非有意な変数:先行降雨AE14, AE30
P(B|A) ≈ P(B)
P(A|B) ≈ P(A)

研究地域では先行降雨は地すべり発生と相関が低い。
降雨の激しさ(総量、継続時間、強度)とともに地すべり確率が増加。
ただし、極端な値では確率が減少する傾向(サンプルサイズの問題と定義バイアスによる)。

2次元ベイズ解析の結果
特定のI-D値で確率が急激に増加。これはシステムの状態の根本的変化、物理的閾値の存在を示唆。不確実性を含めた警報システムの構築に有用。

考察

1. ベイズ手法の利点
統計的厳密性：事前確率と周辺確率を考慮:尤度P(B|A)だけでなく、降雨の頻度P(B)と地すべり発生率P(A)を組み込む。
認知バイアスの回避:人間の直感的判断は事前確率を無視しがち(「80%の地すべりがI>50で発生」≠「I>50で地すべり確率80%」)。
不確実性の定量化:0〜1の連続値で確率を表現、決定論的手法の曖昧性(「閾値超過時に何が起こるか?」)を解消、信頼区間により推定の信頼性を評価

2. 先行降雨の非有意性
一般に細粒土壌では先行降雨が斜面安定性に重要とされるが、本研究では相関が低い。研究地域の地すべりは主にイベント降雨に対する急速な水文応答が支配的。60%の地すべりが降雨終了時またはその直後に発生。
深層地すべりでも、長期的な要因よりも短期的な降雨イベントが引き金となる。

3. 物理的閾値の示唆
I-D平面の特定領域で地すべり確率が急激に上昇。これはシステムの状態変化(安定→不安定)を示唆。
確率的手法であっても、背後に物理的メカニズム(臨界間隙水圧など)が存在。確率分布のパターンから、真の物理的閾値の存在を推測可能。

4. 方法論的課題と解決策
多発地すべりの扱い
方針:同一降雨による複数地すべりは1イベントとしてカウント。
理由:P(A) = NA/NR > 1を避け、統計的整合性を保つ。

バイアスへの対応
問題:誘因降雨は地すべり発生日で打ち切られるが、非誘因降雨は降雨終了まで継続。
影響:極端な長継続時間での確率推定に影響。
緩和:信頼区間の提示、データの60%が降雨終了時に発生することで影響を軽減。

Landslide発生降雨をベイズ手法で評価した少し古い文献です。基本的には、国内で研究されている降雨数ベースの評価と変わりません。
国内では短期雨量と土壌雨量指数の２軸で評価されます。スネーク曲線を書いてピークで評価したりするのですが、原点から離れるにつれてレアな雨（P(B,C) ≒０）となり、確率が極端に高くなります。これ、補正したいですよね。

文献：Contact modelling approach to simulate water flow in a rock fracture

Contact modelling approach to simulate water flow in a rock fracture under normal loading – Modelling report of Task 10.2.2. Task 10 of SKB Task Force GWFTS – Validation approaches for groundwater flow and transport modelling with discrete features – SKB.com

AI要約

背景
未開口の自然岩盤亀裂に対し、法線荷重を変化させたときの流量を「ブラインド予測」する課題。
亀裂を有する岩試料で水理‐力学試験を実施後、亀裂を開口し両面を高精度スキャン、その後再装着して再試験。
目的：単一亀裂での荷重‐変形‐流れの関係と流路のチャンネリング、ならびに実務的な検証手法の検討。

手法
使用コード：COMSOL Multiphysics。
流れ：低レイノルズ数の定常流とし、Navier–Stokes式から慣性項を除いたStokes流として 3D 解析。
岩盤変形：線形弾性・等方・連続体とし、亀裂面の接触を「拡張ラグランジュ法」による接触条件で表現。

モデル群：
Analytical Model：平行平板間流れ（立方則）でスキャン精度が流量に与える影響を評価。
Parallel Plate Model：3D平行平板の数値解で解析解との整合とメッシュ条件を確認。
Small Surface Model：小さな実測亀裂面（Task 10.2.1）を用い、上面のx,y,z変位とz軸回転の影響をDSDによる感度解析。
Nearing Contact Model：上下岩体を初期ギャップから押し付ける接触モデル（1 MPa荷重）。
Departing Contact Model：初期的に重なりを与えた状態から上盤を押し上げる接触モデル（0,1,2,4,6,8 MPa）。
接触オフセット：数値的な「めり込み」を避けるため、上下面の最小間隔を強制（0.05–0.1 mm）。この仮定が開口と流量を増大。

結果
Analytical / Parallel Plate：開口 = スキャン精度0.035 mmに対応する流量は約0.1 ml/sで、実測流量と同程度。実測から逆算した等価平行平板開口はいずれもスキャン精度以下。
Small Surface Model：x,y,z変位はいずれも統計的に有意で、とくにz方向変位が流量に最も強く影響。z軸回転はほぼ無影響。
Nearing / Departing Contact：接触解析から得た変形後空隙でStokes流解析を実施。

いずれの荷重・接触アプローチでも、予測流量は実測より数桁大きい。接触オフセットやスキャン精度、初期位置誤差、補間誤差、弾性仮定を考慮しても、この差を説明できない。
von Mises応力は局所的に一軸圧縮強度を超え、現実には局所破壊や塑性変形が起こる可能性が高いことを示唆。

考察
スキャン精度0.035 mmでは、実測レベルの微小流量を信頼性高く再現するには不足。逆算開口が常に精度以下であるため、幾何データに根本的な限界。
上下亀裂面の水平方向移動（x,y）とより精密なインターロッキング、ならびに岩の塑性変形・破壊を無視しているため、実際より開口が大きく評価され、流量を過大評価している可能性が高い。
接触オフセットは数値上必要だが、同時にモデルに系統的バイアス（過大開口）を導入する。
既知の不確実性を考慮しても、モデルと実験の乖離は説明できず、「未知の支配的要因」が欠落していると結論。
接触モデル＋詳細幾何に基づく手法は、本課題のような極めて小さい流量の予測には実務的でなく、より小さい試料や高精度スキャン、あるいは別の概念モデルが必要と提言。
併せて、複数現象・複数検証を通じたモデルの「妥当性ランクR」を定義する検証評価指標を提案し、今後のモデル群の統合的評価に用いることを示唆。

0.035mmで精度が不足とは、驚きです。いえ、このような制度が出せるのも驚きなのですが。地層処分ですので細かい。
使われた 3 つの接触法（ペナルティ、拡張ラグランジュ、Nitsche）、どれも知りませんでした。拡張ラグランジュが最もオーバーラップ（食い込み）が小さいものの計算コストが高いそうです。オーバーラップを完全にはゼロにするには非常に細かいメッシュが必要で、現実的なメッシュでは 0.01–0.1 mm 程度の食い込みが残ると報告されていますが、それでも完全には抑えきれなかったとのこと。扱うスケールが全く異なります。

文献：Site scale modelling of groundwater evolution

Site scale modelling of groundwater evolution at SFR. A modelling application with iDP – SKB.com

AI要約

背景
スウェーデン核燃料・廃棄物管理会社（SKB）は、Forsmarkにある低・中レベル放射性廃棄物の最終処分場SFRの拡張を計画している。地質学的な処分の候補地の信頼性ある評価には、地下水の流動（hydrogeology）と地下水化学（hydrochemistry）を統合したモデル構築が必要である。しかし、数値的・概念的な制約や計算負荷のため、多くの従来研究ではこれらの重要なプロセスを個別に扱ってきた。本研究では、SFR周辺の地圏の進化を記述する3次元反応輸送モデルの構築と検証を目的とし、地下水流動と化学反応を統合したモデル化に取り組んだ 。

手法
本研究では、DarcyToolsとPFLOTRANという2つの計算コードを連結する「iDP（interface DarcyTools-PFLOTRAN）」を活用した。DarcyToolsは複雑な地層中の流れ・物質輸送を、PFLOTRANは高性能計算環境での流れ・反応輸送をそれぞれ扱う。DarcyToolsで得られた透水性・圧力場等をPFLOTRAN形式に変換し、反応輸送シミュレーションを行う。地圏を等価連続多孔質媒体（ECPM）として上流から下流までを模擬し、rock matrix領域は「multi-continuum approach（Lichtner法）」により反応性を持たせている。
モデルの検証のため、1D、2D、多孔質岩石行列を想定したモデルや解析解を利用。化学反応は前研究FASTREACTとの比較のためベースケース（方解石＋ヘマタイト平衡）とバリアントケース（方解石＋非晶質硫化鉄(FeS(am))平衡）の2通りを設定。3500年間（2500〜6000AD）の沿岸線移動および気候変動・地下水組成の変化も考慮した。

結果
3D反応輸送計算により、岩石マトリクス領域を含むモデルと含まないモデル（FASTREACTと同水準の設定）の比較が可能となった。
塩化物（Cl）の分布は、希釈・海岸線移動による水侵入の影響をよく再現し、岩石行列の拡散効果による急激な化学変化の緩和が見られた。
pH値は両ケースとも7.2〜8で安定。pe値、還元状態も3500年の計算期間で安定して保持された。
Fe2+等については、岩石マトリクスでの biotite 溶解による供給があり、反応的な鉱物（ヘマタイトあるいはFeS(am)）の存在条件によって、濃度推移が異なった。
主要カチオンの動態、水理パターン、時空間分布についても詳細な可視化を3Dで実現した。

Fe^2+の濃度進展には岩石マトリクスの反応性の影響が顕著であった。バリアントケース（FeS(am)平衡）ではFe^2+の濃度はベースケース（ヘマタイト平衡）に比べて低下し、FeS(am)の沈殿が効果的なFe^2+の吸収源＝シンクとなっていることが示された。 

Cl^-のような保存的成分（反応しない成分）は水収支および拡散・移流の指標として有用であり、モデルでの到達遅延や空間分布が確認された

考察
岩石マトリクス（そのパラメータ）の設定によって、地下水組成の急激な変化（例：希釈効果、流動路の変化）が拡散・反応により緩和されることが示された。
従来の1DモデルやFASTREACTによる結果と、今回の3Dモデル（岩石マトリクスを加味したもの）の比較で、岩石マトリクスを考慮することの重要性が強調された。
モデルパラメータ（特に岩石マトリクスの空隙率、表面積等）が地質環境ごとに大きく異なるため、今後は現場データによるパラメータ同定が重要になる。
開発した3Dモデル（iDP）は、今後より複雑・広範なスケールでの安全評価や、放射性核種の挙動・バリア相互作用のモデル化が期待される。

亀裂（fracture）の移流分散と、岩石マトリクス拡散を解いているようです。岩石中の拡散パラメータを設定するのは難しいでしょうね。
反応計算の地球化学的パラメータは前研究FASTREACT（Román-Ross et al. 2014）のモデルに準拠とのこと。どのようなものでしょうか。

1次元移流分散（有限差分法、Python）

GitHubで水理公式集例題集2024 のコードを眺めていた時に、地下水の不飽和浸透や移流分散計算がUPされているのを見かけました。

EXCEL版 と Fortran版 があり、懐かしいなあと思いつつ「例題 1.24 水平一次元帯水層を対象とした溶質輸送過程の計算（数値解）」の Fortran コードを Python に移植してみました。

で、結果が全く同じだったことには少し驚かされました。

以下、コードに沿ったFTCS（Forward-Time Centered-Space）による離散化とスクリプトです。手元に例題集が無いため、表現は異なるかもしれません。スクリプトではスライスを利用することで計算が1行となり、見やすくなりました。実際の地下では流速や流量が一定ではないので、ここまで簡素にはできませんが、まずはここからなのでしょう。

1. 支配方程式

∂c/∂t + v·∂c/∂x = D·∂²c/∂x² ……(1)

ここで、

c : 濃度 (mg/L)

t : 時間 (day)

x : 距離 (m)

v : 移流速度 (m/day)

D : 分散係数 (m²/day)

2. 差分近似

2.1 時間微分の前進差分

時刻 t における時間微分を前進差分で近似すると、

∂c/∂t ≈ (c_i^(t+1) - c_i^t) / Δt ……(2)

ここで、

c_i^t : 格子点 x_i、時刻 t における濃度

Δt : 時間ステップ

2.2 空間2次微分の中心差分（拡散項）

格子点 x_i における空間2次微分を中心差分で近似すると、

∂²c/∂x² ≈ (c_(i+1)^t - 2·c_i^t + c_(i-1)^t) / Δx² ……(3)

ここで、

Δx : 空間格子間隔

2.3 空間1次微分の後退差分（移流項）

v > 0（下流方向の流れ）を仮定し、上流差分（upwind差分）を適用すると、

∂c/∂x ≈ (c_i^t - c_(i-1)^t) / Δx ……(4)

3. 差分方程式の導出

式(1)に式(2)、(3)、(4)を代入すると、

(c_i^(t+1) - c_i^t) / Δt + v·(c_i^t - c_(i-1)^t) / Δx

= D·(c_(i+1)^t - 2·c_i^t + c_(i-1)^t) / Δx² ……(5)

式(5)の両辺に Δt を掛け、c_i^(t+1) について解くと、

c_i^(t+1) = c_i^t

+ (D·Δt/Δx²)·(c_(i+1)^t - 2·c_i^t + c_(i-1)^t)

- (v·Δt/Δx)·(c_i^t - c_(i-1)^t) ……(6)

ここで、

r_diff = D·Δt / Δx² （拡散数）

r_adv = v·Δt / Δx （クーラン数）

を用いると、式(6)は以下となる。

c_i^(t+1) = c_i^t

+ r_diff·(c_(i+1)^t - 2·c_i^t + c_(i-1)^t)

- r_adv·(c_i^t - c_(i-1)^t) ……(7)

import numpy as np
import csv

# +++++++++++++++++++++++++++++++++++++++++++++++
# 1次元移流分散方程式の数値解析
# 有限差分法：FTCS（Forward-Time Centered-Space）
# ∂c/∂t + v ∂c/∂x = D ∂²c/∂x²
# +++++++++++++++++++++++++++++++++++++++++++++++

def reidai082():
    # ----- パラメータ設定 -----
    xmax  = 30.0        # 解析領域の長さ (m)
    dx    = 1.0         # 差分格子間隔 (m)
    dx2 = dx ** 2     # 差分格子間隔の二乗
    imax = int(xmax / dx) + 1    # 下流境界の格子点番号

    xo  = 10.0   # 濃度観測地点（m）
    ii1 = int(xo / dx)   # Python 0-based インデックス

    vel   = 0.1         # 実流速 (m/day)
    al    = 1.0         # 縦分散長 (m)
    tau   = 1.0         # 屈曲率 (-)
    dm    = 1.0e-5      # 分子拡散係数 (cm^2/s)
    # 分散係数 D (m2/day)
    # dm を m2/day に換算するために 1 cm2/s = 8.64 m2/day
    D     = al * abs(vel) + tau * dm * 8.64

    bc    = 1.0         # 上流境界の濃度 (mg/L)

    tend  = 300.0       # 計算時間 (day)
    dt    = 1.0         # 時間ステップ (day)
    nmax  = int(tend / dt)   # 計算時間のステップ数

    tjikan = 1.0        # 結果をファイルに出力する時間間隔（day）
    iout1 = int(tjikan / dt)   # 結果をファイルに出力する時間間隔のステップ数

    # ----- グリッド設定 -----
    x    = np.linspace(0, xmax, imax)

    # ----- 配列の用意 -----
    c  = np.zeros(imax)  # 新しい時刻の濃度
    co = np.zeros(imax)  # １ステップ前の濃度

    # ----- 係数の用意 -----
    r_diff = D * dt / dx2
    r_adv  = vel * dt / dx

    # 初期条件：上流境界だけ bc
    co[0] = bc

    # ----- 出力ファイル準備 -----
    out_filename = 'output.csv'
    with open(out_filename, mode='w', newline='') as f:
        writer = csv.writer(f)
        writer.writerow(['elapsed-time(days)', 'concentration(mg/L)'])
        # t=0 のデータ
        writer.writerow([f"{0.0:.1f}", f"{co[ii1]:.7f}"])

        # ----- 時間ループ -----
        time = 0.0
        iout  = 0
        for n in range(1, nmax+1):
            time += dt
            iout += 1

            # ---- 境界条件 ----
            co[0]      = bc      # 上流
            co[-1]     = 0.0     # 下流

            # ---- 内部格子の更新 (FTCS) ----
            # c[i] = co[i] + D*dt/dx^2*(co[i+1]-2*co[i]+co[i-1])
            #                  - vel*dt/dx*(co[i]-co[i-1])
            # スライスで一括計算
            # co[1:-1]：インデックス 1～N−2 の内部点
            # co[2:] ：インデックス 2～N−1
            # co[:-2] ：インデックス 0～N−3

            c[1:-1] = (co[1:-1]
                       + r_diff * (co[2:] - 2*co[1:-1] + co[:-2])
                       - r_adv  * (co[1:-1] - co[:-2]))

            # 更新
            co[:] = c[:]

            # ---- 出力 ----
            if iout == iout1:
                print(f"{time:6.1f} days")
                writer.writerow([f"{time:.1f}", f"{co[ii1]:.7f}"])
                iout = 0

LiDARによる土石流観測手法の比較

土研資料とGeomorphology文献の違いをAIさんに整理してもらいました。

基本情報

項目	土木研究所資料	Geomorphology論文
文献名	二次元レーザスキャナによる土石流の流量観測手法	Monitoring debris flow dynamics: insights from 4D-LiDAR observations in Ohya landslide, Central Japan
発行年	2023年11月	2025年
発行元/掲載誌	国立研究開発法人土木研究所資料 第4445号	Geomorphology 482 (2025) 109800
著者	今森直紀、清水武志、池島剛、伊藤誠記	Tatsuki Kaneko, Fumitoshi Imaizumi, Tomoya Osada, Saleh Yousefi, Shoki Takayama
観測地	桜島有村川流域	大谷崩れ(Ohya landslide)、静岡県中部

使用機器

項目	土木研究所資料	Geomorphology論文
レーザスキャナ種類	二次元レーザスキャナ(2D-LiDAR)	4D-LiDAR(3D+時間)
機種	北陽電機 UXM-30LAH-EWA	Livox Horizon (Livox Corp.)
波長	905nm近赤外線レーザ光	記載なし
ステップ角	0.125度	記載なし
測定点数	1520点(190度走査)	記載なし
測距範囲	0.1~30m	記載なし
サンプリングレート	20Hz(50ms毎)	0.1秒間隔(10Hz相当)
測定分解能	1mm	記載なし

観測システム

項目	土木研究所資料	Geomorphology論文
走査次元	1軸走査(断面計測)	3次元走査
設置箇所数	1箇所	2箇所(Site U: 上流、Site D: 下流)
併用機器	非接触型流速計、ビデオカメラ	ビデオカメラ(Sony HDR-CX470、60fps)
電源	UPS経由の商用電源	バッテリー+ソーラーパネル
データ記録装置	ノートPC	Raspberry Pi(オンボードコンピュータ)
起動方式	常時稼働	ワイヤーセンサーによる自動起動
データ保存間隔	連続記録	1時間間隔のファイル保存

計測・解析手法の比較

項目	土木研究所資料	Geomorphology論文
計測対象	流下断面積(水通し断面)	3次元表面形状、縦断・横断形状
座標系	極座標→平面直角座標変換	極座標→直交座標変換(x:鉛直、y:横断、z:流下方向)
解析ソフトウェア	Python 3.9(独自プログラム)	CloudCompare (v2.13)、QGIS (v3.16)
ノイズ除去	1秒間の中央値抽出(雨滴除去)	記載なし
DEM作成	記載なし	複数のグリッドサイズでDEM作成
地形指標	流下断面積のみ	勾配、粗度、これらの平均値・標準偏差

流量算出手法

項目	土木研究所資料	Geomorphology論文
流速計測	非接触型流速計(表面流速×3/5)	ビデオカメラ+LiDAR(直接計測)
流下断面積	LiDARから直接算出(台形公式)	LiDARから3次元的に算出
流量計算式	Qs = U × Ad	詳細な記載なし
時間分解能	1秒(中央値処理後)	0.1秒

観測対象とする土石流の特徴

項目	土木研究所資料	Geomorphology論文
流域面積	約1.6km²	約0.37km²(市ノ沢)
チャンネル特性	砂防堰堤水通し(固定床)	自然渓流(移動床・固定床混在)
河床勾配	記載なし	上流部:2537.3°、下流部:1520°
土石流タイプ分類	記載なし	完全飽和流、部分飽和流
観測イベント数	2022年に4回発生中2回計測	2023年8月3日の1イベント詳細解析

データ処理・出力の比較

項目	土木研究所資料	Geomorphology論文
出力データ形式	CSV(時刻、平均水位、流下断面積、平均流速、流量)	点群データ、DEM、各種地形指標
可視化	ハイドログラフ、断面図、アニメーション(mp4)	縦断・横断プロファイル、地形指標グラフ
プログラム言語	Python 3.9	記載なし(CloudCompareとQGIS使用)
公開データ	プログラムソースコード+計算データ(DVD-R添付)	データ公開の記載なし
ライセンス	プログラムはCC BY-SA 4.0	論文はCC BYライセンス

研究成果・新知見の比較

項目	土木研究所資料	Geomorphology論文
主な成果	・流下断面積の定量的・高時間分解能計測 ・夜間・豪雨時の計測可能性 ・従来手法との比較検証	・完全飽和流と部分飽和流の形態差異 ・表面粗度と乱流の関係 ・堆積過程における逆勾配地形形成
流速に関する知見	表面流速×3/5を平均流速として使用	先端部より後続部の表面流速が高い場合あり
断面形状	台形断面として計算	凸型断面(部分飽和流)の観測
堆積特性	記載なし	先端部の堆積が後続流の移動性を低下させる
表面形状指標	記載なし	勾配・粗度の平均と標準偏差が流動特性を反映

技術的課題と今後の展望

項目	土木研究所資料	Geomorphology論文
データ取得の安定性	・HDD→SSD換装による改善提案 ・温度上昇対策 ・リモート監視の必要性	・霧・豪雨時の検出範囲制限 ・3m以内の点群歪み ・観測継続性の確保
流速計測の課題	・1点のみの計測限界 ・3D-LiDARへの発展可能性	・0.1秒間隔での粒子追跡の限界 ・粒径との相関の課題
維持管理	・光学窓の火山灰付着対策 ・定期的な清掃(約2ヶ月毎)	・ワイヤーセンサーの設置・維持 ・ソーラーパネル+バッテリーの管理
今後の発展	・ROSを用いた制御プログラム開発 ・降雨連動起動システム	・自然渓流での更なるデータ蓄積 ・複数地点での同時観測