國立臺灣大學理學院地質科學系 碩士論文
Department of Geosciences College of Science
National Taiwan University Master Thesis
台灣中央山脈地區地震所引起淺層地殼速度與靜態體積應變 擾動之相關性
Coseismic Velocity Reduction Correlated with Volumetric Strain Change Induced by Recent Large Earthquakes in the Central
Range of Taiwan 吳欣玫
Sin-Mei Wu
指導教授：洪淑蕙 博士 Advisor: Shu-Huei Hung, Ph.D.
中華民國 103 年 7 月
誌謝
首先非常感謝我的指導教授洪淑蕙老師，研究所兩年間無論是在研究上、課 業上都是非常重要的引領人以及知識的寶庫。老師在教學上非常的認真，對於學 生的疑問都會耐心的指導；在研究上老師讓我有自己思考學習的空間並適時的給 予重要的建議與協助，甚至於生活中的困難也都盡量給予幫助。感謝老師， 這兩 年間我覺得自己學到了非常多知識，並且對於未來有好的規劃與期許。也很感謝 郭本垣老師，因為大學時期參加中研院暑期生計畫，老師的細心指導與特立獨行 的教學風格，開始對於研究題目漸漸產生興趣，是我選擇繼續就讀研究所不可或 缺的契機。在研究時期梁文宗老師也幫助了我非常多，除了資料的提供，不管在 研究上或是任何的問題老師都非常熱心且有耐心的撥出時間跟我討論，真的是十 足欽佩老師對待學生的熱忱與付出！感謝龔源成老師在教學或研究上都有非常獨 到的見解，讓我體會到研究上可能有不同的發展或是應該要去彌補的地方，也很 欣賞老師獨特的個性以及教學風格。感謝李憲宗老師在研究上提供了非常重要的 幫助以及對於仔細的批改碩士論文與提出了重要的見解。感謝喬凌雲老師對於論 文認真的批閱以及口試時提出了重要的疑問讓我能再重新思考研究內容，得到更 進一步的心得。
由於本身能力不足，在研究過程中常常遇到許多問題與困難，在此要好好感 謝欣穎學姐對於我的任何疑問都非常積極且耐心的幫我尋找答案，以及柯彥廷學 長和侑頻學姐總是能熱心認真的給予我幫助和提供研究上或生活上非常有用的資 訊。謝謝珮如學姐總是辛苦的擔當助教跟我們一起討論作業或是扛起課程領導人 的角色。謝謝陳映年學長與莊育菱學姐總是很熱心提供研究上的資訊以及討論，
謝謝詹忠翰學長以及趙韋安和黃信樺學長在研究上提供幫助，也謝謝尤水輝學長 提供很多幫助以及情義相挺的咖啡，謝謝梓殷學姐除了有耐心回答我的問題也常 常提起有趣的話題讓實驗室多一點色彩，最後也謝謝張毓軒學長在電腦硬體部分 給了我非常大的幫助。
再來要感謝一起在研究生活中努力的同學們，謝謝陳力維在研究上與生活上 都幫助了我非常多，也總是很熱心的分享研究資訊以及在程式上或是研究內容進 行討論。唐楚欣、李詩婷妳們是換環境到台大後一起陪伴奮鬥的朋友，謝謝妳們 不管是研究上或生活上總是有事互相幫忙，休閒時一起玩樂，是這兩年來重要的 精神支助，也謝謝小龜我們是同實驗室一起熬夜互相聊天的好夥伴！再來也謝謝 羅翊菁、林姿綺學妹們，實驗室也因妳們而豐富。
最後也是最重要的，要謝謝我的家人。謝謝你們盡全力的協助我使我無後顧 之憂的進行研究，謝謝你們支持我朝我自己想走的路前進，是我最不可或缺的精
中文摘要
地震引起破裂帶鄰近區域岩石強度下降或是孔隙率增加以及應力和應變狀態 的改變，造成區域性彈性介質受到擾動並進一步影響淺層地殼同震速度的變化。
為了探討台灣中央山脈地區 2010 年 M 6.4 甲仙地震以及 2013 年 M 6.2 南投序 列地震斷層破裂引起的地殼應變與岩石強度改變的潛在性原因，我們使用中央氣 象局地震觀測網 (CWBSN) 提供的短周期地震儀以及台灣寬頻地震網 (BATS) 提 供的寬頻地震儀連續資料，並利用周遭噪訊交互相關函數 (CCF) 來重建經驗格林 函數 (EGF)，透過比較疊加長時間與短時間的經驗格林函數表面波到時後尾波的 走時偏移來得到地震前後地殼速度隨時間變化的關係。
由頻率範圍 0.1-0.8Hz 經驗格林函數計算過後的結果顯示，甲仙與南投地震前 後均有明顯的速度改變，與測站區域位置進行比較，甲仙地震有較明顯同震速度 下降區域為震央西方數公里處餘震活動頻繁地區，同時也為逆衝斷層上盤區域；
南投地震明顯同震下降區域則為穿越震源之斷層上盤處，亦是餘震活動頻繁區 域。
透過地震鄰近區域一維速度構造加以計算經驗格林函數主要周期 3 - 5 秒表面 波尾波的敏感深度結果顯示，甲仙與南投地震速度變化結果敏感深度皆在 10 公里 之內，暗示所量測到的同震速度變化皆由淺部構造或介質受到地震造成的擾動影 響所致。為了探討彈性介質與速度改變的關聯性，我們利用有限斷層滑移分佈資 料計算同震靜態體積在空間上的分佈並與速度變化分佈進行比較，發現在淺部地 殼同震速度下降區域與靜態體積伸張區域位置結果一致。我們推論量測到的同震 速度變化與地震應力造成之淺部地殼靜態體積應變改變結果相關。
關鍵字：經驗格林函數、交互相關函數、同震速度下降、靜態體積應變
Abstract
Coseismic velocity reduction has been considered to be mediated by perturbations of stress and strain conditions in the crust and/or increased permeability/porosity of fractured rocks within the quake-damaged zones. To investigate potential changes and causes in crustal strains and rock properties of areas impacted by the earthquake ruptures of the two recent large events, Jiasian and Nantou, (Mw > 6 and focal depth > 20 km) occurring in the south Central Range of Taiwan, we construct empirical Green’s functions (EGFs) from cross-correlation functions (CCFs) of continuous ambient noise between available station pairs near the epicenters from the short-period Central Weather Bureau Seismic Network (CWBSN) and the Broadband Array in Taiwan for Seismology (BATS). The temporal variations in seismic velocity perturbations are estimated by measuring the relative time shifts of late-arriving coda waves between short-term and long-term stacked EGFs.
The resulting EGFs at 0.1-0.9 Hz show the statistically significant coseismic velocity reduction immediately after both the events. The velocity drop is detected most pronouncedly from the pairs with the interstation paths traversing through the hanging-wall block of the ruptured fault. The sensitivity of surface wave coda arrivals to shear wave speed in the dominant period range of 3-5 s is confined within the depth of 10 km, where the crust mostly experienced the coseismic dilatational strain change induced by the slip distribution from the finite-fault models. Compared with the coseismic slip distribution from GPS data and finite-fault inversion, peak ground velocity, and slip-induced volumetric strain, we suggest the coseismic velocity reduction associated with these two events is plausibly caused by the induced dilatational strain in the shallow part of the crust
目錄
口試委員會審定書 ...i
誌謝 ... ii
中文摘要 ... iii
Abstract ...iv
第一章 緒論... 1
1.1 地殼介質特性隨時間變化 ... 1
1.1.1 地震前後地殼介質變化 ... 1
1.1.2 地熱活動地殼介質變化 ... 2
1.1.3 其他地殼介質變化 ... 3
1.2 周遭噪訊與尾波特性 ... 5
1.3 研究的地震 ... 7
1.3.1 2010 甲仙地震 ... 7
1.3.2 2013 南投序列地震 ... 8
1.4 研究問題與探討 ... 10
第二章 研究方法與資料分析 ... 11
2.1 理論背景 ... 11
2.1.1 周遭噪訊交互相關函數 ... 11
2.1.2 尾波干涉技術. ... 13
2.2 資料收集及處理 ... 18
2.2.1 儀器與連續資料 ... 18
2.2.2 經驗格林函數重建 ... 18
2.2.3 參考與現時格林函數重建 ... 19
2.2.4 尾波時間窗的選擇 ... 24
2.3 量測相對走時偏移方法 ... 27
2.3.1 移動視窗交叉頻譜法 ... 27
2.3.2 拉張法 ... 29
第三章 研究結果... 33
3.1 甲仙地震前後速度變化 ... 34
3.2 南投地震前後速度變化 ... 38
3.2.1 3 月 27 日南投地震結果 ... 38
3.2.2 6 月 2 日南投地震結果 ... 41
3.2.3 南投序列地震速度擾動結果數據比較 ... 44
第四章 討論... 48
4.1 頻率與影響深度 ... 48
4.2 同震靜態體積應變 ... 50
4.2.1 靜態體積應變與速度變化關聯 ... 51
4.2.2 同震靜態體積應變與速度變化分佈 ... 52
4.2.3 地震造成地殼非線性反映 ... 53
4.3 GPS 同震地表位移與強地動觀測 ... 56
4.4 南投 3 月與 6 月地震比較 ... 60
第五章 結論... 63
圖目錄
圖 1-1、 日本 Tohoku 地震前後速度及非均向性變化 ... 3
圖 1-2、 台灣台東地震前後淺地殼速度變化 ... 4
圖 1-3、 法國留尼旺島富爾奈斯火山活動期間速度變化與波形退相干性 ... 4
圖 1-4、 2009-2010 墨西哥地區慢地震事件期間非火山長微震活動度與不同周期的 速度變化比較圖 ... 5
圖 1-5、 噪訊源分布影響經驗格林函數正負時間軸對稱性差異示意圖 ... 7
圖 1-6、 甲仙地震有限斷層逆推的滑移量分布 ... 8
圖 1-7、 南投地震 GPS 同震水平位移與餘震分布 ... 9
圖 1-8、 短週期的次要微震所引起的台灣地區周遭噪訊分佈強度. ... 10
圖 2-1、 二維周遭噪訊交互相關函數來源分佈與能量關係示意圖 ... 13
圖 2-2、 聲波多重散射模擬示意圖 ... 17
圖 2-3、 聲學實驗測試不同溫度下花崗岩樣本的速度變化 ... 17
圖 2-4、 研究使用測站與地震後一個月餘震分布圖 ... 21
圖 2-5、 S13 感應器與 STS-1 感應器速度型儀器響應圖 ... 22
圖 2-6、 單日連續資料切割視窗示意圖 ... 23
圖 2-7、 經驗格林函數疊加天數與長時間參考格林函數相關係數關係圖 ... 23
圖 2-8、 CGF 和 RGF 波形圖 ... 24
圖 2-9、 選取尾波時間窗方法示意圖 ... 25
圖 2-10、 選取尾波時間窗示意圖 ... 25
圖 2-11、 移動視窗交頻法量測走時偏差操過程 ... 31
圖 2-12、 拉張法量測走時偏差操作過程 ... 32
圖 3-1、 經驗格林函數疊加方法 ... 33
圖 3-2、 SSD-WTP 站兩種疊加方法量測結果比較 ... 34
圖 3-3、 甲仙地震 SSD-WTP 站與 CHN1-SSD 站結果 ... 35
圖 3-4、 甲仙地震 MWCSM 結果 ... 36
圖 3-5、 甲仙地震拉張法量測結果與區域地表構造的分佈 ... 36
圖 3-6、 甲仙地震 SSD-WTP 測站對震前與震 MWCSM 量測結果 ... 37
圖 3-7、 南投地震期間 3 月速度變化結果 ... 38
圖 3-8、 南投 3 月地震 MWCSM 結果 ... 39
圖 3-9、 南投 3 月地震拉張法量測結果在區域地表構造上的分佈 ... 39
圖 3-10、 南投 3 月地震 SSLB-WUSB 測站對震前與震後 MWCSM 量測結果 ... 40
圖 3-11、 南投地震期間 6 月速度變化結果 ... 41
圖 3-12、 南投 6 月地震 MWCSM 結果 ... 42
圖 3-13、 南投 6 月地震拉張法量測結果在區域地表構造上的分佈 ... 42
圖 3-14、 南投 6 月地震 TDCB-WUSB 測站對震前與震後 MWCSM 量測結果 ... 43
圖 4-1、 甲仙地震區域基態表面波相速度受剪力波速度擾動一維敏感度算核隨深 度變化與一維速度模型圖 ... 49
圖 4-2、 南投地震區域基態表面波相速度受剪力波速度擾動一維敏感度算核隨深 度變化與一維速度模型圖 ... 50 圖 4-3、 甲仙地震側向深度 5 公里同震靜態體積應變與地震前後一個月平均速度
圖 4-5、 南投 3 月與 6 月地震側向深度 5 公里的同震靜態體積應變與地震前後一
個月平均速度變化分布圖 ... 55
圖 4-6、 甲仙地震 PGA 和 PGV 在地表的分布圖 ... 57
圖 4-7、 甲仙地震 GPS 逆推得到同震滑移量投影至地表的分佈 ... 58
圖 4-8、 甲仙地震 PGA 量值與在地表的分佈 ... 58
圖 4-9、 3/27 南投地震 PGA 量值與在地表的分佈 ... 59
圖 4-10、 6/2 南投地震 PGA 量值在地表的分佈 ... 59
圖 4-11、 同震靜態體積應變於 SSLB-CHGB 測站對沿深度方向的剖面 ... 61
圖 4-12、 南投 3 月與 6 月地震前一個月平均速度擾動直方圖 ... 62
表目錄
表 1-1、 甲仙與南投地震相關資訊. ... 9
表 2-1、 甲仙與南投地震資料與測站相關資訊. ... 19
表 2-2、 甲仙地震使用測站位置資訊. ... 20
表 2-3、 南投地震使用測站位置資訊. ... 20
表 2-4、 建置格林函數參數 ... 24
表 2-5、 南投地震使用尾波時間窗與測站間距 ... 26
表 3-1、 南投 3 月地震前後一個月平均速度差值 ... 45
表 3-2、 南投 6 月地震前後一個月平均速度差值 ... 46
表 3-3、 南投 3 月與 6 月地震前後一個月拉張法平均速度差值比較 ... 47
第一章 緒論
1.1 地殼介質特性隨時間變化
過去的地球物理研究受限於天然地震測站資料的分布，主要著重於研究空間 上的靜態地球內部結構。隨著技術以及理論研究發展的進步，開始有利用重複地 震來觀測地殼速度隨時間變化 [Poupinet et al., 1984; Schaff and Beroza, 2004] 以及 計算內核自轉速率變化 [Song et al., 1996; Zhang et al., 2005]，或以人工炸測實驗來 研究動態地球內部速度構造 [Nishimura et al., 2005]。
隨著近年周遭噪訊干涉法理論的發展，可重建反映兩測站間底下速度構造的 經驗格林函數，因其不受限於地震時空分布與震源機制，只需引入兩測站連續資 料，使得地球物理研究不再只受限於空間，更開啟了時間上的研究與發展。透過 比較重建後短期與長期經驗格林函數，可以得到短期經驗格林函數在時間軸上的 訊息，常用來監測火山活動期間的速度變化、地震前後速度與非均向性等介質特 性的改變，或是地熱或油氣層等構造的變化 [Brenguier et al., 2008a; 2008b; Nakata and Snieder, 2012; Takagi and Okada, 2012; Wegler et al., 2009]
1.1.1 地震前後地殼介質變化
地殼速度隨時間變化反映著內部介質受到應力溫度或是其他事件如地震以及 火山活動影響造成地殼彈性介質的改變。Nakata and Snieder [2012] 利用日本 KiK 地震網的地表與井下站地震波形進行解迴旋 (deconvolution) 得到在兩站間傳遞的 的 S 波波場，並藉由測量 2011 𝑀𝑤 9.0 地震發生前後 S 波的相對到時差分析震後 速度以及非均向性 (anisotropy) 的變化。結果顯示震後平均速度下降 6 %；非均向
性係數則上升60 %，並在三個月之後有明顯的恢復情形 (圖 1-1)。Yu and Hung [2012] 利用量測周遭噪訊經驗格林函數的表面波尾波走時偏移來量測 2006 年 𝑀𝑤 6.1 台東地震速度變化，發現地震引起的速度變化與在地震在破裂過程中，互 相垂直的兩個斷層交錯處產生極大滑移錯動而產生近地表淺層地殼的破壞有關 (圖 1-2)。
1.1.2 地熱活動地殼介質變化
地熱活動區域的淺部地殼，受到來自更深地底下的高溫液體或岩漿沿裂隙侵入，
造成區域性的介質特性改變。一般地熱或火山活動監測有採集溶水後特有分子的 火山氣體濃度分析、地表溫度梯度和熱流值的感測、衛星定位系統連續紀錄地表 變形量以及即時偵測火山微震活動和特有高頻水滴狀的地震波訊號。Maeda [2010]
利用測站連續紀錄的自我相關函數來觀測九州區域淺部地殼速度變化，發現慢速 異常區域與地熱和群震活動區域相關，推論觀測結果是由於張裂縫隙和流體入侵 造成速度下降所致 。 Obermann et al. [2013] 利用尾波干涉技術 (Coda Wave Interferometry, CWI) 來觀測法國留尼旺島上富爾奈斯火山活動期間速度變化與尾 波波形退相干性 (Decoherency) 的關聯，並對於不同時期火山噴發進行監測與定 位 (圖 1-3)。
1.1.3 其他地殼介質變化
Rivet et al. [2014] 利用周遭噪訊交互相關函數中表面波的尾波部分量測其到 時的偏移，發現在墨西哥格雷羅州 2009-2010 年期間的慢地震事件中地殼的速度變 化與非火山長微震活動具時空上的相關性，透過不同頻段的測量結果解析速度變 化在深度方向的分佈情形 (圖 1-4)。其他也有研究淺層速度隨季節性的變化，顯示 可能與降雨量以及所引起地下水位變化有關 [Sens-Schönfelder and Wegler, 2006;
Meier et al., 2010; Nakata and Snieder, 2012; Yu and Hung, 2012]。
圖 1-1、日本 Tohoku 地震前後速度及非均向性變化。上圖黑色資料點為地震之均 向速度，藍色線條為相鄰九個資料點平均速度；下圖黑色交叉資料點為非均向性 係數 (由快慢軸速度求得)，藍色線條為九個資料點平均。綠色直線為 Tohoku 地震 發生日期，紅色直線與灰色區域分別為平均速度與一個標準差範圍 [取自 Nakata and Snieder, 2012]。
圖 1-2、 台灣台東地震前後淺地殼速度變化，由疊加五日與三年之高頻單站自我 相關函數測量所得。上部區塊為速度隨時間的變化；下部為相關係數。兩垂直黑 色直線為同震速度變化影響期間，放大標示於右上角 [取自 Yu and Hung, 2012]。
圖 1-3、 法國留尼旺島富爾奈斯火山活動期間的地殼速度變化 (a) 與尾波波形退 相干性 (b)。不同顏色代表不同測站交互相關函數量測的結果，標示在 (c)，灰色 條帶區塊代表十月以及十二月火山噴發期間。(c) 研究區域地形圖以及測站相對位
(c)
圖 1-4、 2009-2010 墨西哥地區慢地震事件期間非火山長微震活動度與不同周期的 速度變化比較圖。(a)(b)(c) 藍色線條分別代表周期 5-6.8 秒、13-17.9 秒、19-26 秒 並分別以 5.9 秒、15.4 秒以及 22.5 秒為核心周期的速度隨時間變化；紅色線條為 頻率 2-7 赫茲非火山長微震能量隨時間變化情形。(d)為不同周期下兩者之間的交 互相關係數 [取自 Rivet et al., 2014]。
1.2 周遭噪訊與尾波特性
周遭噪訊干涉法理論基礎基於在有限空間的線性系統下，波場的傳遞經隨機不 均勻的散射以及邊界的反射後達到能量均分的擴散波場，此時將兩點的波場紀錄 做交互相關函數運算 (Cross-Correlation Function, CCF)，可重建兩點之間的經驗格 林函數 [Lobkis and Weaver, 2001; Weaver, 2005; Campillo, 2003]。
因為噪訊不受限於地震的分布時間或是震源機制的影響，近年來被廣泛地運用。
最 常 見的 是透過 頻 散分 析來進 行 周遭噪 訊表 面波層析成 像法 (ambient noise surface wave tomography) [Shapiro and Campillo, 2004; Shapiro et al., 2005; Yao et al., 2006; Bensen et al., 2007; Lin et al., 2008]。除了擷取經驗格林函數的相位資訊，噪 訊研究開始引入振幅的資訊。Lin et al. [2011] 量測 10 到 18 秒雷利波振幅能量在
空間上的衰減值，發現與內華達州和懷俄明州區域地震產生的波形所量測的衰減 量值一致。Lawrence and Prieto [2011] 利用 8 到 32 秒噪訊經驗格林函數在頻率域 間的一致性 (Coherency)，經由貝塞爾函數 (Bessel function) 來估算美洲西部沿海 空間上的衰減係數 (Attenuation Coefficient)，並進行層析成像的逆推。其結果與傳 統利用地震訊號的雙站法 (Two station method) 所求得衰減構造層析影像一致，暗 示周遭噪訊在振幅上也提供了可靠的資訊。Stehly et al. [2006] 利用加州地區測站 的垂直分量噪訊分析經驗格林函數雷利波時，發現 5 至 10 秒周期下有明顯穩定的 峰值，且在正負時間軸能量不均，推測是由主要來自於西邊海岸線的海浪拍打與 在 該 週 期 頻 段 微 震 活 躍 度 達 到 峰 值 有 關 ， 即 所 謂 的 次 要 微 地 動 (Secondary microseism) 活動有關 (圖 1-5)。
理論上，只要訊號符合能量均分的擴散波場條件就可以完整的重建在交相關函 數正負時間軸對原點呈現鏡面對稱的經驗格林函數。然而實際上地表並非理想的 均勻散射場，主要周遭噪訊源往往是來自於海浪或洋流與海岸海底地形間的互動 以及大氣的擾動 [Kedar and Webb., 2005; Shapiro, 2005; Stehly et al., 2006]，而這些 現象通常會隨季節產生變化，因此常造成測站網周圍噪訊源在方位角上的分佈不 均，其所重建的經驗格林函數的表面波於正負時間軸上的振幅與相位不同，對估 算隨時間變化的微小速度擾動易造成極大的誤差，因此本研究還應用了尾波干涉 技術以解決噪訊源分布不均的問題。
尾波不同於直達波可清楚定義其粒子振動方向及傳遞路徑，通常持續長時間的 振盪，無法用傳統波線理論來定義其路徑，可視為在不均質的速度構造中多次散
2001; Nakahara and Margerin, 2011]，雖然失去了震源的訊息但卻保留了傳播介質的 訊息，尤其對於介質中微弱變化更為敏感。常被用於實驗監測當中岩石物理性質 的變化或是地殼構造微弱的擾動 [Grêt et al., 2004]，以及火山活動的監測等 [Brenguier et al., 2011; Chaput et al., 2012; Obermann et al., 2013; Ballmer et al., 2013]。
圖 1-5、 噪訊源分布影響經驗格林函數正負時間軸對稱性差異示意圖。(a)為噪訊 源均勻分布，其經驗格林函數於正負時間軸呈現對稱。(b)與(c)為不對稱噪訊源，
交互相關後正負時間軸不對稱 [取自 Stehly et al., 2006]。
1.3 研究的地震 1.3.1 2010 甲仙地震
2010 年 3 月 4 日上午 8 時 18 分台北時間，於高雄縣桃源鄉和茂林鄉的交界處 發生芮氏規模 6.4 地震。中央氣象局 (Central Weather Bureau) 發布之地震深度約 為 22.6 公里，震央位置為北緯 22.97 度，東經 120.71 度。甲仙地震引起劇烈的搖 晃，造成台灣中南部 96 人受傷，震度在嘉南地區達到 6 級；南投與東部地區震度 4 級；北部區域則為 2 級。甲仙地震震源較深，並未出露至地表，屬於目前未知的 盲斷層。Lee et al. [2013] 利用有限斷層模型 (Finite fault slip model) 根據斷層面解 逆推計算出甲仙地震斷層破裂過程以及滑移量值分佈情形 (圖 1-6)。
1.3.2 2013 南投序列地震
2013 年 3 月 27 日上午 10 時 3 分台北時間，於南投縣仁愛鄉發生芮氏規模 6.2 地震。中央氣象局發布地震於北緯 23.9 度，東經 121.05 度，深度約為 19.4 公里。
同年 6 月 2 日下午 1 時 43 分台北時間於南投縣魚池鄉發生了芮氏規模 6.5 的地震，
震央於北緯 23.86 度東經 120.97 度，深度約為 14.5 公里。兩個地震震央距離相隔 不到 10 公里，皆造成台灣中部區域嚴重的損傷，根據中央氣象局地震測報中心解 釋，地震是由距離車籠埔斷層以東約 40 公里處的盲斷層所致，斷層型態皆為接近 南北走向之東傾逆衝斷層。南投地震主震區域位置與 1999 年集集地震主要餘震區 域大致相符，且 6 月地震深度較 3 月淺，有同一斷層往淺部破裂趨勢，因此引起 了廣泛的討論。南投序列地震 GPS 水平位移與深度的分佈見圖 1-7。
圖 1-6、 圖左為甲仙地震有限斷層逆推的滑移量分布，粉紅色星形為震央，黑 色箭頭代表每一子斷層滑移量，藍色圓圈代表規模 1.2 至 5.7 之間的餘震分布。
圖右為甲仙地震三維滑移量分布示意圖，星形代表中央氣象局公布的震源位置 [取自 Lee et al., 2013]。
圖 1-7、 上圖為 3 月 27 日 (藍色箭頭) 與 6 月 2 日南投地震 (紅色箭頭) GPS 同震水平位移 (以澎湖白沙站為基準)。紅色海灘球為南投地震主震震源機制；
黑色為集集地震餘震與其他規模大於 5.5 之震源機制。藍色與橘色圓點分別代 表集集地震前後的地震活動。淺藍與淺紫色代表 3 月與 6 月南投地震的餘震。
下圖為深度剖面圖
，
圖示同上 [取自 Chuang et al., 2013]。表 1-1、 甲仙與南投地震相關資訊。
Earthquake Date Depth Strike Dip Rake
Jiasian 2010/3/4 22.6 324 39 67
Nantou 2013/3/27 19.4 344 23 62
2013/6/2 14.5 30 37 106
1.4 研究問題與探討
由於台灣自然噪訊源分布不均，主要來自於西側的台灣海峽 [Chen et al., 2011]
(圖 1-8)，使得經驗格林函數表面波在正負時間軸不對稱，因此本研究透過尾波干 涉技術以及周遭噪訊資料的引入，利用多重散射來達到能量均分的擴散場條件，
來進行動態地殼速度隨著時間微小擾動的研究，並進一步討論介質擾動的機制與 速度的關係。
圖 1-8、 由短週期的次要微地動 (SPSM) 所引起的台灣地區周遭噪訊分佈強度，
左圖為 Z-Z 方向；右圖為 T-T 方向。箭頭長度與兩測站間交互相關函數振幅相對 所有測站對的平均振幅值的強度成正比，平均振幅值的比例尺於右下角顯示 [取自 Chen et al., 2011]。
第二章 研究方法與資料分析
2.1 理論背景
2.1.1 周遭噪訊交互相關函數
周遭噪訊理論基礎的始源是 Callen and Welton 於 1951 年證實的漲落耗散定理 (Fluctuation-Dissipation Theorem)，應用在噪訊研究立基於其線性響應理論。簡而 言之，在近乎平衡且線性的系統下受到微小短暫的擾動，在此非平衡過程中反應 會越趨於平衡，亦即可以平衡狀態之下所受的擾動性質表示。
以下進行公式推導說明有限空間的線性系統下，波場的傳遞經隨機不均勻的散 射以及邊界的反射後達到能量均分 (equipartition) 的擴散波場 (diffuse wave field)，
此時將兩點的波場紀錄做交互相關函數運算，可重建兩點之間的經驗格林函數 [Lobkis and Weaver, 2001 & Snieder and Wapenaar, 2010]。
在有限空間中的線性系統下，不考慮非彈性衰減，其彈性波場位移 u 可用常 態振盪模式展開 (normal mode expansion)：
𝑢(𝒓, 𝑡) = ∑ 𝑢𝑚𝜔(𝒓)
𝑚 𝑚 (𝑎𝑚𝑠𝑖𝑛𝜔𝑚𝑡 + 𝑏𝑚𝑐𝑜𝑠𝜔𝑚𝑡) , (2.1.1 式) 其中 𝑟 為空間中任一點位置向量； 𝑢(𝒓, 𝑡) 為在 𝑟 位置時間 𝑡 的位移向量
(displacement)；𝑎𝑚 , 𝑏𝑚為模式係數 (modal excitation function)；𝑢𝑚(𝒓) 為第 𝑚 個模 式(mode)所對應的特徵函數(eigenfunctions)；𝜔𝑚為其特徵角頻率(eigenfrequencies)
將位移場對時間做一次微分可得速度場：
𝑣(𝒓, 𝑡) = ∑ 𝑢𝑚 𝑚(𝒓)(𝑎𝑚𝑐𝑜𝑠𝜔𝑚𝑡 − 𝑏𝑚𝑠𝑖𝑛𝜔𝑚𝑡) . (2.1.2 式) 兩點之間的格林函數相當於在𝒓𝑨施予單一脈衝 (impulse) 於 𝒓𝑩得到的理論波場，
可表示成：
𝐺(𝒓𝑨, 𝒓𝑩, 𝑡) = ∑ 𝑢𝑚(𝒓𝑨𝜔)𝑢𝑚(𝒓𝑩)
𝑚 𝑚 𝑠𝑖𝑛(𝜔𝑚𝑡) 𝐻(𝑡), (2.1.3 式) 其中 𝐻(𝑡)為階梯函數(Heaviside function)，當時間為負時其值為0，時間為正時為1 對時間做一次微分可得速度場格林函數：
𝐺(𝒗)(𝒓𝑨, 𝒓𝑩, 𝑡) = ∑ 𝑢𝑚 𝑚(𝒓𝑨)𝑢𝑚(𝒓𝑩)𝑐𝑜𝑠(𝜔𝑚𝑡) 𝐻(𝑡). (2.1.4 式) 另外，由 A、B 兩測站兩速度波場進行交互相關運算得到的經驗格林函數表示成：
𝐶𝐴𝐵(𝒗)(𝜏) = 〈𝑣(𝒓𝑨, 𝑡)𝑣(𝒓𝑩, 𝑡 + 𝜏)〉 , (2.1.5 式) 其中〈 〉 代表整體平均數 (ensemble average)。將 (2.1.2 式) 代入 (2.1.5 式)，經過展
開運算整理成
𝐶𝐴𝐵(𝒗)(𝜏) = ∑𝑛,𝑚𝑢𝑛(𝒓𝑨)𝑢𝑚(𝒓𝑩)×
{〈𝑎𝑛𝑎𝑚〉𝑐𝑜𝑠𝜔𝑛𝑡𝑐𝑜𝑠𝜔𝑚(𝑡 + 𝜏) − 〈𝑎𝑛𝑏𝑚〉𝑐𝑜𝑠𝜔𝑛𝑡𝑠𝑖𝑛𝜔𝑚(𝑡 + 𝜏) −
〈𝑏𝑛𝑎𝑚〉𝑠𝑖𝑛𝜔𝑛𝑡𝑐𝑜𝑠𝜔𝑚(𝑡 + 𝜏) + 〈𝑏𝑛𝑏𝑚〉𝑠𝑖𝑛𝜔𝑛𝑡𝑠𝑖𝑛𝜔𝑚(𝑡 + 𝜏)} =
𝐹2(𝜔) ∑ 𝑢𝑚 𝑚(𝒓𝑨)𝑢𝑚(𝒓𝑩)𝑐𝑜𝑠(𝜔𝑚𝜏), (2.1.6 式) 其中 𝐹(𝜔) 為波場功率密度 (power density)。
在一個能量均分的擴散波場下，各 mode 所對應的模式係數呈現是以零為平均 值，彼此無相關性的隨機變數分佈，亦即：
〈𝑎𝑛𝑎𝑚〉 = 𝐹2𝛿𝑛𝑚，〈𝑏𝑛𝑏𝑚〉 = 𝐹2𝛿𝑛𝑚，〈𝑎𝑛𝑏𝑚〉 = 0 . 對照以上算式，(2.1.6 式) 可簡化成：
𝐶𝐴𝐵(𝒗)(𝜏) = 𝐹2{𝐺(𝒗)(𝒓𝑨, 𝒓𝑩, 𝜏) + 𝐺(𝒗)(𝒓𝑨, 𝒓𝑩, −𝜏)} , (2.1.7 式) 其中 𝜏 代表在交互相關函數計算時，B 點波場相對 A 點時間往前或往後推移的量。
𝜏 > 0 可視為在 A 點施予脈衝於 B 接收反應；𝜏 < 0 則視為 B 點施力 A 點接收。
(Interference) 效應，造成主要建設性干涉方位角範圍集中在以兩記錄點連線之間 0° 及 180° 附近，也就是說重建後的經驗格林函數主要能量以兩站連線延伸角度噪 訊源產生的波場紀錄影響最劇 [Snieder and Wapenaar, 2010] (圖 2-1)。
圖 2-1、 在二維的開放系統下周遭噪訊交互相關函數來源分佈與能量關係示意圖，
(a)紅點為噪訊源分佈，以半徑 𝑟𝑠 以及方位角 ∅𝑠 表示， 𝑋𝐴 和 𝑋𝐵 為接收測站。
(b) 𝑋𝐴和 𝑋𝐵 之間進行交互相關運算後時間遲滯與方位角的關係，綠色虛線區域標 示了主要建設性干涉範圍。(c)交互相關後 ±0.6 𝑠 內總和的經驗格林函數，可分為 正時間軸經驗格林函數 𝐺(𝑟𝐴, 𝑟𝐵, 𝑡) 以及負時間軸經驗格林函數 𝐺(𝑟𝐴, 𝑟𝐵, −𝑡) [取自 Snieder and Wapenaar, 2010]。
2.1.2 尾波干涉技術
地震波形的尾波是由一連串在介質中重複散射所架構而成的能量。由於對於 同個區域內介質不斷重複的取樣，使得尾波對於介質中的微小擾動非常敏感 (圖 2-2)。Snieder et al. [2002] 進行花崗岩實驗以溫度為操作變因來模擬聲波的多重散 射並觀察波形與速度變化。發現尾波波形在溫度較高 (50℃) 比溫度較低 (45℃) 到 時來得晚，並且有隨著散射時間越久累積走時差越大的情形 (圖 2-3)。Aki and Chouet [1975] 利用區域地震來量測時空上尾波的衰減係數 (Q) 變化，進一步探討
地殼中異質性的分佈情形。Jin and Aki [1986] 量測區域小地震尾波衰減係數隨時 間的變化，發現 1976 年唐山大地震以及 1975 年海城地震之前衰減係數皆達到異 常低值。之後許多研究也利用尾波衰減係數來監測近地表應力場隨時間的改變，
但都僅擷取尾波的振幅並沒有使用相位的資訊。 Lauterborn et al. [1995] 利用人工 干涉儀的實驗，得到較精準的多重散射後尾波的振幅以及相位的資訊。尾波的干 涉技術將介質視為自然干涉儀，隨著波場在受擾動範圍內散射距離越長，其累積 的走時差異增大的特性，廣泛應用於監測地殼介質的微小變化或是實驗中岩石物 理性質的變化等 [Poupient et al., 1984; Grêt et al., 2005; Snieder et al., 2002;
Brenguier et al., 2011; Rivet et al., 2014]。以下根據 Snieder et al. [2002] 針對尾波干 涉技術應用於監測時間序列上速度的變化進行理論公式推導:
假設未受擾動的波場紀錄 𝑢𝑢𝑛𝑝以所有可能路徑 𝑃 由費曼路徑總和 (Feynman
path summation) 函數來表示：
𝑢𝑢𝑛𝑝(𝑡) = ∑ 𝐴𝑃 𝑃𝑆(𝑡 − 𝑡𝑃) , (2.1.8 式) 其中費曼路徑定義為通過一系列的散射點；𝑡𝑃為經過路徑 𝑃 的到時；𝐴𝑃 為相對應 的振幅；𝑆(𝑡) 為跟震源相關的波場能量。(2.1.8 式) 並不能闡述波行路徑的散射點 狀態，僅能代表所有可能路徑總和疊加的總和效應。
當散射點或是震源位置受到極小的擾動且波場行進路徑距離遠小於平均自由 路徑 (mean free path) 時，則幾何擴散與散射強度的效應可被忽略，反映於波形上 的變化由每個路徑上的到時變化 𝑇𝑃 所主導：
𝑢𝑝𝑒𝑟(𝑡) = ∑ 𝐴𝑃 𝑃𝑆(𝑡 − 𝑡𝑃− 𝑇𝑃) . (2.1.9 式)
𝑅(𝑡,𝑇)(𝑡𝑠) ≡ ∫𝑡−𝑇𝑡+𝑇𝑢𝑢𝑛𝑝(𝑡′)𝑢𝑝𝑒𝑟(𝑡′+𝑡𝑠)𝑑𝑡′
√∫𝑡−𝑇𝑡+𝑇𝑢𝑢𝑛𝑝2(𝑡′)𝑑𝑡′∫𝑡−𝑇𝑡+𝑇𝑢𝑝𝑒𝑟2(𝑡′)𝑑𝑡′ , (2.1.10 式) 其中積分計算的時間段是根據路徑在波場被擾動前的到時 𝑡 為中心，設定前後增加
減少因擾動造成的到時變化量的時間段，即 2𝑇 的寬度；𝑡𝑠 則代表擾動波場相對非 擾動波場的遲滯時間。將 (2.1.8 式) 及 (2.1.9 式) 代入，則導出表不同路徑交互相 關的雙總和 (double summation), ∑𝑃𝑃′ 。若表示成矩陣運算形式，可視為代表一致 路徑的交相關性對角線項 (diagonal terms) (𝑃 = 𝑃′)，非對角線項為不同路徑的交互 相關性項 (cross terms) (𝑃 ≠ 𝑃′)。假設不同路徑的交互相關性為不一致 (incoherent) 的，其積分平均後對整體並無貢獻，可被忽略，則(2.1.10 式)可改寫為：
𝑅(𝑡,𝑇)(𝑡𝑠) ≈∑𝑃(𝑡,𝑇)∑ 𝐴𝑃2𝐴𝐶(𝜏𝑃−𝑡𝑠)
𝑃2𝐶(0)
𝑃(𝑡,𝑇) , (2.1.11 式) 其中 ∑𝑃(𝑡,𝑇) 代表所有具相同路徑的擾動和未擾動波場在有限時間窗交互相關函 數的總和。𝐶(𝜏) 為震源訊號的自我相關函數，表示為：
𝐶(𝜏) ≡ ∫ 𝑆(𝑡−∞∞ ′+ 𝜏)𝑆(𝑡′)𝑑𝑡′. (2.1.12 式) 若震源遲滯時間 𝜏 遠小於波場的主要週期時，可以二階的泰勒展開式
(second-order Taylor expansion) 表示𝐶(𝜏) ≡ 𝐶(0)(1 −12𝜔̅2𝜏2)，其中 𝜔̅2 代表時間窗 內到達的所有散射波相的頻率均方值。因此(2.1.11 式)可替換為：
𝑅(𝑡,𝑇)(𝑡𝑠) = 1 −12𝜔̅2〈(𝜏 − 𝑡𝑠)2〉(𝑡,𝑇) , (2.1.13 式) 其中〈… 〉(𝑡,𝑇)代表所有路徑到時的遲滯擾動量在時間窗內範圍 (𝑡 − 𝑇, 𝑡 + 𝑇) 的平均 值。當 𝑡𝑠 = 〈𝜏〉(𝑡,𝑇) 時 𝑅(𝑡,𝑇)(𝑡𝑠) 可達到最大值為：
𝑅𝑚𝑎𝑥(𝑡,𝑇) = 1 −12𝜔̅2𝜎𝜏2 , (2.1.14 式) 其中 𝜎𝜏2 為到時擾動的變異度 (variance)。因此可以藉由擾動前後的波形來估算到
時擾動的平均值以及變異度。
若假設擾動前後平均路徑長度沒有改變，則平均到時擾動 〈𝜏〉(𝑡,𝑇) = 0。受擾動 後散射點位置改變以均方根距離 𝛿 表示， 𝑛 為散射點的數量，則 𝑛 =𝑣𝑡𝑙，t 為波在 散射介質中行進的時間，𝑙 為散射點間的距離，則到時擾動的變異度可以下列表 示：
𝜎𝜏2 = 2𝛿𝑣𝑙2𝑡
∗ , (2.1.15 式) 其中 𝑙∗ 為散射點間傳遞的平均自由路徑 (transport mean free path)，代表傳遞過程 中已經失去震源以及方向資訊的距離。將 (2.1.14 式) 及 (2.1.15 式) 代入可得
𝛿2 = (1 − 𝑅𝑚𝑎𝑥(𝑡,𝑇))𝜔̅𝑣𝑙2∗𝑡 . (2.1.16 式) 對於固定的散射點來說，不同型態的擾動造成速度的改變量為一常數的值 𝛿𝑣，
則波傳遞的平均到時擾動 〈𝜏〉(𝑡,𝑇) = − (𝛿𝑣𝑣) 𝑡，只與波到時 𝑡 有關，與任一路徑是非
相關的，所以當時間窗足夠小時，𝜎𝜏 = 0。速度變化與交互相關函數後得到的最大 時間遲滯關係為：
𝛿𝑣𝑣 = −〈𝜏〉(𝑡,𝑇)𝑡 . (2.1.17 式)
圖 2-2、 聲波多重散射模擬示意圖。圖左實心與空心圓形分別代表擾動前後散射 點位置；星號代表聲源；三角形代表接收器。圖右實線與虛線分別代表擾動前後 的波形 [取自 Snieder et al., 2002]。
圖 2-3、 聲學實驗測試不同溫度下花崗岩樣本的速度變化。紅色與藍色分別代表 溫度在 50℃ 與 45℃ 記錄的波形。右上角圖示代表直達波波形，右下為尾波波形。
可以看出尾波波形走時偏差較直達波明顯，而且隨著到時有累積走時偏差的趨勢。
[取自 Snieder et al., 2002]。
2.2 資料收集及處理
2.2.1 儀器與連續資料
資料部分針對甲仙地震與南投地震分別選用震央距 60 公里內，由中央氣象局 即時觀測網 CWBSN (Central Weather Bureau Seismic Network) 短周期地震儀(Short Period Seismometer) 與台灣寬頻地震網 BATS (Broadband Array in Taiwan for Seismology) 寬頻地震儀 (Broadband Seismometers) 所記錄的連續資料 (圖 2-4) (表 2-1)。其中短周期地震儀使用 S13 速度型感應器，對於 0.1 Hz 以下的訊號較不 敏感；寬頻地震儀使用 STS-1、STS-2 等感應器，敏感範圍頻寬較寬，儀器響應如 (圖 2-5)。取樣點皆為每秒 100 點，分別使用涵蓋發震時間前後將近一年的資料量，
甲仙與南投地震使用測站詳細資訊分別見表 2-2、表 2-3。
2.2.2 經驗格林函數重建
首先將垂直分量的原始連續資料，經過移除時間序列上的平均值、線性趨勢 線以及儀器響應後得到位移的連續記錄。
為了達到較均勻噪訊源的環境場，將每天的連續資料以 600 秒為單位一次移 動 200 秒分別切割成數百個重疊的小視窗，再進行交互相關函數的運算 (圖 2-6)，
最後將每個視窗的經驗格林函數疊加平均可得單位為一日的交互相關函數。在進 行交互相關運算時，為了抑制地震訊號或是測站資料有突值的影響，在頻率域為 了使每個視窗頻譜分佈呈現白頻譜 (white spectrum)，我們不改變相位將振幅正歸 化 (whitening)。以及在時間域計算一整天的波形振幅的標準差，將超過 1.2 倍標準
2.2.3 參考與現時格林函數重建
周遭噪訊法是利用環境中許多隨機微弱的能量在空間中重複反射散射，經由 建設或破壞性干涉的作用之後，利用交互相關運算自看似隨機、無規律的波場記 錄中提取出有意義的訊號。不同於傳統地震學使用地震能量訊號，周遭噪訊能量 相對較低，因此在重建經驗格林函數之後必須再經由疊加的方式將區域性一致的 訊號彰顯出來。
我們將所有資料天數的經驗格林函數疊加起來作為反映長期平均兩站間底下 速度構造的參考格林函數 (Reference Green’s Function, RGF)，並為了能更穩定的測 得短時間格林函數的微小變化，我們計算疊加不同天數的 CCF 與長時間疊加整年 的 RGF 的相關係數 (圖 2-7)，在穩定波形與解析度之間取得最適當的疊加天數。
最後我們選擇疊加 30 天做為之後用來量測的的現時格林函數 (Current Green’s Function, CGF) (圖 2-8)。
表 2-1、 甲仙與南投地震資料與測站相關資訊。
EQ Network Instrument Sampling Rate Time
Jiasian CWBSN S13 100 sps 2009-9/2010-6
Nantou BATS BB 100 sps 2013-1/2013/10
表 2-2、 甲仙地震使用測站位置資訊。
測站名 經度(°𝐄) 緯度(°𝐍) 高程(m) 震央距(km) 天數
ALS 120.8133 23.5082 2413.4 60.4974 358 SGS 120.5906 23.0804 277.5 17.2891 367 SSD 120.6401 22.7443 148.3 25.9879 284 STY 120.7657 23.1607 639.7 21.8631 254 WTP 120.6223 23.2436 560.3 31.5844 343
CHN1 120.5284 23.185 360 30.12 368
CHN4 120.5939 23.3512 205 43.8318 350
TWG 121.0799 22.8176 195 41.51 365 TWL 120.5022 23.2638 590 38.8576 333 ELD 121.0251 23.187 1040 40.224 347
表 2-3、 南投地震使用測站位置資訊。
測站名 經度(°𝐄) 緯度(°𝐍) 高程(m) 震央距(km) 天數
SSLB 120.9540 23.7876 516 15.8229 312
TDCB 121.1583 24.2527 1308.4 40.5646 289
WUSB 121.1175 23.9919 1522.4 12.2732 298
YULB 121.2971 23.3925 376.1 56.9386 313
VWDT 121.1412 23.7537 2559 18.6674 309
CHGB 121.1740 24.0602 1846.5 21.7604 294
HGSD 121.4239 23.4921 104 59.0841 318
圖 2-4、 研究使用測站與地震後一個月餘震分布圖。深紅色與深藍色三角形分別 代表甲仙與南投地震使用的測站；橘色與淺藍色圓點代表地震後震央鄰近區域一 個月內的地震分佈；綠色星形代表震央位置。
圖 2-5、S13 感應器 (左) 與 STS-1 感應器 (右) 的速度型儀器響應圖。上圖為振幅 響應下圖為相位響應。S13 自然頻率為 1 Hz；STS-1 自然頻率為 0.0028 Hz。
圖 2-6、 上圖為單日連續資料切割視窗示意圖，每個視窗以 600 秒為單位並一次 移動 200 秒進行交互相關運算。下圖為單日所有視窗交互相關函數疊加平均後的 結果。
圖 2-7、 疊加不同天數現時格林函數 (CGF) 與長時間疊加之參考格林函數 (RGF) 相關係數關係圖。
圖 2-8、疊加 30 日的現時格林函數 (紅) 與長時間疊加之參考格林函數 (藍) 波形 圖。
表 2-4、 建置格林函數參數。
EQ Type Window Overlap Threshold Whitening Frequency
Jiasian CCF 600 s 67% S.D.(1.2) Yes 0.1-0.9 Hz
Nantou CCF 600 s 67% S.D.(1.2) Yes 0.1-0.8 Hz
2.2.4 尾波時間窗的選擇
尾波被視為在不均質速度構造中多次散射的結果，多重路徑效應使得尾波僅 隱含著介質本身的特性。在尾波時間窗選擇中先將正負時間軸經驗格林函數平均，
並以直達表面波後第一個最大振幅到時為起始時間，經由所選擇的尾波時間窗長 度範圍，以一秒為單位依序往後移動時間窗的起始時間，將時間窗內截取的波形 與相對應的時間窗得參考波形進行相關係數運算，最後可得到兼具穩定性與正確 性的尾波時間窗 (圖 2-9)。
甲仙地震根據 Wegler et al. [2009] 利用尾波干涉量測日本 𝑀𝑤6.6 新瀉地震前 後速度變化得到的經驗尾波時間窗長度 100 秒為主要使用的尾波長度。由於南投
為分界分別選取不同長度之尾波時間窗 (50 公里以上為 80 秒；以下 60 秒) (圖 2-10)，
相關資訊如表 2-5。
圖 2-9、 選取尾波時間窗方法示意圖。下圖紅色與藍色分別代表 RGF 與 CGF 時 間軸正方向波形；上圖紫色線條代表尾波時間窗的範圍，中間黑色圓點代表此尾 波時間窗所有 CGF 與 RGF 相對應的時間窗相關係數值。
圖 2-10、 選取尾波時間窗長度示意圖。上圖 SSLB-WUSB 站間距約為 28 公里，
使用 60 秒尾波時間窗長度；下圖 SSLB-YULB 站間距約為 56 公里，使用長度為 80 秒。
表 2-5、 南投地震使用尾波時間窗與測站間距。
測站名 尾波(起始) 尾波(結束) 相關係數 測站間距
CHGB-HGSD 85.71 165.71 0.7331206 70.36800385
CHGB-SSLB 34.7 94.7 0.8879341 37.58110428
CHGB-TDCB 19.83 79.83 0.9231488 21.29185867
CHGB-WUSB 17.62 77.62 0.8639932 9.49539471
CHGB-YULB 38.13 118.13 0.8029519 74.97919464
HGSD-SSLB 59.38 139.38 0.8998209 62.37511826
HGSD-TDCB 77.17 157.17 0.7284148 90.54434204
HGSD-WUSB 40.68 120.68 0.8522858 66.55834198
SSLB-TDCB 22.37 102.37 0.9239273 55.45619202
SSLB-WUSB 19.14 79.14 0.9177147 28.08648682
SSLB-YULB 22.96 102.96 0.9180021 56.02643967
TDCB-WUSB 23.56 83.56 0.8967490 29.08964348
TDCB-YULB 46.42 126.42 0.8524396 96.19794464
WUSB-YULB 31.38 111.38 0.9459124 68.84261322
CHGB-VWDT 20.94 80.94 0.9123160 34.09152222
HGSD-VWDT 28.32 88.32 0.8675210 44.80438232
SSLB-VWDT 13.94 73.94 0.9292399 19.43283081
TDCB-VWDT 31.33 111.33 0.8612059 55.18714142
VWDT-WUSB 27.22 87.22 0.8865389 26.47696877
VWDT-YULB 18.89 78.89 0.9093367 43.04622650
2.3 量測相對走時偏移方法
假設地殼介質速度擾動量 (𝛿𝑣𝑣) 均勻分布(homogeneous)，隨著波散射距離越長，
相對走時變化量 (𝛿𝜏) 就越大。也就是尾波走時偏差會隨著到時有線性增加的趨勢，
其負斜率值等同於速度變化的比值，𝛿𝜏𝜏 = − (𝛿𝑣𝑣) (詳見 2.1.2 節)。
經由給定不同參數的尾波時間窗長度與位置，量測長時間的 RGF 與短時間 CGF 的尾波走時偏差值，再根據尾波本身的到時 τ 做線性回歸分析即可得到地殼 速度隨時間變化的關係。
量 測 兩 波 形 訊 號 相 對 走 時 差 或 相 位 差 最 常 使 用 的 方 法 為 交 互 相 關 法 (Cross-correlation) ， 此 方 法 在 訊 號 處 理 上 是 用 來 量 測 兩 者 波 形 的 相 似 程 度 (Similarity)，透過不斷改變其中之一的訊號偏移時間 𝛿𝜏，來計算兩者波形的相關係 數 (Correlation coefficient)。當相關係數達到最大值，其所對應到的偏移時間即為 相對走時差。本研究使用移動視窗交叉頻譜法與拉張法，做為計算速度擾動量的 方法，將在以下分別介紹。
2.3.1 移動視窗交叉頻譜法
Poupinet et al. [1984] 提 出 移 動 式 窗 交 叉 頻 譜 法 (Moving Window Cross Spectral Method, 以下簡稱 MWCSM) ，用來計算重複地震體波或尾波的走時偏差 值，將各頻率含量納入考量並計算誤差，以改善量測時的準確度。
假設 𝑥(𝑡) 與 𝑦(𝑡) 分別代表 CGF 和 RGF 所選取時間窗的尾波訊號，則相對偏 移時間之交互相關函數 (Cross-Correlation Function, CCF) 為：
𝐶𝑥𝑦(𝛿𝜏) = ∫ 𝑥(𝑡)𝑦(𝑡 + 𝛿𝜏)𝑑𝑡 . (2.3.1 式) 將 𝑥(𝑡) 與 𝑦(𝑡) 經傅利葉轉換公式 (Fourier Transform) 至頻率域之頻譜分別
以 𝑋(𝑓) 與 𝑌(𝑓) 表示，則交叉頻譜 (cross-spectrum) 為：
𝑆𝑥𝑦(𝑓) = 𝑋(𝑓)𝑌∗(𝑓) . (2.3.2 式) 而 𝑆𝑥𝑦(𝑓) 在特定頻率𝑓𝑖的相位 𝜑𝑖 與相對偏移時間 𝛿𝜏 的關係可表示成：
𝜑𝑖 = 2𝜋(𝛿𝜏)𝑓𝑖 , (2.3.3 式) 將 𝜑𝑖 對 𝑓𝑖 的變化以通過原點的直線作線性迴歸所得的斜率與 2𝜋 的比值則為 兩訊號之相對走時偏移 𝛿𝜏。在計算斜率時，考慮到不同頻率下波形的一致性造成 不同誤差範圍值，進而影響估算的走時偏差(圖 2-11a)。因此在迴歸斜率時以兩者 波 形 在 不 同 頻 率 𝑓𝑖 下 的 一 致 性 (Coherency) 做 為 所 對 應 相 位 值 𝜑𝑖 的 權 重 𝑤𝑖 [Brenguier et al., 2008a]：
𝑤𝑖 = |𝑆𝑥𝑦1−𝐶(𝑓𝑖)|𝐶𝑖2
𝑖2 , (2.3.4 式) 其中𝐶𝑖 =|𝑆𝑥𝑦(𝑓𝑖)|
√𝑆𝑥𝑥𝑆𝑦𝑦 為波形在不同頻率 𝑓𝑖 下的一致性，為交叉頻譜平滑之後的振幅
|𝑆𝑥𝑦(𝑓𝑖)| 與兩者自我頻譜乘積的平方根 (Auto spectrum) 比值。兩訊號走時偏移的 誤差值為：
𝐸𝛿𝜏 = √∑𝑀𝑖=1𝑤𝑖2(𝜑𝑖−𝑎𝑓𝑖)2
∑𝑀𝑖=1𝑤𝑖2𝑓𝑖2 , (2.3.5 式) 其中 𝑎 為線性迴歸交互頻譜相位 𝜑𝑖 對頻率 𝑓𝑖 關係之直線斜率 2𝜋(𝛿𝜏)；𝑀為線性迴 歸時涵蓋的頻率數量。
MWCSM 是以移動視窗的方式來量測 CGF 與 RGF 的相對走時偏移，根據 (2.3.3 式)分別估算正負軸上每個移動視窗內尾波的走時偏差值 𝛿𝜏，再根據每個移 動視窗所估算之走時偏差值給予權重值 𝐸𝛿𝜏 進行線性迴歸，求得通過零點的直線斜
其負值即為速度擾動量 𝛿𝑣𝑣 (圖 2-11b)。其速度擾動值的誤差 𝐸𝛿𝜏/𝜏 為
𝐸𝛿𝜏/𝜏 = √∑𝜎𝛿𝜏2𝜏
𝑖2
𝑁𝑖=1 , 𝜎𝛿𝜏2 = ∑𝑁𝑖=1(𝛿𝜏𝑖− 𝛿𝜏̅̅̅̅)𝑖 2/𝑁 , 𝛿𝜏̅̅̅̅ =𝑖 𝛿𝜏𝜏 × 𝜏𝑖 ,
其中 𝑁 為移動視窗個數， 𝜎𝛿𝜏 為走時偏移 𝛿𝜏 的標準差， 𝛿𝜏𝑖 與 𝛿𝜏̅̅̅̅ 分別代表第 𝑖 個𝑖 移動視窗的走時偏移量和迴歸斜率的預測值，標準差為兩者平方差總和平均後的 平方根。
2.3.2 拉張法
拉張法 (Stretch Method) 是透過輸入速度擾動量範圍，以線性內插方式拉長或 壓縮時間軸以改變 RGF 波形後，與短時間現時 CGF 做交互相關函數，
直至兩者達到最佳擬合程度，也就是相關係數最大值 [Sens Schönfelder and Wegler, 2006] (圖 2-12)。
進行測量前，須將正負時間軸訊號平均，長期平均之 RGF 代表介質速度擾動 前狀態的格林函數，其尾波訊號以 𝑦(𝑡) 表示；當介質速度發生均勻微小的擾動量 𝛿𝑣𝑣， 造成短時間 CGF 尾波 𝑥(𝑡) 走時偏移，預期改變的尾波訊號可以寫成 [Duputel et al., 2009]：
𝑥(𝑡) = 𝑦 (𝑡 (1 +𝛿𝑣𝑣)) . (2.3.6 式)
當速度擾動為正時介質速度變快，CGF 相對 RGF 較快則時間軸被拉長；反之 速度變慢，CGF 時間軸被壓縮。
由於噪訊紀錄為離散訊號，在計算拉長或壓縮 RGF 和 CGF 的相關係數時必須 重新內插資料點，使每一點振幅所對應的時間點一致。若 RGF 到時為 𝜏𝑖 的尾波受 速度擾動 𝛿𝑣𝑣 產生走時偏移至新的時間點 𝜏𝑗 = 𝜏𝑖(1 +𝛿𝑣𝑣)，則被拉長或壓縮的 RGF 在未擾動前 𝜏𝑖 時間點的振幅由其在時間上相鄰兩點 𝑗 + 1 和 𝑗 (時間分別為 𝜏𝑗 和
𝜏𝑗+1 ) 所對應的振幅線性內插求得：
x(𝜏𝑖) = 𝑦𝑗+ (𝑦𝑗+1− 𝑦𝑗)/𝛿𝑡[𝜏𝑖(1 +𝛿𝑣𝑣) − 𝜏𝑗] . (2.3.7 式) 其中 𝛿𝑡 為採樣時間間隔 (0.01s)，每一組波形以 0.001 % 為單位，微小改變的速 度變化量 𝛿𝑣𝑣 所產生拉張或壓縮後的波形進行 600 次運算，測試之速度擾動範圍為
−0.3 % ~ 0.3 % 。
圖 2-11、使用合成之 CGF(黑)及 RGF (紅)，預設速度擾動量為−1%時示範移動視 窗交叉頻譜法量測走時偏差操作過程。(a)由上至下再由左至右分別為 12-28 秒 (標 示於(b)的藍色虛線區間內) 兩波形在不同頻率的一致性 (coherency，𝑐𝑖)；於頻率 域對相位 𝜑𝑖 做線性回歸得到相對走時差(直線斜率)；兩波形交互頻譜 𝑆𝑥𝑦 對不同頻 率振幅含量；利用波形在頻率含量的一致性以及交互頻譜所計算的權重 𝑤𝑖，四圖 單位皆無單位因次，橫軸單位為赫茲(Hz)。(b)上圖放大標示 CGF 和 RGF 尾波訊號，
左側與右側分別為負與正時間軸。中圖為完整包含直達波經驗格林函數波形。下 圖為以 17 秒視窗長度進行移動視窗量測之走時偏差，利用最小平方差進行線性回 歸得到速度變化，其負斜率值相當於速度變化量，由圖中可得相對速度擾動接近 −1 % [余, 2011]。
圖 2-12、 拉張法量測走時偏差操作過程。使用合成之 CGF (黑) 及 RGF (紅)，
測試用拉張法量測的相對走時偏移和速度擾動量。左圖為在不同速度擾動量下對 RGF 的時間軸拉張或壓縮之後的波形和 CGF 的比較。右圖為對應的波形相關係數，
在速度擾動量為 −1 % 時相關係數至最大值 [余, 2011]。
EQ Time
EQ-reference
Start time-reference
第三章 研究結果
本研究利用移動視窗交叉頻譜法與拉張法量測短時間與長時間經驗格林函數 沿著時間軸的速度擾動，為了使量測結果兼具穩定及可信度，必須對經驗格林函 數進行疊加以獲取穩定波形。隨著疊加天數遞增，對於日期的解析度亦隨之下降，
於是經過交互相關函數的分析 (本文 2.2.3 節)，本研究最後選擇疊加 30 天經驗格 林函數波形並將參考日期以一天為單位移動，來得到單日的解析度。
在疊加經驗格林函數的過程中，我們使用兩種疊加方式，藉以用不同的觀點 來討論其量測結果的物理量 (圖 3-1)。首先為沿著時間軸疊加，讓經驗格林函數順 著時間軸進行疊加，適合用來監測較長期間的速度擾動，或是週期性變化以及與 其他地體活動事件的相關性；再來是以事件發生日期為參考點疊加，主要是要探 討此事件對週遭影響的程度以及在時間上影響的長度。圖 3-2 比較了兩種疊加方法 應用在甲仙地震的結果，可看出兩種方法皆觀測到地震造成時間軸上的速度擾動，
照時間疊加得到較平滑的數據結果， 依地震疊加則得到由地震影響明確的相對速 度擾動值。由於本研究探討地震活動造成鄰近區域淺部地殼的速度擾動，後面章 節討論以依地震參考日期疊加的成果為主。
圖 3-1、 經驗格林函數疊加方法，上圖為依事件日期為參考點疊加，下圖為順時 間軸疊加。
圖 3-2、 SSD-WTP 站兩種疊加方法量測結果比較。上圖為照時間疊加；下圖為照 地震發生日期疊加，黑色虛線為甲仙地震日期，藍色資料點為移動視窗交叉頻譜 法結果；紅色為拉張法結果。。
3.1 甲仙地震前後速度變化
經由 MWCSM 以及拉張法量測，甲仙地震前後有明顯的同震速度下降 (圖 3-3)。
同時我們計算了地震前後一個月的平均速度擾動差值，用以代表地震前後穩定的 速度變化並賦予在統計學上的意義。由圖 3-4 可看出在統計上甲仙地震前後有明顯 的速度降低，其最大值可至將近 −0.1 %。量測到速度降低的測站對主要穿越於震 央西方數公里處，以區域構造位置與速度變化的分佈圖來看 (圖 3-5)，此區位於東 傾逆衝斷層上盤位置以及餘震活動頻繁區域，暗示著甲仙地震期間可能有明顯的 應力環境改變或是受到地震造成構造上的破壞。
Yu and Hung [2012] 利用周遭噪訊的高頻單站自我相關函數 (Auto-correlation function, ACF) 來量測 2006 年台東地震前後速度擾動，發現地震事件後有明顯的 同震速度下降以及地震發生數天前經驗格林函數波形的相關係數值 (Correlation Coefficient, CC) 開始有明顯降低的情形；Obermann et al. [2013] 利用周遭噪訊交 互 相 關 函 數 來 監 測 火 山 活 動 期 間 有 明 顯 的 速 度 擾 動 以 及 波 形 的 退 相 干 性 (Decoherency)。因此我們在使用拉張法量測時，也同時針對波形的一致性進行比 較與探討，由圖 3-3 可看出 SSD-WTP 測站對在甲仙地震影響前一個月，交互相關 係數值開始有些微的降低情形。本研究利用雙站 CCF 的結果是根據 30 天疊加平均 的結果，相對 Yu and Hung [2012] 單站三天疊加結果解析度較差，可能可以解釋 波形退相干性相對台東地震前較不明顯；CHN1-SSD 則顯示震前波形的一致性有 明顯的降低。圖 3-6 為甲仙地震前後 2 日移動視窗交叉頻譜法的量測結果，展示主 震之後尾波到時遲滯現象，並且在 CCF 正負到時時間軸成反對稱以及地震前後線 性回歸的斜率變化。
圖 3-3、 甲仙地震 SSD-WTP 站 (上圖) 與 CHN1-SSD 站 (下圖) 結果。藍色資料 點為移動視窗交叉頻譜法結果；紅色為拉張法結果；黑色資料點為拉張法交互相 關係數值；垂直虛線為主震發生日期。
圖 3-4、 甲仙地震 MWCSM 結果。垂直粗黑線表示主震發生日期；垂直虛線表示 地震前後一個月日期；紅色直線表示地震前後一個月的平均速度擾動；灰色區域 為地震前後一個月速度標準差範圍；數據點顏色表示 MWCSM 計算線性回歸方均 根誤差量值；右上角數字以誤差值為權重計算的前後一個月平均速度擾動以及標 準差範圍。
Drop (%)
圖 3-5、 甲仙地震拉張法量測結果與區域地表構造的分佈。連線部分為地震前後 一個月平均速度變化量值，暖色系代表速度下降；冷色系為速度上升，綠色圓點 為主震後一個月的地震分佈。
圖 3-6、 2010/3/4 甲仙地震 SSD-WTP 測站對 3/2 震前 (上圖) 與震 3/6 後 (下圖) MWCSM 量測結果，圖說參見 (2.3.1 節)。
3.2 南投地震前後速度變化
3.2.1 3 月 27 日南投地震結果
3 月南投地震也量測到明顯的同震速度下降 (圖 3-7)，在統計上計算地震前後 一個月平均速度擾動差值也觀察到震後平均速度的降低 (圖 3-8)。我們將兩種量測 方法得到的數值進行整理比較 (表 3-1)，並利用兩方法量測的一致性作為結果穩定 性的標準。將結果與地表構造分佈圖比較，發現主要同震速度降低測站對大致穿 越震央鄰近區域，同時也是東傾逆衝斷層上盤以及餘震活動區域 (圖 3-9)。
為了確認量測結果的可信度，同時比較拉張法量測的交互相關係數值 (圖 3-7) 以及 MWCSM 在地震前後 1 日的計算過程 (圖 3-10)。測量時拉張法波形的相關係 數值幾乎皆大於 0.6，並且在 3 月地震影響期間有些微的下降。圖 3-10 展示 MWCSM 量測結果，顯示 SSLB-WUSB 測站對線性迴歸的誤差值較小，以及在地 震前後明顯的斜率值變化。
圖 3-8、南投地震 MWCSM 結果，以 3 月主震日期做為疊加參考日期。
圖 3-9、 3/27 南投地震拉張法量測結果在區域地表構造上的分佈。
Velocity Drop (%)
圖 3-10、 2013/3/27 南投地震 SSLB-WUSB 測站對 3/26 震前 (上圖) 與 3/28 震後
3.2.2 6 月 2 日南投地震結果
6 月南投地震也量測到明顯的同震速度下降 (圖 3-11)，圖 3-12 結果也顯示震 後平均速度的降低。將兩方法量測的一致性作為結果穩定性的標準進行整理比較 (表 3-2)，發現同震速度擾動結果與 3 月地震雖不完全相同，但主要同震速度降低 測站對大致穿越震央鄰近區域，同時也是東傾逆衝斷層上盤以及餘震活動區域 (圖 3-13)，暗示著震源機制相近兩地震造成區域淺部地殼彈性介質影響的相似性。
為了確認量測結果的可信度，同時也比較拉張法量測的交互相關係數值 (圖 3-11) 以及 MWCSM 在震前與震後計算過程 (圖 3-14)，結果顯示 TDCB-WUSB 在 地震前後量測的誤差較小；圖 3-11 可看到 6 月地震影響期間相關係數值有明顯的 降低，並且約在震前半個月就有下降的趨勢。
圖 3-11、 南投地震期間速度變化結果，以 6 月主震日期做為疊加參考日期，黑色 虛線為 3 月與 6 月主震發生日期。
圖 3-12、南投地震 MWCSM 結果，以 6 月主震日期做為疊加參考日期。
city Drop (%)
圖 3-14、 2013/6/2 南投地震 TDCB-WUSB 測站對 5/31 震前 (上圖) 與 6/4 震後 (下 圖) MWCSM 量測結果。
