2026年11月1日日曜日

文献:GeoDualSPHysics: a high-performance SPH solver for large deformation modelling of geomaterials with two-way coupling to multi-body systems

GeoDualSPHysics: a high-performance SPH solver for large deformation modelling of geomaterials with two-way coupling to multi-body systems

AI要約

1. 背景(Background)

  •  従来の課題: 地盤工学や防災工学における地滑り、土石流、斜面崩壊などの現象は、土砂が液体のように振る舞う「極大変形(大流動)」を伴います。従来の有限要素法(FEM)などの格子基盤の手法では、メッシュが過度に歪むため、このような大変形を安定して計算することが困難でした。
  • 粒子法(SPH)の台頭と限界: メッシュを用いない「粒子法(SPH)」は大変形問題に適していますが、計算コストが非常に高く、数百万〜数千万粒子に及ぶ広大な実地形の大規模シミュレーションを現実的な時間で処理することが困難でした。
  •  開発の動機: 既存の非常に強力なオープンソースGPU加速型流体ソルバー「DualSPHysics」の並列計算フレームワークを基盤とし、地盤固有の力学挙動(弾塑性)や構造物との複雑な相互作用を、数億粒子スケールで超高速に計算できる実用的な地盤用ソルバーの確立が必要とされていました。

2. 手法(Methodology)

  • •連続体力学のラグランジュ定式化:地盤材料を連続体として捉え、質量保存則および運動量保存則をSPHの粒子離散化形式で解きます。カーネル関数には、粒子の異常な凝集(クランピング不安定性)を防ぎ安定性に優れる「5次Wendland(C2)カーネル」が採用されています。
  • 地盤の構成モデル(Constitutive Models):材料の応力状態を評価するため、弾塑性モデル(例:Drucker-Prager Yield Criterion など)を組み込んでいます。毎タイムステップで、まず弾性試行応力を計算し、降伏条件を超えた場合には「リターンマッピングアルゴリズム(Return Mapping Algorithm)」を用いて、応力を降伏曲面上に引き戻すプラスチックひずみ計算を行います。
  •  ノイズ除去応力処理技術(Noise-free Stress Treatment):従来のSPHで課題となっていた「応力場の数値的なノイズや振動」を抑制・安定化させ、高精度に応力分布を解像する特別な平滑化・拡散技術を導入しています。
  • 改良型境界条件(Modified Dynamic Boundary Conditions / mDBC)の拡張:固体壁や地盤境界での粒子配置の乱れを補正するため、境界粒子の外側に仮想的なゴーストノードを投影し、周囲の材料粒子から変数を1次あるいは0次の補間(Shepard補間)で定義する手法を地盤用に拡張・適用しています。
  • マルチボディシステム(MBS)との双方向結合:物理演算エンジン Project Chrono と双方向(Two-way)でカップリングしています。SPH側で構造物(強体)粒子にかかる地盤・流体からの力を集計してChronoに渡し、Chrono側で構造物の並進・回転運動や衝突、拘束条件を計算して再びSPH側へ戻すというプロセスを、時間積分ステップごとに高同期でループさせます。
  •  2相流(水–土骨格)フレームワークの統合:間隙水圧(pw)と土骨格の変位(u)を連成して解く水・土砂の2相流モデルが組み込まれており、雨水の浸透による有効応力の低下や地盤内の圧密現象(Terzaghiの圧密など)を機械的・水理学的境界条件のもとで同時に扱えるよう拡張されています。

3. 結果(Results)

  • 圧倒的な計算パフォーマンス(GPU加速):CUDA並列化スキームの最適化により、単一のGPU(Single GPU)であっても最大で1億粒子オーダーという極めて高解像度かつ大規模な地盤変形計算を実行できることが実証されました。
  • ベンチマークによる検証:1次元圧密試験(Terzaghi理論解との比較)や、球体サンプルの水理負荷試験(Cryer効果の再現)において、計算された間隙水圧の推移や地盤の沈下量が理論値・参照解と高い精度で一致することが確認されました。
  • 実規模の崩壊シミュレーション:5メートルおよび8メートル規模のデリケートな粘土質斜面(Sensitive-clay slopes)の崩壊実験では、雨水の浸透に伴って斜面が段階的に後退崩壊(Retrogressive failure)を起こし、最終的に大流動へと至る複雑なシーケンスを破綻なく再現することに成功しました。

4. 考察(Discussion)

  • 本手法の優位性:GeoDualSPHysicsは、これまで「大規模すぎて計算が終わらない」、あるいは「メッシュの歪みで途中で止まる」とされていた、「豪雨・浸透 → 斜面崩壊 → 構造物への衝突・破壊」という一連の連鎖的災害(Cascading disasters)を一つのシステム内で一気通貫かつ超高速にシミュレートできる点が極めて強力であると評価されています。MBS(Project Chrono)との結合は、防護壁や杭などの防災構造物の耐力評価において非常に実用的です。
  • 今後の課題と限界:一方で、2相流シミュレーションにおける「初期応力状態(Prestress)の設定」や、非常に透水係数が低いケースにおける「間隙水圧の収束性と計算安定性の担保」にはまだ改善の余地が残されています。また、人工的な数値粘性(Artificial Viscosity)や応力拡散(Stress Diffusion)などのパラメータ設定に一部依存する箇所があり、これらを物理的な土質試験データとどのように厳密にキャリブレーション(校正)していくかが、今後の予測精度向上の鍵として議論されています。

土水2層のSPHです。DualSPHysicsベースです。
私が作成したものに比べ、圧密、不飽和透水係数の実装で、こちらの方が進んでいます。良いですね。弾塑性とビンガム流体の切り替えは、実務で使うには調整が必要そう。私も個々は割り切りましたが、ひと段落付けば他の研究成果も取り入れるべきでしょうね。
バグを見つけて修正する速度がAIの進化で劇的に速くなりました。また、機能の追加も同様です。このようにAIを利用してコードを作っていくと、結局同じようなコードができてしまいそうです。不飽和や圧密も容易に実装できると思います。ま、過去の手法の実装や論文に出ている内容は容易に再現できるようになったので、その先を考えられるために時間をとれるようになったと考えるべきなのでしょう。

AI と英語学習 その2

英語でのポスター発表を終えました。

技術面では、オーストリアやノルウェーの観測事例やハザードに対する考え方を学ぶことができました。日本も負けずについていかねば、と改めて感じました。

また、「AIで作ったな」とわかるポスターやスライドが増えていたのも印象的でした。使い勝手が良いので多くの人が利用しているのでしょうが、発表から個性が失われていくさみしさも。ま、私も助長してしまっているのですが。

一方、英語の面で痛感したのは、自身の雑談力の低さです。いくつかの国の方々と話をして、日本のことを説明できない。食事の席では出汁の説明ができず、日本の文化についてもうまく話せませんでした。練習していた英語のアウトプットといえば、技術英語ばかりだったからでしょう。純粋な英会話の練習も必要でした。

スマホの音声モードでAIを使えば会話になります。が、これだけでは上達しないと思います。まずは Input からでしょう。
雑談でよく聞かれる話題について、英語で短く話せる説明をAIに用意してもらうのが効率的かもしれません。食事、温泉、お箸のマナー、季節の行事、地方の話、そして自分の趣味。それらを自分の言葉で言えるまで練習しましょうか。

技術以外のことを考えられるようになったのは、余裕が出てきたのか、それとも歳のせいか。本業の成果を年々厳しく求められていますので、次の機会を得られるかはわかりません。が、プライベートとして気長にいきましょう。

2026年9月21日月曜日

AI と英語学習

AI の進化により、国外の多くの論文を短時間で要約、把握できるようになりました。

とても便利なのですが、苦手な英語を読む必要がなくなったので、より一層、苦手になりました。計算機を使えば計算能力が落ちると考えた昔の人みたいですが、実感しているので笑ってもいられません。
どうにかしないといけない、と思い AIに頼ることにしました。AI を使う障害は AI を利用して克服しよう、という発想です。
まず作ったのは、Readerです。PDFを読み込み、1文づつEnterキーで読み進めます。同時に音声読み上げを実施してもらうので誤って覚えていた発音のチェックにもなりますし、和約、構文、単語の意味を表示しているので、わからない文、単語、イディオムが出てきたらそれらを見る、聴く、という使い方ができます。和訳等のデフォルト表示は無料枠を利用しているため回数に制限があったり精度は低かったりするのですが、都度調べるよりは思考が途切れないので助かります。特に音声は飽きなくて良いですね。

難しい文章が出てきた場合は有料枠の AI 翻訳&解説で詳しく調べます。ま、回数制限を回避するためにも基本はこちらをデフォにする使い方になるのかな。DeepL や ChatGPT の API も使えるようにしてもらいましたが、そちらは試していません。バグも出てきていますので、使いながら改善していきましょう。


もう一つ。近々、英語でのポスター発表を控えています。AI にスキルを使って作らせたポスターと原稿を与え、二十数時間の発表練習教材を作ってもらいました。ずいぶん前に書いた論文なので中身を忘れかけていましたが、10時間ほどこなす内に徐々に思い出してきました。
作ってもらったのは、主要な単語、内容のまとめ、挨拶、説明(2種)、質疑応答です。後半は Claude のチャットに1時間で実施する内容を毎回渡して、音声モードで会話する形にしました。専門的な内容でも練習ができるのが良い点。反応に一呼吸入る場合があるのと、時々暴走するのが悪い点、でしょうか。もう少し AI にフィードバックしながら練習内容を詰める必要がありそうですが、ひとまず教材はできました。
ま、20時間程度で話せるようになるわけではないのですが、今までの学習よりは効率的だとは思います。


考え方なんでしょうね。
AI が何かの障害になるのであれば、別の有利になるような使い方を考える。そういった発想も今後は頻繁に必要になってくるのでしょう。

文献:Causality-based Feature Selection: Methods and Evaluations

DOI: 10.1145/3409382
URL: https://doi.org/10.1145/3409382
書誌: Yu, K., Liu, L., & Li, J. (2020). Causality-based Feature Selection: Methods and Evaluations. ACM Computing Surveys, 53(5), 1–36.

AI要約

【背景】

  • 古典的な特徴量選択は、特徴量とクラス変数の相関(共起)に基づく。因果関係は捉えない。相関だけで選ばれた特徴(例: 喫煙者の黄色い指と肺がん)は、予測に効いても機構の説明にならない。データの分布が変わると、予測も崩れる。
  • 因果特徴量、特にMarkov boundary(MB: クラス変数の親・子・配偶者)には理論保証がある。忠実性等の仮定の下で、MBは最大予測力を持つ最小の特徴集合になる。解釈可能で頑健なモデルにつながる。
  • 多数の因果特徴量選択アルゴリズムが提案されてきた。しかし、包括的なレビューは存在しなかった。比較の基盤となるオープンソース実装も存在しなかった。

【手法】

因果ベース特徴量選択の包括的レビューと大規模実証評価である。まず、ベイジアンネットワーク・Markov blanket/Markov boundary・忠実性・因果十分性といった基礎概念を定義する。対象は、クラス変数のMBまたはPC(親子: 直接原因と直接結果)をデータから学習する枠組みである。アルゴリズムは2系統に分類される。制約ベースは、条件付き独立性検定を用いる。この系統は、さらに2型に分かれる。IAMB系は、MBを一括推定する(同時学習型)。MMMB・HITON-MBは、PCと配偶者を分けて探す(分割統治型)。スコアベースは、SLL・S2TMBなどBN構造スコアを用いる。MB内で原因と結果を区別するアルゴリズムも概観される。次に、代表的アルゴリズムを実装したオープンソースパッケージCausalFSを開発した(C言語、GitHubで公開)。評価には2種のデータを用いた。一つは、ベンチマークBNから生成した合成データである。もう一つは、特徴量数57〜100,000の8つの実データセット(spambase、madelon、gisette、arcene、dexter、dorothea、infant、bankruptcy)である。選択された特徴集合で、NB・KNN・SVMの分類器を学習した。予測精度・コンパクト性(選択特徴数)・AUC・計算時間を比較した。

【結果】

  • クラス変数のPC(親子)集合による分類性能は、MBに劣らなかった。PC学習は、MB学習よりはるかに効率的だった。実用上はPC学習が現実的と確認された。
  • 標本サイズが大きい場合、同時MB学習型(IAMB系等)は他のアルゴリズムより有意に高速だった。予測精度も競合した。FBEDが最速だった。
  • スコアベースのMB/PC学習は、MBが大きいと計算コストが法外になった(例: spambaseでMBサイズ30超のSLL・S2TMB)。制約ベースに対する明確な優位を示さなかった。
  • 合成データでは、後退戦略や対称性チェックが偽陽性の除去を通じて精度を大きく改善した。一方、実データでは、偽陽性を許容する方が分類に有利な場合があった。理由は、ノイズと小標本である。
  • クラス不均衡データ(infant、bankruptcy、dorothea)では、全アルゴリズムのAUCが大きく低下した。既存手法は、不均衡データに対処できていない。

【考察】

  • MBの理論保証を持つ因果特徴量選択は、解釈可能で頑健な機械学習への有望な方向とされる。CausalFSが、手法比較と新規開発の共通基盤になる。
  • 忠実性・因果十分性の仮定は、実データでしばしば破れる。欠測・ノイズを含む低品質データ、ストリーミングデータ、不均衡クラスへの対応が未解決課題として残る。
  • 相関ベースの予測変数は、分布シフトの下で予測を誤りうる(黄色い指の例)。本論文は、この予測と因果の乖離を明確に概念整理した。


文献:DirectLiNGAM: A Direct Method for Learning a Linear Non-Gaussian Structural Equation Model

書誌: Shimizu, S., Inazumi, T., Sogawa, Y., et al. (2011). DirectLiNGAM: A Direct Method for Learning a Linear Non-Gaussian Structural Equation Model. JMLR, 12, 1225–1248.

AI要約

【背景】
  • 共分散構造のみを使う従来の構造方程式モデルやベイジアンネットワークの推定は、構造に関する事前知識なしには、因果順序と接続強度を含む完全な構造を識別できない。
  • 非ガウス性を利用するLiNGAM(Shimizu et al., 2006)では、完全な構造が識別可能である。しかし既存のICAベースの推定(ICA-LiNGAM)は反復探索に依存する。初期値や収束基準しだいで局所最適に陥る。有限ステップでの正解収束も保証されない。
  • ICA-LiNGAMの置換ステップは尺度不変でない。そのため、変数のスケールの違いで誤った順序を返しうる。

【手法】
対象は線形非ガウス非巡回モデル(LiNGAM)である。このモデルはx=Bx+eと書ける。外生擾乱eは非ガウスで相互に独立とする。潜在交絡はないと仮定する。提案手法DirectLiNGAMは、このモデルの因果順序を反復探索なしに直接推定する。核は、Darmois-Skitovitchの定理から導いた補題である。この補題によれば、変数x_jが外生であることと、他の各変数をx_jに単回帰した残差のすべてがx_jと独立であることは同値である。そこで、カーネル法による相互情報量の推定(Bach-Jordanのカーネル独立性尺度)を用いる。残差との独立性が最も高い変数を外生と判定し、因果順序の先頭に置く。次に、その変数の影響を最小二乗回帰で残りの全変数から除去する。残差ベクトルにも再びLiNGAMが成り立つ。残差の因果順序は、元の変数の因果順序と一致する。したがって、この操作の繰り返しだけで、変数の個数と同じ固定ステップ数で因果順序全体が決まる。順序の確定後、接続強度行列Bは最小二乗などの通常の回帰で推定する。さらに、「x_iからx_jへ有向道がある/ない」を成分とする事前知識行列を組み込む拡張も与えた。この拡張は、外生候補の絞り込みと残差計算の省略により、精度と計算時間を改善する。評価では、人工データ(変数10〜100、標本500〜2000、疎・密ネットワーク×18種の非ガウス分布)でICA-LiNGAMと比較した。物理(二重振り子)と社会学(地位達成モデル)の実データにも適用した。実データではPC・GESとも比較した。

【結果】
  • モデル仮定が満たされ標本サイズが無限のとき、DirectLiNGAMは変数数と同じ固定ステップ数での正解収束を保証する。初期値や収束基準などのアルゴリズム的パラメータは不要である。
  • 人工データの全条件(疎・密×次元10〜100×n=500〜2000)で、真のBとの誤差(Frobenius距離)はICA-LiNGAMより小さかった(例: 疎・100次元・n=2000で2.35対36.52)。
  • 事前知識(50%)の利用で、精度はさらに向上した。計算時間も短縮された(密・100次元・n=2000で17.24時間→19.80分)。
  • 計算量はO(np^3M^2+p^4M^3)で、ICA-LiNGAMより大きい。ただし収束保証と既知の計算量をもつ。
  • 実データでは、領域知識と整合する結果を得た。二重振り子では、「角速度は角度で決まる」という知識と整合する有意な辺を検出した(ICA-LiNGAMは有意な辺なし、PC/GESは一部の辺の向きを誤る)。社会学データでは、推定された辺のほぼすべてが領域知識と整合した(誤りは息子の学歴→父の学歴の1辺のみ)。

【考察】
  • 非ガウス性を使えば、辺の向きまで識別できる。このLiNGAMの理論的利点を、収束保証つきの実用アルゴリズムとして実現した。
  • カーネル独立性尺度の計算が、実行時間のボトルネックである。仮定(線形性・非ガウス性・潜在交絡なし)の破れの検定、計算効率の改善、潜在交絡がある場合への拡張が、今後の課題と明記されている。社会学データでの1辺の誤りは、仮定の破れの兆候と解釈されている。
  • 事前知識行列により、「有向道の有無」を制約として与えられる。部分的な領域知識の注入で、精度と速度の両方が改善する。

文献:Metalearners for estimating heterogeneous treatment effects using machine learning

DOI: 10.1073/pnas.1804597116
URL: https://doi.org/10.1073/pnas.1804597116
書誌: Künzel, S. R., Sekhon, J. S., Bickel, P. J., & Yu, B. (2019). Metalearners for estimating heterogeneous treatment effects using machine learning. PNAS, 116(10), 4156–4165.

AI要約

【背景】

  • 平均処置効果(ATE)を超えて、属性で条件付けた効果=条件付き平均処置効果(CATE)を推定したい需要が高まっていた。しかし、教師あり学習アルゴリズムはCATEを直接推定するようには設計されていない。
  • 既存の代表的なメタ学習器には、T-learnerとS-learnerがある。前者は、処置群と対照群を別々に回帰して差をとる。後者は、処置の有無を通常の特徴として1本のモデルで回帰する。どちらも、群サイズが不均衡な場合やCATEが応答関数より単純な場合に非効率になる。
  • 同一の対象では、「処置あり」「処置なし」の両方の結果を観測できない。そのため、個別の処置効果は直接には識別できない(因果推論の根本問題)。

【手法】

任意の教師あり学習アルゴリズムをCATE推定器へ変換する「メタ学習器」の枠組みを整理した。その上で、新手法X-learnerを提案した。X-learnerは3段階からなる。第1段では、対照群・処置群それぞれの応答関数(共変量から結果を予測する関数)を任意の基学習器で推定する。第2段では、「代入処置効果」を作る。処置群では、観測結果から対照側モデルの予測値を引く。対照群では、処置側モデルの予測値から観測結果を引く。この代入処置効果を目的変数として、各群でCATEを回帰する。第3段では、2つのCATE推定を重み関数(傾向スコアの推定値を推奨)で加重平均する。この構成の利点は2つある。一方の群の情報を、他方の群の推定に活かせる。CATEの線形性・疎性・平滑性といった構造を、第2段の回帰で直接利用できる。理論面では、応答関数の収束レートとCATEの収束レートを分けて解析した。実装では、ランダムフォレスト(honest RF)とBART(ベイズ的な回帰木のアンサンブル)を基学習器に用いた。広範なシミュレーションに加え、政治学の2つのフィールド実験(社会的圧力による投票率向上実験・戸別訪問による偏見低減実験)を再解析した。実装はソフトウェアパッケージhteとして公開されている。CATEのブートストラップ信頼区間も備える。

【結果】

  • シミュレーションでは、どのメタ学習器も一様には最良でなかった。ただし、X-learnerは最悪になることが無かった。群サイズが不均衡な場合や、CATEに線形・疎などの正則性がある場合に特に良好だった。
  • 理論では、CATEが線形で対照群が十分大きいとき、X-learnerが処置群サイズnについてパラメトリックレートn^(-1)を達成することを証明した(Theorem 2)。CATEに正則性を仮定しない場合は、T-learnerと同じミニマックス最適レートを達成する。
  • 投票率向上実験(対照191,243人・処置38,218人、ATE 0.081)の再解析では、過去の投票が1回のみの層に有意な負の反動(backlash)を検出した。効果は投票履歴3回の層で最大だった。この異質性は標的化に有用である。
  • 偏見低減実験(最終標本n=501、ATE 0.22)では、原著者が線形モデルで見つけられなかった効果の異質性をX-RFが検出した。S-learnerは推定を0へ縮小し、T-learnerは散らばりが最大で、X-learnerは中間だった。
  • 実データに基づくシミュレーションでは、小標本でX-RFのRMSEがS-RF・T-RFより小さかった。

【考察】

  • X-learnerの本質的な強みは、「CATEが応答関数より滑らかで単純」という実データで頻繁に見られる状況を利用できる点にある。逆に、CATEがほぼ0ならS-learnerが有利になる。応答関数が群間で全く異なるなら、T-learnerが有利になる。この使い分けの指針を与えた。
  • 基学習器の選択も精度を大きく左右する。小データ・大域的構造にはBART、大データ・局所的構造にはRFが向く。因果フォレストの性能は、RFを基学習器とするメタ学習器と同程度だった。
  • CATEの信頼区間はブートストラップに依拠する。形式的な保証のある区間の構成は、今後の課題と明記されている。


文献:Spatial Causal Inference in the Presence of Unmeasured Confounding and Interference

DOI: 10.48550/arXiv.2303.08218
URL: https://arxiv.org/abs/2303.08218
書誌: Papadogeorgou, G., & Samanta, S. (2023). Spatial Causal Inference in the Presence of Unmeasured Confounding and Interference. arXiv:2303.08218.

AI要約

【背景】

  • 空間データの因果推論では、空間的干渉と未計測の空間交絡が絡み合う。一方だけを調べると、他方をそれと誤認する恐れがある。誤認は誤った結論につながる。両者を同時に扱う方法論は乏しかった。
  • 曝露(効果を知りたい原因側の変数)自体が、共通の空間トレンドという固有の空間依存を持つ場合がある。この依存は、交絡が全く無い場合でさえ、独立観測を前提とする標準的な推定を無効にし得る。
  • 動機となる応用は、米国の発電所からのSO2排出が高齢者の心血管死亡率へ与える効果の推定である(2,440郡。局所と近傍50km圏の曝露)。

【手法】

空間データ用の因果ダイアグラム(spatial causal graphs)を導入した。まずペア構造の空間データで、局所効果と干渉効果をポテンシャルアウトカムにより定義した。未計測共変量Uと曝露Zの固有の空間依存は、潜在的な共通トレンドとして表した。直接空間交絡・間接空間交絡・干渉とその組合せからなる6つのDAGシナリオを解析した。解析の結果、次を理論的に示した。空間交絡と干渉は互いに化ける。曝露が空間的だと標準解析は無効になり得る。したがって局所と近傍の共変量を同時に調整し、局所効果と干渉効果を同時に推定すべきである。次に、相互接続した空間ネットワークでは、隣接行列と近傍平均曝露により両効果を定義した。推定はベイズ因果推論(欠測ポテンシャルアウトカムの事後分布からの補定)の枠組みで行う。未計測交絡がある場合、アウトカムモデルに空間ランダム効果を入れるだけでは不十分である。曝露モデルの同時指定が必要になる。著者らは(U,Z)にCAR型精度行列を持つ多変量正規の同時分布を仮定した。線形のポテンシャルアウトカムモデルの下、リング状グラフでの識別定理を証明した(Theorem 1)。未計測空間交絡が存在しても、標本数無限大の極限で局所効果・干渉効果を含む全パラメータが識別可能である。評価はシミュレーションとSO2排出の実データ解析(MCMC・B-スプラインの曝露反応曲線・感度分析)で行った。

【結果】

  • シミュレーション(6シナリオ×3標本数×500データセット)では、未計測空間交絡があると、測定共変量のみを調整したOLSが大きく偏った。95%区間の被覆率は局所効果で0〜0.3%まで崩壊した。
  • 提案法は同じ交絡下でバイアスを大幅に縮小した。全シナリオで被覆率をほぼ名目(約90〜99%)に維持した。
  • 交絡が無いシナリオでは、提案法はほぼ不偏だった。OLSに対する効率の損失もわずかだった。
  • 実データ(SO2→心血管死亡)では、OLSと提案法の曝露反応曲線が互いの信用区間の内側にあった。仮定した形の未計測空間交絡は、OLSを偏らせるほど強くなかった。ただし空間相関パラメータの事後平均は−0.26だった(95%CI −0.41〜−0.10)。曝露と局所的に相関する未計測変数の存在自体は示唆された。
  • 曝露を対数線形で入れる感度分析では、未計測空間交絡の証拠がやや現れた。提案法の効果推定はOLSより正方向にシフトした。

【考察】

  • 空間交絡と干渉は互いに誤認され得る。このため実務では、局所と近傍の共変量調整、および局所効果と干渉効果の同時推定を必ず行うべきである。
  • 未計測空間交絡は完全な欠落変数だけから生じるのではない。測定誤差や関数形の誤指定からも生じ得る。
  • 限界: アウトカム変数自体の固有の空間依存は未対応である(実データでもその兆候がみられた)。識別定理はリング状グラフと線形モデルの設定で証明されている。一般の空間構造は今後の課題である。近傍曝露の地理的近接による定義は、大気輸送の実態と乖離し得る。



文献:The Central Role of the Propensity Score in Observational Studies for Causal Effects

DOI: 10.1093/biomet/70.1.41
URL: https://doi.org/10.1093/biomet/70.1.41
書誌: Rosenbaum, P. R., & Rubin, D. B. (1983). The Central Role of the Propensity Score in Observational Studies for Causal Effects. Biometrika, 70(1), 41–55.

AI要約

【背景】

  • 無作為化のない観察研究では、処置群と対照群で共変量の分布が系統的に異なる。そのため結果の単純比較は交絡バイアスを受ける。
  • 多数の共変量で直接に層別やマッチングを行うと、セルは指数的に細分化される。比較相手のいないセルが続出し、調整は実行不能になる。
  • 多次元の共変量を、交絡調整に十分な情報を保ったまま1次元に縮約するスコアの理論的基盤が求められていた。

【手法】

本論文は理論論文である。潜在結果枠組みの下で議論を進める。この枠組みでは、各個体の処置下と非処置下の結果のうち一方しか観測できない。まず傾向スコアe(x)を定義する。これは、共変量xを与えたときに処置を受ける条件付き確率である。より一般に「バランシングスコアb(x)」を定義する。これは、b(x)で条件付けると共変量xと処置指示zが独立になる関数である。次に両者の関係を定理群として整理する。さらに、因果効果の識別条件として「強く無視できる処置割当(strongly ignorable treatment assignment)」を導入する。この条件は2つの要求からなる。第1に、xを与えたとき潜在結果の組(r1, r0)と処置zが条件付き独立である。第2に、すべてのxで0<e(x)<1が成り立つ。後半の条件は、のちにoverlap(共通サポート)条件と呼ばれる。この条件の下で、同じスコア値をもつ処置群と対照群の結果の差の期待値は、その点での平均処置効果に一致する。応用として、傾向スコアによるペアマッチング・層別(サブクラス化)・共分散調整の3手法を統一的に正当化する。その際、単一共変量の5層別に関するCochran(1968)の結果を、多次元共変量の場合へ拡張する形で接続する。実務では傾向スコアは未知である。そのため、ロジスティック回帰などによる推定スコアの利用も論じている。

【結果】

  • 傾向スコアe(x)はバランシングスコアである。すなわち、e(x)で条件付けると共変量xと処置zは独立になる。
  • b(x)がバランシングスコアであるための必要十分条件は、e(x)がb(x)の関数になることである。すなわち傾向スコアは「最も粗い」バランシングスコアである。これにより、多次元共変量の調整を1次元に縮約できる。
  • 処置割当がxを与えて強く無視できるならば、任意のバランシングスコアを与えても強く無視できる。このとき、同一スコア値での処置群と対照群の結果の差の期待値は、その点での平均処置効果に一致する。
  • この理論により、傾向スコアによるペアマッチング・層別・共分散調整が、因果効果の不偏推定として統一的に正当化される。
  • 傾向スコアの5層別は、測定された共変量による偏りの約9割を除去すると期待できる。

【考察】

  • 傾向スコアが釣り合わせるのは、測定された共変量の分布のみである。未測定の交絡には対処できない。
  • 割当の無視可能性は検証不能な仮定である。感度分析などによる補完が必要になる。
  • 0<e(x)<1(overlap)は識別の必須要件である。処置確率が0や1に近い領域では、「共変量が同じで処置だけが違う」比較相手が実質的に存在しない。その領域の効果推定は、モデルによる外挿にならざるを得ない。



文献:Principles of confounder selection

DOI: 10.1007/s10654-019-00494-6
URL: https://doi.org/10.1007/s10654-019-00494-6
書誌: VanderWeele, T. J. (2019). Principles of confounder selection. European Journal of Epidemiology, 34(3), 211–219.

AI要約

【背景】

  • 観察研究の因果推論では、調整すべき共変量集合の選択が決定的である。Pearlのバックドア基準は、全共変量間の因果構造(完全なDAG)の知識を要する。この知識は、実際には得られないことが多い。
  • 実務でよく使われる「処置前の全変数を調整する」基準は、コライダー調整によるM-biasを招きうる。「曝露と結果の共通原因のみを調整する」基準は、本来バイアスを消せる変数の調整漏れを招きうる。
  • 統計的な変数選択(10%変化基準、前進/後退選択など)は、交絡因子と媒介因子を区別できない。実質的知識の代替にはならない。

【手法】

各共変量について「曝露の原因か」「アウトカムの原因か」だけが分かると想定する。この完全なDAGより弱い知識の下で、実用的な交絡選択基準を提案する。まず、反実仮想の記法で交絡制御の条件(Ya⊥A|C)を定義する。次に、2つの通説的基準の失敗例をDAGで示す。処置前基準の下では、コライダーである処置前変数Lの調整がバックドアパスを開く(M-bias)。共通原因基準は、未測定共通原因の代理変数Cの調整を見落とす。そのうえで、disjunctive cause criterion(曝露またはアウトカムの原因である共変量をすべて調整する基準)を導入する。理論的根拠は次の性質である(VanderWeele & Shpitser 2011)。測定共変量のいずれかの部分集合で交絡制御が可能なら、この基準で選んだ集合でも必ず制御が可能になる。さらに、実用上の2つの修正を加える。第一に、操作変数と分かっている変数を除外する。除外の目的は、未測定交絡下でのバイアス増幅(Z-bias)の回避である。第二に、未測定共通原因の代理(proxy)変数を追加する。加えて、共変量の測定タイミングの問題を論じる。BMIと運動の例では、同時測定は交絡と媒介を区別できない。処置より前の時点の測定値の調整が望ましい。最後に、統計的選択法(前進/後退選択、change-in-estimate、高次元傾向スコアHDPS、TMLE等の機械学習)を体系的に評価する。評価の観点は、交絡制御の妥当性・選択後推論・非崩壊性である。

【結果】

  • 「修正disjunctive cause criterion」を定式化した。すなわち、曝露またはアウトカム(または双方)の原因である共変量をすべて調整する。既知の操作変数は除外する。未測定共通原因の代理変数は追加する。
  • disjunctive cause criterionは次の性質を持つ。測定共変量のどこかに交絡制御に十分な部分集合が存在すれば、この基準の選択集合も必ず十分になる。処置前基準と共通原因基準は、この性質を持たない。
  • 操作変数の調整は、未測定交絡の下でバイアスをむしろ増幅しうる(Z-bias)。このため、判明している操作変数は調整から除外すべきと整理した。
  • change-in-estimate基準は、オッズ比・ハザード比の非崩壊性の下で不適切である。また、データ駆動の変数選択後にそのまま当てはめた信頼区間は無効になる。
  • 統計的選択を使う場合も、初期共変量集合そのものが交絡制御に十分であることが前提になる。したがって、実質的知識は省略できないと結論した。

【考察】

  • ・完全なDAGの知識なしでも、理論的裏付けのある調整集合を設計できる。ただし、各共変量が曝露/アウトカムの原因か否かという実質的知識は依然として必要である。
  • ・TMLE等の原理的な機械学習ベースの選択は有望とされる。ただし、小標本での性能は未解明とされる。



文献:Double/Debiased Machine Learning for Treatment and Structural Parameters

DOI: 10.1111/ectj.12097
URL: https://doi.org/10.1111/ectj.12097
書誌: Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., & Robins, J. (2018). Double/Debiased Machine Learning for Treatment and Structural Parameters. The Econometrics Journal, 21(1), C1–C68.

AI要約

【背景】

  • 高次元の局外(ニュイサンス)パラメータの存在下で、低次元の因果パラメータθ0を推定したい。しかし、局外パラメータ空間の複雑さが標本サイズとともに増える現代的設定では、古典的セミパラメトリック理論の前提(Donsker性など)が破綻する。
  • 機械学習で局外関数を推定して素朴に推定式へ代入すると、正則化バイアスと過学習がθ0の推定量に一次で流れ込む。その結果、√N一致性が失われる。
  • lasso・ランダムフォレスト等の任意の機械学習手法を局外推定に使いながら、θ0の妥当な信頼区間を得られる一般理論が必要だった。

【手法】

本論文は、部分線形回帰モデル Y=Dθ0+g0(X)+U、D=m0(X)+V を主導例とする。まず、素朴なプラグイン推定量の誤差を分解し、g0の学習バイアスに由来する項の発散を示す。この問題を除去する要素として、2つが提案される。第一は、Neyman直交スコアの使用である。これは、局外パラメータの推定誤差への感度が一次でゼロになるモーメント条件を指す。具体的には、処置Dから調整変数Xの効果を回帰m0で取り除いた残差と、結果Yからg0の寄与を取り除いた残差の双方を使って回帰する。この「二重の予測」が、double machine learning(DML)の名の由来である。第二は、cross-fitting(交差適合)である。標本を分割し、局外関数を補助標本で学習し、主標本でθ0を推定する。役割を入れ替えて平均すると、過学習による自己バイアスを避けつつ全標本の効率を回復できる。理論条件は、局外推定量の収束レートに関する弱い要請のみである。このため、lasso・ridge・ランダムフォレスト・ブースティング木・深層ニューラルネットやそれらのアンサンブルを利用できる。この一般理論は、部分線形回帰、部分線形IVモデル、非交絡下のATE・ATT、IV設定のLATEの推定に適用された。さらに、3つの実証例(ペンシルベニア再就職ボーナス実験、401(k)の資産効果、Acemoglu et al. (2001)の制度と経済成長のIV)が使用法を例示する。

【結果】

  • DML推定量は、真値のN^(−1/2)近傍に集中する。近似的に不偏で、正規分布に従う。妥当な信頼区間も構成できる。これらの性質が、弱い条件下で証明された。
  • 数値実験では、非直交の素朴な機械学習推定量の分布が真値から大きく偏った。正規近似も破綻した。直交化したDML推定量は、ほぼ不偏だった。正規近似も良く成立した。
  • ボーナス実験の再分析では、失業期間(対数)へのATEがどの機械学習手法でも−0.073〜−0.085となった。いずれも、5%水準で有意だった。
  • 401(k)適格性の純金融資産へのATEは、無調整で$19,559(SE 1413)だった。柔軟な調整により、推定値は大幅に減衰した。適格性を操作変数とする加入のLATEは、全手法で正かつ有意($8,944〜$11,764)だった。
  • lasso・回帰木・フォレスト・ブースティング・ニューラルネット・アンサンブルの推定値は、広く整合した。この結果は、直交化の理論的予測と一致する。

【考察】

  • Neyman直交化とcross-fittingという2つの単純な要素が、任意の高性能な機械学習器を因果パラメータ推定へ安全に組み込む汎用手続きを与える。
  • 局外関数の推定には、収束レート条件が必要である。極端に複雑な局外関数では、推定は依然として失敗しうる。
  • 401(k)の例では、5分割の交差適合が2分割より標準誤差を小さくした。ただし、分割数と精度の一般的関係はないと注意されている。



2026年7月6日月曜日

AI と"考える"技術者

AI を使うことで、コーディングの壁がなくなりました。
これまで土水問わず多様なシミュレーションコードを触ってきましたが、これらは全てAIが組めるようになりました。プログラムの組み方を知らなくても、AI に教えてもらいながらシミュレーションができる時代になってしまいました。この点で、シミュレーションを専門にする技術者のアドバンテージはなくなりつつあります。生き残るには、これまでにないアイデアを出し続けられる人か、AI の誤りを正せる人か。そう考えると、若い方々をどのように育てるべきか悩むところです。効率化のためには AI に頼るべきですが、生き残るためには AI に頼らず時間をかけて考え、失敗を繰り返しながら判断力をつけて頂きたい。いずれにせよ、AI によって一定レベルの技術者の需要が失われるのは確実でしょうから、それを見据えて今何をするか、何をしたいかを考えてもらう必要があるでしょう。

建設コンサルタントが委託されている設計業務も、ある程度 AI ができるレベルになっています。各種基準書は AI にやさしい形で提供され始めましたし、点群、平面図、用地、各種制約条件等を与えて判断ポイントを教えると、どの会社でも 1, 2 年で実用化できる環境は整っています。もしかすると、国が設計スキル群を公開するかもしれません。そうなると、設計者の数が今ほど必要なくなります。より短時間で誰でもできる仕事になるので、設計単価も安くなるでしょう。そもそも、誰でも一定品質で設計できるように基準やマニュアルが整備されてきた背景があるため、画一化、統一化と相性の良いAIこそ、この理念の行きつく答えなのかもしれません。ま、コンサル業界がそこまで見ていたら、失業を伴う方向に本気で取り組むことはないと思いますが。

そう考えると、重要となるのは測量や調査です。人はもちろん、AIであっても地下の状態はわかりません。いかに精度よく、効率的に3Dデータを現場から持って帰られるか?地下の状態を的確に把握することができるか?近い将来、地質踏査やボーリング調査がこれまで以上に重要視されると思います。ま、さらに先はそれらも自動化されるのでしょうが。

そしてもう一つ。現状の AI は国産ではありません。先日の Mythos や Fable 5 をめぐる一件では、そのサプライチェーンリスクを改めて突きつけられました。頼りすぎるのも危ういが、頼らなければ後れを取る。AI を使うことを前提としながらも、それに依存しきらない考える技術者をどう育てるか。これが、短期的には避けて通れない課題なのだと思います。


2026年7月5日日曜日

3次元すべり面推定→安定計算(逆算)→(順算)LSM作成 その2

3次元すべり面推定→安定計算(逆算)→(順算)LSM作成の続きです。https://phreeqc.blogspot.com/2025/01/3lsm.html

この流れを全てWEBアプリ経由で動かせるようにしていたのですが、この週末に中身を整えました。

SLBLでは、崩壊前DEMに崩壊後DEMも加え、推定面が削剥部より浅ければ最大深にアンカーする処理を追加。 C やすべり面深度の統計量書き出しも追加しました(こちらはまだOpenにしていません)。

安定計算 は bug fixのみ、Region3D は引張/圧縮領域の分布図追加、高速化、bug fixを行いました。Win と Linux で Python library の挙動が異なる箇所はAIが見つけてくれました。ありがたい。

14000箇所の崩壊地のすべり面推定→逆算→順算が半日かからない程度。実用的な時間となりました。
もう少し精度面で手を入れたいと思います。また余裕がある時に手を入れましょう。

スマホアプリでH/V

M5Stack Tough の振動計をアレンジして、H/Vを画面表示するプログラムを作成しました。

Tough の画面上のスタートボタンを押して計測を開始し、計測が終わったらHVSRを表示。ピーク周波数を数値で表示するようにしました(ついでにS波速度を仮定し、1/4波長則で換算深度も表示させました)。これならToughを地面に置いてボタンを押すだけで、どなたでも相対的な硬軟、堆積物の深さがわかる、はずです。
理屈はあっているのですが、ADXL355で微動を評価できるのかは半信半疑。先日の文献でも欲しい3Hz以下がイマイチなようでしたので、使えたらラッキー程度の思想です。どこかで試さなくては。

スマホからも Bluetooth 接続で Tough を操作できるようにしました。スマホ単独でも測定可能ですが、文献を見るとスタックが必要とのこと(時間をかけると評価できるとの報告に驚きですが)。

Tough を制御しデータを取り込むことで、長時間のスペクトル計算が可能です。GNSSによる位置情報をデータに付与して保存し、地図上に表示できるようにもしました。



最近のスマホは処理能力が上がっており、数分間の振動データのFFTも瞬時に処理できました。数多くのデータも余裕で保存できます。位置情報をつけてクラウドに保存する基盤も整っています。そこにAIが入ることで、コーディングの壁がなくなりました。もう、アイデア次第で大抵のことはできそうです。

この程度のことは独創的とは言えません。今後はより尖ったアイデアが必要になってきそうです。

2026年7月2日木曜日

Claude

Claude で Fabale5 を再び使えるようになりました。

早速、過去のコード群を修正。あきらめかけていた資料収集を深堀。で、1時間でクレジットを消費。強いけど、高い。

「デフォルトでリモートコントロールを有効にする」を ON にすると、Code に 別PC のプロジェクトが表示されるようになりました。他PCの計算をリモートで制御できるようになっています。Linux版 Desktop が出たので、Windows , Linux へスマホからでも、指示ができるようになりました。知りませんでした。

Claude Science も発表されました。何ができるのでしょうか。


ひとつづつ、試していきましょう。


2026年6月21日日曜日

文献:Use of MEMS Accelerometer ADXL355 in Microtremor Surveys | IEEE Conference Publication

AI要約
1. 背景
地震調査においてMEMSセンサの活用が期待されています。
従来のMEMSセンサは0.1〜2 Hzの低周波領域における感度が専用ジオフォンに及ばないという課題がありました。本研究では、従来より低ノイズ(22.5 µg/√Hz)なADXL355を用い、微動探査、特にH/V法(中村メソッド)への適用可能性を検証しました。

2. 手法:機器特性の評価と差し引き
機器固有の特性(ノイズ)を正確に把握し、測定データから除去するために以下の手順が取られました。
  • 特性評価(自己ノイズの記録): 微動測定に先立ち、振動のない「ゼロ運動レベル」での計測を、道路から離れた建物の地下(コンクリート基礎上)で70時間にわたって実施しました。
  • 多項式近似: 記録されたノイズスペクトルは滑らかな形状を示したため、これを3次多項式で近似し、センサ固有のホワイトノイズ曲線として定義しました。
  • ノイズの差し引き: 実際の現場で計測された微動スペクトルから、解析の第一段階として、この近似された多項式曲線を差し引くことで、機器特性の影響を排除した純粋な振動成分を抽出しました。
3. 結果
微動の特定: フィルタリング後のスペクトルにより、地盤の自然振動数(約6 Hz)と人工的な振動(8 Hz以上)を明確に区別して特定することができました。

4. 考察
MEMSセンサの周波数応答は、低周波領域において感度が物理的に制限されているため、微弱な深部からの共振を正確に捉えきれない可能性があります。本研究では、0.05〜3 Hzの低周波領域については精度が不十分(不確実)であると判断されました。この帯域でもいくつかの弱いピークが観測されましたが、センサの感度限界に近いことから、信頼性に欠けると判断されました。最終的なH/Vスペクトル図(Fig. 9)において、0.05〜3 Hzの範囲は破線で示されており、この領域の正確な評価には専用の3成分ジオフォンを用いた追加調査が必要であると結論付けられます。

ADXL355とマイコンでH/Vを算出する端末を作成し、スマホから操作できるようにはしました。MEMSですのでノイズの大きさから半分あきらめつつ作っていたのですが、意外とそれらしいH/Vが取れるので調べてみることに。引っかかった文献がこれでした。3Hz以下が難しいというのはそこそこ実体験に合うので、ある程度正しい情報なのだろうと思います。

文献:A Novel High-Frequency Landslide Monitoring Device Based on MEMS Sensors and Real-Time Early Warning Method

A Novel High-Frequency Landslide Monitoring Device Based on MEMS Sensors and Real-Time Early Warning Method

AI要約

1. 背景

地滑りは世界的に甚大な被害をもたらす地質災害であり、特に中国のような山岳地帯では、居住地近くで発生する「小規模ながら大きな被害」を引き起こす突発的な地滑りが深刻な脅威となっています。既存の監視技術(GNSS、InSAR、地中傾斜計など)は一定の成果を上げているものの、「高コスト」「設置の複雑さ」「データ処理の遅延」といった課題があります。特に、衛星測位(GNSS)はサンプリングレートが通常1Hzと低く、崩壊直前の急激な加速フェーズ(数秒〜数分)を詳細に捉えることが困難でした。そのため、低コストで設置が容易、かつミリ秒単位の動態を捉えられる高頻度なリアルタイム監視システムの構築が求められていました。

2. 手法

本研究で開発された GeoMAS (Geological MEMS Alert System) は、以下の技術的特徴を備えています。

  • システム構成: 斜面に設置する「監視ユニット」と、居住区に設置する「アラームユニット」の2ユニット構成です。これにより、物理的な距離(最大5kmのLoRa通信)を利用して避難時間を稼ぐ設計となっています。
  • 高頻度サンプリング: 3軸MEMS加速度計を用い、100〜1000Hzの極めて高い頻度でデータを収集します。これにより、目視では不可能なミリ秒単位の微細な振動や前兆を記録できます。
  • エッジコンピューティングによる高度な信号処理: ウェーブレット閾値フィルタリング(sym8基底)とカルマンフィルターを組み合わせ、環境ノイズを抑制しつつ、加速度から速度・変位を算出する際の積分ドリフトを補正します。
  • パラメータ算出: 加速度、変位に加え、重力成分の空間分解による傾斜角も同時に算出します。
  • 早期警戒アルゴリズム:
    運動強度因子 (y): 加速度の大きさと、異常イベントの発生頻度(時間間隔)を組み合わせた独自の指標を算出します。
    変化率 (y ′ ): y の値を4次多項式でフィッティングし、その微分から変化のトレンドを把握します。
    決定マトリックス: y と ∣y’ ∣ を閾値と比較し、4段階(警報、警告、警戒、安全)のレベルを判定します。また、複数デバイスの情報を統合したリスク指数 (F) による広域評価も行います。
  • 省電力・通信設計: 監視ユニットは20mW以下の低消費電力で動作し、ソーラーパネルとバッテリーで自律運用が可能です。通信は4GとLoRa(自組織ネットワーク)のデュアルモードを採用し、緊急時でも安定した警報送信を維持します。

3. 結果

屋外の大型実験サイトにおいて、人工降雨による斜面崩壊シミュレーションを行い、性能を検証しました。

  • 前兆の早期検知: 目視で亀裂が確認される前の段階で、システムは運動強度因子の異常を検知し、「警戒 (Vigilance)」レベルを通知しました。
  • 警報の発令: 背面に引張亀裂が発生したのとほぼ同時(崩壊の9秒前)に、最高レベルの「警報 (Alarm)」を自動発令しました。
  • 避難時間の確保: 前兆信号の検知から最終的な崩壊まで、システム全体で13秒前に警告を行うことに成功しました。実際の運用では、斜面から居住区までの土砂到達時間を加味することで、さらに長い避難ウィンドウが確保されるとされています。

4. 考察

既存技術との相補性: 本システムは、InSARやGNSSが苦手とする「崩壊直前の過渡的な加速」を捉える能力に長けています。一方で、数ヶ月単位の緩やかな変位の追跡にはGNSSが適しており、これらを統合した「空・宇宙・地上」の共同監視ネットワークの構築が理想的です。
 加速度積分による変位誤差の累積(ドリフト)や、環境変化(雨量など)に応じた閾値の動的調整が今後の課題です。
今後は機械学習を用いた閾値の最適化や、風力・振動発電などのさらなるエネルギー回収技術の導入が検討されています。

崩れる直前1分に焦点を当て、開発されたようです。個人的にはそのコンセプトよりもハードの方が気になりました。ソーラーパネルもLoRaも普及している技術ですが、それをコンパクトにうまくまとめられています。私も手元の機器でこの組み合わせを考えていましたので、参考にさせていただきましょう。

2026年6月20日土曜日

文献:Groundwater pH buffering in carbonate aquifers exposed to high-calcium ash backfill


AI要約
1. 背景
廃止された浅い鉱山跡は、地盤沈下を引き起こし、インフラや公共の安全を脅かす重大なリスクとなります。この対策として、石炭や油母頁岩(オイルシェール)の燃焼副産物である高カルシウム灰( CaO 含有率20%以上の灰)をバックフィルとして地下空洞に注入し、地盤を安定化させる手法が取られています。高カルシウム灰は自己硬化性を持つため、ポルトランドセメントの使用量を抑え、コストと二酸化炭素排出量を削減できる利点があります。
しかし、室内試験ではこれらの灰がpH 12.5 を超える強アルカリ性の浸出水を発生させることが確認されており、地下水環境への悪影響が懸念されてきました。
一方で、世界の陸地の約15%を占める炭酸塩岩は、重炭酸塩(HCO₃⁻)と二酸化炭素(CO₂)の系による強い天然の緩衝能力を持っています。本研究は、実際の炭酸塩帯水層において、この緩衝能力が強アルカリ性のバックフィル浸出水をどの程度抑制できるかを検証することを目的としています。

2. 手法(詳細)
本研究では、室内試験、地球化学モデリング、および2年間にわたる長期的な現場実証実験を組み合わせた統合的なアプローチが採用されました。
  • 現場実験の設計: エストニア北東部の、上部被覆層が10m未満と薄い廃止されたオイルシェール鉱山が選定されました。
  • バックフィル注入: 遊離石灰( CaO )を約 30wt% 含むオイルシェール飛灰を使用しました。これを水灰比 0.65〜0.70 で混合してスラリー化し、計74m3を地下の坑道に注入しました。これにより、10mのテスト区間に平均厚さ約1.6mのバックフィル体が形成されました。
  • モニタリング体制: 注入地点から約10mの距離に、上流側(バックグラウンド用:DH1, DH2)と下流側(影響評価用:DH3, DH4)の計4つの観測井を設置しました。
  • 連続測定: pH、水温、水位を自動センサーにより連続的に記録しました。
  • 定期的サンプリング: 定期的に地下水を採取し、主要イオンや重金属(As, Ba, Cd, Cr, Pb, Se)の濃度を分析しました。
  • 地球化学モデリング: ソフトウェア「The Geochemist's Workbench」を用い、バックフィルから発生したアルカリ性孔隙水が帯水層へ移動する際の化学反応を1D反応輸送計算でシミュレートしました。
  • モデルの設定: 浅い帯水層であることを考慮し、大気中のCO₂と平衡状態にある「開放系」を仮定しました。地下水の移動速度は10 m/日、初期孔隙率は10%に設定されています。
  • 室内浸出試験: 規格(EN 12457-3)に基づき、「新鮮な灰」と「28日間養生して硬化したバックフィル」の両方に対して試験が行われました。なお、硬化した材料は試験前に4mm未満の粒子サイズに粉砕されました。規格(DIN EN 12457-3)によれば、高固形分かつ特定の粒子サイズ以下の材料を対象とした、固液比 2 L/kg および 8 L/kg の2段階バッチ試験です。
  • 得られた成分溶出挙動をモデリングの境界条件として使用しました。
3. 結果
  • pHの変化: 室内試験ではpH 12.0〜12.7 の強アルカリ性が確認されましたが、現場の観測井(10m地点)では有意なpHの上昇は見られませんでした。モデリングの結果とも一致し、アルカリ度はバックフィル近傍で急速に減衰することが示されました。
  • 主要イオンの挙動: 最も顕著な変化は硫酸塩(SO₄²⁻)濃度の増加でした。一方で、カルシウム(Ca²⁺)や重炭酸塩(HCO₃⁻)は、後述する方解石の沈殿に伴い減少する傾向がモデルで予測されました。
  • 孔隙率の減少: モデリングにより、アルカリ性浸出水と地下水の反応で方解石(Calcite)が沈殿し、孔隙率が30日間で10%から約7.5%へと低下することが予測されました。
  • 重金属: 重金属類の系統的な濃度上昇は確認されませんでした。これはバックフィル注入によるpH上昇が限定的であったことも要因の一つと考えられます。
4. 考察
本研究により、炭酸塩帯水層には高カルシウム灰由来の強アルカリを効果的に中和する能力があることが実証されました。
  • 緩衝メカニズム: 地下水に溶解しているCO₂から生成される炭酸が、浸出水の水酸化物イオン(OH⁻)を消費し、方解石として沈殿させることでpHを背景値付近まで低下させます。このプロセスは、大気とのガス交換が容易な浅い帯水層で特に有効です。
  • 水理的障壁の形成: 方解石の沈殿による孔隙率の低下は、バックフィル周辺の透水性を下げ、浸出水の拡散を抑制する「水理的障壁」として機能する可能性があります。
  • 真の環境懸念: 炭酸塩帯水層において、高カルシウム灰の使用による主要な課題はpHの上昇ではなく、硫酸塩(SO₄²⁻)の増加であることが明らかになりました。飲用水源として利用される地域では、pHよりも硫酸塩に焦点を当てたモニタリングと対策が不可欠です。

20年くらい前でしょうか、同じような室内実験結果をシミュレーションの境界条件にしたことを思い出しました。国内ですので指針もなく、理屈が誤っていないことだけが拠り所だったのですが、この論文を見て安心しました。このような規格を知っていたら、もう少し説明性も上がったのかもしれません。

生成AIとオジサンのシナジー

後輩君「表面波の分散曲線を得られない」とのこと。

データの一部を見せてもらいました。
まずはコリレーション。そしてCMP重合。で、位相シフト。Claude に指示したところ、何の問題もなく結果が出ました。彼も Claude を使用して試行錯誤していたのにたどり着けていません。

この差は何なのか?と考えた結果、経験か?となりました。
知識はあっても、不慣れであれば指示が大雑把になるかもしれません。求める答えにたどり着くには、AIの提案に対して細部の判断、方向性のずれの早期修正が必要です。手順、見たい周波数帯、測点の選定範囲、間隔など、基礎知識として点で持っていても、線でつながっていないと適切な指示ができないのでしょう。

以前は基礎知識を習得後、自分でプログラムを組むか市販ソフトの使い方を覚える必要がありました。が、今はそこまで必要ありません。順に何をしたいか、何に気を付けないといけないか、誤っていたら修正させるなど、口頭でそのまま AI に渡せる能力があれば十分です。若い人たちに対しては、ソフトの指導は必要なく、考え方の指導がより重要になっています。もしかすると、オジサンに優しい時代が来たのかもしれません。

基礎知識、ノウハウ、応用力がより重要になっています。オジサンの経験を Skill として記録し展開することで、停滞していた技術の伝承もある程度進むのかもしれません。


2026年6月14日日曜日

生成AIのSkillと研究

Fable5 が使用停止になっている Claude。

UBUNTU では Code を CLI で利用しています。
Windows では GUI?。こちらには Chat, Cowork, Code が含まれています。Chat, Cowork は アカウント共通で、別 PC でも(Chatは)スマホでも同じセッションを扱えます。カスタマイズした skill も共通です。

が、Code の Skill は PC 毎。カスタマイズした Skill とは別です。
plugin も含めて整理すると3系統です。

① アカウント共通のSkill(クラウド型)— Chat / Cowork
保存先: Anthropic アカウント(クラウド)
同期範囲: ログインした全端末で共通。Skill は自分のアカウント専用(プライベート)。Enterprise はオーナーが組織全体へ配布可。
追加方法: アプリ内 Customize > Skills →「+」→ Browse skills / Create skill(Create skill 内で ZIP アップロード)
実行環境: コード実行コンテナ(要:Code execution 有効化 / Settings → Capabilities)
対応形式: 自己完結した単体 Skill のみ。ZIP はルートに Skill フォルダごと入れ、中に SKILL.md。

② PC別のSkill(ローカル型)— Code(Claude Code)
保存先: このPCのファイルシステム。-g 指定時は %USERPROFILE%.claude\skills\、-g なしはそのプロジェクトの ..claude\skills
同期範囲: PCのみ(端末ごとに個別インストール。アカウント同期なし)
追加方法: ターミナルで npx skills add ... -a claude-code(-g=全体 / なし=そのフォルダ(プロジェクト)限定)
実行環境: ローカル

③ Plugin
中身: Skill+エージェント+コマンド+フックの束
保存先: PCの %USERPROFILE%.claude\plugins\(ローカル)
同期範囲: PCのみ(アカウント同期は公式に無し)
追加方法: Claude Code 内で /plugin marketplace add owner/repo → /plugin install 名前@マーケット名
実行環境: ローカル(Claude Code)

node が入っていない PC では、Code に頼めば両方ともいれてくれます。
npx skills add https://github.com/vercel-labs/skills --skill find-skills -a claude-code -g
npx skills add https://github.com/anthropics/skills --skill frontend-design -a claude-code -g
npx skills add remotion-dev/skills --skill remotion-best-practices -a claude-code -g

論文作成に使えそうな Plugin, Skill では、以下が引っ掛かりました。
npx skills add dsebastien/ai-skill-scholar --skill  "*" -a claude-code -g
npx skills add dsebastien/ai-skill-arxiv --skill "*" -a claude-code -g

https://github.com/tam07pb915/ethical-academic-writing
これはclaude.aiでの利用が推奨されています。①クラウド型 Chat / Cowork ZIP 用です。
ダウンロード → claude.ai の Customize > Skills → Create skill でアップロード。
中身は Skill なので、②としても可能。
npx skills add https://github.com/tam07pb915/ethical-academic-writing --skill ethical-academic-writing -a claude-code -g

以下はフックの自動実行が将来的に不安だったので入れていません。が、プラグインの入れ方として残しておきます。
https://github.com/Galaxy-Dawn/claude-scholar/blob/main/README.ja-JP.md

インストール
/plugin marketplace add Galaxy-Dawn/claude-scholar
/plugin install claude-scholar@claude-scholar
/reload-plugins

削除
/plugin uninstall claude-scholar@claude-scholar
/plugin marketplace remove claude-scholar


研究にAIを使わないという選択肢はもうありません。有効な使い方、プロセスが提供されていますし、若い方にとっては標準ツールの一つになっているでしょう。
どう使うかが問われる中で、自分の使い方を模索しつつ、同時に対応も求められている──そのような状況です。




2026年6月9日火曜日

FGAM2

地形も入れたFGAM1よりも、雨だけFGAM2の方が汎化性能が良さそうなので、そちらでコードを動かしてみました。

が、得られた降雨応答局面 F(t, X)がおかしい。崩壊直前より33日前の方が効きが良い。実際は逆です。最初に確率を稼いでおいて、段々減じるような非物理的モデルになっています。



論文では「近い日が効いている」とのことなので、逆転しています。でも、得られた確率はあっている。うーん。
ふと、空間CVの一つを見てるからか!と気づきました。全データ使えばどうか?と思い、予測していた能登も学習データ(紀伊+宮城)に混ぜてみました。
結果はOK。論文通り、近い日も効きました。

mgcvのように省けそうなルーチンが含まれて手数が多くなっているのですが、その割に物理モデルとまでは言えない結果。これなら従来の実行雨量の方がシンプルでは?と思い試してみました。
実効雨量を説明変数にしたロジスティック回帰と、FGAM2 を、同じ学習/検証分割で比較。


半減期が短いほど能登で良い(3日 0.846 > … > 60日 0.770 > 単純累積 0.747)。=最近の雨ほど効く。これは能登の実態に合っています。
FGBM2が合うイベントではなかったのでしょう。

2026年6月7日日曜日

文献:On the use of rainfall time series for regional landslide prediction by means of functional regression

On the use of rainfall time series for regional landslide prediction by means of functional regression - ScienceDirect

コード/データ: GitHub - Vitorecacho/Functional-Regression-Landslide: This is a comprehensive R script designed for Landslide Susceptibility Modeling. It compares Functional Generalized Additive Models (FGAM) against standard Generalized Additive Models (GAM). · GitHub— R言語、`refund` / `mgcv` パッケージ使用

時空間という視点は以前私が寝かせた案の一つ。どのように扱っているのかな?と見れば日本のデータ。これは見なければ、ということで早速AIに投げました。

AI要約

 1. 背景
 降雨誘発型地すべりの早期警戒は、50年以上にわたり降雨閾値に依存してきた。降雨の強度・継続時間・累積量がある臨界値を超えると地すべり発生確率が急増するという経験的・統計的アプローチで、直感的で一般市民にも伝わりやすく世界的に成功してきた。
 降雨閾値は地形特性とは独立に推定され、後段で「掛け合わせ」によって統合されることが多い。しかし統計的には、確率の単純な積で同時確率を求められるのは変数が互いに独立な場合のみ。実際には地形は降雨の時空間分布に強く影響する(植生による蒸発散など)。したがって降雨と感受性は独立と仮定できず、両者の事後的統合は統計的に問題がある。
これを受けて、降雨を感受性モデルに直接組み込む時空間モデリング(space-time modeling)が登場。データを「データキューブ」(平面=空間、高さ=時間)として捉え、地すべりの有無を応答変数とする。
これまでの時空間モデルのほぼすべてが、降雨を時間窓ごとの累積値(スカラー値)に集約してしまい、時系列本来の豊かな情報を失っている。一方、連続した降雨信号をそのまま扱う試みでは、スカラー手法を大きく上回る性能が報告された。
深層学習(LSTM・CNN)の課題:①ブラックボックス性、②大量のラベル付きデータが必要(正確な発生時刻を持つインベントリは希少)、③学習された時間特徴に物理的意味づけができない。
関数回帰(functional regression)の利点:降雨時系列上に滑らかな係数関数 β(t) を明示的に与えるため、「どの時点の降雨が最も効くか」を直接解釈でき、早期警戒システム設計に有用。
Moreno et al. (2025) は関数回帰を試したが、地形+降雨を入れた場合、長い時間窓ではスカラー版と関数版の性能差が消えた。著者らは「これは時系列がわずか6日と短すぎたため」と仮説を立て、より長い時間窓で再検証することを本研究の目的とした。

 2. 手法
 2.1 対象データ(日本の3つの地すべりインベントリ)
国土地理院(GSI)の資料を基に、降雨誘発型地すべりの3事例を使用。

| 地域 | 誘発イベント | 時期 | 地すべりポリゴン数 |
| 紀伊半島 | 台風Talas | 2011年9月 | 1,902 |
| 能登半島 | 豪雨 | 2024年9月 | 1,540 |
| 宮城・福島 | 台風Hagibis | 2019年10月 | 847(土石流・表層崩壊のみ) |

2.2 データキューブの構築
マッピング単位: グリッドではなく斜面ユニット(Slope Units, SU)を採用(`r.slopeunits` + QGISのSZ-pluginで最適化)。
総観測数 64,250 SU(Kii 49,306 / Noto 6,056 / Miyagi 8,888)。不安定SUの割合はそれぞれ3%, 10%, 4%。
静的共変量(5つ): 傾斜(Slope)、断面曲率(Profile curvature)、斜面方位を分解した Eastness / Northness、岩相(Lithology:堆積岩・付加体・火成岩・変成岩の4分類)。すべてGSI DEM等から導出。
動的(関数)予測子:CHIRPS-GEFS による日降水量(予報データ、空間解像度5km)。発生日とその33日前まで=計34日の降雨時系列を、先行降雨(準備的)+誘発降雨を含む形で構築。

2.3 GAM と FGAM
GLM → GAM → FGAM という拡張系列
GLM(一般化線形モデル): 予測子と地すべり発生のlog-odds(対数オッズ)の間に線形関係を仮定する。この線形性が制約となる。
GAM(一般化加法モデル, Generalized Additive Model): 応答変数の期待値を、予測子の滑らかな(非線形の)関数の和で表せるようにしてGLMを一般化する。この滑らかな関数にはスプライン(spline)を用いる。加法構造(additive structure)のため、各共変量の寄与を個別に分離・解釈できる(機械学習のブラックボックス性と対照的)。
指数型分布族の応答を扱え、本研究では地すべり有無を二値(Bernoulli分布)して扱いロジットリンク関数 logit(p) = log(p/(1−p)) で線形予測子を [0,1] の確率に変換する。

FGAM(関数一般化加法モデル, Functional GAM): GAMをさらに拡張し、スプラインのような滑らかな関数項に加えて、関数項(functional term)=予測子そのものが時系列(関数)であるものを組み込めるようにしたもの(McLean et al., 2014)。
降雨時系列に対して回帰でき、地すべり過程に内在する非線形な時間ダイナミクスを表現できる。技術的には各時点での寄与を「ヒートマップ」として捉え、s次元×t次元のパラメータ空間になる。関数項は係数関数(無限次元)と降雨データ関数の L²内積で定義される。無限次元のままでは推定不能なので、K個の基底関数(B-spline)に展開して K次元の有限問題に変換する。

Y(s,t) ~ Bernoulli(p_s,t)

logit(p_s,t) = α + β1·eastness_s + β2·northness_s + γ·litho_s
              + f1(slope_s)          ← 1次元スプライン(基底7個)
              + f2(profcurvature_s)  ← 1次元スプライン(基底5個)
              + f3(precipitation_s,t) ← 2次元スプライン(各次元4基底)=関数項

降雨の関数項は、降雨曲線と各基底関数の積を積分してスカラー共変量 Zik を事前計算し、それを滑らか項としてモデルに投入する(**Penalized Functional Regression, PFR**)。

 2.4 推定
REML(制限付き最尤法)で平滑化パラメータを推定(Rの refund + mgcv バックエンド)。GCV(一般化交差検証)は過小平滑化・局所解の問題があるため不採用。
過学習を防ぐため、係数関数 β(t) の「うねり(wiggliness、2階微分の2乗積分)」に基づくペナルティを課す(平滑化パラメータ λ で適合度と滑らかさのトレードオフを調整)。
基底数 K は有効自由度(EDF)が K より十分小さいことを基準に選択。

 2.5 ベンチマーク(4モデルの比較設計)

| モデル | 降雨の扱い | 地形共変量 |
| FGAM1 | 関数(時系列) | あり |
| FGAM2 | 関数(時系列) | なし(降雨のみ) |
| GAM1 | スカラー(累積値) | あり |
| GAM2 | スカラー(累積値) | なし |

検証戦略: 時間窓を1日ずつ増やして34日まで性能変化を観察。
ランダム交差検証(訓練比率5%〜90%、各100回ブートストラップ)。
空間交差検証(leave-one-inventory-out):2地域で訓練・残り1地域で検証(転移性能を評価)。
評価指標:AUC(カットオフ非依存)、混同行列・TP/TN率(Youden指数でカットオフ設定)、Precision-Recall。

3. 結果
3.1 共変量の効果(FGAM1)
傾斜: ほぼシグモイド状。〜10°は低確率、10〜35°でほぼ線形に確率増加、その後変曲(最大傾斜で確率約0.3)。表層崩壊・土石流という地すべりタイプと整合(緩斜面は安定、急峻すぎると土層が乗らない)。
断面曲率: ほぼ平坦地形でピーク(確率〜0.55)、凸型・凹型では急減。
岩相: 付加体のみが有意に不安定と関連。
方位: 東向き・南向き斜面で不安定性が高まる。
降雨の関数効果: 33日前(day −33)はどんな強度でも発生確率ほぼゼロ。発生日に近づくと確率が急増し、ピークは発生の1〜2日前。3事例とも約3日続いた暴風雨で誘発されたことと完全に整合。

3.2 時間窓に対する性能
GAM2が一貫して最低性能(AUC開始点〜0.73)。
GAM1は地形予測子のおかげで〜0.80から開始。
関数モデル(FGAM1/FGAM2)はスカラー版を上回り、かつ挙動が安定(滑らかに性能向上、特にday −16以降)。GAMは時間窓に対し不規則変動。
GAM1とGAM2の性能差は大きいが、FGAM1とFGAM2の差は小さい。

3.3 データ量に対する安定性
訓練データ割合を変えたときのモデル安定性ベンチマーク。AUCは補集合(検証サブセット)で計算し、手順を100回ブートストラップして箱ひげ図を作成。
FGAMは一貫してGAMを上回り、少ない訓練データで早く性能の漸近線に到達。ランダム交差検証でFGAMはAUC > 0.81(acceptable〜good)を安定維持、GAMは0.81未満。
ランダムCVは過大評価しがち(Anderssen et al., 2006)。隣接SUは地形が類似し降雨も5km格子を共有するため、訓練と検証に近接SUが混在して空間的自己相関で結果が楽観的。 FGAM2は地すべり検出(TP/P, 平均83.8%)でFGAM1(80.7%)よりわずかに高いが、安定斜面の識別(TN/N)では76.2% vs 69.9%でFGAM1が優位=降雨のみのFGAM2はFN(見逃し)を多く出す。

3.4 空間交差検証(leave-one-inventory-out)
| 訓練 | テスト | モデル | AUC | 評価 |
| Miyagi & Noto | Kii | FGAM1 | 0.749 | Moderate |
| Kii & Miyagi | Noto | FGAM1 | 0.681 | Acceptable |
| Noto & Kii | Miyagi | FGAM1 | 0.613 | Acceptable |
| Miyagi & Noto | Kii | FGAM2 | 0.746 | Moderate |
| Kii & Miyagi | Noto | FGAM2 | 0.757 | Moderate |
| Noto & Kii | Miyagi | FGAM2 | 0.723 | Moderate |

FGAM1のAUCは地域により0.61〜0.75(ランダムCVより低い=地域間の一般化は難しい)。FGAM2は3地域すべてで安定して良好(0.72〜0.76)。地形相互作用を入れたFGAM1は日本の多様な地質で性能が落ちることがある。

4. 考察(Discussion)
4.1 支持する論点(強み)
多変量フレームの統計的正当性:感受性評価と降雨閾値を別々に行うのは数値的に妥当でない。多変量回帰は共変量の相互影響を同時推定できる。
解釈性: 統計モデルは回帰係数や確率を直接見ることでモデル自体の理解が得られる。機械学習のLIME・SHAPは事後的(post-hoc)な近似にすぎず、モデル内部の指標ではない。
コンカービティ(concurvity)検証: 加法FGAMが地形(静的)と降雨(動的)の効果を適切に分離できることを確認(降雨0.61/0.35、傾斜0.40等、いずれも0.8閾値以下)。
降雨閾値の概念を継承: 降雨を発生日から遡って累積する手法は従来の降雨閾値と類似で、過去50年の早期警戒の知見を関数回帰に直接統合できる。
運用上の意義: 誤りの種類が重要。FN(見逃し)は人命に直結、FP(空振り)は信頼低下を招く。
予報降雨を使用: あえて観測でなく予報降雨を使い、前向き(forward-looking)に検証。早期警戒システムの要件を満たす設計。

4.2 反対する論点・限界
時間窓の問題: 性能はday −23付近で漸近(GAMはday −29でピーク)。適切な時間窓の定義が未解決の大きな課題。
地すべりタイプ未分類: 3事例は主に表層崩壊・土石流(≒単一タイプ)。タイプ別早期警戒(浅い/深い地すべりの区別)ができず、降雨時系列の地形学的解釈にも制約。深い・大規模な地すべりはより長い時系列が必要と推測(浸透に時間がかかるため)。
データ要件の厳しさ: 高頻度降雨データと正確な発生時刻を持つインベントリが必須。多くのインベントリは時刻が月・季節単位で曖昧。時間的不確実性は係数関数 β(t) にノイズを注入し、特に発生直前の重要な数日で降雨効果を過小評価する。データ希少地域では従来の静的感受性マップ・経験的閾値の継続使用が依然適切。
空間転移の難しさ: Leave-one-out検証(未学習の別地域+別イベントへの転移)で性能が低下(FGAM1で0.61〜0.75)。日本の地質的多様性により、地形相互作用を含むFGAM1は学習データと条件が異なる地域で性能が落ち、降雨のみのFGAM2の方が転移には安定(0.72〜0.76)。

  検証の整理:①ランダムCV=同一イベント内でのSU内挿(リーク込みで楽観的)、②空間CV=未学習の別地域への転移(0.6〜0.76)。
社会・運用的側面: 降雨閾値が今も主流なのはシンプルさゆえ。「○○mm降ったら警戒」は専門知識不要で誰でも理解できるが、関数回帰の確率はそうではない。多くの人に理解される警報ほど従われやすく、結果的にリスク低減につながる。性能だけが最良の基準とは限らず、行政機関・住民との協働による警報設計が必要。

4.3 結論
連続降雨時系列を用いた関数回帰(FGAM)は、スカラー集約モデルより一貫して優れた地すべり予測を与えることを、日本の3事例・多様な交差検証で実証。総合的にはPrecision-RecallでFGAM1が最良。
研究コミュニティに対し、スカラー降雨集約値の慣習的使用を批判的に見直すよう提言。
地域→他地域の転移は性能が落ち(0.6〜0.76)、全国共通の単一式は難しい。転移設定では降雨のみFGAM2が比較的健闘。検証例はまだ少なく、地形・降雨・地すべりタイプを変えた系統的検証が必要。

今後の課題: ①地域差を説明する変数・階層/転移モデルでの汎化、②タイプ別かつ発生時刻精度の高いインベントリでの検証、③降雨の高解像度化・地域別の時間窓最適化、④発生確率(感受性)に加え規模(intensity)・完全なハザードへの拡張、⑤確率出力を行政・住民が受け入れる警報へ翻訳する社会実装。


雨についてはRのライブラリを利用しながら33日前からの崩壊確率を出されてています。実効雨量や近年の機械学習ではもっと前の雨を使いますが、3つのイベントは台風のような短期指標があればよいので、33日とされたのかもしれません。目的が長期なのであれば、違ったイベントで見たいですね。
ランダムCVはリークの可能性があるというのはよく見かけます。汎化の妥当な見積りは空間CV側でしょう。そのFGAM1,2の結果を見ると、地形・地質を加えた効果が相対的に小さい=地形・地質を交えた汎化が難しい→降雨だけの方(FGAM2)がまだ良い、でしょうか。結論を得るには、もう少し地域を増やす必要があるのでしょう。
各地域は1イベントのみのため、「同一地域の“次の”豪雨を当てる(時間的転移)」は本研究では検証されていません。その場合はもっと性能が上がると期待します。

ちょっと雨量の関数化がややこしいのでコードを見てみましょうか。


2026年6月5日金曜日

SPHのGPU化

SPH(粒子法)シミュレーションコードの GPU 化に取り掛かりました。

まずは、土のみ3Dコードから。テストには 斜面崩壊(粒子数77万弱、200,00step)を使用しました。実行環境は Dtransu と同じです。
GPU: NVIDIA RTX 4000 Ada(20GB、CC 8.9)
CPU: AMD EPYC 9754(128 コア)

いくつかの試行を得て、こちらは GPU 版が CPU 版に勝てない結果となりました。
OpenMP  32 / 64 / 128 スレッド 94.4 / 60.3 / 48.0s 
OpenACC 75〜84s 
CUDA Fortran + Thrust(set_box のみ)  41.7〜55.3s  

計算コードの変数を倍精度で計算していたため、FP64律速かつメモリ帯域律速となり、搭載GPU にとって二重に不利、CPU に有利でした。特に処理時間の約45%を占める make_interaction とkernel_correction は、GPU 側で改善できませんでした。粒子のセル順ソートは GPU では悪化、コアレス化は微妙でした。
  
セル分割処理だけはGPUで効果がありました。この処理は整数演算 が主体で、FP64律速にも帯域律速にも縛られなかったのでしょう。
  • セル分割(set_box):空間を格子(セル)に区切り、「どの粒子がどのセルにいるか」の対応表(名簿)を作る処理、近傍探索の前準備です。GPU 化では内部のソートに Thrust の sort_by_key (CUB の基数ソート) を使い、インデックス(cell_id と粒子番号)だけを並べ替えて名簿を作ります(粒子データ本体は動かさない)。これが GPU で効いた処理です。
  • 粒子のセル順ソート(reorder_particles):粒子データ本体そのものを、セルの順番に並べ替えておく「席替え」(50 step 毎)。近くにある粒子をメモリ上でも近くに置き、キャッシュ効率を上げる狙い。set_box内のソートと違い、倍精度の粒子配列を丸ごと移動するため、GPU では悪化しました (CPU ではキャッシュ効果で改善)。
  • メモリアクセスのコアレス化(pair_layout):粒子ペアのデータを GPU が読みやすい並び順に変える。GPU は隣り合うスレッドが連続したメモリを読むと速い(=コアレス化)ので、その並びに合わせる。
セル分割+ソート+コアレス化は DualSPHysics の講義で教えていただいたのを覚えていました。当時はメモリ上の最適配置程度で大きな影響があるとは思わず、ちょっとしたテクニック程度の認識でした。ほかにも、DualSPHysicsでは dx/dw など過去の値に依存しない変数は FP32、累積する変数は FP64 でドリフトを避ける工夫がなされています。

土のみコードは両者を倍精度で計算しており、一部を単精度にすると過去の結果(降伏、流体化の判定)の完全再現が困難。単精度化は周囲の応力均衡を変化させ、隣接粒子の降伏タイミングを次々とずらしていく連鎖を引き起こします。これは避けたいので、今回は速度を犠牲にして見送り。 

今回はここまで。次は、残したFP32化に手を付けるか、良いGPUに入れ替えるか、でしょうか。A100 / H100 クラスの GPU であれば、CPUを超えるかもしれません。

ソフトウェアサポートのカタチ

様々なソフトウェアサポートを利用しています。

近年、電話サポートが少なくなり、メールやチャットが主流になりました。どこも人手不足、コスト改善を迫られているのでしょう。

質は様々。電話サポートは意思疎通が簡単で、短期に問題を解決できます。ユーザーの限られている国内のソフトウェアメーカーや、サポート専門会社、外資系でも特殊なソフトだから実施できるのかもしれません。

ESRI Japan さんは、昔(10年以上前)はレスポンスが早く質も高かった。が、今は AI と1日1回の壁打ちをしているような感覚に陥ります。メールでの返事は1日1回、なので意思疎通ができるまで時間がかかる、帰ってくる内容の質が悪い、問題が解決しない、しかもお高い。AI に聞く方が速くて確実になってきました。
https://phreeqc.blogspot.com/2024/03/survey123-web.html

大手さんほど AI を利用する方向に進むでしょう。が、ユーザーとしては対話がBEST。
立場が逆であればと我が身に引き比べながら、お客様と対話し問題解決を図らなくてはならないとあらためて感じています。

2026年6月1日月曜日

Dtransu のGPU化 その3

テスト環境:CPUが高性能、GPUは中(FP64は低)性能でアンバランス。
GPU: NVIDIA RTX 4000 Ada(20GB、CC 8.9)
CPU: AMD EPYC 9754(128 コア)

1. 密度流 (v2 sample3.dtr, KAN3=2 密度連成) 

流れ:GPU + 粒子追跡:OpenMP  のハイブリッドとしました。
流れ系 (setelm ~48% + solpcg ~28%の self-time) が構造的に並列で GPU 理想形。ここを GPU カーネル化して大幅短縮です。

粒子追跡 MOVE1 連鎖 (~20%) はバケット探索の IF/GOTO  分岐が多く、GPU だと並列効率が落ちました。これは OpenMP スレッドに流す形が BESTでした。

2. 表面流+移流分散 

 一番速かったのは ompacc で、密度流と同じ戦略(流れ:GPU + 粒子:OMP)。全構成中で最大のを記録しました。
surface_flow.f90 (D8 地表水ルーティング) の self-time が ~2.3% しかないため、にここは GPU 化不要、host 常駐で問題なし。残りは密度流と同じハイブリッドという構成でした。


当初、なぜ FP64 が弱い GPU でも速くなったのか不思議でした。AIに聞くと、「典型的な arithmetic intensity(バイトあたりFLOP数)が低い 処理だったから」+「16tだから」とのこと。この辺りをAIに頼るようではまだまだかな。

  • 1要素読むごとの計算量はわずか。演算器はメモリ待ちで遊んでいるので足を引っ張らない。
  • GPUは数千〜数万スレッドを同時に走らせ、メモリレイテンシをスレッド切替で隠蔽。結果、自分の~360 GB/sをほぼ飽和できた。
  • CPU 側の並列効率が悪く、460 GB/sを使い切れず~16t相当で頭打ち。

純 GPU (acc)  が負けるのは、ハードの性能差が大きなこともありますが、分岐の多い粒子追跡が GPU の thread divergence に弱いのも一因。ここをテコ入れするとGPUの方が速くなるかもしれません。なお、途中でSPH の Box Search 手法を取り入れてみましたが、劇的な効果は見られませんでした。
このような計算をしていると、FP64に強いGPUだとどうなるのか、試してみたくなります。

改変コードの公開は配布元が許可されていません。時代の流れでコーディング問題はほぼ解決したので、いずれ配布元でもGPU対応版やMPI版、表流水連携版を公開されると思います。比較は出るまで待ちましょう。

2026年5月22日金曜日

Dtransu のGPU化 その2

Dtransu Ver.2 が手元にあったので、こちらもGPU化に着手。

テストは付属の密度流です。
もともと、OMP対応済みだったので、GPU化のみでした。が、何度かテストするうちに、再現性に劣ることがわかりました。

128tまで実施しましたが、16tで頭打ち。

v1 でも密度流をかけると、同じように再現性を確保できませんでした。テストするモノです。
密度流で PCG の流れが変わるのと、粒子追跡の OMP 化部分が引っかかっているようでした。最終的には v1 で固定しましたが、OMP+GPU の速度は v2 よりやや低下しました。
最初は v2 GPU+OMP でパラスタし、最後に v1 で提出、でしょうか。


2026年5月18日月曜日

Dtransu のGPU化

Dtransu のGPU化を実施。

AIさんがプロファイルをとってくれるので、非常に効率的でした。

RTX4000を使っていたのでFP64性能はイマイチ。それでも、2.5倍速くなりました。次は計算用途も考慮して選択しましょう。

ずっと優先度低で眠らせていた課題。今年の頭に外した途端、解決しました。残りは不連続体+連続体ですが、SPHのGPU化が先ですね。

優先度中:機械学習のスキル増強
優先度低:流体+個体(不連続体+連続体)+振動
優先度低:Dtransu の MPI/GPU 対応
優先度低:地表流+地下水+移流拡散


2026年5月10日日曜日

文献:High-resolution raindrop counting via instantaneous frequency sensing

High-resolution raindrop counting via instantaneous frequency sensing on hydrophobic elastic membranes | PLOS One

AI要約

1. 背景
雨粒の数やサイズを測定する装置であるディスドロメータは、気象観測において重要ですが、従来の装置には課題がありました。高精度な標準機器(JWDなど)は高価で、交流電源や専門知識を必要とするため、広範囲な設置が困難です。
一方で、安価な圧電素子を用いたセンサーは「振幅閾値方式」を採用しており、一度衝撃を検知すると一定時間計算を停止(ロックアウト)するため、時間解像度が低く、連続して衝突する雨粒を見逃したり、小さな雨粒を過小評価したりするという欠点がありました。
本研究は、市販のマイクと疎水性弾性膜を組み合わせ、低コストながらこれらの限界を克服する新しい検知手法を提案しています。

2. 手法
本手法は、振幅ではなく「瞬時周波数」の変化に着目している点が革新的です。

  • ハードウェア構成: 一般的なPVCパイプの片端に、安価なプラスチック膜(PETやHDPEなど)を張り、ドラムのような構造にしています。この膜に雨粒が衝突した際の音をマイクで集音します。
  • 物理的原理: 雨粒が疎水性弾性膜に衝突すると、膜の固有振動(低周波)とは別に、12–18 kHzの非常に短い高周波過渡振動が発生することを発見しました。この高周波信号は雨粒の衝突に特有のものであり、低周波の背景ノイズと分離しやすい特性があります。
  • Acoustic Feature Model (AFM): 音声データから、信号の局所的なピーク(山)と谷の位置を抽出し、ノイズを除去した上で「局所的な周波数推定値」へと変換します。
  • 機械学習 (1D-CNN): AFMで得られた特徴量を、1次元畳み込みニューラルネットワーク(CNN)に入力します。このネットワークは、風などのノイズと雨粒の衝突を判別し、衝突の確率を算出します。
  • 軽量設計: リソースの限られたマイコンやスマートフォンで動作(TinyML)させることを目的としており、モデルのパラメータ数は1000個未満と極めて軽量です。

3. 結果

  • 高い時間解像度: 従来の手法が30ms程度の解像度だったのに対し、本手法は約10ミリ秒の時間解像度を達成しました。
  • 精度: 評価実験において、予測されたイベントの80%~90%が実際の衝突時刻から10ms以内に収まっており、高い精度で雨粒を数えられることが示されました。
  • ノイズ耐性: 振幅を計算に使用しないため、低周波の振幅ノイズ(風など)に対して頑健であり、高い感度を維持したまま背景ノイズの影響を抑えることに成功しました。
  • モデルの有効性: DensenetやResidualブロックを用いたCNNモデルは、基本的なモデルよりもノイズ抑制能力が高いことが確認されました。

4. 考察
本研究の最大の特徴は、材料の物理的な減衰に頼るのではなく、信号処理と機械学習によって時間解像度を向上させた点にあります。

  • 利点と応用: 非常に低コストな材料で構築できるため、教育現場や発展途上国での気象観測ネットワーク、短期間のキャンペーン調査などに適しています。
  • 現在の課題: 強風による高周波ノイズが誤検知の原因となることがあり、さらなるノイズ除去技術や風防の検討が必要です。また、一度の衝突を二重に数えてしまうアルゴリズム上の課題も指摘されています。
  • 今後の展望: 今回は雨粒の「カウント」に焦点を当てていますが、将来的には個々の雨粒の運動エネルギーの推定や、平均的な粒径分布に基づいた降雨強度の算出へと拡張することが期待されています。
  • 未解明の物理: なぜ軟らかな疎水性膜への衝突で高周波振動が発生するのかという微視的な物理メカニズムは完全には解明されておらず、さらなる研究の余地を残しています。

塩ビパイプに膜を張って雨音を計測するだけという、これ以上ない実用的なアイデアです。
パイプのサイズや膜の張力は低周波の共振には影響しますが、検知に用いる高周波過渡振動(12–18kHz)には影響しないとのこと。また、 PETやHDPEといった異なる材質でも同様の信号が観察されています。身近な材料かつDIYによる物理的な製作誤差に対しても極めて高い堅牢性を持っている点が素晴らしい。
機械学習モデルも最初からエッジAI向きに作られています。今後の雨量推定がどうなるか、期待を込めて待ちましょう。


コード整備 その3

おまけです。

 ・振動計
200Hzサンプリング以上にするとデータ保存時にバグっていたのですが、AIに見せると問題点を指摘してくれました。浅いコピーをやめて深いコピー&数を増やすことで問題を回避。500Hzサンプリングでも問題なく保存できるようになりました。

・3次元安定解析
これもAIがバグを発見してくれました。学習量の多いであろうPythonで書いていたためか、割と早く修正が完了しました。正しい部分も怪しそうなところは指摘することもありましたが、そこに注意しておけば、コードの照査にかなり有効です。

・RegionGrow3D
MATLABコードをPythonに変換。近い将来できるようになると思っていましたが、もう既に可能だったとは。
1回目の変換で結果はピクセルベースで86%合致。もう少し詰める必要はありそうですが、なかなか賢い。

全体を通してみると以下のような感触。
◎既存コードを他の言語に変換したり修正したりする。
◎コードの照査
◎WEBサイトの解析
〇新たな実装
〇機械学習の自動化

まだまだ性能が上がっていくのでしょう。楽しみです。


2026年5月9日土曜日

AIがAIを鍛える

GWの間、コア写真から柱状図を作成する機械学習モデルの訓練を、ほぼAIに任せきりにしていました。先日記載したとおり、ほとんど人手を加えず、AIの提案を承認・却下するだけ。最後は、気が向いたときにRQD等のスコアの伸びを確認する程度でした。手をかけなくてもスコアが徐々に伸びていくのは、精神的に楽。連休も終わりましたので、ひとまずここで一区切りとします。

昨年まではモデルをどう構築するか、比較対象としてどれだけのモデルを用意するか、最適化をどう進めるかなどを考えながら自ら組んでいました。これらの作業も、今年からはほぼAIに任せられそうです。相手は人ではありませんから、提案方針に何度ダメ出しをしても気を使う必要はありません。納得がいくまで試行錯誤させることができます。文献収集力や実装の速さは人を上回り、初歩的なミスに目を瞑れば、24時間無休で働いてくれる部下ができたようなものです。

一方で、誰でも機械学習モデルを組める状況になったことで、AI活用そのものを差別化要因とするのは難しくなるでしょう。土木分野でも、画像や動画を使った異常検出や点検、数値データを使った予測といった難しいとまでは言えない作業の自動化があります。これらに投資してきた企業が案件を獲得してきた側面もあります。しかしAIの進化により、誰もが短期間で同等のモデルを構築できる時代になってしまったため、こうした先行投資・利益回収型のビジネスモデルは終焉を迎えつつあるように感じます。短期的には実績がある分、有利な状況が続くとは思いますが。

これからは、より一層アイデアと新規性が求められる時代になるでしょう。発想の豊かな技術者こそが重宝されるはずです。そのためには、基礎力を疎かにしないことが欠かせません。若い方々には、便利なAIを活用しつつ、時間をかけてこそ身につく基礎力を養ってほしいと思います。

2026年5月6日水曜日

文献:A multi-layer SPH method for generic water–soil dynamic coupling problems

A multi-layer SPH method for generic water–soil dynamic coupling problems. Part I: Revisit, theory, and validation - ScienceDirect

AI要約

1. 背景
水と土砂の動的な相互作用は、地滑りや土石流といった自然現象から、防波堤やダムの浸透流といった工学的問題まで、幅広く存在します。これらの現象を数値シミュレーションで扱うには、土砂の大きな変形、水面の自由表面流、そして両相間の複雑な相互作用を同時に扱う必要があります。従来の格子法(有限要素法など)では、大きな変形によるメッシュの歪みや自由表面の追跡が困難でした。一方、SPH法(粒子法)はメッシュフリーであるため、これらの問題に適しています。しかし、従来のマルチレイヤーSPHモデルの多くは、空隙率(porosity)の時間・空間的な変化を正確に考慮していなかったり、流体の体積保存に欠陥があったりするなどの課題がありました。
本研究は、混合体理論に基づき、これらの課題を克服した包括的な3次元数値フレームワークを提案することを目的としています。

2. 手法
本研究で提案された手法の核心は、混合体理論(Mixture Theory)マルチレイヤー(多層)粒子で表現するSPHの実装にあります。
数学的定式化(ADFとIDF): 以下の2つのモデルが提案されました。

  • 仮密度ベース定式化 (ADF): 混合体単位体積あたりの質量(仮密度)を用いる手法。
  • 真密度ベース定式化 (IDF): 各相の実際の材料密度(真密度)を用いる手法。

特にIDFは、空隙率の空間的・時間的な変化を厳密に考慮できるよう、本研究で新たに質量保存則が導出されました。

支配方程式と構成則:

  • 水相: 弱圧縮性ニュートン流体としてモデル化され、状態方程式(EOS)を用いて圧力を算出します。
  • 土砂相: 弾塑性構成則(Drucker-Prager降伏基準)を用い、土砂骨格の変形と有効応力を扱います。
  • 相互作用力: 浮力(buoyancy force)と粘性抗力(viscous drag force)を考慮します。抗力には、層流から乱流まで対応可能なErgunの式に基づいた2次形式が採用されています。

数値的な工夫:

  • 粒子体積の補正: 流体粒子が土砂領域に出入りする際、空隙率の変化に応じて粒子の体積を適切に調整することで、流体の体積保存を確実にします。
  • 補正シェパードフィルタ (CSDF): 粒子の一貫性を確保し、自由表面付近のノイズを低減しつつ滑らかな圧力分布を得るための新しいフィルタ手法が導入されました。
  • 境界条件: 壁面における粒子の貫通防止と、自由滑り/ノンスリップ条件を適切に扱うための一般化境界粒子法が実装されました。

3. 結果
提案されたモデルの妥当性は、複数の検証ケースを通じて確認されました。

  • 体積保存の検証: 水が土砂領域に流れ込む単純な崩壊試験において、粒子体積を補正する手法(VA)が正確な水深を再現できることを示しました。一方で、補正を行わない従来の手法では、流体体積が保存されず不正確な結果となりました。
  • U字管内の浸透流: 空間的に空隙率が変化する土砂を通る水の流れをシミュレーションし、解析解と非常によく一致することを確認しました。また、CSDFの導入により、長時間のシミュレーションでも自由表面の乱れが抑えられることが実証されました。
  • 重力下での沈下: 水中に沈められた土砂に重力が加わる動的な過程を扱い、水圧と土砂の有効応力が最終的に理論的な静水圧分布に収束することを確認しました。
  • 水中地滑りの再現: 実験スケールの水中地滑りをシミュレートし、土砂の堆積形状や、地滑りによって誘発される波(impulsive wave)の形状が実験データと良好に一致しました。

4. 考察
本研究の結果から、以下の点が考察されています。

  • ADFとIDFの比較: 両モデルは数学的には等価ですが、SPHの実装においてはIDFの方が数値的に安定し、精度が高いことが判明しました。これは、ADFが空隙率の変化に対して敏感すぎ、圧力計算にノイズが乗りやすいためです。
  • 統一的なフレームワーク: 提案手法は、純粋な水領域、乾燥した土砂領域、およびそれらが混合した飽和領域を、単一の方程式系で統一的に扱える点が強力です。
  • 今後の課題: 3次元シミュレーションは計算負荷が高いため、GPUによる加速が不可欠です。また、水中地滑りの実験で見られた土砂の急激な体積膨張(細粒分の拡散など)を完全に再現するには、単なる構成則を超えた、水と土砂のより高度な界面混合モデルが必要であることも示唆されました。

密度を変更しないと、計算が複雑にならず理解が楽。IDFがADFに比べて安定であることも実感しました。ただ、地中で体積を膨張させるというのは初期配置が面倒。私は粒子分割・統合で対応しました。発表が2022年ですから、改良案は出ているのかもしれません。


文献:3D semantic mapping of surface geological features

3D semantic mapping of surface geological features - ScienceDirect 

AI要約

1. 背景
地学研究(地形学、災害評価、環境モニタリングなど)において、3D空間上のデータに「岩」「断層」などの意味(セマンティクス)を割り当てるセマンティックマッピングは非常に重要です。しかし、以下の課題が大きな障壁となっていました。

  • アノテーションの困難さ: 数百万点におよぶ大規模な3D点群データに、手作業で正確なラベル(注釈)を付けるのは極めて困難で時間がかかります。
  • トレーニングデータの不足: 地質学分野に特化した、機械学習用の汎用的な3Dデータセットが不足しています。
  • 既存モデルの限界: コンピュータビジョン分野の既存3Dセグメンテーションモデルは、構造化されていない屋外の広大な地質データには最適化されていませんでした。
本研究は、これら「学習データの不足」と「大規模データの処理」という2つの大きな課題を解決することを目指しています。

2. 手法:独自の統合アルゴリズム「SegMo3D」
提案手法の核となるのは、既存の大規模視覚モデル(LVM)セグメンテーション・モザイキング・アルゴリズムの統合です。

  • 2Dセグメンテーション(意味の抽出): ドローン画像に対し、学習不要(ゼロショット)で動作するSAM2を適用し、画像内のオブジェクトを自動で切り出します。さらに、Grounding DINOを併用することで、テキスト(例:「rock(岩)」)による検索・抽出を可能にしました。
  • 3D-2D投影(空間の紐付け): SfM(Structure-from-Motion)を用いて、写真から3D点群を生成し、カメラの位置・向きを推定します。この情報を使い、3Dの点を2D画像上に投影することで、どの点群がどの画像ピクセル(およびセグメンテーション)に対応するかを計算します。この際、手前の物体に隠れた点を除去する「デプスフィルタリング」を適用します。
  • セグメンテーション・モザイキング(新規性の核): 異なる角度から撮影された多数の画像による認識結果を、3D空間で一つに繋ぎ合わせます。
  • 投影フィルタリング済み3D IoU: 異なる視点からの認識結果をマージする際、3D点を一度別の画像平面に再投影して重なりを確認する独自の基準を導入し、視点違いによるノイズを排除しました。
  • 確率的投票戦略: 個々の画像の認識結果をそのまま信じるのではなく、複数の画像からの「証拠」を確率的に集計して最終的なラベルを決定します。これにより、一部の画像で認識ミスがあっても、全体として正確な判断が可能になります。
  • 効率化設計: 計算時間が画像数に対して線形に増えるように設計されており、数千万点規模の膨大なデータも現実的な時間で処理できます。

3. 結果
合成データと実世界のUAV(ドローン)データの両方で検証が行われました。

  • 精度の向上: 合成データ(Kubric)を用いたテストでは、既存の最先端手法(SAM3D、SAMPro3D)と比較して、誤検出(false positives)を大幅に減らし、高い平均精度(mAP)を記録しました。
  • 実世界データへの適応: アリゾナ州の岩壁データ(約1,250万点)において、自立した岩の抽出に成功しました。競合手法のSAM3Dが大規模データに対応できず、データを間引く(ダウンサンプリング)必要があったのに対し、SegMo3Dは高解像度のまま処理を完遂しました。
  • ケーススタディ: 抽出した岩の形状から「高さと幅の比率(脆弱性の指標)」を自動計算し、地震発生時に崩落しやすい岩(PBR:不安定な岩)の脆弱性分布マップを作成することに成功しました。

4. 考察

  • 実用性: 専門的な学習データを用意することなく、既存のSfMワークフロー(Agisoft MetashapeやOpenDroneMapなど)に組み込んで、誰でも高精度な3Dセマンティックマップを作成できる道を開きました。
  • 技術的な進歩: 従来の「下から順に繋ぎ合わせる(エラーが蓄積しやすい)」方式から、複数の視点を統合する「確率的決定」方式へと転換したことで、複雑な地質環境での安定性が飛躍的に向上しました。

今後の課題:

  • メモリ消費: 大規模な関連付けデータを保持するために、非常に大きなRAM容量(実験では約40GB)を必要とします。
  • SfM精度への依存: カメラ位置の推定精度が低いと、空間的なズレが生じます。
  • 2Dモデルの限界: 最終的な精度はSAM2などの2D認識能力に依存するため、地質学に特化した2D学習モデルが登場すれば、さらに精度が向上すると期待されます。

UAV画像による点検にも使えそうです。
voting による視点間のセグメンテーション統合という方法は、見たことがありません。素晴らしいですね。

資料では、不安定な岩(PBR: Precariously Balanced Rocks)が転倒するために必要な最小のPGA(PGA overturning)は下式とされています。

PGA  =1.3g・tan(W/H)

g: 重力加速度
W/H: 幅と高さの比(幅を高さで割ったもの)

この関係性は、地震時のハザード分析や落石リスクの評価に直接活用できるとのことですが、どの程度役立つかは未知。スクリーニングに使う程度でしょうか?落石源としての独立性指標も入ると、より実務で使いやすくなるでしょう。


Dtransu + 地表流 +移流分散

Dtransu に地表流を組み込みました。

以前(といっても10年も前になるのですが)、地表の水頭Pを排除しないように改造していました。それを地表流として平面上で散らすだけです。散らす方法はそろっていますので、あとは組むだけでした。それもAIが頑張ってくれるようになったので、重い腰を上げることに。
https://phreeqc.blogspot.com/2020/05/kinematic-wave-diffusion-wave.html

組んでみると、それなりに谷部には水が流れます。が、思ったほど高くなりません。30mメッシュで 1㎞×1㎞ のテストケースを作っていたのですが、もっと谷部を細かく分割しないと高くなりません。また、表層も薄くしないとPが精度よく出ません。つまり、地下水のみの計算よりもメッシュを細かく切る部分が発生するということです。ま、当然なのですが。
計算時間もかかります。地下水の流れは遅いので dt を飛ばせるのですが、地表流はその間にサブステップを回す必要があります。
濃度では、同じPでも地下水からの湧き出し、降雨による希釈を分離する必要がありました。こちらは組んでいる最中に気が付きました。

組んでは見たものの、思ったよりも計算時間がかかるので実用的でないなあと。また寝かせることになりそうです。

2026年5月4日月曜日

地理空間MCP Server

地理空間情報:地理空間MCP Server - MLIT Geospatial MCP Server - (α版) - 国土交通省

国交省さんが公開されたMCPサーバーです。第一印象は「対応が早い」でした。河川砂防技術基準の学習用データセット公開に続き、これです。今後、さらに展開されるのでしょう。

動画では Claude Code を利用されています。これ、かなり使い勝手が良いので私も利用しています。もうChat に戻れません。

私の勤めている会社は Chat 止まりです。MCPサーバーにアクセスするためのツールはホワイトリストに入っていないため利用できません。が、そのうち対応してくれるでしょう。たぶん。
Office365  もMCP 経由で操作できるようになりましたので、早く業務メールも処理できるようにして欲しいのですが。

コード整備 その2

・ボーリングコアから柱状図作成
以前は手抜きで生成AIに丸投げ。そりゃダメよね、という結果でした。
今回は画像認識、評価を正攻法で積み上げる仕様で着手。古典的な画像処理からSAM2, Dino v2, YOLO8など、何を使うか選択させるところから実装、比較、評価、次の方針提案までをAIさんにお任せ。こちらは出てきた結果を見て方針変更の必要性を適時与えるだけ。手がかかりません。
コア、サンプリング箇所くらいは教師ナシでもXMLから判断できるだろうと思いきや、意外と難しい。仕方がないので少ない教師データを与えるとで適切に学習が進みました。教師データ作成のためのセグメント指定アプリもAIさんに作っていただきました。ありがたや。
まだセグメンテーションが完璧ではなく、その後の評価もまだまだですので方針を変える必要ありかもしれません。が、自動で試行を繰り返してくれるのは(時間がかかっても)手がかからないのでうれしい限りです。

・PCからエッジAIへ
PCで作成した巨大モデルをセンサー側に載せようと思い、方針作成をAIさんに指示。実機が手元になかったのテストまではできていませんが、数時間でプログラム作成までたどり着きました。しかも、上記の柱状図作成と並行で進行。部下を2人得たような気分です。

・過去コードのバグフィックス
過去に書いたコードを読ませてバグを調査。それなりに出てきます。もともと、これが目的の一つで有償AIを利用しているのですから、ありがたいと思うべきなのでしょうが。ま、修正もAIさんにお願いするので、それほど手はかからないか。

AIを利用したコーディングには、ある程度の知識やノウハウを有していた方が効率的です。その点、これまでのコーディング経験は役立つだろうと思います。若い方はノウハウなしで作っているので、また違った感性が育つのでしょう。
1年後、どのようにAIが変化して、技術者が対応するようになっているのか、楽しみです。

2026年5月2日土曜日

コード整備

先月、地方への出張が続いたためか、出不精がより落ち着いて家にこもっています。時間のあるうちに、いくつかのコードとマニュアルを整備することにしました。

・機械学習モデルの利用サンプル
お客様のご依頼で、Pythonを触ったことのない方でも使えるようにブラウザ上で判定させるサンプルモデルを作成。Python 環境を Docker file にまとめ、OSに関わらず動くようにしました。

・水文水質データベースのスクレーパー
国土交通省 水文水質データベースから、以下のデータをCSVファイルとして取得できるようにしました。
  • 観測所一覧 雨量・水位流量・ダム堰
  • 時間雨量 指定観測所の任意期間の時間雨量 [mm/h]
  • ダム諸量 指定ダムの任意期間の流域平均雨量・貯水量・流入量・放流量・貯水率
  • 水位・流量 指定観測所の任意期間の水位 [m]・流量 [m³/s]
このスクレーパーは、後輩さんからの依頼です。以前、Python で似たようなスクリプトを作成していましたが、Pythonを知らない方でも簡単に使えるようにC++へ移植・拡張。すべてブラウザのフォーム操作のみで完結するGUIも作成し、コマンドプロンプトでなくとも動かせるようにしました。
リクエスト間隔を2秒に固定しサーバー負荷への配慮しつつ、取得済みデータはローカルにキャッシュして二度取得しないようにしました。どの程度のアクセス頻度がOUTなのか判断できなかったため、公開は控えています。

・SPH 土ー水連成
Drag Force を考慮した2相2層でのコードを整備。もともと書いていたコードに対し、以下の点を変更しました。 
  • IDF(Peng 2022)への移植 — 水相を intrinsic density rho_tilde で記述し、自由水-飽和土界面の数値不整合を解消。ベースライン定式化を置換し、以後の全機能の土台に。
  • 境界条件の刷新 — mDBC(English 2022) へ統一(Shepard ghost + hydrostatic 補正)。
  • 数値安定化レイヤ追加 — Tait γ=7において、Riemann 圧力拡散(Zhang 2017)、δ-SPH 拡散、Bonet-Lok カーネル補正、HHT-α 積分を追加。
最も効いた改編はIDF 移植でした。
上記以外にもAIに多様な文献を集めさせて実装し試しています。便利な時代になったのですが、それら棄却した内容に対する理解が追い付いていません。復習が必要です。

つづきます。

2026年4月27日月曜日

AI コーディング

AIの助けを借りてコーディングし始めたころも「すごい」と思いましたが、今はAI主導でコーディングできるようになっています。あらためて「すごい」。これから何度驚かされるのでしょうか。

1年ほど寝かせていたコードですが、2週間ほどで問題点を見つけ、概ねそれらしく動くようになりました。その間、次々に露呈する問題に対して対策を施した文献を探し、実装し、試行したのはAIさんです。1年前と比べても、文献収集力、試行錯誤の効率が段違いです。仕事で使わない手はありません。

この時代、学生さんから見たら、AIを使ってもプログラムを書けないオジサンは老害に見えるかもしれませんね。経験ベースのオジサンたちが古いやり方を踏襲しようとすると、それだけでハラスメント?になるかもしれません。時代は繰り返す、でしょうか。昔も今も、このような構造は変わらないのでしょう。

さあ、コーディング問題は片付いたも同然です。寝かせていた子たちを起こしていきましょう。

文献:Surface-wave dispersion spectrum inversion method applied to Love and Rayleigh waves

Surface-wave dispersion spectrum inversion method applied to Love and Rayleigh waves recorded by distributed acoustic sensingDAS surface-wave inversion | Geophysics | GeoScienceWorld

AI要約

1. 背景(問題点)

DASの特性と制約: DASは既存の光ファイバー網を利用して高密度な空間サンプリングを低コストで実現できる技術ですが、標準的な直線状ケーブルでは主にケーブル軸方向の動的ひずみを記録します。
波の分離の困難さ: 軸方向の成分のみを記録するという特性上、レイリー波とラブ波を分離することが困難です。
従来手法の限界: 従来の表面波解析(分散曲線解析)では、震源とケーブルを一直線上に配置する「インライン」観測で得られる純粋なレイリー波(またはラブ波)のみを扱うのが一般的でした。
クロスラインの問題: 震源がケーブルの延長線上にない「クロスライン」配置では、レイリー波とラブ波が混在して記録されるため、標準的な手法では正確なモデル化ができませんでした。これを解決するために特殊なケーブルやアレイを用いると、導入コストが大幅に増大するという課題がありました。

2. 手法

本研究では、分散曲線を抽出する代わりに分散スペクトル(周波数-位相速度領域)を直接インバージョンする手法を提案しています。
データ処理(観測スペクトルの構築)
スラントスタック法: 多チャンネルの記録をフーリエ変換で周波数領域に変換し、位相速度を仮定した位相シフトを与えて重ね合わせることで、分散スペクトルを構築します。振幅の変動を抑えるための正規化も行われます。

フォワードモデリング(理論スペクトルの生成)
モード集約法(Modal Summation Method): 推定された地下構造モデル(Vs,Vp,ρ、層の厚さ)に基づき、各受振器位置における放射方向(radial)と横方向(transverse)の変位記録を算出します。
DAS記録への変換: 計算された変位を、震源と受振器の方位角に基づいてケーブル軸方向の変位へと線形結合します。
ひずみとゲージ長の考慮: 得られた変位記録に1次中心差分を適用することで、DASが記録する「動的ひずみ」へと変換します。この差分をとる距離が実際のDASのゲージ長に対応します。
スペクトル化: この理論波形に対しても観測データと同様のスラントスタック法を適用し、理論的な分散スペクトルを生成します。

最適化
遺伝的アルゴリズム(GA): 観測スペクトルと理論スペクトルの差(二乗平均平方根誤差)を最小化するように、遺伝的アルゴリズムを用いて最適なS波速度(Vs)モデルを全域的に探索します。

3. 結果

数値シミュレーション: クロスライン配置(混合波)のデータからでも、インライン配置と同等の精度でVs構造を復元できることが確認されました。
低速度層(LVL)の検出: ラブ波に含まれる低周波成分が追加の制約を与えるため、レイリー波単独(インライン)の解析よりも地下深部や低速度層を正確に復元できることが示されました。
実データへの適用: ネバダ州の地熱地帯(PoroTomoプロジェクト)のアクティブ震源および環境ノイズ(NCF)データに適用した結果、異なるデータセットや観測配置の間で一貫性のある地下構造モデルが得られました。

4. 考察

観測設計の柔軟性: 本手法により、震源を必ずしもケーブルの直線上に配置する必要がなくなり、既存の通信用光ファイバー(ダークファイバー)を活用した調査が容易になります。
データ利用効率の向上: 従来は解析が困難だった混在信号を有効活用できるため、利用可能なデータ量が増大し、より高解像度な地下イメージの構築が可能になります。
パラメータ感度:密度の不確実性が解析結果に与える影響は5%未満と極めて限定的です。
P波速度(Vp) については、未固結堆積物では影響が小さいですが、固結した岩盤では誤差が最大20%に達する場合があるため、必要に応じてVp も同時にインバージョンに含めるべきです。
ロバスト性: スラントスタック法はインコヒーレントなノイズを抑制する効果があり、10〜20%程度のノイズが含まれるデータや、10〜15度程度の地形勾配がある場合でも、浅部構造については安定した結果が得られます。

分散曲線を抽出する必要がなく、また Love, Rayleighを気にしなくて良いので、叩く位置をファイバーから外せます。微動計の上を叩いているようなもので、すぐに飽和しちゃいますからね。
この程度ならすぐに組めそうです。

2026年3月28日土曜日

地上開度 その2

久しぶりに地上開度を求めることになりました。

以前の記録を見ると、15年ぶりです。その間、利用することはありましたが、自分の手で作ることはありませんでした。

ArcGIS pro v3.6 を確認しましたが、未だ作成機能は実装されていません。国外製品に実装されないということは、国内限定の指標なのでしょう。

今回は python でArc 用のツールボックスを作って計算してみました。
で、上記 QGIS v3.40  (アドイン:Processing Saga NextGen Provider) との結果の比較段階で、相違に気づきました。

SAGA の Morphometric Protection Index は仰角(水平からの上向きの角度 [radian])です。地上開度(鉛直からの角度)に変換すると、水平:MPI 0度→開度90度となります。水平から下向きの見下ろし角は無視されているため、痩せ尾根や山頂などでも最大90度、平野や滑走路のような広い水平場所でも90度となります。
一方で、地上開度の論文では水平から下向きまで考慮されており、90度以上の値をとります。

NoData 周辺の取り扱いも注意が必要です。特に池や湖。今回は1探索方向の全セルが NoData の場合は90度に設定し、その後に8方向平均をとりました。その方向の探索をやめるとか、8方向平均をとらずに NoData 扱いにするとか、いくつかの考え方があると思います。どのような処理をされたデータなのか、確認が必要です。

15年前はこれらの点に気づいていませんでした。もらった開度データに関しても、どの範囲をとっていたかまで見ていませんでした。今回のように論文を見て実装するなど、手を動かさないとダメでした。
基本、開度を自動で計算できるツールは存在していませんので、頂くデータに関しては、作成手法とデータの範囲に留意しましょう。

*********************************************
20260329 追記
地上開度を計算するArcGIS  pro 用の Python toolbox (.pyt) を GitHub の Public リポジトリに置きました。

2026年3月23日月曜日

文献:Real-time lithology identification while drilling based on drilling parameters analysis with machine learning

Real-time lithology identification while drilling based on drilling parameters analysis with machine learning | Geomechanics and Geophysics for Geo-Energy and Geo-Resources | Springer Nature Link

AI要約

1. 背景
石炭鉱山における災害(斜面崩落、ガス突出など)の防止やガス抽出の効率化には、正確な地層(岩相)情報が不可欠です。従来の技術では、技術者が掘り屑(カッティングス)の色や形状を目視で確認していましたが、精度が低く、個人の経験に依存するという課題がありました。近年、振動、音響、画像解析を用いたAI手法も研究されていますが、地下の過酷な環境やノイズの干渉により、実用化には至っていませんでした。特に、「ドリルストリングの摩擦抵抗」が掘削データに与える影響や、膨大なデータのリアルタイム処理が、機械学習を導入する上での大きな障壁となっていました。

2. 手法
研究チームは、以下のハードウェア、前処理技術、機械学習アルゴリズムを組み合わせたシステムを構築しました。

掘削機(スマート掘削リグ): MWD(Measurement While Drilling)システムを搭載した自走式スマート掘削リグを開発しました。姿勢制御、ロッド操作、遠隔制御の各サブシステムを備え、掘削データを自動的に記録します。また、地下のリグと地上のワークステーションを専用ネットワークで結ぶ遠隔データ処理プラットフォームにより、大規模なデータのリアルタイム解析を可能にしました。

前処理:
高速データ検索: 掘削ステータスを14のコード(掘削、ロッド接続など)に分類し、膨大なデータから有効な掘削情報を即座に抽出します。
摩擦検知と補正: ドリルストリングの座屈(バックリング)安定性解析を行い、座屈が起きない条件下でビットが岩石を破壊していない時の負荷を計測して、摩擦抵抗を自動検知します。
データクリーニング: 部品の起動・停止時に発生する「衝撃」などの無効な信号をフィルタリングして除外します。

機械学習手法: 複数のアルゴリズムを組み合わせる「ソフト投票型アンサンブル学習」を採用しました。個別の識別器として、SVM(Quadraticカーネル)、決定木(Medium Tree)、KNN(Fine KNN)、ニューラルネットワーク(Bilayered NN)の4つを統合し、それぞれの予測確率の重み付き平均から最終的な結果を出力します。

説明変数(入力パラメータ): 補正後の推力(Thrust)、トルク(Torque)、および貫入速度(ROP)、回転数(RPM)の4種類を使用します。

目的変数: 石炭 or 岩石(2値分類)

識別結果の継続性を分析(10セットごとのデータ群で判定)することで、さらに精度を高める手法を導入しています。

3. 結果
識別精度: アンサンブル学習モデルは、テストデータに対して98.79%という極めて高い岩相識別精度を達成しました。
掘削効率の向上: 識別された岩相に合わせて掘削パラメータ(回転数や推力)を最適化した結果、従来のシステムと比較して純掘削効率が18.12%向上しました。
省人化: 従来3人で行っていた作業を1人に削減できる可能性を示し、コスト削減と安全性の向上に寄与しました。

4. 考察
本研究の鍵は、ドリルストリングの摩擦をリアルタイムで検知し、掘削データを補正したことにあります。これにより、地層情報とパラメータの相関性が大幅に強化されました。また、地上と地下をシームレスに結ぶデータプラットフォームが実用性を支えています。
 一方で、課題も残されています。現在のモデルは単一のデータソースに基づいているため、地質条件が大きく異なる地域への適用には限界があります。今後は、掘り屑の画像や掘削振動など、複数の情報を統合した「マルチソース情報」による学習を行うことで、汎用性と精度のさらなる向上が期待されています。

TBM, ドリルジャンボなどと同様に、ノンコアボーリングでもトルクや速度をロギング可能であれば、このような使い方ができるのでしょう。地質を知りたい場合はボアホールカメラを後で挿入すれば良いだけです(実際、数m~10mピッチでノンコアボーリングを実施し、ボアホールカメラから柱状図を作成した経験があります)。ボーリングの数を必要とする場合には、有効な手段でしょう。

機械学習よりも掘削機の方が主要な成果です。既にデータさえあれば分析ができる環境が整っていますので、あとは良いデータを集めるだけですから。


文献:Automated geotechnical logging of core box images using machine learning

Automated geotechnical logging of core box images using machine learning: CIM Journal: Vol 16 , No 1 - Get Access

AI要約

背景
地質工学的な不確実性を低減するためには、信頼性の高いデータの量を増やすことが不可欠です。しかし、初期段階のプロジェクトでは、高品質な地質工学用ボアホールの数は限られており、大部分を占める探査用ドリルのデータは品質が不安定なRock Quality Designation (RQD)などに限定されています。また、探査用の実物コアは、分析のために分割されたり、経年劣化や廃棄によって利用できなくなることが多く、コア写真のみが唯一の信頼できる情報源となることが多々あります。従来のコア写真の再ロギングは、時間がかかるだけでなく、ロガー(記録者)の主観やバイアスに影響されやすいという課題がありました。これらの課題を解決し、大量の既存コア写真を客観的かつ効率的にデータ化するため、機械学習を用いたワークフローが開発されました。

手法
機械学習手法: ディープラーニングの一種であるMask_RCNN(畳み込みニューラルネットワーク)が使用されています。これは、画素を自動的に分類・グループ化して特徴マップを作成し、分類を行うモデルです。

前処理:50〜100枚程度のコア画像に対し、各カテゴリーを表すセグメントをポリゴンで囲んでラベル付けを行います。データセットはトレーニングセットと検証セットに分割され、モデルの評価と調整が行われます。
コアの列を特定し、画像を直線化(リニアライズ)して深度を参照するための別のモデルも利用します。

説明変数(入力データ): コア箱のカラー写真です。一貫した照明と色バランスで撮影された写真が、トレーニングに適しています。

目的変数(分類指標): 以下の3つのシステムが開発されました。
Core Damage Index (CDI): コア径に対する砕屑物の平均サイズに基づき、岩盤の状態を5つのカテゴリーに分類します。

  • クラス1: 砕屑物のサイズがコア径の2倍より大きい(> 2x)。従来のRQD(岩質指標)に相当する良好な状態です。
  • クラス2: 砕屑物のサイズがコア径の1〜2倍(1 - 2x)。SCR(固体コア採取率)に相当します。
  • クラス3: 砕屑物のサイズがコア径の0.5〜1倍(0.5 - 1x)。
  • クラス4: 砕屑物のサイズがコア径の0.25〜0.5倍(0.25 - 0.5x)。
  • クラス5: 砕屑物のサイズがコア径の0.25倍未満(< 0.25x)。これには角礫(breccia)、砂、ガウジ(gouge)などが含まれます

Pseudo-RQD(擬似RQD): コア画像上の割れ目と破砕帯(ラブルゾーン)をマッピングして計算されます。自然な不連続面(割れ目)と掘削時などの機械的な破壊を区別しないため、従来のRQDと区別するために「Pseudo」と呼ばれています。
Break Frequency(割れ目頻度): Pseudo-RQDと同様に、割れ目の数から推定されます。

ポストプロセッシング: 5cm未満の間隔のフィルタリング、同じクラスの連続する間隔の結合などのクレンジング処理を経て、以下を表形式でエクスポートします。

  • ホール名(ボアホールID)
  • 開始深度・終了深度(インターバル)
  • 割り当てられたクラス(CDIクラスなど)
  • 該当するコア写真への参照

Leapfrog Geoなどのソフトウェアにこれらのデータをインポートすることで、以下のような高度なモデリングが可能になります。

  • 地質工学ドメインモデルの構築: 抽出されたCDI(コア損傷指数)やPseudo-RQDデータを、地質、構造、岩盤、水理地質の各コンポーネントモデルを構築するためのインプットとして使用します。
  • 断層帯の詳細なモデリング: CDIは5cm単位という非常に高い解像度でデータを生成できるため、断層の核(コア)と損傷帯(ダメージゾーン)を明確に区別し、断層の幅を正確に3D空間で評価できます。
  • 岩盤の変動性の可視化: 従来の1.5m間隔のRQDでは捉えきれなかった岩盤内の細かな品質変化を、3D空間上で連続的に表現できるようになります。

結果
高解像度なデータの取得: CDIを用いることで、従来の1.5m間隔のRQDロギングと比較して4倍以上の高解像度で地質工学的特徴を捉えることが可能になりました。

客観性の確保: ロガーの主観を排除し、砕屑物サイズという客観的な基準で分類することで、一貫したデータセットを構築できました。

損傷の識別: 掘削時の機械的な損傷と、断層核や断層損傷帯のような地質学的な損傷を区別しやすくなり、構造モデリングや水理地質モデリングに貢献するデータが得られました。

考察
客観性の重要性: 機械学習は、主観的な入力が増えるほど、信頼性を維持するために必要な学習データ量が指数関数的に増加します。そのため、RQDのような主観的な判断(自然な割れ目か機械的な破壊かの区別など)を伴う指標は、自動化には不向きであると指摘されています。

デジタル化の限界: コア写真からは岩石の強度や、画像に現れない微細な割れ目の表面状態などを特定することはできません。したがって、この自動ロギングは経験豊富な技術者による詳細なロギングを完全に置き換えるものではなく、既存の膨大なデータを補完するためのツールとして位置づけられます。

モデルの汎用性: 特定の堆積物や条件に過剰適合(オーバーフィッティング)することを避けるため、分類システムを一般化し、客観的な基準(コア径に対する砕屑物サイズなど)を採用することが、信頼性の高いMLモデル構築の鍵となります。

まず思ったのは、オペレータさんによってコアの品質が異なるため、実務に適用するには壁があるだろうな、という点でした。毎回腕の良い、費用の高いオペさんに掘ってもらうことはできないのでしょう。時々、ボロボロのコアに出会うことがありました。
品質を保ちつつ利益を出すというのは企業にとって当たり前のことですが、担当者レベルになるとExecutivesからの圧力もあり、利益が優先される場合があります。これに盾突いて品質を優先する技術者は(程度にもよりますが)、まず民間企業では偉くなれません。結果、利益や株主優先の企業が出てくるのでしょう。

3次元モデラーに入力できるようにしている点は、さすがです。疑似RQDなど、やや品質が落ちたとしても、多数の孔のデータを自動マッピングできるのは基図としてはありがたいと思います。AIに解釈させた結果を提示させるのではなく、画像上の客観的事実のみをAIにより高解像度で提示させたら、後は技術者の仕事です。おかしな点は技術者がログをまとめるなり、モデリングするなりの過程で気づくと思います。

多くのボーリングを整理する必要に迫られた場合、なかなか有向かもしれません。

2026年3月22日日曜日

文献:Machine learning for drill core image analysis

Machine learning for drill core image analysis: A review - ScienceDirect

AI要約

1. はじめに:鉱物探査におけるドリルコアの重要性と機械学習の台頭
世界的な原材料需要の増加、特に再生可能エネルギーへの移行に伴う戦略的原材料の確保が急務となる中、鉱物探査プロセスの効率化が求められています。ドリルコアは、地下の地質状態を物理的に保存した最も直接的な記録であり、その解析(ロギング)は探査の根幹をなします。
近年、このドリルコア写真をデジタルデータとして機械学習(ML)で解析する試みが急速に進化しています。ドリルコア写真は、低コストかつ迅速に撮影可能で、岩相、変質、鉱物、層序、地盤工学的特徴といった多様な視覚的指標を保存しているため、MLの入力データとして非常に魅力的です。

2. 機械学習が取り組む3つの主要な地質学的タスク

① 岩相予測(Lithology Prediction)
岩石の種類(砂岩、大理石、玄武岩など)を分類するタスクで、最も研究が進んでいます。
手法と精度: 畳み込みニューラルネットワーク(CNN)が主流であり、ResNetやGoogLeNetなどのアーキテクチャを用いた研究では、90%から96.7%という極めて高い正解率が報告されています。
技術的アプローチ: 多くの場合、ImageNetで学習済みのモデルをドリルコアデータで微調整する「転移学習」が用いられます。
課題: 砂岩と石灰岩のように視覚的に酷似した岩石の判別には限界があり、人間による化学テスト(塩酸反応など)を代替するには至っていません。

② 地盤工学的解析(RQD/FS推定)
岩盤の品質指標(RQD)や亀裂間隔(FS)を推定するタスクです。
成果: モデルによる推定誤差は手動測定に対して±2.5〜3%以内に収まり、処理速度は人間の約1800倍に達する可能性があるとされています。
汎用性: 亀裂や破砕帯といった構造的特徴は、岩石の種類よりも普遍的なパターンを持つため、最も自動化や汎用化に適した領域と考えられています。

③ 鉱物予測(Mineralogical Prediction)
画像テクスチャ(質感)から特定の鉱物の含有量や品位を推定します。
手法: 従来のテクスチャ記述子(LBPやGLCMなど)を用いる手法と、CNNを用いる手法があります。
複雑性: 鉱物学的性質は微視的なスケールで決定されることが多いため、画像というマクロな視覚情報から推測するこのタスクは、3つの中で最も抽象度が高く困難です。

3. 技術的・実務的な主要課題と議論
研究レベルでの高精度とは裏腹に、実務への完全導入にはいくつかの大きな壁が存在します。

汎用化の壁(Generalization): ある現場で学習したモデルを別の現場や地域に適用すると、精度が著しく低下します。これは、同じ岩相でも地域によって見え方が異なるためです。
データの一貫性と主観性: 地質学者が作成する「正解ラベル」自体に主観が含まれており、専門家間でのラベル不一致(label non-agreement)がモデルの学習を阻害しています。
データの不均衡: 自然界では特定の岩相や希少鉱物は少なく、データが偏っているため、稀なケースをAIが学習しにくいという本質的な問題があります。
スケールの不一致: AIが数ミリ単位で予測を行うのに対し、地質学的な解釈はメートル単位の連続性を重視します。このギャップを埋めるための「予測結果の平滑化(smoothing)」などの後処理が重要視されています。

4. 将来展望:自動化から「コンテキスト認識型アシスタント」へ
ヒューマン・イン・ザ・ループ(HITL): AIが予測を行い、不確実な箇所のみを人間が修正し、その結果をAIが再学習する「地質学者・イン・ザ・ループ」のサイクルが提案されています。これにより、現場固有の文脈(コンテキスト)を効率的にモデルに教え込むことができます。
コンテキスト認識: 空間的な地質モデルや層序学的な知識をMLの制約条件として統合することで、「古い岩石の上に若い岩石があるのは不自然」といった地質学的妥当性を持った予測が可能になります。
ベンチマークの必要性: 異なるモデルを公平に比較するため、多様な地質環境を網羅した「標準的なデータセット群(ベンチマーク・スイート)」の構築が提案されています。

結論として、機械学習は地質学者の仕事を奪うものではなく、解釈の矛盾を指摘したり、膨大なルーチンワークを高速化したりする「高度な助手」へと進化していくことが期待されています。

ルレオ工科大学(スウェーデン)の研究者を中心とし、スウェーデンの大手鉱業企業であるBoliden Mines の実務家が加わった、学術界と産業界の多角的な専門家チームによって執筆されています。レビューなので新たな知見や目新しい手法がないのは仕方ありませんが、それでもHITLの方針はやや残念。いえ、私も(手を抜いて)失敗しましたし、現状では納得なのですが、この壁を破る必要があります。ま、要素ごとに解決して組み合わせる、いずれ丸っとログを吐かせるなどの手法は時間の問題でしょう。
https://phreeqc.blogspot.com/2025/04/blog-post.html

さらに、TBMのようにボーリングマシンでもデータをとれるような仕組みを考え、そのデータを柱状図作成用の機械学習モデルに取り込み、3次元地質モデルを推定させる段階までは、近い将来可能でしょう。要素技術としては既に存在するため、難しくはないと思います。

国内ではこのようなコアに対する機械学習モデルやデータセットを公開しようといった動きを聞きませんし、機械学習自体を使えない技術者たちが Executives になりつつあるため、当面は遅れたままでしょう。ビジネスチャンスではあるので、フットワークの軽い企業は文献に紹介されているような要素技術を商品として宣伝を始めるでしょう。そういった企業がやがて地質の機械学習に関する老舗となり、時代が追い付いた段階でに優位に立てるのかもしれません。


2026年3月21日土曜日

文献:Automated classification of drill core photographs

Automated classification of drill core photographs: A complete workflow for weakly supervised training - ScienceDirect

AI要約

1. 背景
鉱物探査において、地質学者がコアを視覚的に検査することは伝統的な手法ですが、分類が部分的で主観的になりやすく、専門家によって結果が異なるという課題があります。近年、深層学習(CNNなど)を用いた自動予測が注目されていますが、学習に必要な大量のラベル付きデータを作成するには多大な時間と労力が必要です。本研究は、最小限のラベル付けで効率的に岩石を分類できる完全なワークフローを構築することを目的としています。

2. 手法
本手法は、以下の5つの手順で構成されています。

前処理 (Preprocessing)
レンズ歪み補正: 広角レンズによる「樽型歪み」を補正し、その後の直線検出の精度を高めます。
コア箱の分離: 木製のコア箱から岩石コア部分を分離するため、Cannyエッジ検出やモルフォロジー演算(膨張、収縮、オープニングなど)を用います。
細分化 (Slicing): 分離されたコア画像を、解析に適した標準的なサイズ(300×60ピクセル、実寸約18cm長)に分割します。

機械学習手法

自己組織化マップ(SOM)によるクラスタリングとラベル付け
SOMは、高次元のデータを視覚的に理解しやすい低次元(通常は2次元)の格子状マップに変換する「トポロジー保持型」の無向ニューラルネットワークです。
仕組みと構造: 抽出された96個の特徴量(色記述子48個、テクスチャ記述子48個)を入力として使用します。これを10×10(合計100個)のニューロン格子に投影します。
トポロジーの保存: 色や質感が似ている画像は、マップ上でも物理的に近いユニットに配置されます。この性質により、数万枚の画像を1枚ずつ確認するのではなく、似た画像が集まった「クラスター」ごとに効率よく処理できます。
効率的なラベル付け: 地質学者は各ノードから代表的な画像(例:各ノードからランダムに選ばれた5枚)を確認し、そのクラスター全体に対して「漂白」や「片麻岩」といった地質学的なラベルを割り当てます。これにより、数時間という短時間で、数千から数万枚規模の学習用データセットを作成することが可能になります。

深層学習モデル:ベースラインとVGG16
ベースラインモデル (Baseline Model):
本データセットに合わせて開発されたシンプルな3層の畳み込みニューラルネットワーク(CNN)です。
構造: 3つの畳み込みブロック(サイズ32, 32, 64の3×3カーネル、ReLU活性化関数)で構成され、各層の後にバッチ正規化とマックスプーリングが適用されます。最後にソフトマックス活性化関数を持つ全結合層が配置されています。
VGG16: ImageNetデータセットで事前学習された、13層の畳み込み層と5層のマックスプーリング層を持つ非常に深く複雑なアーキテクチャです。
手法: 転移学習(Transfer Learning)のアプローチを採用し、最後の全結合層のみを今回の岩石データで訓練しました。
比較の結果: いずれの地質ドメインにおいても、複雑なVGG16よりもシンプルなベースラインモデルの方が高い精度を示しました。

半教師あり学習(Mean Teacher)
ラベル付きデータが不足している場合に、ラベルのない大量の画像を有効活用して学習の安定性と精度を高める手法です。
「生徒」と「教師」の2モデル構造: 同じ構造を持つ「生徒(Student)」モデルと「教師(Teacher)」モデルを使用します。
一貫性コスト (Consistency Cost): 同じ画像に対して、それぞれ異なるランダムなノイズ(輝度変化や反転など)を加えて両方のモデルに入力します。「滑らかさの仮定(Smoothness Assumption)」に基づき、小さな摂動(ノイズ)を加えても予測結果は変わらないはずであるため、2つのモデルの予測結果の距離(平均二乗誤差:MSE)を一貫性損失として計算し、これを最小化するように学習します。
重みの更新: 生徒モデルは通常の学習(損失関数の最小化)を行いますが、教師モデルの重みは、生徒モデルの過去の重みを指数移動平均(EMA)で反映させることで、より安定した「smoothed(平滑化された)」モデルとして機能します。
利点: 特にサンプル数が極端に少ないクラス(例:ダイアテクサイト)や、複雑な結晶質岩(基盤岩)の分類において、精度の向上に寄与しました。

説明変数 (Explanatory Variables / Features)
色情報: 3Dカラーヒストグラム(RGB空間)から抽出した512個の記述子を、主成分分析(PCA)で48個に圧縮して使用します。
テクスチャ情報: ガボールフィルタ(Gabor filters)を用いて、24通りの畳み込み演算から平均と標準偏差を算出し、計48個の特徴量を取得します。

目的変数 (Objective Variables / Target)
地質学的ドメインごとに以下のクラスを設定しています。
盆地(堆積岩): 漂白(Bleached)、部分除去(Partial removal)、ヘマタイト化(Hematized)、非コア(Non-core)。
基盤岩(結晶質岩): 片麻岩(Gneiss)、メタテキサイト(Metatexite)、ダイアテクサイト(Diatexite)、ペグマタイト(Pegmatite)、非コア。

評価関数 (Evaluation Function)
教師あり学習: カテゴリカル・クロスエントロピーを損失関数として使用します。
半教師あり学習: 教師あり損失(クロスエントロピー)と、教師モデルと生徒モデルの予測の一致度を示す「一貫性損失(MSE:平均二乗誤差)」の合計を最小化します。

3. 結果
盆地サブセット: ベースラインモデルを用いた教師あり学習で、0.932という最も高い分類精度を達成しました。
基盤岩サブセット: 半教師あり学習(Mean Teacher)を採用したベースラインモデルが0.751の精度で最高の結果となりました。
いずれのケースにおいても、複雑なVGG16よりシンプルなベースラインモデルの方が高い精度を示しました。

4. 考察
精度の差: 基盤岩の精度が盆地より低いのは、岩相の複雑さ、特定のクラス(ダイアテクサイト)のサンプル不足(40サンプル)、および変成度の変化が連続的でクラス間の境界が曖昧であるためです。
バイアスの軽減: SOMを用いたクラスタリングにより、地質学者の主観による「オペレーター・バイアス」を一貫したアルゴリズムで補完し、再現性の高い分類が可能になります。

今後の展望: 照明条件の統一や、Segment Anything Model (SAM)のような高度な画像分割手法を導入することで、さらなる精度の向上が期待されます。

まだコアを分けるという段階で、柱状図を作ることはできません。ラベリングの手間をなくすためにクラスタリングを利用するという主張はよく見ますが、経験的には正確ではない場合が多くあります。いずれも発展途中なのでしょうが、いずれ解決すると思います。


文献:Direct inversion of surface wave dispersion

Direct inversion of surface wave dispersion for three-dimensional shallow crustal structure based on ray tracing: methodology and application | Geophysical Journal International | Oxford Academic

AI要約

背景
表面波トモグラフィーは、地球の地殻や上部マントルの構造を研究する上で重要な役割を果たしてきました。従来の一般的な手法は、まず観測された走時データから2次元の位相速度(または群速度)マップを作成し、次に各地点の分散データのインバージョンにより1次元の横波(S波)速度プロファイルを求めるという、2段階のプロセスを踏みます。しかし、この手法や従来の直接反転法では、多くの場合、大圏コース(Great Circle)伝播を仮定しており、複雑な構造を持つ地殻浅部などでは、波の伝播経路の歪みが無視できず、推定結果にバイアスが生じる可能性がありました。一方で、波形全体を用いるフルウェーブフォームインバージョン(FWI)やアジョイントトモグラフィーは高精度ですが、特に高周波数帯域において計算コストが非常に高いという課題があります。

手法
本論文で提案された手法は、中間的な速度マップを作成せず、分散データから直接3次元のS波速度構造を求める単一ステップのインバージョン法です。主な特徴は以下の通りです。

Fast Marching Method( FMM) によるレイトレーシング: 各周期において、FMMを用いて震源・受信点間の走時と波の伝播経路を計算します。これにより、不均質メディアにおける大圏コースからの経路の逸脱を考慮することが可能になります。

反復的な更新: 反転プロセスにおいて、新しく得られた3次元速度モデルに基づいて、波の経路とデータ感度行列(カーネル)を逐次更新します。

ウェーブレット基底を用いた疎性制約付き反転: モデルパラメータをウェーブレット領域で表現します。ウェーブレット変換の多重解像度特性を利用し、データ制約が強い地域では高解像度を、制約が弱い地域では長波長の特徴を保持するようにします。

L1ノルム正則化: ウェーブレット係数がまばら(スパース)であると仮定し、L1 ノルム正則化を適用します。これにより、従来の滑らかさを強制する正則化よりも、データの分布に適合した適応的な解像度が得られます。

計算アルゴリズム: L1ノルムベースの最小化問題を解くために、反復再重合最小二乗法(Iteratively Reweighted Least Squares: IRLS)とLSQRアルゴリズムを組み合わせて使用しています。

物性間の関係: P波速度(α)と密度(ρ)は、地殻に関する経験的な関係式(Brocher, 2005)を用いてS波速度(β)に関連付けられており、これらによる分散への影響も感度行列に含まれています。

結果
台湾の台北盆地における、環境ノイズ相関法から抽出された短周期(0.5〜3.0秒)のレイリー波位相速度データにこの手法を適用しました。

計算の収束: 反復的な反転により、走時残差の標準偏差は1.44秒から0.91秒に減少し、残差の平均はほぼゼロになりました。

速度構造の解明: 盆地内の低速度域と山地(大屯火山群など)の高速度域が明瞭に描き出されました。これらの結果は地表の地質構造と良好に一致しています。

従来手法との比較: 従来の2段階反転法の結果と概ね一致しましたが、データ網羅性が低い地域では、本手法の方がアーティファクトが少なく、より安定した結果が得られました。

解像度テスト: チェッカーボードテストにより、データ制約が良い地域では水平・垂直方向に高い解像度を持つことが確認されました。

考察
提案手法は、計算効率と精度のバランスに優れており、特に構造が複雑な地殻浅部の調査に有効です。

利点: 2次元速度マップの構築が不要であり、3次元参照モデルの組み込みが容易です。また、ウェーブレット基底の採用により、データの密度に応じてモデルの細かさが自動的に調整されます。

限界と仮定: 現状では高周波近似(レイ理論)に基づいており、滑らかな不均質構造を想定しています。

今後の展望: 有限頻度効果(Finite-frequency effects)を取り入れた感度カーネルの導入や、実体波(P波・S波)の走時データとの共同反転(Joint Inversion)への拡張が期待されます。これにより、地殻浅部におけるP波・S波両方の速度構造をより正確に制約できる可能性があります。

こちらも、数式を追うのに時間がかかりました。パスを計算する点や、1次元のインバージョンを介さないという点が新鮮です。



2026年3月18日水曜日

文献:3D Wave‐Equation Dispersion Inversion of Distributed Acoustic Sensing Data


AI要約
3D波動方程式分散インバージョン(3D-DAS-WD)法の解析手順
3D-DAS-WD法は、DASひずみデータに直接基づく3D波動方程式分散インバージョン手法であり、ひずみ成分を速度に変換する際の数値誤差を回避し、低SNR(信号対雑音比)のDASデータにおけるFull Waveform Inversion(FWI)の収束問題に対処することを目的としています。この手法は、浅い地表付近の地質構造を高分解能で3Dイメージングする能力を示しており、地下の空洞などの特徴を効果的に特定できます。

1. DASの基本原理とデータ取得
• DASの測定原理: 従来の地震計が粒子速度を測定するのに対し、DASは光ファイバーのゲージ長に沿ったひずみ速度を測定します。

2. 初期S波速度モデルの作成
• 解析の出発点として、初期S波速度モデルを準備します。このモデルは、反復的なインバージョンプロセスを通じて改善されていきます。

3. 予測ひずみ速度DASデータの順モデリング
• 作成した初期S波速度モデルを用いて、forward modeling を実行し、予測されるひずみ速度DASデータを生成します。この研究では、スペクトル要素法に基づく弾性モデリングコードSPECFEM3Dが使用されています。
• 計算される速度成分は、式(14)および(15)を使ってDASひずみ速度成分に変換されます。

4. 観測データと予測データの分散曲線の導出
• MASW として分散曲線を作成します。
• 分散曲線は、linear Radon transform を用いて計算された分散スペクトルの最大値に基づいて抽出されます。
• この研究では、各ショットから10度間隔で分散曲線が抽出されます。

5. 目的関数の定義と勾配の計算
• 目的関数(Objective Function): 3D-DAS-WD法は、観測された波数κobs (θ,ω)と予測されたDASデータ表面波分散曲線κ(θ,ω)の二乗差の合計を最小化します。目的関数Jは以下のように定義されます: 
J = (1/2) ∑θ ∑ω (κ(θ,ω) − κobs (θ,ω))^2
ここで、θは方位角、ωは周波数です。

• S波速度に対する勾配: S波速度s(x)に対する目的関数の勾配δ(x)は、以下の式で計算されます: 
δ(x) = ∂J / ∂s(x) = ∑θ ∑ω Δκ(θ,ω) (∂κ(θ,ω) / ∂s(x))

• Fréchet 微分の導出: Fréchet 微分 (∂κ(θ,ω) / ∂s(x)) は、予測データと観測データのk-ω領域における相関を示す接続関数(connection function)Φを用いて導出されます。これは、随伴状態法(adjoint state method)を使用して計算されます。最終的にS波速度 s を λ (ラメの第一定数) と μ (剛性率) に関連付けることで、Fréchet 微分が導出されます
∂κ(θ,ω) / ∂s(x)
 = − (1/A) R{∫ ∂D(g,ω)/∂s(x) D̂*obs(g,θ,ω) dg} 
 = − (1/A) R{∫⟨GR ∂u(g,ω)/∂s(x), D̂obs(g,θ,ω)⟩ dg}
= − 4ρs (∂κ(m) / ∂λ) + 2ρs (∂κ(m) / ∂μ) 

・∂κ(m) / ∂λ、∂κ(m) / ∂μは、 波場u(m)(合成地震計データ)とバックプロパゲーションされた残差u'(D̂∗obs(g,θ,ω)をソースとする)の内積の積分として表される。

6. 速度モデルの更新
• ステップ長の決定: safeguarded backtracking line search method を用いて、モデル更新のためのステップ長αが決定されます。
• モデルの反復更新: S波速度モデルは、計算された勾配とステップ長αを用いて、steepest-descent formula により反復的に更新されます。 ここで、

7. 収束判定と最終S波速度モデルの出力
• ステップ3から6のプロセスは、目的関数が指定された誤差閾値以下になるか、または最大反復回数に達するまで繰り返されます。
• 収束条件が満たされた後、最終的なS波速度インバージョンモデルが出力されます。

さすがにAI要約のみでは理解できず、数式を追うのに時間がかかりました。
DAS(1軸センサー)を利用して、どのようにレイリー波を取り出すか、取り出した後、どのようにインバージョンをかけると精度が落ちにくいか、という内容でした。地震分野の知識の蓄積を感じる文献です。
3次元解析を目的としてファイバーをクロスで張るなら、通常のピックを使った方が楽なように感じましたが。ま、一つの方法として覚えておきましょう。