2021年5月23日日曜日

緊急事態宣言とバイアス

緊急事態宣言が各地で出ています。

個人レベルで対策可能なことは皆さん以前から実施されているでしょう。そういう点では宣言が出ても日常生活に変化がないかもしれません。

勤め先では、リモートワークに対応する支社と対面を続ける支社で分かれています。top executive の 考え方の違いです(ほとんどが前者ですが)。 

「マスクは意味がない」と言われていた方も、心境の変化かマスクをつけられる機会が多くなりました。ほとんどの方が常時マスクを外さないからでしょう。

一方、味覚異常、発熱、咳などの症状を有する方が PCR 陰性ということで休まず出社されていたり、それを静観する executives であったり。ここにきて正常性バイアスを振り払うのは大変なようです。

リスク管理や危機管理に対する個人の考え方や能力が露呈しています。よく見ておきましょう。


地下構造のモニタリング

波動を利用した地下構造のモニタリングにかかわる総論的資料2件。導入部分主体ですが基礎知識として有用な個所を ( ..)φ。

中原恒「地震波干渉法 その 2 応用」地震波干渉法による地下構造のモニタリング

実際には常時微動源の分布が等方的ではないことが多く,その場合は常時微動の相互相関はグリーン関数に対応しない.しかし,その場合でも,常時微動源が時間的に安定しているならば,常時微動の相互相関を用いて地下構造のモニタリングを行うことが可能であることが指摘されている[Hadziioannouetal.(2009)].
波形から地下構造の変化を検出する手法は,すでに相似地震の解析のために提案されていたコーダ波干渉法[Poupinetetal. (1984),Sniederetal. (2002),Snieder (2006)]がよく使われる.
常時微動だけではなく,自然地震を用いることも可能であり,こちらの場合は常時というわけにはいかないが,自然地震が発生するたびにイメージングが可能となり,様々な時期の波動場の相互相関関数や自己相関関数を比較することにより,その期間内に地下構造の変化があったのかどうかを監視できる.
大・中地震の震源域では,2000年鳥取県西部地震に伴う変化を相互相関関数により検出したSawazakietal. (2009) や2011年東北地方太平洋沖地震に伴う変化を自己相関関数により検出したNakahara (2015) がある.

澤崎郁「大地震に伴う地下構造変化とその回復過程について」

地下構造変化に関する研究は,高密度かつ高精度な地震観測網の整備[例えば,Okadaetal. (2004)]と,コーダ波干渉法や地震波干渉法に代表される解析手法の高度化[例えば,Sniederetal. (2002),Curtisetal. (2006)]に伴い,この十数年の間に飛躍的な進歩を遂げた.
大地震に伴う地震波速度の変化は,古くは断層運動に伴う断層面近傍での静的応力変化を反映するものとして解釈されることが多かった[例えば,早川(1951),Poupinetetal. (1984),Lietal. (1998)]が,近年は強震動に伴う地盤浅部の損傷が速度変化の主要因として有力視されている[例えば,Rubinstein andBeroza (2005),PengandBen-Zion(2006),Sawazakietal. (2006),Wuetal. (2009a),Takagietal.(2012)].
強震動により地盤の物性が変化して弱震時とは異なる地盤増幅特性を示す現象は,地盤の非線形応答として知られているが,強震動が終息した後も地盤の物性が元に戻らない場合がある[例えば,AguirreandIrikura(1997),PavlenkoandIrikura(2002)].
強震動によりいったん低下した地震波速度は,数か月から数年をかけて経過時間の対数に比例しながら回復する[例えば,Rubinstein andBeroza (2004a),Sawazakietal. (2009)].
この対数比例型の速度回復は,岩石試料に振動を加えた後の固有振動数の変化を測る実験においても報告されており,岩石の物性が加振前の状態に向かいゆっくりと回復していくこの現象は,スローダイナミクス(slow dynamics)として知られている[例えば,TenCateetal. (2000)].

2点間の伝播速度を(Active なり Passive なり)何らかの手法で求める、位相 or 群速度をモニタリングする、といった2段階の流れ。ambient noise だと長期間のデータが必要になるので時間の解像度が落ちます。ちなみに、Hi-net のデータを1年分DLしようとすると、1か所あたり丸1日かかります。そういう意味でも、余震のコーダ波を利用するのが現実的でしょうか。


2021年5月16日日曜日

passive 手法

passive 手法のテスト。

*** CC-FJpy ******************************
https://github.com/ColinLii/CC-FJpy
CC-FJpy: A Python Package for seismic ambient noise cross-correlation and the frequency-Bessel transform method
・複数の観測点間の相互相関関数をプロット。
・このコードに限りませんが、基本、長期データが必要です。F-netだと申し込み時に制限(1GB以下、リクエスト数(容量)制限)があるため、難しい。

*** noisepy ******************************
https://noise-python.readthedocs.io/en/latest/
・複数の観測点間の相互相関関数をプロット
・SAC を ASDF に変換して読み込み。
・noise_module が import できないエラー(S0B_to_ASDF.py)。

*** AANT *********************************
http://penglipk.github.io/AANT/
Adaptive ambient noise tomography and its application to the Garlock Fault, southern California
・sac 形式のデータを読み込み、速度を推定し、速度場を逆算してくれます(トモグラフィー)。
・cutsac のソースが含まれておらず中止。
・理想的なツールなのですが、更新が止まっています。非常に惜しい。

*** BIDO *********************************
https://staff.aist.go.jp/ikuo-chou/BIDO/2.0/bidodl.html
・時間領域と周波数領域の結果を比較する目的で実施。
・2.0 では 水平動も利用されていましたが、 3.0 では鉛直動のみ扱うとのこと。
・Ubuntu の binary が入っていましたが、動かなかったのでコンパイルしました。
・15分程度のデータであれば動作。速度が異常。
・1時間程度のデータは読み込み時にエラー。ディメンジョン修正箇所がわからず。

*** FMST *********************************
http://rses.anu.edu.au/~nick/surftomo.html
・トモグラフィーのみのようでした。分散曲線を出した後のステップです。
・速度を事前に用意することで、それらを満たす速度場を逆算。
・マニュアルを読みましたが、今回は手を出さず。

*** seismic-noise-tomography *************
https://github.com/bgoutorbe/seismic-noise-tomography
Rayleigh wave group velocities at periods of 6–23 s across Brazil from ambient noise tomography
・Python2.7でした。ライブラリがそろっていた 3.8 環境で動かし始めましたが、変換が面倒になり断念。

*** amb_noise_processing *****************
http://hestia.lgs.jussieu.fr/~boschil/downloads.html
Two-receiver measurements of phase velocity: cross-validationof ambient-noise and earthquake-based observations
・これもPython2.7でした。変換が print 文のみでしたので、3.8 にて実施。


この他、MATLAB を利用した GUI 対応コードもいくつかありました。さすがに購入しようとは思いませんでしたが、仕事で必要になれば選択肢に入るでしょう。

私の環境では、最後の amb_noise_processing のみ分散曲線まで確認できました。詳細を追いましょう。

************************************
20210516追加
「Two-receiver measurements of phase velocity: cross-validationof ambient-noise and earthquake-based observations」より、passive 手法で共通と思われる部分。

  • Microseisms are most energetic in the period band between 5 and 30 s. At the shorter period end of this range, surface waves from teleseismic earthquakes are difficult to observe. Local to regional seismicity is needed to obtain phase velocities from earthquake data in the frequency range between about 3 and 10 s (Endrun et al. 2004).
  • 1年のデータが必要(1年以上は不要)
  • Records are filtered in the frequency band of interest (2–200s). サンプリング周波数1Hz以下⇒容量削減、計算負荷軽減。
  • The data is then cut into overlapping time windows (e.g. half an hour with 50 per cent overlap), cosine-tapered at both ends.

 

2021年5月15日土曜日

CPS 3.30

 先週入れた  CPS 3.30 のコマンドです。

 ひとまず、動作確認ができたモノから記載。

*** 震源(dist)追加 *********************************************
$ gsac
GSAC> r *.U
GSAC> ch evla -4.52 evlo 101.374 evdp 34
GSAC> wh
GSAC> quit          

*** plot *********************************************************
GSAC> sort up dist
GSAC> p

*** mft *********************************************************
$ sacmft96 -f 'filename'.sac -PMIN 4.0 -PMAX 100.0 -a0 100.0 -A -VMIN 2.0 -VMAX 5.0 -U cm -R -S
$ plotnps -K < MFT96.PLT > temp.ps

*** pom ********************************************************
$ ls *.sac > cmdfil
$ sacpom96 -C cmdfil -PMIN 4.0 -PMAX 100.0 -nray 100 -A -VMIN 2.0 -VMAX 5.0 -R -S
$ rm cmdfil
$ plotnps -K -R < POM96.PLT > temp.ps

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

(理解できていませんが)active 手法はひとまずクリア。あとは passive。

2021年5月8日土曜日

表面波の伝播 その2

先日の文献で紹介されていたツールです。

Computer Programs in Seismology 3.30
http://www.eas.slu.edu/eqc/eqccps.html


MFT が提供されていましたので、試しに使ってみました。

Cygwin 環境を作るのはイマサラなので、Linux 環境を選択。SSH や WSL2 だと GUI がダメなので、手元の laptop を Win10 との dual boot に変更(途中、BIOS を触って Win の boot loader が立ち上がるようになってしまいましたが、grub を選択して F10 で書き込めばOKでした(機種によって書き込み方が違うので、また迷わないように( ..)φ))。

これで Ubuntu 20.04 が入ったので、次はCPS。インストール手順は以下の通り。

$ sudo apt install make gcc gfortran

隠しファイル表示
/home/.profile にパス追加
PATH=$PATH:/home/CPS/PROGRAMS.330/bin
export PATH

このまま make すると ヘッダーファイルがないといわれます。追加。
Xlib.hの追加
$ sudo apt install libx11-dev
curses.hの追加
$ sudo apt-get install libncurses5-dev

$ ./Setup LINUX6440
$ ./C

これでOK。で、起動。
MFT は GUI を備えていましたので、まずこれから触ってみました(do_mft)。
残念ながら、Hi-net は読み込めず。F-net は読めました。同じ防災科研さんのツールで SAC 形式に変換していたのですが、前者についてはヘッダーを読めなかったようです。gsac では Hi-net も読み込み & プロット可能でしたので、do_mft の問題でしょうか。

続いて同じく GUI の do_pom。こちらは以下が出典。
Analysis of dispersive waves by wave field transformation
https://pubs.geoscienceworld.org/geophysics/article-abstract/46/6/869/68551/analysis-of-dispersive-by-wave-field
相変わらず理屈はついていけないのですが、結果は MFT よりマシ。先日の位相速度と同程度でした。うーん。

なかなか前に進みません。やはり、素人がすぐに結果を出せるほど甘くない内容です。

**************************************
20210511追記 その1

追加すべきライブラリは TOP ページに書かれていました。
http://www.eas.slu.edu/eqc/eqccps.html

Ubuntu Linux

You can use the program download page, or open a terminal window and enter the following commands:

sudo apt-get install make
sudo apt-get install gcc
sudo apt-get install gfortran
sudo apt-get install xorg-dev
sudo apt-get install libncurses5-dev
sudo apt-get install gawk

the following programs are useful:

sudo apt-get install gv
sudo apt-get install ImageMagick


追記 その2
F-net も読めないモノがありました。長いと読めないのでしょうか?

20210514追記
do_mft のバックで動いている sacmft96 などでは読めました。do_mft、do_pom のファイル認識に不具合があるようです。

 
20210626追記
Kimiyuki Asano et al.(2017)Surface wave group velocity in the Osaka sedimentary basin, Japan, estimated using ambient noise cross-correlation functions

The group delay time of the surface wave at each period was estimated using the multiple flter analysis technique (Dziewonski et al. 1969). 

MFTから求められる群速度。まだ理論は理解できていません。
http://www.eas.slu.edu/eqc/eqc_cps/TUTORIAL/LessonA/index.html
http://www.eas.slu.edu/eqc/eqc_cps/TUTORIAL/EMPIRICAL_GREEN/MFT.pdf

2021年5月7日金曜日

表面波の伝播

Hi-net、F-net から表面波を取り出せないか試行錯誤していました。

先日見た時間領域フィルタの文献では、遠地地震で綺麗に表面波が現れていました。
これはS波より遅れて到達する点、継続時間が分~時間単位になる点で、見た目だけでも同定は容易でしょう。が、震源が観測局に近接している場合や微小地震では判定困難。地震分野では、この表面波速度をどのように同定しているのでしょうか?

調べてみると、ありました。

まず、passive な手法。
微動+相互相関を利用した表面波位相速度の推定です。この手法は図書にも載っていました。

・林田ほか「中京堆積盆地における表面波群速度の推定-Hi-net連続地震観測記録を用いた地震干渉法に基づく検証-」
https://www.jstage.jst.go.jp/article/zisin/66/4/66_127/_article/-char/ja/
・Bruno Goutorbe et al.(2015) Rayleigh wave group velocities at periods of 6–23 s across Brazil from ambient noise tomography
https://academic.oup.com/gji/article/203/2/869/578005 

passive な手法は、土木分野でもハイブリッド表面波探査として利用されています。アレー探査から着想されたのでしょう。また、地震後の剛性回復チェックにも利用されているようです。
このような微動でレーリー波位相速度を求めておき、地盤構成を推定しておけば、地震が起こった際に各周波数で目立つ波がレーリ波かどうかチェックできるかもしれません。

次に、 active な手法。地震動を利用します。土木分野でよく利用されるアクティブ探査と同様の考え方です。違いは1イベントかつランダムな受振点配置。震央距離を利用するようですが、どのように群速度を計算するのでしょうか?

・Corchete V. et al.(2007) Shear wave velocity structure of the Sinai Peninsula from Rayleigh wave analysis
https://link.springer.com/article/10.1007/s10712-007-9027-6
・Dziewonski A. et al.(1969) A technique for the analysis of transient seismic signals.
https://pubs.geoscienceworld.org/ssa/bssa/article-abstract/59/1/427/101646/A-technique-for-the-analysis-of-transient-seismic

前者は以下の2種を組み合わせているようですが、頭がついていきません。
・Multiple Filter Technique (MFT)
・Time Variable Filtering (TVF)

GitHub 上で両者を扱うソースを順に見てみましたが、どれも一長一短。記載通り Python2.7 の環境を作っても動かなかったり、fortran だったり(最初はインタラクティブ希望)、面倒だったり。

ひとまず、弾性波探査のP波のピックアップと同じように、各周波数帯の波形から位相速度を取り出してみることにしました。
が、観測点が100、対象周波数が20あれば、2000個の波形に対してマニュアルピッキングになります。これは非効率。自動でできないかと obspy の Trigger/Picker をいくつか試したものの、どれもイマイチ(そういえば、表面波探査の解析では OYO さんのソフトで自動処理していましたが、どのようにピックアップしていたのでしょう?)。
結局、record section にして直線を引くことで平均位相速度を見出し、それを中心周波数に対してプロットしました。エンベロープを書けば、群速度も得られると思いますが、ひとまずここまで。

手を動かすほど、理解していなかった点、学ぶべき点が見えてきました。
機会があればプロに指導を仰ぎたいですね。

*************************************
20210508追記

現代地球科学入門シリーズ6「地震学」
p95-97 「6.6 表面波分散曲線の基づく推定」
・最終的に実施した手法
・鉛直成分でレーリー波、トラバース成分でラブ波を測定
p273-274 「21.2 常時微動の相互相関関数解析に基づく地震波伝播速度の推定
・passive
・地震直後の速度低下とその後の回復


2021年5月1日土曜日

C言語の基礎

GW 中に、C言語を眺めておこうと考えていました。

英語同様に、辞書があれば読めるものの全く書けない、という状況です。必要に迫られることが少なく、積極的に利用する機会はありませんでした。が、書ける能力があれば楽だったかも、ということは数回ありました。将来?に備えた弱点補強です。

取り急ぎ作りたいモノもないので、いくつかのサイトを見て1週間程度で完了できそうな以下を実施。

一週間で身につくC言語の基本
http://c-lang.sevendays-study.com/index.html

基礎編はほぼ復習になりましたが、応用編は知らないことばかり。あわせてメモしておきましょう。

基礎編
・void とは、関数の戻り値が無いことを意味する。
・Visaul Studio では、セキュリティの関係で scanf を使えない(scanf_s ならOK)。
 プロジェクトを右クリック → プロパティ → c/c++ → SDL チェックを「いいえ」に変更
 ※画面に直接入力する仕様は少ない(が、たまにあると忘れる)のでメモ。
・7日目、ファイル分割の御作法。
 #ifndef、#define、#endif マクロの意味、2重インクルードの防止、global 変数に対する extern修飾子。

応用編
2日目 変数のアドレス
・通常の変数:int a    値a    アドレス&a
・ポインタ変数:int* p    値*p    アドレスp
    p = &a:    pに、aのアドレスを代入    *pはaと同じものになる
    *p = 300:    *pに300を代入
3日目 ポインタと配列
・ポインタへのアドレス代入:p1 = &ar1[0]
    配列変数    アドレス    ポインタ    ポインタ変数の値
    ar1[0]    &ar1[0]    p1        *p1
    ar1[1]    &ar1[1]    p1+1    *(p1+1)
・ポインタ変数p1に配列d[3] = { 0.2 , 0.4 , 0.6 }を代入:配列変数dの先頭アドレス&d[0]の値を代入するのと同じ処理
・動的メモリ確保 malloc()~解放 free()を使う場合、stdlib.hをインクルード
5日目 構造体
・構造体の成分の変数「メンバ」:(構造体変数名).(メンバ)
・データ渡し⇒スタック領域の圧迫、データコピー(無駄な処理)発生。⇒ポインタ渡しが通常

初心者がまず躓くといわれているポインタ。shallow・deep copy の関係と似ているようです。概念は理解できましたが、慣れるまでに時間がかかるでしょう。

丁寧な御説明で分かりやすく、2日で終わりました。感謝。
今回はココまで。