2025年11月13日木曜日

mce=off

数か月前から splash screen で再起動がかかる Ubuntu。

GRUB_CMDLINE_LINUX_DEFAULT に mce=off を追記したら、すんなり起動しました。

何だろう?

2025年11月6日木曜日

薄片画像 + ViT

PoreViT: Automated pore typing in carbonate rocks using vision transformers and neighborhood features - ScienceDirect

AI要約 

背景

  • 炭酸塩岩の孔隙(ポア)をその生成過程や物理的特性に基づいて分類することは、炭酸塩岩の微細構造と物理特性の関係を理解する上で重要である。
  • 従来の孔隙分類は主に手作業で行われ、非効率で主観的、かつ大規模解析には不向きである。
  • 既存の自動化手法は、孔のサイズや形状といった単純な特徴に依存しており、遺伝的クラス(生成過程に基づく分類)を正確に識別できていない。
  • CNN(畳み込みニューラルネットワーク)は画像分類に強みがあるが、局所的特徴の抽出に優れる一方で、画像全体のグローバルな文脈情報を捉えるのが苦手であるため、複雑な孔隙タイプの分類には限界がある。
  • そこで、Vision Transformer(ViT)を用い、局所的特徴とグローバルな文脈情報を融合し、孔隙の近隣情報(隣接する孔の情報)も活用することで分類精度の向上を目指す。

手法

  • PoreViTモデルを提案。これはViTをベースにしたモデルで、薄片画像中のマクロ孔隙をLucia分類(interparticle, touching vug, separate vug)に分類する。
  • モデルの特徴は、ViTの特徴量にGlobal Token Addition(GTA)層を加え、CNNから抽出した空間的特徴と融合するFeature Fusionブロックを持つ点。
  • 入力データは、単一の孔だけでなく、その周囲の孔やテクスチャ情報を含む近隣情報を含めてモデルに与えることで、局所的な孔系のトポロジーを認識可能にしている。
  • データセットは25枚の高解像度薄片スキャンから4115のラベルを取得し、これを用いてモデルを訓練・評価。
  • 従来のCNNモデル(DenseNet121)と比較し、精度(precision)と再現率(recall)で約4%の絶対的改善を達成。

結果

  • テストセットでPoreViTは正確度93.6%、precision 0.92、recall 0.93、F1 0.92を達成。最良CNN(DenseNet121)に対し精度・再現率で絶対+4.0%、相対約+4.5%の改善を示した。
  • クラス別ではinterparticleが最良(Precision 0.95、Recall 0.93、F1 0.94)、誤分類は主にseparate vugで発生。
  • 高解像度薄片1枚(2650孔)での全面推論は約1.5分で完了し、精度93%、Precision 0.931、Recall 0.930、F1 0.927を達成、処理スループットは1枚約2.8分(前処理1.3分+推論/後処理1.5分)で、1日160–170枚のバッチ処理が可能。
  • 近隣情報の活用が分類精度向上に大きく寄与している。
  • ViTの自己注意機構により、画像全体の長距離依存関係を捉えられるため、CNNの局所的受容野の制約を克服できている。
  • これにより、複雑な炭酸塩岩の孔隙タイプ分類において、より正確でスケーラブルな自動化が可能となった。

考察

  • CNNは局所的特徴の抽出に優れるが、グローバルな文脈理解が不足し、複雑な孔隙分類には不十分であることが示された。
  • ViTは画像をパッチに分割し、それぞれをトークンとして扱い、自己注意機構で長距離の関係性を直接モデル化できるため、孔隙の空間的・構造的関係をより深く理解できる。
  • 近隣孔の情報を含めることで、単一孔の形状だけでなく、その周囲の孔との関係性も考慮した分類が可能となり、実際の地質学的解釈に近い分析が実現。
  • 本研究は、炭酸塩岩の微細孔隙構造の高スループットかつ定量的な解析を可能にし、石油工学や地質学における炭酸塩岩評価の効率化と精度向上に貢献する。
  • 限界として、入力画像の高解像度・セグメンテーション品質・注釈品質への依存、224×224・パッチ分割による連続関係の断片化、ViTのデータ要求性、学習分布外リソファシスへの一般化不確実性、クラス不均衡の影響などが指摘される。
  • 将来展望として、Choquette & Prayの高次遺伝的分類への拡張(クラス数増でデータ収集が課題)、3DマイクロCTへの適用による真の接続性評価、Swinなどの多尺度ViTや自己教師あり事前学習、ハイブリッド拡張、予測不確実性の定量化などが提案されている。

 薄片画像 と 機械学習は相性が良いでしょう。地質屋さんなら鉱物の同定に利用するといった発想に至るようですが、特許もあるようなので注意。

2025年9月13日土曜日

3次スプライン関数で崩壊面推定 その2

Assessing the landslide failure surface depth and volume: A new spline interpolation method - ScienceDirect

GitHubにコードあり、と書かれていましたが、本日時点で ReadMe と License しか公開されていません。スプラインを使うだけの簡単な手法なので、書いてみました。

import numpy as np
from scipy.interpolate import CubicSpline
from tqdm import tqdm

def spline_failure_surface_pointwise(
    dem: np.ndarray,
    landslide_mask: np.ndarray,
    fix_points=None,
    max_iter=1000,
    tol=0.01,
    slope_angle_thresh=1
):

    failure_surface = dem.copy()

    fixed_mask = np.zeros_like(dem, dtype=bool)
    fixed_depths = np.zeros_like(dem)
    if fix_points:
        for (i, j), val in fix_points.items():
            fixed_mask[i, j] = True
            fixed_depths[i, j] = val

    nrows, ncols = dem.shape

    print('dem.shape:', dem.shape)
    print('landslide_mask.shape:', landslide_mask.shape)

    for it in tqdm(range(max_iter)):
        prev_surface = failure_surface.copy()

        for i in range(nrows):
            for j in range(ncols):
                if not landslide_mask[i,j]:
                    continue

                if fixed_mask[i,j]:
                    failure_surface[i,j] = fixed_depths[i,j]
                    continue

                j_idx_list = [jj for jj in range(ncols)
                             if landslide_mask[i,jj] and jj != j]
                if len(j_idx_list) >= 2:
                    row_vals = [failure_surface[i,jj] if not fixed_mask[i,jj] else fixed_depths[i,jj]
                                for jj in j_idx_list]
                    row_spline = CubicSpline(j_idx_list, row_vals, bc_type='not-a-knot')
                    zintp_i = row_spline(j)
                else:
                    zintp_i = failure_surface[i,j]

                i_idx_list = [ii for ii in range(nrows)
                             if landslide_mask[ii,j] and ii != i]
                if len(i_idx_list) >= 2:
                    col_vals = [failure_surface[ii,j] if not fixed_mask[ii,j] else fixed_depths[ii,j]
                                for ii in i_idx_list]
                    col_spline = CubicSpline(i_idx_list, col_vals, bc_type='not-a-knot')
                    zintp_j = col_spline(i)
                else:
                    zintp_j = failure_surface[i,j]

                # 両方向平均
                zavg = (zintp_i + zintp_j) / 2

                # 文献式:
                if zavg < failure_surface[i, j]:
                    failure_surface[i, j] = zavg

        # 収束判定(標高変化量)
        change = np.abs(failure_surface - prev_surface).max()
        if change < tol:
            print(f"Converged at iteration {it}")
            break

    return failure_surface

制御点をある関数にあてはめて残りを推定するというのはGEORAMA等と同じ考え方。GEORAMAはGUI(CAD)を利用しているので制御点を与えやすいのがメリット。

適当な地形を作って動かしてみたのが以下。


結果はイマイチ。イタレーションを3000回にしたのですが、まだ収束していません。これでも数時間かかったのですが、並列化しても半分ぐらいでしょう。DEMの解像度が高い、あるいは崩壊範囲が広いと時間のかかる手法だということがわかりました。こうなると、スプラインよりもGEORAMAの方が優秀、SLBLの方が簡素。関数を変更したら良いのでしょうね。https://www.blogger.com/blog/post/edit/269984068930159765/2223718948702813275

ま、ボーリングによる崩壊面深度を反映する場合の簡易手法ということで、覚えておきましょう。

3次スプライン関数で崩壊面推定

Assessing the landslide failure surface depth and volume: A new spline interpolation method - ScienceDirect

AI要約

背景
土砂崩れの滑落面の深さや形状、移動した土砂の体積の正確な推定は、災害評価、ボーリング計画、対策のために不可欠である。土砂崩れの滑落面が完全に露出している例は少なく、データ不足の地域では正確な評価が課題となる。これに応じて、コスト効率のよいスプライン補間法を用いた新しい手法が提案された。

手法
データ準備:
Hope SlideについてはLiDARデータから取得した8mDEM、Downie Slideについては30mDEMを使用。
スプライン補間法:
滑落面を3Dで再構築するために3次スプライン補間を使用。滑落面の深さや体積を推定するために反復計算を行い、極端な境界条件(谷先、フランク、冠部など)を利用して滑落面の深さや形状を最適化し、滑落面の「最浅」「中間」、「最深」の状態を算出する。既存データ(ボーリングデータ、露出面など)を反映させることで結果の精度を向上させることが可能。
推定処理のステップ:

  1. Landslide の境界を手動で定義し、その内部を分析領域として設定。
  2. スプライン関数はx方向(東西軸)およびy方向(南北軸)に設定。必要に応じて座標系を回転させることで、データを地すべり方向に合わせる。
  3. 対象グリッドノードを一時的に削除し、その周囲の残りデータに基づいて該当部分を再計算する。削除したノードの補間値をスプライン計算に戻し、またその部分を全体の滑落面に統合する。
  4. 各グリッドノード(i, j)に関して、反復回数 ( t ) 毎に生成されたデータ(例えば標高値 )を前回の結果と比較。推定した滑落面が収束するまでプロセスを反復。
  5. 確率密度関数(KDE: Kernel Density Estimation)を用いて、滑落面の深さに関する確率を算出。

結果
Hope Slideでは最深の滑落面での体積を65百万m³。既存のドナー指標(Donati, 2019)で示された体積48.4百万m³よりも約25%多い。
露出した滑落面領域をモデルに追加すると、体積が61百万m³となり、精度が向上。
Downie Slideでは最深の滑落面での体積は880百万m³。既存データに基づく過去の推定体積(1000億m³)との比較では、情報不足の下でも信頼性の高い結果を提供可能であった。ボーリング深度を追加してモデルを改良した結果、体積は889百万m³となり、精度が向上した。

考察
過去の標準的手法(カットアンドフィル法やGISベースの手法)よりも結果が正確であり、特に滑落面が非露出である場合の精度が向上。
手法は追加データ(ボーリング位置、露出面など)を組み込む柔軟性を持ち、精度をさらに高める可能性がある。
本手法が全てのケースに適用可能ではない点や、結果を解釈する際の注意点(例えば、境界条件の適切な設定)があることに留意する必要がある。

確率分布の表示に惹かれましたが、提案手法(スプライン補間とその反復計算)の過程で得られた値をもとに統計処理(カーネル密度推定等)を行って得られる、手法依存の確率分布でした。手法を適用できることが前提での値であり、実際のすべり面が存在する「物理的な確率」を直接表すものではありません。ま、これはどの手法にも言えるものですけれど。

そうなると、国内ではどの手法(関数)がどのようなタイプの崩壊や地すべりに適用できるか?という整理が必要です。これがないと、3次元形状を推定して3次元の安定解析を進めるにしても、その根拠はこれまで通り「経験」止まりで科学には至りませんし、説明性に劣ります。確率分布も意味のない計算になります。

これは、地質屋さんの出番ですね。


2025年9月6日土曜日

崩土の移動距離

Empirical Estimation of Landslide Runout Distance Using Geometrical Approximations in the Colombian North–East Andean Region

AI要約

背景
地滑りは、部分的な斜面の静的平衡喪失による地質的災害であり、インフラ・環境・経済への重大な影響を及ぼすことがある。特にコロンビア北東アンデス地域では、20世紀初頭以降約1.3百万米ドルの損害が報告されており、道路、建物、農地などが影響を受けている。地滑りの移動距離(ランアウト距離, LRD)はリスク評価において重要であり、特に幾何学的近似を用いた経験的モデルは広域的な研究や重要区域の特定に役立つ。コロンビア北東アンデス地域は地滑りが多発するエリアであり、地質や地形、土地利用などが複雑に絡み合う環境である。

手法
データ収集:
サンプルは、コロンビア地質調査局の地滑りインベントリからランダム抽出されたデータセットを基に、無人航空機(UAV)を使用して地形を詳細に調査。計49の地滑りが分析対象に選ばれた。
幾何学的モデルと推計式:
地滑りの特性(深さ、幅、長さ)、全高低差、高低差に基づく傾斜角などを計測。これらを用いて経験的な幾何学的モデルを適用し、ランアウト距離を推定。
既存の経験的モデル(e.g., Finlay et al. 1999, Hunter and Fell 2003)を検証し、分析されたデータセットに最適なモデルを選定した。

結果
データ分析:
地滑りランアウト距離は、土地のカバー(作物地、水域など)や地質的特性(e.g., Bocas層での砂岩)に依存することが確認された。
傾斜角度や地滑り体積とランアウト距離の間に統計的な相関があることが示された。特にランアウト距離と傾斜角においては、相関係数R²=0.32である。
モデルパフォーマンス:
Finlay et al. (1999)やHunter and Fell (2003)モデルなど既存のモデルを検証した結果、特定のモデルが他のモデルよりも良好な精度でランアウト距離を推定できることがわかった。
誤差率(RE)はモデル間で異なり、推計結果が過小評価される場合が多かったが、一部モデルは過小評価を抑えつつ信頼性の高い予測を提供した。

考察
傾斜角や地形特性はランアウト距離(LRD)の評価において鍵となる要素であるが、既存モデルによる推計の精度には依然として課題がある。誤った過小評価は、工学的実務において深刻なリスクを引き起こす可能性があるため、過大評価が望ましいとされた。
高度落差(H)に依存しないLRD推定方法の提案が行われており、これにより評価プロセスが単純化される可能性がある。
将来的には、地滑りリスク評価や高リスク地域の特定における幾何学的近似モデルの改善が必要であり、対象地域全体にわたる実証的データの収集拡張が重要である。


機械学習ではなく、重回帰。
特徴量に土砂量が使われていますが、推定式から求められています。それでも役に立つといった点は、特徴量エンジニアリングで利用できそうです。試してみましょう。

ポリゴンの崩壊域と流動/堆積域の自動区分

Towards automatic delineation of landslide source and runout - ScienceDirect

AI要約

1. 課題の背景と重要性
従来の地滑りインベントリでは、発生源と流動域が区別されず、一体として記録されることが多いため、予測モデルの精度や体積推定、ハザード評価に悪影響を及ぼしていた。発生源と流動域は力学的特性が異なるため、これらを分離することで、地滑りの発生メカニズムや伝播、堆積パターンなどをより深く理解し、より信頼性の高いハザード・リスク評価、体積推定、運動学的シミュレーションが可能になる。

2. アプローチの概要
モデルは、地滑りのトポロジー情報(形状の複雑さ)と形態情報(物理的な属性)を組み合わせて使用する。これらの特性を特徴量として、機械学習モデル(Random Forest)を訓練し、発生源領域を識別する。

3. 入力データ(特徴量)
モデルへの入力となる特徴量は、主に以下の2種類。

トポロジー特性:
代数トポロジーの手法を用いて、地滑りの地形における「穴(holes)」や「連結成分(connected components)」といった複雑な形状を数値化する。具体的には、パーシステントホモロジー(persistent homology)の概念に基づき、平均寿命(Average Lifetime)、点の数、ベッチ曲線(Betti-Curve)、ヴァッサースタイン振幅(Wasserstein Amplitude)、ボトルネック振幅(Bottleneck Amplitude)などの指標が計算される。これらは地滑りのコンパクトさ、伝播ゾーンの屈曲度、斜面の勾配の変化などを示唆する。最も重要なトポロジー特性には、穴のボトルネック振幅(BA_H)や連結成分の平均寿命(AL_C)が含まれる。

形態特性:
地滑りの物理的属性を示す要素で、ヘイムのエネルギールール(Heim energy line)に基づいた移動距離(travel distance)、移動角度(travel angle)、移動高さ(travel height)、および推定速度(velocity)などが含まれる。これらは地滑りの幾何学から導出され、運動学的指標として機能する。重要な形態特性には、移動距離、移動角度、移動高さ、平均速度が含まれる。これらの特性の中から、相互相関の高いものや特徴量重要度の低いものは排除され、モデルの予測精度を高める上位10の特徴量が特定された。

4. 回帰問題としての定式化
この研究では、発生源領域を地滑り全体の伝播長に対する発生源領域の伝播長の比率として識別する。この比率は、スカープ領域を囲むバウンディングボックスの伝播長と、地滑り全体を囲むバウンディングボックスの伝播長を用いて算出される。回帰モデルの主な目標は、予測された比率と実際の比率との間の誤差を最小化することとなる。

5. モデルの訓練と評価
訓練データ:
ドミニカ、ネパール、トルコ(バルトゥン)、イタリア(ベッルーノ)、日本(新潟)などの地域からの、約30,000の地滑りサンプルを使用してモデルを訓練および評価した。これらのデータセットの一部には、すでに発生源と流動域が分離された「グラウンドトゥルース」が含まれており、モデルのテストと予測精度評価に利用した。

訓練手法:
10 k-fold交差検定を用いた。

精度:
モデルの平均偏差指標は15%未満で、標準偏差は約8%であった。ドミニカで平均10.37%、ネパールで平均7.48%(最も良い結果)、トルコで平均12.30%の偏差を示した。

6. モデルの堅牢性
このモデルは、さまざまな地理的設定やトリガーの種類(降雨、地震、歴史的イベント)の地滑りに対して、発生源領域を高い精度で特定できることが示された。また、スライドやフローといった地滑りの動きの種類の違いに関わらず、発生源領域を効果的に識別できることも示された。

まさかのトポロジー。以前、読んだけど利用していないなと思いつつ、見返すと6年前でした。https://phreeqc.blogspot.com/2020/01/blog-post_24.html

日本のデータも含まれています。が、利用は他国の方。うーん。

崩壊土砂量推定式 V=αA^γ

崩壊生産土砂量推定式のパラメータ<i>α </i>値・<i>γ </i>値に及ぼす崩壊深や地質の影響 -2014年8月豪雨で発生した丹波地域と広島地域の土砂災害の事例-

AI要約

1. 背景
斜面崩壊地から生産される土砂は土石流となり、下流域に甚大な災害を引き起こすため、崩壊生産土砂量の早期かつ精度良い推定が求められている。崩壊生産土砂量(V)は、崩壊面積(A)と崩壊深の積で概ね求められ、V=αA^γ の累乗式で近似できることが知られている。これまで、様々な国で崩壊生産土砂量の推定手法が研究され、多くの研究者がα値・γ値を報告してきた。しかし、既存研究にはいくつかの問題点が存在する。

  • 計測値の不正確さの可能性: 崩壊深や崩壊面積の値が先行研究からの引用であり、計測者の技量や解析精度の粗さによるばらつきがある。
  • サンプル偏り: 特定の地域や規模の大きな崩壊地のデータがα値・γ値の計算に影響を与えている可能性がある。
  • 地質分類の欠如: 誘因や素因が異なっていても同一のα値・γ値を使用する研究がある一方で、地質によって崩壊の形態や規模が異なることが指摘されており、地質を考慮する必要がある。
本研究は、LP地形データを用いて崩壊深や地質が崩壊生産土砂量推定式のパラメータα値・γ値に及ぼす影響を明らかにすることを目的とし、既存研究のα値・γ値との比較も行う。

2. 手法
本研究では、2014年8月豪雨で土砂災害が発生した兵庫県丹波市と広島県広島市を調査地域とした。

崩壊地の抽出と土砂量・崩壊深の計算:
災害前後の1m分解能LP地形データを使用し、ArcGIS Pro2.8を用いて標高値の変化量を計算した。
崩壊地は、斜面勾配が20度以上かつ外縁部に厚さ0.5m以上の侵食が発生した領域を対象とし、空中写真で確認して認定した。
崩壊地ポリゴンごとに、侵食部の面積(Ae)と堆積部の面積(Ad)を計算し、崩壊面積Aは両者を合算した「A = Ad + Ae」で算出。
崩壊生産土砂量Vは、侵食量(Ve)と堆積量(Vd)を差分した「V = Ve - Vd」で求めた。
崩壊深Dは、崩壊地ポリゴン内における侵食部の標高値の変化量の平均値として算出。

α値とγ値の計算:
崩壊面積Aと崩壊生産土砂量Vから、累乗式V=αA^γでフィッティングすることにより、地質別にα値とγ値を計算。
対象とした地質は、丹波地域の「丹波帯頁岩砂岩」と「超丹波帯砂岩」、広島地域の「花崗岩類」「流紋岩類」「苅田層泥岩」「玖珂層泥岩」の計6種類。

3. 結果
崩壊面積Aと平均崩壊深Dの関係:
崩壊面積Aと平均崩壊深Dは、D=mA^nで近似すると、ばらつきがあるものの、地質に関係なく概ね正の相関がみられた。近似直線の傾きを示すn値は0.03〜0.36の正の範囲にあり、切片を示すm値は0.20〜1.27の範囲であった 。m値とn値の間には決定係数R2値が0.94の強い負の相関が確認された 。

地質別の崩壊生産土砂量とα値、γ値:
崩壊面積Aと崩壊生産土砂量Vの関係は、V=αAγで近似すると、全ての地質で決定係数R2値が0.996以上と非常に高く、この式で概ねVを推定できることが示された。傾きを示すγ値は0.94〜1.49の正の範囲にあり、切片を示すα値は0.08〜2.02の範囲であった。地質ごとの特徴として、丹波帯頁岩砂岩と苅田層泥岩はγ値が1.34〜1.49と他の地質よりも大きな傾向が見られた。崩壊数の合計は丹波地域が777箇所、広島地域が381箇所。崩壊面積の中央値は約150〜340m2の範囲にあり、丹波帯頁岩砂岩、超丹波帯砂岩、玖珂層砂岩で大きくな傾向があった。崩壊深の中央値は超丹波帯砂岩が1.6m以上と大きい一方で、他の地質は約1.3〜1.5mの範囲であった。

4. 考察
α値とγ値の関係性:
α値とγ値の間には負の相関があることが確認された。
崩壊面積に対する崩壊深の増加割合を示すn値が大きくなると、V=αA^γの傾きに相当するγ値が増大する傾向が見られました。これは、崩壊生産土砂量が崩壊面積と崩壊深の積で表されることと整合する結果であった。α値は極めて小規模な崩壊地の土砂量を意味するため、γ値に比べて地域や地質による差が若干不明瞭になることがある。

地質とα値・γ値の関係:
火成岩類(本研究の花崗岩類、流紋岩類、および既存研究の斑れい岩)は、概ねα値が大きく(0.76〜2.02)、γ値が小さい(0.94〜1.06)傾向に分布する。これは、これらの地質が崩壊面積に対する崩壊深の増加割合が小さい崩壊形状の分布特性を持つためと考えられる。堆積岩類は、α値が0.08〜1.90、γ値が0.97〜1.49と広範囲に分布。これは、地質ごとに崩壊形状の分布特性が多様であることに起因すると考えられる。変成岩(既存研究のホルンフェルス、黒雲母片麻岩など)は、α値が小さく(0.03)、γ値がかなり大きい(1.55)という特徴を示した。これは、崩壊面積に対する崩壊深の増加割合が大きい崩壊形状の分布特性を示唆している。

既存研究との比較:
Guzzetti et al.(2009)やLarsen et al.(2010)が提案するα値・γ値は、本研究の値と比較して、α値が小さくγ値が大きい傾向にある。ただし、丹波帯頁岩砂岩のように大規模な崩壊地を含む場合、本研究のα値・γ値はこれらの既存研究の値に近くなることがわかった。既存研究のα値・γ値を用いて土砂量を推定した場合、Larsen et al.(2010)の岩盤崩壊の場合を除き、本研究で計測された実際の土砂量よりも過小評価される傾向があった。これは、既存研究のα値が小さいため、本研究の調査対象である中小規模の崩壊地が卓越する状況では、崩壊生産土砂量を過小評価してしまうためと考えられる。一方で、Larsen et al.(2010)の岩盤崩壊のα値・γ値を用いた場合、本研究とは崩壊深の規模が大きく異なるため、土砂量が約2倍に過大評価される結果となった。

5. 今後の展望
本研究で扱った地質の種類は限られているものの、今後さらに計測データを増やすことで、α値・γ値を地質ごとにある程度体系的に整理できると考えている。これにより、土砂災害発生時の砂防基本計画などにおいて、崩壊生産土砂量推定式の活用を目指し、地質ごとのα値・γ値のデータを蓄積していくことが期待される。

地質によらず、この形の式を使えそうだという点、地質によって係数が変化する点、大規模か小規模かにより切片αの重要性が変わる点など、V=αA^γ に着目して研究された結果です。降雨、地震によらず式は適用できるが、地域は考慮しないとダメ。というか、不確定な内容を地域(係数)に押し込んでいるのかもしれません。
データを蓄積しないと予測には使えませんし、そもそも正確なポリゴン、崩壊面の推定、土砂量の算出がないと正しい係数を得られません。手間がかかる研究です。


2025年9月4日木曜日

胆振東部地震の土砂量推定式


AI要約
背景
地震による土壌地すべり(EQIL)の体積推定は、造山過程や表層変動の理解に不可欠である。特に2018年北海道胆振東部地震(M6.6)では、多数の浅い土壌地すべりが発生し、その体積推定には高精度なスケーリング関係式が必要とされている。従来の研究では面積と体積の経験的なべき乗則 ( V = α A^γ ) が用いられているが、このスケーリング指数 γ の不確実性が大きく、地域や土質によって異なることも指摘されている。また、高解像度DEM(デジタル標高モデル)の普及により正確な体積推定が可能になってきたものの、多くはフィールド調査や低解像度データに依存していたため誤差も大きかった。本研究では北海道胆振東部地震で発生した1719件の浅い土壌地すべりを対象に、高解像度LiDAR DEMを用いて新たな面積-体積関係式を構築し、不確実性評価も行うことを目的としている。

手法
対象地域:2018年北海道胆振東部地震による約500 km²範囲内。
データ:1m空間分解能の事前・事後LiDAR DEM(事前29.41 km²=校正領域、事後512 km²=検証領域)。
地すべりマッピング:既存3つの異なる土地滑りインベントリ(GSIなど)と、新規作成した校正領域内1719件についてDEM差分法で負変化部分から自動抽出し、人手で写真照合・修正。
フィールド検証:50箇所訪問し深さ測定など現場確認。
面積・長さ・幅算出はGISツール利用。
体積計算はセル単位でDEM差分値×セル面積を集計。
面積-体積関係式( V = α A^γ )  のパラメータα, γ を最小二乗法等で求める。
他文献から得られた既存スケーリング指数との比較、不確実性評価も実施。
地殻隆起量との比較検討にはGNSS観測値から補間した垂直変位マップ利用し総隆起量算出。

結果
新規構築した面積-体積関係式ではスケーリング指数γ=1.15となった。この値は文献中報告されているγ=1.0~1.9範囲内だが、日本国内他研究より若干低めか類似水準。
北海道胆振地域全流域から侵食された総堆砂量は6400万~7200万m³と推定された。これは同時期GNSSデータから求めた共震隆起量約555万m³より遥かに大きく、大規模地震および豪雨による浸食作用が局所的隆起効果を上回っていること示唆する。
複数インベントリ間および既存文献との比較分析では、材料特性(土質)、地域特有条件(火山灰層厚さ等)、解析方法違いなどによってγ値やα値にはばらつきあり、不確実性要因となっていた。

考察
本研究成果は、日本特有とも言える厚い火山灰層下に形成された軟弱な土壌条件下で発生する浅層型EQILについて、新しい高精度かつ現場検証済みスケール則を提示しており、その適用範囲拡大への貢献が期待できる。また、大規模EQILイベント時には浸食作用による堆砂量増加が共震隆起効果以上となり得るため、このバランス把握は造山運動や表層プロセス理解にも重要だと論じている。不確実性要因として材料種別ごとの異なるγ設定必要性やDEM精度向上効果にも触れており、更なる多地点多イベント解析への展開可能性も示唆している。
面積-体積関係式( V = α A^γ ) が土質によって異なるよ、という内容です。胆振東部ではγ=1.15でした。
他国の方々が胆振東部を分析されています。どのようなつながりでしょうか?








崩壊土砂量の推定


AI要約
背景
降雨による浅層地すべりの体積(VLDR Volume of Landslides due to Rainfall)の定量化は、リスク管理、緊急対応、土木設計、経済評価、環境保護において重要である。地すべり体積の推定は、地形変化の把握やハザードマップの更新に寄与し、土壌の移動規模を示すことでリスクゾーンの分類(低・中・高)を可能にする。従来の地すべり感受性モデルは発生確率のゾーニングに重点を置くが、体積の推定には十分でない。過去の研究では、リモートセンシングや航空写真、降雨データを用いて体積推定が試みられてきたが、地形・植生・気候要因の複合的影響や非線形性を十分に考慮していない場合が多い。特に、韓国のような地形が複雑で気候変動の影響を受けやすい地域では、より精度の高い体積予測モデルの構築が求められている。

手法
本研究では、韓国を対象に、降雨(誘発因子)と地形・土壌・環境(素因因子)を説明変数とし、浅層地すべり体積を予測するために9種類の機械学習モデルを比較検討した。使用したモデルは、OLS(最小二乗法)、RF(ランダムフォレスト)、SVM(サポートベクターマシン)、EGB(極端勾配ブースティング)、GLM(一般化線形モデル)、DT(決定木)、DNN(深層ニューラルネットワーク)、KNN(k近傍法)、RR(リッジ回帰)である。これらのモデルは、物理ベースモデルと経験的モデルのパラメータを組み合わせ、非線形性を含む複雑な特徴を捉えることを目指す。モデル評価指標には決定係数(R²)、平均絶対誤差(MAE)、二乗平均平方根誤差(RMSE)、平均絶対パーセント誤差(MAPE)、対称平均絶対パーセント誤差(SMAPE)を用いた。データセットは2011~2012年に韓国各地で収集された455件の地すべり記録で、走行長、幅、深さ、体積、地形構成、植生、降雨履歴などの特徴量を含む。特徴量の多重共線性を排除し、モデルの過学習を防ぐために適切な前処理を行った。

結果
モデルの予測性能はモデル間で大きく異なり、OLSやRRは低い決定係数(R²=約0.27~0.30)で最も性能が劣った。一方、EGB、DNN、RFは高い性能を示し、特にEGBはトレーニングセットでR²=0.96、テストセットでR²=0.88を達成した。RFもトレーニング・テストセットでのR²差が小さく安定した予測を示した。DNNはSMAPEが低く、誤差が小さい予測を行った。特徴量の重要度分析では、斜面長(slope length)が最も重要な予測因子であり、次いで降雨量(特に最大1時間降雨量)、斜面角度、斜面方位、標高が重要であった。土壌深さや森林密度の寄与は小さかった。火災履歴は散布図では大きな地すべり体積と関連が示唆されたが、モデルの重要度では有意ではなかった。これらの結果は、地形的特徴と降雨強度が浅層地すべり体積の主要な決定因子であることを示す。

考察
本研究は、物理的・経験的モデルのパラメータを統合し、非線形性を考慮した機械学習モデルによって浅層地すべり体積の予測精度を向上させた。特にEGBやRFは過学習を抑えつつ高精度な予測を実現し、地域特性に基づく体積推定に有効である。斜面長や降雨強度の重要性は、地すべりの発生メカニズムと整合的であり、これらの因子を重視した防災対策が有効と考えられる。火災履歴や土壌深さの影響が小さいことは、地域特性やデータのばらつきによる可能性があるため、さらなるデータ収集とモデル改良が望ましい。気候変動による降雨パターンの変化を踏まえ、継続的なモニタリングとモデルの適応が必要である。今後は他地域への適用可能性や、より多様な環境因子の統合が課題である。

機械学習で崩壊土砂量を推定する内容。相関係数が高すぎるような気がしますが、ま、一例ということで。
斜面長が効くと言われると、そうなのかなとも思います。では、斜面長を作る際にどの程度の面積、凹凸まで許容するか?工夫の必要がありそうです。

Visual C++ 2008 Redist SP1 インストール時エラー

Win11 にソフトをインストール後、起動させようとすると強制終了してしまう問題が発生しました。
以前より見て見ぬふりをしていましたが、Ver.Upしても使えないと限界。対応することにしました。

イベントビューワーでは VC++ の再配布ファイルがないとのこと。

Microsoft.VC90.MFC,processorArchitecture="x86",publicKeyToken="1fc8b3b9a1e18e3b",type="win32",version="9.0.21022.8" が見つかりませんでした。 詳細な診断を行うには sxstrace.exe を実行してください。

sxstrace.exe trace -logfile:sxstrace.etl
(問題のプログラムを実行してエラーを再現)
sxstrace.exe parse -logfile:sxstrace.etl -outfile:sxstrace.txt

吐かれたTXTを見ると、Microsoft.VC90.MFC, version="9.0.21022.8" がないという、同じ文言。古いのでサポートは切れてますが、入れましょう。ということで SP1 をインストール。が、またもエラー。

Visual C++ 2008 Redistributableのインストールで問題発生 - エラー1935
https://www.reddit.com/r/techsupport/comments/qsv2gf/trouble_installing_visual_c_2008_redistributable/?tl=ja

ファイアーウォールに引っかかっていたという初歩的な内容。もっと早く調べるべきでした。

2025年9月3日水曜日

Docking Station WD22

 Dell の Docking Station WD22 に外付け HDD を USB-A で接続すると、本体再起動後にマウントが外れて認識されない症状に遭遇しました。

本体に USB接続していた時は見られなかった症状のため、ハードウェア側の不具合を疑い試行錯誤していたのですが、結局はソフトウェア側(OS側)の設定の影響でした。 

Win11で外付けHDDを休止(スリープ)させないようにすればOK。

  1. スタートメニュー → 「すべてのアプリ」 → 「Windowsツール」 → 「コントロールパネル」
  2. 「ハードウェアとサウンド」→「電源オプション」を選択
  3. 現在使用中の電源プラン右側の「プラン設定の変更」をクリック
  4. 「詳細な電源設定の変更」をクリックし、「ハードディスク」→「次の時間が経過後ハードディスクの電源を切る」の値を「なし(0分)」に設定

これで解決。外部ディスプレーが時々切断される問題も同時に解決しました。よくわかりませんが復帰させるときに Dock が頑張っていたのでしょう。

2025年8月27日水曜日

conda チャンネル書き換え

 miniconda を入れたら、defults チャンネルを削除できなくなっていました。

$ conda config --remove channels defaults


~/miniconda3/.condarc を書き換えるだけでした。

channels:
  - conda-forge


これで conda-forge がデフォになりました。

2025年8月26日火曜日

分位点回帰のピンボール損失

ピンボール損失の意味

分位点回帰では、目的変数の条件付き \(\tau\) 分位点を推定するために、 ピンボール損失(チェック関数)を用いる。二乗誤差が平均を推定するのに対し、 ピンボール損失は任意の分位点 \(\tau \in (0,1)\) に対応する。

残差を \(r_i = y_i - F_i\)(実測−予測)とすると、ピンボール損失は

\[ L(y_i, F_i) = \rho_\tau(r_i) = \begin{cases} \tau \, r_i, & r_i \ge 0 \quad \\[6pt] (\tau - 1)\, r_i = (1-\tau)\,|r_i|, & r_i < 0 \quad \end{cases} \]

非対称な重み付けにより、アンダーとオーバーを異なる強さで罰する。 例えば \(\tau=0.9\) では、下側(予測が小さすぎ)の誤差に強いペナルティがかかり、 予測は上に引き上げられて 90% 分位点に一致しやすくなる。

ピンボール損失(チェック関数)の具体例

\(\tau = 0.8\)、実測 \(y_i=5.0\) に対し、予測 \(F_i\) をいくつか試す。残差 \(r_i=y_i-F_i\) と損失 \(\rho_\tau(r_i)\) は、

予測 \(F_i\) 残差 \(r_i=y_i-F_i\) 判定 ピンボール損失 \(\rho_{0.8}(r_i)\)
4.0 +1.0 アンダー (\(r_i\ge 0\)) \(0.8 \times 1.0 = 0.8\)
5.0 0.0 一致 0
6.0 -1.0 オーバー (\(r_i<0\)) \((0.8-1)\times(-1.0)=0.2\)

同じ誤差量でも、\(\tau=0.8\) ではアンダー(+1.0)の損失 0.8 が、 オーバー(-1.0)の損失 0.2 より大きい=非対称性が働く。

最小化問題の構成

与えられたデータ \(y_i\) に対し、予測値 \(F_i\) を選ぶことで、損失の合計を最小化する。

\[ F^\star = \operatorname{arg\,min}_{F} \sum_{i=1}^n \rho_\tau(y_i - F_i) \]

分位点の定義

ある分布において、\(\tau\) 分位点 \(Q_\tau\) は、以下を満たす。

\[ P(Y \leq Q_\tau) = \tau \]

これは、観測値 \(Y\) に対し分位点 \(Q_\tau\) 以下となる確率が \(\tau\) であることを意味する。

ピンボール損失と分位点の関係性の直感的理解

ピンボール損失は、予測値 \(F_i\) が \(\tau\) 分位点 \(Q_\tau\) から外れた場合に、外れの方向に応じて異なるペナルティを課すことで、分位点の推定を誘導する。 具体的には、以下の性質が働く。

  • 例:\(\tau = 0.8\) を推定したい状況を考える。もし予測値 \(F_i\) が真の 80% 分位点 \(Q_{0.8}\) よりも1だけ小さい場合、\( y_i > F_i \) となるデータ点(残差が正)に対するペナルティは \( 0.8 \) であり、\( y_i \leq F_i \) となるデータ点(残差が負)に対するペナルティ \( 0.2 \) よりも4倍大きくなる。このため、\(F_i\) はより大きな値に引き上げられる傾向がある。
  • 例:逆に、予測値 \(F_i\) が真の 80% 分位点 \(Q_{0.8}\) よりも1だけ大きい場合、\( y_i \leq F_i \) となるデータ点(残差が負)に対するペナルティは \( 0.2 \) であり、\( y_i > F_i \) となるデータ点(残差が正)に対するペナルティ \( 0.8 \) よりも小さくなる。このため、\(F_i\) はより小さな値に引き下げられる傾向がある。

このように、ピンボール損失を最小化する過程で、データ全体として、予測値 \(F_i\) が \(\tau\) 分位点\(Q_\tau\)に近づくように調整される。

2025年8月25日月曜日

GBDT 二値分類の勾配計算例

前提条件

  • ラベル \( y_i \in \{0,1\} \)
  • モデルのスコア(ロジット値)を \( F_i \) とし、予測確率はシグモイド関数で変換:
\[ p_i = \sigma(F_i) = \frac{1}{1 + e^{-F_i}} \]
  • 損失関数(Log Loss):
  • \[ L(y_i, p_i) = -[ y_i \log p_i + (1-y_i) \log (1-p_i) ] \]


損失関数の微分(勾配)

損失関数の勾配 \( g_i \) は L を \( F_i \) で微分して求める。チェーンルールを使うと、

\[ \frac{\partial L}{\partial F_i} = \frac{\partial L}{\partial p_i} \cdot \frac{\partial p_i}{\partial F_i} \]

各項は、

\[ \frac{\partial L}{\partial p_i} = -\frac{y_i}{p_i} + \frac{1-y_i}{1-p_i} \] \[ \frac{\partial p_i}{\partial F_i} = p_i (1 - p_i) \]

よって、

\[ \begin{align*} g_i & = \frac{\partial L(y_i, F_i)}{\partial F_i} \\ & = \left( -\frac{y_i}{p_i} + \frac{1-y_i}{1-p_i} \right) \cdot p_i (1 - p_i) \\ & = -y_i (1 - p_i) + (1 - y_i) p_i \\ & = p_i - y_i \end{align*} \]

つまり、勾配 \( g_i \) は「予測確率と正解ラベルの差」となる。
この勾配の負号(\(-g_i\))が疑似残差として使われ、回帰木の目的変数となる。


直観的解釈

ロジット値 \( F \) はシグモイド関数を通して確率 \( p \) に変換される。
\( F \) を正方向に増やすと予測確率 \( p \) も1に近付く。逆に、\( F \) を負方向に増やすと0に近付く。

  • F = 2 → p ≈ 0.88
  • F = 0 → p = 0.5
  • F = -2 → p ≈ 0.12

「\( p - y \)」という勾配が、「今のモデル予測\( p \)が正解\( y \)より高いか低いか」を表し、Fの上げ下げ・量を直接教えてくれる量となる。

  • p = 0.9, y = 1 → g = -0.1(h=0.1)
  • p = 0.1, y = 0 → g = +0.1(h=-0.1)

\( p < y \) なら「Fを上げる」ことで予測確率も上げる
\( p > y \) なら「Fを下げる」ことで予測確率も下げる

______________________________________________________________

計算例:5データ点による勾配ブースティング(2値分類)

前提条件
  • 入力 \( x_1 \sim x_5 \),ラベル \( y_1 \sim y_5 \)
  • \( y = [1, 0, 1, 0, 1] \)
  • 初期ロジット値 \( F_i^{(0)} = 0 \) (全データ)
  • 学習率 \( \eta = 0.1 \) (計算簡単化のため)
Step 1. 初期化
初期ロジット値
\( F_i^{(0)} = 0 \) (i=1〜5)
初期予測確率
\[ p_i^{(0)} = \sigma(F_i^{(0)}) = \frac{1}{1 + e^{0}} = 0.5 \] 
Step 2. 勾配および疑似残差の計算

i \( y_i \) \( p_i^{(0)} \) \( g_i^{(0)} \) 疑似残差
\( -g_i^{(0)} \)
110.5-0.5+0.5
200.5+0.5-0.5
310.5-0.5+0.5
400.5+0.5-0.5
510.5-0.5+0.5

Step 3. 疑似残差を目的変数に回帰木(弱学習器)学習
  • 単純化して、\( y=1 \) のデータには \( h(x)=+0.5 \)、\( y=0 \) のデータには \( h(x)=-0.5 \) を出力する回帰木を仮定(本来は特徴量で分岐)。

Step 4. モデルの更新
\[ F_i^{(1)} = F_i^{(0)} + \eta h_i^{(0)} \]  
i \( h_i^{(0)} \) \( F_i^{(0)} \) \( F_i^{(1)} \)
1+0.500.05
2-0.50-0.05
3+0.500.05
4-0.50-0.05
5+0.500.05

Step 5. 新しい確率 
\[ p_i^{(1)} = \sigma(F_i^{(1)}) = \frac{1}{1 + e^{-F_i^{(1)}}} \] 
  • \( F = 0.05 \) → \( p \approx 0.5125 \)
  • \( F = -0.05 \) → \( p \approx 0.4875 \)

i \( F_i^{(1)} \) \( p_i^{(1)} \)
10.050.5125
2-0.050.4875
30.050.5125
4-0.050.4875
50.050.5125

Step 6. 新しい疑似残差 
\[ g_i^{(1)} = p_i^{(1)} - y_i \]  
i \( y_i \) \( p_i^{(1)} \) \( g_i^{(1)} \) 疑似残差
110.5125-0.4875+0.4875
200.4875+0.4875-0.4875
310.5125-0.4875+0.4875
400.4875+0.4875-0.4875
510.5125-0.4875+0.4875

Step 7. (以降のステップも同様)
  • この新しい疑似残差を目的変数として、再度回帰木(弱学習器)を学習しロジット値を更新。
\[F_i^{(2)} = F_i^{(1)} + \eta h_i^{(1)}\]
  • \(h_i^{(1)} = +0.4875\)(y=1のデータ)、\(-0.4875\)(y=0のデータ) 

i \( h_i^{(1)} \) \( F_i^{(1)} \) \( F_i^{(2)} \)
1+0.48750.050.09875
2-0.4875-0.05-0.09875
3+0.48750.050.09875
4-0.4875-0.05-0.09875
5+0.48750.050.09875

Step 8. 新しい確率の算出
  • \( F = 0.09875 \) の場合:\( p = \frac{1}{1+e^{-0.09875}} \approx 0.5247 \)
  • \( F = -0.09875 \) の場合:\( p = \frac{1}{1+e^{0.09875}} \approx 0.4753 \)

この流れを何度も繰り返すことで、予測確率 p が徐々に正解ラベルに近付く。
学習率が小さいほど、各ステップの修正は少なくなり安定するが、収束には時間がかかる。

2025年8月18日月曜日

土層区分 + RF

A machine learning-based approach for constructing a 3D apparent geological model using multi-resistivity data | Geoscience Letters | Full Text

AI要約

背景
台湾の濁水渓沖積扇(CRAF)は主要な地下水盆であり、高速鉄道が地盤沈下帯を通過するため、地下の地質構造を詳細に理解することが重要とされている。しかし、従来の地下モデルはボーリングデータに大きく依存しており、その高コストと疎な配置により、広範囲の領域で詳細な空間的リソロジー分布を把握することは困難であった。また、先行研究では抵抗率データとボーリング情報を十分に統合した包括的な3D地下モデルの構築には至っておらず、この不足が課題となっていた。

手法
本研究では、濁水渓沖積扇において、垂直電気探査(VES)、過渡期電磁(TEM)、ボーリング孔比抵抗(NBR)を含む複数の比抵抗データを統合し、3D見かけ地質モデル(Apparent Geological Model (AGM) )を構築する包括的なアプローチを提示。まず、異なる手法で取得された比抵抗データの厳密な調和(ハーモナイゼーション)を行い、整合性のある比抵抗値を確保した。この調和の必要性は、同じ研究エリア内での測定にもかかわらず、各データセットの比抵抗値の範囲(上限・下限)に顕著な不一致が見られたためである。本研究では、以下の手順で厳密な調和を行った。

  1. データ範囲の制限とサンプリングレートの調整: まず、各1次元データセット(VES、TEM、NBR)の深度を200mに制限し、Piecewise Cubic Hermite Interpolating Polynomial (PCHIP) メソッドを用いて、1m間隔でサンプリングレートを調整した。 
  2. 次に、Min–Max Re-scalingとして知られる正規化手法を用いて特徴スケーリングを実行した。これは、元のデータ範囲を線形変換し、所定の境界内で比抵抗値の関係を維持するためである。
  3. VESデータ範囲を基準とした変換: 全ての比抵抗データセットは、研究エリア内のVES比抵抗範囲(1.02 Ωmから2512 Ωm)を基準として変換されました。この「データ取得(data retrieval)」と呼ばれるプロセスにより、比抵抗データはまず0-1の範囲に正規化され、その後sci-analysisパッケージを用いてVES比抵抗値に一致するように再スケーリングされました。VESデータが基準として選ばれた理由は、そのデータがNBRおよびTEMデータと比較して比抵抗値の最低および最高の境界(1.02 Ωmおよび2512 Ωm)を両方とも示し、かつ3つの手法の中で最も密に分布している測定値であった点にある。 

次に、従来のソフトウェアの限界を克服するため、Pythonベースのモデリングと動径基底関数補間(RBFI)を用いて3D抵抗率モデルを構築した。これは、複数のソースから得られた点の情報(抵抗率データ)を3次元空間に補間するために用いられた。抵抗率データは、VES、TEM、NBR。スムージング係数(この研究では500に設定)は、「データ点への正確な適合」と「滑らかな補間表面」の間のトレードオフを制御し、ノイズの多いデータセットやスパースなデータセットへの過学習を防ぐ役割を果たす。

その後、ボーリング情報から得られた堆積物タイプ(粘土、砂、砂利)をグラウンドトゥルースデータとして利用し、教師あり機械学習(SML)技術によって3D AGMに変換した。SMLアルゴリズムとしては、ランダムフォレスト(RF)、決定木(DT)、サポートベクターマシン(SVM)、勾配ブースティング(XGBoost)の4種を実装し、混同行列分析、評価指標、ROC曲線を用いて性能評価を行った。

結果と考察
本手法により、データカバレッジは従来の62地点から386地点へと大幅に増加し、空間カバレッジ密度が約84.02%向上した。SMLアルゴリズムの評価では、ランダムフォレスト(RF)が全ての評価指標において0.952という最高の性能を示した。構築された3D AGMは、堆積物タイプと抵抗率の明確な相関を明らかにした。具体的には、粘土層は低比抵抗(59.98 Ωm以下)、砂層は中比抵抗(59.98 Ωm超136.14 Ωm未満)、砂利層は高比抵抗(136.14 Ωm以上)を示した。また、近位扇状地では砂利層が優勢である一方、中間扇状地では主に砂質粘土層が、遠位扇状地(西部の沿岸地域)では粘土質砂が主体であることも判明した。
マルチ比抵抗データを統合し、データカバレッジを大幅に増加させたことにより、地下の詳細な特性をよりきめ細やかに把握することが可能となった。特に、データハーモナイゼーション技術によって、異なる比抵抗データの範囲を整合させたことは、3Dモデリングの精度向上に寄与した。この研究成果は、濁水渓沖積扇だけでなく、他の地域における3D AGM構築にも応用可能であり、地盤沈下メカニズムの理解や地下水管理、資源探査など、多様な地質学的応用において貴重な洞察を提供する。

複数の物理探査の結果をソースとして、機械学習により土層分布を推定しています。地下水の有無に留意でしょう。RBFNで一度訛ったデータを用いているので、過学習にはなり難いのかもしれません。
 近年、他国で多数報告されている3D地質モデル作成時の機械学習の利用ですが、ジョイントインバージョンに頼らなくて済みますし、扱いが容易なので普及するのでは?と注目しています。国内では3次元地質が多用されているとは未だに言えませんし、物理探査も縮小傾向のまま低空飛行を維持、機械学習も未だに普及していない、ということで、また日本が取り残されていくのを残念に思いながら眺めています。

Geochemistry Data + SOM

Correlations among large igneous provinces related to the West Gondwana breakup: A geochemical database reappraisal of Early Cretaceous plumbing systems - ScienceDirect

AI要約

背景
南大西洋初期開口期に関連する大規模火成岩区(LIPs)は主に苦鉄質(火山性および/または深成性)岩石で構成され、通常は短期間の活動パルスで形成される特徴を有する。白亜紀初期には、ウェストゴンドワナ超大陸が広範な伸張性地殻変動を経験し、現在の南米プレートとアフリカプレートの配置が形成された。この地殻変動の初期段階(134 Ma以前)において、パラナマグマティックプロヴィンス(PMP)とエテンデカマグマティックプロヴィンス(EMP)が形成され、これらはウェストゴンドワナの分散以前には連続したLIPを構成していた。また、赤道大西洋縁辺では、白亜紀初期の伸張性地殻変動により北東ブラジルのリフトシステムとベヌエトラフ/ナイジェリアが発達し、赤道大西洋マグマティックプロヴィンス(EQUAMP)に分類される苦鉄質ダイク群およびシル複合体が貫入した。LIPsの特性評価には、地球化学データベース、地質年代データ、および空間的範囲の決定が不可欠であるが、大量のデータを従来の分析手法で扱うことは困難を伴う。PMPおよびEMPに関する地球化学的研究は1980年代から2000年代にかけて主に玄武岩(CFBs)に焦点を当てて行われてきたが、過去20年間でダイク群およびシル複合体に関するデータが顕著に増加した。本研究では、機械学習ベースのツールである自己組織化マップ(SOMs)を用いて、これらの問題を克服し、PMP、EMP、EQUAMPの貫入性ソレイアイト質マグマ活動の解析を行う。

手法
本研究では、PMP、EMP、EQUAMPの地化学データセットを分析するために、Tiベースの岩石学的分類と自己組織化マップ(SOM)アプローチを組み合わせて適用した。SOMは、入力変数が定義する空間におけるベクトル類似性を用いてパターンと関係性を見出す計算技術であり、複雑な入力を持つ多変量空間データセットの統合分析を可能にする。データは、蛍光X線分析(XRF)による主要酸化物の決定、誘導結合プラズマ質量分析(ICP-MS)、原子発光分光分析(AES)、光学発光分光分析(OES)、発光分光分析(ES)、中性子放射化分析(INAA)などの高精度な手法による微量元素および希土類元素の決定から収集された。解析対象から、MgO > 10.5 wt.% の高マグネシウムサンプル、LOI > 3.5 wt.% の変質が疑われるサンプル、およびSiO2 > 65.5 wt.% の珪質・アルカリ性マグマは除外された。SOMの計算にはSiroSOMパッケージが使用され、58行×52列の計3,016ノードのマップサイズが採用された。SOMの出力は、統一距離行列(U-matrix)やコンポーネントプロットを用いて可視化された。また、地質マッピングには、ブラジル地質調査所(CPRM)およびバイーア鉱物調査会社(CBPM)の航空磁気データ、既存の地質図、およびGoogle Earth画像が組み合わされた。磁気データは、還元極性化(RTP)、解析信号振幅(ASA)、および指向性強調フィルターを用いて処理され、ダイク、剪断帯、広域断層の区別が行われた。さらに、GPlates 2.3.0を用いてプレート再構築モデルが操作され、白亜紀初期のマグマ活動と南大西洋関連ホットスポットの古地理学的特徴が視覚化された。

結果
全アルカリ-シリカ(TAS)分類図において、PMP、EMP、EQUAMPのHT(TiO2 > 2 wt.%)スイートは、SiO2 < 56 wt.%、総アルカリ < 6.8 wt.% を示し、宮下線の下にプロットされ、準アルカリ性玄武岩、玄武岩質安山岩、および玄武岩質粗面安山岩に分類された。MgO対TiO2図では、HT群はMgO > 4 wt.% のサンプルにおいてTiO2 = 0.1 MgO + 3.3 の線より上にクラスターを形成し、MgO < 4 wt.% のサンプルではTiO2 = 0.725 MgO の曲線によって他のスイートと区別された。Ti/YパラメーターはEQUAMPソレイアイトにおいて250から325までの連続的なトレンドを示し、高チタンマグマと低チタンマグマの重複が見られた。不適合元素(例:Sr > 400 ppm)が多く、La/Yb > 10 の特徴も示された。進化した岩石(ERs)はSiO2レベルが57~65 wt.%、TiO2 < 2.5 wt.% であり、主に粗面安山岩~粗面岩の領域に限定された。ERはTi/Y対Ti/Zr図においてLTマグマと同様のTi/Y値を示すが、Ti/Zr値は低い。低チタン(LT、TiO2 < 2 wt.%)マグマはHT岩石よりも低いアルカリ/シリカ比を持ち、準アルカリ性玄武岩および玄武岩質安山岩で構成される。PMPとEMPに特徴的なTiO2の遷移的なソレイアイト(TTs)は、TiO2レベルが2.7~1.7 wt.% のサンプルであり、HTとLTの中間的な挙動を示した。

考察
SOMの使用は、専門的な岩石学的知識や既に確立された地球化学的分類を置き換えることを意図するものではなく、むしろデータ集約的な地域比較を行うことを目的としている。SOM手法においては、すべての変数に均等な重みを与えることが可能である。EQUAMPデータベースの検証においては、中央大西洋LIPに属する可能性のあるサンプルを認識し、除去することが必須のステップであった。しかしながら、ナイジェリア地域における初期白亜紀ソレイアイトに関する地球化学的および地質学的情報が不足しているため、赤道大西洋の両側間の正確な相関関係を導き出すことは困難である。

新規性
PMP、EMP、EQUAMPの地化学データセットに対して、Tiベースの岩石学的分類とSOMアプローチを組み合わせて分析を行った点にある。特に、貫入性のダイク群およびシル複合体のデータに焦点を当て、主要酸化物、微量元素、希土類元素、および同位体データ(Sr、Nd、Pb)の定性的な再評価と組み合わせることで、初期南大西洋開口期に関連する苦鉄質ソレイアイト質配管システムを特徴付けた。SOMのような機械学習ベースのツールを用いることで、大量のデータを統合的に分析し、複雑な入力を持つ多変量空間データセットからパターンと関係性を抽出することで、従来の分析では困難であった微細で複雑な関係性の理解を深めることが可能になった点も、本研究の重要な新規性である。

多量の地球化学分析データをグルーピングするのにSiroSOMを使用しています。水質や元素などの分析データをグルーピング、視覚化するのに使えるツールなのでしょう。



水質分析 PCA + SOM


AI要約
背景
北フィンランドにおけるSakatti鉱鉱開発サイトは、複雑な更新世の堆積物と変質した・破砕された基盤岩の領域に位置する。この地域には広範な泥炭地と、地下水、表流水、湿地または泥炭地間の相互作用に関する理解が求められている。具体的には、アーパ泥炭地と呼ばれるパターン化したフェンが存在し、これらは栄養が乏しいが、時折栄養の豊富な泥炭地も見られる。地下水流動の理解は、鉱業活動を計画する上で重要な要素である。

手法
データ収集:
二つの水文化学データセットが使用された。一つはヘルシンキ大学から、もう一つはAA Sakatti Mining Oyから得られた。両者のデータは2015年から2020年の期間として収集された。

データ分析:
  • 主成分分析(PCA): 化学変数を分析し、相関関係を明らかにするために実施された。PCAの結果、主成分はEC、Ca、Mg、HCO₃が含まれる側面と、Al、Cr、Fe、Mn、Co、Ni、Znが含まれる側面の二つに分類された。
  • 自己組織化マップ(SiroSOM): 水質パラメータ間の関連性を理解するために使用され、異なる水のグループが明確にされ、7つの水文化学的クラスタに分類された。
  • SiroSOMはCSIROが開発した教師なしデータマイニング手法で、類似した水質パラメータを持つ水サンプルをクラスタリングする。これは自己組織化マップ(Self-Organizing Maps, SOM)に基づき、最適マッチングユニット(BMU)のベクトルを用いて変数の類似性を2Dマップ上に投影する。
  • 分析に用いたパラメータはpH、電気伝導度(EC)、主要イオン(Na, Ca, K, Mg, Cl, SO₄, HCO₃)、および微量元素(Al, Cr, Co, Ni, Zn)。クラスタは2DシートマップとLeapfrog Geoソフトウェアによる3Dビューで可視化。
地下水流動モデリング:
  • MODFLOW: 事前に構築された3D地質モデルGM2020に基づき、流動のシミュレーションを行った。二つのモデル(PRE1989とPOST2014)を使用し、地域の流れの違いを理解するために、定常状態モデルと過渡状態モデルを計算した。
  • 粒子追跡(MODPATH): 水サンプルの再充填エリアを特定するために、バックス追跡(BW追跡)が行われた。水サンプルの各地点で粒子をトラッキングし、地下水の流動経路を可視化した。
安定同位体分析:
水サンプルにおける異常気象の影響を考慮し、安定同位体のデータを使用して地下水の流入エリアを解析した。

結果
データ分析結果:
SiroSOMにより、7つのクラスタが得られた。
  • SOM 1: 主要な陽イオン、EC、HCO₃、Mn、Feが高く、pHが比較的高い地下水および間隙水(1つの湧水を含む)
  • SOM 2: 主にEC、Mn、Feが高く、CoとZnも高い泥炭地の間隙水(1つの地下水、1つの地表水を含む)
  • SOM 3: Na、K、EC、SO₄、Mn、Niが高く、CoとZnがわずかに高い地下水サンプル 
  • SOM 4: EC、Co、Ni、Zn、Al、Crが高く、MnとFeも高い水(地表水、地下水、泥炭地の間隙水を含む)
  • SOM 5: SO₄がわずかに高く、MnとFeが低い、比較的希釈された地表水 
  • SOM 6: SO₄が高く、pHが比較的高い地表水(1つの地下水、1つの湧水を含む)
  • SOM 7: 最も希釈された地表水(3つの泥炭地の間隙水、1つの湧水を含む)
流動パターン:
水の流動は大規模(>1000 m)、中規模(1000–100 m)、小規模(<100 m)の三つのスケールで確認された。
  • 大規模流路は、泥炭島から始まり、残存した堆積物を通過して、最終的にはキティネン川に至る。
  • 中規模流路は、堆積平野や隣接する指状泥炭地域から地下を経由して川に流れ込む。
  • 小規模流はキティネン川の岸辺やViiankiaapa泥炭地で見られ、地形の変化に関連している。
  • 平均滞留時間: 粒子追跡による平均滞留時間は、大規模流路で約150年、中規模で25年、小規模で6年であった。これにより、流動の特性が明確に示された。
安定同位体データ:
d-excessの分布から、蒸発した水信号(<5)や降水由来の水信号(>8)、混合水信号が確認された。特に、Viiankiaapa泥炭地域では地下水の再充填が確認されており、その水がキティネン川周辺で観察される浅い井戸や湧き水に流入することが示唆された。

水質特性:
大規模流路の水はCa-HCO₃型であることが判明し、SiroSOM分析により、地下水、富栄養な水、希薄な表流水などの水サンプルが異なる水文化学的クラスタに分類された。流れを通じて水の化学的特徴が変化し、特に泥炭層を通過することで多くの化学変数の平均濃度が低下することが示された。

考察
SiroSOMとPCA(主成分分析)の両方の分析結果は、溶解元素濃度が高い3種類の水タイプが存在することを示した。
  • 第1グループ: PCA 1とSOM 1に対応し、EC、Ca、Mg、HCO₃、pHの上昇が特徴
  • 第2グループ: PCA 2とSOM 4、2、3に対応し、Al、Mn、Feが高く、pHが第一のグループ(平均7.2)よりも低い(平均6.3)のが特徴。Cr、Co、Ni、Znといった微量元素も高く、これらは苦鉄質鉱物または苦鉄質硫化物の溶解に関連している。
  • 第3グループ: PCA 3とSOM 6、5、3に対応し、NaとSO₄が高く、K、Cl、Ni上昇
研究成果は、鉱業開発地における地下水のリチャージ領域や流出領域を特定することの重要性を示す。地下水模型と安定同位体分布は、流動経路の検出に適用可能であり、地下水管理や環境モニタリングの計画に役立つ。各水質クラスタとそれらの化学的特徴は、地区特有の水文学的条件に根ざした結果であることが示された。

SiroSOMは商用ツールだそうです。知りませんでした。PCAのみでなく、他の手法と合わせてグループ分割できたという点は実務でも利用価値があります。

2025年8月17日日曜日

GeoChemFoam 適用例

Integrated numerical simulation of reactive bubble columns: A coupled approach with openFOAM and PhreeqC - ScienceDirect

AI要約

1. 背景
気泡塔リアクターは、化学、石油化学、生化学プロセスなど産業界で広く利用される多相流システムである。製造が容易で保守費用が低い利点がある一方で、逆混合が主な欠点とされている。地熱発電所の生産パイプも、水とガス気泡の共存という特性から気泡塔と見なすことができ、その動的な相互作用が重要となる。
しかし、従来の数値シミュレーション手法(Volume of Fluid, Eulerian–Eulerian, Eulerian–Lagrangian)は、それぞれ計算コストや個々の界面の追跡、大規模シミュレーションへの適用性などに課題を抱えていた。特に、気泡塔内の反応をモデル化する際には、簡略化された化学反応セットが用いられることが多く、実世界のHydrogeochemistryにおける複雑な反応を完全に捉えることが困難であった。
HydrogeochemicalモデルであるPhreeqCは、広範な熱力学データベースを持ち、複雑な化学反応(固相、ガス-水相互作用、酸化還元反応など)を扱うことが可能であるが、単独では多相流動や2次元以上の輸送プロセスを扱うことができない。また、既存のOpenFOAMとPhreeqCのカップリング手法(porousMedia4Foam, GeoChemFoam)は、これまでのところEulerian–EulerianやEulerian–Lagrangianといった多相流モデルには対応していない。
本研究は、このような背景から、多相流モデルと地化学モデルを統合し、反応性気泡塔のより包括的で現実的なシミュレーションフレームワークを開発することを目的としている。

2. 手法
本研究では、OpenFOAM (v2006) の「reactingTwoPhaseEulerFoam」ソルバーと、PhreeqCRMコードおよびGeoChemFoam (v-20-4) パッケージを組み合わせたEulerian–Eulerian多相流アプローチを採用した。

  • 流体力学モデル:各相(気体、液体)について、連続の式と運動量保存式を解く。界面力(interfacial force)として、ドラッグ力(Ishii and Zuberモデル)、リフト力(Tomiyamaモデル)、仮想質量、壁潤滑、乱流分散を考慮。
  • 化学種輸送モデル:化学吸着プロセスは輸送方程式を用いて実装され、有効拡散係数には乱流粘性も含まれる。
  • 物質移動:ガス相から液相への物質移動は、2-film theoryに基づいてモデル化。物質移動係数に加え、化学反応の影響を考慮するエンハンスメントファクターが適用される。気泡ダイナミクスが物質移動に与える影響をモデル化するために、動的Sherwood数も利用。
  • 化学反応:液相での化学反応は、OpenFOAM-PhreeqCカップリングを介してPhreeqCRMによって計算される。これにより、簡略化された反応セットではなく、炭酸系に加えてカルシウム-炭酸系、マグネシウム-炭酸系、ナトリウム-炭酸系など、PhreeqCの熱力学データベースで利用可能な広範な水地化学反応を統合的に扱うことが可能。
  • 計算設定:差分方程式の離散化スキームとして、ddtSchemesにEuler、gradSchemesにGauss linear、divSchemesにGauss upwind、laplacianSchemesにGauss linear correctedがを利用。

3. 結果
CO2を高pH溶液に注入するケース(Darmana et al. [19]の確立された試験ケース)を用いて検証。

  • CO2速度と流動特性:気泡塔の軸方向におけるCO2速度のシミュレーション結果は、実験データと良好な一致を示した。ガス体積分率は底部で最も高く、高さとともに減少し、S字型の軌道を示す。気体速度は底部で最も高く、上方へ減少する傾向が見られた。
  • 化学種濃度とpH値:CO2-およびCO3-のシミュレーション濃度は、測定値と比較して低い値を示した。これは、総無機炭素量が実験よりも少ないか、物質移動率が過小評価されている可能性を示唆している。pH値の経時変化では、シミュレーションは実験よりも低下が速く、最終的なpHレベルも低い結果となった。これは、シミュレーションにおける物質移動率の過大評価、または実験と比較して気泡径の差異がある可能性を示唆している。
  • パラメータの影響:見かけのガス速度が増加すると、液体のpHはより速く低下。ガス入口面積が大きいほど、pHの低下が速まることが示された。
実世界例(ミュンヘンの水道水へのガス注入):
  • CO2注入: CO2(aq)濃度は時間の経過とともに増加し、pHは低下。これにより、溶液は方解石とアラゴナイトに対して過飽和状態(SI値が0を超える)となり、平衡状態に達することが示された。
  • 空気注入: CO2ストリッピングによる鉱物沈殿の可能性を評価するために、空気注入シナリオも分析され、種濃度の変化が示された。

4. 考察
本研究で開発されたモデルは、多相流動特性と化学反応の複雑な相互作用を効果的に捉える能力を示した。気泡のダイナミクスが物質移動に与える影響や、物理的・幾何学的パラメータがpH値に与える影響を研究することで、様々な条件下での気泡塔リアクターの挙動理解が深まった。
ただし、シミュレーション結果と実験データの比較では、化学種濃度やpHの経時変化において一部差異が見られた。これは、シミュレーションにおける物質移動率の評価、またはモデルが気泡のブレークアップや合一効果を考慮していないことによる気泡径分布の非現実性などが影響している可能性がある。
PhreeqCとのカップリングにより、従来の簡略化された化学反応セットでは不可能だった、実世界の水溶液における広範な水地化学反応を統合的にモデル化できるようになったことは大きな進歩である。このフレームワークは、地熱発電のようにhydrochemical reactionが技術的な安全性や経済的成功に重要な応用分野において、貴重なツールとなることが期待される。

5. 新規性
研究の主な新規性は以下の点にある。

  • OpenFOAMの流体力学モデル(Eulerian–Eulerian二流体モデル)とPhreeqCの地化学モデルを統合した、新しいカップリングアプローチを提示した。これは、既存のOpenFOAM-PhreeqCカップリングが対応していなかった、多相流モデルとの組み合わせを実現したものである。
  • PhreeqCの広範な熱力学データベースを活用することで、炭酸系だけでなく、カルシウム-炭酸系、マグネシウム-炭酸系、ナトリウム-炭酸系など、実世界の複雑なhydrochemical reactionを包括的にモデルに組み込むことを可能にした。従来の簡略化された反応セットを用いた研究との差別化が図られる。これにより、地熱・電力生産など、hydrochemical reactionが技術的安全性と経済的成功に不可欠な分野への応用可能なフレームワークを提供する。

GeoChemFoam の適用例です。固相との反応は使われていませんが、容易なのでしょうね。

GeoChemFoam

GeoChemFoam: Direct Modelling of Multiphase Reactive Transport in Real Pore Geometries with Equilibrium Reactions | Transport in Porous Media 

AI要約

GeoChemFoamは、OpenFOAMをベースとしたオープンソースのツールボックスであり、多相・多成分反応輸送現象を、実際の細孔構造やマイクロCT画像データ上で直接数値シミュレーションできるソルバーを特徴としている。Geochemical reactionは、流体バルクと固体表面における平衡反応を共に考慮し、流体力学計算はVolume-Of-Fluid(VOF)法、反応計算は米国地質調査所のPhreeqcを使用、統合している。

主な特徴と方法

  • 流体流れ: VOF法による二相(例えば水とガスや油)の界面追跡。インターフェイスの移動やトリプルライン(固体・流体1・流体2の境界)の接触角もモデル化可。
  • 物質輸送: Continuous Species Transfer (CST)法により、多成分の輸送と界面移動を高精度に解く。CST法は界面付近での人工的な物質移動を抑制し、物質拡散、対流、界面移動をシームレスに扱う。
  • 化学反応: 反応はPhreeqcにより逐次的スプリッティング法(transports→reaction)で計算され、バルクおよび表面反応ともに平衡条件を満たす。
  • モデルの検証: 1次元や2次元の単純な系で解析解と比較検証がなされ、さらに3D細孔内CO₂気泡の溶解・微細な砂岩ポア内のCaCl₂溶液注入など、複雑な実例にも適用している。
  • アップスケーリング: 細孔スケールモデルの結果を、より大きなマクロスケール(Darcyスケール)反応輸送モデルへパラメータとしてフィードバック(界面での実効的移動係数や反応速度定数など)。

主な成果と結論

  • 従来モデルの問題点(界面条件の粗雑な扱いによる人工的な物質移動や誤った反応)はCST法により大幅に改善。
  • CO₂バブル溶解やCaCl₂/油システムなど、実在細孔空間での反応輸送を「直接」数値的に再現可能。
  • 細孔スケールでの詳細シミュレーション結果を使い、マクロモデルの「線形転送モデル」や「混合誘起反応速度モデル」にパラメータ化する手法を提案・検証。
  • 表面反応の場合は、単に平衡速度を修正するだけでなく、表面電荷の補正(例:0.95倍)などが精度向上に有効であった。
  • 全ソースコード・データ等は公開されている。

エアを含む流れの計算が得意そうです。porousMedia4Foamとは異なりハイブリッドではないので、微細な構造に特化した計算に向いているのでしょうか。CT画像をシミュレーションに流用できるのはユニークなアプローチです。

淡水pCO₂の決定法

Determining freshwater pCO2 based on geochemical calculation and modelling using PHREEQC - ScienceDirect

AI要約

背景
化石燃料の燃焼により大気中の二酸化炭素(CO₂)濃度が上昇し、これが地球の気候や海洋に影響を及ぼしている。最新の知見によると、大気中のCO₂レベルの上昇はfreshwaterのCO₂分圧(pCO₂)にも影響を与え、これは多くのパラメータによって制御される。また、これまでのところ、現在および過去の淡水pCO₂レベルを測定する標準化された方法は存在しない。

手法
本研究では、PHREEQCという地球化学プログラムを用いて、フィールド、実験室、長期データに基づいたpCO₂の分析方法を記述。PHREEQCは、炭酸の解離定数や、全無機炭素(TIC)、全CO₂濃度(TCO₂)、全アルカリ度(TA)などのさまざまな入力パラメータを使用して、淡水pCO₂を計算・モデル化することができる。このプログラムは、複数の水サンプルの品質と信頼性を評価し、すべての関連する水生複合体の自動計算を可能にする。

結果
実験室データに基づくPHREEQCによるpCO₂計算の結果は、pH、TCO₂、TAといったパラメータに基づくモデル化結果と比較した場合、わずかな乖離が見られた。また、四つの淡水貯水池における長期モニタリングデータでも、PHREEQCによるpCO₂計算は一貫して信頼性が高く、過去のデータとの比較も可能であった。

考察
淡水系におけるpCO₂の変動は多くの環境要因に影響され、これらの要因を正確にモデル化することは容易ではない。PHREEQCを使用することで、過去のpCO₂の変動や気候変動による影響を体系的に評価することができる。特に、有機CO₂種を考慮することが重要であり、今後の研究においてもこのソフトウェアの活用が期待される。

PHREEQCの適用例です。計算は単純で、input例 を先に読んだ方が内容を解釈しやすいと思います。

Temperature °C temp 22.10
pH–value pH 6.406
Density mg kg-1 density 1.00
total CO2 concentration
(TCO2) mg l-1 C(4) # or C 78.66 as CO2 21.45 as C
Alkalinity (TA) mg l-1 Alkalinity 46.34
Calcium mg l-1 Ca 80.09
Chloride mg l-1 Cl 144.46
Magnesium mg l-1 Mg 12.16
Sulphate mg l-1 S(6) 48.06
Potassium mg l-1 K 3.144
Sodium mg l-1 Na 18.62
Oxygen mg l-1 O(0) 37.02
Silicium mg l-1 Si 0.494
Nitrate mg l-1 N(+5) 0.199
Phosphor mg l-1 P 0.150

 EQUILIBRIUM_PHASES Calcite CO2(g)

Hydrogeochemical evolution

Hydrochemical evolution of groundwater in a semi-arid environment verified through natural tracer and geochemical modelling, northwest Australia - ScienceDirect

AI要約

背景
オーストラリア北西部の半乾燥地域に位置するハマーズリー盆地では、採掘作業のために大量の地下水の再配置が必要とされるため、水収支成分の理解が水の管理と保全に不可欠である。地下水の化学組成と進化は、再チャージプロセス、地下水滞留時間、地下の鉱物学によって主に制御されるが、乾燥から半乾燥地域では、降雨量に対する高い蒸発率がこれらのプロセスと複雑に絡み合う。特に、垂直流が浅い地下水流動を支配する小規模で局所的な帯水層における水理化学的進化に関する研究がこれまで不足している。本研究は溶存イオンや水分子の安定同位体などの天然トレーサーを用いて、局所規模の沖積平野帯水層の水文地質モデルを開発することを主要な目的とした。

手法
本研究では、同位体分析と水理化学的手法を組み合わせたアプローチが採用された。具体的には、溶存イオン、水分子の安定同位体(δ2Hおよびδ18O)、塩化物(Cl)、臭化物(Br)、硫黄同位体(δ34S)といった天然トレーサーが使用された。主要および微量イオン濃度とこれらのトレーサーを組み合わせることで、鉱物溶解反応を特定する方針とした。
フォワード地球化学的モデリングには、PHREEQCソフトウェアを用いた。降雨データは、SILOデータベースおよび近隣の2つの気象観測所から取得した。
地下水再チャージ率は、塩化物マスバランス法(R = (Cp / Cr) * P)によって計算した。この手法では、サイクロン水には蒸発係数fRo = 0.995、低降雨イベントにはfRo = 0.96が適用され、様々な降水イベントからの水の混合がシミュレートされた。

結果
観測された塩化物濃度と安定同位体データは、高降雨サイクロンイベントからの再チャージが、高度に蒸発した低降雨イベントからの再チャージと約60:1の比率で混合していることを示した。
地下水は、高いアルカリ度と溶存酸素を特徴とする独特の水理化学的シグネチャを有していた。総溶存固形分(TDS)は淡水から汽水の範囲であるが、ほとんどの地下水はTDSが1000 mg/L未満であった。δ2Hとδ18O濃度は、広範囲のCl濃度にもかかわらず、狭い範囲で変動する傾向が見られた。
主要イオン分布の水理化学的経路モデリングは、地下水が不飽和帯での再チャージ前に降雨の蒸発濃縮によって進化したことを示し、その後、溶存CO2が増加し、炭酸塩鉱物が沈殿したことが確認された。臭化物と塩化物のプロットは、R2=0.93で海水希釈ライン上に濃度が集中し、両イオンが海洋起源(降雨由来)であることを示唆した。しかし、ほとんどの主要イオンは海水希釈ラインよりも高い濃度を示し、水-岩石相互作用が蒸発濃縮以外の主要イオン濃度の増加に寄与していることを示した。
δ34SSO4とSO4/Cl比の負の相関は、比較的枯渇したδ34SSO4値を持つ硫化物が地下水に追加されたことを示唆している。

PHREEQC
1. 化学種の進化 (Fig. 10a)
シミュレーションの各シナリオを通じて、Cl、SO4、pH、方解石飽和指数(SIcalcite)、CO2(gas)のlog10分圧、および方解石沈殿物量の変化が示されている。

  • 初期段階(シナリオI、II): 蒸発濃縮が進行すると、水はアルカリ性のpHを獲得し、方解石が飽和し沈殿。炭酸塩鉱物(方解石、ドロマイト、マグネサイト、セピオライトなど)は、常に最初に沈殿する相であり、地下水進化の重要な分岐点となる。
  • 土壌CO2との平衡化(シナリオIII): 土壌CO2との平衡化により、pHは8.3に低下し、pHと固定されたCO2分圧に応じてアルカリ度とCa濃度が高く調整される。
  • イオン交換(シナリオIV): イオン交換により、主要な陽イオン濃度がさらに変化。観測された中央値地下水濃度と比較して、Naとアルカリ度は5%以内、MgとCa濃度は10%以内で再現。
  • 最終段階(シナリオV): ケイ酸塩と少量のSO4の溶解を考慮することで、シミュレートされた水質(Cl、pH、Ca、アルカリ度、Mg、SO4)は、現場の観測された中央値組成(Fig. 10b)を非常に密接に模倣できた。最終的な水のpHは8.13となり、観測された中央値pH8.2に非常に近い値を示した。

2. 水質進化のメカニズムに関する洞察

  • 蒸発濃縮: 涵養前の不飽和帯における降雨の蒸発濃縮が初期の地下水化学的進化に重要な役割を果たす。
  • CO2の取り込みと炭酸塩沈殿: 有機物の分解による初期の溶解CO2の増加に続き、帯水層内のCO2分圧が大気レベルより10倍高い(pCO2~2.5)ことは、アルミノケイ酸塩の風化によるさらなるアルカリ度の追加を示唆しており、炭酸塩鉱物の沈殿に繋がる。
  • 水-岩石相互作用: 蒸発濃縮された涵養水は、比較的短時間で水-岩石相互作用によって化学組成が変化する。これには、炭酸塩鉱物の平衡化、ケイ酸塩鉱物の加水分解、イオン交換反応が含まれる。
  • 飽和状態: 帯水層の地下水は、アラゴナイト、方解石、ドロマイト、マグネサイト、セピオライトなどの炭酸塩鉱物に対して飽和または過飽和に達したが、石膏や黄鉄鉱に対しては不飽和である。
3. 水文地質学的概念モデルへの貢献
  • 分析的な水化学モデルと数値モデルは、蒸発した涵養水が水-岩石相互作用によって変化するという概念モデルを支持している。この相互作用は比較的短い時間間隔で発生し、水化学は確立される。
  • 主要イオン化学の深度プロファイルは、涵養水が短期間かつ浅い深度で明確な地下水化学を獲得することを示しており、水化学は深度に対して著しく均質であることも示されている。これは、浅層地下水流レジームにおいて垂直流が支配的であり、それが浅層地下水化学の進化に影響を与えているという氾濫原帯水層の水文地質学的概念モデルを支持。

考察
高度に蒸発した低降雨イベント水は、高降雨イベントからの降水と混ざり、下層の帯水層を再チャージするまで土壌プロファイルに留まるというメカニズムが示された。
硫化物の供給源については、高ヒ素黄鉄鉱濃度を特徴とする基盤岩からの黄鉄鉱の酸化である可能性が高い。
帯水層マトリックス中の高い炭素濃度、特に地下水面付近の濃度は、掘削中に観察された炭酸塩沈殿物(2m~4m厚)と関連付けられる。炭酸塩鉱物の沈殿は、地下水進化経路において重要な分岐点となる。炭酸塩沈殿と溶解には、再チャージ水と浅い地下水の蒸発濃縮、土壌中のCO2吸収、異なるPCO2レベルを持つ水の混合、CO2の脱ガスといった複数のプロセスが関与している。

この資料からはPHREEQCが1次元輸送かバッチのどちらを採用したか特定することはできません。が、記述されているモデルの性質はバッチ計算(またはそれに近い逐次的な平衡計算)で実施されたのでしょう。

気候変動に応じたLandslideの増加傾向

Climate change increases the number of landslides at the juncture of the Alpine, Pannonian and Mediterranean regions | Scientific Reports

AI要約

背景

  • 今後数十年間で、降雨の頻度と規模の変化が景観の進化、社会、経済に大きな影響を与えると予想されている。
  • 地すべりは、地球の表面を形成する要因の一つであり、その安定性は降雨に大きく影響される。地球温暖化による降雨事象の頻度と強度の変化は、地すべりの頻度、量、分布を変化させると考えられている。
  • スロベニアは地すべりの影響を受けやすい国であり、アルプス、パンノニア平野、ディナル山脈、地中海地域の接点に位置し、3つの異なる気候帯を持っている。
  • 浅い地すべり(shallow landslides)は主に降雨のピークと降雨強度に、深層地すべり(deep-seated landslides)は長期(月単位または季節単位)の積算降雨と関連する地下水変動に左右されるため、気候変動に対する反応が異なると考えられている。
  • スロベニアでは、国土の3分の1が地すべりに対して非常に脆弱であり、人口の約5分の1が地すべり危険地域に住んでいる。最も一般的な現象は浅い地すべりであり、主に激しい短期間または長期間の降雨イベントによって引き起こされる。

手法
対象地域と期間: スロベニア全域を対象とし、21世紀末までの将来の季節的降雨変動が地すべり発生に与える影響に焦点を当てている。

気候シナリオ: CMIP5地球気候シミュレーションに基づく、中程度の気候シナリオ(RCP4.5)と最悪の気候シナリオ(RCP8.5)を使用。

データ: EURO-CORDEXプロジェクトの6つの地域気候モデル(RCMs)と6つの地球気候モデル(GCMs)の組み合わせから、日降雨データセットを使用。

空間分解能: 元の12.5km分解能のデータを1kmに変換。

期間設定:ベースライン期間: 1981-2010年

将来予測期間: 

  • 1st period (近未来) 2011-2040年
  • 2nd period (世紀半ば) 2041-2070年
  • 3rd period (世紀末) 2071-2100年

分析アルゴリズム:

  • 浅い地すべりの危険地域算出にはスロベニアの国家地すべり予測システムであるMASPREMアルゴリズムを使用。
  • 深層地すべりの代理指標として、mGrovaモデルから導出される地下水涵養(groundwater recharge)を考慮。

評価項目:

  • 将来の季節的な降雨変動が地すべり危険地域にどのように影響するか(浅い地すべり、深層地すべりの誘発要因としての地下水涵養)。
  • 将来予想される地すべり警報の数。
  • 予想される気候変動が土地被覆の変化に与える影響。
  • 地すべり危険度クラス(低、中、高)を考慮し、50パーセンタイル値をプロットしています。

結果
季節ごとの地すべり危険地域(浅い地すべり)の変動:

  • RCP4.5シナリオでは、冬、夏、秋に危険地域が近未来から世紀末にかけて徐々に5%から12%増加。春は世紀半ばより世紀末の方が低くなる傾向がある。(ページ3)
  • RCP8.5シナリオでは、冬に最も大きな変化が予想され、ベースライン期間と比較して世紀半ばから世紀末にかけて最大12%の危険地域増加が見込まれる。春には世紀末までに国土の約10%が危険地域になる可能性がある。夏と秋も世紀半ばで大きな変動が見られる。

地下水涵養(深層地すべり)の変化: 両方のRCPシナリオで、冬に最も顕著な変化が見られ、地下水涵養が6%から8%増加。他の季節では、地下水涵養の変化はわずかであった。

浅い地すべりの将来の警報頻度:

  • RCP8.5シナリオでは、中・高地すべり危険度地域で、すべての季節において警報が大幅に増加。特に夏と秋の高地すべり危険度地域では、世紀半ばから世紀末にかけて200回以上の警報が予想される。
  • RCP4.5シナリオでも警報は増加するが、RCP8.5ほどではない。夏季に警報が最も多く予測される。

警報の空間分布:

  • RCP4.5では、国内の北西部、中央部、東部、北東部で最も多くの警報が発生。警報数は現在から世紀半ば、世紀末にかけて増加し、春と秋に最大の警報が予想される。北西部(アルプス地域)では冬季に最多警報が予想される。
  • RCP8.5でも同様の傾向が見られるが、より顕著であり、地すべり警報はより広い地域に影響を与える。

地表被覆タイプごとの影響:

  • 浅い地すべりでは、冬、春、夏、特に秋にブドウ畑、牧草地、果樹園が主に影響を受ける(20〜40%)。森林、市街地、耕作地への影響は少ない。
  • 深層地すべりでは、冬と秋に牧草地、森林、市街地が主に影響を受ける(20〜40%)。春と夏は地下水位の変化が最も少ないため、影響を受ける地域も少なくなる。
  • 浅い地すべりは春と夏に深層地すべりよりも景観への影響が大きいと予想される。

考察

  • 本研究で用いられたCMIP5の地球気候シミュレーションは、中央ヨーロッパと地中海地域で統計的に有意であることが示されている。IPCCの報告書とも整合性が見られる。
  • スロベニアが属する中央ヨーロッパおよび地中海ヨーロッパでは、夏季の降雨が減少し、激しい嵐や短時間の集中豪雨といった極端な事象が増加している。
  • 地すべり危険地域は冬に最も多く、次に夏に多いと予想されるのは、冬の降雨量増加と夏の短時間強雨に起因すると考えられる。
  • MASPREMアルゴリズムは、夏の激しい降雨イベント(短時間で高強度)に強く関連している可能性があり、これにより、冬の降雨増加が見られるにもかかわらず、夏の警報頻度が高いという結果が説明され得る。
  • 研究は気候モデルとMASPREM/mGrovaアルゴリズムの限界を認識しており、特に深層地すべりの誘発が地域の水文地質学的および構造的条件に大きく依存し、MASPREMでは考慮されていない点を指摘。
  • 分析には、シナリオ駆動型気候予測に固有のエラーと、水文学的モデリングの認識論的(Epistemic)不確実性が伴う。しかし、地域水文学的分析はモデリングの不確実性を減らし、広範囲での地すべり活動の変動を推定するのに有用である。
  • 浅い地すべりは主にブドウ畑、牧草地、果樹園で発生し、スロベニアのブドウ畑がテラス状に整備されており、人為的な影響を受けているため、地すべりが発生しやすいことが示唆される。
  • 深層地すべりは主に牧草地、森林、市街地で発生し、その発生は深層地質学的条件、局所的な水文地質学的条件、地盤の地力学的特性に強く影響される。
  • 地域規模での定量的な分析を行い、政策決定者が適応策を設計する上で有用な情報を提供している。
国交省の「気候変動を踏まえた砂防技術検討会令和5年度版とりまとめ」では、梅雨時期の豪雨発生北限拡大、豪雨数~1.5倍、線状降水帯発生数~1.6倍と言われています。
緯度が異なりますが、日本でも災害発生数などは文献に書かれている傾向に近づくのでしょう。


岩盤表層の亀裂発達

 Explorative data analysis from multiparametric monitoring at the Acuto Field Laboratory (Central Italy) for detecting preparatory conditions to rock block instabilities | Italian journal of engineering geology and environment

AI要約

背景
本研究はイタリアのアクートフィールドラボ(Acuto Field Laboratory)で行われ、ロックブロックの不安定性に関連する条件を検出することを目的としている。研究対象の岩壁は、地元の気候条件や自然要因により、徐々に不安定になる可能性がある。特に、季節的な温度変動や降雨、風などの環境要因が岩の変形に与える影響が重視されており、これらの要因が崩落に至る過程の理解に寄与すると考えられる。

手法
アクートフィールドラボにおいて、複数のパラメトリックモニタリングシステムを用いて、環境ストレッサー(気温、降雨、風速など)とそれによる岩の変形を記録。モニタリングは2015年から継続され、様々なセンサーを用いてデータが収集された。特に、岩のブロックとタワーの2つのセクターに焦点を当て、温度、ひずみ、振動、微小地震などのデータを分析した。データ処理にはオープンソースのJupyterノートブックが使用され、データの整合性や正確性が確保されるよう努めた。

結果
得られた結果から、特に温度要因が岩の表面に大きな影響を与えることが示された。データ分析により、温度変化が岩の変形とそれに続く安定性の崩壊に寄与するメカニズムを明らかにした。

温度変化による岩の変形メカニズム
熱膨張と収縮:岩石は温度変化に伴い、熱膨張や収縮を経験する。日中の高温や夜間の低温における温度差が、岩の表面や内部で異なる膨張・収縮を引き起こす。このプロセスにより、岩の内部応力が増大し、微小ひびや変形が進行する可能性がある。

熱疲労: 繰り返される温度変化は岩の疲労を引き起こし、最終的には構造的な弱化を招くことがある。長期間にわたる温度のサイクルは、岩の内部で累積的な応力を生じさせ、これがさらなるひび割れや断裂を誘発する。

水の影響: 温度変化に伴い、降水量などの湿度変化も影響を与える。降雨や雪が岩の亀裂に浸透すると、凍結と融解のサイクルによって亀裂が拡大し、岩の安定性が損なわれる。この過程は特に寒冷地域で顕著である。

微小地震との関連: 岩内の温度変動は、微小地震活動の増加とも関連している。温度変化による物理的な変形が微小地震の発生を誘発し、これがさらに岩の内部構造に影響を及ぼす。これにより、岩破壊の前兆としての微小地震を特定することが可能になる。

安定性の崩壊:温度変化による上述のメカニズムにより、岩は徐々に不安定化し、以下の過程で安定性の崩壊が進行します:

応力の蓄積: 長期的に蓄積された応力が臨界点を超えると、岩の変形が急激に進行し、さらなる亀裂が生じる。

破壊の進展: 初期の微小ひび割れが成長し、連結してより大きな裂け目や崩壊を引き起こす可能性がある。特に、温度変化や降水量のボトルネックが形成されるポイントで、急激な破壊が発生しやすくなる。

外的刺激への脆弱性: 温度変化が進行することで、物理的な強度が低下し、強風や激しい降雨などの外的刺激が加わると、岩の崩壊が誘発されるリスクが高まる。

季節的または日々の温度サイクルが岩の機械的性質に及ぼす影響を理解することは、今後の防災対策において重要である。また、微少地震や音響信号のモニタリングによって、岩の変形の前兆を特定することが可能になることを示唆した。

新規性
岩の不安定性を引き起こす要因としての温度、降雨、風などの環境条件に基づき、岩の変形に対する詳細なモニタリングとデータ分析を行った点にある。特に、複数の異なるセクターにおけるデータを比較分析し、異なるストレッサーが岩に与える影響を包括的に評価したことが、科学的および実用的な価値を持つ。また、最新の技術を用いたデータ収集と解析方法論は、今後の研究にも応用可能な新しいアプローチを提供している。

内容云々よりも、こういった計測を実施していることに、見習うべき点があるのだと思われます。

2025年8月16日土曜日

地すべり調査(ボーリング、微動、モニタリング)

Investigation of intermittent motion mechanisms in large landslides based on in-situ monitoring and microtremor survey - ScienceDirect

AI要約

背景

  • 中国南東部の丘陵地帯では、台風や豪雨による地すべりが頻発しており、甚大な被害をもたらしている。
  • 地すべりのリスクを軽減するためには、地すべり体の地質構造や地下水流動特性を把握することが重要である。
  • 従来のボーリング調査は、空間的な代表性に限界があり、地すべり内部の構造特性の空間的な変動に関する情報が不足している。
  • Yaoshan地すべりでは、2016年の台風により家屋の倒壊や道路の亀裂などの被害が発生し、人命や財産が危機に晒されている。

手法

  • 地すべりの変形と水文特性の関係を調査するために、水平変位、沈下、深層変位の定期的なモニタリングを実施した。
  • 2019年8月には、雨量計、間隙水圧センサー、アレイ式傾斜計、水位計を含むモニタリングシステムを設置し、リアルタイムでの長期モニタリングを行った。
  • 地すべり体の主要な断面における地下構造特性を特定するために、3つの二次元微動探査測線(D1-D3)を主要な滑動帯(Zone II)に設定した。
  • 二次元微動アレイ技術を用いて、二次元せん断波速度プロファイルを作成し、地すべり体の内部構造を詳細に再構築した。
  • 既存のボーリングデータと統合して、滑り面の深度と空間分布を包括的に分析し、岩相層序の区分を行った。

結果

  • 長期的な現場モニタリングにより、Yaoshan地すべりの全体的な三次元変形が可視化され、地すべりを牽引帯、主要滑動帯、抵抗帯の3つのゾーンに区分することができた。
  • せん断波速度プロファイルから、地すべり体の内部構造が明らかになり、低せん断波速度異常帯が潜在的な優先水文経路として特定された。
  • 岩盤層内のボウル状の部分(すべり面下の低速度帯)が、Yaoshan地すべりの動きを促進する上で重要な役割を果たしている可能性が示唆された。
  • 間隙水圧の変動と地すべり移動の相関関係が明らかになり、降雨時に間隙水圧が上昇すると地すべり移動が加速することが示された。

考察

  • Yaoshan地すべりの間欠的な動きは、滑り面の形態、複雑な地下水流動システム、および地質構造的特徴の相互作用によって制御されている。
  • 低透水性の岩盤層内に形成されたボウル状の地形は、浸透した雨水をトラップし、急速に蓄積させることで間隙水圧を上昇させ、地すべりの動きを誘発する可能性がある。
  • 低せん断波速度異常帯は、潜在的な優先水文経路として機能し、地下水の流れを集中させることで地すべりの安定性に影響を与える可能性がある。
  • 季節的な地下水涵養の変化が、滑り面沿いの間隙水圧の周期的な変動を引き起こし、加速と減速を特徴とする間欠的な動きにつながる。

新規性

  • 二次元微動探査法を適用することで、大規模地すべり体の内部構造を詳細に可視化し、従来のボーリング調査では捉えきれない潜在的な優先水文経路やボウル状の地形を特定した点。
  • 現場モニタリングデータと微動探査の結果を組み合わせることで、地すべりの間欠的な動きを制御する要因(滑り面の形態、地下水流動、地質構造)の相互作用を明らかにした点。
  • Yaoshan地すべりの概念的な水文地質モデルを構築し、地すべりのリスク評価や災害予測に貢献する情報を提供した点。

図10を見ると、すべり面下にもダメージゾーン(潜在的なすべり面)が拡大するというようなことでしょうか。そこの水圧上昇が地すべり挙動に関係しているといった、新たな視点での内容です。
アレイ式傾斜計とは何か?と思いましたが、どうやら多段式の孔内傾斜計のようです。昔から国内で使用されているものと変位算出方法が一緒なのかどうかはわかりませんが、検索すると以下のような”array”商品が紹介されていました。


優先順位が低く、長い間読まずに貯めていた文献ですが、もう少しで整理できそうです。AI様様です。大半の文献は読んでも記憶から消えていきますので、何かしらの手法でデジタル化し記録しておくべきなのでしょう。

一方で、入社から数年の若い子らは論文や報告書を書くのに慣れていません。速く書くためには慣れが必要ですし、他人の論文や報告書の多読が有効です。AIの使用に慣れてしまっている世代ですが、あえて論文を読む習慣を身に着けることも必要だと思います。何事もバランスですね。

PS-InSAR 炭鉱地域での適用例

 Land subsidence mapping and monitoring using modified persistent scatterer interferometric synthetic aperture radar in Jharia Coalfield, India | Journal of Earth System Science

AI要約

背景
ジャリア炭鉱地域は、数十年間にわたり地盤沈下の影響を受けている。この研究は、CバンドSARデータを用いて、地下炭鉱火災や地下採掘活動による地表のゆっくりとした変形を調査することを目的としている。

手法
ENVISAT ASARから取得した19枚のSAR画像を使用し、修正PS-InSARを用いて地盤沈下の多時系列分析を行った。
Perissin and Wang (2012)によって提案された修正PS-InSARでは、従来の方法が利用する永久散乱体(PS)だけでなく、部分的に相関する散乱体も永久散乱体として扱うことができるため、データの取得可能性が向上する。これにより、高いデコヒーレンスが発生する地域、特に植生がある農村地域や、散乱体の密度が低いエリアでも効果的にデータを取得し、地表のゆっくりとした変形を検出する能力が向上する。
デコヒーレンスとは、SARデータを使用する際に、散乱体が時間と共にその状態を変化させることによって生じる信号の相関性の喪失を指す。主な要因は以下の通り。

  • 時間的デコヒーレンス: データ取得時の環境条件(例、降水、風など)が変化することで、同じ位置からの信号が異なる状態を示すことにより、干渉パターンが乱れる。
  • 幾何学的デコヒーレンス: 撮影位置や角度の変化により、同じ対象物から得られる信号の位相が異なることで、干渉信号のコヒーレンスが失われる。

成果
主要な研究地点5か所での年間沈下速度は最大29 mm/年、累積沈下量は90 mmと記録された。また、PS-InSARの結果は、実地の沈下場所と一致した。

結論
PS-InSAR技術は、ジャリア炭鉱地域における地表のゆっくりとした変形の検出、モニタリング、マッピングにおいて実用的かつ効果的であることが確認された。この手法は、植生や農村地域における地盤沈下の監視に有用である。地盤沈下の正確な推定には、より多くの画像を分析することで改善される可能性がある。

2020年の文献です。当時の"修正"は今では反映されて表示されていないのかもしれません。プロに会ったら聞いてみましょう。

表面波 スタッキング方法の比較

 Optimal stacking of noise cross-correlation functions | Geophysical Journal International | Oxford Academic

AI要約

背景・目的
地震波形データのノイズクロスコレレーション関数(Noise Cross-correlation Functions; NCF)は、地殻の速度構造や監視、地盤運動解析などに広く使われているが、データ品質のばらつきやノイズの影響により、単純な平均(Linear stacking)では信号品質の劣化が課題となる。本研究は、8種のNCFスタッキング手法の性能を評価し、最適な方法を用途別にガイドする。

評価したスタッキング手法(8種)
Linear stacking(算術平均)
Robust stacking(ロバスト加重平均)
Selective stacking(選抜平均)
Cluster stacking(クラスタ平均、独自新手法)
Phase-Weighted stacking(PWS; 位相加重平均)
Time–frequency Phase-Weighted stacking(tf-PWS)
Nth-root stacking(n乗根的非線形平均)
Adaptive Covariance Filter stacking(ACF)

Linear stacking(算術平均)
すべてのNCFを単純な加算・算術平均として重み1で合成する、最も基本的な手法。すべてのデータが等しく扱われるため、外れ値やノイズの影響を受けやすいが計算は高速。

Robust stacking(ロバスト加重平均)
Pavlis & Vernon (2010)による反復的加重平均法で、リファレンススタックと各NCFの類似度・残差により重み付け。類似度が低いNCFは重みが低くなり、外れ値の影響を抑制する。重みはL1ノルム(内積)とL2ノルム(残差)の両方を用いて定義し、収束まで繰り返す

Selective stacking(選抜平均)
Pearson相関係数などにより「良質なNCF」のみを選抜して平均する手法。設定したしきい値(例:相関係数0.7や0.8)以上のみ重み1、それ未満は重み0として排除することで、質の悪いデータの影響を除去する。基準波形は線形スタックなどを使用可能

Cluster stacking(クラスタ平均)
機械学習のクラスタリング(k-means法など)でNCFを類似したグループに分類し、クラスタ中心(代表波形)同士を振幅比などで加重平均する手法。2クラスタの場合、代表中心の振幅比で重みをつけ、基準値を超える場合は両方加重、それ以外は高振幅側のみ採用する。NCFのパターンや特徴に応じて重みづけ可能な新手法

Phase-Weighted stacking(PWS; 位相加重平均)
Schimmel & Paulssen (1997)の手法で、各サンプルごとに瞬時位相のコヒーレンス(cosine/sineの平均)で重み付けする。コヒーレンスが高いサンプルは強調される非線形スタック。「harshness」パラメータ(本研究は2)で重みの急峻性を調整。外れ値や波形ノイズに強い

Time–frequency Phase-Weighted stacking(tf-PWS)
PWSの時間・周波数拡張版。Stockwell変換(S-Transform)やDOST(Discrete Orthogonal S-Transform)によりNCFを時系列かつ周波数ごとに分解し、両方で位相コヒーレンスによる重み付けを行う。より短周期・非定常なノイズ成分に対応。パラメータ設定により計算時間も異なる

Nth-root stacking(n乗根的非線形平均)
各NCFの絶対値のn乗根をとって合成し、最後に加算結果のn乗に戻す方法。振幅の高い外れ値の影響を抑えつつ信号のコヒーレンスを強調できる。nは根の次数(本研究では通常2乗根)。配列地震学等で利用される非線形フィルタ。

Adaptive Covariance Filter stacking(ACF)
Nakata et al. (2015)による手法。各NCFに対してスペクトル成分の共分散フィルター(適応型偏極化フィルター)を適用し、非コヒーレントなノイズ成分を抑制。その後、線形平均する。ノイズと信号の分離に特化し、地震探査や環境ノイズデータに有効。計算負荷は比較的大きめ

データ・解析対象
米国カスカディア地震観測網(陸・海底両方、複数ネットワーク)
騒音波形をRaw(未処理)とOne-bit(1ビット正規化)で比較
各手法について6つの評価指標:SNR, 位相分散, 収束速度, 時間変化, 振幅減衰, 計算時間

主な結果・知見
原則としてどの方法でも主要な到達波(ballistic/first arrival)は抽出できるが、位相・振幅の保存特性は大きく異なる。
Raw NCFs:Robust stackingが最適(全用途で推奨)
One-bit NCFs:PWSとNth-root stackingを除く全てが実用的
Tomography(速度構造):Linear, Robust, Selective, Cluster, tf-PWS, ACFが有効
Monitoring(地殻変動監視):Linear, Robust, Selective推奨
振幅関連(減衰・地盤運動):Linear, Robust, Selective, Clusterが良い
PWS法はSNRは高いが位相・振幅保存では不安定
Nth-root法はSNRやノイズ耐性は良いが位相情報が劣化する場合あり
Selective/Cluster stacking、ACFは用途によって有効
Robust stackingは全体的にバランスが良い
tf-PWSは計算コストが大きい
距離減衰や一部の目的ではLinear, Robust, Selective, Clusterが安定

実装・応用
PythonパッケージStackMasterをオープンソース公開
時系列データの一般的なスタッキングにも利用可能

線形、あるいはRobust、Selectiveのような単純な方法が多くの用途に適用可といった結論が 面白い点です。NCFにピークが出てこないときは前処理の工夫と共にスタッキング手法の工夫が必要ということでしょうか。

2025年8月15日金曜日

DualSPHysics with Chrono

 nvestigation of subaerial landslide–tsunamis generated by different mass movement types using smoothed particle hydrodynamics - ScienceDirect

AI要約

本研究は、陸上起源の落石・転倒(トップリング)など異なる質量移動タイプ(MMT)が引き起こすサブエリアル地すべり津波(SLT)の差異を、SPH法(DualSPHysics v5.0.4)とProject Chronoの連成により数値的に比較・評価したものである。従来は滑動(sliding)に偏った研究が多く、落下(falling)・転倒(overturning)の波特性は未解明部分が多かった。本研究は計26ケース(長方体・円柱ブロック、岩石想定密度ρs=2700 kg/m3)で、最大波高・振幅、減衰、方向分布などを整理し、無次元パラメータ(Froude数F、相対厚S=s/h、相対長L=l/h、相対幅B=b/h、相対質量M)の関数として新たな経験式を構築した。

土砂(水陸質量)と水の連成方法(数値手法・境界・結合)
連成の基本枠組み:
オープンソースのDualSPHysics v5.0.4(WCSPH)を用いて流体(自由表面水)をSPH粒子で離散化し、反転型の剛体運動はマルチフィジクスエンジンProject Chronoで剛体力学を解く弱結合(loose coupling)。SPH側で算出した流体力をChronoに受け渡し、Chronoが剛体の平動・回転加速度を更新し、これをSPHにフィードバックする時刻歴連成。ChronoはDVI(微分変分不等式)に基づく非平滑マルチボディタイムステッピングで接触・ヒンジ拘束を解く。両ソルバは独立タイムステップで、1 SPHステップ中にChronoが複数サブステップを持ちうる。データ連携はDSPHChronoLibインターフェース。落下型はChrono拘束を用いず重力自由落下。

流体方程式・数値安定化:
ナビエ–ストークスのSPH離散。弱圧縮性(WCSPH)で状態方程式を用い、音速比で圧縮性を制御(c0=71.2 m/s、密度変動<1%を確保)。数値粘性は人工粘性(α=0.01)を採用。

境界条件:
固体・壁はmDBC(modified Dynamic Boundary Condition)を採用。固定境界を静止流体粒子として扱い、ゴースト粒子の線形外挿で密度場を平滑化、壁近傍の圧力安定性と隙間抑制を図る。圧力場安定のため密度拡散項も併用。

計算設定例:
粒子間隔dp=12 mm、CFL=0.2、音速係数17、GPU(RTX A5000)で4秒計算に約40時間(約8300万粒子)。検証では氷密度相当、主試験では岩石を想定して剛体密度ρs=2700 kg/m3を使用。

検証・収束:
大規模実験(氷ブロックの落下・反転)に対し一次波峰a1は±10%、一次波高H1は±15%以内で整合。空中スプラッシュの影響が小さい測点ほど良好。空間収束は反転型で難しさが残るがdp≤12 mmで一次波は概ね良好。時間収束もCFL 0.05–0.20でaMの差は最大約1.3%。

結論
転倒は落下よりも大きいaM, HMを生じ、ただし沿軸減衰は速い。
新しい経験式(発生・伝播)はF, S, B, L, r/h, 方向γの関数で±50%程度の精度帯で整合。
実事例(カピトリオ)での適用は、適用範囲逸脱や地形・運動複合性により誤差が増えるが、初期評価に有効。
滑動との関係式により、MMT間の波規模変換の目安を提示。今後は粒状体、緩斜面、複合運動、適用範囲拡大が課題。

2025年。DualSPHysics で Chrono との連成は以前よりサポートされていますが、このようなことができるのですね。

2025年8月14日木曜日

single-layer multi-phase formulation

 Modelling internal erosion using 2D smoothed particle hydrodynamics (SPH) - ScienceDirect

AI要約

1. 研究の背景と従来の課題
内部浸食現象のモデリングにおいて、従来の連続体モデルは微細粒子の進化や固相変形を十分に考慮できていなかった。また、有限差分法(FDM)や有限要素法(FEM)のようなメッシュベース手法は、大変形を伴う場合にメッシュの歪み問題に直面し、全破壊プロセスの解析が困難である。近年、メッシュフリー手法であるSPHが内部浸食の解析に適用され、ある程度の成功を収めているが、既存のモデルには詳細な説明の不足、検証の欠如、材料強度への浸食による影響の不考慮、数値安定化の課題などが残されていた。

2. 本研究の提案手法
本研究では、以下の点を改善した新しいモデルを提案。

  • 飽和/不飽和条件および浸透誘発浸食、輸送、堆積を記述するためのu-w-p定式化に基づく3相5成分SPHモデル。計算は single-layer multi-phase。
  • 微細粒子含有量と材料強度間の非線形関係を考慮するために修正された「Ma et al. (2022)」の微細粒子依存構成モデル。
  • 浸透-浸食-変形連成プロセスに適応するように拡張された境界処理および拡散アルゴリズム。数値シミュレーションをさらに安定化させるための新しい粘性散逸項。

3. モデルの基礎と構成要素
浸食性多孔質材料を、固相(不浸食性固体骨格と浸食性固体微細粒子)、液相(液体化浸食微細粒子と間隙水)、気相(空気)の5成分からなる多相多成分連続体として扱う。各相/成分の体積分率が定義され、これにより各相の物理量を定量化する 。
支配方程式は、質量保存則と運動量保存則に基づいて定式化され、質量バランス方程式は固形骨格、微細粒子含有量、液体化浸食微細粒子、水密度について記述される。全応力はBishopの有効応力概念を用いて計算する。
構成則は以下のとおり。

  • 浸食/堆積法則: 浸食と堆積のメカニズムを記述するため、「Cividini and Gioda (2004)」による浸食項を修正し、実験結果に合わせた経験的関係を利用。
  • 土の構成モデル: 弾塑性理論とDrucker-Prager降伏基準が適用され、材料強度パラメータ(粘着力と摩擦角)は、ひずみ軟化、吸引圧による硬化、および微細粒子含有量への依存性が考慮される。特に、微細粒子含有量と材料強度の非線形関係を組み込むために「Ma et al. (2022)」のモデルが修正されており、浸食による強度低下が表現されている。
  • 水の構成モデル: 間隙水圧は弱圧縮性流体の状態方程式から計算され、飽和度と吸引圧の関係にはVan Genuchtenモデルが用いられる。透水係数は飽和度とKozeny-Carman式に基づく空隙率の変化に応じて変化する。

2024年の発表です。実務で利用することはないですが、2Dとは言え他国はここまで進んでいます。キーワードに DualSPHysics が載っていますので、いずれ実装されるのでしょう。私は2相2成分で止まっていますが、何とかしないといけないと思わされます。


SSRFEM + SPH

Post-Collapse Evolution of a Rapid Landslide from Sequential Analysis with FE and SPH-Based Models

AI要約

研究対象・地形・地質
Sant’Andrea地すべりは北東イタリア(アルプスのPerarolo di Cadore)に位置し、ボイテ川の左岸に隣接している。崩壊時に河川を堰き止めて下流集落に影響を及ぼす可能性が高い区域。
旧地すべり区域を含む72,000m²の範囲。主な地層は、厚さ30mのデブリ(A/B/C層)と、その下部に石膏・硬石膏を主とする岩盤(D層)。土壌・地質調査や地表変位のモニタリング(RTS等)が継続的に行われている。

FEM(有限要素法)解析
3D-FEM(Midas GTSNX)で現地地形・地質情報をもとにモデル化。地層・水理条件(降雨時の水位など)を詳細に考慮。
FEM解析では、剪断強度低減法(SRM: Strength Reduction Method)を使い、滑落面と不安定体積を評価。
乾燥状態・湿潤(雨後)状態でそれぞれ解析を実施し、「最悪ケース」の水理条件で土壌強度を低減した。

不安定体積(Unstable Volume)の詳細
地表変位分布から「iso-displacement surface(等変位面)」を用いて安定領域と不安定領域を識別。「地表変位の最大値の10%」を閾値として不安定体積を定義。
滑落面・分布は地質・水文条件で大きく変化。乾燥時は深く広範囲、湿潤時は浅く面積が狭い傾向。乾燥時の方が多いが、実際には湿潤時の方が危険度が高い(安全率が下がり引き金となるため)。

10%の根拠
・変位分布と実測データの適合性:モニタリングデータとモデル計算結果の分布が一致するように調整。

  • 7.5%〜15%の範囲を候補とし、その中で安全側かつ現実的な値として10%を採用。
  • 最大せん断ひずみゾーンと閾値により抽出された領域との整合性を確認し、閾値設定の妥当性を支持。
  • 10%はより保守的(安全側)な設定であり、過小評価を防ぐ目的で選定。

SPH(Smoothed Particle Hydrodynamics)による伝播・堆積予測
上記FEMで得た不安定体積をSPHモデルに入力し、土砂流の伝播距離・堆積分布を解析。浅水方程式+Völlmyモデル(粘性・摩擦抵抗の両効果)を採用。特に湿潤時には流動性が高まり、土砂がボイテ川谷でダム状堆積(厚さ15-20m)となった。摩擦角の値(μ=0.45, 0.58)などの設定も安全側・リスク評価に寄与。

感度解析の記載
「摩擦係数(μ)は0.36~0.62、摩擦角ϕとしては20度~32度の範囲で変化させて計算し、土砂の流下距離・堆積の広がりにどの程度影響するかを検証」。具体的には、摩擦角(ϕ)を20〜32度で7パターンに設定し、ξ=500ms^-2で各パターンをシミュレーション。
結果として、「runout distance(流下距離)は摩擦角の値に敏感に反応したが、乱流係数(ξ)の影響は小さかった。

解析結果の詳細
乾燥時は不安定体積が多いが、流動性が低いため堆積量・範囲は限定的。湿潤時は不安定体積はやや少ないが、流動距離(runout)が長く、河川全体を埋める可能性がある。

結論・意義
複数の水理条件下で不安定体積と安全率(FS)を評価。
湿潤時(大雨後)が最も危険度が高く、比較的小さい体積でも流動性増加により下流域のリスクが顕在化。SPHモデルの摩擦係数(μまたは摩擦角ϕ)は、土砂流動・停止位置・堆積形状に対して極めて高い感度を示し、「runout distance」に大きく影響する。

SSRMを利用して不安定領域を求め、SPHで流す2段階の手法です。2021年です。私はSPHのみで解けるように組みましたが、発表するとしても2027年。「摩擦係数(μまたは摩擦角ϕ)は、土砂流動・停止位置・堆積形状に対して極めて高い感度を示す」という結果も同じ。のんびり組んでいたら時代遅れになっちゃいました。 

Windows 11 clamshell mode

Windows 11のラップトップでクラムシェルモードを設定していたのですが、時々切断されます。再度、設定を確認。

1. 電源設定
「コントロールパネル」「ハードウェアとサウンド」「電源オプション」を選択。
左メニューの「電源ボタンの動作の変更」をクリック。
「カバーを閉じたときの動作」で、「電源接続時」を「何もしない」に設定」←これが未設定でした。

2. ディスプレイ設定の調整
デスクトップ上で右クリックして「ディスプレイ設定」を開く。
複数ディスプレイが認識されていることを確認する。
「2のみ表示」「このディスプレイの接続を切断する」で内蔵画面はオフになる。

CAD使用時に頻繁に現れるので GPU ドライバも確認しましたが、最新版。BIOSも最新版。ドック等のファームを最新版にUPし、ひとまず様子見です。

2025年8月13日水曜日

AMD + nvfortran

Dockerが走らなくなっていたPCは再インストールから。

$sudo apt update
$sudo apt install -y nvidia-container-toolkit
$sudo systemctl restart docker

$sudo docker run --gpus all --ipc=host --ulimit memlock=-1 -it --rm -v /media/user/Data/:/workspace/data/ nvcr.io/nvidia/nvhpc:25.5-devel-cuda_multi-ubuntu24.04

$cd workspace/data/src_gpu
$make

makeファイルのオプション:
FLAGS = -O3 \
         -byteswapio \
         -tp=znver4 \
         -Mfma \
         -Mcache_align \
         -Mvect=simd \
         -acc=gpu,multicore \
         -gpu=cc89,mem:managed \
         -mp=allcores,bind \
         -Minfo=accel,mp,inline

-O3: 最適化レベル
-byteswapio: I/O操作でバイト順序を入れ替え
tp=znver4: AMD Zen4アーキテクチャを対象プロセッサとして指定
-Mfma: 融合乗算加算(FMA)命令を有効化
-Mcache_align: キャッシュ利用効率向上のためのデータアライメント
-Mvect=simd: SIMDベクトル化を有効化
-acc=gpu,multicore: GPUとマルチコア加速のためのOpenACCを有効化
-gpu=cc89: Capability 8.9(Hopperアーキテクチャ)を対象
mem:managed: CUDA管理メモリを使用
-mp: OpenMPの並列処理設定
allcores: 利用可能な全CPUコアを使用
bind: スレッドをCPUコアに固定
-Minfo: コンパイラ情報出力
accel: アクセラレータ(GPU)コード生成情報
mp: マルチプロセッシング/OpenMP最適化情報
inline: 関数インライン化情報


Memtest86+

GRUB 操作時に不意に再起動する現象に遭遇。

PC購入時より、計算中にメモリを使いきると再起動する問題があったので、メモリテストを実施。するとメモリテスト中にも再起動が発生。

メモリを抜き差しして、再度テスト。3時間くらい走らせて止めました。
計算で負荷をかけて様子を見ます。


Disable the Nouveau Driver on Ubuntu 24.04

引き続き、お盆の間にPCのメンテナンスです。

後輩君が使用していた計算用の Ubuntu 24.04 v6.14 がカーネルパニックを起こして起動しなくなりました。引っかかていた NVIDIA ドライバを UPDATE したら解決するだろうと思っていたのですが、挙げてみるとたカーネル v6.11 も起動しなくなりました。どうも、Nouveau(オープンソースの NVIDIA ドライバ)と衝突している模様。カーネルの起動オプションに nomodeset を追加して Nouveau ドライバのモードセット(modeset)機能を無効化。

設定方法は以下の通りです。

一時的にnomodesetを使う方法(起動時のみ):
PC起動時にGRUBメニューが表示されたら、Ubuntuの起動エントリを選択し、eキーを押して編集モードに入る。
linux行の中にある quiet splash の部分を探し、その後ろに nomodeset を追加する。
F10キーを押して起動する。

恒久的にnomodesetを設定する方法:
/etc/default/grub
変更前:GRUB_CMDLINE_LINUX_DEFAULT="quiet splash"
変更後:GRUB_CMDLINE_LINUX_DEFAULT="quiet splash nomodeset"
v6.14はまだ起動しませんので、GRUB内の起動順序を入れ替えます。
変更前:GRUB_DEFAULT=0
変更後:GRUB_DEFAULT="1>2"
$ sudo update-grub

これで起動できました。うーん、対処療法。


2025年8月12日火曜日

リモートデスクトップ(GNOME)のパスワード固定 Ubuntu 24.04

22.04でも設定し、その後忘れてしまったリモート時のパス固定。

GNOME Remote Desktop はパスワードをキーリングに保存しますが、起動時にキーリングがロックされているとパスワードが読み込めず、新しいパスワードを自動生成してしまいます。ダミーのキーリングを作成し、パスワードを暗号化しない(空パスワード)状態にすることで、キーリングのロック問題を回避し、パスワードが固定されるようにします。

パスワードを固定する手順

1. Seahorse(パスワードと鍵)をインストール・起動
$ sudo apt install seahorse
$ seahorse

2. 新しいパスワードキーリングを作成
Seahorse の左ペインで「新しいパスワードキーリング」を追加
任意の名前を付けて作成
パスワード入力画面が出たら、リターンキーを2回押して空パスワードで保存
作成したキーリングをデフォルトに設定

3. 作成したキーリングを右クリックし、「デフォルトにする」を選択
リモートデスクトップのパスワードを再設定

4. Seahorse で元の「ログイン」キーリングに戻す
Seahorse を再度開き、作成したダミーキーリングに「GNOME Remote Desktop」関連のパスワードが保存されていることを確認
元の「ログイン」キーリングを右クリックして「デフォルトにする」に戻す

5. その他
オートログイン設定
「設定」-「システム」-「ユーザー」を選択
「自動ログイン」をオン

自動ロック解除(画面ロック無効化)方法
「設定」-「プライバシーとセキュリティー」-「画面ロック」
スイッチをオフ


S波速度構造の推定と機械学習

 Surface wave dispersion curve inversion using mixture density networks | Geophysical Journal International | Oxford Academic

AI要約

目的と背景
浅部S波速度構造の迅速・高精度な推定は耐震・地盤工学などで重要。従来の分散曲線逆解析(線形・大域探索・ベイズ)は、局所解・計算量・不確実性評価・層数/厚さ固定などの課題がある。

提案手法(2ステップML)
ステップ1:分類NNで上部100 mの層数(2–7層+半無限層)を推定。
ステップ2:Mixture Density Network(MDN)で層境界深さ(決定論的出力)と各層S波速度の確率分布(混合ガウス10個)を同時推定。Love・Rayleigh基本モードの分散曲線(1–20 Hz)と各周波数点の不確かさを入力。

学習データと前処理
2–7層モデルで上部100 mをパラメタ化。最上層Vsは100–500 m/sで最小厚さ5 m、以深はVs∈[Vtop+100,1600] m/sで最小厚さ10 m、半無限層はVs∈[max(v)+200,2000] m/s、境界深さ最大85 m。P波速度と密度は経験式で付与(Vp=1.16Vs+1.36、ρ=1.74Vp^0.25)。
gplivemodelでLove/Rayleigh分散曲線を計算し、1–20 Hzで100点に対数再サンプル。0–5%の一様ノイズを加え、同時に不確かさベクトルuを作成し入力に含める。標準化後、学習/検証=90/10に分割。
各層数ごとにMDNを個別学習(出力次元は(層数+2)×10)。平均推定は事後分布の最大値を採用。ハイパーパラメータはHyperopt系で探索。

性能評価(合成データ)
層数分類:未学習4000モデルで正答>60%、±1の誤差が約30%、相関係数r=0.90。中間層数(4–6層)の識別がやや難しい。
MDNのVs・深さ推定:2–4層で高相関(各層r>0.84、深さr>0.80)。5–7層では深部ほど相関低下(Vs 0.62–0.87、深さ0.29–0.67)。半無限層Vsは全体で高相関(r>0.84)。
代表例:3層では狭い事後分布で真値に近い。7層では混合核の多峰性が現れ、曖昧性を示すが大筋で整合。

従来法との比較(合成データ)
比較対象:GEOPSY DINVER(近傍算法NA)とベイズMcMC(Alder et al. 2021 実装)。パラメータ空間は最大7層、Vs=100–2000 m/sに整合。
2層モデル:分類は正解、MDNは深さ・Vsとも真値に近く、推定は1秒程度。NAも良好(約30秒/1万モデル)。McMCも良好だが約2 CPU時間/20万モデル。
5層モデル:MDNは平均値が真値に近く、深さ誤差<5 m、推定は1秒程度。NAは目標適合は良いが真の構造から逸脱(層境界の位置ずれ・深部不解像)、約60秒/3万モデル。McMCは大筋整合、約7 CPU時間/50万モデル。

実データ適用(ミュンヘン)
市内L字配列(辺100 m、4.5 Hzセンサー11台、2時間観測)。HRFKで1–20 HzのLove/Rayleigh分散曲線を抽出。
分類は5層、MDNの結果は5 m付近で速度段差、その後は圧密によると考えられる速度勾配。水位・堆積物境界(GeoPot水文地質3Dモデルの砂→粘土)と整合し、20 m付近の層境界とも概ね対応。

ロバスト性検証
学習外の滑らかな速度勾配(200→1800 m/s)に対して、分類は最大の7層、MDNは構造を再現。順解析の分散曲線は入力に近く、非一意性を示唆。
低周波欠測(1–3 Hzのギャップ)を4種の外挿で補完し、不確かさを大きく設定。分類結果は3–4層に分かれたが、MDNの事後分布を平均すると浅部構造と浅い層境界は安定。30–60 mで差異が出るが最大確率は真値に近い。
深部解像度の限界例として、深部の薄い低速度層(厚さ10–30 m)の感度テストでは、10 mではノイズ下で識別困難、20 mで差がやや現れ、30 mで明瞭。深部ほど最小可探厚が増すことを示唆。

議論・利点と限界
利点:層数を固定せずに分類し、Vsの事後分布と層深さを同時に高速推定。学習後の推論は秒オーダーで、多点・多次元適用に有利。複雑構造に対しNAより高精度、McMCに近い結果を高速に得られる。
限界:深部の解像度低下と非一意性により、層数中間域の分類や深い層境界の推定は不確かさが増大。MDNは事前サンプル型で、観測近傍のサンプル密度が低い場合、McMCより事後分布が広くなる可能性がある。ただし最尤近傍は真値に近い。
実務上の注意:1–20 Hzの広帯域データが望ましい。欠測は外挿+大きな不確かさ付与で平均化戦略が有効。7層以上の極薄層は本帯域では解像困難。
将来展望:H/V比の併用、上位モード導入、説明可能AIの活用が精度・解釈性向上に有望。

 層数を指定することなく、層境界深さ・Vsの確率分布を同時に推定できる、しかも機械学習なので高速というのは従来法に比べて有利で理想的です。データを作って教師データとするのもうまいと思います。学習時の不確かさの作り方は書かれていますが、実データでどう定量化したかは書かれていません。コードを追わないとわからないのかな。


2025年8月11日月曜日

分散曲線の自動ピッキングとドメインシフトへの対応

 Domain Adaptation in Automatic Picking of Phase Velocity Dispersions Based on Deep Learning - Song - 2022 - Journal of Geophysical Research: Solid Earth - Wiley Online Library

AI要約

本研究は、深層学習を用いた位相速度分散曲線の自動ピックアップにおける「サンプル不適合性(sample-inadaptability)」の問題、特に訓練データと異なる地質領域のテストデータでニューラルネットワークの性能が低下する課題に取り組んだ。これは、従来の表面波トモグラフィーが基礎モードのみに依存し、非一意性と精度が低いという問題や、手動ピックアップの時間とエネルギーの消費をディープラーニングで解決しようとする中で生じたものである。

背景1:分散曲線ピッキングにおける時間と労力の課題
地質構造を正確にマッピングするためには、基本モードと上位モードの分散曲線をピッキングするのに時間を要する。従来、位相速度-周波数図からの分散曲線の手作業によるピッキングは、かなりの時間を必要とた。そのため、深層学習を用いた分散曲線の自動ピッキングが試みられていた。

背景2:深層学習モデルの一般化能力の低下(サンプルへの不適応性/ドメインシフト)
深層学習を用いた自動ピッキングにおいて、ある位相速度-周波数図(トレーニングセット)で訓練されたニューラルネットワークが、トレーニングセットとは異なる地質学的領域のテスト図に対しては、その一般化能力や性能が低下するという問題が発生していた。これは「サンプルへの不適応性」と呼ばれ、コンピュータビジョンにおける「ドメインシフト」の問題と同義である。具体的には、東北中国のデータで訓練されたネットワークが、他の地質学的領域からの位相速度-周波数図(例えばアメリカのデータ)に対して性能が低下するという現象が観察された。

この問題に対処するため、本研究では深層学習ベースのピックアップ手法に「ドメイン適応(domain adaptation)」ステップを導入した 。具体的には、画像処理で一般的に使用される「ガンマ変換(gamma transform)」を用いて、位相速度-周波数図の画像のコントラスト(グレーレベル分布)を変更した 。ガンマ変換により、テストデータのグレーレベル分布を訓練データに近づけることで、ニューラルネットワークの汎化性能を向上させた。

分散領域の抽出には Res-Unet++ を使用し、位相速度-周波数図はF-J変換によって取得する。F-J変換は、高次モードの分散曲線を効果的に抽出できる特徴を持っている。最終的な分散曲線の特定には「追求法(pursuit method)」が用いられる。

主な成果と利点:
ガンマ変換により画像のコントラストを増加させると(ガンマ値γを減少させると)、ニューラルネットワークはより多くの分散点を特定した。例えば、合成データにおいて、γ=1からγ=0.7に変化させることで、総計92点の分散点が追加で選択され、その平均精度は99.4%以上であった。
画像の平均グレースケール分布(MGLD)がトレーニングセットの平均MGLDに近いかそれよりも高い場合、高次モードの分散領域がより効果的に抽出された。
この手法は、転移学習や新しいニューラルネットワークの再訓練なしに、訓練されたネットワークの汎化性能を向上させるという利点がある。

課題と注意点:
画像のコントラストを過度に増加させると(γを小さくしすぎると)、ニューラルネットワークがノイズなどのアーティファクトも抽出してしまう可能性があるため、適切なガンマパラメータの選択が重要である。


分散曲線のピッキングはまだ手動です。機械学習で自動化するステップ、そこで問題となるドメインシフトに対応するステップの2手遅れています。
社内の地質部門はさらに遅れています。直営で物理探査を実施していた世代のおじさん達が効率化のもと外注化を進め、技術を伝達せずに退職し、我々の世代もできる人が地質から離れてしまいました。今の若い人たちは将来どうなるのでしょう?私のチームの後輩君の一人には理論もコードも伝えているので、自力で解析できるようになっていますが、地質を知らない子ですからね。
ま、いずれにしてもドメインシフトの問題は機械学習を利用する上で避けて通れない問題です。このような取り組みがあること自体、参考になります。

分散曲線の自動ピッキング その2

Deep learning‐based extraction of surface wave dispersion curves from seismic shot gathers - Chamorro - 2024 - Near Surface Geophysics - Wiley Online Library

AI要約

この研究では、地震ショット収集データから表面波の分散曲線を深層学習を用いて直接抽出する新しい方法を提案した。

背景と課題:
多チャンネル表面波分析(MASW)では、表面波の伝播特性を利用して地盤の情報を得るが、分散曲線の抽出は通常手動で行われ、品質管理が必要であり、多大な労力を要する。

提案手法::
深層学習(residual network)を使用して、地震ショット収集データとそれに関連する分散曲線との直接的なマッピングを確立する。これにより、分散スペクトルの計算が不要となる。

データ生成:
1Dの圧縮波及びせん断波速度、密度モデルを使用して、1000組のショット収集データとその分散曲線を生成した。

結果:
合成データセットと実際のフィールドデータの両方において高品質な分散曲線の予測を実現した。結果は地域の地質と一致し、信頼性の高い速度モデルを提供できた。

不確実性の定量化::
モンテカルロドロップアウトと分布パラメータ推定技術を用いて予測の不確実性を評価し、予測プロセスの品質を確認した。

モデルの可視化:
説明可能なAI(XAI)技術を用いて、ネットワークの学習プロセスを理解し、どのデータが予測に寄与しているのかを可視化した。可視化された結果は、モデルが特定の周波数に対してどのように反応しているのかを示し、解析の透明性を高めることに寄与した。

結論:
分散曲線の自動抽出プロセスを効率化し、フィールドデータに対する分散曲線の予測精度を向上させ、実用的な逆解析の結果を得ることができることを示した。


その1の文献よりは受け入れやすい内容です。データ解析の自動化、省力化に寄与することが期待されます。

分散曲線の自動ピッキング その1

Automatically Extracting Surface‐Wave Group and Phase Velocity Dispersion Curves from Dispersion Spectrograms Using a Convolutional Neural Network | Seismological Research Letters | GeoScienceWorld 

AI要約

研究の背景と目的:
アンビエントノイズトモグラフィー(ANT)は地殻や上部マントルの構造を画像化するために広く用いられているが、その重要なステップの一つが、連続的なアンビエントノイズ記録の相互相関関数から地表波の群速度および位相速度分散曲線を抽出することである。
従来の手法では、手動による抽出は非常に労働集約的で時間がかかる。自動抽出手法も、特に位相速度データにおいては信頼性が低い場合がある。本研究の目的は、これらの課題を解決するため、畳み込みニューラルネットワーク(CNN)を用いた新しい手法「DisperPicker」を提案し、群速度と位相速度の両方の分散曲線を自動的かつ正確に抽出することである。

提案手法「DisperPicker」:
DisperPickerは、CNNと特定の抽出戦略、データ品質管理基準を組み合わせた手法である。CNNの入力は周期-速度ドメインにおけるペアになった群速度および位相速度分散スペクトログラムである。これにより、群速度分散曲線が位相速度分散曲線の抽出を暗黙的にガイドし、位相速度抽出における曖昧さを軽減する。CNNのアーキテクチャは、U-net構造をベースに問題に適応するように簡略化されている。出力としての分散曲線は、ガウス関数を用いて確率画像に変換されている。これにより、CNNの訓練が容易になる。

データセットと訓練:
訓練データには、中国の3つの異なる密な地震アレイ(蘇州市、長寧-昭通、濰坊市)から収集された短周期地表波データが使用された。合計32,404組の群速度および位相速度分散曲線が手動で抽出され、訓練(28,404組)、検証(2000組)、テスト(2000組)に分割された。

結果:
テストセットでの平均損失は1.03であり、過学習がないことが示された。DisperPickerは、従来の自動抽出プログラムAutoPickerと比較して、より多くの分散点を抽出し、群速度と位相速度の両方でより高い精度を示した。特に位相速度の抽出精度が大幅に向上している。DisperPickerで抽出された分散データは、浅層のVSモデルインバージョンにおいても信頼性があることが確認された。

ただし、上記の説明では品質管理基準の具体的な内容や、どのように組み込まれるか(例えば、アルゴリズムや閾値など)については、記載がありません。

2025年8月10日日曜日

表面波探査による地下水位推定

 Physics‐Guided Deep Learning Model for Daily Groundwater Table Maps Estimation Using Passive Surface‐Wave Dispersion - Cunha Teixeira - 2025 - Water Resources Research - Wiley Online Library

AI要約

この研究は、限られた空間的・時間的観測データしか得られない地下水位(GWT)のモニタリングにおける課題に対処するため、DNN(MLP)と受動的な表面波解析(passive-MASW)を組み合わせたアプローチを提案している。
研究サイトでは、鉄道によって誘発される表面波を連続的に記録し、そこから毎日分散曲線(DCs)を測定している。 鉄道線路と平行に5本の測線が配置されており、42個のジオフォンが3m間隔で設置されている。

MLPモデルは、入力としての分散曲線と、出力としての対応する地下水位の深さを用いて訓練された。モデルはセンサーアレイ全体にわたるGWT深度を推定し、日々のGWTマップを生成する。モデルの性能は、訓練ピエゾメーター地点でR²が80%、RMSEが0.03mという高い精度を示した。さらに、訓練データに含まれていない別のテストピエゾメーター地点でも、R²が68%、RMSEが0.03mと良好な結果を示し、空間的な外挿能力があることが実証された。この研究は、ディープラーニングが地震データからGWTマップを推定する上で有効な手段であり、空間的に限られたピエゾメーター情報のみで広範囲の地下水動態を効率的に監視できる実用的かつ効率的な解決策を提供する可能性を示している。

表面波分散の物理::
レイリー波の位相速度(VR)は、地中の周波数(または波長)に応じたS波速度(VS)の構造に依存する。短波長は浅い深度のVSを反映し、長波長は深い深度のV_Sを反映する。この関係は非線形である。

物理モデルの考え方:
地下水位(GWT)の変動が地盤の弾性特性、特にS波速度(VS)に影響を与え、その結果として表面波(レイリー波)の位相速度(VR)の分散特性が変化するという、物理的な関係性を取り扱っている。、地下水位が上昇すると速度は低下し、地下水位が下降すると速度は上昇する。この関係はすべての波長で見られるが、特定の波長でより顕著になる。これは、地下水位の変化(水で飽和しているか否か)によって地盤の剛性が変化し、それがS波速度に影響を与え、最終的に表面波の位相速度に影響するという物理的なメカニズムに基づいている。
Physics-Guided」という言葉は、この地下水ダイナミクスと地震波速度構造の間の物理的な因果関係をモデルの基本的な入力と出力の関連性として捉え、データ駆動型のアプローチの基盤としていることを示している。つまり、地下水位の物理的な変動がどのように地震波の伝播特性に影響を与えるかという知識が、モデルの設計と解釈の指針となっている。

単純に、上載荷重の変化がVsに効いているのでしょう。
あと、Physics-informed と -guided は別物なのですね。もう間違えません。


点群からロックボルトの抽出

 A physics-guided hierarchical deep learning framework for underground rock reinforcement compliance check based on 4D point cloud data - ScienceDirect

AI要約

この論文では、地下鉱山での岩盤補強(ロックボルト)の設置パターンの適合性チェックを自動化・高精度化するために、階層型深層学習フレームワーク(PGHDFramework)の利用を提案している。従来は人手に頼ったため誤差や労力が大きかったが、SLAM-LiDARを使った4次元点群(空間座標+反射強度)データを活用し、手法の汎用性と精度を向上させた。提案手法は、まず物理ベース(LiDAR強度情報)のフィルタリングとクラスタリングでボルト候補を抽出し、その後、PointNet++をベースとした階層型ニューラルネットワーク(4DBDNet)で空間・物理属性を統合した点群セグメンテーションを行う。評価実験では他の先行手法と比べて高い精度・再現率・IoUを達成し、新設坑道でもボルト間隔自動算出による適合チェックを実現できることを示した。今後は粉塵・腐食・視点遮蔽など現場の制約への対応が題として挙げられる。

物理モデルと損失関数について:
本手法の物理的特徴づけの中心は、「強度(LiDAR Intensity)」である。LiDARの強度値は、ターゲット(ボルトと岩盤、金属メッシュ等)の表面反射率ρ、およびレーザーの入射角などの物理パラメータで決定される。論文ではWeitkampのLiDAR方程式を元に、ボルトと非ボルト領域の強度値差が主に表面反射率(材質の違い:金属vs岩盤)に由来する点を物理的に説明している。これに基づき、強度値によるしきい値処理(Intensity Thresholding)が分類に利用されている。
また、Deep Learning部分(PointNet++を用いた4DBDNet)では、強度値(i)を空間座標(x, y, z)とともに標準化し、点群特徴量として多階層で統合的に扱い、ロックボルト識別精度を向上させている。損失関数そのものには直接物理方程式を組み込んでいないが、物理的に導かれた強度しきい値と、それに基づく点群フィルタ/クラスタリング結果を教師データのラベリングに利用し、点群の空間的実測データに物理属性を付与した学習データを生成し、データ駆動型学習を物理的根拠で補強している。

点群の利用方法について詳細:
SLAM-LiDARで取得した点群は(x, y, z, intensity)の4次元データとして保存される。強度しきい値(Intensity Thresholding)で非ボルト点を物理モデルベースでフィルタリング。強度値の違いは主に材質違いなので、しきい値を可視的に調整した。
DBSCANクラスタリングで、密度ベースでボルト点群を抽出。さらにOctreeベースの接続成分解析(CCA)で個別ボルトクラスタを分離。ボルト形状の物理寸法に基づくしきい値処理(dimensional thresholding、buffering)でボルトサイズ・形状に合致するものを抽出。得られたボルトクラスタにラベル付与し、各点について(x, y, z, i, label)の形で教師データ化。値はZ-score標準化し、空間的バラツキと強度の違いを均等に扱えるようにした。
最終的にPointNet++をベースとした4DBDNet(4D Bolt Detection Neural Network)に入力し、多階層・多スケールで空間+物理情報を融合し、点ごとのボルト/非ボルト分類を学習・推論させた。
精度評価は、ボルトレベル(クラスタ単位)、点群レベル(点単位)で実施。精度・再現率・F1スコア・IoUなど多面的指標で先行研究と比較した。

まとめると、本論文の物理モデルは点群の「強度」と「空間座標」を融合し、その違い・相関をラベリング等に活用することで、従来法では解像不十分だった非突出型(Split Set)ボルトの検出と識別精度を大幅に高めている点が最大の特徴である。

損失関数に物理モデルを組み込んでいるのかと期待しましたが、そうではない様子。国内のトンネルではロックボルトは整然と打たれていますので、このようなモデルは必要ないでしょうね。
 

GNN + CAE

 [2107.11529] Graph neural networks for laminar flow prediction around random 2D shapes

書名だけ見て購入し、失敗した図書で引用されていた文献。図書については後日。

AI要約

2次元ランダム形状周りの層流予測問題において、GCNN用いたサロゲートモデルを提案。従来のピクセルベースCNNとは異なり、CFD(数値流体力学)で一般的な三角形メッシュ上で直接グラフ畳み込みを実行可能とすることで、物理シミュレーションと深い親和性を実現した。ランダムな2次元形状まわりの流れ(2,000例)をCFDで計算・データセット化し、これを学習に利用。未知形状に対する速度と圧力場の再構築精度は従来のU-netと比べてパラメータ効率面で優れており、特に境界層付近で高精度であった。GCNNで再計算した流速・圧力場によりドラッグフォース(抗力)も高精度で積分計算でき、GCNNの有用性を示した。

GCNNと解析のリンク方法:
CFD計算で使われる三角形メッシュのノード(頂点)とエッジ(辺)を、グラフのノード・エッジと対応させている。
GCNNはノード間の隣接関係を用いて局所的な特徴抽出を行っている。つまり、CFDでの有限要素法や有限体積法と親和性が非常に高い計算構造である。
メッシュをそのままGNNの入力データとして使うので、CFDで生成された流れデータ(速度・圧力場)をそのままGCNNに投入可能。これにより物理的形状や境界条件、流れの詳細構造を損なうことなく、ニューラルネットモデルによる高速な予測・再構築が可能となっている。

インプット:
・ノード座標(x, y)
・障害物境界情報(どのノードが障害物境界かを表すマスク)
・(場合によって)周囲のメッシュ情報や初期値など

入力は特徴量行列としてGCNN内に与えられる(ノードでは座標やマスク、エッジではノード間の距離など)。

アウトプット:
流速ベクトル u, v(x, y方向速度)、圧力 p

入力された障害物の形と位置などの情報から、GCNNはメッシュノード上の速度・圧力場の空間分布を出力し、その結果がCFDの計算値と合致するよう学習する。予測は未知形状に対しても一般化が可能である。

解析フロー:
ランダムな2D形状をBézier曲線で生成
それぞれの形状に対してCFDで速度・圧力場を計算(三角形メッシュ上に投影)
入力として「三角形メッシュ+座標」を用意
グラフ畳み込みブロックでメッシュ上の情報を伝播・集約・平滑化
各ノードでのu, v, pを出力(流れ場として再構築)
CFD値とのMAEで学習
学習済みモデルで未知形状の流れも高速予測可能

GNNには以前から興味があるのですが、まだ触ったことがありません。感覚的にはどのようなことをしているのか理解できるのですが、具体的に動かすとなると勉強が必要です。
学習用の解析を2000程度実施しても元が取れるというのは、それなりのモノづくり分野ならではでしょう。3次元浸透流では使いづらい、というのは以前にどこかで書き残したと思います。


花崗岩の亀裂、鉱物、引張強度

鉱物分布を考慮したき裂発生・進展解析による圧裂引張試験時の花崗岩の破壊挙動評価 

AI要約

本研究は、花崗岩を対象に鉱物分布を考慮したき裂発生・進展解析(数値解析)により、圧裂引張試験時の破壊挙動を評価したものである。鉱物の分布や含有量が破壊の進展や亀裂生成にどのような影響を及ぼすかを、実験および有限要素法(FEM)によるシミュレーションを組み合わせて解析している。実験データと解析結果を比較し、鉱物分布の違いが破壊挙動に与える影響について議論している。

鉱物含有量・分布:
花崗岩は複数種の鉱物(主に石英、長石、黒雲母など)から構成される。解析に用いた試料は、鉱物分布をデータ化し、各鉱物相ごとに力学特性の違い(ヤング率、破壊強度など)を設定したメッシュモデルを作成している。
REVs(代表的体積要素)ごとに鉱物の割合・配置を反映した力学パラメータを設定し、現実の鉱物分布パターンをモデル化している。

鉱物分布と亀裂の発生・進展との関係:
亀裂(き裂)は、力学的弱部や鉱物境界、脆弱な鉱物が局所的に多く分布する領域から発生しやすい。解析モデルでも鉱物分布に応じて力学的性質が変化するため、亀裂は鉱物境界や低強度鉱物(黒雲母など)近傍から優先的に発生・進展する傾向が出ている。
応力場分布(FEM解析結果)からも、鉱物分布の不均一性が亀裂進展経路に大きく影響し、均質なモデルと比べて破壊開始点や亀裂経路が異なる現象が観測された。
圧裂引張試験の結果、亀裂が主に鉱物の境界や低強度鉱物を通過しやすかった一方で、石英や長石など高強度鉱物で一時的に亀裂の進展方向が変化する様子も見られた。
亀裂発生・進展の予測精度は鉱物の分布データの正確性や力学特性のパラメータ化に依存している。破壊強度や変形特性のばらつきも鉱物分布の違いに起因。

これもピース、かつ引っ張り(圧裂)状態です。亀裂は解析から追いかけているようです。どのように実問題へ発展させていくか、は先の文献と同じで難しいところです。

花崗岩の亀裂、鉱物、圧縮強度

3D analysis of mineralogical and fracture features in granite and impact of biotite distribution on uniaxial compressive strength - ScienceDirect

AIによるまとめです。
最近、特に気になる文献や図書以外はAI任せになりつつあります。ますます英語から離れていきます。

強度と亀裂の基本的な関係:
 花崗岩の鉱物学的異方性(不均一性)と亀裂特性を理解することは、花崗岩の強度と亀裂の進化を評価する上で不可欠。ここでの鉱物学的異方性とは、鉱物の含有量、粒径、粒形、空間分布、粒界特性のばらつきを指し、これが岩石の機械的挙動(機械的強度や亀裂挙動を含む)に影響を及ぼす 。岩石の亀裂発展は以下の8段階。

  1. 亀裂閉鎖
  2. 線形弾性変形
  3. マイクロクラック発生
  4. マイクロクラック成長(安定した亀裂伝播)
  5. マイクロクラックの合体によるマクロクラック発生
  6. マクロクラック成長(不安定な亀裂伝播)
  7. マクロクラック合体
  8. マクロな破壊とポストピーク軟化
荷重の初期段階では粒間亀裂が開始し、より高い応力レベルで粒内亀裂が出現することが確認さた。また、初期荷重段階では、亀裂は黒雲母が豊富な領域に向かって伝播する傾向があることが示された。
亀裂体積比の急激な増加は、試料が不安定な破壊段階に入り、破壊に近づいていることを示唆する 。

鉱物含有量と強度:
先行研究では、鉱物含有量とUCSの間に関係があることが議論されており、特に黒雲母含有量の高さは一般的にUCSの低さと相関があるという負の傾向(逆相関)が報告されている。しかし、本研究の試料では、最も高い黒雲母含有量を持つにもかかわらず、最も高いUCSを示しているものがあり、一般的な傾向に当てはまらないケースとして挙げられる。これは、鉱物の単なる含有量だけでなく、その「分布」が重要であることを示唆している 。

鉱物分布と強度:
鉱物の空間分布は、岩石試料のピーク強度と亀裂発生応力の両方に影響を与える 。本研究では、特に「黒雲母」の分布がUCSに与える影響を分析。「中央領域に黒雲母のような弱い鉱物がより高濃度で存在する花崗岩は、単軸圧縮下で比較的高い強度を示す可能性がある」という仮説が立てられる。

人手を割いていた亀裂や鉱物のセグメンテーションに機械学習を使って省力化するなど、近年の基礎技術を利用した上での報告になっています。(先日聞いたある学会の災害報告会では、国の研究者が未だに古典的な手法の利用を堂々と話されていました。日本も追いついて欲しいものです。)
とはいえ、実務で欲しいのはピースではなくマスとしてのな強度なので、黒雲母のばらつきが対象スケールで不均質と言えない限りは応用の難しい内容です。まだ基礎研究レベルなので、本質的な部分は日本も追いつけるかな?


2025年8月7日木曜日

区間予測

ある時系列データに対し予測区間を求め、異常値が含まれているかどうか判断したいというシチュエーションは多くあります。

例えば、地下水に対する施工の影響。施工前の水位観測データで実行雨量などから簡単なモデルを作成しておき、施工中の水位低下について施工の影響が含まれているかどうか(予測区間を外れているかどうか)を判断したい。このような場合に多用されてきました。

実行雨量から機械学習になっても、考え方は同じです。
区間を作る方法は以下のように複数あります。他にもありますが、分布を仮定せず、早いのは 時系列拡張型の conformal prediction でしょうか。

① ブートストラップ
  • 元データから復元抽出で複数セット作成 → モデル再学習(複数回)で予測分布や区間を構築。
  • 理論背景がシンプルでモデル非依存(どのようなMLでも適用可能)。
  • 分布の仮定をしない分、推定精度はデータ量に強く依存する。
  • 繰返し計算で高コスト/遅い。
② 分位点回帰
  • 目的変数の特定の分位点(例:中央値、95パーセンタイル)を出力する回帰。
  • 分布の仮定不要で非線形構造にも強い。
  • 多数実装GBDT系では実装あり。
  • 分位点ごとに個別学習が必要でコスト増(例:2.5%, 97.5%予測には2モデル)
③  conformal prediction
  • MLモデルでの予測後に別途実施するため、モデル問わず「後付け」可能。
  • 絶対誤差 もしくは正規化絶対誤差を順に並べ、閾値を決めることで区間を予測(回帰) 
  • 誤分類率や区間内率など「理論保証」される。
  • 時系列の取り扱いやデータ分割に工夫が必要。

よくよく考えると、区間を連続して外れると異常となるわけですから、教師ありの異常検知とも言えます。
先日、水位データの異常値を除く方法を問うているプロポがありましたが、こういった方法も使えると思います。

2025年7月26日土曜日

雨量観測所の選定

先日、プロポーザルの特定テーマとして、「がけ崩れ箇所に近接する雨量観測所の選定方法」が問われていました。

降雨確率規模を算出するためにアメダスのような古いデータを使いたい、しかしどこを選んだらよいのかわからない、ということだと思われます。

アメダスでは昭和49年11月1日より観測が始まっています。しかし、観測所が約17㎞間隔で整備されているため、過去の土砂災害に対する降雨トリガー条件を設定する際等に、任意位置での雨量を推定するには、近隣の観測所を何らかの方法で選定し代用せざるを得ません。
このような観測所の選定問題に対し、古典的にはティーセン法を用いたり、周辺観測値の平均を使ったりしている例が多いような気がします。また、「地理的に近接」するとともに、同流域・同斜面等「気象条件が類似」する等の条件を設けて雨量観測所を選定することでも問題を解消できると期待されてきました。QGIS用の自動抽出ツールも提供されていますので、かつては広く受け入れられた方法だったのでしょう。

一方、近年では解析雨量の1㎞メッシュ化により、発生個所での降水量を得やすくなっています。対象箇所を含むメッシュの解析雨量とその周辺観測所を含むメッシュの解析雨量を直接比較し、降雨パターンの似ているメッシュ(観測所)を選定することが可能となっています。つまり、降雨データをベクトルとして扱い、それらの類似度を評価できるので、選択に迷うことがありません。機械学習やデータサイエンスで頻繁に使用される「データ駆動型」のアプローチに基づいた決定方法です。

このように、現在ではデータを根拠として観測局を直接選択することが可能となっています。あえてデータを利用せず、条件や仮定を設けて推定したり、それらの検証を省く必要はありません。が、これをテーマにされているということは、サイエンスとしてではなく、実績のある経験的手法の中から選択したい、ということだったのでしょう。
社会がデータ駆動のさらに次の段階に向かうには、まだまだ時間がかかりそうです。


2025年7月9日水曜日

論文 DASによる震源決定法

Source location of volcanic earthquakes and subsurface characterization using fiber-optic cable and distributed acoustic sensing system | Scientific Reports

震源決定

・2手法の併用
・両者を使用して震源位置を決定し、最も可能性の高い結果を採用。
・解析ではS/N比が4より大きい波形データを使用。
・グリッドサーチでは、すべてのグリッドポイントの到着時間差を事前に計算しておくことにより、ほぼリアルタイムで震源決定可能。

1.到着時間差法
道路のカーブを利用し、L字型アレーとして解析。
入射地震波の遅延性(伝搬方向と入射角)を推測。
30.6〜71.4mの組み合わせを利用。コヒーレンスの高い波を利用
到着時間差は、クロススペクトルの位相差より計算。
クロススペクトルは主に振幅の大きい波から計算されるため、到着時間の差はS波の到着と関連していると仮定した。
波形に12秒の時間ウィンドウを適用。
グリッドサーチを利用。

2. 振幅ソース位置(ASL)法
光ケーブルに沿った地震波の最大振幅を測定し、火山性地震の発生源を特定。
ASL法では、発生源から発生する地震波は、幾何学的な広がりと固有の減衰によって減衰すると仮定。
サイト増幅係数は、地域性構造地震のコーダ波の解析から決定。
火山性地震の最大振幅を光ファイバーケーブルに沿ったすべての測定点で読み取る。
Qを10〜100とすると、結果は大きく変化しない。20を使用。
グリッドサーチの利用。

到着時間差の測定
光ファイバーケーブルは長距離をカバーするため、到着時間差を計算する測定ポイントのペアの間隔距離には多くの選択肢がある。
最大間隔距離をS波の波長(4Hz)の約4分の1に設定。
ケーブルに沿って90m以内に位置する一対の測定点の波形を解析。
波形の相関を調べるために、コヒーレンスが0.5より大きい波形を使用。
ケーブル方向が異なる測定ポイントでの動的ひずみ信号は極性が異なる場合があるため、ケーブルの方向が互いに30°未満で一致している測定ポイントのペアのみを使用。
ケーブルに沿って91.8m以内にある測定ポイントのペアについて、9つの到着時間差を計算した。

次の手順を適用して、どの間隔が適切であったかを判断。
(1)10.2mから91.8mまでの9区間を10.2m刻みでデータを用意し、
(2)最適な位置を取得し、各区間の分散減少を評価。
(3) 分散の減少が10%未満の区間のデータを削除。
(4)その後、ステップ(2)からこのプロセスを繰り返した。
(2) から (4) までのプロセスは、差異削減が10%を超えたときに終了。
30〜70mの間隔を決定。

部位増幅率の推定
直達S波に続くコーダ波は、不均質な構造が地震波を十分に散乱させると、空間的に均一に分布する。
サイト増幅係数は、基準局に対するターゲットサイトにおけるコーダ波の相対振幅より推定。
約50〜200 kmの距離と約10〜100kmの深さで発生するマグニチュード3以上の地殻変動地震を11回利用。
ケーブルに沿った測定点における二乗平均平方根の振幅を計算し、基準点(地震計が仮的に設置されたMP1027)の振幅と比較。
S波走時の2倍以上の経過時間を持つコーダ波を使用。
5秒ごとに0〜50秒の経過時間の相対振幅を推定し、各測定ポイントでの平均振幅を分析のサイト増幅係数として使用した。
相対振幅の標準偏差は非常に小さく、サイトファクターが確実に評価されていることを示す。
サイト増幅率は、光ファイバケーブルに沿って対数スケールで約1のオーダーで変化するため、ASLではこの補正が極めて重要である。


2025年6月25日水曜日

ambient noise cross-correlation の前処理

 後輩君に ambient noise cross-correlation の前処理について聞かれました。

そういえばまとめていなかったなと思い、あらためて文献を整理。まとめてみると、1つの論文に収斂するのがよく分かりました。


No.

Title

Down Sampling

Demean

Detrend

Band Path

交通振動の少ないものを選定

amplitude-versus-offset metric

One-bit Normalization

Running-Absolute-Mean
Normalization

Spectrum Whitening

Correlation

1

岩手県の Hi-net 観測点で観測された常時微動の地震波干渉法解析による群速度の推定

20Hz

 

 

 

 

 

 

2

既設地震計の微小振動記録への地震波干渉法の適用による 農業用ダム地震波伝播特性評価の試み

 

 

0.5Hz-

 

 

 

 

 

 

3

山崎断層におけるトラフィックノイズを用いた地震は干渉法の適用

 

 

 

 

 

 

 

 

 

 

4

常時微動の相互相関関数中の実体波の特徴

 

 

 

0.52.0 Hz

 

 

 

 

 

 

5

地震波干渉法 その 1 歴史的経緯と原理

 

 

 

 

 

 

 

 

 

 

6

地震波干渉法 その 2 応用

 

 

 

360s, 36時間分スタック

7

地震波干渉法による西日本の地殻速度構造(1

 

 

 

 

 

 

 

 

いくつか試した結果、データ長1時間、MaxLag100

8

中京堆積盆地における表面波群速度の推定-Hi-net連続地震観測記録を用いた地震干渉法に基づく検証-

20Hz

 

 

0.052.0 Hz

 

 

 

 

1時間、1年分スタック,MaxLag100

9

福島県の広帯域リニアアレイで観測された常時微動の地震波干渉解析

20Hz

 

 

 

 

 

1sec

 

15分、1ヵ月分スタック(1時間、1日、10日、6カ月と比較)

10

Identification of a nascent tectonic boundary in the San‑in area, southwest Japan, using a 3D wave velocity structure obtained by ambient noise surface wave tomography

 

 

 

 

 

 

 

 

 

10-min segments with a 50% overlap1ヵ月分スタック

11

Global propagation of body waves revealed by cross-correlation analysis of seismic hum

 

 

 

 

 

 

 

 

 

2.8 h segments with a 50% overlap7年分スタック

12

Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements

 

 

 

 

13

Rayleigh wave group velocities at periods of 623 s across Brazil from ambient noise tomography

1Hz

760 or 360 s

 

 

 

 

14

Surface wave group velocity in the Osaka sedimentary basin, Japan, estimated using ambient noise cross-correlation functions

 

 

0.052.0 Hz

 

 

 

10sec

30分、50% overlap1年分スタック

15

Two-receiver measurements of phase velocity cross-validation of ambient-noise and earthquake-based observations

1Hz

 

 

2200s

 

 

 

 

half an hour with 50 per cent overlap

16

Aquifer Monitoring Using Ambient Seismic Noise Recorded With Distributed Acoustic Sensing (DAS) Deployed on Dark Fiber

125Hz

0.002-15Hz

 

 

 

0.5sec

 

1min

17

Distributed Acoustic Sensing for Seismic Monitoring of The Near Surface: A Traffic-Noise Interferometry Case Study

 

 

 

 

 

18

Daily Groundwater Monitoring Using Vehicle-DAS Elastic Full-waveform Inversion

 

3-10Hz

 

 

 

 

 

 


**************************************
20250628 文献追加
18:DASの場合はゲージ長が空間的なlow-pass filter の効果がある。