4.1、 震波干涉法
透過計算均勻散射波場(diffuse field)中任兩點之間連續紀錄的交互相關函數 (cross-correlation function)，可以得到介質的經驗格林函數(e.g., Weaver and Lobkis, 2001, 2002；Shapiro and Camille, 2004；Snieder, 2004)。相當於其中一點處施予脈衝 源(impulse source)，在另一點接收到的脈衝響應(impulse response)(圖 4-1)。而建立 經驗格林函數的過程被稱為震波干涉法(seismic interferometry)。通常利用周遭噪訊 (ambient noises)或是地震尾波(earthquake coda waves)作為經散射場充分散射的波形 來源。前人研究指出，周遭噪訊能量主要來自海浪和海床間耦合(coupling)作用產 生的微地動 (e.g., Longuet‐Higgins, 1950; Stehly et al., 2006 )，周期約為 2-20 秒，並 且主要以表面波的形式傳遞(Shapiro and Campillo, 2004)；而地震尾波為地震直達 波到達之後，在異質性(heterogeneous)介質中多重散射(multiple scattering)而達到能 量均分(equipartition)的狀態(Hennino et al., 2001)。因為尾波散射的路徑範圍廣泛，
平均了來自路徑效應(path effect)和震源輻射型態(radiation pattern)的影響，保留散 射介質中震波的訊息。從 2003 年 Campillo 和 Paul 首度使用地震尾波重建經驗格 林函數後，相關研究至今仍舊方興未艾。
震波干涉的相關研究普遍使用位於地表的測站對(station pair)建立經驗格林函 數以解析地殼內的震波訊息。而近年來的研究以高頻的人為噪訊(e.g., Miyazawa et al., 2008；Lewis and Gerstoft, 2011)或是地震尾波(e.g., Nakata and Snieder, 2012)作為 散射場的資料來源，並將該技術延伸至垂直排列的井下地震儀陣列以解析近地表 的震波特性，如 P 波及 S 波速度，及 S 波震波非均向性等。
由於臺灣地區周遭噪訊主要周期約為 2-8 秒(Chen et al., 2011)，其表面波波長 普遍大於本研究的測站間距 (< 400 m)。因此，我們嘗試以相對高頻(2-8 Hz)的地震 尾波紀錄建立井下-地表測站經驗格林函數，以量測臺灣地區近地表的震波速度和 震波非均向性(圖 4-2)。
圖 4-1、建立兩測站路徑間經驗格林函數示意圖。在散射波場下，看似雜亂的 波形紀錄經過計算兩測站間的交互相關函數，相當於在其中一站施予脈衝源，
而另一站接收到的脈衝響應。(Snieder and Wapenaar, 2010)
圖 4-2、井下-地表測站經驗格林函數正負時間軸示意圖。負時間軸(τ< 0)的
4.2、 理論背景
Lobkis 和 Weaver 在 2001 年以正常模態(normal mode)證明測站間震波紀錄的 交互相關函數等同於經驗格林函數。以下依據 Snieder 和 Wapenaar(2010)簡化其推 導流程。假設一封閉系統中能量均分的均勻散射波場由 m 組正常模態展開(normal mode expansion)表示：
𝑢(𝑟, 𝑡) = ∑
𝑢𝑚𝜔(𝑟)𝑚 𝑚
(𝑎
𝑚𝑠𝑖𝑛 𝜔
𝑚𝑡 + 𝑏
𝑚𝑐𝑜𝑠 𝜔
𝑚𝑡)
……...……… ( 1 )𝑢(𝑟, 𝑡)
：位移場(displacement)紀錄𝜔
𝑚：特徵頻率(eigenfrequncy)𝑢
𝑚(𝑟)
：特徵函數(eigenfunction)T
：時間窗(time window)長度S：模態激發能量(mode’s excitation energy)
𝛿
𝑚𝑛：克羅內克𝛿
函數(Kronecker delta)。𝛿
𝑚𝑛= { 0, m ≠ n
0)速度場經驗格林函數(( 4 )式)相加成正比關係。而正時間軸的速度場經驗格林函 關法(coda wave cross-correlation)重建水平各個分量以及垂直分量的經驗格林函數；
利用地震直達波(P、S waves)將井下紀錄對地表紀錄作解迴旋(de-convolution)，同 樣也可以建立經驗格林函數。最後量測各個分量所得的經驗格林函數的波峰走時，
以估算淺地表的 Vp、Vs 和震波非均向性。以下篇幅以 EGF (empirical Green’s function)代稱經驗格林函數，並以 CCF (cross-correlation function)代稱交互相關函 數。
4.3.1、 資料處理
應等而有所不同，因此必須個別定義每個地震不同測站接收到尾波的起始及結束 時間。由於井下測站背景雜訊較低，我們使用井下測站資料來定義尾波長度。首先，
利用巴特沃斯帶通濾波器(buttterworth band-pass filter)將井下的 E-W、N-S 及 Z 分 量地震紀錄濾波 2-8 Hz，分別處理為水平絕對振幅
H(t)
和垂直絕對振幅Z(t)
： (signal window)，而相鄰時間窗彼此重疊 90%(圖 4-4)。計算每個訊號時間窗和噪訊 時間窗內振幅標準差的比值，定義為訊噪比(signal to noise ratio)：SNR =
𝜎𝑠𝑖𝑛𝑔𝑛𝑎𝑙𝜎𝑛𝑜𝑖𝑠𝑒 ……….………( 11 )
當訊噪比低於 1 時，最後一個訊號時間窗的結尾就作為尾波的結束時間。如 圖 4-5 為 NCUH 測站接收到 2013/6/2 13:43:3.2 南投地震紀錄，選取水平尾波對訊 噪比的時間變化，而圖 4-6 及圖 4-7 為選取水平及垂直分量尾波長度的實際範例。
圖 4-3、本研究所使用 2011 至 2014 年間規模大於 4.0 的地震分佈圖。
圖 4-4、尾波長度選取示意圖。黑線為地震的絕對振幅，黑色直線為訊號時間 窗，綠色直線為噪訊時間窗，淡藍色方塊為尾波的選取時間範圍。
0
Signal window Noise window Coda duration
圖 4-5、選取水平尾波訊號之訊噪比對時間的變化。測站為 NCUH，地震時間 為 2013/6/2 13:43:3.2，震央位於台灣南投，深度 14.5 公里，芮氏規模 6.5。
圖 4-6、選取水平分量尾波長度的實際範例。紅線標示發震時間，藍色方塊標示 尾波選取的範圍。由於井下測站背景雜訊較低，尾波長度以井下測站資料定義。
地震資訊同圖 4-5 圖說。
surface
surface
borehole
borehole HLE
HLN
HL1
HL2
4.3.2、 井下資料方位角修正
建立井下-地表 EGF 之前，必須先校正井下地震儀水平分量資料的方位角。我 們從先前未濾波的速度場水平分量地震資料，將地表資料與井下資料之尾波訊號 提取出來。之後針對井下未修正方位 N-S 分量和 E-W 分量的尾波訊號，以 5˚度的 間隔順時針投影至 0˚度到 355˚度共 72 個方位分量 (0˚定義為正對未修正前井下 的 N-S 分量)。
將井下 72 個方位分量的尾波訊號分別對地表(1)N-S 分量和(2)E-W 分量的尾波 訊號作交互相關運算。流程如下：由於尾波振幅隨時間以對數形式衰減，為了避免 尾波前段直達波主導整體結果，我們把井下及地表水平尾波訊號分割成數個 100 秒
圖 4-7、選取垂直分量尾波長度的實際範例。地震資訊同圖 4-5 圖說。
surface
borehole HLZ
HLZ
的時間窗，相鄰時間窗彼此重疊 90%，並在頻率域進行交互相關運算。確保疊加後 各個時間窗所得 CCF 有相同的功率頻譜密度(power spectrum density)，在頻率域計 算 CCF 時同時進行頻譜白化(spectrum whitening)處理：
𝐶
𝑟𝑟(𝑣),(𝜔) =
|𝑣𝑣𝑟(𝜔)𝑣𝑟,∗(𝜔)𝑟(𝜔)||𝑣𝑟,(𝜔)|………( 12 )
𝐶
𝑟𝑟(𝑣),(𝜔)
：頻率域之交互相關函數，又稱相干性(coherency)𝑣(𝜔)
：速度場紀錄頻譜*
：共軛複數（complex conjugate）最後疊加所有時間窗所得的 CCF 後，以疊加個數平均疊加後 CCF 的振幅。 三次仿樣內插法(cubic spline interpolation)對曲線內插，尋找最大值以確定旋轉角度 正對地表測站 N-S 與 E-W 方向。由於 0˚定義為正對未修正前井下測站 N-S 分量，
因此正對地表測站 N-S 方向的旋轉角度即為方位角修正值；而正對地表測站 E-W
前後有明顯改變。我們以人工方式依照時間分類不同的旋轉角度，並將正對地表測 站 E-W 方向的旋轉角度資料扣除 90˚度，同時納入正對地表測站 N-S 方向的旋轉 角度資料，以中位數(median)作為井下測站整體的方位角修正值。最後依據所得時 間段的方位角修正值校正所有井下資料(圖 4-12)，圖 4-13 為使用 HGSD 測站接收 到 2013/6/2 13:43:3.2 南投地震紀錄，方位角修正前後波形比較。所有測站的方位 角修正值請詳見附錄 A。
圖 4-8、HGSD 測站井下各個方位分量對地表測站(1) N-S 與(2)E-W 分量尾波訊 號的 CCF。左圖結果顯示，當井下水平尾波訊號順時針旋轉至 210˚度左右時，
與地表 N-S 分量尾波訊號的 CCF 在零點處有最大值；右圖顯示井下水平尾波訊 號順時針旋轉300˚度左右時，與地表 E-W 分量尾波訊號的 CCF 在原點處有最 大值。
圖 4-9、正對 HGSD 地表測站 N-S 方向之旋轉角度對 CCF 原點振幅的曲線。
灰線為個別地震所得結果。
圖 4-10、正對 HGSD 地表測站 E-W 方向之旋轉角度對 CCF 原點振幅的曲線。
灰線為個別地震所得結果。
圖 4-11、HGSD 測站井下方位旋轉角度。2011-2014 年選取規模大於 4.0 地震 而得正對地表測站 N-S(黃色三角形)與 E-W 方向(紅色三角形)的井下方位旋轉 角度，於 2012 年 300 日前後有明顯改變。
4.3.3、 地震尾波交互相關法
圖 4-13、HGSD 測站方位角修正前後波形比較。範例為 2013/6/2 13:43:3.2 南投地 震紀錄，濾波 0.1-1 Hz。修正後井下測站波形相似於地表測站波形。疊加 217 筆地震結果，所得的垂直分量井下-地表 EGF。
圖 4-14、使用 CHY 測站垂直分量資料以交互相關法建立的井下-地表測站經驗 格林函數波形。EGF 中波峰走時約為 0.2 秒左右。
Time (sec)
Amplit ude
4.3.4、 地震解迴旋法
除了地震尾波訊號，使用地震直達波同樣可以利用解迴旋(de-convolution)建 立 EGF。當一地震垂直入射井下-地表測站時，在井下測站會接收到到時最早的上 行波(up-going wave)以及經過地表反射的下行波(down-going wave)；而在地表測站 則同時紀錄到上行與下行波，其到時介於井下測站接收到的上行波及下行波到時 之間(圖 4-17)。我們依據 Mehta 等人在 2007 年所提出的解迴旋技術分別建立垂直 及水平分量井下-地表的 EGF。
為了確保地震波垂直入射井下-地表測站，本研究採用 2011-2013 年間測站方 圖 4-16、CHY 測站 115˚度及 25˚度水平經驗格林函數比較。紅色波形為 25˚度 的 EGF，而綠色波形則為 115˚度的 EGF。
Time (sec)
Amplit ude
人工方式挑選地表測站接收到的 P 波和 S 波到時。同時以地表 P 波到時前後 1 秒，
𝐷(𝜔)
：解迴旋函數(de-convolution function)𝑣
𝑠(𝜔)
和𝑣
𝑏(𝜔)
：地表及井下測站速度場紀錄*
：共軛複數（complex conjugate）𝜀
：阻尼係數。避免分母有極小值，本研究採用 0.05。為了增加訊噪比，將個別地震所得的解迴旋函數疊加並將函數振幅以地震個數平 均，建立井下-地表 EGF。
垂直 EGF 負時間軸的波鋒代表於上行 P 波於井下-地表測站間的走時；正時 間軸的波鋒相當於下行 P 波於測站間的走時；同理可知，水平 EGF 正、負時間軸
位角排列，EGF 中波峰走時範圍為 0.7-0.75 秒，在方位角 115˚度時波峰最早到達，
而在25˚度時最晚到達。
圖 4-17、地震垂直入射井下-地表測站示意圖。井下測站會接收到上行波和下 行波，而地表測站則同時紀錄到上行波與下行波，波形相對單純。
surface
borehole
Ś S
Ŝ S
Ŝ Ś
borehole
surface
圖 4-18、所選取 CHY 井下和地表測站垂直速度紀錄。地震時間 2012 11/18 21:02:34.04，震央位於臺灣嘉義，深度 10.26 公里，芮氏規模 2.1。
surface
borehole
HLZHLZ
HLN
HLN
HLE
surface
surface
borehole
圖 4-20、使用 CHY 測站垂直分量資料以解迴旋法建立的井下-地表測站經驗格 林函數波形。EGF 中波峰走時約為 0.2 秒左右。
Time (sec)
Amplit ude
4.3.5、 震波非均向性量測
的非均向性(anisotropic)波速(Alford, 1986; Thomsen, 1988)：𝑉
𝑠(𝜃) = 𝑉
𝑖𝑠𝑜+ 𝑉
𝑎𝑛𝑖cos 2(𝜃 − 𝜑)
………...…………..………..( 16 )度(strength of anisotropy)以及近地表剪力波分離時間(delay time, 𝑇𝑠)：
由於井下-地表測站為垂直分布，地表以上的介質沒有尾波傳遞而未達到均勻 散射波場的條件，所計算的 EGF 在時間的正負兩側並不全然對稱。負時間的 EGF 波場能量來源為地表以下的尾波散射波場；因為地表為強反射面，推論正時間的 EGF 波場能量可能由地表反射的尾波所提供，作用機制尚未明瞭。因此我們以地 震尾波交互相關法建立 EGF 並選擇負時間軸量測震波速度及非均向性。所有測站 的震波非均向性量測結果請詳見附錄 B。