フーリエスペクトル解析

フーリエスペクトル解析 †

（フーリエ）スペクトル解析は，ランダムデータから有用なシグナルを抽出するために，地震記録の分析においても頻繁に使われる基本的な解析手法である．通常の地震観測記録は，揺れの振幅の時間変化を記録するいわゆる時系列データであるが，地盤や建築構造物はある特定の周期（固有周期）の地震動を増幅させる特徴があることから，時系列波形データを周期ごとの振幅，すなわちスペクトルに分解することにより，地盤の増幅特性や構造物の応答特性を明らかにすることができる．また地盤や建物に対する入力としての地震動は，震源断層のせん断すべり運動によって発生し，地表に伝播したものであるが，こうした入力地震動に対するスペクトル解析によって，震源の破壊過程や伝播経路における減衰特性の解明が期待される．

↑

1. フーリエスペクトル †

↑

1.1　フーリエ級数 †

時間の関数 x(t) が区間 $\left[-T/2,T/2\right]$ を基本とする任意の周期の関数であるとき，次式のように三角関数 $\cos \frac{2n\pi t}{T}$ と $\sin \frac{2n\pi t}{T}$ の線形和で級数展開することができる．

$x(t)\sim \frac{a_0}{2}+\sum_{n=1}^\infty \left( a_n\cos\frac{2n\pi t}{T}+b_n\sin\frac{2n\pi t}{T} \right)$

このとき，右辺の級数をフーリエ級数またはフーリエ展開，その係数 a_n, b_n をフーリエ係数と呼ぶ．フーリエ係数は次式により求められる．

$a_n=\frac{2}{T} \int_{-T/2}^{T/2}x(t)\cos\frac{2n\pi t}{T}dt \\ b_n=\frac{2}{T} \int_{-T/2}^{T/2}x(t)\sin\frac{2n\pi t}{T}dt$

なおここで用いられている関数 $\cos\frac{2n\pi t}{T}$ と $\sin\frac{2n\pi t}{T}$ は直交関数列であり，次のような特徴を持つ．

$\begin{align} \int_{-T/2}^{T/2}\cos\frac{2m\pi t}{T}\cos\frac{2n\pi t}{T}dt =\int_{-T/2}^{T/2}\sin\frac{2m\pi t}{T}\sin\frac{2n\pi t}{T}dt &=0 (m \neq n) \\ &=T/2 (m=n) \end{align}$

さらにm=nのときの積分値を1に正規化すると（この場合は係数 $\sqrt{2/T}$ を乗じる）正規直交関数列となる．

三角関数を指数関数に置き換えることにより，次のように複素フーリエ級数が導かれる．

$x(t)\sim \sum_{-\infty}^\infty c_ne^{i2\pi nt/T}\\ c_n=\frac{1}{T} \int_{-T/2}^{T/2}x(t)e^{-i2\pi nt/T}$

↑

1.2 フーリエ積分 †

複素フーリエ級数の周期区間を無限大とする極限を考えたとき，次の複素フーリエ積分が成り立つ．また，関数 x(t) の無限積分が収束し（絶対積分可能であり），全区間で連続，かつ区分的に滑らかであるならば，フーリエ積分はもとの関数 x(t) に一致する（フーリエの積分定理）．

$x(t)=\int_{-\infty}^{\infty}X(\omega)e^{i\omega t}d\omega \\ X(\omega)=\frac{1}{2\pi}\int_{-\infty}^{\infty}x(t)e^{-i\omega t}dt$

このとき， $\omega=2\pi f$ は角周波数，は周波数で周期の逆数である．

↑

2.3 フーリエ振幅スペクトルと位相スペクトル，および位相差分スペクトル †

複素フーリエ係数 c_n やフーリエ成分 $X(\omega)$ は複素数であり，実数 a_n, b_n ，または $a(\omega), b(\omega)$ を用いて次式のように表すことができる．

$c_n=a_n+ib_n=\sqrt{a_n^2+b_n^2} e^{i\theta_n}\\ \theta_n=\tan^{-1}(b_n/a_n)$

$X(\omega)=a(\omega)+ib(\omega)=\sqrt{a(\omega)^2+b(\omega)^2} e^{i\theta(\omega)}\\ \theta(\omega)=\tan^{-1}(b(\omega)/a(\omega))$

このときの絶対値 $|X(\omega)|=\sqrt{a(\omega)^2+b(\omega)^2}$ はフーリエ振幅スペクトルを構成し，角周波数 $\omega$ の波の振幅を表す．また偏角（位相角） $\theta(\omega)$ は位相スペクトルを成す． x(t) が実数であるとき， $|X(\omega)|=|X(-\omega)|$ ， $\theta(\omega)=-\theta(-\omega)$ が成り立ち， x(t) は次式のような実数関数の積分として表現される．

$x(t)=2\int_0^\infty|X(\omega)|\cos(\omega t+\theta(\omega))d\omega$

上式は，任意の時系列データ x(t) が，位相差 $\theta(\omega)$ で重なり合う振幅 $|X(\omega)|$ ，周期 $T=2\pi/\omega$ の余弦関数の群に分解されることを示している．

位相スペクトルに対して，隣り合う周波数の位相角の差（位相差分）の頻度分布を位相差分スペクトルと呼ぶ．位相差分は負の値，すなわち0から $-2\pi$ の間で定義される．地震波形記録のような非定常波の位相差分スペクトルは，元となった時系列データの包絡形状によく似ることが指摘されている．

↑

3 自己相関関数とパワースペクトル †

時間の関数 x(t) に対し，自己相関関数 $C(\tau)$ は次式に示すような $\tau$ 時間離れた関数の積の平均値として定義される．

$C(\tau)=E[x(t)x(t+\tau)]$

定常確率過程の場合，上式は時間平均で置き換えられ，次のように時間ずれ $\tau$ のみの関数となる．

$\begin{align} C(\tau)&=\overline{x(t)x(t+\tau)}\\ &=\lim_{T \to \infty} \frac{1}{T} \int_{-T/2}^{T/2}x(t)x(t+\tau)dt \end{align}$

自己相関関数は偶関数（ $C(\tau)=C(-\tau)$ ）であり， $\tau=0$ で最大値をとる． $C(\tau)$ を C(0) で正規化した $R(\tau)=C(\tau)/C(0)$ を自己相関係数と呼ぶ． R(0)=1 である．

一方，パワースペクトル密度関数 $S(\omega)$ は，角周波数 $\omega$ の波のエネルギーに対する時間平均として，次のように定義される．

$S(\omega)=\lim_{T \to \infty} \frac{1}{T} 2\pi|X(\omega)|^2$

なお， $S(\omega)$ は時間平均をとることからパワースペクトル密度関数と呼ばれるが，単純にパワースペクトルと記述される場合も多い．自己相関関数 $C(\tau)$ とパワースペクトル $S(\omega)$ は次式で示されるように，互いにフーリエ変換と逆フーリエ変換の関係にある(Wiener-Khintchineの公式)．

$C(\tau)=\int_{-\infty}^\infty S(\omega)e^{i\omega\tau}d\omega\\ S(\omega)=\frac{1}{2\pi} \int_{-\infty}^\infty C(\tau)e^{-i\omega\tau}d\tau$

ここで， $\tau=0$ とおくと，次式が得られる．

$\begin{align} C(0)&=\lim_{T \to \infty} \frac{1}{T} \int_{-\infty}^\infty x^2(t)dt \\ &=\int_{-\infty}^\infty S(\omega)d\omega \end{align}$

すなわち，パワースペクトル（密度） $S(\omega)$ は，時系列データの平均パワー $\overline{x^2}$ に対する各周波数成分からの寄与の割合を表している．

↑

4 相互相関関数とクロススペクトル，コヒーレンス †

二つの時間変動の関数， x(t) と y(t) に対し，自己相関関数と同様に時間遅れ $\tau$ を伴う積の平均値として，次式のように相互相関関数 $C_{xy}(\tau)$ と相互相関係数 $R_{xy}(\tau)$ を定義する．

$\begin{align} C_{xy}(\tau)&=\overline{x(t)y(t+\tau)}\\ &=\lim_{T \to \infty} \frac{1}{T} \int_{-T/2}^{T/2}x(t)y(t+\tau)dt \end{align}$

$\begin{align} R_{xy}(\tau)&=\overline{x(t)y(t+\tau)} / \sqrt{\overline{x^2}}\sqrt{\overline{y^2}}\\ &=C_{xy}(\tau) / \sqrt{C_x(0)C_y(0)} \end{align}$

相互相関関数は，システムに対する入出力の関係（例えば入力地震動と建物の応答）や，ランダム変動場の空間相関などを推定する際によく用いられる．

相互相関関数のフーリエ変換として，クロススペクトル $S_{xy}(\omega)$ が定義される．

$C_{xy}(\tau)=\int_{-\infty}^\infty S_{xy}(\omega)e^{i\omega\tau}d\omega\\ S_{xy}(\omega)=\frac{1}{2\pi} \int_{-\infty}^\infty C_{xy}(\tau)e^{-i\omega\tau}d\tau$

さらにクロススペクトル $S_{xy}(\omega)$ は，時系列データ x(t) と y(t) のフーリエスペクトル $X(\omega)$ と $Y(\omega)$ を用いて以下のように書ける．

$S_{xy}(\omega)=\lim_{T \to \infty} \frac{2\pi}{T}X^*(\omega)Y(\omega)$

ここで*は共役関係を表す．

クロススペクトルの二乗を二変数のパワースペクトルで正規化した量，またはその平方根をコヒーレンス $coh^2(\omega)$ ， $coh(\omega)$ と呼ぶ．

$coh^2(\omega)=\frac{|S_{xy}(\omega)|^2}{S_{xx}(\omega)S_{yy}(\omega)}\\ coh(\omega)=\sqrt{coh^2(\omega)}$

コヒーレンスは0から1の間で実数値をとる．またコヒーレンスは，入出力システムにおけるノイズの影響，システムのモデル化の誤差，システム特性自体の時間変動（非線形挙動）を反映する量として，主に解析精度の検証に用いられる．

↑

5 フーリエスペクトル比と伝達関数 †

地盤応答特性の推定や構造物の振動同定においては，対象となる地盤や構造物を地震動の応答システムとみなし，入力記録（基盤地震記録や建物基礎位置の記録）と出力記録（地表記録や上階の記録）からシステムの応答特性を推定する．このときの周波数領域における応答特性関数を一般に伝達関数と呼ぶ．いま，入力側の地震動 x(t) が対象応答システムにほぼ完全に入力する（ノイズや散逸を考慮しなくてもよい）と考えた場合，出力側の観測値 y(t) の推定誤差を最小とする伝達関数 $H_1(\omega)$ は次式で与えられる．

$H_1(\omega)=\frac{S_{xy}(\omega)}{S_{xx}(\omega)}$

このようにして得られた伝達関数 $H_1(\omega)$ と入力地震動 $X(\omega)$ から推定される出力と，実際に観測される出力 $Y(\omega)$ とのエネルギー比はコヒーレンス $coh^2(\omega)$ に一致する．すなわちコヒーレンスは同定された応答システムの信頼性を測る目安となる．

$coh^2(\omega)=\frac{|S_{xy}(\omega)|^2}{S_{xx}(\omega)S_{yy}(\omega)}=\frac{|H_1(\omega)|^2S_{xx}(\omega)}{S_{yy}(\omega)}$

一方，地盤応答の伝達関数の推定では，入力となる地中地震波に表面波などの散逸成分が多く含まれていると考えられることから，以下に示すフーリエスペクトル比 $H_R(\omega)$ を伝達関数として用いることが多い.

$H_R(\omega)=\sqrt{\frac{S_{yy}(\omega)}{S_{xx}(\omega)}}$

なお H_1 と H_R の間には以下の関係が成り立つ．

$H_R(\omega)=\frac{|H_1(\omega)|}{\sqrt{coh^2(\omega)}}$

コヒーレンスは常に1以下の正の実数になるので， H_1 により与えられる伝達関数推定値は，フーリエスペクトル比 H_R に等しいか，より小さな値を示す．

↑

6 FFTとスペクトルの平滑化 †

ここまで，フーリエスペクトルとそれに関わるいくつかの概念を提示してきた．実際の観測記録は時間長が有限であり，また計算機処理の際には離散値を扱うので， $\infty$ までの計算は行われない．現在，地震動のフーリエ解析では高速フーリエ変換FFT(Fast Fourier Transform)と呼ばれる計算手法がもっとも一般的に用いられている．FFTのプログラムは多くのコンピュータライブラリに含まれており，また下に載せた参考文献にもFortranプログラムが掲載されているので，詳細はそちらを参照していただくこととし，ここでは使い方の注意点のみ以下に簡潔に記す．

時系列データ個数は個（nは整数）に固定する．データ数が個に達しない場合は，になるまで後続に0振幅を追加する．
時系列データには不連続点を作らないようにする．不連続点があるとGibbs現象と呼ばれるフーリエ級数和の振動・飛び出しが生じる．特にデータの最初と最後が有限の振幅を持たないようにする．
データのサンプリング周期を $\Delta t$ ，データ数をとしたとき，周波数のサンプリング $\Delta f$ は $1/(N\Delta t)$ となる．このとき，N/2番目の周波数 $1/(2\Delta t)$ が検出しうる最大の周波数（Nyquist周波数）であり，元のデータにこれよりも高周波の成分が入っているとエイリアジングと呼ばれる誤差が生じる．このため，地震観測においては，デジタルデータにNyquist周波数を超える成分が混入しないよう，事前にローパスフィルターを掛けておく必要がある．

実地震データの解析から得られたフーリエ振幅スペクトルは振幅の凹凸が大きく，スペクトル比を取った場合などにはさらに凹凸の度合いが増幅する．こうした問題を避けるために，スペクトルの平滑化をおこなう．平滑化の操作は，周波数領域で一定のバンド幅をもつウィンドウ関数を乗じて移動平均を行うことでなされる．平滑化でよく用いられるウィンドウ関数としてはParzenウィンドウ，あるいはそのデジタル表記に相当するHanningウィンドウなどがある．Parzenウィンドウ W(f) は次式で表される．

$\begin{align} W(f)&=\frac{3}{4}u\left(\frac{\sin A}{A}\right)^4\\ A&=\frac{\pi uf}{2} \end{align}$

一方，Hanningウィンドウは中心周波数の振幅を0.5倍し，その前後， $f\pm \Delta f$ の振幅を0.25倍して重みつき平均を取る操作に相当する．この操作をn回繰り返すことにより，nが十分に大きい場合に得られる重み分布はParzenウィンドウに近づく．ParzenウィンドウとHanningウィンドウのバンド幅（単位Hz）は，それぞれパラメータと繰り返し回数nを用いて次式で表される．