羽佐田葉子 2007年3月24日 アクロス研究会＠静岡大学

Presentation on theme: "羽佐田葉子 2007年3月24日 アクロス研究会＠静岡大学"— Presentation transcript:

1 羽佐田葉子 2007年3月24日 アクロス研究会＠静岡大学
アクロス解析入門 羽佐田葉子 2007年3月24日 アクロス研究会＠静岡大学

2 アクロスデータ解析必須項目 離散フーリエ変換 …デジタル信号処理 アクロス信号とは ノイズと誤差の統計学 …確率・統計
離散フーリエ変換 …デジタル信号処理 複素数の表現 アクロス信号とは ノイズと誤差の統計学 …確率・統計 テンソル伝達関数 …線形代数学 まだあるかも･･･

3 アクロスデータ解析の流れ（１） 時間領域データ 送信信号で除算 フーリエ変換 誤差つき伝達関数 周波数領域データ 正逆回転合成 信号成分
取り出し ノイズレベル 推定 座標変換 重みつきスタッキング テンソル伝達関数 誤差つき周波数領域データ

4 アクロスデータ解析の流れ（２） テンソル伝達関数 存否イベント解析 周波数窓 波の到着イベント フーリエ逆変換 時間領域波形

5 離散フーリエ変換 Discrete Fourier transform
時間領域 time domain xn フーリエ変換 フーリエ逆変換 周波数領域 frequency domain Xk = Ck + i Sk 時間 t 周波数 f xn = Sk { Ck cos (2p fk tn) + Sk sin (2p fk tn) } 時間領域データxnをcosとsinの重ね合わせで表現したときの 各周波数成分の係数を複素数の形で表したのがXk

6 複素数の表現 z = a + i b = r exp(i q) ←オイラーの公式 実数部分 a = Re[z] real(z)
虚数部分 b = Im[z] imag(z) 絶対値 r = | z | abs(z) 偏角 q = Arg[z] angle(z) ↑Matlab表記 a b z q r Re Im

7 離散フーリエ変換 xn = Sk { Ck cos (2p fk tn) + Sk sin (2p fk tn) } 実数部分がcosの係数
データ長の中に整数周期入る周波数の波で表現 それ以外の周波数は整数周期入る波の合成で表される 時間 t …

8 離散フーリエ変換 時間領域 x ( t ) 周波数領域 X ( f ) フーリエ変換 フーリエ逆変換 時間 t 周波数 f
Dt Df = 1/(NDt) NDt 1/Dt フーリエ変換　　Xk = S xn exp (-2pi fk tn) / N　　　X=fft(x)/N フーリエ逆変換　xn = S Xk exp (2pi fk tn)　　　　　 x=ifft(X)*N 　　n = 1, …, N 　k = 1, …, N 　　tn = (n-1)Dt 　fk = (k-1)Df 　　　Df = 1/(NDt) :基本周波数 　　fNyquist = 1/(2Dt) : ナイキスト周波数 　　　ナイキスト周波数より高周波側は、低周波側の複素共役 　　　X2 = X*N X3 = X*N X4 = X*N … ↑Matlab表記

9 いろいろなフーリエ変換対 cosの場合 sinの場合 半端な周期

10 いろいろなフーリエ変換 インパルス

11 アクロスデータ解析の流れ（１） 時間領域データ 送信信号で除算 フーリエ変換 今ココ 誤差つき伝達関数 周波数領域データ 正逆回転合成
信号成分 取り出し ノイズレベル 推定 座標変換 重みつきスタッキング テンソル伝達関数 誤差つき周波数領域データ

12 アクロス観測データの離散フーリエ変換 例）100Hzサンプリング、200秒間のデータ x1, … , x20000 Dt = 0.01sec
フーリエ変換すると、 　　X1, … , X Df = 1/200 = 0.005Hz ナイキスト周波数は　1/(2×0.01) = 50Hz 50Hzよりも高周波の成分は0～50Hzの ところに現れる　→　aliasing （折り返し歪） この場合、アクロス信号は0~50Hzの範囲で 0.005Hzの倍数の周波数でなければならない！

13 アクロス信号（弾性波FM回転の場合） 例）周波数を10~20Hzの間で 変調させる。変調周期は50秒
時間 (s) 50sec 回転周波数 10Hz 20Hz 200sec 例）周波数を10~20Hzの間で 変調させる。変調周期は50秒 →アクロス信号は 搬送波周波数を基準として 1/50=0.02Hz間隔に現れる 振幅 フーリエ変換 搬送波周波数 振幅 0.02Hz ノイズ チャンネル 0.005Hz 周波数 f 周波数 (Hz)

14 アクロスデータ解析の流れ（１） 時間領域データ 送信信号で除算 フーリエ変換 誤差つき伝達関数 周波数領域データ 正逆回転合成 信号成分
取り出し ノイズレベル 推定 座標変換 今ココ 次ココ 重みつきスタッキング テンソル伝達関数 誤差つき周波数領域データ

15 ノイズと誤差の統計学 正規白色雑音 xn~N(0,s2)
白色雑音のフーリエ変換 　　実数部分と虚数部分が 　　それぞれ独立に同じ分布に従う Re[Xk]~N(0,s2/(2N)), Im[Xk]~N(0,s2/(2N)) 標準偏差は s/(2N)1/2、データ長の-1/2乗に比例　→スタッキングの効果

16 ノイズと誤差の統計学 アクロスデータに含まれる誤差の推定 実際の自然のノイズレベルは周波数依存 推定した誤差をスタッキングの重みに使用
ノイズチャンネルからノイズレベルを推定 ノイズを正規白色雑音と仮定すると、 ノイズチャンネルの振幅の２乗平均から 誤差の分散が推定できる 実際の自然のノイズレベルは周波数依存 適当な周波数範囲を区切って、その中では分散が一定と仮定して推定 推定した誤差をスタッキングの重みに使用

17 アクロスデータ解析の流れ（１） 時間領域データ 送信信号で除算 フーリエ変換 誤差つき伝達関数 周波数領域データ 正逆回転合成 信号成分
取り出し ノイズレベル 推定 座標変換 重みつきスタッキング テンソル伝達関数 誤差つき周波数領域データ 今ココ

18 テンソル伝達関数 伝達関数 transfer function とは… 線形システムの特性を表す関数
ある周波数の信号が、そのシステムによってどのように変化するか 絶対値 | H(w) | が振幅の増幅率を、 偏角Arg[H(w)]が位相の進みを表す 線形システム 入力 X(w) 出力 Y(w) 伝達関数 H(w) = Y(w)/X(w)

19 テンソル伝達関数 アクロスで求まる伝達関数 弾性波アクロスの場合 受信信号を送信信号で割り、地下の伝播特性を求める
入力：震源の力ベクトル(N) 出力：地面の変位ベクトル(m) 正逆回転の合成、座標変換により求める テンソル伝達関数 与えた力の向き、 大きさに対する 変位の向き、大きさが 分かる Ur HrR HrT HrZ FR Ut HtR HtT HtZ FT Uz HzR HzT HzZ FZ =

20 アクロスデータ解析の流れ（２） 今ココ テンソル伝達関数 存否イベント解析 周波数窓 波の到着イベント フーリエ逆変換 時間領域波形

21 伝達関数を時間領域で見る 伝達関数（周波数領域） ＝入力が全て１の場合の出力 フーリエ逆変換で時間領域に変換すると？
　　　　＝入力が全て１の場合の出力 フーリエ逆変換で時間領域に変換すると？ 振幅が１の全ての周波数のcosを足し合わせると時刻0のインパルス 全ての周波数の出力を足し合わせると、時刻0にインパルスを入力したときの出力 時刻0のインパルスがどんなふうに伝わるか？ 線形システム 入力 X(w) 出力 Y(w) 伝達関数 H(w) = Y(w)/X(w)

22 伝達関数を時間領域で見る 時間領域で見ることで、 別の経路を通ってきた波を 分離できる
現実には全ての周波数を 観測できないので インパルスにはならない なるべく広い周波数範囲を 観測すればシャープなパルスが見える 存否イベント解析を使うと 狭い周波数範囲でも それなりにパルスを分離可能 震源 時間(s) 10km 反射波 20km 30km P波 S波

23 アクロスデータ解析必須項目 離散フーリエ変換 …デジタル信号処理 アクロス信号とは ノイズと誤差の統計学 …確率・統計
離散フーリエ変換 …デジタル信号処理 複素数の表現 アクロス信号とは ノイズと誤差の統計学 …確率・統計 テンソル伝達関数 …線形代数学 まだあるかも･･･