2020年12月31日木曜日

やり残し事項 2020

GEE コミュニティのチュートリアルを再度実施してみようと手を動かしました。
https://phreeqc.blogspot.com/2020/12/gee-sentinel-1.html

1か月でPart3まで増えていました。仮説検定を用いた判定がメイン。今回は任意の場所でも1~3まで走りました。が、結果はイマイチ。
これまで、InSAR、Change Detection ともに、良い結果を得られたことがありません。イマイチなんですよね。やはり、有償ソフトの出番でしょうか。この辺りは来年の短期目標にしましょう。

職種が変わったこともあり、地質関連の作業がなくなりました。帯磁率計などはやり残したままで終わりました。

あらためてやり残し事項を整理しておきましょう。

優先度高:粒子法コード のGPGPU対応
優先度高:SAR・画像の Change Detection
優先度低:PSInSAR・・・700万で購入可
優先度低:DtransuのGPU対応・・・100万で購入可
優先度低:地表流+地下水・・・GSFLOW、HydroGeoSphere
優先度中:OpenMI・・・PHREEQCとDtransuの連成等

中期目標4年目は惨敗。新たな目標を考えましょう。

2020年12月30日水曜日

DataFrames.jl vs. pandas

 Python は私のようなノンプロが手を出しやすい簡素な言語です。

一方、ライブラリへの依存が高いため、それらの互換性や version に応じて複数の仮想環境の構築を迫られることがあります。将来的には、作成していたコードがライブラリのVer.UP や開発中止で修正を余儀なくされることもあるでしょう(前者は既に経験済みです)。
個人的には 精度や速度が必須となる科学技術計算 を Fortran 、データ整理や可視化を Python、といったような使い分けをしています。現状維持で良いともいえるのですが、そろそろ代替え言語を試しておくべきと思い始めました(リスク管理です)。

先日、ある方から Julia を勧められました。Python と同じくらい容易な書き方で、計算が速く、コンパイルができるとのこと。
それは良い。

ということで、試してみました。
少し古いのですが、LTS の v1.0.5 + win10 (i5のラップトップ)です。インストールはプレコンパイル版を入れるだけ。そのままでは Jupyter での利用は不可でしたが、何度か試行し、動くようになりました。
今回の比較は csv 読み込みからの DataFrame 操作。頻繁に使う機能です(Juria と Python の比較ではなく、DataFrames.jl と pandas の比較になりました)。
結果は以下の通り。

csv読み込み df=DataFrame(CSV.File(open(read, 1170万行 84.8s
df6=DataFrame(CSV.File(open(read, 6万行 0.50s
df=pd.read_csv 1170万行 51.4s
df6=pd.read_csv 6万行 0.14s
カラム名変更 rename!(df,   0.21s
df.rename(columns= ,inplace=True)   0.02s
inner_join dfj=innerjoin(df, df6 290万行 15.5s
dfj=pd.merge(df, df6 290万行 5.8s
型変換 dfj._=Date.(dfj._,"yyyy/mm/dd")   86.7s
dfj['__']=pd.to_datetime(dfj['__'])   1.4s
期間抽出 df0=dfj[dfj[:,"__"] .< _, :] 190万行 4.6s
df0=dfj[dfj['_’] < _] 190万行 3.6s
集約 df0_st=stack(df0, 380万行 5.6s
df0_pv=pd.pivot_table(df0,values= 190万行 1.4s

pandas の圧勝でしたね。強い。

Julia でもパッケージを使うので Python と同じような問題は潜在的に有しているのでしょう。
イタレーションでは Julia の方が速くなるそうですので、扱う作業内容により言語を選択する形になるでしょうか?(今回の結果を見る限り、当面の出番はないかもしれませんが)
https://julialang.org/benchmarks/ 

引き続きその動向に注目しましょう。

2020年12月29日火曜日

常微分方程式

振動などの物理現象を解く際には、常微分方程式の解法が利用されています。

が、既に学んだことすら記憶にありません。
冬休み中に再読。以下にメモ。学ぶべき時に励んでおけばよかった。

*******************************

未知関数 y(x) と x および導関数 y'(x),y"(x),・・・yn(x) からなる方程式を n階常微分方程式という。
F(x,y(x),y'(x),y"(x),・・・yn(x))=0
このとき、y(x) を常微分方程式の解という。(方程式を解く=式を求める)

常微分方程式の分類
・線形・非線形
・定数係数・変数係数
・同次・非同次

登坂宣好「微分方程式の解法と応用」より常微分方程式の例
線形‐変数係数‐非同次形
\begin{align*}\frac{dy}{dx}\left(x\right)+x^2y\left(x\right)=e^x\end{align*}
線形‐定数係数‐非同次形
\begin{align*}\frac{dy}{dx}\left(x\right)+3y\left(x\right)=e^x\end{align*}
線形‐変数係数‐同次形
\begin{align*}\frac{dy}{dx}\left(x\right)+x^2y\left(x\right)=0\end{align*}
線形‐定数係数‐同次形
\begin{align*}\frac{dy}{dx}\left(x\right)+3y\left(x\right)=0\end{align*}
非線形‐非同次形
\begin{align*}y\left(x\right)\frac{dy}{dx}\left(x\right)+3y\left(x\right)=x^3\end{align*}

解法(最近、使用した解法のみ記載)
線形-1階-定数係数-非同次形
\begin{align*}\displaystyle\frac{dy}{dx}\left(x\right)=f\left(x\right)\end{align*}
直接積分形:両辺をxで積分(両辺にdxを乗ずる形にして積分)
\begin{align*}f\left(x\right)=x : y\left(x\right)=\int xdx=\frac{1}{2}x^2+c\end{align*}\begin{align*}f\left(x\right)=sinx:y\left(x\right)=\int sinxdx=-cosx+c\end{align*}
線形-1階-定数係数-同次形
\begin{align*}\displaystyle\frac{dy}{dx}\left(x\right)=ay\left(x\right)\end{align*}
指数関数解:
\begin{align*}\frac{d}{dx}e^{\lambda x}=\lambda e^{\lambda x}=ae^{\lambda x}\end{align*}\begin{align*}\left(\lambda-a\right)e^{\lambda x}=0\rightarrow\lambda=a\end{align*}
解として\begin{align*}e^{ax}\end{align*}を有する。定数倍した\begin{align*}Ce^{ax}\end{align*}も解となる。(一般解)

線形-n階-定数係数-同次形
\begin{align*}\frac{d^ny}{dx^n}\left(x\right)+a_1\frac{d^{n-1}y}{dx^{n-1}}\left(x\right)+\ldots+a_ny\left(x\right)=0\end{align*}左辺に\begin{align*}y=C_1y_1+C_2y_2+\ldots+C_ny_n\end{align*}を代入すると\begin{align*}\frac{d^n}{dx^n}\left(C_1y_1+C_2y_2+\ldots\right)+a_1\frac{d^{n-1}}{dx^{n-1}}\left(C_1y_1+C_2y_2+\ldots\right)+a_n\left(C_1y_1+C_2y_2+\ldots\right)\end{align*}\begin{align*}{=C}_1\left({\frac{d^n}{dx^n}y}_1+a_1\frac{d^{n-1}}{dx^{n-1}}y_1+\ldots+a_ny_1\right)+C_2\left({\frac{d^n}{dx^n}y}_2+a_1\frac{d^{n-1}}{dx^{n-1}}y_2+\ldots+a_ny_2\right)+\ldots\end{align*}\begin{align*}{=C}_1・0+C_2・0+\ldots C_n・0=0\end{align*}よって\begin{align*}y=C_1y_1+C_2y_2+\ldots+C_ny_n\end{align*}は一般解となる。
解の形が\begin{align*}y=Ce^{\lambda x}\end{align*}のとき、
\begin{align*}C\ \left(\lambda^n+a_1\lambda^{n-1}+\ldots+a_n\right)e^{\lambda x}=0\end{align*}
C=0のとき、y=0は自明な解。\begin{align*}\lambda^n+a_1\lambda^{n-1}+a_2\lambda^{n-2}+\ldots+a_n=0 :特性方程式\end{align*}のとき、非自明な解\begin{align*}y=C_1e^{α_1x}+C_2e^{α_2x}+⋯+C_ne^{α_nx}\end{align*}が存在。

非線形(特別な形式)
変数分離型解法:両辺をg(y)で割ってから積分(g(y)≠0)
\begin{align*}\displaystyle\frac{dy}{dx}\left(x\right)=f\left(x\right)g(y\left(x\right))\end{align*}\begin{align*}\frac{1}{g\left(y\left(x\right)\right)}dy=f\left(x\right)dx\end{align*}\begin{align*}\int\frac{1}{g\left(y\left(x\right)\right)}dy=\int f\left(x\right)dx\end{align*}


2020年12月27日日曜日

Mann-Whitney's U test

Mann-Whitney's U test を実施。

scipy.stats.mannwhitneyu を利用していたのですが、ランキングのつけ方が書かれていません。Σ の計算に n(n+1)/2 が含まれていますので、1位、2位、2位、4位のような同順位のつけ方はまずい。そこは正しく計算してくれているのだろうと思いながら、初めて使うので確認。

ソースでは rankdata が使用されています。

ranked = rankdata(np.concatenate((x, y)))

https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rankdata.html?highlight=rank#scipy.stats.rankdata

>>> from scipy.stats import rankdata
>>> rankdata([0, 2, 3, 2])
array([ 1. ,  2.5,  4. ,  2.5])
>>> rankdata([0, 2, 3, 2], method='min')
array([ 1,  2,  4,  2])
>>> rankdata([0, 2, 3, 2], method='max')
array([ 1,  3,  4,  3])
>>> rankdata([0, 2, 3, 2], method='dense')
array([ 1,  2,  3,  2])
>>> rankdata([0, 2, 3, 2], method='ordinal')
array([ 1,  2,  4,  3]) 

はい、一番上の平均 (default is ‘average’) が利用されていました。OKです。

比較する母集団の平均、標準偏差、サンプル数の違いはどのように影響するのでしょうか?簡単な正規分布のプラス側を取り出してチェックしてみました(プロットは正規化済み)。シード値は設けず、何度かテストを繰り返して両側・有意水準5%にかかわる傾向を見てみました。

x1=np.random.normal(1, 10, 1000)
x2=np.random.normal(1, 10, 1000)
x1=x1[x1 > 0]
x2=x2[x2 > 0]   

当然ですが、棄却されない分布です(稀に棄却されるのが気になりますが)。


2群が同じ平均・標準偏差で、一方のサンプルサイズを1000から100000に増やしたケース:検定結果に大きく影響せず(棄却され難い)。



2群が同じ平均・サンプルサイズで、一方の標準偏差を10から20に増やしたケース:検定結果に影響(棄却され易い)。



2群が同じ標準偏差・サンプルサイズで、一方の平均を1から2に増やしたケース:検定結果に大きく影響せず(棄却され難い)。

 

2群が同じ標準偏差・サンプルサイズで、一方の平均を1から3に増やしたケース:検定結果に影響(棄却され易い)。



2群が同じ標準偏差で、一方の平均を1から3、サンプルサイズを1000から100000に増やしたケース:検定結果に影響(棄却され易い)。


ランダムサンプリングをしているため、棄却・棄却されないという結果は固定になりません。それを繰り返し実行しただけですが、比較的容易に傾向が見えてきました。
平均2、3は見た目が微妙ですね。やや大きいかな?というような違いでも「有意差がないとは言えない」という傾向を出せるようです。人によってはやや大きい、あるいは同じように見える、などと意見が分かれるような分布ですが、ある考え方に従い再現性のある答えを出してくれるのはありがたい。

最近では、土研交流研究員の方が景観に配慮した護岸に関する類似性の評価に、この検定を利用されていたように思います。アイデア次第でいろいろな分野に適用できそうです。

好きではなかった統計的手法ですが、その合理性をもっと利用すべきでしょうとも考えるようになってきました。


2020年12月22日火曜日

shpfile in ZIP

先日、防災科研さんの地すべり地形分布図をプロットしようとして、ちょっとした問題に引っかかりました。

1.ZIP ファイルに5つの SHP ファイルが納められている。
2.ZIP ファイルは150個。
3.ZIP ファイルから、目的の SHP のみを解凍せずに取り出したい。
4.取り出した SHP を結合し、1つの SHP として保存したい。

QGIS で試みるも、複数の SHP からその一部を機械的に取り出せません。
ArcGIS Pro 2.5 では ZIP を扱えません。技があるのかとサポートに問い合わせましたが、やはり標準では対応しておらず、開発&有償になるとのこと(余談ですが、サポートさんのレスが以前に比べずいぶん遅くなりました)。

それほど大仰なモノでもないでしょう。
ということで、簡素にできたのがコチラ。

import os
import glob
import pandas as pd
import geopandas as gpd

path = './input/*.zip' #folder 'input'
zfiles = glob.glob(path) #zip files in the folder
gdf_300=gpd.GeoDataFrame()

for zfile in zfiles:
    fname=os.path.splitext(os.path.basename(zfile))[0]
    gdf_temp=gpd.read_file('zip://'+str(zfile)+'!ls300_'+fname+'.shp')
    gdf_300=pd.concat([gdf_300, gdf_temp])

gdf_300=gdf_300.reset_index(drop=True)
gdf_300.to_file('./output/300.shp', encoding='utf-8')

GIS のカジュアルユーザーだからでしょうか、Python との併用が合理的に感じます。 


2020年12月17日木曜日

5次メッシュを3次メッシュに集約

AVS30(5次メッシュ)を 3次メッシュ内で平均化したいなあと思案。

全国だと EXCEL では行数が足りないので、正攻法の ArcGISで処理。

時間がかかるので Python でも追っかけ処理。
shpを読んだ後の主な処理はこれだけ。

#AVS30=0(水域)を削除
gdf=gdf[gdf['5th_CODE'] != 0]
#3次メッシュコード作成
gdf['3rd_CODE']=gdf['5th_CODE'].astype(np.int64)//100
#3次メッシュにおける平均を求める
df=gdf.groupby('3rd_CODE').mean()

空間演算が必要ないので、Python の方が早く終わりました。
便利。

2020年12月14日月曜日

scipy.optimize.curve_fit

任意関数の最小二乗法は EXCEL でも Pythonでも可能です。

精度、速度、記載の単純さで scipy.optimize.curve_fit の圧勝です。初期値によっては一発で求まらないこともありますが、 EXCELソルバーに比べると安定かつ高速で解を求めてくれます。
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html

変分推論は難しいですね。収束させるのにあたりのつけ方が必要になるようです。で、結局は最小二乗法に戻ってしまったり。慣れが必要なのでしょう。

 

2020年12月13日日曜日

DataFrame のアダマール積

Pandas の DataFrame でよく使う便利な乗算。
間違えやすいので、( ..)φメモメモ

**************************************

2020年12月7日月曜日

地下水 と MCMC

ようやくMCMC を概ね理解し動かせるようになったので、地下水モデル10章の不確実性に関する部分へ少しだけ立ち戻ることに。
https://phreeqc.blogspot.com/2019/07/4.html

MCMC関連の文献を探してみると、簡単な計算をされている事例がありました。第2著者以降は日本の方ですね。
Julien Boulange, Hirozumi Watanabe, Shinpei Akai (2017) A Markov Chain Monte Carlo technique for parameter estimation and inference in pesticide fate and transport modeling

  • 気候、水収支、土壌、農薬特性などの40以上のパラメータから、農薬濃度の予測精度に大きな影響を与えることが報告されている以下の4パラメータを選択
  • 農薬溶解速度(kdiss)、水田土壌中の農薬の一次分解速度(kbio)、農薬の脱離速度(kdes)、農薬分配係数(kd)
  • これらの事前確率分布を作成(最大・最小間で一様分布)
  • 尤度は正規分布(E(θ)は観測値と予測値の平均2乗誤差)
  • MCMC で事後確率分布を作成

尤度を計算する前に採択されたパラメータでシミュを回す必要があります。大きなモデルを扱う実務向きではないですね。

実務で扱う地下水モデルで不確実性を検討するには、PEST一択なのでしょうか?確かに、市販のパッケージソフトには組み込まれており、使い勝手も良かった印象はあります。加えて、逆解析を古くから扱っているせいか、Bayes も良く出てきます。海外では当たり前なのかな?

もう少し、文献等を読んでみましょう。

2020年12月6日日曜日

change detection

Ohki et al.(2020)
Landslide detection in mountainous forest areas using polarimetry and interferometric coherence
https://earth-planets-space.springeropen.com/articles/10.1186/s40623-020-01191-5


ALOS-2 のデータ (Lバンド、空間解像度約6m)を用いた、土砂移動箇所検出の検討事例です。北海道胆振東部地震(2018年9月)と九州北部豪雨(2017年7月)を対象としています。
提案手法は change detection の範疇ですが、画像データを利用していません。波のデータのみです。そして最後はお決りの機械学習です(といっても決定木です)。

特徴

  • PolSAR、InSAR、DEM を組み合わせた解析

結果

  • すべての特徴量(InSAR、PolSAR、DEM)を用いた場合、kappa>0.6程度の精度。
  • 正確な検出には、少なくとも二重偏波(HHとHV)が必要。特に大きな入射角では四重偏波を推奨。
  • 精度向上には、より高い空間分解能、より高い時間周波数、より多くの方向からの観測が必要である。

留意点

  • Hサイト(北海道)は、試験地周辺の火山の火砕流堆積物に覆われている。そのため、胆振地震では浅い土砂崩壊が多く、比較的緩やかな斜面(<30°)で発生したと考えられている。Hサイトを用いて構築したモデルは,急峻な地形での多発したFサイトへの適用性が悪い。地質の違いに留意する必要あり。


最近の国交省さんの動向を見ていますと、今後、土砂災害分野における各種衛星データの利用は拡大して行くようですね。関連マニュアルも複数出ていますし、ALOS-3へ向けてのシンポジウムでも以下のような報告がされています。
https://www.pco-prime.com/alos-3_alos-4_sympo/pdf/1505.pdf 

機械学習を用いた change detection の精度は高いといえない状況ですが、いずれ向上し、災害時の人手不足問題も(change detection に関しては)回避できるでしょう。
検出アルゴリズムの構築以外は特に難しい作業はないので、備えておきましょう。

2020年12月5日土曜日

スパースモデリング

EXCELをよく使っていた頃、データに多項式をフィッティングする場面で何次式を使うかは感覚的でした。

低次だと合いませんし、高次だとノイズまで合わせてしまいます(いわゆる過学習)。ちょうどよいところを探す方法は、、、古くからありましたね。

LASSO
https://scikit-learn.org/stable/modules/linear_model.html#lasso
\begin{align*}
\min_{w}{ \frac{1}{2n_{\text{samples}}} ||y - X w||_2 ^ 2 + \alpha ||w||_1}\end{align*}

この損失関数はモデルの適合性を表す誤差項と、複雑さを表すペナルティー項で構成されています。αは求めたいパラメータwを制御する上位のパラメータ(ハイパーパラメータ)。αを大きくするほど大きなL1ノルムに対してペナルティーがかかる仕組みです。するとL1ノルムは小さくならざるを得なくなり、スパースな解が導かれる(が適合性は低くなりやすい)という仕組み。逆にαを小さくするほどL1ノルムを大きく取れるようになり、適合性が向上しやすい(が複雑で、過学習に陥りやすい)。αをどのように調整・採用するかがキモです。

その判断基準は以下の通り。

・情報量基準(AIC、BICなど)
・cross‐validation

情報量基準の利用では客観性に優れるものの感覚的にずれていると感じることがあるでしょう。一方、cv は理屈でなく現実的。ベースを統計にしているか、機械学習にしているかなどで好き嫌いが分かれそうです。

スパースを仮定した場合、非ゼロのwを増やすと組み込まれるXも高次まで多く含まれて適合性は良くなる(が複雑になる)といった表現が理想です。が、L0ノルムを扱うと組み合わせ(計算コスト)が膨大になるとのこと。そこでα*L1ノルムで代用しているそうです。
ちなみに、L1、L2 は機械学習でも多用されていましたね。回帰だと以下の通りです。

・L0正則化・・・非ゼロノルムの個数
・L1正則化・・・L1ノルム利用、LASSO回帰
・L2正則化・・・L2ノルム利用、Ridge回帰
・L1+L2正則化・・・Elastic Net

scikit-learn を利用すれば、これらのモデリングは容易。ありがたいですね。

 ==========================================
ラグランジュの未定乗数法を適用した場合の図や説明も web 上には多くあります。そちらの方が直感的で理解し易いでしょう。

関連内容はコチラ↓
https://phreeqc.blogspot.com/2019/06/blog-post_19.html
https://phreeqc.blogspot.com/2020/02/blog-post.html


2020年12月3日木曜日

GEE+ Sentinel-1

SAR データを利用した Change Detection を整理。

以前利用していたDocker は -vオプションを使用しなければ動きました。機能が追加されているようです。
https://phreeqc.blogspot.com/2019/04/sar-change-detection-on-google-earth.html

GEE コミュニティのチュートリアルにもありました。Part1は動きましたが、Part2の後半はダメでした。
https://developers.google.com/earth-engine/tutorials/community/detecting-changes-in-sentinel-1-imagery-pt-2

いずれも、GEE+ Sentinel-1。解析エンジンもデータも海外製です。データは政府事業の Tellus 提供より揃っています。
https://phreeqc.blogspot.com/2019/02/google-earth-engine.html
国内の複数個所で Change Detection を試行しましたが、動作も問題ありません。きちんと結果は出てきます。
扱っているモノがモノなので、まだ精度面では物足りません。が、手軽に加工できるところはありがたいですね。将来的にはアルゴリズムが向上し、精度面でも耐え得るモノになるでしょう。

Tellus でも ALOS2 のデータが公開され始めました。が、未だ直近1年分くらいは公開されていません。過去分も、2~3年前くらいまで。直近のデータがないのは辛い。
しかも、開発環境を割り当ててもらえないと、自由に加工できない状況は継続中。https://phreeqc.blogspot.com/2020/09/tellusar.html

まだまだ海外製には追いつけない感が漂っています。頑張って欲しいですね。


2020年11月23日月曜日

MCMC

MCMCを実装するライブラリはいくつかあるようです。

PyMC3:Theano
PyMC4:TensorFlow
Pyro:PyTorch
NumPyro:JAX

手元の環境で動いたのはPyroのみ。環境構築に関してはシビアなようです。

web上のサンプルを見てみますと、書き方はほぼPyMC3と同じでした。動かすのは簡単(traceplot がないのは残念。seaborn で書きました)。

が、写経し終わると理解できていない箇所が露わに。

 y_model = pyro.sample('y_model', dist.Normal(a*x + b, sd), obs=y)

ここの obs の計算内容がわかりません。ここがキモなのはわかっているのですが、マニュアルに詳細が書かれていませんし、ソースも追えません。尤度×事前確率にどのようにかかわっているのでしょう?
おそらく、理解している方々が必要として使われているので、学び方が逆なのでしょうね。

ま、他の部分は理解できていることが分かりました。メモしておきましょう。

 **************************************
20201123 MCMC メモ
 
データ駆動科学

  • データを出発点と考える=因果律を遡る
  • 従来:変形係数を確定(モデルを確定)→応力、ひずみがノイズを含み確率的分布(→安全率を含めた設計)
  • データ駆動:応力、ひずみが確定→変形係数が確率的分布(逆解析的発想)
  • データx,y(ひずみ、応力)が得られたとき、パラメータa(変形係数)が取り得る確率分布(事後確率)を求めたい。
  • 必要なのは、①パラメータaを与えたときにyが取り得る確率分布(尤度)、②aの確立分布(事前確立)

取得データD[x,y]

  • 確定x:ひずみx
  • 確定y:応力y


予測分布モデル (例)物理モデル+ノイズ
 物理モデル y=ax+b

  • 未確定a:平均μa,標準偏差σaとして設定(事前確率P(a)=N(μa,σa2)を仮定)
  • 未確定b:平均μb,標準偏差σbとして設定(事前確率P(b)=N(μb,σb2)を仮定)
  • 未確定σn:一様分布(U([n])を仮定)
    ←いずれもデータDを見て取りうる範囲を仮定

 ノイズ

  • 重畳ノイズを正規分布として設定(尤度P(y_model)=N(ax+b,σn2)を仮定


予測分布モデルから θ[a:変形係数、b:切片、sd:標準偏差]の事後確率 p(θ|D)を求める。

  • p(θ|D) = p(D|θ)p(θ)/p(D) ∝ p(D|θ)p(θ)
  • 事後確率 ∝ 尤度×事前確率
  • 事前確率、尤度に正規分布を仮定した場合、事後分布を解析的に求めることが可能(共役事前分布)。
  • しかし、一般的には計算が困難。→MCMC法などで帰納的に事後確率(最大)を求める。


マルコフ連鎖モンテカルロ法(MCMC法)

  • モンテカルロ法:乱数を使ったシミュレーション手法
  • マルコフ連鎖:現在の状態が、前時刻の状態のみに依存するモデル→メトロポリス法などの利用
  1. θ0を仮定
  2. 乱数εを加え、θ1を作成。
  3. 確率密度の比rを計算r=p(D|θ1)p(θ1)/ p(D|θ0)p(θ0)
    ノイズを正規分布と仮定した場合、p(D|θ) ∝ exp(-NE(θ)/σ2)
    あるθ0、θ1における exp(-NE(θ)/σ2)p(θ) を求める
    ←obs=D を利用 ※点での推定
  4. 確率密度の比rを評価
    r>1でθ1を採用
    r<=1かつ乱数U(1,0)<rでθ1を採用、U>=rでx0を採用
  5. 1~4の繰り返し←exp(-NE(θ)/σ2)p(θ)の分布形状を推定

 

他手法との違い

  • 最小二乗法:残差の平方和 E(θ) を最小化
  • 最尤推定:a,b,σについて導関数=0の位置(パラメータの値)を求める(global minimum とは限らない)
  • ベイズ推定:パラメータを確率的に扱う→分布(形)を求める 
******************************************
20201129
MAP推定:パラメータ値の同定のみの場合(パラメータの分布を求めない)
  • 作成した予測分布モデル+観測データ(obs)に対し、MAP(maximum a posterior :最大事後確率)推定でパラメータを同定
  • pyro.condition で同定パラメータをモデルに反映
 
20201202
変分推論:分布も推定
  • svi = SVI(model, guide, optimizer, loss=Trace_ELBO())でモデル構造を作成
  • step method でデータDをモデルに渡して調整←obsの役目
  • predictive で予測分布を作成(モンテカルロ近似)
概ね、obs の役割について理解しました。
数式は理解できていません。残念ながら後回しです。
 
20201204
数式もOK(おそらく)。
これで MCMC(メトロポリス法利用)は理解できたと思います。

2020年11月22日日曜日

PyMC

PyMCを使いたくて苦闘。

MCMCを理解する過程で必要となったのですが、これが動きません。

PyMC3 が動いたと教えていただいたライブラリVer.の組み合わせを真似し、仮想環境を構築したのですが、ダメでした。PIP を混ぜてPyMC4にもしてみましたが、ダメ。Anaconda を最新版に入れなおしてもダメ、gcc 関連のパスを通してもダメ。お手上げでした。

よく引っかかったのは Theano。これ、Deep Learning にも流用できる数値計算ライブラリとして有名だったのですが、既に開発が終わって時が経っています。
ちなみに、PyMC4 は TF を利用。TFも2年ほど前はよく利用されていましたが、今は PyTorch に首位の座を奪われているように感じます。

バックエンドが変更になるたびにコードが動かなくなるのは避けたいですね。Docker しかないかな?

2020年11月17日火曜日

DO CONCURRENT with GPUs

Accelerating Fortran DO CONCURRENT with GPUs and the NVIDIA HPC SDK
https://developer.nvidia.com/blog/accelerating-fortran-do-concurrent-with-gpus-and-the-nvidia-hpc-sdk/

昨日のブログです。
OpenACC のデータ転送が面倒で困っていたところ、この記事を見かけました。

 All data movement between host memory and GPU device memory is performed implicitly and automatically under the control of CUDA Unified Memory.

これでACCと同程度のパフォーマンスが出たらラッキー。頻繁にメモリへのアクセスがあるので過剰な期待はしていないのですが、いくらか早くなれば嬉しい。早く試したい!

Docker はまだ 20.9 (SDK 単体も)。公開まであと1か月くらいでしょうか?

**************************************

20201217追記

先日、Docker イメージが公開されました。
早速試してみたところ、-stdpar=multicore では OpenMP と同じ結果、速度でした。
が、-stdpar=gpu では発散。うまく計算しません。
コンパイル情報を見てみると、gpu では DO CONCURRENT 内の DO LOOP も自動並列化で GPU に載せられていました。これかな?

2020年11月12日木曜日

NVIDIA HPC SDK

OpenACC を利用しようと考えていました。

ACC の実装は比較的手軽だそうですが、場合によっては CUDA に並ぶことがあるとのこと(以前、チャレンジしたときは失敗しました)。
私はプロではないので、C & CUDA によるゴリゴリチューニングはできません。で、OpenMP からの ACC で、あわよくば狙いです。

OpenACC を gcc がサポートしているようなので使おうかな?などと考えながら、いつもの PGI Community Edition をチェック。以前のVer.と変わらないなあ、と製品版の価格をチェック。すると、以下の表示。

PGI Professional 製品の販売は終了しました。(詳細)
( ゚Д゚)?!
3月で販売終了していたとのこと。知りませんでした。
PGI Compilers&Tools は NVIDIA ブランドの NVIDIA HPC SDK ソフトウェアとして、新しく生まれ変わりました。
NVIDIA HPC SDK
https://developer.nvidia.com/hpc-sdk

ついに PGI の冠が取れたのね。
プロファイラーも以前のまま含まれています。しかも、GPU Cloud で配布されています。
https://ngc.nvidia.com/catalog/containers/nvidia:nvhpc

はい、素晴らしい。さすが NVIDIA。
もうこれは、Linux + Docker + HPC SDK 一択でしょう。

最新のアーキテクチャに対応したコンパイラーが無償で配布されるのも素晴らしいのですが、これ、Win版の Fortran コンパイラーとしても使われだすのではないでしょうか?

俄然、やる気になりました。使ってみましょう。

**************************************
20201113追記

使ってみました。
dockerイメージを走らせただけ(nvidia-docker を使わなくなっていました)。その上で何の問題もなくコンパイルできました。
計算中に GPU を使っています。kernels 構文を差し込んだだけですので、遅いのですが。
チューニングには時間を取られそうですが、制御は拍子抜けするくらい簡単でした。ありがたいですね。

2020年11月3日火曜日

SPH 基本5:長所と短所

SPHスキームの長所と短所

圧縮性(SPH法)
長所:陽解法のため逐次代入計算により時間更新を行う。計算が簡単。
短所:粒子の重なりや抜けが生じやすい。結果として数値的不安定性を招きやすい。
→安定させるためにタイムステップを短くする必要がある。(計算負荷増大)
→人口粘性の導入が必要。ただし、適度に設定しないと訛りすぎる。

弱圧縮性(WSPH法)
長所:陽解法。疑似非圧縮性流体の計算が可能。
短所:計算の安定性のために0.5%程度の圧縮性を許容することが必要。
→体積保存性に欠陥

非圧縮性(ISPH法)
・SPH法の微分演算子とMPS法の計算アルゴリズムの組み合わせ。逆も可。半陰解法。
長所:非物理的な人工粘性を導入する必要がなく安定した計算が実行可能。
短所:圧力場の解に激しい攪乱を伴う
→インターバル平均で攪乱は消失。物理ベースCGなどで利用可能。
→防波堤に作用する衝撃波圧など瞬時の値を扱う問題には使えない。
※陽解法、半陰解法のいずれにおいても、安定性向上に有効なのが粒子再配列(高精度粒子法として分類される)
・XSPH
・Particle Shifting:J=-v∇C
・Optimized Particle Shifting

 

SPH 基本4:圧力場の振動

弱圧縮性SPH数値スキームでは、圧力場がスプリアス振動を示します。

状態方程式を用いた場合、密度の変化が非常に小さくても、圧力の変化は指数関数的に生じます。それらは、物理的にはありえないような巨大な圧力を生じさせます。

対策
Density filtering: Shepard Filter (0th order), Moving Least Squares (1st order)
人工拡散項の追加:δ-SPH, Fourtakas et al. (2019)
Riemann solver:Godunov SPH

δ-SPH
WSPH法の数値的不安定(引張不安定)の制御手法として、δ-SPH法が近年よく用いられる。δ-SPH法では密度勾配をTaylor級数の一次精度で算定することにより拡散項の評価精度を向上し、自由表面の跳躍を抑制している。
後藤仁志「粒子法 連続体・混相流・粒状体のための計算科学」より
Molteni, D. and A. Colagrossi. (2009)A simple procedure to improve the pressure evaluation in hydrodynamic context using the SPH.
S.Marrone et al.(2011) δ-SPH model for simulating violent impact flows

Fourtakas et al. (2019)
Local uniform stencil (LUST) boundary condition for arbitrary 3-D boundaries in parallel smoothed particle hydrodynamics (SPH) models

Godunov SPH
Riemann問題の解を適用し、必要最低限の粘性を自動的に導入する手法。あらたなパラメータ設定が不要。
Inutsuka et al. (2002) Reformulation of Smoothed Particle Hydrodynamics with Riemann Solver

2020年11月2日月曜日

SPH 基本3:精度

SPH 空間補間が、平滑化長を h とすると、 𝑂(ℎ2) に等しい精度で収束。

1Dの場合、

\begin{align*}
\left\langle A\left(r\right)\right\rangle=\int_{\mathrm{\Omega}}{A\left(r^\prime\right)W\left(r-r^\prime\right)dr^\prime}\end{align*}
A(r’) を r 周りでテーラー展開すると
\begin{align*}
A\left(r^\prime\right)\cong A\left(r\right)+\frac{\partial A\left(r\right)}{\partial r}\left(r-r^\prime\right)+\frac{1}{2}\frac{\partial^2A\left(r\right)}{\partial r^2}\left(r-r^\prime\right)^2+・・・\end{align*}
SPHでは

\begin{align*}
\left\langle A\left(r^\prime\right)\right\rangle=A\left(r\right)\int_{\mathrm{\Omega}}{W\left(r-r^\prime\right)dr^\prime}\end{align*}

\begin{align*}+\frac{\partial A\left(r\right)}{\partial r}\int_{\mathrm{\Omega}}\left(r-r^\prime\right)W\left(r-r^\prime\right)dr^\prime\end{align*}

\begin{align*}+\frac{1}{2}\frac{\partial^2A\left(r\right)}{\partial r}\int_{\mathrm{\Omega}}{\left(r-r^\prime\right)^2W\left(r-r^\prime\right)}dr^\prime\end{align*}

\begin{align*}+・・・\end{align*}

 1.正規化:影響域で積分すると1になる。

\begin{align*}
\int_{\mathrm{\Omega}}{W\left(r-r^\prime\right)dr^\prime=1}\end{align*}
 2.偶関数:中心に対し点対称

\begin{align*}
W\left(r-r^\prime\right)=W\left(r^\prime-r\right)\end{align*}

\begin{align*}
\int_{\mathrm{\Omega}}\left(r-r^\prime\right)W\left(r-r^\prime\right)dr^\prime\end{align*}

\begin{align*}=\int_{\mathrm{\Omega}}{\left(r-r^\prime\right)W\left(r^\prime-r\right)dr^\prime}\end{align*}

\begin{align*}=-\int_{\mathrm{\Omega}}\left(r-r^\prime\right)W\left(r-r^\prime\right)dr^\prime\end{align*}

\begin{align*}=0\end{align*}

3.\begin{align*}q=\frac{r-r^\prime}{h}, w\left(r-r^\prime\right)=\frac{\alpha}{h}f\left(q\right), 𝑓: ℝ+ とすると\end{align*}

\begin{align*}\alpha\int_{\mathrm{\Omega}_q} f\left(q\right)dq=1\end{align*}

\begin{align*}dr=dq・h\end{align*}

\begin{align*}
\alpha\int_{\mathrm{\Omega}_q}{\left(r-r^\prime\right)^2f\left(q\right)dq}=\alpha h^2\int_{\mathrm{\Omega}_q}{q^2f\left(q\right)dq}\end{align*}
1, 2, 3より、

\begin{align*}
\left\langle A\left(r^\prime\right)\right\rangle=A\left(r\right)+\frac{1}{2}\frac{\partial^2A\left(r\right)}{\partial r}\alpha h^2\int_{\mathrm{\Omega}_q}{q^2f\left(q\right)dq}+・・・\end{align*}

\begin{align*}=A\left(r\right)+O\left(h^2\right)\end{align*}


2020年11月1日日曜日

SPH 基本2:カーネル

 カーネル関数の基本的な性質

  • 正規化:影響域で積分すると1になる。
  • compact support:影響域が有限。区間の端では関数値が0になる。
  • 非負性:kernel は定義域内では非負である。
  • 単調減少性:中心から離れるに従って単調に減少する。(近接している粒子ほど強く影響しあう)。
  • Dirac のδ関数への収束性:smoothing lengthが0となる極限で Kernel が関数 f(x) 自体に一致。
  • 偶関数:中心に対し点対称
  • 平滑性:kernel とその微分が連続かつ滑らかな関数形状を有する。(計算が安定) 
後藤仁志「粒子法 連続体・混相流・粒状体のための計算科学」より

2020年10月31日土曜日

SPH 基本1:補間式

汎用スカラー関数fのSPH補間式(連続・離散レベルの両方) 𝑓: ℝ3→ℝ

SPH法では、f(x)はx近傍の領域Ωにおける積分として記載される。
位置 x での物理量f(x)はディラックのデルタ関数を使うと、

\begin{align*}
f({x})=\iiint_{\mathrm{\Omega}}{f(\mathbf{x}\prime)\delta(\left|{x}-{x}\prime\right|)dV}\end{align*}

h→0 でW(x,h)→δ(x) となるような関数を使って近似すると、
\begin{align*}
\langle{f({x})}\rangle=\iiint_{\mathrm{\Omega}}{f(\mathbf{x}\prime)W(\left|{x}-{x}\prime\right|,h)dV}\end{align*}

これを離散化すると、
\begin{align*}
\langle{f({x_i})}\rangle=\sum_j{f(\mathbf{x_j})W(\left|{x_i}-{x_j}\right|,h)V_j}\end{align*}

これは、中心点xiの物理量を周辺の離散点xjの物理量に重みを乗じて合計した値とみなせる。


2020年10月27日火曜日

オンラインセミナー

コロナ禍により、セミナーが一斉にオンライン化してます。

今のところ、メリットしかありません。享受しています。

・場所を問わない。移動中でも参加できる。
・疑問に思ったら、即チャットで質問可能。
・施設も紙の配布資料もなくなるためか、無料のコースが多い。
・後で録画を見直せる。
海外のセミナーも受講できる!

先週、イタリアのパルマ大学の方が SPH 関連のオンラインセミナーを開かれました。以前から受けてみたいなあと考えていたのですが、欧州は遠い。さすがに私費では難しい。しかも適当な英語しか話せない。ハードルが高かった。
そのような私がオンライン開催を見つけた時にはすぐに申し込んでいました。SPH 関連の専門用語なら聞き取れるだろうし、数式を追えばついていけるだろうと。イタリアなのでネイティブでないだろうから聞き取りやすいだろうと。まったくのポジティブシンキングです。

申し込んだ後、なぜかしばらく忘れていました。さらに案内メールも見逃していました。で、気づいたのが当日。1日目の録画のURLが送られてきて開催日初日だったことに気づきました。
次の日から出張と重なっていたのですが、道中に1日目の動画を確認。残る4日も、ホテルで受講。幸い時差の関係で開始時刻が夕方の16時~17時でした。最初の1時間くらいは受講できない日もありましたが、毎回15分くらい遅れて始まるので、それもありがたい。後で動画を見ておくことで、十分にフォローできました(最後の方は数式すら半分以上わからなかったのですが)。
宿題も土日をかけて提出。録画が役に立ちました。
内容的に申し分なく、しかも無料で受講させていただいたことには感謝しかありません。

海外セミナー+チャット+録画は、私にとって最強の組み合わせでした。大きなハードルが、一気に解消されました。
欲張りすぎたのか、頭がパンパンな感じ。でも、今日も国内のセミナーに申し込んでしまいました。

オンラインセミナー、ありがたいの一言です。

****************************************
20201031追記

final examination に合格し、3ECTSをいただきました。Grazie!


2020年10月18日日曜日

equifinality problem

昨日の図書では、数値モデルに同位体をトレーサーとして融合する長所・短所が示されていました。
長所:モデルの信頼性向上
短所:校正すべき未知パラメータの増加、固有の不確定性(δ値の時空間的な不均質性など)

今まで、短所には目を向けず長所ばかりを見ていました。が、ある程度その見方も許容されるようです。一般には同位体データを含めて校正・検証されたモデルの方が信頼性が高いと言え、その融合に期待されているとありましたので。

この長所ですが、具体的には等結果性問題への対応ということになります。
等結果性(equifinality):異なる構造・パラメータの組み合わせが同じ結果を示す
浸透流解析では、支配方程式で右辺の各項の値が異なっても、合計が左辺の値として一致することと同義です。

これはわかりやすい話です。
一般的な3次元浸透流解析では、地質構造(ジオメトリ)を1つ決定した後、水理パラメータの調整で水位や流量等の再現を試みます。これは、ジオメトリの調整・再作成に大きなコストを必要とするため、それを固定した後、相対的に安価なパラスタにて再現を試みているということです。つまり、構造の不確実性をパラメータに一手に担わせて調整していることになります(「地下水モデル」にも書かれていましたので、世界的な流儀でしょう)。
本来は文献調査や地質調査結果から複数の3次元地質構造の可能性を地質屋さんは推定します(通常は、その中の可能性の高いと考える構造を地質屋さんが1つ選択・提示しているにすぎません)。ところが、どのジオメトリでも浸透流解析では観測値を概ね再現できます。パラメータ側でその差を調整できるのです。つまり、同じ結果を再現できるモデル(ジオメトリとパラメータの組み合わせ)が複数誕生するのです。これが等結果性。
しかしながら、将来発生するイベント(トンネル掘削など)を入力した場合、どのモデルも同じ予測値を示すとは限りません。再現計算完了時点では、どれが実現象に最も近いモデルなのかわからないのです。
そこで、それらのモデルにδ18Oなどを乗せて流し、観測位置で得られたδ値を再現できるかチェックすれば、モデルの良否が判明します。時空間的に変化しない疑似トレーサーでも可能でしょう(時空間的に変化しやイオンを選択する場合には、Geochemical Modeling との融合が必要です)。これが等結果性への対応です。

「多角的に校正・検証したモデルは、より正確な予測値を出すことが期待できる」ということです。

環境同位体の利用方法

とある大学の書店に立ち寄りました。

理学部のないキャンパスなのか、地質系の図書はありませんでした。が、環境系のコーナーに以下の図書を見つけました(こういった配置、時々あります)。

山中勤「環境同位体による水循環トレーシング」

δ18O などの利用法は、文献からの知識や経験的に身に着けたように思い返されます。ですので、この図書の5章「地下水・湧水・温泉」の区分は系統的とまでは言えないものの、頭の整理に役立ちました。 

********************************************
1.涵養標高の推定:降水のδ値は標高が高くなるほど低くなる

2.流動系区分:δ値の高低で平面的な流動系区分

3.混合割合評価:
涵養源の寄与率評価:δダイアグラム上で河川水・降水・田面水などの混在割合を区分、GIS上でプロット
温泉水:
δダイアグラム上で温泉水・天水(近隣の地下水や河川水)・非天水などの混在割合を区分

4.地下水流動モデルとの融合:
直接利用:δ2H, δ18O をトレーサーとして使用
間接利用:Particle Tracking から算出した河川水寄与率を、安定同位体やイオンの混合計算結果と比較

*******************************************

特に浸透流モデルのブラッシュアップにて、移流分散と粒子追跡の利用を同じくくりにできるか考えようとしたことがありませんでした(頭が固い)。
確かに、このような説明は国内で見ません。聞けば地下水のプロ?でも知らない方がいらっしゃるようです。土木分野で普及していないだけかもしれませんが。

他にも、経験的に知っていたことが明文化されており、あらためて知識として習得できたように思える個所がありました。新たな図書として、手元に置いておきましょう。


2020年10月14日水曜日

微動計3:波形のリアルタイム表示

微動計の設定に目途がつきました。

あとは測定しながら波形を確認したい。
これは簡単だろうと思って手を付けたのですが、ハマりました。

PC と微動計を無線で接続し UDP で吐き出されるデータを受信するだけなのですが、表示ソフトの画面に出てきません。ping は通っていますし、ポートもソフトがつかんでいるのですが画面に動きなし。ソフトが悪いのか?とPythonで試すも改善せず。有線ではどうか?これもダメ。
他の PC では?と、OS をクリーンインストールしたばかりの PC を用意しましたが、こちらもダメ。
念のため、もう1台の PC で試すと、
あっさりできました。原因不明。うーん。

その1台で波形を見ながら軽く衝撃を与えてみます。サンプリング周波数を変えながら、最適な測定条件をチェックします。やはり、確認しながらの作業は効率的です。

ある程度見通しが立ったので、次回、プレテストへ。
綺麗な波形が取れることを期待しています。


2020年10月13日火曜日

微動計2:オーバーサンプリング

微動計内では、これらのフィルターをいつ適用しているのか?

キーワードはアンチエイリアシング、オーバーサンプリング、デシメーション。
メーカーさんに伺うと、製品内の処理手順は以下の旭化成エレクトロニクスさんの解説と同じ手順でした。
https://www.akm.com/jp/ja/technology/technical-tutorial/basic-knowledge-adc/oversampling-advantage-1/
https://www.akm.com/jp/ja/technology/technical-tutorial/basic-knowledge-adc/oversampling-advantage-2/

オーバーサンプリングしておけば、測定対象の周波数帯に含まれるエイリアシングの影響を小さくすることができる。
その後、測定対象よりも高い周波数帯を、再度、デジタルフィルタ(LPF)でカット・デシメーションすればOK。
https://www.toyo.co.jp/mecha/casestudy/detail/id=12894

実際、複数の微動計でサンプリング周波数を変え、同じ衝撃をはかると異なる波形が得られます。低周波数に抑えて不思議な結果に悩むなら、最初から高周波にして物事を単純化しておくべきなのでしょうね。特に私のような入門者にとっては。

どちらかというと、安価で安定したセンサーを作る目的で発達した技術のようです。
これ以上は踏み込みませんが、奥が深そうな世界です。

2020年10月11日日曜日

微動計1: 信号処理

微動計を用いた衝撃の測定に呼ばれました。

今まで他の方が頑張っていました。が、得られたデータがイマイチだったということが判明。で、地盤を対象に微動計を使っていた私も御一緒することに。

どうも、微動計に備わっていたフィルタの周波数特性をうまく活用できていなかったのが一因のようでした。
その選定から見直そうとしたのですが、私も応用が利くほど詳しくありません。調べていくうちにわからない用語が出てきました。標準モード、群遅延特性、デシメーション。
使うツールの詳細を理解していないと、正しい調査ができないのは共通。いままでプロの指導でうまくいっていただけのことであり、場を変えると途端にダメになるのは同じ。
ちょうど週末で時間があり、この機会に理解することにしました。

用語自体は信号処理分野で広く使われているようです。いくつかの図書を見ましたが、基礎中の基礎といった感じでした。
わかりやすかったのが図書よりも以下のサイト。確かに、音響ととらえると理解しやすい(ホールの設計者はこのようなことまで考えていらっしゃるのですね!)。
以下、解釈した内容を( ..)φメモメモ。


**************************************
ヤマハサウンドシステム
短期集中連載 超解説FIR!第1回~第3回
https://www.yamaha-ss.co.jp/published-articles/journals-04.html

1.用語
1-1.周波数特性・位相特性・群遅延特性
入力したインパルス信号の応答を周波数領域で見る(フーリエ変換)と各周波数のエネルギー応答がわかる。これが周波数特性。
各周波数が時間的にどのくらい遅れて到達しているかを見たものが位相特性。
位相を周波数で微分したものが群遅延特性。
周波数とその位相の変化量(遅れ)が比例していれば群遅延は一定(直線位相特性)。  
ex.
1Hz→1波長(2π、1秒)
2Hz→2波長(4π、1秒)
3Hz→3波長(6π、1秒)

1-2.インパルス応答、伝達関数
システム(ブラックボックス)にインパルスを入力した時に出力された信号がインパルス応答。
インパルスは全周波数で同じエネルギー量持っているため、インパルス応答を見ただけでブラックボックスの特性=伝達関数が分かる。
インパルス応答を周波数の観点で見ると周波数特性が分かり、周波数ごとの到達時間を調べれば位相特性が分かる。
インパルス応答、伝達関数、周波数特性、位相特性等は、データの見方を変えただけ。

2.特性を作る(リバーブエフェクトなど)
2-1.FIR フィルタ
FIR の Tap にある乗算器の係数はフィルターのインパルス応答そのもの。
係数操作でインパルス応答を人工的に加工することが可能。IIR のように最小位相として位相が変化するインパルス応答などを作ることも可能。
目的とする周波数カーブを描いてそれを逆フーリエ変換すれば、カーブのインパルス応答が得られる。
そのインパルス応答を乗算器の係数に入れることで目的の周波数カーブのフィルターを作ることができる。
FIR フィルタは入力信号(サンプルデータ)とインパルス応答波形の2つを畳み込んで新たなデータを作り出す。つまり「畳み込みをするための演算器」。

2-2.FFT と窓関数
入力信号も IR も共に FFT して周波数領域のデータに変換。
畳み込み積分は周波数領域に変換すると単純な掛け算として非常に高速の計算が可能。
その結果を逆フーリエ変換(IFFT)して元の時間領域のデータに戻す。
ただし、FFT を行なうため処理するデータを切り出す際にサイドローブ(雑音)が発生。
窓関数でサイドローブの軽減処理をしているが、軽減しようとすると反対に周波数分解能が悪くなるという二律背反の関係がある。
周波数分解能:   〇     △      × 
サイドローブ軽減: ×     △     〇
窓関数:    ハミング  ハニング  ブラックマン

**************************************

インパルスにδ関数が使われているのは当然なのに新鮮。つながりますね。 

 次は、これらの処理順序について。


文献 その2

ためがちな文献。

最近は 土石流や SPH 関連をよく読んでいたので、それ以外は放置でした。
以下、気になったモノと箇所を書きとめておきます。

**************************************
K-NET・KiK-netのPS検層記録に基づくVs・VPおよび深さの関係
物理探査/ 72 巻 (2019)
Vs-1000m/s 以下の領域では、Vp=500m/s or Vp=1500m/sに大別される。
土質によらない。不飽和、飽和による。
地下水の既知、未知に応じた回帰式を提案。
ポアソン比:K-NET:0.45以上(深度20mまで)、KiK-net:0.33(ばらつき大)

レイリー波の分散曲線の近似計算法の提案と地下構造推定への応用
土木学会論文集/1997 巻 (1997) 577 号
多層水平地盤構造における基本モードレイリ-波の位相速度の分散曲線をS波速度と層厚のみのパラメータで近似。

深層学習による崩壊・非崩壊地の自動判読手法の開発
日本地すべり学会誌/56 巻 (2019) 5 号
SonyNNC使用。深層崩壊地のタイル画像作成・利用

地すべり運動を磁気特性から探る-リングせん断試験の供試体を用いた帯磁率および残留磁化実験による検討-
日本地すべり学会誌/56 巻 (2019) 5 号
リングせん断試験により供試体の磁性粒子が円周方向に配列。磁気特性からすべり運動の方向を再現できる可能性がある。

天然ダムの決壊による洪水流下の予測と対策
砂防学会誌/45 巻 (1992-1993) 1 号
天然ダムの決壊時刻を推定する簡易推定図
単位幅ピーク流量と単位幅貯水面積との関係による単位幅ピーク流量の簡易予測図
(1)天然ダム決壊時のピーク流量の支配的な指標はダム高さと河床勾配である。
(2)天然ダムが高いほど単位幅ピーク流量が大。河床勾配が大きいほど単位幅ピーク流量小。
 
ネット上で引っかかった文献
塩化ナトリウムの地下への浸透による水質変化に関する研究
凍結防止剤の散布
高速道路からの距離と塩化物イオンの関係:150m付近まで基準超過
安山岩と花崗岩領域の水質変化の違い
 
昨今の災害事例から見た道路斜面防災の在り方と留意点
治水3法から土砂災害防止法までの変遷・・・基本的にその対象が人家等の保全対象を含むもの。道路には対象から外れており、自衛手段が求めらる。
切土のり面の崩壊
①供用後5年以内
②過去に変状を来したのり面
③急速施工
④発生しやすい地質
⑤80%以上が深度5m、平均N値23未満。
⑥背後にリニアメント(断層、素因)
盛土のり面の崩壊
①盛土を通るリニアメント
②ボトルネック地形
③道路を挟む上流側にレベルバンクや集水地形
④崩壊の半数は供用後5年以内。20年以降でやや増える。

 

2020年10月9日金曜日

geochemical modeling と地下水の流路

Isotope hydrology and hydrogeochemical modeling of Troodos Fractured Aquifer, Cyprus: the development of hydrogeological descriptions of observed water types
https://www.sciencedirect.com/science/article/pii/S0883292720302729

安定同位体と geochemical modeling を組み合わせ、 地下水の涵養源と流路(水質の進化経路)について推察された報告です。

geochemical modeling は PHREEQC を利用。地下水の起源、流路、岩石との反応プロセスを示されています。他の複雑な帯水層システムにも適用可能とされていますが、ま、妥当でしょう。
浸透流計算を外しているため、PHAST よりお手軽です。この程度なら国内でも稀に受け入れてもらえるかもしれません。

残念ながら、国内にはこのような手法を扱う研究者がほとんどいらっしゃいません。研究者のお墨付きがなかったり、基準を作成してもらえないと、なかなかお客様も首を縦に振っていただけません。
そのためか技術者側もフレームワーク自体を理解しようとせず、「水質を調べて区分しました」「地下水年代を調べて区分しました」「〇%づつ混じっている可能性があります」程度で終わっているケースをよく見ます(こちらも残念です)。ま、利益に直結するルートが開拓されていないにもかかわらず、それに時間と頭を使うのはビジネスにとって省くべき無駄というのも理解できます。海外では図書レベルなのですが。

PHREEQC にハマってから10年以上が過ぎました。次の10年も同じようなものなのかもしれません。どこかに生かせる職業はないでしょうか?


2020年10月6日火曜日

SPlisHSPlasH

SPlisHSPlasH
https://github.com/InteractiveComputerGraphics/SPlisHSPlasH

先日の PySPH よりも式が豊富です。
しかも GUI で入れ替え可能。素晴らしい。

手元のソースの状態方程式や密度の時間微分の計算式を何度か入れ替えて試算していたのですが、これがあればずいぶん楽になるのでは?

などと考えていましたが、甘かった。
example はきちんと動くものの、自分でモデルを作る方法がわかりません。複雑な粒子配置はどのようにすれば実現できるのでしょうか?うーん。

「物理ベースシミュ」と説明されていますし、「泡」への対応も済んでいるようです。Maya や Blender のプラグインも用意されています。CG分野ですね。
この分野では、やはりリアルタイム計算を目指すのでしょう。速度面への対応は以下の通り。確かに、速い。
・neighborhood search on CPU or GPU
・supports vectorization using AVX

うーん。このまま寝かせるには勿体ない。
どこかに取っ掛かりがないでしょうか?


2020年10月4日日曜日

PySPH

PySPH を見てみました。
https://github.com/pypr/pysph

使用されている式は以下の通り。境界条件がイマイチはっきりしないのですが、それ以外は一通り(固体まで)そろっています。複数の式が出典とともに並べられていますので、文献を見比べたり探さなくて済みます。これだけでも参考になりますね。
https://pysph.readthedocs.io/en/latest/reference/equations.html

サンプルを動かすだけなら問題ないのですが、新たなジオメトリを作る方法がわかりません。以下のチュートを動かして少し理解できた程度でした。
https://github.com/pypr/pysph/tree/master/docs/tutorial

ブロック感覚で式を組み合わせる方式なので手軽です。生かしたいですね。

 

2020年9月25日金曜日

SPH と状態方程式

SPH のソースを書き換えようとして、ふと気づきました。

粒子法では計算中に水面を取得するのが面倒?

計算中にどの粒子のどのあたりを自由表面とすべきかがわかりません。なので他の手法のように水深から容易に静水圧を出せません。いえ、post 処理で使っている方法があるので可能ではあるのですが、計算量が多くなります。
で、どうしているかというと、状態方程式から出すようです。

例えば、
「Modelling multi-phase liquid-sediment scour and resuspension induced by rapid flows using Smoothed Particle Hydrodynamics (SPH) accelerated with a Graphics Processing Unit (GPU)」
「Smoothed Particle Hydrodynamics: A Meshfree Particle Method」ソース:p408, 409, 400, 本文:137

P=B((ρ/ρ0)^γ-1)
P=c^2ρ

簡単な流儀ですが、慣れるまで時間がかかりそうです。

これまで MPS の解説に触れる機会はあったものの、SPH は多くなかったように思われます。良書が国内になく、また過去の経緯もあり個人的には国内の SPH に良い印象を持っていませんでした。
が、先月から文献や洋書を読み進めていたところ、 比較的新しい研究が多いことに気づかされました。実装も他の手法に比べ容易ですし、プレはCAD + DualSPHysics + Python で済みます。これ、使えるのでは?と思うようになってきました。海外の動向に目を向けず、食わず嫌いだったのかもしれません。

あとは経験でしょうか。
しばらく手を動かしてみましょう。

2020年9月23日水曜日

Flux

D. Violeau "Fluid Mechanics and the SPH Method"より。

数式の展開はなんとかついていけますが、説明がイマイチわからない。でも綺麗なのでメモ( ..)φ。

∂B/∂t + gradB・u=(1/ρ)div(Jb gradB) + Sb

Variable    Flux  Source
B  Jb gradB    Sb
u σ    0
ek σ・u     -f
eint  Jt gradT f
e σ・u+JtgradT    0

応力テンソルが運動量フラックスに、 σ・u がエネルギーフラックスに。

ここから σ を書き換えてナビエ・ストークスへ。結局、圧力をどう計算するかだけの違いだけです。そう見ると、地表流や河床変動計算の流体部分、津波、土砂移動、さらには固体まで、すべて同じ数式に見えてきます。
固体も連続体力学のスキームをあてはめるだけだったとは。多くの文献で書かれていますので、古くからある考え方なのでしょう。

読んでみるものです。

 

2020年9月20日日曜日

Windows Remote Desktop + OpenGL

Windows の Remote Desktopで OpenGL を利用するソフトを起動したい場合があります。

といっても、普段からすべてのソフトに対して OpenGL の利用をチェックしているわけではありません。大抵はリモート接続後に起動しないとわかって初めて気づきます。で、起動しなくて困ると(情けない)。 

 調べてみると解決方法はあるようですね。
https://developer.nvidia.com/designworks
https://blog.techlab-xe.net/rdp-opengl-app/ 

今日は後者で対応しました。接続し直すと起動していたので一安心。ヨカッタ。

普段、把握しているソフトは立ち上げてから離席することで対応していました。が、NVIDIA さんの方法なら根本的に解決しそうです。

連休明けに試してみましょう。

**************************************

20201004
前者はインストールするだけ。こちらの方が簡単でした。


2020年9月19日土曜日

インターンシップ

インターンシップに参加されていた学生さんが帰られました。

夏休みを利用する学生さんが多く、先月からのリストを見るとかなりの人数になっていました。単位の出る1週間程度の方が多いそうですね。知りませんでした。

先月、私も対応しました。
これまで学んできた知識と実地がつながると、学生さんは目が覚めたようでした。アハ体験ですね。最初の希望分野とは異なる方向だったのですが、今後の選択肢が増えたようでした。

社会人にとっては毎年のことでも、学生さんにとっては生涯一度、人生を左右する経験になるかもしれません。積極的に多くの会社に出向いて社員を活用してほしいですね。知ることでさらなる勉強の必要性に気付かれるかもしれません。まだ学んでいない世界の方がはるかに大きいのです。眼前には大きな可能性が広がっているのです。

インターンシップ、結果的に良い選択へのサポートになればと考えます。


2020年9月12日土曜日

トポロジー

トポロジーに高次元を縮約して図化する研究があると話をされたお客様がいらっしゃいました。

記憶にあるのは超立方体や3次元トーラス、ドーナツ=コーヒーカップの印象程度。いえ、ドーナツ見せられてコーヒーカップと一緒でしょって言われても、たいていの方は腑に落ちないでしょう。

あらためて見直さないといけないなと思いつつ、本屋に行っても素人には難しい図書ばかり。その中でもわかりやすそうで読み始めたのが以下の図書。
市原一裕「低次元の幾何からポアンカレ予想へ」

トポロジー:位相幾何学
多様体:ある点の近傍ではn次元のユークリッド空間のある点の近傍と同相となるような、ユークリッド空間内の図形:円は1次元多様体。

トポロジーを利用した設計というのもあるそうです。
私の分野で利用する機会はなかったのですが、単に気づいていないだけかもしれません。情報は仕入れておきましょう。


2020年9月10日木曜日

SPH の図書

集めていた SPH の文献を読み終わりました。

今回は、DualSPHysics のビンガム流体、個体などを扱った文献を中心に10本程度。
最初はあちこち巡りながら必要部分のみを読んでいたのですが、最終的にはどれも一通り読み終わることになりました。部分的に異なる数式が使われているものの、基本は同じです。バリエーションが少なく、FEMほど複雑ではないので、同じような書き方になるのでしょう。

‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐

連続の式:∂ρ/∂tを離散化

運動方程式∂v/∂tを離散化

 σ 固体:弾塑性のフレームワーク
 (降伏関数、Dマトリックス、ひずみ速度テンソル、スピンテンソルなど)
 σ 液体:流体のフレームワーク
 (粘性、HBPモデル、浸透力(ダルシー)、浮遊砂など)
 状態方程式
 人工粘性
 レイリー減衰
 XSPH

境界条件

‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐

役立ったのは、 連続体力学の知識(あと、土質力学、線形代数など)。
見た目粒子ですが、解くのは連続体ですから当然です。が、つい流体のイメージが先行し、先日まで実施していた河床変動計算のイメージに寄っていました。いえ、NSは共通なんですけどね。

今日、1か月前に頼んでいた図書が英国から到着。
「Smoothed Particle Hydrodynamics: A Meshfree Particle Method」
楽しみにしていたのですが、待っている間に主たる内容をほぼ理解していました。高かったのですが。
ま、手触りが良いし、ソースコードがついていたので良しとしましょう。

****************************** 

 20200912追記

2020年9月6日日曜日

空間の離散化

浸透流解析でメッシュと精度の関係が何かに書かれていたはずと、他支店のプロ?に聞いてみました。

結果は御存知ないようでした。

資料を探して出てきたのが、昨年購入した「地下水モデル」の5章でした。勘違いしていた点もあるのでメモ( ..)φ

**************************************

水平方向の節点間隔
・節点間隔に対する感度解析を行うことは滅多にない
・河川流出3λ、井戸流出4λ、節点間隔λ以下(理想は0.1λ)
・構造格子:1.5倍ルール

鉛直方向の接点間隔
・ほとんどの状況では、鉛直方向の節点間隔が水平に比べてかなり大きくても、必ずしも数値誤差が生じるわけではない
→異方性により水平方向に流れが向くため。

・鉛直方向の間隔が水平方向の√(Kh/Kv)倍で等方性メッシュと同じ精度の数値解となる
→Kh/Kv=100であれば、水平方向の10倍にできる(未固結堆積物の異方性比は2~100)
→透水係数のアップスケールによる異方性比の表現(スケーリングよりは難透水層を明示的に組み込む方が望ましい)

・透水係数の鉛直方向の変化が重要な場合や鉛直方向の漏水性を扱う場合は3つ以上のレイヤーが必要
・粒子追跡が含まれる場合、上下2層の透水性の大きな差を表現するには3つ以上の遷移レイヤーが必要
・水頭から鉛直流動を計算するのみであればレイヤーの追加は不要
→鉛直透水係数のパラスタによる調整で可

傾斜
・10度以下なら水平として計算可
・10度以上は透水係数テンソル

パラメータの不確実性
・透水係数の幅は12オーダー(特に火成岩、火山岩の亀裂性―非亀裂性)
・キャリブレーションによる改善(9章)、不確実性解析(10章)


2020年9月3日木曜日

TelluSAR

 Tellus に 新しい アドインが実装されました。

TelluSAR
https://www.tellusxdp.com/market/tool_detail/tellus-product/294

干渉画像を作成してくれます。

コンセプトが良いですね。
TellusOS 上でアドインを開くと、見ている場所の SAR データを選んでくれます。一つ選択すると、干渉画像の作成可能な別の日のデータを選出してくれます。これ、地味に面倒でしたのでありがたい機能です。

残念ながら、サーバーが弱く処理の順番待ちが数時間発生していますし、5つの結果までしか保有できません。DLして削除しようとするも動きません。結果も「補正されていない」そうですが、どの補正がなされていないのかは明示されていません。
開発環境は順番待ちで数週間音沙汰ナシ。うーん、使いづらい。(政府の事業が終わったら、これで有料になるのでしょうか?)

ま、それでも今回のアドイン実装や LEVEL1.1 の公開は予想以上でした。
Tellus では以前から SAR データを公開しようとされていましたが、一部地域の画像が公開されていただけでしたし、コチラもこの程度までだろうとあきらめていました。が、先日見た限りでは範囲が広がり、2018~2019年までのアーカイブデータが公開されていました。このまま、数時間前までのデータが揃えば、災害時でも役に立つでしょう(いえ、順番待ちかな)。

今後のデータ公開に期待しましょう。


2020年8月9日日曜日

DualSPHysics 5.0

  DualSPHysics が ver.5.0 に up しました。


今回の注目はコチラ↓

3.17 Rheology models and non-Newtonian formulations

https://github.com/DualSPHysics/DualSPHysics/wiki/3.-SPH-formulation#317-rheology-models-and-non-newtonian-formulations


以前の土砂流下モデルを引っ張り出してテスト中ですが、結果が実現象に近づきつつあることを確認できました。

まだ、式に対応するパラメータ(m, n)の設定をまだ見つけられないなど、計算に使っている式と設定値の対応を完全に把握できていません。

もう少し時間が必要です。


2020年8月7日金曜日

SNAP 7.0

esa の SNAP が Ver.7 にアップされていました。

以下、SNAP v7.0 + ALOS2 (SM, L1.1) での InSAR の手順です。
1- Open product (import ALOS-2 CEOS file)
2- Coregistration
3- Interferogram Formation
4- Topographic Phase Removal (SRTM 1sec)
5- Coherence Estimation
6- Goldstein Phase Filtering
7- Multilook
8- Snaphu export and import (Phase unwrapping)
1- では解凍後のデータを import します。
zip を open すると、以下の通り歯抜け画像になります。
このあたりは v6.0 から変わっていません。

また、読めるのは StripMap (高解像度モード) のみといった報告もあります。私の手元には SM1 しかないので確認できませんが、その可能性はあるでしょうね。
ま、流行の災害対応では L1.5 の Geotiff 利用となっていますので、SNAP は必要ないのですが。

2- 以降は SNAP を使った「処理」です。手元のデータで 6- まで実施しました。
が、干渉しません。
2時期が離れすぎたか?などと思いつつ coregistration の結果を RGB 表示してみました。結果、2時期で尾根が重なっておらず、大きなズレが生じています。

検索してみますと、2-の段階で失敗することが多いようですね。
Radar --> Coregistration --> Coregistration. As the “Initial Offset Method” was set “product geolocation”. 
これに Nearest Neighbor を加えて試したところ、2時期の画像がほぼ一致しました。V6.0 時もこれが原因だったのかもしれません。時間は短縮され17分。マシンスペックの差でしょうか。

最終的に地震時や火山地域のような広域の綺麗な干渉縞は得られませんでしたが、局所的な縞は得られました。まだ手順が足りないのでしょう。

**************************************
20200808追記
同様の状態がココにも。
ENVI SARscape、高価ですが優秀みたいですね。


2020年8月5日水曜日

河床変動計算終了

TELEMAC を一通り使い終わりました。河床変動や下流端の設定で悩みましたが、きちんと動くようになりました。

Fortran での1次元河床変動計算(土石流)のカスタマイズも終了。いくつかのモデルで計算が終わりました。

やっと終わったと思ったら、SNAP と DualSPHysics がアップデートされていました。

うーん。
ひとつづつ確認しましょう。

2020年7月23日木曜日

リニア中央新幹線 有識者会議

先日、御客様との話の中で、静岡のリニア問題に関する話題が出ました。

これ、詳細な議事録が公開されています。以前より興味があり目を通していたのですが、お客様の中にも気になる方がいらっしゃるようですね。

リニア中央新幹線静岡工区 有識者会議

議事録を読んだ際に感心したのが以下の点。(第2回議事録)
  • 解析モデルの不確実性ををふまえた上でのリスク管理システムを作っていただきたいと住民側が要望している。
  • 通常は住民側が安全・危険の2択で質問・要望するのに対し、エンジニア側がリスク管理の視点で説明する。しかし、ここでは逆転している(住民側がリスク管理、リスクコミュニケーションの向上をJR東海さんにお願いしている)。
住民側として説明をうけるはずの静岡県が総監技術士の視点で発言しているのに対し、JR東海さんが主に専門技術士の視点で説明しているような異例の流れとして読めます。JR東海さんの地下水解析の委託先は国際航業さんとのこと。リスク対応方針を決める立場ではない担当技術者の心中を察します。
モデルの不確実性については、「地下水モデル」にも書かれていました(備忘録としてまとめた内容はどこかに置いたかな?)。理解できていない部分がまだ残っていますが、大事な内容です。この不確実性の結果も見据えて静岡県にアドバイスされているブレインがいらっしゃれば、その方は優秀です。見習わねば。


もう一点。使用するアルゴリズムへの委員からの指摘です。(第3回議事録)
  • 委員からは、そもそもの解析ソフトももっと最新のもの、国際的なものがあるので、そちらでやって欲しいとの意見を散々申し上げているが、JR東海からは、高橋の方法だとか国鉄時代のものを使用すると言われてしまうので、そうなると私達としても結果を検証できない。
第1回議事録に「水収支モデル(TOWNBY)」の記載があります。これが「国鉄時代」のアルゴリズムです。TOWNBY 以外の新たな手法を要望される方がいらっしゃるのを初めて知りました。この方の仰る国際的なモノとはどれでしょうね。
広域モデルである程度の精度で予測し、その結果を個別の詳細モデルの境界条件として利用し精度向上を図るという2段階手法は常套手段だと考えます。それを100mメッシュ(感覚では広域モデル相当)のTOWNBY のみで説明を通そうとされるのは、確かに乱暴かもしれません。それ故、副知事の仰る不確実性を明記し、それを踏まえたリスク管理が求められるのでしょう(厳しいですが、理解します)。
私なら、詳細モデルとして何を使うでしょうか?おそらく、知り合いの先生のアルゴリズム(浸透流+地表流NS)を利用させてもらえないか、まず相談するでしょうね。ダメなら GETFLOWS でしょうか。高価ですし、自分で修正できないのが難点ですが。
不確実性は統計手法が良いでしょう。浸透流の予測結果に対して実施したことはありませんが、z分布、t分布などを仮定した検定のような手法を当てはめることができそうです。


目指すは専門技術・リスク管理技術の両輪での解決。総監の視点が必要です。
私レベルだと一人で解決に導くことは難しいでしょう。が、技術者としては腕の見せ所です。JR東海さんが合意に向かう道筋を見つけられるよう、今後を期待しています。


2020年7月19日日曜日

クロネッカーのδ、ディラックのδ関数

紛らわしい名称ですが、内容は異なります。


クロネッカーのδ(単位行列など)

δij=1 (i=j)
    =0 (i≠j)


ディラックのδ関数(SPHカーネルなど)

δ(x-x’)=∞ (x=x’)
   =0 (x≠x’)

Change Detection

航空写真や衛星画像を用いた classification や change detection は古くからあります。

それ故、手法は多数あり、目的や要求精度などによって自由に選択できる状況です。手元の洋書には、change detection のみで8つの手法が紹介されています。

以前に拝見した国総研さんのレポートに、以下の案内がありました。
合成開口レーダ(SAR)画像による土砂災害判読の手引きの作成

関連資料はコチラ↓
2偏波SAR 画像による大規模崩壊及び河道閉塞箇所の判読調査手法(案)
災害時における衛星画像等の活用を促進

2時期の強度画像を用いた目視による崩壊個所の検出作業です。change detection の一つなので、要求される技術レベルは(目検索であれば)低い方です。入社当時、LANDSATデータを画像化してPhotoshopでRGB合成していた、そのような作業と同レベルですね。
昨年の台風19号のように短時間かつマンパワーを必要とする広域災害が今後も発生すると想定されるため、各自治体で分散実施できるように手引きにまとめられたのかもしれません。※今月の SAR は現行体制のままのようです。

残念ながら PSInSAR のような、いくつかの分野に展開できる作業は含まれていません。相対的に高度であることや、現行体制の諸事情?により展開されないのでしょう。で、マンパワーの必要な「作業」のみを外に出す形になるのだと想像しています。

現在、Detection では使用するソースの多様化、高次元化により、コンペ等で機械学習による精度向上が競われているところです。近い将来、マンパワーが機械学習に取って代わるでしょう。民間も売れる技術を整備しなければならないのですが、「SAR って何?」という executives が占める会社では現状、「作業」すら難しい状況です。
ま、ここは他社さんに任せましょう。


2020年7月13日月曜日

極小アレイ + H/V

山本他「岩手県久慈市における微動アレイ探査による地盤振動特性の評価」物理探査 2019 年 72 巻 p. 8-16
極小微動アレイ(3Hz~、深度10m程度の分散曲線)のみでなく、H/V(ピーク周期)も満足する速度構造を推定。


極小アレイは手軽な反面、探査震度が浅いと感じていたのは、私だけではなかったようです。それを補完する目的で(理論)H/V を併用するという考え方は、現場作業が増えず理想的です。
が、理論 H/V の計算手順がわかりません。ミディアムレスポンスって何?状態。

で、プロに聞きました。

ありがたいことに、900ページ超のテキストを御貸し頂きました。また、ソフトも紹介いただきました。

残念ながら、このテキストの導入部(連続体力学部分)以外は全く理解できません。とっかかりすら掴めません。ソフトがあるので実務への適用は困らないのですが、導出くらいは理解したいところです。うーん。学生の頃にもっと勉強しておくべきでした。後悔。

ま、次の機会にプロと相談しながら計算してみましょう。

2020年7月8日水曜日

睡眠

災害は待ってくれません。

今年の災害シーズンは始まったばかりなのですが、いきなり夜半から出勤。

一仕事終えてようやく帰宅。人のためだからできるのかも。

今日は寝ましょう。

************************
20200709
隣の席の方も徹夜でした。2徹にならないようお客様の方が配慮してくださったようです。あれ?

20200713
ブーストが切れるとガクッときます。肉体的にも精神的にも。
年ですね~。


2020年7月6日月曜日

コロナ禍と災害

毎年どこかに必ず訪れる災害。

被災せずとも心が痛みます。

今回の熊本県の災害でも、SAR利用での浸水範囲、土砂崩れ位置が公表されていました。
今日は航測会社さんが、可視画像による判読結果を公開されていました。
もう、一連の流れです。

今年はコロナ禍における災害発生が早くから想定されており、避難所設営時の対応が求められていました。危機管理マニュアルに沿って被災者に負担なく実施できたと願うところです。

2020年7月4日土曜日

TELEMAC-2D + SISYPHE

Finally, riverbed change calculation.
Coupling SISYPHE with TELEMAC-2D (v8p1r1).

Prepare 2 'cas' and 2 'cli' files for sis and t2d respectively. If you set a boundary condition as unsteady sediment discharge, it will be recommended to prepare 2 'liq' files for better interpretability.

The setting that bothered me till the end was 'cli'.
*I've only tried bed-lord, so the following is limited to that. 

**********************************
cli (boundary condition file )
**********************************
First, set the boundary conditions. The user manual (v8p1) contains incorrect expressions . You must enter a non-zero value when setting the time series sediment inflow.
Time-series of sediment discharge
The SISYPHE’s boundary condition file must contain the flags as shown below:
4 5 5 0.0 0.0 0.0 0.0 4 0.0 0.0 0.0 565 1
4 5 5 1.0 0.0 0.0 0.0 4 0.0 0.0 0.0 565 1
Regardless of 1.0 and so on, the time series data in the liq file will take precedence. 
The keyword PRESCRIBED SOLID DISCHARGES must be also included in the steering file,with an arbitrary (but not 0.0) value.

If you set solid discharge =1.0 m3/s (without porosity)  in the 'liq' file and NON COHESIVE BED POROSITY = 0.4, printouts information about solid is 1.66666 m3/s (with porosity) at the time step. In manual,
When the keyword PRESCRIBED SOLID DISCHARGES is used, the mass balance provided in the listing printouts information accounts for the pores = Qb/(1-λ), with λ the porosity.
 
*********************************
liquid file(from the forum)
*********************************
liquid Input file for bed load or suspended load
T QG(1) QG(2) CG(3) CG(1,1) CG(1,2) CG(2,1) CG(2,2) CG(3,1) CG(3,2)
(1) is the number of the boundary

How to IMPOSE a time-series of BED-LOAD discharge ?
I managed to run a simulation with Liquid Boundary File for Sisyphe (version V8P0)
As you did, it's necessary to put QG instead of Q2BOR. Remenber that the unit is m3/s (without porosity).
 
*********************************
rigid bed
*********************************
In manual, 
Non-erodible bed everywhere
Replace in noerod.f the line of code: 
CALL OV(’X=Y+C ’,ZR,ZF,ZF,-100.D0,NPOIN)
with
\texttt{CALL OV(’X=Y+C ’,ZR,ZF,ZF,0.D0,NPOIN)}
 But source shows below
CALL OV('X=Y+C   ', X=ZR, Y=ZF, C=-100.D0, DIM1=NPOIN)
So, I fixed it like ↓
CALL OV('X=Y+C   ', X=ZR, Y=ZF, C=0.D0, DIM1=NPOIN)
I'm checking it now.


2020年7月3日金曜日

TELEMAC parallel

TELEMAC のインストール時に、mpich2 が含まれていました。

「OpenMP ではなくて MPI か」で素通りしました。というのも、MPI の概要は理解していましたが、具体的な利用手順、コマンドを一切知りませんでした。

せっかくなので telemac2D を parallel で動かそうと思ったのですが、そのコマンドがわかりません。マニュアルを見ても具体例が書かれていません。で、Forum を検索。発見。

まずは初回にユーザーと pass を登録。
mpiexec -register

次に実行。
CONFIGNAME '-c wing64mpi' を加えただけです。CONFIGNAME って、config ファイルに書かれているコンパイラー識別コードのようなモノのことだったのですね。 
telemac2d.py -c wing64mpi ***.cas

cas ファイル内で4コア指定
PARALLEL PROCESSORS = 4

これで parallel 版も動きました。single で14.5時間かかった自由表面流れのモデルが、4コアで9.5時間でした。早くなっています(といっても、2次元計算にしては遅いと感じています)。

Forum を見る限り、領域分割を行っているようです。今、詳細を追いかけようとは思っていませんが、いずれ理解したいですね。

2020年7月2日木曜日

TELEMAC keywords

TELEMAC のマニュアルが薄くてラッキーと感じていたのは最初だけでした。

離散化は FEM or FVM ですが、まず、その切り替えにかかわるキーワードの設定が書かれていません。Forumで検索して見つけたのですが、このような説明不足と感じるキーワードが進むにつれて多く出てきました。都度、調べて修正し、計算を実施。

設定を進めていくうちに、エラーが減り、計算速度が向上し始めました。最終的には、初期モデルに河床変動計算を加えた非定常計算が 5倍速となりました。恐ろしい。
(time series of sediment discharge の設定は時間がかかりました。マニュアルが誤っており、Forum にも具体的な設定方法が上がっていません。これについては後日、まとめておきましょう。)

以下、備忘録的に控えていたリンクや説明を記載しておきます。

*********************************
SVN revision
*********************************
CollabNet Subversion Repository: opentelemac - Revision 15180: /tags
..
v6p0r8/
v6p1/
v6p2r1/
v6p3r0/
v6p3r1/
v6p3r2/
v7p0r0/
v7p0r1/
v7p1r0/
v7p1r1/
v7p2r0/
v7p2r1/
v7p2r2/
v7p2r3/
v7p3r0/
v7p3r1/
v8p0r0/
v8p0r1/
v8p0r2/
v8p0r3/
v8p1r0/
v8p1r1/

*********************************
データ、ソフト
*********************************
TELEMAC-2D tutorial 

Graphical interface

Blue Kenue

Fudaa

SALOME

*********************************
Blue Kenue v3.3.4 操作
*********************************
メッシュ(slf)ファイルの地形標高を指定するキーワードは大文字で「BOTTOM」
×Bottom・・・認識せず、標高0mで計算されます。
〇BOTTOM

アニメーション表示
左のツリーから見たい結果を右クリックして Animate にチェック。

TXT書き出し
左のツリーから DEPTH 等を選択。
書き出したい時間を表示した後に tool bar から file - save copy as で XYZ ファイル保存

断面
open line  作成 → 右クリックから resample → tools - map object で結果を割り当て→ tool bar から Map Object で結果を当てはめ。 → file - save copy as で distance (xy) 保存。

水と土砂の境界条件設定位置、数を合わせること。
土砂流入がなくても、流入0で位置・数を合わせないと、計算時にエラーが発生。

*********************************
keywords(Forum からの引用)
*********************************
COMPATIBLE COMPUTATION OF FLUXES

FREE SURFACE GRADIENT COMPATIBILITY
Theoritically, this option is relevant with both options (1: primitive equations and 2: wave equation). The key idea behind, is to break the rigorous compatibility between water depth and velocity components which verify continuity equation in order to damp spurious oscillations that can appear in areas with high gradient of topography.
For Telemac-2D, this option is only implemented with wave equation. This is logical while this option is more fast and stable than option 1 and thus is oscillations appears with this option they will appear and will grow more and more with option 1.

MASS-LUMPING ON H
TREATMENT OF NEGATIVE DEPTH
IMPLICATION FOR DEPTH
IMPLICATION FOR VELOCITIES
MASS-LUMPING ON H : this changes the matrix stemming from the derivative in time of the depth into a diagonal. In some sense moves from finite elements to finite volumes. Effects : stabilises, smoothes. In steady state just stabilises, no effect on results.

TREATMENT OF NEGATIVE DEPTHS : the fastest algorithms in Telemac allow some slightly negative depths. To get strictly positive depths TREATMENT OF NEGATIVE DEPTHS = 2 recomputes the continuity equation and yields an equation solved with machine accuracy and positive depths. 
IMPLICITATION FOR DEPTH (or VELOCITY) : takes depth or velocity explicitely (time n) or implicitely (time n+1) in various terms of the equations. Implicit is more stable but smoothes, semi-implicit (0.5) is more accurate but the limit for stability. Important only for waves, has no effect on stedy state flows.

×FRICTION ANGLE OF SEDIMENT・・・マニュアルの誤記
FRICTION ANGLE OF THE SEDIMENT

×VARIABLES FOR GRAPHIC PRINTOUTS ='U,V,S,H,B,E,M,R,QS*'
VARIABLES FOR GRAPHIC PRINTOUTS =U,V,S,H,B,E,M,R,QS*


*********************************
FEM, FVM(Forum からの引用)
*********************************
how to choose the "finite volume scheme"
EQUATIONS = 'SAINT-VENANT VF'

**********************************
restart
**********************************
COMPUTATION CONTINUED        = YES
PREVIOUS COMPUTATION FILE    ='********.slf'
INITIAL TIME SET TO ZERO         = NO

1時間毎に記録し、1日後(24番目)の結果を利用したい場合
RECORD NUMBER FOR RESTART  = 24



**********************************
20200702修正
モデル内の計算時間を短くしていたのを忘れていました。
設定しすぎると、逆にタイムステップを乗り越えられない場合が出てきました。まだまだ触らないとわからないですね。

20200703追記
リスタート機能に関するキーワードを追加。
SISYPHE のkeywords を移動。↓
parallel はコチラにまとめました↓