2024年11月4日月曜日

Landslide Susceptibility Map using ML その5

Full article: Near real-time spatial prediction of earthquake-induced landslides: A novel interpretable self-supervised learning method

こちらは transformer を利用されています。
pre-training 後に fine-tuning を実施したところ、他の手法より AUC がよかったよ、という報告です。

transformer ではグローバルな発生データを活用する pre-training が可能であり、データ量の多さを活かして高度な特徴学習ができそうです。それをローカルの発生・非発生データで fine-tuning することで、未知のデータに対する汎化性能を保ちつつ、ローカルな特性を捉えたモデルを構築できます。ローカルのために世界のデータを利用するという報告は"その2"で書き残しました。が、この場合は XGBoost ですので pre-training の概念がないですし、学習時に利用するにしても相応の非発生データが必要になります。

利用する特徴量を決めておいて、世界でデータを整備しておけば、事前学習済みデータとして配布・利用できそうです。幸か不幸か国内はこれからですので、特徴量に使える国内データの整備が進むとありがたいですね。


2024年11月3日日曜日

Landslide Susceptibility Map using ML その4

Full article: An integrated neural network method for landslide susceptibility assessment based on time-series InSAR deformation dynamic features 

時系列 DInSAR を特徴量として使用されています。変動量はSBASでもチェックをされているようです。SARというと、つい地震前後の差分をイメージしてしまいますが、地震前の変動量を利用することは言われて初めて気づいた重要なポイントですね。

24 stages of time-series InSAR cumulative deformation information are taken every 96 days per quarter. 

国内で Landslide Sasceptibility Map を産総研さんの研究以外で聞いたことがありません。それを作成するための機械学習も土木分野では浸透していません。あっても SVM とか Random Forest などクラシカルな手法が使われているように感じます。が、この分野の研究では複雑な特徴に対応するため、または精度向上のために DNN 等が利用されています。

Based on traditional linear statistical analysis, machine learning methods stand out by virtue of their ability to examine large amounts of data independently. Machine learning methods, such as random forest (RF) (Dou et al. Citation2019) and logistic regression (Zhang et al. Citation2019), have been extensively used in LSA. However, under the requirements of complex scenes or high precision, traditional machine learning algorithms cannot meet actual demand (He et al. Citation2021a; Zhao et al. Citation2022). Building on the neural networks present in machine learning, the neural network method effectively predicts complex nonlinear dynamic systems. It has been widely and successfully introduced into the field of LSA, including the convolutional neural network (CNN) (Wang, Fang, and Hong Citation2019; Gao et al. Citation2023a), recurrent neural network (RNN) and deep belief network (DBN) (Chen et al. Citation2020). The convolutional layer of CNN can extract multidimensional features from the input images and has good performance (Hakim et al. Citation2022). Gated recurrent unit (GRU) network of RNN variant has good performance in processing sequence data (Zhao et al. Citation2022). With the complexity of the environment, when faced with a limited sample, the ensemble learning model is also widely used in the LSA. For example, Wang et al. (Citation2022) conducted the LSA based on the XGBoost ensemble learning model. Lv et al. (Citation2022) combined CNN, DBN and ResNet models with the ensemble learning techniques of Stacking, Bagging and Boosting to generate the LSA.

日本は2ステップ遅れている状況です。せめて LSM 作成のためのデータは整備していただきたいものです。そうすると土木に携わっていない機械学習エンジニアが参加しやすく、多様な目的のマップが作成されるようになるでしょう。その時々の、最新のアーキテクチャで。

2024年10月27日日曜日

Landslide Susceptibility Map using ML その3

 昨日、産総研さんの第41回地質調査総合センターシンポジウムを拝聴しました。

地質に関するCIM やデジタル化の話、九州のハザードマップ作成についての報告がありました。

ハザードマップは産総研の方の報告でした。その中には LSM の話もありました。産総研では LSM を「地すべり感受性マップ」と訳されているようです(が、この訳は一般的ではないとおっしゃっていました)。スライドには「崩れやすさマップ」という記載もありましたので、今後はわかりやすい後者に落ち着くのかもしれません。

LSM 作成の話がメインのように思えましたが、表題はハザードマップで、しかも使用されているインベントリは被災を伴う災害データベースであり、作成されたマップはリスクマップに近いようでした。「感受性とハザードはどう違うのか」と質問された方がいらっしゃいましたが、「いつ、どこで、まで答えるのがハザード」というような回答でした。また、災害データを使われた理由については、「被災の有無にかかわらない崩壊データも一部で使っている。ないものは使えない」といった回答でした。
それをどうするか、は研究者ではなく技術者の仕事、でしょうか?Googlig では産総研さんからのリスクマップ作製に関する委託業務がいくつか引っかかりますので、初めからLSMでなくリスクマップを作成する目的だったのか、単に予算の問題だったのかもしれません。
九州地域における斜面災害リスク評価主題図改良業務:SVMを利用
九州北部地域の地質情報解析によるリスクマップ(案)の作成業務:災害特性に着目


国内では崩壊データを整備する仕組みがありません(国が災害データを収集する仕組みはあります)。災害データを使うと、民家に近いところが危ない、危険区域や危険個所に指定されている範囲が危ない、という結果になります。これは当然で、あちこち崩壊しても被災した箇所しか国には報告されない、報告外の箇所は未崩壊として誤って扱われてしまいやすいというのが理由です。そのインベントリを教師データに使用すれば、属性としての地質よりも区域等が効いてきます。機械学習モデルの中で地質の重要性が薄れてしまうのです(地すべりはそれでも地質が効くのでかなり重要なのでしょう)。

では、LSMを作成するのに災害データしか入手できない場合はどうするか?

①航空写真、衛星写真等から作成する。
産総研さんの「一部」というのはこれかもしれません。時間と費用がかかるものの、現状では最も精度の高い方法と言えるでしょう。写真等から抽出するための画像処理や機械学習が進化すれば解決すると思います。が、現状では最後に人手が必要です。もう少し待ちましょう。

②文献で利用されているデータを利用する。
昨日の文献通りです。ローカルの評価にデータが不足すれば、世界中から集めてくればよいのです。文献個々のサプリメントデータを集めるのも有効です。機械学習にかける際には多少の工夫が必要になると思いますが、データ不足には効果的です。

③LSMをあきらめてリスクマップを作成する。
この場合、LSMを作成したうえで、ある地域に台風〇〇号並みの豪雨が来たらどうなるか、その際に被災しやすい場所はどこか?といったような手順は踏めません。被災関係なく豪雨や地震で崩れそうな場所を予測することも正確にはできないのですが、被災しやすい箇所は推定できます。

産総研さんの報告からは、①と③のいずれの方法をどの程度採用しているのかは明確ではありませんでした。聞く限りでは作成されたハザードマップ?リスクマップ?よりは収集された属性データを公開していただく方が、後々利用目的に応じて柔軟なマップ作りができるように感じました。

2024年10月25日金曜日

Landslide Susceptibility Map using ML その2

2019年~2024年まで、機械学習を利用したLSM作成に関する文献を20編程度集めました。

  • 文献で使われている手法は、RF や SVM からスタートし、DNN が現れて、2022年ごろから XGB や LightGBM が加わっていました。
  • 2022年までは ML の結果のみで LSM を作成していますが、その後は何かと組み合わせるパターンが加わっています。
  • imbalance data への対応は under samplimg が主体。1:1~1:2が良い結果を導いているようです。
  • 地形スケールの影響を検討している文献もちらほら。30m程度が多いように感じるものの、何が良いかはケースによる、といった文献がありました。
  • ブロック全体からポイントを作るか、中心から作るか。滑落崖からつくるか移動体から作るか。微妙な差でしたが滑落崖全体から作るのが良かったようです。

気になった文献のメモ( ..)φ
データ構成のフローチャート Fig. 6 がわかりやすいし、全体構成が参考になります。
Important considerations in machine learning-based landslide susceptibility assessment under future climate conditions | Acta Geotechnica

  • undersampling (one-sided selection, k-means clustering, gridded hyperspace even sampling, random sampling)の比較。グリッド型ハイパースペースサンプリング手法が有効
  • XGBoostモデルの外挿予測の確認。

 Even though the case study presented in this work is implemented for the state of California, using only California-based datasets for model training may pose several challenges. For instance, the number of positive data samples may be insufficient, leading to a severely imbalanced dataset. In addition, given the changing climate conditions, developing a robust model for future landslide susceptibility assessment requires a diverse dataset not limited to the historical California data to prevent issues with extrapolation due to the impacts of climate change. To address these challenges, a global dataset is used here for training to ensure adequate data samples and reduce extrapolation errors. 

Research indicates that tree-based models, particularly XGBoost [24], are robust to correlated features, and variables do not necessarily need to be removed due to high correlation. However, to reduce model dimensionality, only one out of each two highly correlated variables is kept in the model.

機械学習をかじっていると転移学習が有効であることを知っているので、この文献の似た考え方は受け入れやすいでしょう。データを利用させてもらえると、日本の LSM は作成しやすくなります。

2024年10月22日火曜日

機械学習による水位予測 その2

 ”water level machine learning” + Googlimg で多数引っかかる内容です。

頭から読んでみたところ、地すべりのSusceptibility Map の作成の文献とは傾向が異なっており、アルゴリズムの比較が多く目につきました。データセットによって予測性能の良いアルゴリズムが異なる点はよく知られていますが、見事にバラバラ。それらのアンサンブルが少ないのは機械学習分野の流れと異なりますが、浸透するにはもう少し年月が必要なのでしょう。

内容的には玉石混合ですが、アルゴリズムの比較という点では有用ですので、引っかかったものから書き残しておきます。

Exploring machine learning algorithms for accurate water level forecasting in Muda river, Malaysia - PubMed

  • DNN, LSTM, XGBoost
  • 1日前の水位を利用する場合はDNNがBEST。
  • 7日間の予想ではLSTMが良好。
  • データの量と質にに大きく依存する。

Deep Machine Learning-Based Water Level Prediction Model for Colombo Flood Detention Area

  • DNN, LSTM
  • Daily rainfall, Daily evaporation, minimum daily temperature, maximum daily temperature, daily relative humidity at daytime/nighttime, and daily average wind speed
  • LSTMの方が良好

Prediction of Water Level Using Machine Learning and Deep Learning Techniques | Iranian Journal of Science and Technology, Transactions of Civil Engineering

  • Random Forest, XGBoost, RNN, BiLSTM, CONV1D-BiLSTM
  • XGBoostがBEST

Water level prediction using various machine learning algorithms: a case study of Durian Tunggal river, Malaysia

  • 線形回帰 (LR)、相互作用回帰 (IR)、ロバスト回帰 (RR)、ステップワイズ回帰 (SR)、サポートベクター回帰 (SVR)、ブーストツリーアンサンブル回帰 (BOOSTER)、バッグドツリーアンサンブル回帰 (BAGER)、XGBoost、ツリー回帰 (TR)、ガウス過程回帰 (GPR)
  • 29年間の日降水量を使用。
  • 自己相関関数によるラグタイムを考慮した4つの入力データ(シナリオ)で試行
  • 過学習対応のため、Lossをearly stoppingに使用
  • 異なるフォールド(3、5、7、9)を使用
  • GPRがBEST
  • 精度を高めるには長いラグデータが必要
  • 過学習の問題は通常、小さなデータセットを使用してモデルを開発するときに発生するが、1990年から2019年までの日次データの長さはそのような問題を回避するのに許容できる長さであった
  • 不確実性分析を実施

Modeling the fluctuations of groundwater level by employing ensemble deep learning techniques

  • DL, アンサンブルDL
  • 4 つの井戸の地下水位から残りの井戸の地下水位の予測
  • 最大 20 日前の地下水位記録の時系列を入力として5つの井戸の翌日の地下水位を予測
  • 合計レコード数 276、2017 年 10 月 20 日から 2018 年 5 月 1 日までトレーニングセット (70%、194 レコード) 、2018 年 5 月 2 日から 2018 年 7 月 22 日までテストセット(30%、82 レコード)
  • EDLの方が良好

2024年10月21日月曜日

機械学習による水位予測

2016年に DNN を利用した水位予測を試行していました。
https://phreeqc.blogspot.com/2016/11/using-deep-learning-4.html

それから8年が経過し、ようやく地下水チームからも機械学習による予測に興味を持つ方が現れ始めました。
8年の間に、"その3"で書き残していたように、いくつかのラグを有する特徴量を利用しないとうまく機能しないことがわかっています。また、決定木系の回帰では、渇水年など外挿となる予測が困難なことも体感しました(コチラの文献を見ると理解しやすいと思います)。

国内では Signate で自治体主催の水位予測コンペが開催されたり、関連文献も見かけるようになりました。が、国外のほうが圧倒的に進んでいると感じています。国内では物理ベースのモデル構築が好まれて、海外よりもより保守的だったのかもしれません。物理モデルの構築にはいくつかの仮定が必要であり、問題を包括的に理解しシンプルに表現する力量が必要です(昔の人は偉かった)。この力量が問われてきた世界では、仮定が必要なく、説明性も低い機械学習を学問として受け入れるのに抵抗があったのかもしれません。が、大量のデータを容易に得られる現代では、そこから新たな知見を発見するというデータ駆動型のアプローチも重要なんですよね。

関連する内容が、以下に書かれています。
Full article: Water level prediction using various machine learning algorithms: a case study of Durian Tunggal river, Malaysia

For example, in Malaysia, a physical-based model was developed to assess one river's floodplain and water level (Mohamad et al., Citation2014). In order to build such a physical model, there was the need to collect massive data and information on top of the cost involved in building such a physical model. Over time, numerical models overcame the limitations of the physics-based models. For instance, a numerical model was developed by (Wu et al., Citation2014) to forecast water levels at the Yangtze River. However, a study conducted by (Guan et al., Citation2013) reported that many errors were encountered during the development of a numerical model. Therefore, despite the noticeable improvements in the numerical models, they still have limitations, such as the need to mimic some of the physical phenomena to improve their accuracy and reliability.
Recently, data-driven techniques were shown to have overcome traditional models’ drawbacks and proved to be more accurate in modeling complex engineering problems. 

ま、データがたくさん必要なのは機械学習も同じです(いえ、それ以上)。が、仮定の必要がなく、ハイパーパラメータの調整さえすればあとは自動でモデルを構築してくれる、ほぼリアルタイムで予測結果を算出できる、しかも予測性能が高くなるとなれば、使わない手はないでしょう。合理性の観点から国外で進んだのかもしれません。

続く。


2024年10月20日日曜日

水位予測に使用する特徴量

地すべり速度、水位の予測性能を Random Forest と重回帰で比較した文献です。
A comparative study of random forests and multiple linear regression in the prediction of landslide velocity | Landslides

RFの回帰の仕組みとしては Fig8 が直感的に理解しやすいと思います。この木が多数あって平均をとるのがRF、重みを含めているのが勾配ブースティングのイメージです。線形補間と異なり階段状の予測モデルになるので、ある程度深い木 or 木の数が必要、それで過学習にならないためには多数のデータや特徴量が必要、ということになります。

sckit-learn の cheat sheet では、50以下でデータを集めなおし、10万以下で線形SVRなど、さらに必要なら非線形やアンサンブル手法(RF、勾配ブースティングなど)という流れになっています。
12. Choosing the right estimator — scikit-learn 1.5.2 documentation

この文献の良いところは入力データを明示しているところです。地すべりの速度を予測する前に水位を予測するのですが、それに使用した75の特徴量をラグ付きで示されています。文献によってはどの程度のラグを使ったかが明確にされていない場合もあるので、表1のように整理して表示されるとありがたいです(理解できない内容もありますが)。

2024年10月18日金曜日

Python から NASへアクセス

Python から NAS へアクセスする必要が出てきました。

Windowsであれば問題にならないのですが、Ubuntu だとマウントが必要です。
CIFSを利用します。コマンド打ちの方が圧倒的にシンプルです。

$ sudo apt install cifs-utils
$ sudo mount -t cifs //[nas_ip]/[nas_share] /mnt/[name] -o user="user_name",pass="nas_password"

どうしても Python スクリプト内で実施する必要がある場合は以下。

import subprocess
#マウントディレクトリを作成。
mount_point = "/mnt/[name]"
subprocess.run(["sudo", "mkdir", "-p", mount_point],
                input=[sudo_password].encode(),
                stdout=subprocess.PIPE,
                stderr=subprocess.PIPE)
#CIFS共有をマウント
mount_cmd = ["sudo", "mount", "-t", "cifs",
              f"//{nas_ip}/{nas_share}",
              mount_point,
              "-o",
              f"user={user_name}, pass={nas_password}"]
subprocess.run(mount_cmd,
               input=[sudo_password].encode(),
               stdout=subprocess.PIPE,
               stderr=subprocess.PIPE)


2024年10月14日月曜日

Blind Source Separation

ブラインド信号源分離(BBS)は、2つのソースが混在した観測記録(波形)から、元のソースを分離しようという試みです。

音声であれば同じ媒体(空気)を通すので経路中の波形の変化は無視できそうですが、地震波の場合は2つの観測点まで同じ地盤を通るわけではないので、その変化を無視できそうにありません。その場合でも、どの程度使えそうか試してみました。

・ウェーブレット変換
・独立成分分析(ICA)
・非負値行列因子分解(NMF)

ウェーブレット変換は、簡単なsin波の合成でも振幅が大きく異なるとダメでした。
ICAはOKです。ソースの1つとして同じ地震波形を利用した場合、2つ目のソースの振幅が大きく異なっても奇麗に抽出できました。が、自然界にそのような都合の良いモデル波形はなく、実際の2カ所で観測された波形を使うとダメでした。NMFも同じくダメでした。

なかなか難しいようです。


2024年10月11日金曜日

ヘテロダイン検波

DAS(ヘテロダインOTDR)を理解するため、順に可視化してみました(Thanks! GPT)。
誤っていたらいずれ気づくでしょう。その際に修正しましょう。

まずは import。

import numpy as np
from scipy.signal import hilbert, butter, filtfilt
import matplotlib.pyplot as plt
 

レーザーを作成します。

# サンプリング周波数
fs = 5000  # サンプリング周波数
t = np.linspace(0, 1, fs, endpoint=False)  # 1秒間の時間軸

# 信号の周波数
f_signal = 1000  # 1000Hzの振動信号
signal = np.sin(2 * np.pi * f_signal * t)

# 周波数シフト
carrier_freq = 1010
carrier_wave = np.sin(2 * np.pi * carrier_freq * t)

# 後方散乱光
# 光ファイバーから得られた2ショット
# 位相が少しずれている
scattered_signal1 = np.sin(2*np.pi*(carrier_freq)*t
                   + 2*np.pi*0.1)
scattered_signal2 = np.sin(2*np.pi*(carrier_freq)*t
                   + 2*np.pi*0.2)
                   
plt.plot(t[:100],signal[:100],label="signal")
plt.plot(t[:100],carrier_wave[:100],label="carrier",c="green")
plt.plot(t[:100],scattered_signal1[:100],label="scatteres1")
plt.plot(t[:100],scattered_signal2[:100],label="scatteres2")
plt.legend()
 

戻ってきた光を干渉させます。
sin波の積和の公式より、高周波成分と低周波成分が生まれます。

 # 干渉信号の生成
interfered_signal1 = scattered_signal1 * signal
interfered_signal2 = scattered_signal2 * signal
plt.plot(t[:100],interfered_signal1[:100])
plt.plot(t[:100],interfered_signal2[:100])
 

光の高周波数は検出できませんが、干渉後の低周波成分は検出可能です。

 # Butterworthフィルタ(ローパス)
cutoff = 50
# 正規化周波数
Wn = cutoff / (f_signal / 2)
# Butterworthフィルタの設計
b, a = butter(4, Wn, 'low')
filtered_signal1 = filtfilt(b, a, interfered_signal1)
filtered_signal2 = filtfilt(b, a, interfered_signal2)
plt.plot(t,filtered_signal1)
plt.plot(t,filtered_signal2)
 

低周波側の信号から位相を取り出し、2ショット間の位相差を算出します。
ここはプロのノウハウがあるそうです。図書にも見当たらないので推定です。

 # ヒルベルト変換を用いた瞬時位相の抽出
analytic_signal1 = hilbert(filtered_signal1)
instantaneous_phase1 = np.unwrap(np.angle(analytic_signal1))
analytic_signal2 = hilbert(filtered_signal2)
instantaneous_phase2 = np.unwrap(np.angle(analytic_signal2))
# 瞬時周波数の計算
instantaneous_frequency1 = np.diff(instantaneous_phase1)
                                   / (2.0 * np.pi) * fs
instantaneous_frequency2 = np.diff(instantaneous_phase2)
                                   / (2.0 * np.pi) * fs

plt.figure(figsize=(10, 10))
# プロット:位相
plt.subplot(311)
plt.plot(t, instantaneous_phase1)
plt.plot(t, instantaneous_phase2)
plt.title('Extracted Phase from Interfered Signal')
plt.xlabel('Time [s]')
plt.ylabel('Phase [radians]')
plt.grid(True)

plt.subplot(312)
plt.plot(t, instantaneous_phase2-instantaneous_phase1)
plt.title('Extracted Phase from Interfered Signal')
plt.xlabel('Time [s]')
plt.ylabel('Phase [radians]')
plt.grid(True)

# プロット:周波数
plt.subplot(313)
# t変数の長さをinstantaneous_frequencyに合わせる
plt.plot(t[:-1], instantaneous_frequency1)
plt.plot(t[:-1], instantaneous_frequency2)
plt.plot(t[:-1], instantaneous_frequency2
                  -instantaneous_frequency1)
plt.title('Instantaneous Frequency')
plt.xlabel('Time [s]')
plt.ylabel('Frequency [Hz]')
plt.ylim(-20,20)
plt.grid(True)
plt.tight_layout()
 


2024年9月29日日曜日

Complex Hilbert PCA

DASでの震源決定に複素主成分分析を利用されている事例を拝聴しました。
確かに、密な観測点ですので相互相関は得られそうですし、ノイズも落としやすいと思います。

複素主成分分析についてはよく知らなかったので調べてみました。
時系列データに対し、ヒルベルト変換をかませてからPCAにかけるだけのようですね。手順としては以下の通りです。 (by GPTさん)

1.時間領域の波形読み込み

2. 信号の標準化
標準化(mean = 0, standard deviation = 1)することで相互相関の計算が比較的均一になります。
signals_mod = [(s - np.mean(s)) / np.std(s) for s in signals]

3. ヒルベルト変換
ヒルベルト変換を使って各信号の解析信号(analytic signals)を計算。ヒルベルト変換は信号を複素数表現に変換し、瞬時位相や包絡線を得るために使用されます。
analytic_signals = [scipy.signal.hilbert(signal) for signal in signals_mod]

4.エルミート共役行列の計算
解析信号のエルミート共役(Hermitian transpose)を計算。エルミート共役は、行列の転置を取ってから各要素の複素共役を取る操作です。
analytic_signals_h =np.conjugate(analytic_signals.T)

5. 相互相関行列の計算
解析信号の相互相関行列を計算。複素数表現を使用することで、信号の振幅変化だけでなく、位相変化も考慮した相互相関を計算することができます。
corr_matrix = np.dot(analytic_signals, analytic_signals_h) * (1/fs)
1/fs倍することで正規化しています

6.主成分分析
相互相関行列の固有値(eigenvalues)と固有ベクトル(eigenvectors)を計算。
重要な情報を保持しつつデータの次元を削減します。小さな固有値に対応する成分は、ノイズである可能性が高く、これらの成分を無視することでデータのノイズ除去が可能です。
eigenvalues, eigenvectors = np.linalg.eigh(corr_matrix)
principal_components = np.dot(eigenvectors.T, analytic_signals)

ここから主成分をいくつか抜き出して、ノイズを落とした信号を再構築すれば、前処理は終わりです。再度、相互相関でΔtを推定するなり、位相差を計算するなりして利用するのでしょう。

***********************************************
20241007追記

複素行列Aのエルミート転置イロイロ。
A.H
np.swapaxes(np.conjugate(A),1,0)
np.conjugate(A).T
np.conjugate(A.T)

.Hが簡単です。

2024年9月27日金曜日

Landslide Susceptibility Map using ML

2020年頃から機械学習が Landslide Susceptibility Map の作成に用いられるようになりました。今日では安定計算と同じような基本ツールの一つになっています。日本の土木分野は、かなり出遅れています。

初期の文献では、ある事象を単一モデルで学習し、予測マッピングするだけでした。その後、インバランスデータの取り扱いやアルゴリズムの選択などの文献が現れています。参考になるのは以下。1:2がBESTだった例です。
Comparative study on landslide susceptibility mapping based on unbalanced sample ratio | Scientific Reports (nature.com)

今日では既知の地震に対してのみならず、新たな地震に対しても利用できるような汎用性を持たせる作り方が報告されています。新たな地震に対する危険度や損害の程度、それがどこで発生するかを予測したいからでしょう。過去の地震をトレース、予測するだけだと予防や対策につながりませんので。

また、MLと別の分析結果を組み合わせる例が増えているようです。予測された危険度毎に建築物の数を集計するなど、家屋への被害の程度を簡易的に表現する単純な方法でも現段階では有用だと思われます。(建築物の数を考慮すると Risk Map になります。)
Enhancing co-seismic landslide susceptibility, building exposure, and risk analysis through machine learning | Scientific Reports (nature.com)

MLとKDEとのアダマール積、InSARの結果などを組み合わせて評価している例もあります。残念ながら、これらは汎用性に劣るため流行らないでしょう。先のReg3Dは後で組み合わせるより先に計算しておいてMLに投入する方が良いでしょう。どの範囲、というのはDEMだけだと難しいですから、確率をメッシュに与えた方が当たりそうな気がします。

MLと何を組み合わせるか、アイデアはたくさんあると思います。基本ツールの結果を組み合わせ、本来解くべき問題に早々に挑めると、少しは前進するのではないでしょうか。


2024年9月26日木曜日

Landslide Susceptibility Map using 3D LEM

RegionGrow3D: A Deterministic Analysis for Characterizing Discrete Three‐Dimensional Landslide Source Areas on a Regional Scale - Mathews - 2024 - Journal of Geophysical Research: Earth Surface - Wiley Online Library

In this study, the RegionGrow3D (RG3D) model is developed to broadly simulate the area, volume, and location of landslides on a regional scale (≥1,000 km2) using 3D, limit-equilibrium (LE)-based slope stability modeling. Furthermore, RG3D is incorporated into a susceptibility framework that quantifies landsliding uncertainty using a distribution of soil shear strengths and their associated probabilities, back-calculated from inventoried landslides using 3D LE-based landslide forensics.

動画もあります。
Landslide susceptibility map generated using RegionGrow3D (youtube.com) 

インベントリがあれば、そこから逆算してΦを出す、その確率密度の分布から累積分布を書いて susceptibility とみなす、その確率でどの大きさの地すべりがどこに分布するかを推定する。このような土質力学ベースの計算が可能な Matlab ツールをUSGSが提供しています。

文献では、すべり層厚を決めるのに隆起量などをパラメータとした別のシミュレーションを利用していますが、この辺りはインベントリの平均層厚でも構わないでしょう。

すべる・すべらないのバイナリを拡張し確率を持たせた点、thresholdに応じて地すべりの大きさを扱える点が魅力です。日本は完全に置いてけぼりです。


2024年9月23日月曜日

Susceptibility, Hazard, Risk

4~5年前から、海外の文献で “Landslide Susceptibility Map” をよく見かけるようになりました。

”感受性マップ”と直訳するとわかりにくいのですが、危険性や事象が起こりそうな確率を表す地図を意味しています。日本で使われるハザードマップが近いと思います。
が、Susceptibility と Hazard は使い方が異なります。前者は空間に対し、後者は時空間に適用されます。さらに、Hazard が発生した場合の影響度を考慮するとRiskになります(KY活動からリスクアセスメントに移った頃を思い出せば、容易に理解できると思います)。
このあたり、検索してみるとすぐに出てきます。例えば以下。

DOI: 10.5772/intechopen.100504
2. Definition and concepts
・Susceptibility refers to the probability of occurrence of an event within a selected type during a given location.
・Hazard refers to the probability of occurrence of an event within a selected type and magnitude during a given location within a reference period.
・Risk refers to the expected losses or damage by events during a given region, which are the products of susceptibility, hazard, and elements in peril.
・Vulnerability means the degree of loss to a given element of the set of elements in peril resulting from the occurrence of natural phenomena of a given magnitude. It is expressed on a scale from 0 (no damage) to 1 (total damage). 
・Elements at risk is potentially vulnerable of properties, population, and economic activities including public services in peril during a given area.

用語を正しく用いないと伝わりません。忘れないようにしましょう。

2024年9月22日日曜日

hDVS

An Introduction to Distributed Optical Fibre Sensors

先の洋書よりもハードウェアや原理について詳しく書かれています。先日使用した hDVS について、大まかな流れは理解できました(細かい点はわからないことばかりで先は長い)。

以下、流れの抜粋です。( ..)φ

狭帯域レーザーの出力は、Optical local oscilator (OLO)パスとパルス形成パスに分割される。

前者は角周波数 ω0 で放射され、ELo・exp(jω0t) の形の時間依存電界を持つ。

後者はAcoustic optical module (AOM) によってω0 + Δω のプローブパルスに周波数シフトされる。

サーキュレータ等を使用して、プローブパルスをセンシングファイバーに発射し、後方散乱光を OLO と混合する別の結合器を介してバランス型受信機に導く。

検出器に降りかかる電界は次のように表される。

Eωt = EL0・exp(jω0t)+Eb(t)・exp{j[(ω0 + Δω)t + Φ(t)]},  (3.29)

光検出器は表面に到達する光パワーに応じて電流を供給するので、光電流は次のように表される。

iphot = Rd{EL0^2+Eb^2+2EL0Eb・exp[j(Δωt+Φ(t))]}  (3.30)

バンドパスフィルタを使用して中間周波数IF=Δω のみを選択することで、OLO の電界と後方散乱信号の積であるヘテロダイン項ihet(つまり、式 3.30 の 3 番目の項) が得られる。これは、受信器の出力を支配する唯一の光学的寄与である。

ihet = Rd{2EL0Eb・exp[j(Δωt+Φ(t))]}  (3.31)

この信号を二乗することにより、後方散乱信号と直接比較できる波形が得られる。つまり、自由空間の各ポイントから返される後方散乱パワーに比例する。

ヘテロダイン OTDR では、バンドパスフィルタ適用時に後方散乱信号の振幅が検出され、位相情報が破棄される。対照的に hDVS では、後方散乱の位相は、位相情報を保持しながらも電子機器で処理できる周波数であるダウンシフト IF を介して測定される。(プローブパルスの周波数が50MHzシフトアップされている場合、ヘテロダイン後方散乱信号の周期は、約2mのファイバーに相当する。)

処理の手順
1.ファイバーに沿った距離として位相を検出。
2.レーザーのショット間の位相差を算出。
3.選択した距離間隔(ゲージ長)にて位相差を空間微分。
※ゲージ長全体にわたる位相の変化をレーザーショット数の関数として記録することも可能。これにより、ひずみの時間微分(ひずみ速度)が得られる。

これらの処理後、ファイバーに沿った距離と遅い時間(つまり、プローブパルスの連続に見合った時間スケール)の関数として、アンラップされた位相の画像が利用可能になる。


2024年9月15日日曜日

Distributed Acoustic Sensing in Geophysics

Distributed Acoustic Sensing in Geophysics: Methods and Applications 


昨年購入し、何度か読むにつれ徐々に頭に入ってきました。基本原理についても最初に記載があります。


1章

<A(z)> = 1/(A0·Fs)·τ(z)⊗[v(z+L/2)−v(z−L/2)] = 1/(A0·Fs)·τ(z)⊗ε˙(z)·L
A0 = λ/(4π·neff·Kε) = 115nm

The elongation corresponding to ΔΦ = 1rad is A0 = 115nm, calculated for λ = 1550, neff = 1.468  and Kε = 0.73, which has been measured for conventional fiber. 

The DAS signal is a convolution of pulse shape  with a measured field which is the spatial difference in fiber elongation speed of points separated by a gauge length.


2章

Φ(z,t) = atan(Q/I)
ε(z,t) = Φ(z,t)·λ/(4π·neff·Kε)/L

Gauge Length
・ゲージ大 → 帯域小, パルス大 → SNR大
・ゲージ小 → 地震信号の忠実度は高まるが、SNRが低下する。
・必要なデータ帯域幅のゲージ長を確認する必要がある。
・良好な信号を得るために、最適なゲージ長、パルス幅を選択する。

Fading
・Phase unwrapping algorithm
・フェーディングは DAS 取得の自然な特徴。常にハードウェアとソフトウェアの両方で対処する必要あり。
・複数の光周波数を使用する IU は、本質的にフェーディングの低減に優れている。
・取得後の処理でフェーディングに対処可。

Common-Mode Noise
・Stack
・周辺の環境音によって発生。
・簡単な取得後の処理で、そのほとんどを除去可能。

Spacial Calibration of Channels
・Tap test
・チャネルの位置を決定することは、結果として得られる精度にとって重要。
・既知の位置をキャリブレーションポイントとして使用し、インターリーブチャネルをそこから補間する。

その他
・可能であれば、マルチモード ファイバーではなくシングルモードファイバーが望ましい。
・接着されたファイバーケーブルを使用すると、より優れたSNRデータが得られる。回収可能なファイバー配置方法を使用して目的に適したデータを取得することも可能。
・光ファイバー接続を清潔に保つ必要あり。
・DAS データでタイミング情報が確実に保持されるよう、GPS への接続が必要。
・ファイバーの入射角応答がジオフォンの応答と異なる場合は、まずその影響を判断し、次にジオフォンのような応答を除去する。
・DASのSNRはジオフォンからのデータよりも本質的に低くなる。しかし、高密度で取得可能であるため大幅に改善可。


2024年9月13日金曜日

DAS その2

1年ほど前にこんなことを書いていたのですが、先週、DAS (Distributed Acoustic Sensing) に触れる機会がありました。

と言っても、私は作業員で路上を叩くだけ(Tap Test)。後輩君がメインで進めるお手伝いした形です。

まだ機器の構造と理論を理解できていないのですが、理解したいと思わせる結果を見せてくれます。橋の揺れ、ダンプの走行など数㎞の揺れを表示してくれるのが面白い。測定画面を見ているだけでも時間が過ぎるくらい、久々に興味を惹かれる機器でした。

今年の物理探査学会講演会でも1セッション設けられるほど注目されています。早く理解しましょう。以下、文献より備忘録です。

Seismic Monitoring With Distributed Acoustic Sensing From the Near-Surface to the Deep Oceans

φ = k·n·2z
Δφ = 4π/λ·Δn·Δz + 4π/λ·n·ε·Δz = 4π/λ·Δneff·Δz
Δneff = Δn + n·ε
ε˙(z) = [v(z+L/2)-v(z-L/2)]/L

where k is the wavenumber (k = 2π/λ, with λ the light central wavelength), n is the refractive index and 2z is the two-way length that the backscattered light travels. 


2段階認証

 google の2段階認証アプリに異常が発生し、他のアプリに切り替えました。

その間、約1週間。ひとまず会社とGoogleの2アカウントは復帰できましたが、それ以外は不具合の発生しているアプリの中で眠っています。ここにもようやくたどり着きました。

GitHub がまだ眠っています。ひとまずローカルからアップはできるのでしばらく待つつもりです。

怖いですね。

2024年8月18日日曜日

core MP135

core MP135 を購入。
Python コードを動かすまでの手順です。

balena-eacharを管理者モードで起動
microSDにイメージを焼く

USBケーブル接続後、シリアル通信
ID: root
Pass: root
$ passwd root

ユーザー追加
$ adduser "user name"
sudo group に追加
$ gpasswd -a "user name" sudo

root パーティション拡張
$ cd /usr/local/m5stack
$ bash resize_mmc.sh 

LANに接続後、SSH接続
IP確認
$ ifconfig

・自動起動させるプログラムをWinSCPでコピー。
・python3 で稼働確認。
・vimでrc.local を変更し、自動起動設定。


double-layer two-phase formulation その3

double-layer two-phase SPHの実装が概ね完了しました。

実装を通し、理論と実際を理解できました。
現状、2点の問題があるようです。

1つ目は seepage force (drag force) の計算式。

これまでの文献にも、DualSPHysicsでも同じ式が使われています。が、透水係数の大きな場合にしか実用的ではないようです。1e-4㎝/sオーダーだと、γ/kが大きくなりすぎて、少しの速度差が大きな拘束力を生みます。土塊がすべろうとしても、地下水がその動きに追随するまで拘束されるような挙動も見受けられました。

2つ目は密度変化。

土の間隙内の水が地表に浸出した場合、水相の密度は大きく変化しますが、この表現が困難。いえ、pressureの式で参照している基準密度を変更すれば容易に反映できるのですが、粒子が重なり合い、挙動がやや不安定でした。初期配置も面倒で、実用的ではないですね。水収支はともかく、土の動きだけを正しく反映したい、ということであれば計算の軽い single-layer two-phase formulation でよいと思います。

土相側でも層毎に密度を変えている場合、その境界で密度が訛るのは現状では仕方がないと考えています。が、この密度変化の問題は将来的に何とか解決されそうな気がします。たとえば、土層別にフラグを付けて、個別に計算するとか。それが実現象を正しく反映するようになれば、より実務に浸透していくのでしょう。

今年前半より取り組んできた SPH コードの修正と実装は、これで一旦終了です。
結果的には単相での計算はOK、2相は「定性的評価」まで、土木分野の実務で求められる諸問題には適用が難しい場合が多いと判断されます。が、もう少しだと思います。

将来に期待しましょう。

2024年8月13日火曜日

double-layer two-phase formulation その2

振り返ると、double-layer での実装に取り掛かってから1か月以上かかっています。

以下を参考に始めたのですが、試行錯誤の結果、最終的にはPersianSPHの文献と同じような形に落ち着きました。
A coupled fluid-solid SPH approach to modelling flow through deformable porous media - ScienceDirect


2次元、3次元共に土、水の相互作用までは動くのですが、地下水が地上に出た場合の密度変化をうまく表現できません。相互作用の検証も終わっていませんので、まだ時間がかかりそうです。

SPHでの2層2相表現は定性的にとどまると伺ったことがあります。それほど新しいものでもないので、国内でも実装されている方はいらっしゃいます。その先生がおっしゃるのですから、そうなのかもしれません。

まだ1か月。もう少し追いかけてみましょう。

2024年8月6日火曜日

システム修復

片道4時間、日帰りで、観測システムを修復してきました。

Linuxに不慣れだった後輩君が設定していたため細部の設定に漏れがあったこと、私の管理ミスでこれを見逃していたこと、夏場の高温でハードの一部に不具合が発生していたこと、最終的には停電がトリガーとなったこと等、いくつかの要因が重なり問題となりました。で、現場まで赴くことに。

原因を突き止めてしまえばその場で対応できる、比較的簡単なソフト側の対策で済んだのですが、熱だけはハード対策が必要そうでした。が、これは高価な製品を購入しないと対応できません。うーん。


システム修復後、ついでに DNS を設定。GUIからなぜか反映されていなかったので、コマンドで変更。

$ sudo vi /etc/systemd/resolved.conf

[Resolve]
DNS=8.8.8.8

$ sudo systemctl restart systemd-resolved.service


apt update は成功しますが、upgrade で失敗します。古いリストが残ったままなのでしょう。リストを全消ししてから再度 update & upgrade。upgrade には時間がかかりました。

$ sudo rm -r /var/lib/apt/lists/*


ひとまず周辺環境を含め目についた問題は全てつぶしておきました。今回は停電により、隠れていた問題も同時に顕在化しました。最初に見逃さないようにしないと、このように不要な対応が発生します。気を付けましょう。


2024年7月27日土曜日

MEMS Rain Gauge

IOTデバイスを用いて、振動から雨量を推定するプログラムを公開しました。https://github.com/T40O0/ADXL355_SPI_M5_SD_RainGauge/tree/main

約1年間観測して、改良を施し、安定して測定できるようになりました。

ポイントはデジタルフィルターを利用しないことでした。雨滴による衝撃を観測しているので、高周波まで扱えるほうが精度よく雨量を推定できることがわかりました。周波数領域は扱えませんが、ただただ時間軸で衝撃をカウントするという単純な方法が推定には良いようです。

転倒マスとは異なり、構造的に火山灰等で詰まることがなく、メンテナンスの楽な点が長所でしょう。火山灰を被っても、雨で洗い流されます。
安価な点も魅力です。個人で作るなら材料費2万円弱、従来の1/5~1/10の価格ですので、多点に設置できます。近い場所でも雨量は異なりますので、各地すべりブロックに1台、などという計画も可能です。
さらに、1滴の振動から観測できますから、時間方向の解像度が転倒マス型に比べ格段に向上します。

欠点は消費電力。電池のみで数か月持てば理想的なのですが、ソーラー+バッテリーか商用電源が必要です。モバイルバッテリーでも1週間しか持ちません。
また、振動のノイズに弱いことも短所でしょう。雨が降っていなくても、加速度を観測すると微小な雨量としてカウントします。設置場所毎にノイズを除去するための閾値の設定が必要です。

まだまだ改良の余地がありますが、ひとまずここで公開。良いアイデアが浮かべば追加していきましょう。

SPHの水圧計算

前回から約1か月かかりましたが、2相2層での水圧計算に目途が立ちました。

https://www.sciencedirect.com/science/article/pii/S002076831730286X

上記の78-82式を実装することで解けるようにはなったのですが、まだ土層を出た後の水の密度変化までは安定して解けません。最初の含水率のまま解くことはできるのですが、空中でそれは誤りですし、勝手に密度を上昇させることもできません。圧力が下がって自然に密度が上昇する、というように計算させるとすぐに発散します。そのあたり、この論文では明記されていません。土層の計算だけ追うのであれば、1層での計算が良いと思いますが、まだ何か見落としているようです。

ひとまず、次の drag force を安定して解けるか見てみましょう。

2024年7月3日水曜日

double-layer two-phase formulation

SPHにおいて、two-phase を single-layer で実装するか、double-layerにするか?

前者は計算が速い!
でも、水の流れは表現できません。(工夫すればできるのかもしれませんが、single ではなくなります。)

後者は計算が遅い!しかも圧倒的に。
でも、水の流れは表現できます。

考えた結果、もともと計算してみたかった「崩土から水の抜ける表現」が可能な後者を選択しました。既存のプログラムを容易に変更できそうな点も good です。

2相 を double-layer で解く文献はたくさんあります。が、なぜか式が統一されていません。以下が最も詳しく書かれていましたので、これを参考に実装します。
A coupled fluid-solid SPH approach to modelling flow through deformable porous media - ScienceDirect

式を整理しながら不飽和の表現を single-layer の文献からいただきました。これで、飽和・不飽和浸透と力の関係を解けるはずです。

コツコツ、組んでみましょう。

**********************************
20240707追記
地下水の上下で飽和度(水の単体)を変更すると、境界で爆発してしまいます。計算過程を追えば当然なのですが、土層境界よりも挙動が激しいようでした。飽和度を低くするには文献のように水粒子をまばらに配置するしかなさそうです。
似た現象ですが、地表流がある直下で飽和度が過剰に高くなります。透水係数が実際よりも高く計算され、水が入りやすくなります。なかなか難しい。
SPHでは密度変化がないか小さい現象(弱圧縮)を扱うのが前提なので、単体の異なる土層とか飽和‐不飽和‐地表流などを扱うのが難しいのでしょうね。

2024年6月26日水曜日

single-layer two-phase formulation

SPH で流体とソリッドをカップリングさせる方法は比較的古くから提案されています。

以前、PersianSPHを触っていた際に文献でおおむね理解していたのですが、実装目的であらためて詳細を確認しました。

基本は土質力学レベルです。それを丁寧に組み合わせると two-phase formulation の出来上がりです。SPH の場合はこれを2レイヤーで計算します。ここで倍の計算量。さらに相互作用を反映するため、土粒子の位置での水粒子の動き、水粒子位置での土粒子の動きを計算します。単純に4倍ではありませんが、計算量は飛躍的に増えます。

これを解決した、という文献がコチラ⇓
Two-phase fully-coupled smoothed particle hydrodynamics (SPH) model for unsaturated soils and its application to rainfall-induced slope collapse - ScienceDirect
A general smoothed particle hydrodynamics (SPH) formulation for coupled liquid flow and solid deformation in porous media - ScienceDirect
1レイヤーで2相の計算が可能に工夫されています。不飽和浸透も取り扱えるので一通りの計算ができそうです。が、境界条件が面倒。いえ、実装されたものを使用するのであればDtransu を利用するのと変わらないでしょうが、汎用性を持たせて実装しようとすると気が遠くなります。

著者には DualSPHysics の開発者が含まれています。GPU対応で開発されているのでしょう。実装されるまで待ちましょうか。

2024年6月23日日曜日

火山地帯の観測機器

先週は火山につけた観測機器のメンテナンスでした。

いくつか新しい機器を設置し、古い機器を入れ替えて帰社したのですが、あらゆるモノが錆びていました。火山ガスと雨でしょうね。防食、大事です。

PCは火山灰を吸い込み異音が発生。分解、清掃しても、治るかどうか。寿命を短くしたのは間違いなさそうです。

こういった環境では、安価な機器を使い捨て&ばらまきにするか、お金をかけて防食し長期運用&集中させるかのどちらかでしょう。費用が1~2桁違いますので、重要度に応じて使い分けるのが良いのでしょう。


2024年6月16日日曜日

SPHのベンチマーク

組み終わったSPHで静水圧のチェックは終わりました。残る弾塑性のチェックをどうしようか、と調べてみましたが、適当なベンチマークはなさそうでした。

流体のベンチマークは多く見かけます。が、弾塑性のものは少ないようです。FEMと同様に、室内3軸圧縮試験であったり、支持力公式との比較をいくつかの文献で見かけますが、ベンチマークとして共有されるには至っていないようです。

3軸は少し面倒ですので、支持力を試してみました。
が、文献に多くあるように、解析値のほうがやや大きく出る結果になります。うーん。ま、FEMと同程度に大きな値になっているので、あっているともいえるのですが。

新しい機能も組み込みつつ、もう少し探してみましょう。

2024年6月8日土曜日

EOS

equation of state

WCSPH で圧力を計算する際に用いられる式です。Monahan らが提案した Tait 型だそうですが、これ、簡単に見えて扱いが難しい。

係数を決めるのに最大流速を考慮しないといけないのですが、流速は計算しないとわかりません。ま、感覚で大体この程度、圧縮率は〇%程度をネラッテ音速ハコノクライ、みたいな決め方をします。DualSPHysics に使われている multi-phase (liquid-sediment)の文献では、音速はドメイン内の最大流速の10倍以上が利用されています。これは、マッハ数が0.1なので、概ね0.5%までの圧縮を許容するということでしょう。
https://research.manchester.ac.uk/en/publications/modelling-multi-phase-liquid-sediment-scour-and-resuspension-indu

では、個体が弾性から塑性になり、さらに非ニュートン流体として流動するような現象を追う場合はどうするか?

弾性、塑性では変形係数とポアソン比から音速を求められますので、それぞれ別の値を利用するという方法が考えられます。流体は水のように扱うことも可能です。
が、それぞれ別の音速を使うと計算が破綻しやすくなります。試算では塑性の計算で圧力が低くなり、塑性領域が広くなって流動化しやすくなるような結果を得ました。あと、2Dでは動くが3Dだと進まなくなるとか(2Dの計算って、問題やバグが内在しても露見しにくいことを、最近強く感じています)。

結論としては、弾性として求めた音速をそのまま利用するか、徐々に切り替えるか、なのでしょう。文献には出てこないのでノウハウの部類なのでしょうが、皆さんどうされているのか知りたいところです。


2024年5月30日木曜日

弾塑性SPH その2

新たに頂いたコードを参考に、以下の組み合わせで計算ができるように改造しました。

・弾性
・弾塑性
・弾性‐流体
・弾塑性‐流体
・流体

これらのほとんどは日本語の図書でも紹介されています。個体と流体を同じように実装できるのはSPHならではでしょう。

弾塑性(一部流体)の挙動はコチラ。端から土圧により崩壊し始めます。

少しc・φをあげると、崩れにくくなります。その場合、まず角が楔のように崩れ、次に長辺の2/3が崩れ、最後に残りが崩れるような時間差が見られました。ちなみに、前回のアップしているイメージは「弾性‐流体」の計算結果で、全体がペシャっと潰れるような挙動を示しています。これだけ自由な表現ができると、多様な自然現象を再現することが可能でしょう(予測は相変わらず困難ですが)。

まだ、いくつか疑問点が残っています。資料を探して早く完成させましょう。

2024年5月26日日曜日

壁粒子の圧力

SPHで壁際の粒子が壁に張り付く現象を、後回しにしていました。
ある程度動くようになったので取り組んでみることに。

まずは mDBC で様子見。
効果はありますが、残念ながら完全ではありません。いくつかの粒子は壁に残り、計算時間も余計にかかりました。

負圧が効いているのは想像がつくのですが、他の修正方法を思い付きません。
そこで googling。
すると、「まさにコレ!」といった回答を見つけました。条件も全く同じです。
Trying my chances here; anyone knows why particles stick to the wall in my SPH simulation? - Specific Domains / Modelling & Simulations - Julia Programming Language (julialang.org)

Here the simplest fix is simply to put the boundary particle density to reference density, if it ever gets below that reference value.
In this way a wall can never produce suction, but only “push”.

当たり前すぎて一般的には解説されていないのでしょうね。初心者には思いつかないノウハウです。
密度でなく圧力を変更。数行の追加で計算時間もかからず、簡単に解決できました。


2024年5月18日土曜日

弾塑性 SPH

2020年の秋から始めた弾塑性SPHコードの修正が完了しました。

基本コードを頂いてから寝かせていたため実質1年程度ですが、SPHを知ってからは14年目です。長い。
何とかSPHの離散化や各種補正方法、境界条件について基本的なしくみを理解したと言える証になるでしょう。いや、よかった。 

いくつかのハマり処があったのですが、式の実装時の静水圧Pと静水圧応力σの符号の整理が効いたと思われます。多くの文献を引用する際に、これらの書き方の流儀を統一してから実装に落とし込む必要があったのでしょう。

粒子が凝集する現象を hourglass,  zero energy mode だと勘違いしたり、他にも似たような勘違いで多くの寄り道をしましたが、振り返るとそれらも良い勉強になったと思います。
https://phreeqc.blogspot.com/2022/05/blog-post.html

早くVerification を済ませたいですね。ベンチマークは何が良いのでしょうか?
それが済んだら高速化。倍精度を多く使っていますので、GPUではなくスパコン利用(MPI並列)のほうが現実的でしょうね。

まだまだ先は長い。


2024年5月6日月曜日

相互相関関数の計算

西田究のホームページ (u-tokyo.ac.jp)
解説 (u-tokyo.ac.jp)

地震動に関する様々なシミュレーション結果をお手軽に見ることができる、素晴らしいサイトです。特に、地震波干渉法で相互相関関数が出来上がっていく過程を見たことがなかったので、非常に興味深く見ていました。

  初期
数十分後

相互相関関数を求める際は、時間領域(直接法)か周波数領域(FFT)を利用します。配列の要素数が多ければ周波数領域での計算が早いでしょう。Python の scipy.signal.correlate では手法を指定できますし、指定しなくとも早い方を自動で選択してくれます(Ver.1.13.0)。 numpy.correlate は直接法のみです(Ver.1.26)。

以下、計算例です。設定次第なのですが、以下だと答えは①≠②=③
①直接法

sig1 = np.std(data1) sig2 = np.std(data2) data1 = data1-data1.mean() data2 = data2-data2.mean() cross_correlation = np.correlate(data1, data2, mode='same'                                 ) /sig1/sig2/N

②FFT1

f1 = np.fft.fft(data1) f2 = np.fft.fft(data2) cross_spectrum = f1 * np.conj(f2) cross_correlation = np.real(np.fft.ifft(cross_spectrum)                             ) /sig1/sig2/N

③FFT2

f1 = np.fft.rfft(data1) f2 = np.fft.rfft(data2) cross_spectrum2 = f1 * np.conj(f2) cross_correlation = np.fft.irfft(cross_spectrum2                                 ) /sig1/sig2/N

0をパディング + シフトで①=②となります。

①直接法

padded_data1 = np.pad(data1, (0, N), 'constant')
padded_data2 = np.pad(data2, (0, N), 'constant')
cross_correlation = np.correlate(padded_data1,
                 padded_data2,
                 mode='same' ) /sig1/sig2/N

②FFT1

f1 = np.fft.fft(padded_data1) f2 = np.fft.fft(padded_data2) cross_spectrum = f1 * np.conj(f2)
cross_correlation = np.real(np.fft.ifft(cross_spectrum))
cross_correlation = np.fft.fftshift(cross_correlation ) /sig1/sig2/N