やり残し事項 2019

今年は研究職を半分担当しました。

内業が主体となり、時間的にも体力的にも楽になりました。報酬の対象となる勤務時間だけを見ると確実に短縮されています。働き方改革と銘打って時短を推奨されている御時世ですが、素直な後輩たちの動向を見ていると将来への危機感を抱かざるを得ません。

人のことを言えないのですが、やり残し事項は2年連続ほぼそのまま。
かわりではないですが、3Dレーザー関連を扱えるようになりました。python を利用したデータ分析、機械学習関連のテクニックもさらに身に付けることができました。


以下、やり残し事項です。昨年とほぼ同内容ですが、一応書いておきましょう。

道具
帯磁率計
ガンマ線測定器（携帯型）

技術
動的解析（耐震、液状化）
斜面設計

コード・ソフト
優先度低：DtransuのCUDA化・・・100万で購入可。
優先度低：PSIｎSAR・・・700万で購入可
優先度中：地表流＋地下水・・・GSFLOW、HydroGeoSphere
優先度中：GMS・・・J-SHISの取り込みまではできました。いつか続きを。
優先度中：OpenMI・・・PHREEQCとDtransuの連成等


中期目標3年目はもがき中。あと2年。間に合うかな。

的中率、見逃し率

砂防分野では「的中率」や「見逃し率」などの用語をよく聞きます。

独自に進化したのか、機械学習分野と同じ定義でも異なる用語が用いられています。毎回、変換に戸惑います。

先日、以下の論文に砂防分野での定義の記載がありました。書き留めておきましょう。
堤ほか「ストリームチューブによる地形分割を基にした表層崩壊解析手法」2019, 砂防学会誌, 72, 2, 343, 3-13

TP:A
FN:B
FP:C
TN:D

的中率=(A+D)/（A+B+C+D）
捕捉率=A/(A+B)
空振り率=C/(C+D)
見逃し率=B/(A+B)

正答率：accuracy rate = 的中率
再現率:recallのうち、感度:sensitivity = 捕捉率
　　　　　　　　　特異度：pecificity = 空振り率
適合率：precisionの利用で、1-陰性的中率=見逃し率

imbalanced data のままだと、accuracy は機能しません。個々の評価指標を押さえておく必要があるのですが、重視すべきもののみ砂防分野で再定義されているということでしょう。

***********************************
20200304追記
「的中率」よりは「適中率」をよく見ますね（私は後者を使っています）。
また、上記の文献と気象庁では、一部の用語の定義が異なっています。ややこしや。
https://www.data.jma.go.jp/fcd/yoho/kensho/explanation.html

		予 　　報
		降水あり	降水なし	計
実況	降水あり	A	B	N1 = (A+B)
	降水なし	C	D	N2 = (C+D)
	計	M1 = (A+C)	M2 = (B+D)	N = (A+B+C+D)

• 適中率（降水の有無の適中率）（％）＝（A＋D）÷N×100
• 「降水あり」予報の適中率（％）　　＝A÷M1×100　（一致率に同じ）
• 「降水なし」予報の適中率（％）　　＝D÷M2×100
• 見逃し率（％）　　　　　　　　　　＝B÷N×100 
• 空振り率（％）　　　　　　　　　　＝C÷N×100
• 捕捉率（％）　　　　　　　　　　　＝A÷N1×100
• 一致率（％）　　　　　　　　　　　＝A÷M1×100　（「降水あり」予報の適中率に同じ）

機械学習モデルの解釈可能性

機械学習モデルの解釈可能性に関する、素晴らしいスライドと kaggle マイクロコースです。スライドの１、マイクロコースの２が昨日の特徴量評価にあたります。

スライド
https://speakerdeck.com/line_developers/machine-learning-and-interpretability

1. どの特徴量が重要か
　Feature Importance
2. 各特徴量が予測にどう影響するか
　Partial Dependence
　Surrogate Model
3. ある予測結果に対して特徴量がどう寄与するか
　LIME, SHAP
　Grad-CAM, Grad-CAM++

kaggle
https://www.kaggle.com/learn/machine-learning-explainability

1. Use Cases for Model Insights
　Why and when do you need insights?
2. Permutation Importance
　What features does your model think are important?
3.Partial Plots
　How does each feature affect your predictions?
4.SHAP Values
　Understand individual predictions
5.Advanced Uses of SHAP Values
　Aggregate SHAP values for even more detailed model insights


両者に出てくるSHAP、発想は面白いですね。単純な数式でモデルの推定結果を模擬し、各特徴量にかかる係数で定量的に結果への寄与度を表現しています。これは決定木ベース以外の説明できなかったモデルに対する解釈性の要望を反映したものでしょう。

が、正解を模擬しないのはなぜでしょう。あくまでモデルの解釈性を加えることが目的だからでしょうか。それとも理解できていないのでしょうか？
精度を犠牲にしつつ解釈性を加えた簡易モデルの追加。そこまでして「ブラックボックス」を避ける必要はないと思います。最後に判断するのはヒトなので。

このあたり、機械学習分野の今後の動向が注目されます。


********************************
20200104追記
技術評論社「Kaggleで勝つデータ分析の技術」
6.2特徴量選択および特徴量の重要度
6.2.2 特徴量の重要度を用いる方法
・モデルから出力される重要度：RF、GBDT
・他の方法：Permutation Importance、null importance、boruta、特徴量の機械的大量生成

特徴量選択

機械学習の特徴量選択について。

機械学習にかける前に有効な特徴量を選択するのですが、その手法は以下のように区分されています。https://www.oreilly.co.jp/books/9784873117980/

1. filter method、単変量統計: 特徴量個々の分離性・相関係数・分散などを利用。
2. wrapper method、反復特徴量選択: 再帰的特徴量削減(Recursive Feature Elimination, RFE)、Boruta など。feature_importance 利用。計算コスト高。

決定木ベースの手法による feature_importance を利用すれば、モデル作成と重要度の評価が同時に実施できます。３つ目の手法としてモデルベース特徴量選択と呼ぶそうです。

分類問題では、決定木ベースのモデルに依存した選択手法が主体でしょう。ライブラリも提供されており、環境は整っています。

逆に、決定木ベース以外のモデルでは feature_importance 属性を利用できませんので、RFE はダメ、モデルベース特徴量選択も当然ダメです（XGBoost や LightGBM 等のモデルが分類問題に対し高精度に機能することを示してきた中で、あえて他の手法に積極的に移行する動機はないのですが）。

決定木ベース以外では、モデルに依存しない指標を利用するようですね。例えば、Permutation importance（組もうかと思ったら、wrapper がありました）。

これは素人でも容易に思いつく手法なのですが、高コスト。完走させるために軽い手法を選択しても結果が良くないと、その中での指標の低下量、重要度を比較するだけになります（木を見て森を見ず状態）。

いずれにしても、何らかの評価・解釈を行って重要度の低い特徴量を破棄し、モデル全体の計算コスト低下、予測精度向上に努める必要があります。

特徴量の評価・選択は前処理に続き必要な作業だと、あらためて認識させられました。

********************************
20200104追記
技術評論社「Kaggleで勝つデータ分析の技術」でも同じ3区分でした。
6.2特徴量選択および特徴量の重要度
6.2.1 単変量統計を用いる方法
6.2.2 特徴量の重要度を用いる方法
6.2.3 反復して探索する方法

手法に関しては以下の記載がありました。
「なお、scikit-learn のドキュメントにはモデルの選び方として図4.8がありますが、データが十分にあるテーブルデータの回帰や分類のタスクでは、GBDT でうまくいくことが多く、このように細かに条件で分岐させる必要はあまりないでしょう。」p230

SVM

延長戦突入。

今日はサポートベクターマシン（support vector machine, SVM）。SVM はデータ数が少ない場合に用いられる手法で、昔よく利用されていました。が、今回は大量のデータを処理する必要があります。力技でこなそうと、scikit-Learn で組んでみました。が、やはり遅い。

SVCではなく、SVRだとどうか？と試しましたが、終わりません。

GPU版がないかな？と思い RAPIDS の中を覗いてみると、cuML の中に SVC がありました。が、fit()を使うと GTX1060 3GB のメモリに載りません（想定内）。scikit-Learn の cross validation に放り投げてみるも連携しません。

デモデータだと計算速度が40倍も違ったので、GPU は利用したいところです。

RTX8000 48GB なら動くでしょう。汎用性はなくなりますが、このマシンのみで動かすことを前提に、組んでみましょうか。

********************************
20200104追記
結局、他のライブラリとの連携の関係で、scikit-Learn を選択。年末から調整を続けていますが、最終モデルは40時間以上経過した今でもまだ完走しません。遅すぎ。特徴量を減らさないとダメということはわかりました。
技術評論社「Kaggleで勝つデータ分析の技術」では以下の記載がありました。SVM の使われない理由が明記された書き物を見たのは初めてです。

本章では、テーブルデータを扱う分析コンペにおいて使われるモデルとして以下を紹介します。・・・・・・ なお、教師あり学習のモデルとして良く紹介されるモデルのうち、サポートベクターマシンは精度や計算速度が見劣りするため、あまり使われることはありません。

********************************
20200109追記
O'REILLY「Pythonで始める機械学習」にも、以下のような記述がありました。

• サンプル数が大きくなるとうまく機能しない（10万サンプルぐらいまではうまく機能する）
• 10万サンプルぐらいになると、実行時間やメモリ使用量の面で難しくなる。
• 注意深い前処理とパラメータ調整が必要。→決定木ベースが良く用いられる理由の一つ。
• 検証が難しい。予測された理由を理解・説明することが困難。

 Ｃ、gamma の影響についてはｐ98に図解あり。

不飽和透水試験

不飽和透水試験機のデモに立ち会いました。

試験機自体は古くからあるのですが、小型化、高速化したと宣伝されてた機種です。小型化はまだしも、機器の改良だけで高速化？という点が気になっていました。

実際に見せていただくと、測定自体は確かに速い。が、不飽和透水試験で時間がかかるのは測定前の浸透。今回は肝心のプレ浸透が済まされており、試験だけの営業でした。
何時間浸透させたのかと聞いてみると、ほとんどしていないとのこと。つまり、測定が速い=高速化と宣伝されているようでした。うーん。

ま、プレ浸透を長くすれば飽和透水係数を得られるわけでもありません。が、雑に扱うのは論外。短時間のプレ浸透で飽和とみなせる理屈はないでしょう。
https://phreeqc.blogspot.com/2019/09/blog-post_17.html

残念ながら私の勘違いだったのですが、コンセプトは良いと思われます。これからも開発が続くことでしょう。一つぐらいは手元に置いておきたい機材でした。

正弦波

正弦波
=================================

振幅u0、角速度ω のとき、

u(t)=u0sin(ωt)

これはx=0での変位。
x=x位置での変位は、ちょっとしたアイデアで求まる。

x=0からx=xまで波が到達するのにかかる時間はx/c秒。
(x/c秒前の原点の状態を求めれば良い）

u(t,x)=u0sin(ω(t-x/c))

T=2πr/rω=2π/ω
λ=c/f
f=1/T
k=2π/λ　より

u(t,x)=u0sin(ωt-2πx/Tfλ)
　     =u0sin(ωt-2πx/λ)
　     =u0sin(ωt-kx)

RIVeR

PIV、STIV で河川の表面流速を比較しました。

以前、Python で組んでいたPIV＋STIV。これに動画をよませるだよませるだけです。
https://phreeqc.blogspot.com/2019/10/piv-stiv.html
https://phreeqc.blogspot.com/2019/10/stiv.html

今回は流速を実際に計測していますので、答え合わせができます。
PIV は静かな流れだと流速を掴みづらいようです。パラメータを調整しましたが、目立つ浮遊物がない限り、実際よりも低速と判定されます（というよりも認識しません）。
STIVは安定しています。浮遊物がなくても安定して流速を把握できます。
両者とも正しい答えを出すのですが、安定性で STIV 優位でしょうか。結果的には土研さんや国交省さんと同じ選択となりました。
http://www.jsece.or.jp/event/conf/abstract/2014/pdf/P2-56.pdf

表面流速が得られても、流量を出すにはひと手間必要です。複数のラインで平均流速を出したのち、断面形状と位置をあわせてから、それぞれの流速が受け持つ範囲を決めて流量を出す必要があります。このあたり、1から組むのは面倒。
ソフトがないかな？と探してみると、ありました。またしてもUSGSが噛んでいるのでしょうか。

RIVeR
https://www.sciencedirect.com/science/article/pii/S0098300417307045?via%3Dihub#!
http://riverdischarge.blogspot.com/

使ってみたところ、優秀。計算はもちろん、静止画への変換から幾何補正まで一通りの機能が備わっています。Unshake と言って、画像のブレを補正してくれる機能まであるので、スマホ撮影でも対応できそうです。しかも、PIV 版と STIV 版の両方が備わっています。
残念ながら、PIV で一部の流速を得られないのは同じ。これはソフトに起因するものではなく、手法に起因するものでしょう。一方、STIV では安定した結果を得られます。良いですね。

PIV で指定横断に直行方向の流速が得られる点もgood。Python では u,v 成分から指定横断位置での直交方向の流速を出す必要がありました。これを自動で算出してくれます。

残る課題は高さの把握。
洪水時等に幾何補正を為せるだけの高さ情報が必要になります。これを効率的に得るためにはどのようにすればよいかアイデアが必要があるでしょう。この点では、水路のような形状だと楽なのですが。

ひとまず、理屈を理解し、ソフトを準備しました。
残る課題をクリアし、本番に備えましょう。

ネットワークカメラ

ここ２か月で、３台のネットワークカメラを設置。

カメラの世界、進んでいますね。

http のみならず様々なプロトコルに対応しています。容易にデータを入手・転載できるほか、警報メールなども出せるようになっています。当然、発報にはトリガーが必要で、音の異常、動きの検出などを判定するためのアルゴリズムが入っています。精度は別として、技術レベルとしては高くないですからね。
https://phreeqc.blogspot.com/2019/10/blog-post_27.html

ハード面では、夜も見えるように近赤外領域を撮れるカメラも普及しているようです。併設するライトまで「近赤外LED」です。これは、自動車にも使われていますので、カメラの世界では浸透しているのでしょう。
高解像度化や省電力化も進んでいるようで、ライトを合わせて数Wの消費で済むように工夫されているものもありました。

知らない間に、ハードやソフトが変化しています。
さらに未来の状況は容易に想像できます。それを考えると、ここで脱落するわけには参りません。
発展させる立場にないことはやや残念ですが、開発者に敬意をもちつつ、積極的に利用させていただきましょう。

Fourier Transform in Python

先週末から触っていたフーリエ変換。結局、Python を使いました。

有限複素フーリエ係数Ck。Nで割る前の値が出ます。

Ck=np.fft.rfft(data)/N

周波数。ここまではいつも通り。

freq=np.fft.rfftfreq(N,d=1./fs)

躓いたのがランニングスペクトル（スペクトログラムという方が正解でしょうか）。scipyでの計算結果が EXCEL の計算結果と一致せず、最初は何の値が吐き出されているのかわかりませんでした。

from scipy import signal

freqs, length, Sx = signal.spectrogram(data, 

                                       fs=fs, 

                                       nperseg=segment,

                                       noverlap=segment-1,

                                   #raw-data

                                       window=('tukey', 0.), 

                                   # Remove linear trend along axis from data.

                                        detrend = False,

                                  # power spectral density(V**2/Hz)

                                  # or power spectrum (V**2) 

　　　　　　　　　　    scaling='spectrum', 

　　　　　　　　　   #mode='magnitude'
                                    )

ｗindowを削り、傾き補正を削り、スペクトルを指定すると、結果が一致。デフォルトではこれらが指定されていたため、明記しないと外れませんでした。

得られる値は Power。継続時間 T=ndt は乗じられていません。

この T や N が乗じられているか否かは式を見ないとわからないのですが、明記されていない場合があります（今回もそうでした）。

プロに聞けば、それぞれの流儀があるそうです。地震分野でも「私は乗じない」と言われる方がいらっしゃるそうで、決まりはないとのこと。ま、Nの場合は逆算する際に気を付けておけば良いだけですので、あまり気にする必要はないのでしょう。

EXCELでのミスも見つかりました。Power の最初と最後は2倍しない、など。数式を見直すと、確かにそのようになっていました。

初心者は手を動かさないと正解にたどり着けません。EXCEL で整理し、Python で答え合わせをしました。これで次は大丈夫でしょう。

Fourier Transform in EXCEL

この週末、EXCEL でのフーリエ変換を試していました。

これまで FFT は Python を使用しており、仕事でも Python を使うつもりでいたのですが、今回は EXCEL。理由は「後輩君が EXCEL で始めたから」。

私は空間、彼は時間を扱っており、別の作業です。が、DFT 部分は共通。彼が EXCEL で作り始めた（が正解にたどり着けていない）ので、今回はそれをフォローできる体制を整えておくことにしました。

EXCELでは、以下の3つの方法が考えられます。
1．分析ツールの「フーリエ解析」を使用する。
2．関数を使用する。
3．VBAで組む。

まず、分析ツールのフーリエ解析。試したところ、一番手軽でした。
当然、後輩君はこれを使用しています。が、 FFT なのでデータ数は2の累乗に限定されます。たとえパディングでかわしたとしても、採用されている計算式が明記されていないため、結果をデータ数で割るべきかどうかわかりません。故に振幅も出せません。
（VBA での結果と比較すると、データ数で割っていない値が表示されていました。）

関数では複素数を扱うのがやや面倒です。四則演算のみでもIM○○といった関数を使う必要に迫られます。さらに、k次の答えを求めるためには、k列必要です。これでは非現実的。

VBA では複素有限フーリエ変換を組もうとしました。が、一旦躓きました。VBAでも複素数を扱いづらい。WorksheetFunction でもエラーが出るので原因究明しかけました。が、結論としてはその時間をとらずに有限フーリエ変換に変更しました。今回はデータ数が700以下と少ないので、FFTでない方が効率的なのかもしれません。
ここまで来れば簡単。k、mの2重ループで有限フーリエ係数を出して、それから振幅、フーリエ振幅、Power、複素振幅を計算します。

計算後、スペクトルを図化する際にどの振幅を使うか迷いました。
調べてみると、フーリエ振幅ですね。今まで、スペクトル表示の際はただの振幅を利用しているものと思っていましたが、基本はフーリエ振幅を使っているようです（ランニングスペクトルではPowerが一般的）。慣習のようですね。分野にもよるのでしょう。


EXCEL では VBA が楽でした。分析ツールの出す答えが何を示すのかも分かったので、今後はコチラも使えます。
彼が週明けに自力で正解にたどり着いていたなら、それに越したことはありません。どうなっているでしょうか。
ま、頭の整理ができましたので、手を動かして良かったと思います。

*******************
20191217
結局、私は Python を利用しました。これが楽でした。

TREND-POINT

地上レーザーのデータを可視化しています。

サーフェス化して断面を切っていますが、同時に航空LPのデータも切っています。大体、同じ形状になります。両者の精度、思ったより良いですね。

レーザーのデータから地表を抽出するには、木や草を除去する必要があります。また、座標の変換も必要な場合があるでしょう。これ、ReCap ではできません。

プロに聞いたら、ライセンスに空きのあるソフトがありました。
TREND-POINT
https://const.fukuicompu.co.jp/products/trendpoint/

Leica のデータも読めました。
が、樹木等を除去するのはかなり困難。プロもトライ＆エラーで解決するそうです。航測会社さんの処理はとても速いのですが、きっと独自のノウハウがあるのでしょう。
これはあきらめました。

座標変換は簡単。
図面から公共座標を拾って点群に与えるだけ。それを計算させると独立座標から変換してくれます。この機能だけでも使えそうです。

プロに聞いたところ、点群を平面図に変換できるソフトはないようです。が、いずれできるでしょうね。フルサーフェス化したり、ソリッド化したり。

新たなハードが出ると新たなデータが生まれ、それを効率よく処理するソフトが整備されます。昔からそうですが、近年は浸透が速く、多様化しているように思えます。
頑張って、ついて行きましょう。

RTC 360

再度、3次元地上レーザースキャナー。

ミスを繰り返さないよう、準備期間を設けて実戦！
というほど気負う必要もないですね。注意点を頭に入れておけば、簡単です。

今回は 前回と同じ Leica RTC360 を使用。

• 距離：中密度で 130m 飛びます。現場では、概ね 100m 程度飛んでいるのは確認できました。
• 座標：前回、ターゲットを入れないと座標付けに手間取ることが分かったため、全箇所ターゲット設置 ＋ VRS 測定。

取得した点群は iPad で合成し目視確認。今のところ、順調です。

明日も落ち着いて作業しましょう。

スマホで GIS

GIS データを現場で確認する手段を探しました。

複数のデータを重ねているのですが、それらを個別に、現場毎に事前に印刷するのは面倒です。online が当たり前の時代、データを クラウド上に up すればスマホで確認くらいはできるでしょう、と調べ始めました。

まずは QGIS。
これはオフライン&Android しか方法がありませんでした。データをスマホにコピーし、アプリ「QField for QGIS」で閲覧するのみ。

次に ArcGIS。
Pro でも Map でも良いのですが、オンラインでグループ内共有し、アプリ「Collector for ArcGIS」で閲覧します。が、非常に複雑な手順です。更新した場合も手間。サポートさんに聞いても時間がかかるばかりで、こちらは現場に間に合いませんでした。

どちらが良いか？というと、QGISでしょうか。
今回は携帯電波の入らない現場が複数あり、オフラインでデータを閲覧できるQField for QGIS に助けられました（地理院タイルはダメでしたが）。バッテリーは1日持続するようになりましたが、エリア外は広く残っています。オフラインありきで準備しないといけないですね。


********************************
20191201追記
地理院タイルはキャッシュを使えました。
キャッシュ＋立ち入り禁止区域＋スマホGPS で位置関係が分かりました。便利。

プロの仕事？ その2

地上レーザースキャナーを初めて使いました。

機材所有部署から事前にレクチャーを受け準備していたのですが、事情により現場の前日に別機種に変更。再度1時間ほどレクチャーを受けました。
テスト中、手元では70m程の点群しか表示されていなかったのですが、射程は倍だそうで、映っていない点群も解析時には使用できるとのことでした。

で、現場。
今まで3次元地質モデルを作成する際には、地表面を作るのがネックになっていました。
測量図面から3次元に起こすのは手間。UAV ＋ SfM が楽なのですが、自分では積極的に操作したくはない。結果的に他部署にお願いする形になっていました。
ですが、このスキャナーであれば墜落の心配はなく、数分でデータを取得してくれます。何度か場所を変える必要は生じますが、GPSも付いており、簡単・楽です。

帰社後にデータを確認。
が、70m以遠のデータが出てきません。スキャナーを所有していた部署にも確認したのですが、「そのようなことはないはず」との回答。弱りました。

結論から言うと、データは取れていませんでした。いくつかのモードの中で、点群を最も密に取得できるモードを選択していたのですが、この場合は65mまでと’マニュアル’に記載がありました。そう、マニュアルに。
読まずに出かけました。反省。

で、現場はやり直し。
これがプロの仕事？と問われると恥ずかしくなるような内容です。ベテランですし、変化に追従する方法をよく知っていましたが、実施しておらず失敗。お粗末です。

次は失敗しません。十分に準備していきましょう。

プロの仕事？

ドライブレコーダーの電圧感知機能がイマイチ。

購入したカー用品店にて電圧を計って不具合を説明をしたところ、本体交換となりました。が、交換しても直らない。で、2回目は配線ごと交換となりました。配線は自分で行っていたのですが、無償で取外し＆取付け＆今からできる、ということでしたので若いプロにお任せしました。

が、いつまで待っても作業中。それほど時間がかかる作業ではないのにどうしたのだろうと思いながらピットに行ってみると、「作業は終わったがドアロックができなくなった」とのこと。面白い。なぜドラレコを外しただけでロックができなくなるのかわかりません。

3人がかりで2時間ほど調べていただき、やっと原因がわかりました。室内イルミのヒューズが飛んでいました。サービスマニュアルを見ない限り、イルミとロックが連動しているとは容易に思いつかないでしょう（ま、ロックしたらイルミが切れるので逆から考えると繋がりますが）。

なぜ飛んだのだろうか？と帰ってきて配線部分を確認してみると、わかりました。取付の際にカプラー手前の細線からエレクトロタップで常時電源を分岐させていたのですが、その先のギボシや2股部分（変更・拡張の可能性を考えていました）で配線を外さずにタップごと、根元からすべて外されていました。これ、アースを取っていたフレーム近くにあるので、触れさせると一発でショートします。
調べてみるとこのカプラー、フットランプ用でした（怪我の功名。翌日、フットランプを付けました）。

エレクトロタップ外してしまったら次の取付けはどうするのでしょう（ヒューズボックスでACCを取れないのですが）。しかも、外すときに本体側の線を切断するリスクもあったろうに。というか、作業前にバッテリー外してなかったんですね。

若いプロがミスをして、ベテランのプロがフォローする。どこかで見たような構図です。が、アマがわかるようなケアレスミスを起こさないよう、事前に教育しておいていただきたいところです。

これ、ブーメランで帰ってきます。

地すべりと地質

若い方に、興味のある研究テーマを聞いてみました。

「加水ハロイサイト」をキーワードに挙げられていました。近年の地震災害でマスコミに取り上げられた鉱物です。

以前は層状粘土鉱物が注目されていました。それらが整理・オーソライズされるに至らず、別の鉱物へ新展開を見出すようになった感がありますね。長い間この業界にいたにもかかわらず、何の成果も出せなかったことには、責任の一端を感じます。

地すべり災害は地質（岩種）に関係するといわます（がけ崩れ災や土石流災も岩種に関係しているといわれますが、これは疑問です）。
新鮮でも軟らかい岩種、（先天的、後天的にかかわらず）亀裂の発達した岩種に特徴的なのは感覚で理解できます。それらの粘土鉱物の力学的性質に着目するのも自然でしょう。
が、地質屋さんとしてはここまで。鉱物と力学的考察について統一見解が出されるまでには至りませんでした。

より難しくしているのは、地形的・水理的な力学バランスと構成鉱物の力学的性質が重なった結果を、一方だけで解釈しようとしているからでしょう。いえ、着手時は２属性のうち1属性に着目するとしておきながら、いつの間にかそれだけで binary classification として着地させたいロジックのすり替え、焦りがあったのかもしれません。

ま、一方だけでもデータが整理されるのは、後進のためにはありがたいことです。が、現段階では、統計や機械学習の視点ではデータが足りません（この程度の問題では機械学習に頼るまでもないでしょうが）。

データが揃えば、自ずと解けるレベルの問題でしょう。
ハロイサイトでも良いのですが、それだけに着目せず、力学的条件もあわせて整理しておけば、いつか後進が解いてくれることでしょう。

BASEGRAIN

BASEGRAIN
https://basement.ethz.ch/download/tools/basegrain.html

ダム堤体のリップラップ材や河床材料の調査に使えるソフトです。
一通り動かしてみました。動作のおかしなところもありましたが、ほぼ問題ないですね。便利です。

ソフトの存在を知らなかったのはもちろん、このようなニーズがあるとも知りませんでした。UAVを利用して大量のデータ取得、そのデータを利用した粒径解析が行われているようです。

ただ表面の粒径だけなので、河床変動シミュ用の河床材料調査には向いてないのでしょうか？いえ、広域での傾向がわかるので検証等に向いているのでしょうか？わかりません。
学ばねば。

近年、いろいろな分野で広域かつ大量のデータを取得・解析することができるようになりました。扱う側も意識改革を迫られています。
変化に追随するのは大変です。が、そこも技術者として楽しむべきところでしょう。

QGIS Desktop 3.8.3 with GRASS 7.6.1

ArcGIS の入っている PC を他の方が利用。

QGIS があるから良いか、と思いながら作業開始。
QGIS では GRASS の機能を使用できます。SAGA も入っており組み合わせると概ね主要機能が揃うというのはありがたいですね。

が、今回はハマりました。
Ver3.8では、GRASS が外されて起動するのがデフォのようです。いえ、GRASS のメニューは出てくるのですが、使用するとエラーが発生します。
どうも、起動時に「QGIS Desktop 3.8.3 with GRASS 7.6.1」を選択しないとダメなようです。3.4の頃はデフォで動いていたのですが。
https://gis.stackexchange.com/questions/313724/grass-not-working-with-qgis-3-6

安定性やサポートでは Arc に勝てません。が、無償利用できるのはありがたい。感謝ですね。

動態検知

千葉県での水害・災害が今秋3回目。

千葉に限りませんが、人命を失う惨事を繰り返さないように、早期に逃げてもらうきっかけを作れないか？と思案中。その過程で、ここ数日は動画の中から「異常」をピックアップ・周知する手法を模索中でした。

まず思いつくのは異常検知。
時間軸で主成分ベクトルを求める場合のように、全期間通して学習させる必要はありません。計算コストが小さくて済みます。でも、オンタイムでできるか不安。

どうしましょう？などと考えながら風呂に入っていると、簡単な方法を思いつきました。動態検知です。いえ、手法は前から知っていましたが、適用性をまったく考えていませんでした。人や車など、自動運転への適用に頭が向かっていましたので。

これ、OpenCV だと簡単。example に少し手を加えると完成するレベルです。しかも公開されているライブカメラでもラップトップで処理が追いつく程高速。休みの間、データを取り続けましたが難なく動きました。使えます。
動態検知した面積を積分し時系列に直して移動平均をとると、トリガーを設定できそうでした。こうなると異常検知メールも出せるようにしよう、などと手を出しかけましたが、解除設定等が面倒になったため、テスト送信のみで終了。今回はココまでとしました。

この手法、課題も多くあります。単純な前データとの差分ではダメで、加重平均をとって係数（半減期のようなもの）を調整したり、動きとみなす閾値を決めたりする必要があります。しかも、カメラ設置場所毎に調整が必要です。まだ実用段階には程遠い（広域でなければ簡単かも）。

この動態検知、OpenCV が動けばよいので、ラズパイでもできてしまいます。というか、ドライブレコーダーの駐車監視がこの系統なのでしょう。いくつかの機種には動態検知入っていますので（OpenCVとは限りませんが）。
そう考えると、すでに社会に広く普及している技術であり、土木分野が遅れているだけ、気づかなかっただけ、ということになります（もっと柔らかい頭が欲しいですね）。

他分野のプロが介入すれば、土木分野の多くの課題がすぐに解決されるような気がします。自動運転のように市場は大きくないので、なかなか踏み込んではいただけないでしょう。が、人命がかかっています。CSR がてら興味を持たれる方はいませんかね。

無知は罪 2

計算用に3台集めたPC。

全てメモリを 64GB にアップ。内蔵 HDD として 4TB 追加。これでギリのスペック。

それを聞いた他支店の技術者？が、手元の PCを「64GBにしたい（してほしい）」と言い出しました。よくよく聞いてみると、必要な作業にどのようなスペックが必要か理解せずに PC を買ったようでした。新しい PC = 速い の理屈？らしいです。実装されているメモリやコア数も御存知ないようでした。

調べてみると、実装 16GB。
計算ではメモリを 8GB しか使ってません。使いたい Arc は 32bit 版。メモリを増やすよりは 64bit 版の Arc に移行して、推奨GPUを付けた方が効果あるでしょう（それを調べて PC を発注すべきですが）。
一方、今後予定している計算に必要なメモリをプロに確認してもらったら、まさかの64GB越えでした。マザーは 128 まで実装できるようなので新 PC も無駄にはならないですが、複数台 PC を確保しておいた方が現実的だったでしょう。

これ、あるあるです。
必要なスペックを調べずに、ただただ「良いやつ」が欲しいといわれる方。
メモリ、HDD、GPU、ソフトを理解せずに無駄遣いをしているにもかかわらず、他の数万円の消費には厳しい方。傍から見るとバランス悪い人。

無知は罪です。

===========================
20191029追記
数GBのデータ保存用に、TBのHDDを買いたいとのこと。
もう少し調べてほしいです。

建築屋さんのこだわり

面白い階段を見かけました。

建築屋さんのこだわりを感じる階段です。
花崗岩中の暗色包有物を表現されたのでしょうか（一部の石は違うように見えましたので色だけ？）。

いずれにしても、どこまでの方が気付いていらっしゃるでしょうね。

絵でみる水文観測

「絵でみる水文観測」監修：国土交通省　中部地方整備局

基本的内容が絵で描かれており超初心者には良いと思います。古い内容ですが、基本は変わっていないため問題ないでしょう。

雨量計の設置に関して備忘録です。
===============================
設置位置
・障害物の高さの4倍以上離すことが望ましい。
・避ける箇所：尾根、谷、傾斜地、ビルの屋上、ビルの下、屋根の上
電線
・保護パイプに入った対候性に優れた2芯の電線
・途中でつながない

地震波と技術者

地震波に関して技術者と話をする機会があります。

微動や強震動を扱う物理屋さん。この方たちは周波数領域の中で生きていらっしゃるようで（本人談）、スペクトル、固有周期、卓越周期に関する話題が多くなります。
斜面防災に関しては専門外であるものの、周波数領域での議論がなされていないと「聞くに値しない」と思われるようです。

建築屋さんや橋梁屋さんは、地盤の固有周期や構造物の減衰を当たり前のように扱われています。こちらは基準に含まれていますので、思い入れというよりは設計作業の一つと捉えられているようです。
斜面防災の話に興味を持たれますが「専門でないためわからない」という一歩引いたスタンスで聞かれることが多いように思われます。

地質屋さん。断層や地質構造に関する話題が主体ですね。弾性波や表面波などを扱われていますが、波の議論はほぼされません。
斜面防災では被災要因（素因）として地質に着目されるようです。この話がないと「片手落ち」と感じるようです。

いろいろな立場の方がいらっしゃいますし、新たな視点を得られて面白い。
恵まれた環境です。

履歴減衰と再載荷

履歴減衰の話を聞いていて、一つ繋がりました。

片側に応力を与え続け塑性化を進行させた後、除荷・載荷を行うと、履歴ループが原点からズレます。逆方向に応力を加えず、もう一度片側に応力をかけ続けた場合、見たことのある形になります。これ、プレッシャーメーター試験や圧密試験での再載荷のループと同形状です。

地震の履歴減衰は Voigt モデル・sin波等から入って楕円の方程式に導かれるため、疑問を挟む余地がありませんでした。コンクリートの履歴減衰も、この説明をふまえるとイメージしやすいと思います。
が、物理的過程は上記とほぼ同じ。これがつながらなかったのは、やはり視点が固まっていたのでしょう。

かといって、プレッシャーメーター試験や圧密試験での除荷・載荷ループがダッシュポットでモデル化できるのか？と言われると、違う気がします。確かに、透水係数と時間がダッシュポットの役割を果たしているのかもしれません。が、弾塑性の全応力解析でもこのループは再現できます。土の場合、何がこの再載荷ループを作るのでしょうか？
https://phreeqc.blogspot.com/2011/12/4.html

ひとまず、イメージは繋がりました。
わからない範囲も絞ることができました。
もう少しのような気がします。

台風19号

台風19号でお亡くなりになられた方が日々増えています。ご冥福をお祈りします。

通過後の週明けには斜面防災の後輩君から関係者に速報が入っていました。メディアでは洪水被害が主に取り上げられ斜面崩壊は目立ちませんし、だいち防災WEBポータルの SAR の解析結果でもイマイチ見つけにくい。速報でも現状を把握するにはありがたい状況でした。

その後、慰問を兼ねて営業さんが情報収集に走りましたが、洪水跡に阻まれ立ち入ることもできない箇所が多かったようです。自治体の担当者はそれ以上に走り回っておられ、電話すら繋がらないところもあったとのこと。
1週間経ちましたが、まだ全容は把握できていないようです。

毎年起こる災害。
技術者として、何とかお役に立ちたいものです。

Amazon の危機管理？

9月末に発生した Amazon の個人情報漏洩。

「申し訳ない。もれちゃった。」
という内容の連絡が来て、被害を被ったことに気づきました。メール内容があっさりしていましたので何度か問い合わせをし、ようやく実態を掴めました。
今回、公的なアナウンスはないようです。

登録サイトを広げない、登録先毎にアドレスを変える、カードを登録しないなど、各企業の情報漏洩に対し被害を抑える対策をとっていました。が、いざ当たるとやはり残念。

仕事でも、情報漏洩対策はとっています。が、漏洩を前提とした危機管理手順は講じておかなくてはなりません。
少なくとも、「申し訳ない。もれちゃった。」「修正したからもう大丈夫だよ」だけでは危機管理とは言えません。
力関係もありますが、Amazon の対応にも学ぶべきことがあると気づかされました。

orange

orange
https://orange.biolab.si/download/

Anaconda に含まれるまれるツールです。前処理だけかと思いましたが、いくつかの機械学習を実装していました。
一通り試してみましたが、

便利です。

手元に大きなデータセットがなかったので、実務に耐え得るかどうかはわかりません。が、簡単な特徴量分析には使えます。

近年、AI といったバズワードのもと、各種ツールがそろってきました。いずれ programing なしでも結果を得られる便利な汎用ツールが揃うでしょう。

今は期待が過剰に進んでいる時期です。が、できること、できないことが分かった後、一つの学ぶべきツールとなるのは必然と思われます。

過剰な期待は捨てつつ、変化には追随しましょう。

学ぶべきこと

運動後に若い方たちと夕食。

子供のような世代の方たちでしたが、思った以上に多様でした。

凱旋門賞にネットで賭けている子
本省に転職する子
FXで溶かしている子
副業で儲けている子
全員キャッシュレス。

でも、悩んでいることは私の若いころと同じ。
環境は変わりましたが、変わらぬところもあるようです。
若い方たちに学び、持っているものは伝える。

そのような機会が多くなれば良いと思います。

STIV

STIV の文献はコチラ。

藤田2003「時空間画像を利用した河川表面波紋の移流速度計測」
https://ci.nii.ac.jp/naid/10018851453
藤田2009「STIVによる劣悪な撮影条件での河川洪水流計測」
https://ci.nii.ac.jp/naid/10025309107/

違いはコヒーレンスを閾値として使うか重みとして使うかぐらいですね。
読んでいくうちに実装が簡単そうに思えてきました。動画読み込み部分は PIV で作っていましたので、その続きで組んでいくことに。

実装は思った通り簡単で、3～4時間ほどで終わりました。が、Verification ができません。比較するソフトがないので、微妙な答えだとあっているのか確かめようがないのです。これ、重力探査の時もそうでした。
https://phreeqc.blogspot.com/2019/07/2.html

あと、時空間画像ができた時点で、見た目で勾配を決めてしまった方が良いかも。人の感覚は鋭いものです。綺麗に勾配を出せますし、処理も速い。
そうなるとPython のように CUI でなく、GUI で処理したいですね。勾配を決める以外にも、測線を決めたり、既知の測量点を画像から指定する際には GUI が効率的でしょう。実務で使うなら、市販ソフトを購入した方が良いでしょうね。

また微妙な結果になってしまいましたが、ひとまず理論は理解できました。

PIV と STIV

PIV を触っていたのが7年前。
https://phreeqc.blogspot.com/2012/04/piv.html

このPIV、今は OpenPIV の example だけで十分実装できるレベル。画像を2枚指定して、比較する大きさを決め、相互相関を取り、閾値以下のベクトルを表示するだけ。1時間ほどで実装できました。
ソースとしての動画の扱いや、その幾何補正は OpenCV。こちらも1時間ほどで組めました。簡単でした。

平成26年版の河川砂防技術基準では、主要な流量観測手法に PIV が含まれています。が、雨滴等に弱いためか、土研さんは STIV 等を勧めています。この STIV、論文発表が2003年なので古い技術です。
https://www.pwri.go.jp/team/hydro_eng/manual.htm
どちらも速度を出すために、時空間を扱います。
PIV は空間寄りの計算法ですが、STIV は時間寄りです。1次元のライン上で、輝度の時間変化を取り出し連結します。縦軸に空間（距離）、横軸に時間を取るように連結させると、その画像の輝度の傾きで速度を出せます。理解し易い発想です。

その後、平成29年には国土交通省の「水文観測業務規程」が改正され、「その他」として画像による河川流量の観測方法が加わりました。
http://www1.river.go.jp/
https://www.pwri.go.jp/jpn/about/pr/event/2018/1011/pdf/kouen4.pdf

手法は何であれ「できない」は避けたいところ。いつか STIV も実装しようかな。

*************************************
20191222追記

土研さんの「流量観測の高度化マニュアル（高水流量観測編）Ver1.2」D1より。

『河川砂防技術基準 調査編、第 2 章 第 4 節－3、平成 26 年 4 月』には、画像処理型流速計測法として LSPIV が示されている。しかし、LSPIV はパラメータ設定や撮影動画像の分解能に計測精度が依存するため、近年ではほとんど用いられていない。
一方、LSPIV に替わる画像処理型流速計測法として STIV、PTV として浮子を活用するFloat-PTV が実用化されつつあるため、本ガイドラインでは LSPIV は除外している。

https://www.pwri.go.jp/team/hydro_eng/manual/manual_D1.pdf

Neural Network Console

今回も SONY さん。

Neural Network Console
https://dl.sony.com/

Prediction One が全自動だとすると、コチラは半自動。プログラミング不要ですが、層構成を D&D で作る必要があります。逆に言うと、GUI でネットワークを構成できるので、各層の意味を理解していたら非常に簡単。LSTM、CNN も、GUIで組めます。

ハイパーパラメータ探索のみならず、層構成もベイズ最適化等を用いて自動最適化し、パレート最適解を示してくれます。
残念ながら、手持ちのデータで3日以上回しましたが、期待したほど広範囲を探索してくれませんでした。が、有用な機能なので、もう少し触る必要ありでしょう。

GPU も使えます。複雑な設定は必要ありません。
結果の PowerPoint への書き出しも可能。かゆいところに手が届いています。

このように設計は素晴らしいのですが、まだまだバグを潰せていないようです。1週間ほど使用しただけなのですが、training 段階で走ったデータが evaluation で走らない、PowerPoint 書き出しで学習曲線の一部しか載らない、混同行列が表示されない、等。エラーには再現性あり、なしがあり、クセを掴みにくいところが悩ましいですね。ま、バグ取りがてらの無償公開なのでしょう。有償になるかもしれませんが、次のVer.に期待しましょう。

今のところ、手持ちのデータでは最後まで走らないので、自動最適化で把握した構成を Keras 等で再現・調整するといった使い方になるでしょうか。ま、使い勝手の良い探索ツールが一つ増えたという認識で良いでしょう。

現場作業の基本

週末に観測機器を設置。

既存のネットワークに追加するだけなのですが、単管や電気配線を新たに配置するため、後輩君の練習がてらお手伝いとして参加。

設計主体の後輩君だったため、現場作業は不慣れ。
インパクト、電工ペンチの使い方や番線の締め方まで、ほぼ最初からでした。

番線の締め方やモノレールの使い方などは現場作業の基本です。が、現場主体の地質屋さんでも、御存知ない場合があります。おそらく、先輩から受け継がなかったのでしょう（動画でも確認できる時代ですが）。

ま、今回は「次はできる」と言わしめたので、目的は達成。
若い方も活躍できるように、基本的なことは学べる環境ができると良いと思います。

Prediction One

SONY さんの Automatic ML「Prediction One」
https://predictionone.sony.biz/

H2O の Driverless AI （最近、国内販売の代理店ができたようです）とコンセプトは同じです。

まずはKaggle の Credit Card Fraud Detection でテスト。
このデータは不均衡で、陽性0.2％です。ダメだろうなあと思いつつ前処理なしで喰わせてみました。
結果は上々。AUCが0.981。驚き。
残念なのは速度。GPUを使いませんし、CPUのコアも1つしか使用していません。メモリもかなり余っていました。

次に、手持ちのデータでテスト。同様に0.1%の Imbalanced data です。
結果はエラー。「予測ターゲットのユニーク数が少なすぎます。」と出ます。交差検証必須を外しても、学習・予測データを最初から分けてもダメ。わずか0.1％の差か？データの並びか？

前処理済みのデータではどうか？と思い over sampling 済みのデータでテスト。
が、これもダメ。60万行×40列なので仕様上は問題ないはずですが、学習を始めると1時間程度で停止します。

残念ながら前処理で python を使いますし、大きなデータを扱わないわけにはいかないので、実務には向かないようです。が、そこそこ精度は出そうな点、評価に資するための十分なレポートを作成してくれる点、ドキュメントに基礎知識等の解説が丁寧に解説されている点、日本語を扱える点から、入門者には良いツールでしょう。

国産 Automatic ML。Driverless AI レベルまで進化すれば  great machine learning expert の不在は重要でなくなるかも（前処理に強い方は必須です）。
https://phreeqc.blogspot.com/2019/09/rules-of-machine-learning.html

SONY さん、応援します。

全国の地下水流動経路

近年、GETFLOWS の利用例をよく見かけます。

昨年、地下水学会誌でも「水循環基本計画の下での地下水に関する取り組み」が特集として組まれていました。こうなると地表水、地下水の連成が必須になります。できないだけは避けたいのですが、まだまだ実務で自由に使えるほど実装できていません。

と、見ないふりをしている間に、このようなサイトが公開されていました。
https://www.getc.co.jp/webmap/

trial ということは、有償版があるのでしょう。
うーん。素晴らしいと思う反面、コード非公開のため盲信してよいのかわかりません。相変わらず微妙なのですが、一部ではネームバリュー含め利用価値はあるでしょう。
うーん。

Rules of Machine Learning

Rules of Machine Learning
http://martin.zinkevich.org/rules_of_ml/rules_of_ml.pdf

To make great products:
 do machine learning like the great engineer you are, not like the great machine learning expert you aren’t.Most of the problems you will face are, in fact, engineering problems. Even with all the resources of a great machine learning expert, most of the gains come from great features, not great machine learning algorithms.

ないものねだりでしょうか。
機械学習のプロも傍に欲しいと常々思っています。より良いモデルができるのでは？時間を短縮できるのでは？云々。
そのようなときに目にした文章です。有名なのか、全訳されている方もいらっしゃいました。ありがたい。

不飽和地盤を対象とした現場透水試験

2年前に不飽和地盤を対象とした現場透水試験法が、地盤工学会で基準化されました。

それまで、現場透水試験といえば、飽和帯を対象としたものが主体でした。
過年度の報告書を借りてくると、不飽和帯に対し飽和帯の試験を適用して、透水係数を求めている例を見かけることも屡々。「長時間注入しているから飽和したと見なす」、「流量が一定になったから飽和帯とみなす」、などという期待のもとに実施されていたり、何も記載せず淡々と飽和帯の式を流用している例も多く見ました。背景には、不飽和の式が基準にないから飽和の式を使うといったマニュアル依存が強くあったようです（基準外の論文ベースの不飽和試験を実施している例は見ませんでした）。

「長時間注入しているから飽和したと見なす」等の期待について、論ずるに値しないことを示すデータは、以下の報告にあります。

西垣ほか「不飽和地盤を対象とした現場透水試験法に関する課題の抽出と改良に関する考察」

https://www.jstage.jst.go.jp/article/jagh/60/3/60_289/_article/-char/ja

• 初期飽和度により求められる透水係数が異なる。これは、不飽和帯に注入することで間隙空気が封入されるため。
• CO2ガスの事前注入が有効。

大きくは上記2点の成果が中心の報告です。他にも、示唆に富む報告が数多く含まれており、個人的には「当たり」の資料でした。

過去の研究をもとに、普遍化・共有化できる部分のみがマニュアル化・基準化されているとも言えます。その観点ではマニュアルや基準の適用範囲に多様な自然をすべて合わせることはできません。適用外の自然をどのように調べるか？が地質屋の力量と言われたこともあります（素晴らしい施主でした！）。

今まで力量を問われ対応できなかった方々も、せっかく基準化されたのですから、使えばよいと思います。不飽和地盤を対象とした試験、今後は広く使われるようになることに期待しましょう。

土石流と水温

土石流発生のタイミングを水温で予測しようとする研究です。

【研究成果】土石流発生リスクを地下水の温度で予測

https://www.hiroshima-u.ac.jp/news/53312

降雨がある程度強くなると、地下水位が上昇するために、地盤内でも下から上に水分が広がっていきますが、通常は地表から地下に降りていく雨水によって地表面の熱が伝わっているために、地下水の温度は上昇します。

ところが、平成30年7月豪雨では地下水位上昇時に水温が急激に低下することが観測されました(図)。

これは、地下深いところにある岩盤内の水圧が高まって、岩盤から表層の風化層に冷たい水が供給されて地下水位が上昇したことを示しており、豪雨時に基盤から表層斜面に水が供給されることを裏付ける貴重なデータとなりました。

この研究、ビジネスに展開するには次の3点をクリアしなければならないでしょう。
・非発生斜面との比較（非発生斜面では温度低下が認められない？鈍い？）。
・流域の全斜面はもちろん、全国に温度計を設置するのは不可能。代替案が必要。
・短期雨量に起因する土石流発生。

1点目。研究の根源でありビジネス以前の問題です。今後、他の谷での観測結果を待ちたいところです。
2点目。文献では1つの谷に対し5箇所で観測を行っていました。5箇所でも少ないのに、他の斜面、全国展開となるとお手上げでしょう。が、契約民家裏の谷のみをビジネスの対象とする、土砂災害警戒情報にあわせてメッシュ対応にする、数値実験や SWI を利用するなど、代替え案は考えられます。
3点目への対応は厳しいでしょう。先日の岡山県北部で発生した土石流は拾えないでしょう。

いずれにしても、まだこれからの研究。
今後のデータ蓄積・分析に期待しましょう。

土石流と降雨パターン

学主導の岡山県での土石流調査に同行しました。

継続3時間程度の短期雨量が効いた災害です。土石流は長期雨量で発生しやすい≠短期雨量は効かない、です。

現場を見ると、今まで見たことのある土石流現場とは感覚的に異なっていました。発生土砂量や分布範囲はやや少なく思える程度なのですが、粒径が小さい。渓流の洗い流しも少なかったと思われます。地質や傾斜の違いがあるのかもしれませんが、降雨パターンの違いの方が大きいように感じました。

帰路の途中、先生が「GRIBデータの入手が遅れる」と話されていました。購入されてはいないようです。学なら京大の整理したデータが使えると思いますが、欲しいデータは含まれていないのかもしれません。
http://database.rish.kyoto-u.ac.jp/arch/jmadata/
ArcGIS ならそれをアニメーション表示できます。
https://blog.esrij.com/arcgisblog/wp-content/uploads/2016/08/ArcGISforMeteorology.pdf
Arc を購入されていないようでしたら、Git で公開されているコードを利用する方法もあります。便利な時代になりました。
https://github.com/catdance124/synthetic_radar_GPV


今回、発生土砂量、流動範囲、粒径についても予測モデルに入れるべき特徴量だと実感しました。今後、検討する機会があれば意識して入れてみましょう。

大規模・長期停電

千葉県の停電が4日目に。

近年、災害発生時の情報伝達にスマホは必須となっています。しかし、今回のような大規模かつ長期の停電、通信障害では全く役に立ちません。スマホへの過度の期待が浮き彫りとなりました。

キャッシュレスは以前より災害時の虚弱性を指摘されていました。停電が1週間続けば、現金でも危ういかもしれません。

電線を地中に埋設すればここまで大規模になっていなかったと思われます。が、復旧は長引くでしょう。

今後も多くの課題が浮かび上がりそうです。

しかし、今はそれどころではありません。
現場の方々のプレッシャーや負担は計り知れません。今は作業員を含めた現地の方々の安全衛生環境の維持と早期復旧を願わざるを得ません。

建設産業と消費税増税

9月に変更契約になる仕事。

あれ？消費税はどうなるの？
ということで、営業部門に確認。以下を紹介してくれました。

建設産業における消費税の転嫁対策について
http://www.mlit.go.jp/totikensangyo/const/totikensangyo_const_tk1_000063.html

ポイント①どの時点で課税されるのか？契約日ではなく、「引渡し日」時点の税率が適用されます
ポイント②経過措置とは？ 消費税率引上げの半年より前に締結した契約は、旧税率が適用されます
○消費税率１０％適用に係る指定日・・・平成３１年４月１日

なるほど。引渡しが10月を超える場合、
・4月以降の契約は10%
・3月以前契約の繰り越しは、変更分のみ10%
ですね。
各自治体もこれに倣い、工事・コンサルタント業務の両方に適用されているようです。

これ、下請契約にも触れられています。
元請が4月契約10月納期であったとしても、下請契約の工期が9月末であった場合、8%対象になるのでしょうか？その理屈だと、9月末で一旦契約を完了する元請が多く出てきそうですが。また、減額変更の場合はどうなるの？
ま、そのあたりはきちんと整理されているのでしょう。

物品購入で10月納期でも、8%の会社があったり10%の会社があったり。
適正な処理ができているのか不安になりますね。営業・経理さんに頼りっきりになりそうです。

地下水の年代測定を省力化する採水法

先月末、農研機構さんが以下の手法を公開されました。

「省力的な採水法による六フッ化硫黄を指標とした地下水の年代測定」
http://www.naro.affrc.go.jp/publicity_report/press/laboratory/nire/131962.html
http://www.naro.affrc.go.jp/publicity_report/publication/pamphlet/tech-pamph/131331.html

以前、プロに教えて頂いた手法とほぼ同じですので、農研機構さん独自の手法というわけではありません。空気に触れさせない簡単な方法を考えると、同じところに落ち着くのでしょう。

資料の中で、定量的に評価されている点がポイントでしょう。2地点のみですが、SF6を使うにしてはやや古い時代ですので、一般化できているものと期待したいところです。

即時残余耐震性能判定システム

ニュースを見ていますと、住宅に安価な地震計をつけて損傷を評価するシステムが紹介されました。

住宅の１階と２階に各1台（計2台）の地震計を付けています。振動が閾値を超えると地震と判定し、同時にTVに損傷度が出るようです。通常の波形も出るので、欲しくなりました。
ただ、どのように評価するのかまでは紹介がありませんでした。損傷があると固有振動数に変化が出るのでしょうが、どこまで変化すれば「危険」と判断するのでしょう？

気になって調べたところ、引っ掛かりました。
東京大学地震研究所  楠 浩一「即時残余耐震性能判定システム」
http://www.eri.u-tokyo.ac.jp/wp-content/uploads/2018/10/ERI-nl-plus_No29-web-A4.pdf

「ウェーブレット変換を用いることで、安価な加速度計の測定値でも2階積分してノイズを除去し、安定して変形量を算出できるようになり、判定システムが現実味を帯びてきました」

なるほど。すばらしい。
加速度計なんですね。2台の差分を積分して2階の相対変位にまでもっていく。F=maで力を出して、その勾配を変形係数とみなし、降伏・破壊を判定する。というところでしょうか。いえ、単純な力の出し方だと勾配は折れないですね。力は１階の加速度のみで出せば良いか。ずれるかな？
地震検知や2階積分のアルゴリズムは既にあるので、この程度なら電子工作でもつくれそうです。

理屈だと固有振動数の変化やスマホアプリでも判定可能です。が、住宅メーカーの協力がないと難しいでしょう。というか、メーカーさんが付加価値として作られるかもしれませんね。

高いハードがなくても、アイデア次第で、既存の技術で新しいビジネスに挑戦できる。すばらしい例だと思います。
私に足りない想像力です。脱帽。

OpenACC 成功？

寝かせていたコードを、起こしてみました。これの続きです。
https://phreeqc.blogspot.com/2018/04/openacc.html

・PGI Community Edition 19.4

あっさり、計算が進みました（もう慣れました）。

結果を見ると、精度は許容範囲。
ただし、計算時間は全く短縮していません。GPUは使っていますが、計算に使ってくれているのでしょうか？

Intel のコンパイラーはよくできているとあらためて関心。
また寝かせましょう。

******************************************
20190925追記
プロと話をしていると、NVIDIA さんは数値解析分野で思ったほど力を発揮できなかったため、機械学習分野にシフトしていった、とのこと。そういえば、以前、GPGPU用にチューニングされたライブラリを使わないと速くならないとアドバイスを受けていました。
そうなのかなあ？

現場管理

昔、他社さんと一緒に働いていた頃の話。

「いかに現場に行かないで仕事を済ませるかが、利益向上につながってしまいます」

「どこの会社も同じですよ」

私が勤めている会社では、現場に行く頻度・時間は個人の裁量です。現場稼働中、私は基本、毎日現場にいますが、他の方はそうでもないようです。

昨日、下請けさんからの話

「○○会社さんが言ってましたよ。○○さんが現場にこなさすぎると。」

「・・・そうだと思います。」

「○○会社さんは初めてだし若い方が多いので、もっと大事にしないと。」

「・・・ありがとうございます。」

今年度から私の所属が変更され、今後、現場に出ることはほぼなくなりました。が、最後となろう現場にて、ありがたい助言をいただきました。

ただ、それを社内で伝えたとしても彼らなりに言い分があり、改善はなされないでしょう。経験上、これまで問題が発生していないことも大きな判断材料になっているのだと思われます。将来は別として。

コンプライアンス、バイアス、リスク管理。最低限の知識があれば、すぐに気づく内容です。該当部署のリーダーは総監取得にも苦労されているようなので、まだ先でしょうか？

危ういですねえ。

残念なオジサマ

講習会に来られていた少し残念なオジサマ。

小人数であったためか、講師が話されている途中で質問を投げかけます。熱心に質問・議論をされる姿勢は良いことなのですが、そのレベルの低さ、頻度、口調には閉口。もっと勉強してから来られるべきでしたね。
度々、進行が止まるため周りの方も迷惑そうでしたが（受講者がしびれを切らして回答される場面もありましたが）、一番前の席だったためか、元からそのような性格なのか、最後までお気付きになりませんでした。

残念ながら、私にも少し身に覚えがあります。
社内会議でレベルの低い技術報告が出てきたときには指摘すべきかどうか迷います。executives が（知ってか知らずか）指摘しないのだからこのままで良いのでは？などとも思うのですが、意見することはあります。が、たいていの場合は話になりません（基礎を疎かにしていなければ、最初から誤った報告は出されないでしょう）。で、結果的には多くの人が共有すべきでない時間を作ってしまいます。さらに空気を読む必要があるのでしょう。

若手なら育てる義務がありますが、オジサマを育てる義務はありません。
無知は罪。自分はそうならないよう、毎日コツコツ積み上げましょう。

SAR 緊急観測オーダーと課題

だいち防災 web。
久しぶりに訪れると、公開内容が増えたように感じました。気のせいでしょうか？
https://jaxa-dis.maps.arcgis.com/home/index.html

危機管理のための SAR 緊急観測オーダーが普及しつつあるようです。このサイトを見る限り、お盆前後の台風10号では、紀伊半島や中四国等にオーダーが出されていたようです。
今回の佐賀県豪雨では、28日昼に確認した時点で強度画像の合成結果と手動判定結果が公開されていました。が、使っているのは28日０時撮影のデータ。判定された浸水域は既に過去のものとなっており、危機管理の視点では情報が古く使えないでしょう。図らずも、次の課題（リアルタイム把握）が見えてきました。

ま、せっかくの衛星ですので、このような目的でも利用頻度は上げて頂きたいですね。

pooling

配列を pooling する方法です。

CNN ではなく、ただ単に配列を最大値で pooling したいと考え、その方法を探しました。
最初は numpy や scipy にあるだろうと簡単に考えていたのですが、見つけることができず。

最終的には、tensorflow の max_pooling を利用しました。このような感じ↓
******************************************

In [1]:

import numpy as np
from scipy.ndimage.filters import maximum_filter
import tensorflow as tf

In [2]:

x1=np.array([[0,1,2,3],[8,9,10,11],[4,5,6,7],[15,14,13,12]])
x1

Out[2]:

array([[ 0,  1,  2,  3],
       [ 8,  9, 10, 11],
       [ 4,  5,  6,  7],
       [15, 14, 13, 12]])

In [3]:

x = x1.reshape([1, 4, 4, 1])
x

Out[3]:

array([[[[ 0],
         [ 1],
         [ 2],
         [ 3]],

        [[ 8],
         [ 9],
         [10],
         [11]],

        [[ 4],
         [ 5],
         [ 6],
         [ 7]],

        [[15],
         [14],
         [13],
         [12]]]])

In [4]:

x.shape

Out[4]:

(1, 4, 4, 1)

In [5]:

maximum_filter(x1,size=[2,2])

Out[5]:

array([[ 0,  1,  2,  3],
       [ 8,  9, 10, 11],
       [ 8,  9, 10, 11],
       [15, 15, 14, 13]])

In [7]:

pool=tf.nn.max_pool(x, ksize=[1, 2, 2,1], strides=[1, 2, 2, 1], padding='SAME')
se = tf.Session()
pool = se.run([pool])
np.array(pool).reshape([1,2,2])

Out[7]:

array([[[ 9, 11],
        [15, 13]]])

メモリオーバー

興味本位でいくつかを対象に FFT や PCA をかけています。

これが興味深い。想像してなかったことや、「考えると当前」といったような結果が浮かび上がります。先の 1/f もそうですね。

この過程で、頻繁に遭遇するのがメモリーオーバー。
48GB 積んでいますが、簡単にオーバーします。GPUならなおさら。前処理としてスライスしたり、pooling したりするのですが、ダメですね。ビッグデータどころか1時間のデータも扱えません。

簡単なのは物量作戦。メモリが足りなければ積めば良い、です。が、ある程度で頭打ちになるでしょう。何とか、効率よい処理方法を考えないと。FFT はできそうですが、PCAは思いつきません。
うーん、どうしましょう。

1/f ゆらぎ

動画を時間軸で FFT。

予想通り、スペクトルにすると低周波が高く、高周波が低くなる形状です。微動では見慣れない形状ですので面白いなあと思いつつ、当初の着目点に戻って整理を進めました。

昨夜、文献を眺めていますと、この形状についての説明がありました。

寺西ほか「動画像の時系列周波数特性によるゆらぎ解析」法政大学情報メディア教育研究センター研究報告 19, 149-153, 2006

周波数の低下とともにパワースペクトラムが増加 するような信号の中で、パワースペクトラムの振幅が周波数に対して反比例する信号を「１/ｆゆらぎ」と呼ぶ。

これでしたか！と目が覚めた気分。
昔、「１/ｆゆらぎ」が商品名についていたのを覚えています。が、意味も違いも分かっていませんでした（現在、そのような商品はありませんので、多くの方が違いを感じなかったのでしょう）。ここにきて、周波数との逆相関=「1/f」だと、ようやく理解しました。

これ、自然に多くあるとのことですが、本当でしょうか？
少なくとも、微動にはありません。
雨音は単位時間であれば高周波が卓越しているように思われます。が、音の場合、ランニングスペクトルで評価するのではなく、単位時間を代表する1つの数値を時間軸方向に並べ、そのスペクトルを取る必要があるのでしょう（時間軸での「ゆらぎ」ですから）。そうすると「1/f」に近づくのかもしれません。
画像だと、グレースケール変換、畳み込み、固有ベクトルなどが代表値として考えられます。風は風圧で良いでしょう。が、音は音圧ではダメでしょう（音色の方が重要な気がします）。知りたいですね。そのあたりは昔研究されているでしょう。探せばいくらでも出てきそうです。

まだまだ面白い会遇があるものです。
色々手を出してみるものですね。

numpy での処理

python で3次元配列を扱うことが立て続けにありました。

これらを空間方向、時間軸方向に扱う処理方法が求められます。が、最初は効率的な方法がわかりませんでした（そもそも、3次元配列に整形する段階で躓きました）。

意外と web 上に情報がなかったので、最終的に採用した numpy での処理方法を書き残しておきます。
******************************************

In [1]:

import numpy as np

In [2]:

x1=np.array([[1,1,1],[2,2,2],[3,3,3]])
x2=np.array([[11,11,11],[22,22,22],[33,33,33]])
x3=np.array([[111,111,111],[222,222,222],[333,333,333]])

In [3]:

x=np.array([x1,x2,x3])
np.shape(x)

Out[3]:

(3, 3, 3)

In [4]:

x

Out[4]:

array([[[  1,   1,   1],
        [  2,   2,   2],
        [  3,   3,   3]],

       [[ 11,  11,  11],
        [ 22,  22,  22],
        [ 33,  33,  33]],

       [[111, 111, 111],
        [222, 222, 222],
        [333, 333, 333]]])

In [5]:

x[:,0,0]

Out[5]:

array([  1,  11, 111])

In [6]:

xx=[]
xx.append(x1.copy())
xx.append(x2.copy())
xx.append(x3.copy())
xx=np.array(xx)
np.shape(xx)

Out[6]:

(3, 3, 3)

In [7]:

xx

Out[7]:

array([[[  1,   1,   1],
        [  2,   2,   2],
        [  3,   3,   3]],

       [[ 11,  11,  11],
        [ 22,  22,  22],
        [ 33,  33,  33]],

       [[111, 111, 111],
        [222, 222, 222],
        [333, 333, 333]]])

In [8]:

xx[:,0,0]

Out[8]:

array([  1,  11, 111])

In [9]:

np.mean(xx,axis=0)

Out[9]:

array([[ 41.,  41.,  41.],
       [ 82.,  82.,  82.],
       [123., 123., 123.]])

In [10]:

np.mean(xx,axis=1)

Out[10]:

array([[  2.,   2.,   2.],
       [ 22.,  22.,  22.],
       [222., 222., 222.]])

In [11]:

np.mean(xx,axis=2)

Out[11]:

array([[  1.,   2.,   3.],
       [ 11.,  22.,  33.],
       [111., 222., 333.]])

In [12]:

np.mean(xx,axis=(1,2))

Out[12]:

array([  2.,  22., 222.])

In [13]:

np.mean(xx,axis=(0,1))

Out[13]:

array([82., 82., 82.])

衛星データ

オープンな衛星データを扱えるサイト。AWS でも扱っているのは知りませんでした。

https://www.sentinel-hub.com/
https://earthengine.google.com/
https://registry.opendata.aws/

Tellus でも扱っていますが、日本の首都圏主体です。せめてALOS 程度は全国整備してほしいところです。

衛星データと防災

SAR分析チャレンジ
https://sardac.jp/

4月から5月にかけて開催されていた、情報通信研究機構 機構の啓蒙活動。基本的なコードを GitHub 上で配布し、それを組み合わせる＋新たに組むことで、SAR データと防災に関するアイデアをリンクさせるハンズオンだったようです。
公開された成果を見ると、皆さん面白いアイデアを出されていました。小さな自治体を対象としたアイデア主体で、これまで緊急観測の視点から広域利用を考えていた（そして宇宙村から外れている限り現実的でないと距離を取っていた）私にとって、新たな視点からの、驚かされる成果が多くありました。これからの利活用が期待されるところです。

SARの活用は古くから検討されています。が、SIP 等で提案されているモニタリング関連の技術レベルは数年前から変化なし。ソフトを買えば容易に達成できるレベルですが、高額なためか流行っていません。
コンペ以外にも、上記のような啓蒙活動を通しユーザーを増やすことは、技術レベルを進める一歩かもしれません。
https://www.jst.go.jp/sip/k07_kadai_01.html
http://pixel.eri.u-tokyo.ac.jp/sarws2013/sarws2013.html

一方、ArcGIS のブログでは、可視画像を使った土砂災害発生個所の抽出例が載っています。これは航空写真だけでなく、衛星の光学画像でも使えるでしょう。LANDSAT での抽出は昔から知られていますので、データ配布だけでなく結果の見やすいサイトもできています。
https://blog.esrij.com/2018/09/20/post-31278/
https://www.jstage.jst.go.jp/article/journalofmmij/132/6/132_96/_article/-char/ja/
https://apps.sentinel-hub.com/sentinel-playground/　もしくは　https://www.sentinel-hub.com/explore/eobrowser

衛星データやDEMなど、身近に使えるデータはたくさんあります。'おいしい'データは制限がかかっていますが、それでもアイデア、プログラミングの腕、やる気があれば、まだまだビジネスチャンスの埋もれている可能性があります。「SAR分析チャレンジ」によって、あらためて認識させられました。

k-空間

東京電機大学出版会「リモートセンシングのための合成開口レーダーの基礎」p325より

g(x) のフーリエ変換 G(f) の f
1/距離：空間周波数（spatial frequency）
1/時間：周波数（frequency）

空間周波数 fの代わりに k=2πf がしばしば用いられる。
k-空間（またはフーリエ空間（Fourier domain/space)）：k=2πf
kは電磁波等の波数(k=2π/λ)と混同しないように。

先日の重力探査の文献では「kは波数」等と明記されていました。一方、Fourier domain の wave number と表記された文献もあったように思います。もう一度読みましょうか。

波の知識

微動の整理を進める中で、必要な知識も一緒に整理しています。

有限フーリエ近似、有限フーリエ変換とその複素表現、フーリエ積分、周波数領域・時間領域でのフィルタや積分、窓関数、畳み込みなど。電気回路や画像処理にもつながっており、必要に迫られている技術です。これらを実装しながら、理解していない箇所をつぶしています（全く分からないところも残っていますが）。

お盆ですがプロが出社してましたので、いくつか質問もできました。おかげで解決できた問題もあります。ありがたい環境です。

一つ山を越えると、次の山が見えてきます。
このまま進みましょう。

ObsPy その2

ObsPyで計算した変位が、何かおかしいと悩んでいました。

他のソフトで計算した結果と大きく異なります。EXCELで計算した値ともあいません。速度も微妙に違っています。
加速度はあってたはず・・・と思い再確認したところ、ch２とch3 が異なっていました。15分のうち、ある1分間のデータだけがおかしい。
Mergeの仕方がまずかったか？と、その部分だけを読み込みましたが、おかしな値が再現されます。では係数の乗算がまずかったか？いえ、読み込み直後の生データの段階でおかしくなっていました。
処理ではなくて、ライブラリの想定するWin データではなかったということでしょう。プロ曰く、Winデータには方言があるそうですので。

地震の世界も案外狭いようで、その折々、必要なソフトを必要な方が作られているようです。いくつかのソフトや頂いたコードでは、今回欲しい結果に至るまでに2種以上の組み合わせが必要になりました。ObsPy＋Python なら、一気通貫できると思ったのですが、ダメでした。他のソフトで処理した波形やスペクトル比を TXT で読み込み、図化するだけなら対応できますが。うーん、これでは今までと一緒です。
結局、既存ソフトの機能追加をプロに依頼しました。

結果は出せませんでしたが、手を動かすことで理解は深まりました。それだけでも価値はあったと考えます。
わかっていたつもりでしたが、まだまだ初心者。あらためて分かった「わかっていないこと」について、今後身に着けていきましょう。

ObsPy

微動の結果を整理しています。

専用のソフトはあるのですが、自分でも中身を追って理解を深めようと思い、この週末に組んでいました。使ったのは最近のお気に入り、python。既にfortran のコードがあるので、それを片手に組んでみました。

まずは win データの読み込みから。
ライブラリがあるでしょう、と探してみたら、ありました。

ObsPy
https://docs.obspy.org/

微動ではなく、地震波を処理するためのフレームワークのようです。公開サーバーからデータを入手し、分析・図化することが容易にできます。進んでいますね。

obspy を使って、1分の win データを読み込むと、データが1個の stream として認識されました。時間やサンプリング周波数などのヘッダ付きデータです。各チャンネルはその下位の trace に入り、data で読みだせました。
フォルダ内の全 win データを merge で連続データに集約し、plot で図化。integral 2回で加速度を変位に。あとは matplot で図化。particle motion の出来上がりです（この機能はなぜか見当たりませんでした）。簡単。

integral の有効数字が気になるのと、窓関数の使い方がわからなかったのですが、それは今後。
いずれにしても、便利な世の中になりました。

******************************************
20190807追記
winデータの読み込みに、問題がありました。残念。

地すべりと微動

地すべりと微動に関する文献を読んでいました。

集めた中で最も古い文献が1973年です。意外と昔から研究されてきたようですね。港湾基準に取り込まれてから流行りだしたのかな？と思っていましたが、知りませんでした。

（気になった）主な内容は、以下の通りです。
・地すべりブロック内では、3~4Hz付近にて増幅率の増す報告が複数あり。
・地すべりブロック外でも同様の固有振動数を持つが、増幅はやや小さい。
・2~5Hzの particle motion では、表層のクラック直交方向の揺れが大きく見える。
・地質毎に卓越周波数が異なる。
　シラス:2Hz、花崗岩:6～10Hz、四万十層群:15～20Hz。
・表面波探査のS波速度と1/4波長則にて、地すべり深度を推定。誤差15%程度。

残念ながら報告例が少なく、微動が地すべりに十分活用されてきたとは言えないようです。盛土では動的解析を入れて微動を再現している例もありましたが、地すべりではまだまだ。平地で測る方が楽ですし、興味を持たれる方も少ないのでしょう。

以下、読んだ論文の概要です。
******************************************
島
常時微動と地すべり地への応用
地すべり　1973 年 9 巻 4 号 p. 1-8
・N-S、E-W方向とも大体同じような傾向。卓越周期ははっきりしない。→下層と上層の速度比はあまり大きくないか、あるいは表層の粘性係数が大きいか、また両者が一緒に原因しているかどちらかである。
・常時微動の観測をふつう平地でおこなうと、100mぐらいはなれたところで大きく変化しない。ところが、地すべり地をふくめて山地になると測点間の関係がかなりくずれる。このことはボーリング結果をみても同じであるから、やはり地下構造の変化を反映しているものとみられる

和田ほか
地すべり地の Crack 群の雑微動に対する影響
地震 1973 年 26 巻 4 号 p. 316-325
・亀ノ瀬地すべり
・2~5Hzのピーク：地すべりCrackに直交。crackが振動特性に強く影響。
・低周波領域（0.5~0.7Hz）：地すべり地全域で一様な振動特性。深部構造を表している。
・10Hzのピーク：断続的な振幅の強弱。人工ノイズ
・地すべり地のような地盤構造の複雑な地域では、単一の卓越周期のみに注目し、 垂直方向のみの地盤構造と対比するのは危険。
・地盤構造のディメンジョンに相当する個々の周波数領域についてそれぞれ考察することが必要
・地域下では潜在している crack 群 の発見と地すべり運動との関連性等を提供する可能性あり。

島
長野県鬼無里村地すべり地における常時微動と表層の振動特性
地すべり　1975 年 12 巻 1 号 p. 34-42
・粘性係数が大きいほど（Q値が小さいほど）減衰大
・地すべり崩落地に近いところでは粘性係数ξ=5×106～107CGS程度(Q値1～2程度）
・直接地すべりの影響をうけないところではξ=5×105CGS程度(Q値20程度）
・地すべりの発生からある程度時間が経過すると、地盤の弾性的性質は元の安定した状態に復元するものと思われる。表層の粘性係数は小さく（Q値は大きく）なり、スペクトルからみて明瞭卓越周期がみとめられる。

泉谷ほか
奈良尾地すべり地における常時微動特性
地すべり　1978 年 15 巻 3 号 p. 17-22
・仮定1：常時微動は基盤からのSH波垂直入射に対する地盤の応答である.
・仮定2:常時微動の源は時間的、 空間的にランダムに分布しており、 基盤より入射する波のy、 z成分のスペクトルは相等しい.
・仮定3:地すべり地のクラック群による土塊の異方性に着目して常時微動記録を解析することが、地すべりの状態を知る一つの手段となり得る。
・地表層上部を異方性体で近似した地盤モデル
・クラック群の到達している平均的な深さと、 クラックの混み具合とを推定できる解析手法を見出した.

泉谷ほか
常時微動測定による地すべり地盤調査の一手法
土木学会論文報告集 1981 年 309 号 p. 159-162
・「奈良尾地すべり地における常時微動特性」の改善
・常時微動によって推定された地下構造と、ボーリング調査等が対応するかを調べるにとどまっている.

島
常時微動特性から推定される味大豆地すべりの発生機構
地すべり　1981 年 18 巻 1 号 p. 15-25_1
・これまで、地すべり付近では振幅の大きくなる振動数はあらわれないで、振動数とともに振幅が小さくなるか、白色雑音のように振動数に関係せず振幅がほぼ一定になるようなスペクトルを示すことが多い。
・地すべり地下方では3～5Hz付近の振動数でスペクトル振幅が卓越する傾向を示す。
・味大豆地すべりも同傾向。
・地すべり地のように地形や地下構造が複雑なところでは、水平・鉛直成分の組み合わせによる楕円軌跡は水平軸に対して傾斜していることが多い。
・レイリー型の表面波であるなら、傾斜角は層構造の傾斜とみることができる。
・振動数ごとに傾斜角が異なっているときは、高振動数においてはごく地表の構造の傾斜に対応し、低振動数の場合は深い構造の傾斜に対応する。
・常時微動のスペクトルから地下構造を推定する場合、成層構造の中に低速度層を仮定してみると都合がよい。
・振動方向の分布にあずかる振動数と地下構造とを対比させると、低速度層より上部にある比較的硬い層が地すべり運動に寄与していることが分った。

川邉ほか
地すべりに及ぼす地震動の影響
地すべり　1983 年 36 巻 2 号 p. 5-16
・静岡県由比地すべり地において地震および常時微動の観測を行った。
・深度別二組の地震計の記録からスペクトル比を求め、 それをもとに重複反射理論を使って地表層の密度・剛性率・Q値の諸定数を推定した。
・固有周期で顕著なものはA点0.28秒(3.6Hz、増幅率6～10倍）、B点0.10秒 (10Hz、増幅率は5倍程度）
・0.1秒付近の振幅の増幅には深さ10m程度までの表層が関与。 0.3秒付近の振幅の増幅には、より深い幅が関与しているようである。
・推定された定数と重複反射理論を用いて、 基盤から地表層にある地震波が入射した場合に地表層内に発生する加速度 ・勢断応力を計算し、 土質試験より得 られたせん断強度と比較した。
・A点における最大加速度は、地表で最大で約3倍に増幅。最大勢断応力は境界面で最大約2.5kg/cm2となった。
・B点の各深さにおける最大せん断応力と、土質試験によるそれぞれの深さでの一面せん断強度を比較すると、 入力地震波の最大加速度が477ga1以上で表層内に破壊面
ができる。
・求められた加速度・勢断応力を由比の4個所の地すべり斜面の安全解析に導入した結果、 安全率が1になるのは傾斜約25度の斜面で50～110gal、約15度では100～220ga1のときであった。

秀島ほか
地すべり斜面における常時微動観測と一考察
開発土木研究所月報第457号　1991
・地すべり地域内の各成分の増幅は、地すべり域外のものと比べていくぶん大きい
・上下動成分は、地すべり域内および域外とも水乎2成分と比ぺて増幅は小さい
・地すぺり域内外でも、スペクトルの卓越振動数に相違はほとんど認められない。洪積層（第二種地盤）の基盤の振動特性そのものが強く反映された結果と考えられる。
・ハンドパスフィルクーは、既応の解析例をもとに2-5Hz。N-S方向(すべり頂部のクラックに直交する方向）が卓越した振動。

川邉
斜面表層の振動特性と不安定化
日本地すべり学会誌　2005年 42 巻 2 号 112-114
・地質毎に卓越した周波数
　シラス:2Hz
　花崗岩:6～10Hz
　四万十層群:15～20Hz
・表層部の固有周期、地質や土質、層厚などによって決定される地盤の振動特性を反映。
・下層から入射する地震動の卓越周期との関係から、表層部での安定性をある程度説明することが可能。

森ほか
微小地震観測による地すべり土塊の三次元形状と地震応答特性の評価
士木学会論文集Al(構造・地震工学）、 Vol.68、 No. 4(地震工学論文集第31-b巻）、 1_395-1_406、2012
・当地すべり地の微動は振幅が低く、電気ノイズに埋もれてしまうことが多い．微動探査は有効に用いることができない．
・微小地震観測は、地すべり地における卓越振動数の評価に有効。
・地すべりブロックの外側の安定した地山を基準にしたときのブロック内側の移動土塊部分の水平動スペクトル比（H/H比）は地震応答特性の評価には有効である．
・表面波探査等によるS波速度と合わせて地すべり深度を推定可能。誤差15%程度（H=Vs/4f、f=4H/Vs）

芝崎ほか
複数深度での地震動観測結果に基づく地すべり土塊の固有周期
日本地すべり学会誌　2016 年 53 巻 6 号 p. 227-234
・譲原地すべりおよび由比地すべりに設置された３深度（想定すべり面よりも下位の基盤岩層、想定すべり面の直上部および地表付近）の地震動観測を行い、高速フーリエ変換 によりスペクトルの増幅率を求め、地すべり土塊の固有周期について検討を行った。
・すべり面直上部付近の最大加速度は基盤岩層の最大加速度の1.0～1.9倍の値を、地表付近のそれは１．１～3.6倍の値を示した。
・いずれの地すべりにおいても、周期が概ね1～2秒を超える加速度フーリエ振幅スペクトルは基盤岩層と想定すべり面直上部、地表付近でほぼ同じ値を示した。一方、概ね1～2秒より短い周期では、基盤岩層に比べて想定すべり面より直上部および地表付近で大きな値を示す特徴が見られた。
・譲原地すべりにおける地すべり土塊の水平方向の固有周期は0.22～0.25秒で上下方向のそれは0.10～0.13秒、由比地すべりでは水平方向が0.37秒で上下方向が0.20~0.22秒と推定された。
・固有周期と加速度フーリエ振幅スペクトルの最大値を示す周期は、譲原地すべりでは異なる値を示した。
・由比地すべりでは、2009年の駿河湾沖地震において、深度40mと深度1mの水平方向の加速度フーリエ振幅スペクトルが最大値を示す周期と固有周期は同じ値を示した。
・由比地すべりで地中変位が発生しなかった原因として加速度が地すべりの変位を生じさせるよりも小さかったことが考えられた。

相関性と分離性

sklearn.preprocessing に LabelEncoder があります。

これ、categorical な変数を数値に変換・分類してくれます。sklearn を使わずとも、その前に既に整理されている場合もあるでしょう。馴染みのある分類方法です。

LabelBinarizer や OneHotEncoder もあります。
当初は馴染みのない上、スパースなのにメモリを多く使うため、その利点がわかりませんでした。が、気づくと便利。データの前処理で target との相関分析を行う際、並びを考えなくて良くなります。e>ａ>c>bの順で相関性があった場合、これらで機械的に番号を振ると良い相関を得られません。

機械学習にかけるうえでは、相関性だけでなく、分離性についても着目すべきでしょう。random forest は相関性重視ですが、LightGBMは分離性も見ています。扱う手法が何を重視しているかによって、与えるべきデータ、前処理が異なります。難しい。

2値分類の場合、KDE 等でピークの分離性を見てみます。が、数百もあると目検索では辛い。最大ピークだけ抜き出して距離を測る方法をとっていますが、まだ納得できていません。第2、第3ピークも考慮したいのですが、まだ実装できていません。プロはどうなさっているのでしょう。
https://stackoverflow.com/questions/16255622/peak-of-the-kernel-density-estimation
http://ibis.t.u-tokyo.ac.jp/suzuki/lecture/2015/dataanalysis/L9.pdf

データを理解するためには、ツールと表現方法を実装する必要があります。
もっと多くの知見を得たいものです。

相互相関関数

先輩から「これ、よかったよ」と勧められた図書です。

ディジタル画像処理編集委員会「ディジタル画像処理 [改訂新版]」

全く興味がなかったのですが、お借りして良かったですね。
フィルタ処理や領域処理等、各種画像処理の考え方が包括的に掲載されています。画像を xy 2次元＋RGB 3次元の配列と考えると、線形代数等で扱われている処理がそのまま使えます。回転テンソル、固有値解析など。ベイズ統計やフーリエ変換も出てきました。
画像処理には何らかの数値処理を施しているとは思っていましたが、求める結果を得ることが先で、その中身まで深く考えたことはありませんでした。まさか領域分割の一つにカーネル密度関数とその極大位置が利用されているとは。興味を持たないというのは怖い。

特に興味を惹かれたのが、相互相関関数を利用した画像マッチング。ちょうど、「SPAC係数は自己相関でなく相互相関でない？」などと思案中でしたので、驚き。

残念ながらこの図書には実装できるまでの実用的な記載はありません（ソフトが処理してくれますから必要ないのです）。python だと 2D の相互相関関数もあるかな？と思って調べると、ありました。まさに画像マッチング。

scipy.signal.correlate2d
https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.correlate2d.html
コチラ↓はノーマライズしてます。地球統計学でも標準偏差で割っていた似た式がありました。
https://stackoverflow.com/questions/1819124/image-comparison-algorithm

correlate2d には mode が3種 （full, valid, same）あるようです。どのような結果になるのでしょうか？

このような結果でした↓

a=np.array([[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]])
b=np.array([[1,1,1],[1,1,1],[1,1,1]])
print(a)
print(b)

[[ 1  2  3  4]
 [ 5  6  7  8]
 [ 9 10 11 12]
 [13 14 15 16]]
[[1 1 1]
 [1 1 1]
 [1 1 1]]

cor2 = signal.correlate2d(a, b, mode='full')
print(cor2)

[[ 1  3  6  9  7  4]
 [ 6 14 24 30 22 12]
 [15 33 54 63 45 24]
 [27 57 90 99 69 36]
 [22 46 72 78 54 28]
 [13 27 42 45 31 16]]

cor2 = signal.correlate2d(a, b, mode='valid')
print(cor2)

[[54 63]
 [90 99]]

cor2 = signal.correlate2d(a, b, mode='same')
print(cor2)

[[14 24 30 22]
 [33 54 63 45]
 [57 90 99 69]
 [46 72 78 54]]

左上から順に重ね合わせる作業は3つとも同じ。返すarray の範囲が異なっているだけです。same が使えそうですね。
個人的には、最大値で該当箇所を見つけられる以上に、異なる次元のベクトルを比較して類似度を判定できる性質に、相互相関関数の価値を見出せそうです。

興味がないと割り切ってしまえば終わります。なんでも手を出してみるものですね。
先輩に感謝。

事前確率

ベイズ統計について初心者用の図書を購入しました。

松原望「やさしく知りたい先端科学シリーズ1 ベイズ統計学」

本当に優しく、2日で読み終わりました。で、感想。

「私はベイジアンだったのか」

物事が起こるたびに、期待値を変更します。それが経験則となります。その過程をベイズ更新として肯定されているので、馴染みやすかったですね。

以前、機械学習による胃がんの検出例についての疑問を書き残しています。

https://phreeqc.blogspot.com/2019/05/ai.html

これと同様の内容を、ベイズ統計の考え方で説明されている個所がありました。わかりやすかったですね。事前確率の重要性を尤度と同レベルで述べられています。

この場合、事前確率が非常に小さい、ということを見落としてしまうと、80%の精度の検査で陽性、ということに驚いてしまうかもしれません。メディアなどでもAの事前確率を隠してBの尤度（条件付き確率）だけを提示して印象を操作しようとすることがあります。これからの世の中では事前確率やベイズの定理を意識しながら、冷静に提示された数字や情報、データを読む力が肝要となってくるでしょう。

雨が降って、災害が発生するかしないか。尤度を90%程度に引き上げることは比較的容易です。問題なのは上記の例と同じくらい低い事前確率。尤度を９9％近くまで引き上げないと、「命の危険があります」と繰り返してもオオカミ少年状態。この数%が難しい。
医学分野ではどの程度でスクリーニングに「有用」と判断されているのか知りたいところです。

以前からデータを読む力は重要であったに違いありません。近年では、量・種類共に膨大なデータを加工する能力と共に、その力が必要に迫られているのだと考えます。

Kuzushiji Recognition

Kuzushiji Recognition
https://www.kaggle.com/c/kuzushiji-recognition

先日始まった kaggle のコンペです。日本の人文学オープンデータ共同利用センターなどがホストをされているようです。

このような人文学の研究団体を初めて知りました。いえ、当然あるのでしょうけど、縁がありませんでした。
HPを見てみましたが、興味の引かれるデータを多く公開されています。試してみようか？のみならず、少し読んでみようかと心動かされてしまいますね。
考えてみると、こういった読めなくなった文字を機械学習で認識させる、翻訳するというのは王道です。

土木分野でも、このようなコンペを開催すれば、物によっては比較的安く、良いアルゴリズムを得ることができると思います。プロポよりも効率的でしょう（経産省さんは昨年実施されていました）。

知らない分野でも、知らない間に機械学習の利用が進んでいます。
井の中で目標を立てることなく、外にも目を向けましょう。

地下水モデル その4

10章は不確実性。

私に欠けている知見であり意気込んで読み始めたのですが、正直よくわからないところがありました。（備忘録以上のまとめ量になったので掲載を控えます。）
理解するには手を動かさないとダメなようです。

MCMC、ベイズ統計は、機械学習でも出てきます。最近、本屋でもよく見かけます。データ同化は古くからあるようですが、まだ理解できていません。次はこのあたりの知識を仕入れないといけないでしょうか？

11章は報告書。
この章以外にも、裁判を念頭に置いた地下水モデリングについて触れられています。立ち位置含め、技術者が気を付けないといけないところでしょう。
シビルアクションを見てみましょうか。

以下、備忘録です。
************************************
11.2 モデリングの報告書
11.2.4.1 水文地質環境
・流向を模式的に矢印で示した図はモデルの周囲境界と領域内における一般的な流行を表すために用いることができる。
11.2.4.2 システムの特性
・パラメータの記述（パラメータを求めた手法、推定範囲、空間分布）
11.2.4.2 概念モデル
・画像（模式図・図表）により概念モデルを表現
11.2.5 キャリブレーション

11.4 査読
・査読チェックリスト
・対立する関係者との争点にならない部分では合意を求め、重要な争点を明らかにすることができれば、訴訟は合理化。モデリングはその枠組みを提供する。
・モデリング過程では、利害関係者と意思決定者との意思疎通が重要であるとの意思疎通が高まりつつある。

地下水モデル その3

9章はモデルキャリブレーションです。

4年前の時点で PEST の使用は当たり前、手動の試行錯誤のみはダメ、と断言されています。10章の不確実性にかかわる箇所でもあり、必須なのでしょう。そうなると、使用するソフト制限されます。

ペナルティ項の追加はタンクモデル（SCE-UA）の最適化でも使われていましたね。機械学習でも L1, L2正則化として用いられていました。目的は異なりますが、数式にすると同じです。このあたりは古くからある理論が実務で利用されているということでしょう。

V&Vは数年前に聞き始めましたが、ぼーっとしているうちに、古くなっていたようです。恐ろしい。

あと、この章だけではないですが、予算の配分についても若干の記載があります。日本産の教科書にはない点です。有益です。

全体として地下水にかかわる技術者ならば、心得ておくべき内容が多く記載されていると考えます。あたりですね。

地下水モデル その2

地表流と地下水の連成に興味を惹かれ、6章を読みました。主に MODFLOW で扱われている内容です。
以下、備忘録です。

************************************
6.2 揚水井および注入井
・3つのモデル化
　・ユーザーが層間での揚水・注水量を割り振る：単純すぎて実際の問題にはほとんど適用できない。
　・鉛直透水係数を高くする：許容範囲内の解が得られるが、計算は不安定。乱流・井戸損失は表現できない。
　・MODFLOW井戸パッケージMAW、MNW

BOX6.井戸節点周辺における接点間隔の設定指針
・有効井戸半径re=0.208a→a=4.81～6.66rw

6.3 湧き出し・吸い込み
・フラックス（L/T）で与えるか？堆積流動速度（L^3/T）で与えるか？コードによる。

6.4 排水および湧水
・トンネル排水等：水頭依存境界
・流出のコンダクタンス（L^2/T）をキャリブレーションで推定するのが通常。

6.5 河川
・計算流出量から河川水位を計算
　SFR1パッケージ：マニング式
　SFR2パッケージ：河道特性をセル毎に指定可
・地下水モデルに河川水位の計算を含むと、非線形方程式が加わることで解の安定性に影響。
・解析対象から隔たった地物に対し簡便な河川パッケージ等を適用し、近傍の河川でSFRを利用するというのが通例。

Box6.2
・地表流と不飽和流コードの連成：Hydrogeosphere→時間・空間ステップを共に細かくする必要あり。長期の計算では大きな計算資源が必要。
・計算時間の短縮→GSFLOW:地表水+1次元不飽和浸透+地下水モデル
Box6.3
・地表水：少数地点の河川ハイドロ→時間的に密、空間的に粗なキャリブ
・地下水：時間的に粗、空間的に密なキャリブ
・地下水からの湧出（飽和余剰地表流）の影響大。

地下水モデル その1

先月発売になった「地下水モデル」の気になる箇所を、約1か月かけてコツコツ読んでいました。

地下水モデル―実践的シミュレーションの基礎― 第2版
https://www.kyoritsu-pub.co.jp/bookdetail/9784320047365

和訳初版が1994年ですので、25年振りの改稿です。実際は2015年版の和訳ですので、既に4年前の情報になっています。が、学ぶべきことは多くありました。

第1版との大きな違いは、「キャリブレーション」と「不確実性評価」。この2点が新たに追加されています。後者は私にとってピンポイントで不足、必要に迫られていた内容です。また、前者も十分とは言えません。
実務的には PEST の使用が前提でしょう。PEST は9年前に知り、便利なツールだと感じていましたが、それっきり。時代が進み、何とか実施しないといけないレベルまできていることは、はっきり自覚できました。

以下、1章の備忘録です。
************************************
1.2 モデルとは
・物理モデル：タンク、カラム
・数理モデル
　・データ駆動型：ブラックボックス、ニューラルネット
　・プロセス型：確率論的モデル、決定論的モデル

1.3 モデルの目的
・不確実性が必ず伴うことを強調するため、「予測」よりも「予報」という用語を好んで用いる。
・気象「予報」は確率（例えば、降水確率）によって記述される。

1.4 モデルの限界
・地下水モデルは現実の単純化→近似に制約を受ける。
・実際の複雑な自然を表現仕切ることのできるモデルは存在しない。
・不確実性を評価し、報告しなければならない。

1.5 モデルのバイアス
・依頼主から金銭的授受を受けているモデル作成者が中立を保ちバイアスから逃げられるのか？
・中立性を保ち、プロとしての信頼性を維持することが重要。

1.7 よくあるモデリングの誤り
・プロジェクトの半分の時間と予算をキャリブレーションに当てることを推奨する。
・不確実性解析の実施→予想より長い時間がかかる。モデリングの初期段階に立ち返る必要あり。
・モデル報告に十分な時間を取る。読むに耐え得る包括的な報告書を作成→再検証に必要。

モデルビルダー ArcGIS10.6

1週間ほどArcGISに取り組んでいました。

ポリゴン同士の重なった面積を集計するだけだったのですが、いくつかの条件を満たす必要があり、意外と時間がかかりました（初心者ということが最大の原因ですが）。

その中で知ったのがモデルビルダーと Python。
モデルビルダーで組んだモデルを Python スクリプトとして export できました。QGIS でも Python をサポートしていましたが、Arcでも同様（どちらが先かは知りません）。
繰り返し処理の部分は手直しが必要です、が、形にはなっています。

疑問点はニーズ。どのようなニーズがあるのでしょうか？
モデルビルダーでできることをあらためて Python で組もうとは思わないですし、その他の必要性を感じるほどのヘビーユーザーでもありません。

ヘビーユーザーの使用例を知りたいですね。まだまだ「既知の未知」があるようです。早く必要に迫られるレベルまでたどり着きたいものです。

重力勾配テンソル その2

続きです。

空間周波数fの定義は「回/m」。
対象範囲から最大・最小波長λを出して、最小・最大周波数fに変換すれば良いだけ。先の STACK OVERFLOW に python コードの例が載っていましたが、もっと簡潔に書けます。
周波数がわからなくても、波長から角波数k=2π/λを出せますし、波数だけならフーリエ変換時に０から順に出せるので、計算すら必要ないかもしれません。

曖昧な点は、実装しながら確認することに。

まずは文献に沿ったデータを用意。
次に「新・地震動のスペクトル解析入門」のソースを確認。これでも十分に短いのですが、python では numpy で1行、np.fft.fftn だけでOK。ありがたい。2D だと fft →転置→ fft →転置と同じ結果になったので、この順で内部処理しているだけかもしれません。

実装自体は容易で、2晩ほどでできました。

が、結果が合いません。似たような重力勾配や水平微分の分布になるのですが、オーダーの異なる場合があるのと、ky の符号が逆になっているように見えます。
角波数でダメならただの波数を試したり、ナイキスト周波数以上の領域を全て考慮したり、しなかったり。組み合わせによっては、それっぽい結果になるのですが、数値を記載している文献がないため確証を得られません。3日ほど考えましたが、最終的にはあきらめて寝かすことに。詳細に書かれた文献が、出てくるかもしれません。
もっと簡単に結果を出せると思っていたのですが、残念。

理解にはまだ時間が必要なようです。
ま、重力探査の現状と自分のレベルが分かっただけでも、良しとしましょう。

******************************************
20190813追記
fftn の並びに対応する fftfreq で使う符号、波数のあたりが怪しいようです。修正すると、それらしくなりました。が、検証できないのでここまで。もう少し寝かせましょう。

重力勾配テンソル

講習で知った、重力異常値から重力勾配テンソルを導く方法について、調べていました。

理論はコチラ。
Kevin L. Mickus, Juan Homero Hinojosa
The complete gravity gradient tensor derived from the vertical component of gravity: a Fourier transform technique
https://www.sciencedirect.com/science/article/pii/S0926985101000313

重力異常を重力ポテンシャルの2階微分した重力勾配の形に直しているだけなのですが、フーリエ領域で波数を使った乗算の形にしており、水平成分も算出できるように工夫されています。これを9成分の重力勾配テンソルに整理し逆変換すれば、空間領域の重力勾配テンソルを求められるという流れ。おそらく重力偏差の測定・解釈から思いつかれたのでしょう。

シンプルな流れなのですが、理解に時間を要しました（符号の誤りもありましたし）。
特にわからなかったのが波数の考え方。今回は時間領域への変換ではなく、フーリエ（空間）領域？への変換。ここで躓きました。
2次元平面から波数を求めるイメージはできても、具体的な数値にする方法を複数思い付き、定められません。
調べてみると、同じように悩まれている方がいらっしゃいました。
https://stackoverflow.com/questions/7161417/how-to-calculate-wavenumber-domain-coordinates-from-a-2d-fft

これ、案外メジャーでした。画像処理分野で。
2次元平面での周波数を「空間周波数」と呼ぶそうです。そういえば、画像処理でもハイパスフィルターとかありましたね。画像を扱う際は見た目さえ良くなればOKなので、具体的な周波数には意識を向けていませんでした。が、内部ではまったく同じ計算をしています。うーん繋がる。

続く。

物理探査

先週、物理探査学会の講習会に参加しました。

若い方を対象にしているのかな？と思いながらの参加でしたが、意外とオジ（イ）サマも多く来られていました。
内容は基礎的なことから最新の内容まで。参加するまではもったいないかな？とも考えていましたが、講義についていけないところもあり、個人的には充実した内容の講習会でした。

全体としては、取得データ数の増大、3次元解析の充実が印象に残りました。データの大量取得、大量処理が可能となっている現在、いつまでも2次元にとどまる必要はないのでしょう。3次元解析、当たり前にできると言えるようになりたいものです。

また、処理能力の向上に関しても知らないことが多くありました。多重反射をある程度除去できるようになっているとか、重力探査でソースの走向傾斜を推定できるようになっているとか（これ、ほぼ理解できるようになるまで、3日ほど論文を読み返しました）。

年を取ると「知っている」と勘違いすることが多くなります。気を付けてはいるのですが、なかなか本当に「知らない」ことを自覚するに至りません。今回は良い機会でした。
調べて身につけましょう。

以下、個人的な備忘録です。要チェック！
*****************************************************

・レーダと電磁探査（ハンドブック図9.1）
　・マクスウェル方程式、タンデルタ

・空中電磁探査
　・逆解析で深度決定。
　・表皮深度を使う簡易法では、誤差が5倍程度。

・電気探査
　・4端子法により接地抵抗をオミット。

・SAR
　・運が良ければ2時間ほどで取得可。
　・品質証明が肝要。

・微動アレイ探査
・SPAC係数
　1.分子：クロススペクトルS12
　2.分母：パワスペクトルS11
　3.ρ12=real[S12/(√S11√S22)]
　4.ρ12aveとρ13ave、ρ14aveの算術平均を出す

・位相速度の求め方
　・スパック係数と位相速度を、第一種0次のベッセル関数J0を介して関係づける。
　・J0：アレーサイズから高周波数が決まる。低周波数側は急激に落ち始めるところ。

DETERMINATION OF SOIL SHEAR MODULE AT DEPTHS BY IN-SITU VIBRATORY TECHNIQUES, ARMY ENGINEER WATERWAYS EXPERIMENT STATION VICKSBURG MISS

・重力探査
・半自動解釈手法

• 固有値・固有ベクトルを用いた解析
  ・産総研データ
  　→フーリエ変換
  　→積分→微分→Gz、Gx,Gyがそろう
  　→半自動解釈で断層傾斜角

• 高密度体(基盤岩)の方向に重力偏差テンソルの最大固有ベクトルが向くことを利用し、断層 傾斜角を推定
• 日本でも、地熱地域を中心に、重力偏差探査(重力ポテンシャルの3次元空間微分)が実施されてきている [空中探査]
• インバージョンよりも短時間で重力異常や磁気異常の異常源を推定する半自動解釈手法は、広範囲の大雑把な構造を知る解析に向いている
• JOGMECのHPから申請。データ取得。

コサイン類似度

コサイン類似度が出てきました。

以前、自然言語処理の図書で出てきたように思います。が、無関係でしたので読み飛ばしていました。

今回、あらためて数式を見ると、分子がベクトルの内積。分母がノルム。ノルムで割って単位ベクトルにした後、内積を取っているイメージ。それが cosθ になるからコサイン類似度。方向が同じなら似ている（cosθ=1.）ということでした。

これ、高校か大学1年生レベルの数学ですよね。ただの「2ベクトルのなす角」です（と言いながら、頭に入っていなかったのですが）。

まだまだ基礎力不足。頑張りましょう。

データの理解

最近、データサイエンスの本質は「データを理解すること」にあるように思えてなりません。

機械学習がいくら優れていても、与えるデータの質が悪ければ、それなりの結果にしかなりません。質が悪ければ良いデータに加工して、または質の良いデータを加えて、良いデータのみを選別して、与えるデータを見極める。「前処理」「特徴量エンジニアリング」と呼ばれる過程です。データを解釈する感覚と、多くの処理手法を知り置く必要があります。

これ、地味に難しいですよね。
データを見て、分離性や相関性を見極め、機械学習に与える。Random Forest では予想通りの重要度を示すのですが、LightGBM では真逆の重要度として出てくることも。
最初に多くの機械学習手法にかけて、重要な特徴量を選別する方が良いのかもしれません。プロはどうされているのでしょう？

機械学習に限らないのですが、「データをよく見て理解する」というのは重要なのに難しい。機械学習では、それについても一部で自動化が進みつつあるようですが、まだまだ人の感覚・判断も重要でしょう。感覚を掴み、磨きたいですね。

Print Word & EXCEL Files

この時期に珍しく、後輩君が遅くまで頑張っていました。

どうも、複数のフォルダに分かれた大量の Word, EXCEL ファイルの印刷を命じられたようです。

これは機械の得意分野。機械にまかせるべき作業です。

方法はいくつかあるでしょう。
今回はテスト感覚で Pythonを選択。COM 等を使用したコードはネットに転がっていましたので、それを利用。

動かしてみると、意外と速い。
昔、EXCEL＋VBA で同じようなマクロを組んだことがありましたが、それよりも速いように感じました。

10年後、Python が使い続けられているかは疑問です。が、今、何らかの言語＋アルゴリズムを身に着けることは技術者にとって必須です（変化に柔軟に対応することは、それ以上に重要です）。
常に手を動かし、今も、将来も困らないようなレベルに身を置くことができるよう、行動しましょう。

再現期間、確率降水量の計算

カーネル法では、カーネル×重みを足し合わせて既知の値にフィッティングしました。

これは、予想される形状が複雑な場合に有効です。
では、予想される形状が単純な場合、一つのカーネルを調整するだけでフィッティングできるかもしれません。コレ、曲線近似、曲線回帰に属するようです。

数種類のカーネル関数で試行し、もっとも適合したものを選択する、というのは常套手段。例えば、コチラ↓。

国交省気象庁「確率降水量の推定方法」
https://www.data.jma.go.jp/cpdinfo/riskmap/cal_qt.html

(1) グンベル分布
(2) 一般化極値（GEV）分布
(3) 平方根指数型最大値分布
(4) 対数ピアソンⅢ型分布
(5) 対数正規分布

これらの関数を降雨量にあてはめ、もっともフィットしたものを選択する、というのが大まかな計算の流れ。
後半のSLSCというのが経験的というのは、今回初めて知りました。情けないことですが、ツール任せの部分であり、残差等で比較しているものと勘違いしていました。しっかり読まなくては。

実務で何度か計算していますが、回数は片手で足りるほど。上記の反省 & 理解を深めるために、前半の計算を追ってみました。使ったのは Python。

計算や図化は簡単だったのですが、水文統計ユーティリティーでの結果と差異が生じます。試算ではヒストグラムを作る時のビン数によりフィッティング結果が変わるのですが、ユーティリティーではどのビン数で計算したのかが不明（後のヒストグラム表示でビン数を変えても分布は変化しないので、内部である程度決めた値を使っているのかもしれません）。フィッティング後の係数も不明な点は、以前から疑問視していた通り。

また、確率紙についてもこれまで真剣に見ておらず、対数軸と誤認していました。反省。
計算を追う過程で先輩の古い（50年以上前の）教科書をお借りし、ようやく理解しました（Rには確率紙のライブラリがあったのですが、Python にはなかったので計算・図化しました）。
が、この図で経験的にフィットの良し悪しを判別する必要があるのでしょうか？
ヒストグラムを作った段階で何らかの最適化手法によるフィッティングをかけているので、残差の小さいカーネル関数を選択すれば良いのではないでしょうか？不都合があるのでしょうか？
確立紙プロットは確認程度でも良いですよね。確かに、50年以上前だと確率紙にプロットするしかなかったのかもしれません。が、今は容易に計算できます。

得られた知識はあったものの、新たな疑問も増えました。
解決するのは、まだ先のようです。

***********************************
20190626追記
河川砂防技術基準に以下の内容があります。

解析対象水文資料を用いて候補モデルの母数を求める際には、標本の大きさに応じて適切な推定法を用いるなどの手法があり、積率法、L 積率法、最尤法等の手法が用いられている。なお、小標本（標本サイズ＜30）については、L 積率法がよく用いられている。 

なぜか Win ＋ Anaconda が動かないので追認できませんが、おそらく、これですね。

***********************************
20190627追記
L積率法でパラメーターを固定すると、ユーティティーの値に近づきました。
ただ、先輩たちも昔からこの手の検算はされてきたようで「ユーティリティと答えが一致しない」というのは知られていたようです。いくつか理由が考えられるそうですが。
検算するなら、水理公式集が良さそうです。

カーネル法とニューラルネット

「カーネル密度推定」と「カーネル法」。言葉は似ていますが異なります。
\begin{align*}\hat{f}(x) = \sum_{i=1}^n w_i K(x_i,x)\end{align*}
なぜ重みwを付けたのか？
フィッテイングしたいからです。既知の値に。

この辺が分かりやすく、詳しいです。↓
https://qiita.com/wsuzume/items/09a59036c8944fd563ff

既知の値と計算値との残差＋正則化項を最小にする重みを計算し、決定します。ニューラルネットに似ています。

ニューラルネットでカーネルに2次元のガウス関数を用いている例が身近にありますね。これです。
国土交通省河川局砂防部「国土交通省河川局砂防部と気象庁予報部の連携による土砂災害警戒避難基準雨量の設定手法（案）」平成１７年６月
https://www.mlit.go.jp/river/shishin_guideline/sabo/dsk_tebiki_h1706.pdf

統計、避けては通れないようです。

カーネル密度推定（KDE）多次元

p次元のカーネル密度推定式
https://web.as.uky.edu/statistics/users/pbreheny/621/F10/notes/10-28.pdf
\begin{align*}\hat{f}(x) = \frac{1}{n} \sum_{i=1}^n\prod_{j=1}^p \frac{1}{h_j}K\bigg(\frac{x_j-X_{ij}}{h_j}\bigg)\end{align*}
Pythonだと、例えばこれだけ。便利。

sns.kdeplot(x, y, shade=True, shade_lowest=False)

https://seaborn.pydata.org/generated/seaborn.kdeplot.html

カーネル密度推定（KDE）直感的理解

KDEの直感的な解釈です。
厳密でなくとも、概略イメージを先に持つことができると、後に数式がすっと頭に入りますし忘れ難いと思います。以下、個人的な解釈です。

「家の裏山が崩れた」時、ある観測値が0.5だったとします。

観測値0.5であれば、崩れる確率が高いと事後解釈できるかもしれません。が、0.4だとどうか？0.2だとどうか？わかりません。実は、1.0の方が崩壊しやすいのかもしれません。

ひとまずデータが集まるまでは、0.5周辺で崩壊確率がかなり高く、離れるほど低くなると捉えましょう。たとえば、以下のように。

この形は「正規分布（ガウス分布）」。
山のような形でなく、どのような形（三角でも矩形でも）を選んで良いのですが、モデル化では連続（なめらか）かつ数式（関数）として扱える方が便利。
縦軸は確率密度と呼ばれています。横軸の特定範囲でたせば（面積）、その範囲で発生する確率になります（これが確率密度関数の定義。全て足せば1.0）。なお、バンド幅h=0.05を考慮した形にしていますが、そこはひとまず気にしないことにします。

しばらくすると、また崩壊が発生。この時の観測値が0.06。
この時、崩壊を導く値（ピーク）が２つあるのでは、と考えた場合、こうなります。

次の観測値が0.45、0.95。
重なっているところがあるので、単純に足してデータ数４（とバンド幅0.05）で割って、ならしておきましょう。破線で表示される密度分布で、3つのピークができました。

観測を続けて足し合わせると、このような形になりました。0.05、0.5、0.9付近で崩壊する確率が高いと判断できます。

このような尤もらしい破線を推定する作業をカーネル密度推定（KDE）と呼びます。
上記の場合、カーネルとしてガウス分布を採用し、不明な母集団の関数を推定したという作業になります（正しく推定できているかどうかは別問題）。

そういえば、先日読んでいたベイズ統計の教科書にも似たような記述がありました。
カーネルを使った推定手法は、他にもいくつかあるようです。それらは理解できていませんが、必要になった際に覚えましょう。

カーネル密度推定（KDE）数式

カーネル密度推定を試行。

手元に数式がなかったので、ひとまずネットで検索。数式は微妙に違いましたが、だいたいこんな感じ（それで良いのか？）

一次元のカーネル密度推定式
\begin{align*}\hat{f}(x) = \frac{1}{nh} \sum_{i=1}^n K\bigg(\frac{x-x_i}{h}\bigg)\end{align*}

n：データ数
h：バンド幅
x：計算位置
xi：i番目のデータ（例えばi番目の観測値）
K()：カーネル関数

x位置でのK()を足してnhで割っています。

カーネル関数K()の例：正規分布（ガウス分布）の確率密度関数
\begin{align*}\ K(\nu)=\frac{1}{\sqrt{2\pi \sigma^2}}\exp{\left\{-\frac{(\nu-\mu)^2}{2\sigma^2}\right\}}
\end{align*}

EXCELだと、手間。面倒。

そもそも、EXCELでは約105万行までのデータしか扱えません。1000万データなど、開くこともできません。

Pythonだと、簡単。

sns.kdeplot(X,bw=0.05)

seabornだけでなく、他にも複数のライブラリが KDE に対応しています。
ひとまず、700万行×13列のデータは一気に処理できました。便利。

エラー修復 (GEORAMA 2019)

CTC さんの GEORAMA 2019を使用している際に遭遇する問題です。

ユーザー側で対処に迫られる（かなり厳しい）事項のみ、経験的対処法を備忘録として残しておきます（サポートさんに伺っても、GEORAMA が原因ではないという見解のようです）。

Q1:図面を存できません。
コマンド: _QSAVE Object(s) added to DB too late during save. Handles 32C1AD to 32C1AF
32C1AD: type AcDbBlockTableRecord, owner 1 (AcDbBlockTable)
32C1AE: type AcDbBlockBegin, owner 32C1AD (AcDbBlockTableRecord)
32C1AF: type AcDbBlockEnd, owner 32C1AD (AcDbBlockTableRecord)

A1: GEORAMA で境界関連図を再作画しましょう。
保存できるようになります。
昔はこのエラーはなかったのですが、最近よく遭遇します。


Q2：dwgを開いても、GEORAMAパネルにデータが表示されません。

A2：GEORAMAのデータが（なぜか）飛んでしまいます。2019からではなく、過去のバージョンからシバシバ遭遇します。
対処法1：悲しいですが「GEORAMAデータはいつか失われるもの」という認識で、バックアップをこまめに取るというのが経験的な対処法。
対処法２：設定のやり直し。
直近のバックアップがない場合、作り直すしか選択肢は残っていません。
が、初期データ登録から順番に追っていくと、地質テーブル、境界テーブルは容易に復活します。ボーリングデータリストは「ボーリングデータチェック」で認識してくれます。
ここまでで一度保存し、ファイルを閉じて再度図面を開くと、「鉛直断面図」が復活します。残念ながら境界ポイントは認識してくれません。一度削除し、再設定が必要です。

機械学習とCAE

機械学習とCAEの融合が進んでいるようです。
http://www.cae21.org/kansai_cae/konwakai/kansai_67/kansai67_annai.shtml

表題をみても、まったくわかりません。
想像もできません。知りたいのですが、雲の上。

うーん。ついて行きたい。