2021年1月31日日曜日

確率密度関数のフィッテイング

確率密度関数のデータへのフィッテイングは SciPy の各分布に実装されています。

最尤推定法のようですが、一度に多数の関数を比較するようにはなっていません。
何かないか?と探してみると、すぐにヒット。皆さん、考えることは同じようです。
https://stackoverflow.com/questions/6620471/fitting-empirical-distribution-to-theoretical-ones-with-scipy-python
distfit は動きませんでした。が、tmthydvnprt & Melquíades Ochoa のコードが使いやすいですね。最小二乗法を用いた総当たり検索です。

手元のデータに当てはめると、数時間かかっても計算し続けていましたので、関数を減らして再実施。無事に最良の関数を選んでくれました。ありがたい。

以前、確率降水量の計算でヒストグラムを作った際に、残差の小さい関数を選択すれば良いのではないかと考えていました。
https://phreeqc.blogspot.com/2019/06/blog-post_20.html

ひとまず、どのようなデータでも計算できるようになりました。

あとは最適化手法でしょうか。

2021年1月26日火曜日

デルタ関数のフーリエ変換

 ω=2πf=2π/T・・・角振動数ωは振動数fの2π倍
\begin{align*}cos{\omega t},\ \sin{\omega t}・・・ω=1(=cost,\ \sin{t})\end{align*}で周期T=2πの周期運動
\begin{align*}x=x_0cos{\left(\omega t+\alpha\right)}\end{align*}調和運動:()の中身が t の一次式

 

自由振動の一般解
\begin{align*}m\ddot{x}+kx=0\end{align*}
【表現1】解の形を\begin{align*}x=Xe^{\lambda t}, λ=iω\end{align*}とおいて展開
\begin{align*}x=X_1e^{i\omega t}+X_2e^{-i\omega t}\end{align*}2つの任意定数を含んだ解=一般解
【表現2】オイラーの公式を用いて書き換え
\begin{align*}x=\left(X_1+X_2\right)cos{\omega t}+i\left(X_1-X_2\right)sin{\omega t}\end{align*}
\begin{align*}=acos{\omega t}+bsin{\omega t}\end{align*}
【表現3】三角関数の合成の公式を用いて書き換え
\begin{align*}x=a_0cos{\left(\omega t+\alpha\right)}\end{align*}

 

周期関数のフーリエ級数
【表現2】の利用
\begin{align*}f\left(t\right)=\frac{1}{2}a_0+a_1cos{\omega t}+a_2cos{2\omega t+\ldots}\end{align*}\begin{align*}+b1sinωt+b2sin2ωt+⋯\end{align*}
オイラーの公式を利用して表現を変えると
\begin{align*}cos{\theta}=\frac{e^{i\theta}+e^{-i\theta}}{2},\ sin{\theta}=-i\frac{e^{i\theta}-e^{-i\theta}}{2}\end{align*}
\begin{align*}f\left(t\right)=\frac{1}{2}a_0+\sum_{r=1}^{\infty}\left(a_r\frac{e^{ir\omega t}+e^{-ir\omega t}}{2}{-jb}_r\frac{e^{ir\omega t}-e^{-ir\omega t}}{2}\right)\end{align*}
複素フーリエ級数
\begin{align*}f\left(t\right)=\sum_{r=‐∞}^{∞}{F_re^{ir\omega t}}\end{align*}
複素フーリエ係数
\begin{align*}F_r=\frac{1}{T}\int_{-T/2}^{T/2}{f\left(t\right)e^{-ir\omega t}dt}\end{align*}

非周期関数のフーリエ積分
\begin{align*}f\left(t\right)=\int_{-\infty}^{\infty}{F\left(\omega\right)e^{i\omega t}d\omega}\end{align*}
フーリエ変換・・・フーリエ級数のフーリエ係数に相当
\begin{align*}F\left(\omega\right)=\frac{1}{2\pi}\int_{-\infty}^{\infty}{f\left(t\right)e^{-i\omega t}dt}\end{align*}

 

デルタ関数のフーリエ変換
t=0の時、e^0=1より
\begin{align*}F\left(\omega\right)=\frac{1}{2\pi}\int_{-\infty}^{\infty}{\delta\left(t\right)e^{-i\omega t}dt}=\frac{1}{2\pi}\end{align*}




2021年1月23日土曜日

Windows コマンド その2

ちょうど1年前、windows コマンドを利用したファイル操作について書き残していました。
https://phreeqc.blogspot.com/2020/01/windows.html

今回は対象のフォルダ内のファイルをリスト化しておき、fort 文でコピーや圧縮を行う書き方です。前回と異なるのは、操作前後でファイル名を変えている点。

リストを参照しファイルをコピー

rem @echo off
set listfile=aaa.txt
set input_path="\\bbb\ccc\"
set output_path="D:\eee\fff\"
FOR /F %%I IN (%listfile%) DO (
    copy %input_path%pre_%%I.csv %output_path%post_%%I.csv
)
echo complete!
pause

 

リストを参照し、ファイル毎に個別圧縮(7-zip 利用)

rem @echo off
set ZIP_PATH="C:\Program Files\7-Zip\7z.exe"
set listfile=aaa.txt
set input_path="\\bbb\ccc"
set output_path="D:\eee\fff"

FOR /F %%I IN (%listfile%) DO (
    %ZIP_PATH% a -tzip %output_path%\%%I.zip %input_PATH%\%%I.csv -mx=4
)
echo complete!
pause

 

2021年1月18日月曜日

地形の曲率

西田ほか「数値地形モデルに基づく地震時山腹崩壊斜面の地形解析」1997

いくつかの国総研資料にて、この文献による「平均曲率」が利用されています。文献でいう4方向ラプラシアンでは斜度により数値が変化するため、変化しない平均曲率が良いとのことでした。

式自体は汎用ソフトに実装されている式と似ています。が、5mDEMの利用で計算メッシュの距離を10m~50m とすると、カーネルが 3x3 では収まりません。これ、GISソフトで対応しているものがないのでは?

まず、ArcGIS。この文献で否定されている式のみの実装なのでダメですね。
SAGA には複数の曲率計算方法が実装されていますが、該当するものはないようです。
https://sourceforge.net/p/saga-gis/code/ci/master/tree/saga-gis/src/tools/terrain_analysis/ta_morphometry/Morphometry.cpp
GISではないものの、Surfer にもいくつか曲率の計算が実装されています。ウェーブレットで9x9のユーザー定義フィルターを使えましたので、もしかするとできるのかもしれません。が、手元になかったので確認できませんでした。
http://surferhelp.goldensoftware.com/gridops/profile_curvature.htm

数字の並びにフィルタをかけていくのは画像処理と同じ作業です。OpenCV であれば、どのような間隔・カーネルでも対応可能です(しかも計算が速い)。いつもの Python だと、手順は以下の通りです。

1. GDAL で GeoTIFF 読み込み
2. NumPy でカーネル作成
3. OpenCV で 2D フィルタをかけて ∂h/∂x ほかを計算
4. H を計算
5. GDAL で GeoTIFF 書き出し

簡単、と思って手を動かし始めたのですが、詰まりました。
いえ、すべて正しく組めて速く走るのですが、できた GeoTIFF を ArcGIS や QGIS に読み込ませると最大・最小の表示が狂います(値自体は正しく読み込まれています)。
一晩考えたものの、原因を特定できませんでした。

今朝、同じ問題で悩まれていた方を発見。これは知りませんでした。
https://gis.stackexchange.com/questions/29064/how-to-get-gdal-to-create-statistics-for-gtiff-in-python

6. 統計量を XNL で出力

最後にこの手順を加えると、あっさり解決。
今まで Geopandas では対応できないラスターの読み書きを避けていましたので、付けが回ってきたのでしょう。

やはり、手を動かすと、何かしらの発見があるものです。

**********************
20220731追記

ArcGIS Pro の エクステンションに追加されたようです。
https://pro.arcgis.com/ja/pro-app/2.8/tool-reference/spatial-analyst/surface-parameters.htm

平均曲率は Minár(2020)とされていますが、その文献の参照先は Gauss(1828)になっていました。

**********************
20220802追記
ラスターでもフィルタの距離を指定できました。
計算結果は±を逆転しても一致しません。

「ローカルサーフェスタイプ」を指定して、サーフェスを作成してから、正規傾斜ライン?=コンター直交方向と接線方向を求めXY軸に定めているようです。谷底ではどうなるのでしょうか。

 

2021年1月11日月曜日

浸透力

M Gholami Korzani et al.
Smoothed Particle Hydrodynamics for investigating hydraulic and mechanical behaviour of an embankment under action of flooding and overburden loads, Computers and Geotechnics 94, 31-45

昨日の文献内で頻繁に参照されていた文献です。
同様の内容でしたが、あっ!と思った部分があり書き残しておきます。

単純な浸透模型で浮力と浸透力について説明されています。その流速はダルシーの法則から求められるのですが、これが重要。式27にあるように、土粒子が移動する場合、その速度(ベクトル)も考慮しないと正しい浸透力が得られません。SPH で解く場合だと当然のように着想するのですが、これまでFEMで解いていた際には考えもしない要素でした。頭が凝り固まっています。

もう一つ、知らなかった式。

However, the Forchheimer equation, i = av + b|v|v, considers the microscopic inertial effect for transient and turbulent flow conditions by adding a quadratic velocity term to Darcy equation.
調べてみると、こちらは日本でも知られているようです。a,b は半経験式のようで、この文献では昨日の Kozeny-Carman Equation をベースとした式が採用されていました。

これらは昨日の内容含め非常にシンプルな説明ですので、実装は難しくないでしょう。いえ、触った感じだと2次元なら著者の c++ コードを使わせてもらえば実務でも対応できるレベルだと思われます。3次元になると実務では厳しい。プレが必須になるでしょうね。そうなると自作になるのかな?

ま、手元にある問題を解決してからですね。

 

2021年1月10日日曜日

Kozeny-Carman Equation

M Gholami Korzani et al.
SPH approach for simulating hydro-mechanical processes with large deformations and variable permeabilities, Acta Geotechnica 13 (2), 303-316
ソース:https://github.com/mghkorzani/persiansph
著者HP:https://korzani.wixsite.com/persiansph/persiansph

シンプルですが、パワフルな内容です。

SPH を利用した土‐水連成の解法が示されています。キーワードは "Kozeny-Carman equation"。

土の SPH を触っていた際、空隙に水を通す方法を考えてみたのですが思いつきませんでした。が、この関係式を使えば簡単。統一的に実装できます。
土が変形→空隙が変化→空隙率に応じ透水係数が変化→水粒子の速度が変化、といった内容です。粒子は離散点であり、物理的な粒子とは異なる点がミソ。

私はこの関係式を知らなかったのですが、公害防止管理者の試験には毎年出題されるくらい有名な式のようです。最初に式を見た際、有効径で透水係数を変化させるのが、いささか乱暴なような気がしました。が、考えてみると層に固有の透水係数を1つ定めて計算する現状の方が乱暴に思えてきます。なかなか良い手法かもしれません。
これに従えば、地表流から浸透流まで、非常にシンプルな支配方程式で記述できます(不飽和は別途追加の必要がありますが)。土の計算を飛ばせば水だけの計算ができるので、簡単なモデルなら3次元も可能かもしれません。
非常に計算負荷が高くなりそうなので実務への適用は2次元どまりかなとも思われますが、将来的には(というかプロなら現段階でも)3次元で問題ないでしょう。GPGPU、何とかしたいですね。

***************************
20200111追記
3年前の状況です。
https://phreeqc.blogspot.com/2017/03/dualsphysics-multi-phase.html

2021年1月5日火曜日

地下水 と MCMC その2

Using Markov Chain Monte Carlo to quantify parameter uncertainty and its effect on predictions of a groundwater flow model - ScienceDirect 

数千回もの浸透流計算が必要となるパラメータ推定手法は、実務向きではありません。https://phreeqc.blogspot.com/2020/12/mcmc.html
この文献では、既存の知見により浸透流計算を省略しています。具体的には、尤度関数内の二乗誤差、誤差分散計算の代替えとして、涵養量と透水係数の比 R/K との相関性を利用します。

overcome this obstacle, 300 realizations of Hassan et al.’s(2002) model were studied and it was found that the terms needed in the likelihood function, namely the sum of squared errors and the error variance are highly correlated to the recharge–conductivity ratio, R/K.

この文献の前に、300 realizations の計算があるのでしょう。それに気づいて、途中で細かく読むのをやめてしまいました。R/K などの相関関係は他の現場で使えたとしても、やはり事前の計算による確認が必要なのでしょうね。

おっ!と思ったり、うーん、と思ったりした文献でした。