不飽和透水試験

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

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

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

ま、プレ浸透を長くすれば飽和透水係数を得られるわけでもありません。が、雑に扱うのは論外。短時間のプレ浸透で飽和とみなせる理屈はないでしょう。
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 のデータも読めました。
が、樹木等を除去するのはかなり困難。プロもトライ＆エラーで解決するそうです。航測会社さんの処理はとても速いのですが、きっと独自のノウハウがあるのでしょう。
これはあきらめました。

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

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

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