2021年11月21日日曜日

差分近似の高精度化

微分の差分近似は、テイラー展開後に目的の項について整理するのみです。

簡単に高次精度化ができるようですが、今まで気にしたことがありませんでした。
反省も込めて、整理です。

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

1次精度‐1階微分
前進差分(テイラー展開)
\begin{align*}f\left(x+\mathrm{\Delta x}\right)=f\left(x\right)+\frac{\mathrm{\Delta x}f^\prime\left(x\right)}{1!}+\frac{{\mathrm{\Delta x}}^2f^{\prime\prime}\left(x\right)}{2!}+\frac{{\mathrm{\Delta x}}^3f^{\prime\prime\prime}\left(x\right)}{3!}+\frac{{\mathrm{\Delta x}}^4f^{(4)}\left(x\right)}{4!}+\frac{{\mathrm{\Delta x}}^5f^{(5)}\left(x\right)}{5!}+\ldots\end{align*} \begin{align*}f^\prime\left(x\right)=\frac{f\left(x+\mathrm{\Delta x}\right)-f\left(x\right)}{\mathrm{\Delta x}}-\frac{{\mathrm{\Delta xf}}^{\prime\prime}\left(x\right)}{2!}-\ldots\end{align*} \begin{align*}=\frac{f\left(x+\mathrm{\Delta x}\right)-f\left(x\right)}{\mathrm{\Delta x}}+O\left(\mathrm{\Delta x}\right)\end{align*}

 後退差分(テイラー展開)
\begin{align*}f\left(x-\mathrm{\Delta x}\right)=f\left(x\right)+\frac{-\mathrm{\Delta x}f^\prime\left(x\right)}{1!}+\frac{{\mathrm{\Delta x}}^2f^{\prime\prime}\left(x\right)}{2!}+\frac{-{\mathrm{\Delta x}}^3f^{\prime\prime\prime}\left(x\right)}{3!}+\frac{{\mathrm{\Delta x}}^4f^{(4)}\left(x\right)}{4!}+\frac{{-\mathrm{\Delta x}}^5f^{(5)}\left(x\right)}{5!}+\ldots\end{align*} \begin{align*}f^\prime\left(x\right)=\frac{f\left(x\right)-f\left(x-\mathrm{\Delta x}\right)}{\mathrm{\Delta x}}+O\left(\mathrm{\Delta x}\right)\end{align*} 

2次精度‐1階微分
中心差分(前進差分‐後退差分)
\begin{align*}f\left(x+\mathrm{\Delta x}\right)-f\left(x-\mathrm{\Delta x}\right)=\frac{2\mathrm{\Delta x}f^\prime\left(x\right)}{1!}+\frac{2{\mathrm{\Delta x}}^3f^{\prime\prime\prime}\left(x\right)}{3!}+\frac{{2\mathrm{\Delta x}}^5f^{(5)}\left(x\right)}{5!}+\ldots\end{align*} \begin{align*}f^\prime\left(x\right)=\frac{f\left(x+\mathrm{\Delta x}\right)-f\left(x-\mathrm{\Delta x}\right)}{2\mathrm{\Delta x}}-\frac{{\mathrm{\Delta x}}^2f^{\prime\prime\prime}\left(x\right)}{3!}-\ldots\end{align*} \begin{align*}=\frac{f\left(x+\mathrm{\Delta x}\right)-f\left(x-\mathrm{\Delta x}\right)}{2\mathrm{\Delta x}}+O\left({\mathrm{\Delta x}}^2\right)\end{align*}

2次精度‐2階微分
中心差分(前進差分+後退差分)
\begin{align*}f\left(x+\mathrm{\Delta x}\right)+f\left(x-\mathrm{\Delta x}\right)=2f\left(x\right)+\frac{{2\mathrm{\Delta x}}^2f^{\prime\prime}\left(x\right)}{2!}+\frac{2{\mathrm{\Delta x}}^4f^{(4)}\left(x\right)}{4!}+\frac{2{\mathrm{\Delta x}}^6f^{(6)}\left(x\right)}{6!}+\ldots\end{align*} \begin{align*}f^{\prime\prime}\left(x\right)=\frac{f\left(x+\mathrm{\Delta x}\right)-2f\left(x\right)+f\left(x-\mathrm{\Delta x}\right)}{{\mathrm{\Delta x}}^2}-\frac{{2\mathrm{\Delta x}}^2f^{(4)}\left(x\right)}{4!}-\ldots\end{align*} \begin{align*}=\frac{f\left(x+\mathrm{\Delta x}\right)-2f\left(x\right)+f\left(x-\mathrm{\Delta x}\right)}{{\mathrm{\Delta x}}^2}+O\left({\mathrm{\Delta x}}^2\right)\end{align*} 

2次精度‐1階微分
中心差分(2Δx)
\begin{align*}f\left(x+2\mathrm{\Delta x}\right)=f\left(x\right)+\frac{2\mathrm{\Delta x}f^\prime\left(x\right)}{1!}+\frac{2^2{\mathrm{\Delta x}}^2f^{\prime\prime}\left(x\right)}{2!}+\frac{2^3{\mathrm{\Delta x}}^3f^{\prime\prime\prime}\left(x\right)}{3!}+\frac{2^4{\mathrm{\Delta x}}^4f^{(4)}\left(x\right)}{4!}+\frac{2^5{\mathrm{\Delta x}}^5f^{(5)}\left(x\right)}{5!}+\ldots\end{align*} \begin{align*}f\left(x-2\mathrm{\Delta x}\right)=f\left(x\right)+\frac{-2\mathrm{\Delta x}f^\prime\left(x\right)}{1!}+\frac{2^2{\mathrm{\Delta x}}^2f^{\prime\prime}\left(x\right)}{2!}+\frac{-2^3{\mathrm{\Delta x}}^3f^{\prime\prime\prime}\left(x\right)}{3!}+\frac{{2^4\mathrm{\Delta x}}^4f^{(4)}\left(x\right)}{4!}+\frac{{-2}^5{\mathrm{\Delta x}}^5f^{(5)}\left(x\right)}{5!}+\ldots\end{align*} \begin{align*}f\left(x+2\mathrm{\Delta x}\right)-f\left(x-2\mathrm{\Delta x}\right)=\frac{2^2\mathrm{\Delta x}f^\prime\left(x\right)}{1!}+\frac{2^4{\mathrm{\Delta x}}^3f^{\prime\prime\prime}\left(x\right)}{3!}+\frac{2^6{\mathrm{\Delta x}}^5f^{(5)}\left(x\right)}{5!}+\ldots\end{align*} \begin{align*}f^\prime\left(x\right)=\frac{f\left(x+2\mathrm{\Delta x}\right)-f\left(x-2\mathrm{\Delta x}\right)}{4\mathrm{\Delta x}}-\frac{{2^2\mathrm{\Delta x}}^2f^{\prime\prime\prime}\left(x\right)}{3!}-\ldots\end{align*} \begin{align*}=\frac{f\left(x+2\mathrm{\Delta x}\right)-f\left(x-2\mathrm{\Delta x}\right)}{4\mathrm{\Delta x}}+O\left({\mathrm{\Delta x}}^2\right)\end{align*} 

4次精度‐1階微分
中心差分(Δx, 2Δx)
\begin{align*}f\left(x+\mathrm{\Delta x}\right)-f\left(x-\mathrm{\Delta x}\right)=\frac{2\mathrm{\Delta x}f^\prime\left(x\right)}{1!}+\frac{2{\mathrm{\Delta x}}^3f^{\prime\prime\prime}\left(x\right)}{3!}+\frac{2{\mathrm{\Delta x}}^5f^{(5)}\left(x\right)}{5!}+\ldots\end{align*} \begin{align*}f\left(x+2\mathrm{\Delta x}\right)-f\left(x-2\mathrm{\Delta x}\right)=\frac{2^2\mathrm{\Delta x}f^\prime\left(x\right)}{1!}+\frac{2^4{\mathrm{\Delta x}}^3f^{\prime\prime\prime}\left(x\right)}{3!}+\frac{2^6{\mathrm{\Delta x}}^5f^{(5)}\left(x\right)}{5!}+\ldots\end{align*} \begin{align*}8\left\{f\left(x+\mathrm{\Delta x}\right)-f\left(x-\mathrm{\Delta x}\right)\right\}-\left\{f\left(x+2\mathrm{\Delta x}\right)-\ f\left(x-2\mathrm{\Delta x}\right)\right\}=12\mathrm{\Delta x}f^\prime\left(x\right)+O\left({\mathrm{\Delta x}}^5\right)\end{align*} \begin{align*}f^\prime\left(x\right)=\frac{8\left\{f\left(x+\mathrm{\Delta x}\right)-f\left(x-\mathrm{\Delta x}\right)\right\}-\left\{f\left(x+2\mathrm{\Delta x}\right)-\ f\left(x-2\mathrm{\Delta x}\right)\right\}}{12\mathrm{\Delta x}}+O\left({\mathrm{\Delta x}}^4\right)\end{align*} 

4次精度‐2階微分
中心差分(Δx, 2Δx)
\begin{align*}f\left(x+\mathrm{\Delta x}\right)+f\left(x-\mathrm{\Delta x}\right)=2f\left(x\right)+\frac{{2\mathrm{\Delta x}}^2f^{\prime\prime}\left(x\right)}{2!}+\frac{2{\mathrm{\Delta x}}^4f^{(4)}\left(x\right)}{4!}+\frac{2{\mathrm{\Delta x}}^6f^{(6)}\left(x\right)}{6!}+\ldots\end{align*} \begin{align*}f\left(x+2\mathrm{\Delta x}\right)+f\left(x-2\mathrm{\Delta x}\right)=2f\left(x\right)+\frac{2^3{\mathrm{\Delta x}}^2f^{\prime\prime}\left(x\right)}{2!}+\frac{2^5{\mathrm{\Delta x}}^4f^{(4)}\left(x\right)}{4!}+\frac{2^7{\mathrm{\Delta x}}^6f^{(6)}\left(x\right)}{6!}+\ldots\end{align*} \begin{align*}16\left\{f\left(x+\mathrm{\Delta x}\right)+f\left(x-\mathrm{\Delta x}\right)\right\}-\left\{f\left(x+2\mathrm{\Delta x}\right)+f\left(x-2\mathrm{\Delta x}\right)\right\}=30f\left(x\right)+\frac{24\mathrm{\Delta x}f^{\prime\prime}\left(x\right)}{2!}+O\left({\mathrm{\Delta x}}^6\right)\end{align*} \begin{align*}f^{\prime\prime}\left(x\right)=\frac{16\left\{f\left(x+\mathrm{\Delta x}\right)+f\left(x-\mathrm{\Delta x}\right)\right\}-30f\left(x\right)-\left\{f\left(x+2\mathrm{\Delta x}\right)+f\left(x-2\mathrm{\Delta x}\right)\right\}}{12{\mathrm{\Delta x}}^2}+O\left({\mathrm{\Delta x}}^4\right)\end{align*} 

 

2021年11月19日金曜日

LSTM の x_train データ

Keras + TensorFlow で LSTM と Seq2Seq。

どちらも与えるデータの形は3次元。
[num_samples, num_times, num_features]

画像を与えるのと同じ感覚ですね。

Sony NNC では num_samples 個の2Dデータとなり、保存するのも削除するのも一苦労でした。が、2次元データを加工して3次元にするのであればハンドリングが容易。

ありがたい。

地下水観測手法と空間スケール

もう一つ、惹かれた文献。

Monitoring ground water storage at mesoscale using seismic noise 30 years of continuous observation and thermo-elastic and hydrological modeling
https://www.nature.com/articles/s41598-017-14468-9

地下水観測手法と対応する空間スケール

small
・traditional monitoring (piezometer networks)
medium (kilometric to decakilometric scale)
・continuous seismic monitoring (ambient seismic noise)
large (>200 000km²) or subsiding areas
・the GRACE gravimetric mission
・InSAR methods

内容

  • Gräfenberg における30年間の地震観測結果(ambient seismic noise)を利用し、地下の力学的変化を検討。
  • 1.2~10秒の周波数帯における速度変化を処理することで、年スケールでのレーリー波速度変化を0.01%オーダーでとらえた。
  • 観測された速度変化は、局所的な帯水層(静水圧)、または地下全体での変化(熱弾性)によって生成される。
  • dv/vの長期的な変化について熱弾性と静水圧のみを考慮するという限定的な近似条件の下で、両者の相対的な貢献度を推定。熱弾性効果を50%、静水圧効果を50%考慮した場合、観測された dv/v との間に最も高い相関係数(0.83)が得られた(Fig.5A)。
  • 観測局グループ間の速度変化(Fig.3)は、水文学的な挙動の空間的なばらつきによるものであり、熱弾性的な影響はより均質であると予想される。
  • 地質学的に安定しているこの地域では、熱弾性と水文学の両方の影響を、年単位の時間スケールと㎞単位の空間スケールでフォワードモデル化可能。
  • 将来的には、深さ方向の分解能向上のため異なる周波数を用いて地下の4次元画像を解析することや、高密度地震アレイを使用して横方向の変動をとらえることが期待される。


GRACEは知りませんでした。
SAR の利用は最近聞いたり読んだりしました。海外では既に浸透しているのかもしれません。国内では山林が多いため、都市部を除き適用が難しい手法でしょうか。
微動を用いた観測手法は国内でも使えそう。というか、世界に誇るべき高密度の観測網が整備されていますので、積極的に研究に利用すべきでしょうね。F-net や Hi‐net のデータを取り入れて速度変化を表示・保存するシステム。技術的にハードルは高くないでしょう。サブダクションの影響がどこまで出てくるか?ん?影響が大きければ、そちら方面に使えるかな?

無償で公開されている地震波形や崩壊データは、近年、海外の研究者にも利用されています。日本の研究者が牽引する形で頑張っていただきたいところです。

 

2021年11月18日木曜日

東北地方太平洋沖地震後の速度低下

ambient noise を利用した cross correlation にて観測局間の表面波速度を求めている事例は海外の文献で多く見かけます。国内でもたまに見ますが、それらはすべて広帯域の低周波を利用しています。

昨日、はじめて Hi-net を利用した文献を見ました。
しかも、内容が興味深い。

Mapping pressurized volcanic fluids from induced crustal seismic velocity drops (science.org)

  • 2011東北地方太平洋沖地震前後で地震波速度変化が認められた(Hi‐net、cross correlation で半年前と比較)。
  • 火山地帯での低下率が大きい。
  • 動的ひずみと動的応力を、PGV(KiK‐net)、レイリー波の平均位相速度(〜3km/s)、平均せん断弾性率(~30×109Pa)より求めた。
  • 地震波による動的ひずみは、static coseismic strain (静的地震時ひずみ?)よりも1~2桁大きい。
  • 速度低下率と動的応力の比を seismic velocity susceptibility (地震速度の感度)と定義。この分布には主要な火山地帯の分布との相関性が認められた。
  • 地震速度の感度は、火山地帯における熱水やマグマ流体の加圧レベルの指標として使用できる 。
  • 火山地帯、特に(静的応力の変化を無視できる)富士山麓での大きな速度低下の主な原因は、東北地方太平洋沖地震による動的応力であると結論づけられる。

使用しているのは Hi-net のデータだけなので、オンタイムでデータを受信し、日々の snapshot を残しておけば、地震の発生毎に全国レベルで速度変化を見比べることが可能かつ容易になります。水位変化や応力変化の把握にも使えるかな?

海外産ツールは多く出ていますが、トモグラフィーまで完走できていません。
https://phreeqc.blogspot.com/2021/05/passive.html
動かせるようにしておかないと。

********** 関連 *********
Static strain and stress changes in eastern Japan due to the 2011 off the Pacific coast of Tohoku Earthquake, as derived from GPS data
https://earth-planets-space.springeropen.com/articles/10.5047/eps.2011.06.049


2021年11月12日金曜日

Rockfall radar

正に ’cool’ という単語が当てはまるシステムです。
https://www.youtube.com/watch?v=o-oVXYkBwgw

使用されているのはスイスの GEOPREVENT 社の商品のようです。
Rockfall radar 
https://www.geopraevent.ch/technologies/rockfall-radar/?lang=en
・all-weather suitability (works in rain, fog, snowfall)
・works day and night

他にも魅力的な商品が多くあります。
Avalanche radar
https://www.geopraevent.ch/technologies/avalanche-radar/?lang=en
これは落石とアルゴリズムを変えているだけでしょうか。

GEORADAR
地上SARは国内でも実績ありですね。親会社?の商品でしょうか。

Deformation camera
https://www.geopraevent.ch/technologies/deformation-camera/?lang=en
カメラ画像から1日の動きをcmオーダーで判定。これは昼のみでしょう。
画像を用いた異常検知では機械学習を利用した文献を見かけます。学習してデプロイしてしまえばオンタイムで判定可能です。精度向上とともに追っかけの商品が出てくるでしょうね。

 
国内のシステムは子供のように見えます。Rockfall radarなど、一度は使ってみたい商品です。

2021年11月7日日曜日

メンテナンス

今週末は計算を回りっぱなしにしておいて、細々したメンテナンスを片づける。

・Microsoft Natural Ergonomic Keyboard 4000 パームレスト張替え
・自転車1のライト、ブレーキシュー、ブレーキワイヤー、Vブレーキパーツ手配
・自転車2のキズ消し
・自動車のフロントガラス 鱗取り

本皮のグレーの端切れを購入し、型を取って切ったところで表裏逆に気付く。ま、良いかとそのまま続行。なかなかの出来栄え。痛んだら、次は茶系にしましょう。

何を消そうと考えていたのか思い出せないほど前に購入していたソフト99「ツヤ仕上げキズ消しセット」。先週、道具を散らかしていたら出てきました。保管していた折り畳み自転車(白)に、拭くだけでは取れない黒いスリ汚れが多数あったことを思い出し、使用。初めて使いましたが簡単でよく落ちました。

来週末は別の自転車のブレーキ回り交換。シフト周りは調整のみでまだ持つかな。そうこうしていると冬タイヤが交換の季節。そのうち、大掃除、年末へ。
頑張ろう。

地表の地震動予測法

技法堂出版「入門数理地震工学」より

地表の地震動予測法は4つに大別される。

(1)経験的方法

  • 観測地震動の最大加速度値等と地震規模、震源距離、地盤条件等の統計処理により地震動の特性値を予測
  • 最大値等の震動特性値を簡単に評価できる
  • 波形は得られない
  • 水平多層弾性体にのみ適用可
  • 逆解析による断層・地層構造の推定可


(2)半経験的方法(経験的・統計的グリーン関数法)

  • 経験的:余震記録が伝播経路・サイト特性を表すと考え、観測記録を用いて大地震の波形を予測
  • 統計的:観測波形の代わりに合成波形を利用
  • 高精度
  • 短周期(2Hz~)


(3)理論的方法(FDM, FEM, 剛性行列法)

  • 震源・伝播経路・サイト特性の全てをモデル化し、波動方程式から波形を予測
  • 震源断層モデル:運動学的断層・動力学的断層
  • 地盤モデル:水平多層弾性体・不正形多層弾性体・不均質不正形多層弾性体
  • 長周期(~2Hz)


(4)ハイブリッド法(半経験的方法+理論的方法)

  • 短周期帯域を統計的グリーン関数法+長周期帯域を理論的方法(例えば剛性行列法)を用いて波形を予測
  • 短周期~長周期をカバー
 

2021年11月5日金曜日

Ubuntu dist-upgrade

Ubuntu は手元に3台あるのですが、2台は 20.04、1台は18.04。

1台だけ環境が異なると支障の出る場合があり、最後の1台も上げましょうと重い腰を上げました。
「Software Updater」で 「upgrade」 をポチ、で、できるはずだったのですがボタンを押しても無反応。

$ sudo apt dist-upgrade
0 upgraded, 0 newly installed, 0 to remove and 1 not upgraded.
cuda-driver が引っ掛かっていました。

$ sudo apt install cuda-drivers
はい、確かにインストール不可→対処↓
https://forums.developer.nvidia.com/t/cuda-install-unmet-dependencies-cuda-depends-cuda-10-0-10-0-130-but-it-is-not-going-to-be-installed/66488/16
$ sudo apt clean
$ sudo apt update
$ sudo apt purge cuda-drivers
$ sudo apt purge nvidia-*
$ sudo apt autoremove
$ sudo apt install cuda-drivers

nvidia-* は迷いましたが、また入れたら良いでしょうと考え enter⏎。
これで完了。

最後に「Software Updater」 で upgrade ボタンを押すと反応しました。


2021年11月4日木曜日

Amplitude of the FTF

lsforce で距離を設定すると土量が少なく推定される原因2つ
https://pubs.geoscienceworld.org/ssa/srl/article-abstract/92/4/2610/596477/lsforce-A-Python-Based-Single-Force-Seismic?redirectedFrom=fulltext

The amplitudes of the FTF produced by lsforce will be lower than their true magnitudes for two principal reasons: the bandlimited nature of the input waveforms and the effects of the Tikhonov regularization scheme.

前者のみでは大きいなあと感じていたところです。
正則化でもそれぞれ影響があるといわれると、確かにそうかもしれないという程度の理解ですが、言われると納得できます。