利用群組地震法探討阿拉斯加隱沒帶的震源特性

(1)

國立臺灣大學理學院海洋研究所碩士論文

Institute of Oceanography College of Science

National Taiwan University Master Thesis

利用群組地震法探討阿拉斯加隱沒帶的震源特性 The Seismic Source Characteristics in the Alaska Subduction Zone Determined by the Cluster Event Method

羅仲宏 Chung-Hung Lo

指導教授：郭本垣博士喬凌雲博士

Advisor: Ban-Yuan Kuo, Ph.D.

Ling-Yun Chiao, Ph.D.

中華民國 107 年 7 月

(2)

(3)

致謝

由衷地感謝幫助我寫出這本論文的所有人。回首兩年前，我初來臺大海洋所，

完全沒有把握自己可以在兩年內順利畢業拿到學位，但是我做到了！我特別要感謝兩位指導教授：郭本垣研究員的氣場一開始讓我震懾，但隨著對您的認識加深，

震懾逐漸轉為敬佩，您總是循循善誘，讓我在迷途中找到正確的道路；喬凌雲教授在海洋所的鼓勵，讓我在 seminar 答辯時，能夠突破重圍、重拾信心。

感謝梁文宗老師在暑期實習對我的啟發並同陳伯飛老師與陳勁吾老師引領我進入地震學的殿堂。感謝柯彥廷老師雖然遠在 Caltech 卻還是不辭辛勞在程式方面提供協助。謝謝曾泰琳老師與黃信樺老師提供許多意見讓我的論文更嚴謹完整。

謝謝龔源成老師、洪淑蕙老師與馬國鳳老師在地震學觀念的教導。我很幸運能夠有許多老師在這條道路上給我幫助，讓我不怕面對困難，收穫滿滿！

感謝在這條路上陪我打拼的眾多親友們。榮均與乃慈常在我不知所措的時候陪伴我，給我繼續走下去的勇氣。遠在集集的如萱幹事姐姐總是用充滿溫度的文字鼓勵我。實驗室的品瑜包辦一切瑣事，讓我可以專心在研究上，不被雜事纏身。

政謙學長、凱勛學長、棟邦、立宸、汶珮是我討論的好夥伴。煒晴帶我理解屬於臺大的強悍、浪漫並身兼吃飯好朋友。全變的晉丞、翔孟、雅琪；地球所的逸威學長、岱潔學姐、志煌學長、冠芙學姐、芳儀學姐、怡君、芷芹；去年畢業的志銘、志傑、至為、芝吟、冠羽；中大天文社的宛庭、立宇、珮欣、鈞凱、佩娟、

哲豪、雲澤；海洋所 R05 的同學們，以上貴人陪伴我走過在研究所的艱辛路途。

暑期生的毓芳、奐鈞、昇均、辰叡、皓哲讓我過足了教學的乾癮。還有許多礙於篇幅未能記錄的朋友們與柴犬迪迪，謝謝你們讓我有個多采多姿的臺大生活。

十萬萬分感謝我的金主爸媽與弟妹，讓我無後顧之憂地完成到目前為止的求學之路。

最後，給翻開本碩士論文的你，我盡量以淺白的方式述說群組地震法與相關的演算法，希望對你的研究有所幫助。

仲宏謹誌於中央研究院

(4)

摘要

斷層在中深部 (50 - 300 km) 高岩壓環境下的破裂機制目前尚有爭論。溫度、

壓力與區域應力場的差別，均會對破裂機制產生不同的影響，進而決定斷層的破裂行為。連續深度剖面的震源參數提供了解不同深度破裂機制的機會。本研究分析地震波頻譜，計算深層地震包含拐角頻率、地震矩、應力降、地震波輻射能量、

地震波幅射效率在內的震源參數，嘗試為深層地震的發生機制提供約束。

本研究利用群組地震法大幅減少逆推過程中的未知參數，有效分拐角頻率與路徑衰減對頻譜的影響。群組地震法的核心概念為建構多重震源-測站對能嚴謹地 定義共同路徑 Q 與各自震源的 fc，由觀測的資料同時決定震源與路徑衰減，以解決在頻譜上拐角頻率與路徑衰減之間的非線性關係。

本研究擷取 IRIS 資料中心 144 個分布於阿拉斯加隱沒帶的地震，規模介於 3 到 5 之間，震源深度 90 到 170 km，並計算中深部地震動態與靜態的震源參數。結果顯示，拐角頻率與地震矩呈現清楚的逆相關，鞏固了逆推架構與結果的可信度。

而拐角頻率、應力降、地震波輻射能量等各個震源參數的計算結果，符合地震理論的自我相似性規範。

阿拉斯加隱沒帶的震源特性具有高應力降 (18.34 ± 1.1 MPa) 與低地震波輻射效率 (0.27) 的特徵。暗示或許在阿拉斯加隱沒帶，深層地震的發生機制主導於熱失控模型。脫水脆化作用釋放出的水可能為熱失控模型的催化劑，而在地震波研究上僅看得出熱失控模型的特徵。本論文認為兩者是相輔相成的關係，共同為中深部地震的發生機制。

關鍵字：群組地震法；鄰域演算法；震源參數；中深部地震破裂機制

(5)

Abstract

The Alaska subduction zone is known as an active subduction zone where earthquakes occur at high rate and widespread along the trench and down dip of the slab.

In recent years, increasing numbers of seismic arrays have been installed in Alaska region that significantly increased aperture of observations. It provides an opportunity to systematically investigate depth-varying seismic source characteristics in subduction zone. The source characteristics can reflect the rupture behaviors of fault, but the differences in rupture behaviors as a function of depth are still notably debated. In this study, we employed a cluster-event method (CEM) to constrain the source parameters as well as along-path attenuation in the Alaska subduction zone. Neighborhood algorithm is applied to solve the nonlinear inverse problem. Using 40 stations from IRIS data management center, we analyzed 144 Alaska local earthquakes spreading over a depth range from surface to several hundred kilometers and a seismic magnitude range from 3 to 5 in 2012~2017. These source parameters are then converted to stress drop and radiated energy at different depths. The fc’s satisfy a self-similar scaling relationship with seismic moment of 𝑓_𝑐 ∝ 𝑀₀⁻³ with a mean stress drop of 18.34 ± 1.10 MPa in Madariaga’s form (Vs model). The lower radiation efficiency and higher stress drop might imply the shear heating instability and dehydration embrittlement as the same important faulting mechanisms for intermediate-depth earthquakes.

Keywords：cluster event method；neighborhood algorithm；source parameter；

faulting mechanisms for intermediate-depth earthquakes

(6)

圖目錄

圖 1.1、fc-t* trade off 說明圖。(a) fc對頻譜影響說明。(b) t*對頻譜影響說明。[Ko et

al., 2012] ... 3

圖 1.2、阿拉斯加區域構造圖。 ... 4

圖 2.1、群組地震法相關描述。(a)群組地震法說明示意圖。[Ko et al., 2012] (b)跨群組事件對 CEM 提供約束。(c)群組地震法與其他方法在逆推自由度的比較。 5 圖 2.2、鄰近震源的波傳路徑示意圖。 ... 6

圖 2.3、模型-資料殘差值分布比較。(a)單一震源-測站對殘差值圖。(b)群組震源- 測站對殘差值圖。 ... 7

圖 2.4、地震頻譜組成示意圖。 ... 8

圖 2.5、fc - t*- Ω0與頻譜的關係。[Shearer P., 2009] ... 9

圖 2.6、地震波形、頻譜展示圖。(a)地震波形圖。(b)訊號頻譜。... 10

圖 2.7、鄰域演算法示意圖。[Sambridge, 1999a] ... 11

圖 2.8、NA 迭代與殘差值的關係。 ... 12

圖 2.9、NA 的收斂情形。 ... 12

圖 3.1、本研究區域 (a) 地震及 (b) 測站分布圖。 ... 13

圖 3.2、資料前處理與資料篩選流程圖。 ... 15

圖 3.3、利用 AIC 自動挑 P 波到時的說明。(a)應用 AIC 法之範例。(b)P 波到時前後各 2 s 的波形與 AIC 值。 ... 16

圖 3.4、波形展示圖。(a)SNR 訊號及噪訊位置圖。(b)P 波及尾波位置圖。 ... 16

圖 3.5、不同濾波頻段對 SNR 的影響。(a)位移波形。(b)(c)(d)為 f1遞增的濾波效果對 SNR 的影響。 ... 17

圖 3.6、利用 SNR 與 PCR 篩選波形的例子。 ... 18

圖 3.7、目視篩選例子。 ... 18

圖 3.8、距離篩選示意圖。 ... 19

圖 3.9、本論文群組組成說明。 ... 20

圖 3.10、對數軸均勻採樣說明。 ... 20

圖 3.11、權重說明。。 ... 21

圖 3.12、Ω₀定義說明，同一個震源-測站對 (λ 搜尋範圍相同)，兩個不同的 X 模型 點。(a)f_c- t*組合： f_c = 5.71，t* = 0.062。 (b)另一組 f_c- t*組合： f_c = 2.33， t* = 0.31 的 Ω₀搜尋過程。 ... 22

圖 4.1、Variance Reduction 之統計圖。 ... 24

圖 4.2、CEM-NA 之逆推範例，範例群組的震央位置。 ... 24

圖 4.3、CEM-NA 之逆推範例，頻譜擬合。 ... 25

圖 4.4、CEM-NA 之逆推範例，全域搜尋的殘差值分布。 ... 26

圖 4.5、跨群組事件的 f_c。 ... 27 圖 4.6、阿拉斯加地區一維物理參數圖。(a)Vp、Vs 隨深度的變化。(b)密度隨深度

(9)

的變化。虛線之間為本論文震源深度的範圍。 ... 28

圖 4.7、阿拉斯加地區近震規模 (ml) 與震矩規模 (mw) 轉換圖。[Ruppert and Hansen, 2010] ... 29

圖 4.8、阿拉斯加地區地震目錄 M₀(catalog M₀) 與頻譜擬合 M₀(average M₀) 的關係圖。 ... 30

圖 4.9、阿拉斯加地區平均 M₀與 f_c的關係圖。 ... 31

圖 4.10、阿拉斯加地區應力降統計圖。 ... 31

圖 4.11、E_R/M₀與 M₀的關係圖。 ... 33

圖 4.12、比例能量的計算值與理論值關係圖。 ... 33

圖 4.13、應力與滑移的能量分布示意圖。[Shearer P., 2009]... 34

圖 4.14、阿拉斯加地區輻射效率圖。(a)統計圖。(b)地震波輻射效率的計算值與理論值關係圖。... 35

圖 5.1、高維度問題之例子。 ... 38

圖 5.2、約束不足問題示意圖。 ... 39

圖 5.3、約束不足問題例子。 ... 40

圖 5.4、NA 不收斂之情形。 ... 40

圖 5.5、ST 逆推方法測試圖。[Ko et al., 2012] ... 41

圖 5.6、ST 近似問題之範例。(a)波形、頻譜擬合圖。(b)全域殘差值分布圖。 (c)CEM-NA 的模型點迭代分布圖。 ... 42

圖 5.7、權重測試圖(W_i,max= 2)。 ... 43

圖 5.8、權重測試圖(W_i,max= 4)。 ... 44

圖 5.9、兩種權重測試頻譜擬合圖。 ... 44

圖 5.10、兩種權重測試，跨群組 f_c圖。 ... 45

圖 5.11、2 種 f_c上限值測試。 ... 46

圖 5.12、隱沒帶的應力降隨深度的變化圖。(a)阿拉斯加 (b)日本。[Ko and Kuo, 2016] ... 47

圖 5.13、阿拉斯加隱沒帶的應力降與比例能量的關係。 ... 48

圖 5.14、隱沒帶比例能量隨深度的變化圖。(a)阿拉斯加 (b)日本。[Ko and Kuo, 2016] ... 49

圖 5.15、比例能量理論與計算值的差別。 ... 50

圖 5.16、(a)阿拉斯加與(b)日本隱沒帶地震波輻射效率隨深度變化圖。[Ko and Kuo, 2016] ... 51

圖 5.17、阿拉斯加的 Q 值趨勢-震源與測站位置俯視圖。 ... 54

圖 5.18、阿拉斯加的 Q 值趨勢-AA’剖面與 BB’剖面。 ... 54

圖 A.1、2 維 Voronoi cell。 ... 62

圖 A.2、(a)馬氏鏈卡通圖(b)機率轉移矩陣 P (c)矩陣 P 的迭代次數結果列表。 .... 63

圖 A.3、加入接受率參數α後的機率轉移矩陣示意圖。 ... 66

(10)

圖 A.5、在 2 維參數空間中 Voronoi cell 的分布情形。 ... 71

圖 A.6、ns、nr 的數值使用示意圖。 ... 72

圖 B.1、窗函數概念圖 ... 73

圖 B.2、不同的窗函數。 ... 74

圖 B.3、主瓣與側瓣示意圖。(a)δ 函數與近似函數。(b)窗函數頻譜取絕對值示意圖。 ... 75

圖 B.4、不同階數 (k) 的橢球窗函數。... 77

圖 C.1、本論文選用的震源與測站分布圖。... 83

表目錄

表 5.1.1、極端群組組成示意表。... 36

表 5.1.2、極端群組組成步驟。 ... 37

表 5.4.1、阿拉斯加的 Q 值趨勢。[Stachnik et al, 2004]... 53

表 C. 1、選用地震資料。 ... 78

表 C. 2、選用測站資訊。 ... 82

(11)

第一章緒論

1.1 中深部地震成因之假說

大部分的地震肇因於斷層錯動。當施加在斷層上的累積應力超越了摩擦力，

斷層便快速滑動，而產生地震。此種作用稱為脆性變形 (brittle failure) [e.g_. Hoek.1968]。當深度增加，壓力也隨之增加，到某一深度後斷層破裂需要其他輔助

機制方能突破岩壓。這個深度通常在 50 ~70 km，在此深度之下的地震，就稱為中深部地震 (Intermediate-Depth Earthquake, IDE)。這些機制目前還未被完全了解，

不過許多假說已被提出 [Green and Houston, 1995 and Prieto et al., 2013]，以下介紹 兩種關於中深部地震發生機制的假說：脫水脆化作用 (dehydration embrittlement) 與熱失控模型 (thermal runaway)。

在隱沒過程中，含水礦物的脫水作用會提高岩石的孔隙壓力，例如玄武岩向榴輝岩的轉變會釋放出水。這個作用使得岩層的正向力減少，進一步降低斷層之 間的摩擦力，創造中深部斷層滑動的契機，這個過程稱為脫水脆化作用 [Green and Houston, 1995]。從實驗室的岩石力學實驗 [e.g., Dobson et al., 2002; Jung et al., 2004]、板塊內地震發生位置與地幔中變質反應位置相符等證據支持這個假設 [Peacock and Wang, 1999; Peacock, 2001; Hacker et al., 2003]。

另一個假說是熱失控模型 [Prieto et al., 2012 ; Prieto et al., 2013 ; Thielmann et al., 2015 ; Thielmann,2017]。在板塊隱沒之前，大量的正斷層發育，在斷層附近的 岩石粒徑受摩擦破壞而變小。海洋岩石圈隱沒之後，這些粒徑小的物質在中深部地區增加斷層滑動的可能性。在高岩壓的環境下，斷層一旦錯動，瞬間產生的高熱會使斷層周圍的岩石融化，斷層之間的摩擦力將因此減低而讓滑動加速，進而產生更多的熱，呈現正向循環。這些融化的岩石會形成假玄武岩玻璃 (pseudotachylyte)。這些假玄武岩質的岩石被擄獲岩帶至地表，在法國的科西嘉島 (the Alpine subduction complex of Corsica) 上被大量發現 [Anderson and Austrheim,

(12)

2006; Deseta et al., 2014]，成為支持此假說的證據。另一個證據是中深部地震的波 幅射效率 (radiation efficiency) 比淺層地震明顯的低 [Prieto et al., 2013]。地震波幅 射效率下降，起因於斷層錯動之後產生的地震波受到瞬間高溫的影響，震波能量因而快速衰減 (thermal dissipation)。

1.2 拐角頻率與時間的制衡關係 (f_c-t* trade-off)

斷層錯動在震源頻譜上最明顯的特徵為拐角頻率 (f_c)。f_c在頻譜中指的是頻率 f > fc之後，其振幅會快速衰減的門檻值。fc與斷層的破裂尺度有關 [Haskell, 1964;

Brune, 1970; Madariaga, 1976; Beresnev 2002]。因此，我們可以透過 f_c一窺斷層的破裂行為。然而，測站接收到的地震波訊號包含波傳遞的路徑效應，而路徑效應可用e^−πft^∗描述，t* 反映介質的衰減特性。在相同 t* 的條件下，越高頻的訊號振幅衰減越多。震源效應與路徑效應在地震波頻譜上表現相似，兩者的成因卻截然 不同：f_c 純粹表現震源破裂行為；t* 反映構造對波傳能量的影響。如果不能有效 分離 fc 與 t* 對頻譜的效應，其結果勢必會嚴重影響震源分析的可靠程度(圖 1.1)。

為了解決 f_c 與 t* 之間的制衡關係 (fc-t* trade-off)，常見的做法有兩種。其一的策略是使用特定工具來消除其中一個效應，如振幅雙站法 (Amplitude ratio method) 使用兩個測站消除震源效應 [e.g. Roth et al., 1999; Shito et al., 2004]；經驗 格林函數法 (Empirical Green’s function) 使用規模一大一小的震源消除路徑效應 [e.g. Mueller, 1985; Viegas et al., 2010]。雖然特定工具能夠減輕制衡關係的影響，

但受限於參考站或參考震源的位置，無法探知大規模區域的物理參數變化。其二 的做法則是Stachnik et al. (2004) 利用逆推法計算阿拉斯加地區全部測站頻譜的 fc， 再將不同測站的 fc 結果取平均，個別測站的 t* 計算使用平均的 fc，然而平均的 fc

是否就代表真實的 fc 亦有爭論。本論文使用Ko et al. (2012) 提出的群組地震法 (cluster event method, CEM)，調查阿拉斯加地區的 fc 與 t* 之間的制衡關係，勾勒出阿拉斯加隱沒帶的中深部震源特性並嘗試為中深部震源的破裂機制提供新的觀測證據。

(13)

不管是使用特定工具或是逆推法，都需要選擇一個參數優先計算，無法在一 次逆推中分離 f_c 與 t*。為數眾多的科學家們一直為減輕兩者間的制衡關係而努力。

近年來，不同的方法陸續推陳出新。如結合逆推法與EGF法的CS-BP法 [Yin et al., 2018] 與考慮震源與測站幾何的CEM法 [Ko et al., 2012]，都為解決 f_c - t* 制衡關係打開新的扉頁。

圖 1.1、fc-t* trade off 說明圖。 (a) fc對頻譜影響說明。fc越大對於頻譜衰減貢獻越 小，而較大的地震通常具有較小的 fc。 (b) t*對頻譜影響說明。t*越大對頻譜衰減 貢獻越大。 (c) fc-t*組合對頻譜的影響範例。利用 Brune 的震源模型，計算兩種不同的參數組合。在有噪訊的環境下，兩個參數組合的頻譜表現相似，對於震源分 析造成極大影響。[圖擷取自 Ko et al., 2012]

1.3 阿拉斯加隱沒帶概述

阿拉斯加隱沒帶位於太平洋火環帶北部，太平洋板塊向北隱沒至北美洲板塊之下。隱沒速率與板塊年齡均呈現由西往東遞減的趨勢：阿留申島弧地區每年隱

(14)

島嶼地區的 60 Ma 至大陸地區的 20 Ma [Goes et al., 2008]。隱沒的太平洋板塊在大 陸地區下方的隱沒厚度約為 45 ~ 55 km [Zhao et al., 1995]。阿拉斯加隱沒系統停止 於阿拉斯加灣 (Gulf of Alaska, 圖 1.2)，在深度 < 50 km 時，其隱沒板塊的傾角極 低，約為 3~5° [Ferris et al., 2003; Ohta et al., 2006]。平緩的隱沒傾角，造就較低溫 的地幔楔 (mantle wedge) 系統，也使得阿拉斯加大陸地區缺少火山弧。依據 Stachnik et al. (2004) 的研究，阿拉斯加的地幔楔較其他隱沒帶低約 150°C。

圖 1.2、阿拉斯加區域構造圖。太平洋板塊約以 55 mm/yr 的速率向西北隱沒至北美 洲板塊下方 [Ferris et al., 2003]。橘線區域為本論文的研究區域。黃線為板塊邊界，

三角形為海溝位置。黑線為阿拉斯加隱沒帶的等深線，間隔為 20 km。(USGS slab 1.0) [Hayes et al., 2012]

(15)

第二章研究方法

2.1 群組地震法 (Cluster-Event Method, CEM)

群組地震法由 Ko et al. (2012) 所提出，以減少逆推自由度的方式減緩 𝑓_𝑐 與 t* 之間的制衡關係 (圖 2.1)，其概念為：

(1) 將位置相鄰的地震視為同一個群組。

(2) 由於震波傳遞路徑相似，群組內地震傳至同一測站的路徑衰減 Q 視為一樣。

(3) 群組內一個地震只能有一個 f_c。

(4) CEM 目前忽略震源方向性 (source directivity) 造成的都卜勒效應。

CEM 建構多重的震源-測站組合能「有效嚴謹地」 (robust) 定義群組內共同路徑𝑡^∗與不同震源的𝑓_c。單一地震不限制只能在一個群組內，這類同時處在 2 個以上的震源稱為跨群組事件 (cross-over event)。群組組成與同一地震的破裂過程無關，

因此代表破裂行為的 fc，理論上不會因震源位在不同群組而改變，因此跨群組事件可約束 CEM。

圖 2.1、群組地震法相關描述。 (a) 群組地震法說明示意圖。鄰近震源到相同測站 共享 Q 值，群組內單一地震僅有一個 fc。 (b) 跨群組事件約束 CEM。紅色星星代表跨群組事件。若在 A、B 兩個個別計算的群組內獲得相似的結果，就可說明結果的嚴謹性。 (c) 群組地震法與其他方法在逆推自由度的比較。以群組內有 3 個地

(16)

CEM 立論基礎為波線 (ray) 因路徑相似而取樣 (sample) 相似的地幔，路徑衰減效應可視為相同。群組內多個震源同時約束路徑效應，從而確保能夠嚴謹的定 義 fc。在此限制下，若研究區域的構造有非常劇烈的變化，CEM 就需特別注意路徑的相似性。以隱沒帶為例 (圖 2.2)：如果震源位置較深 ( > 140 km)，地震波到達海溝後方、海洋岩石圈之上的測站時，其波線會大部分通過隱沒板塊。因為板塊構造的複雜性，鄰近的震源位置也無法保證相似的波傳路徑，所以此類測站接收到的波形須謹慎使用，或者乾脆剔除。規模大於 6 的地震由於破裂過程相對複雜，

亦不適用 CEM 之假設。

圖 2.2、鄰近震源的波傳路徑示意圖。對於群組內的個別地震來說，綠線及藍線的波傳路徑位在板塊內部，不一定是相似的波傳路徑。這種情況可能會不符合 CEM 的假設，逆推結果就有可能存在極大誤差。紅點為震源位置，三角形為測站，黑 線為隱沒板塊 (USGS slab 1.0 [Hayes et al., 2012])，十字為背景地震。

CEM 與傳統方法最大的差異在於，避免考慮先決定震源參數或是路徑衰減的問題，而是以最有效的方式，同時決定波傳介質的衰減效應與地震震源的效應。

而在殘差值的分布上，CEM 假設下的模型解落在等值線封閉區間內，保證了模型解的嚴謹性 (圖 2.3)。

(17)

圖 2.3、模型-資料殘差值分布比較。 (a) 單一震源-測站對殘差值圖。使用 Individual fitting (IF) 計算出來的結果，可觀察到 RMS 最小的區域呈現一個狹長形的帶狀，

並不能很有效的找到一個最佳解。╳為 CEM 計算出來的最佳解，並不會在 IF 計算出來最小值的狹長區域內，這是 CEM 綜合考慮所有震源-測站對的結果。(b) 群組震源-測站對殘差值圖。利用 CEM 搭配鄰域演算法計算出來的結果，僅改變其 中一組 f_c-t* pair 得到的結果。╳落在一個封閉的區間內，可較嚴謹定義模型。

2.2 fc、t*與地震波頻譜估計

地震發生後，在地表測站觀測到的波形，可以拆成四個部分：震源、路徑衰減、近地表場址效應與儀器響應 (圖 2.4)。地震儀記錄到的地震時間序列 A(t)，

可以表示為 (式 2-1)：

A(t) = 𝑆(𝑡) ∗ 𝐸(𝑡) ∗ 𝑅(𝑡) ∗ 𝐼(𝑡) (式 2-1) A(t) 為時間域振幅大小，S(t) 為震源效應，E(t) 為地震波的傳遞路徑效應，也可以想像為波傳介質對震波能量的衰減影響，R(t) 為測站附近的場址效應，I(t) 為儀器響應，符號 ∗ 代表捲積 (convolution)。

頻譜分析相較於時間序列分析在數學上有較簡單的處理方式，如式 2-1 在頻率域中的振幅 A(𝑓)，可寫作 (式 2-2)：

A(𝑓) = 𝑆(𝑓) × 𝐸(𝑓) × 𝑅(𝑓) × 𝐼(𝑓) (式 2-2) 捲積在頻率域中變成乘法，分離各個效應時有明顯的計算優勢。因此，本論文將

(18)

圖 2.4、地震頻譜組成示意圖。[修改自柯彥廷 (2010) 碩士論文]

本論文使用 Brune (1970) 提出的震源振幅頻譜 𝜔⁻² 衰減模型 (式 2-3)：

S(𝑓) = ^M⁰

1+(^𝑓

𝑓c)² (式 2-3)

S(𝑓)表示震源振幅頻譜，M0表示地震矩 (seismic moment)，f_c表示拐角頻率。f_c的物理意義在不同的震源模型有不同的解釋，但均為描述震源本身的破裂性質。在 Haskell 模型中 [Haskell, 1964]，假設斷層面是一個長方形，fc 代表的是上揚時間 (rise time) 與破裂時間 (rupture time) 的效應總合。在 Brune 模型 [Brune, 1970]，

斷層面假設為圓形，f_c與圓形斷層面的半徑有關; Beresnev (2002) 指出 Brune 模型 的 fc與斷層最大的破裂速度有關。

Q 值為對地層本質衰減 (intrinsic attenuation) 的描述。Q 值越大，表示地震波 在傳遞過程中，受到地層非彈性變形而消散的能量越小。E( f ) 表示波傳介質對地 震波能量的衰減貢獻 (式 2-4)，t*反映介質的衰減特性，是沿路徑積分的結果(式 2-5)：

E(𝑓) = e^−π𝑓t^∗ (式 2-4) t^∗ = ∫_source^{𝑟𝑒𝑐𝑒𝑖𝑣𝑒𝑟}_Q(r)^dτ (式 2-5) τ 為沿著波傳方向 (r) 以震波速度 c 傳遞距離 x 的時間。如果 Q 為常數，則 τ 可以積分為走時 T。從e^−π𝑓t^∗的描述中，可知在 t*相同的情況下，頻率 f 越大，能量衰 減越多。

(19)

結合式 2-3 與式 2-4，測站接收到的地震波頻譜可以表示為式 2-6 [Anderson and Hough, 1984]：

A(𝑓) =^𝐶𝑀⁰^𝑒^{−π𝑓𝑡∗}

1+(^𝑓

fc)² ，Ω₀ = 𝐶𝑀₀，𝐶 = ^𝐶^𝑈^𝜙𝜃

4𝜋𝜌𝑅𝑉𝑝3， (式 2-6) 其中 ^M⁰

1+(^𝑓

fc)² 描述震源效應，𝑒^{−π𝑓𝑡}^∗ 代表路徑衰減效應，而 C 包含與頻率無關的因 素，如幾何擴散、P 波速度等，Ω₀則是震波位移振幅頻譜的低頻部分 (圖 2.5)。式 2-6 中的Ω₀、𝑡^∗、f_c為本論文的模型參數。

圖 2.5、f_c、t*、Ω₀與頻譜的關係。斷層錯動之後，若地表有永久位移，位移時間序列會顯示一個躍升，躍升的幅度與滑移量 D 有關，此種波形稱為近場 (near-field)。

距離稍遠的測站，地表不會有永久位移，稱為遠場 (far-field)。遠場頻譜的低頻值 (Ω₀) 與地震矩有關，fc與 t*則共同決定遠場頻譜的形狀。[修改自 Shearer, P. M.

(2009). Introduction to seismology. ]

一般來說，頻譜分析多是使用快速傅立葉轉換 (Fast Fourier Transform, FFT) 的方式，將地震時間序列轉換成頻率域的訊號，轉換基底為正交的正弦函數與餘弦函數。頻譜計算由於有限帶寬與計算設備的限制，無法達到傅立葉轉換的無限高頻要求，在此情形下，頻譜會發生洩漏 (leakage) 的現象。頻譜洩漏導致計算出來的頻譜有劇烈的震盪，增加分析的困難。本論文利用多重視窗頻譜法 (Multi-taper spectrum estimation [Thomson, 1982; Prieto et al., 2007])，使用一系列的離散橢球正

(20)

頻譜洩漏的情形，使得分析的頻譜變得較為平滑 (圖 2.6)。本論文的分析資料為原始波形經由 Multi-taper 法計算而得出的振幅頻譜。(參閱附錄 B)

圖 2.6、地震波形、頻譜展示圖。 (a) 地震波形圖。藍色虛線為 P 波到時，灰底部分為 P 波到時前 1 s 至 P 波到時後 4 s，定為訊號的 5 s 時間窗。 (b) 訊號頻譜。

黑色線為 FFT 的計算結果，紅色線為 Multi-taper 法計算出來的結果。從結果來說，

紅線較黑線平滑。

2.3 鄰域演算法 (Neighborhood Algorithm, NA)

NA 的模型參數空間由 CEM 假設下的震源數與測站數之個數總和決定。每一 個參數 (f_c或是 t*) 都代表參數空間內的一個基底。以單一震源-測站對 (f_c- t* pair) 舉例，即展開一個 2 維的參數空間。這個 2 維平面上的任一點都可能是這個震源- 測站對的模型解。若使用格點搜尋法 (grid search)，土法煉鋼地計算每個模型點與資料之間的殘差值 (misfit)，不僅整體計算速度緩慢，還會消耗大量資源。

為了解決此種計算效率不彰的難題，Sambrige (1999a,b) 提出一種優化的演算 法，利用隨機均勻分布在參數空間內的少數模型點，先計算其與資料的殘差值，

再利用兩兩模型點中垂線切割參數空間，將每個切割出來的區域 (voronoi cell) 內的殘差值由已經取樣的模型點代表 (圖 2.7)。這樣就可以大略知道整個參數空間內的殘差值變化。藉由分辨殘差值的大小，在殘差值小的區域給予更多的採樣點，

以減少無效率的計算。此種先求得參數空間大略殘差值分布 (misfit surface)，再給予不同數量採樣點的作法，便稱為鄰域演算法 (參閱附錄 A)。格點搜尋法依然可以得到相似的結果，只是計算效率較差。

(21)

圖 2.7、鄰域演算法示意圖。 (a)藍色點及紅色點為已採樣的模型點，參數空間被切割為不同大小的區域 (voronoi cell)。灰色區域內的所有模型殘差值由紅點代表。

藉此勾勒出參數空間大略的殘差分布。 (b) 圖(a)中的灰色區域生成 7 個新模型，

並重新切割整個參數空間。由此找出最符合資料的模型解。 [修改自 Sambridge, 1999a]

NA 做為一個有效率的全域搜尋演算法 (圖 2.8)，僅有 2 個控制參數： nr (每次迭代較小殘差值的 nr 個區域) 以及 ns (每次迭代總共要新取樣 ns 個點)。每次迭代生成新取樣點時，將新取樣點分布在殘差值較小的 nr 個區域，避免最佳解被局限在區域最小值 (local minimum) 內。藉由 nr 與 ns 的互相搭配，即可有效率的在廣大的參數空間內找出最符合資料的解。

本論文中，使用 nr = 100、ns = 500 與迭代次數 50 次。NA 的操作步驟如下：

(1) 在參數空間內隨機給予 1000 個初始模型，將參數空間切割成 1000 個區域。

計算初始模型與資料之間的殘差值。

(2) 在殘差值最小的 100 個區域內，平均分配 500 個新模型，即每個區域需要隨機取樣 50 個新模型點。

(3) 將整個參數空間內重新切割成 1000+500×(目前迭代次數)個區域。

(4) 回到步驟 (2)，直到結果收斂或是迭代次數結束。(圖 2.9)

(22)

圖 2.8、NA 迭代與殘差值的關係。此群組共有 4 個測站與 4 個震源，參數空間維度為 8。在迭代 20 次後，殘差值基本上呈現穩定狀態，顯示 NA 的效率。為了確保不同維度也可以達到收斂的情形，本論文將迭代次數定為 50 次。紅線為全區殘差最小值，淺藍線則為每次迭代中模型的殘差最小值，若比目前全區最小值小，

則會取代成為全區最小值。綠線為每次迭代新產生模型的殘差值平均，若趨於一定值，代表 ns 個點均生成在相似的位置，即收斂之意。

圖 2.9、NA 的收斂情形。第一次迭代時採用隨機均勻採樣。第二次迭代後，進行 NA 的計算流程，隨著迭代次數增加，模型解逐漸收斂在一個區域內。

(23)

第三章資料處理與分析

3.1 資料來源

本論文資料來自阿拉斯加地區數個地震觀測網，可在 IRIS (Incorporated Research Institutions for Seismology) 資料中心下載取得 (http://ds.iris.edu/wilber3/find_event)。主要的觀測網如下：Alaska Regional Network, (IRIS 代碼: AK)、Alaska Volcano Observatory (AV)、USArray Transportable Array (TA) 與 National Tsunami Warning Center Alaska Seismic Network (AT) (圖 3.1)。

本論文中的測站取樣率 (sampling rate) 有 2 種：40 Hz 與 50 Hz。Nyquist frequency 為 20 Hz 與 25 Hz。本論文統一設定模型頻譜擬合的最高頻為 20 Hz 的 0.8 倍 [Ko et al., 2012]，即 16 Hz。

圖 3.1、本研究區域 (a) 地震及 (b) 測站分布圖。黃點代表震央位置，方塊代表測站位置，其他圖示同圖 1.2。

(24)

3.2 資料處理流程

本論文對地震波的振幅頻譜進行分析。其步驟可統整為：

(1) 將地震時間序列 (垂直分量)，去除儀器響應後得到位移，並尋找 P 波到時。

(2) 定義訊號、雜訊時間窗；檢視是否通過篩選條件。

(3) 通過篩選條件的波形，利用 Multi-taper 法計算頻譜，並進行人工篩選。

(4) 組成群組。

(5) 藉由 CEM-NA，尋找出最符合資料的模型解。

接下來將針對步驟 (1) 至步驟 (3) 說明 (圖 3.2)。

3.2.1 資料前處理-P 波到時

本論文利用 IRIS 提供的 wilber 3 平台，進行初步地震搜尋。以 IRIS 發布的地震時間為基準，下載發震時間前 2 分鐘至後 5 分鐘，總計 7 分鐘的地震資料。本論文地震資料的經緯度介於56˚N ~ 66˚N 與 142˚W ~160˚W，選用時間為 2012 年 1 月至 2017 年 9 月，地震規模介於 3 ~ 6，震源深度不限 (圖 3.1)。

首先將原始波形垂直分量 (Z 分量) 的平均振幅值 (mean) 與線性趨勢 (linear trend) 移除，並且將波形兩端各 5%的部分給予尖滅函數 (taper)。原始資料去掉儀器響應後變成位移，並加入抗混疊濾波器 (Anti-aliasing fliter, SAC code : freq 0.05 0.1 18 20)。位移經微分成速度，帶通濾波 2 ~ 8 Hz。本論文利用 Akaike information criterion 值 (簡稱為 AIC, 式 3-1, [St-Onge, 2011]) 自動挑出在速度時間序列的 P 波到時，並輔以人工檢視確認。最後將 P 波到時記錄至位移資料中。

AIC(k) = k × log(var(y(1: k))) + (n − k − 1) × (log(var(y(k + 1: n))) (式 3-1) k 會在 1 (起始點)與 n (總點數)之間，log 以 10 為基底，var(y(1: k))表示第一點到第 k 點振幅值的變異數。取 AIC 最大值前的最小值即為 P 波到時 (圖 3.3)。

(25)

圖 3.2、資料前處理與資料篩選流程圖。

下載原始資料發震時間-2~+5 mins

去儀器響應到位移

(詳見 3.2.1 節)

帶通濾波(詳見 3.2.2 節)

(bp co f1 18 n 3 p 2) f1 = 0.1 (Initial)

SNR ≧ 5 f1 ≦ 1 Hz 自動挑選 P 波到時

(詳見 3.2.1 節)

PCR ≧ 1 f1 ≦ 1 Hz 目視確認

P 波到時

SNR < 5 捨棄

f1 > 1 f1 > 1

通過資料篩選可使用資料

SCR < 1

f1 = f1 + 0.1

目視確認頻譜計算的品質

(26)

圖 3.3、利用 AIC 自動挑 P 波到時的說明。(a) 應用 AIC 法之範例。藍線為波形，

紅線是 AIC 值，黑線則為 P 波到時。紅藍線均將最大值定為 1。 (b) P 波到時前後各 2 s 的波形與 AIC 值，綠線為 P 波到時。

3.2.2 資料篩選標準- SNR、PCR

為了盡量降低構造對頻譜的影響，獨立顯著的 P 波是本研究的重要條件。因此，本論文在時間域中設定兩個條件：SNR (Signal Noise Ratio) 與 PCR (P-Coda Ratio)。原始資料經過 3.2.1 節挑出 P 波到時後，將位移時間序列進行 f₁ ~ 18 Hz 帶 通濾波，初始 f1值定為 0.1 Hz。將 P 波到時前 1 s (以下簡稱 P-1 s) 到 P + 4 s 定為 5 s 的訊號時間窗，雜訊時間窗定在 P - 7 到 P - 2 s (圖 3.4a)。

圖 3.4、波形展示圖。 (a)SNR 訊號及噪訊位置圖。灰底部分為訊號，取 P - 1~P + 4 的 5 s 時間窗。紅線圍成的區域為背景雜訊，取 P - 7~P - 2 的 5 s 時間窗。 (b) P 波及尾波位置圖。取 P ~ P + 2 的時間窗為 P 波，P + 2 ~ P + 4 為尾波。

本論文計算兩個時間窗的振幅方均根比值 (RMS) (式 3-2 與式 3-3)，定為SNR (式 3-4)。

RMS_signal = √^∑ ^(y^i,signal⁾²

n

i=1 ，P − 1~P + 4 s (式 3-2)

P P

(27)

RMS_noise= √^∑ ^(y^i,noise⁾²

n i=1

n ，P − 7~P − 2 s (式 3-3) SNR =^RMS^signal

RMS_noise (式 3-4)

yi代表位移振幅值。若SNR > 5，即進入PCR計算；若SNR < 5，則將 f1值加 0.1 ( f₁= 0.1 → f₁= 0.2 )，再進行SNR計算，直到SNR ≧ 5 (圖 3.5)。當 f₁> 1 Hz 時，

若 SNR 仍 < 5 ，則該資料捨棄。

圖 3.5、不同濾波頻段對 SNR 的影響。 (a)位移波形。 (b)(c)(d) 為 f1遞增的濾波 效果對 SNR 的影響。bp co f1 23 n 3 p 2 為使用 3 階的 butterworth 濾波器，進行零 相位帶通濾波 f₁ ~ 23 Hz。SNR > 5 的資料通過篩選。

以上步驟是在噪聲背景中找到顯著 P 波，但並未考慮 P 波後方的尾波 (coda)。

訊號與尾波的振幅比會直接影響 P 波頻譜，因此本論文依照SNR的做法，定義 P ~ P + 2 s 為 P 波時間窗；P + 2 ~ P + 4 s 為 Coda 時間窗 (圖 3.4b)，計算兩個時間窗的振幅方均根比值 (式 3-5 與式 3-6)，定為PCR (式 3-7)。

RMS_p = √^∑ ^(y^i,p⁾²

ni=1

n ，P~P + 2 s (式 3-5)

(28)

PCR = ^RMS^p

RMS_coda (式 3-7)

若PCR ≧ 1，則通過PCR標準。若無，則回到SNR計算，將 f1 + 0.1 直到同時通過SNR與PCR標準 (圖 3.6)。如此篩選過的波形才會進行下一步的 Multi-taper 頻譜計算。最後頻譜經過人工檢視 (圖 3.7)，才會成為可使用的資料。

圖 3.6、利用 SNR 與 PCR 篩選波形的例子。(a)(b) 為同一震源，不同測站的記錄。

(a)與(b)同時通過 SNR 標準(7.10 與 5.17)但(b)未通過 PCR 標準 (0.55)，因此(b)波形被淘汰。

圖 3.7、頻譜目視篩選例子。(a) 好的頻譜例子，頻譜大致呈現 𝜔⁻²衰減。 (b) 地震波可能在板塊內部反彈多次導致頻譜高頻不呈現 𝜔⁻² 衰減。

3.3 地震群組組成

對每個目標地震來說，群組半徑越小，越能保證路徑的相似性，但相對的，

(29)

在群組內的地震就會越少。對於 CEM 來說，群組半徑與群組內震源密度形成制衡 關係 (trade-off)。本節說明 CEM 的群組組成與篩選條件。

通過篩選的資料，才會組成群組。本論文不限定一個震源僅能存在一個群組內。首先選定某個震源當做群組的中心地震，計算震源與測站之間的最小直線距離 R，以中心地震為球心，取 r = 0.05 × R，且 r ≦ 30 km，形成一個以 r 為半徑的球體，球體內的震源視為一個群組 (圖 3.8)。

圖 3.8、距離篩選示意圖。左圖描述群組震源的形成，球體以 r 為半徑，球體內的震源視為一個群組。右圖描述震央與測站的距離篩選，圓形以 2R 為半徑，圓形外的測站不使用[以(X)表示]。星星為震源或震央，三角形為測站。

本論文希望地震波傳路徑不要太長以免帶入過多的構造資訊，因此以中心震央為圓心，2R 當作半徑，2R 最大值為 550 km (震央距 5˚)，若測站與中心震央的距離 > 2R，則該測站不被使用。最終若一個群組內震源數或測站數不足 3 的話，

則該群組不被接受。

假設已經有一組 3 個地震的基本群組，我們會進一步檢視這 3 個震源是否同時被 3 個測站接收到，若無，則該群組捨棄；若有，則再加入 1 個震源，重複檢視是否被同時 3 個以上測站接收到 (圖 3.9)。換句話說，地震與測站的數目不能差異太多。這樣做的好處是群組內地震與測站處於「填滿」的狀態：N 個震源完全擁有相同的測站，而每個測站也擁有相同的震源，這個步驟增加逆推的嚴謹性。

但是缺點是在群組的組成上面不夠靈活，容易侷限在幾個震源密集出現的地區。

(30)

圖 3.9、本論文群組組成說明。左圖為群組示意圖，震源 1 為中心地震，其他震源依距離中心地震之遠近給予編號。右邊標示 Round 1 表示第一階段組成群組：以 3 個震源為基準，測試是否有 3 個以上相同的測站。126，表示 126 這 3 個震源組合沒有相同的 3 個測站。Round 2 將通過 Round 1 的震源組合加入一個震源做檢視，

若有相同的 3 個測站以上，就會進到 Round 3；若無，則接受 Round 1 的群組組成結果。舉例來說，圖中的 12356 表示沒有共同的 3 個以上測站，故接受 1235 的群組組合。右邊標示僅列出部分組合做為展示之用。紅色圈代表群組成立。

3.4 CEM-NA 的殘差值計算

在進行 CEM-NA 逆推時，頻譜擬合的殘差值 (misfit) 是在對數軸上做計算。

若依照一般的點數分布計算，在 10 Hz 後的點密度會遠大於 10 Hz 之前，殘差值將主導於高頻的點數。為了避免此現象，本論文頻譜在對數軸上重新採樣 (圖 3.10)，

確保低頻與高頻對於殘差值的貢獻是一致的。

圖 3.10、對數軸均勻採樣說明，舉出 2 個測站頻譜做為範例。藍點為 Multi-taper 計算出來的頻譜，紅點為均勻採樣之後的頻譜。為了看出差異，藍點振福值為紅點的 10 倍。

(31)

CEM-NA 逆推過程中，本論文將模型擬合最高頻設定為 16 Hz，故將 fc的 NA 搜尋範圍定在 0.1-12 Hz 之間，即限制 f_c的擬合時不能超過 12 Hz； t* 的搜尋範 圍定在 0.001-0.4 s [Ko et al., 2012]。

假設有 3 個震源形成群組並有 4 個測站，其模型參數空間維度即為 7。即 NA 產生每一個模型點 X 裡面會有 7 個參數，X = [𝑓_𝑐,1 𝑓_𝑐,2 𝑓_𝑐,3 𝑡₁^∗ 𝑡₂^∗ 𝑡₃^∗𝑡₄^∗]。Fkj代表第 k 個震源-第 j 個測站模型與資料之間的殘差值。Fkj是 fc、t*與 Ω0的函數 (式 3-8)，

3 個參數會在 CEM-NA 的過程中一起決定：

F(𝑓_𝑐, t^∗, Ω₀)_kj = min (

∑ ω_𝑖×‖^Ω0^𝑒^{−π𝑓𝑖𝑡𝑗}

∗

1+(𝑓𝑖 fc,𝑘)2

−d_kj,𝑖‖ 2 Mi=1

∑^M_i=1ω_𝑖 ) (式 3-8) ω_𝑖表示第 i 個點的權重，由訊號與噪訊的振幅頻譜比值 (δ_i) 決定 (式 3-9 與圖 3-11)，

d_kj,𝑖則為振幅頻譜第 i 個點的值，‖ ‖₂代表 L₂-norm。

δ_i = |log₁₀(f_i,signal) − log₁₀(f_i,noise)|

ω_𝑖 = 0 δ_i< 0.5 ≈ log₁₀(3) ω_𝑖 = 2δ_i 0.5 ≤ δ_i ≤ 2

ω_𝑖 = 4 δ_i > 2 (式 3-9)

圖 3.11、權重說明。本論文使用訊號頻譜與雜訊頻譜在對數軸上的差值δ_i，做為權重基準。

(32)

Ω₀與 fc - t*之間的制衡關係無關，因此本論文在逆推過程中，不將Ω₀放入 NA 的模型參數空間中。本論文計算資料頻譜頻率小於 1 Hz 的振幅平均 (設為 λ)，在 0.05λ 到 50λ 的區間內，以 0.05λ 為間隔，利用式 3-8 計算每一個 λ 倍數值與同一 組 f_c,k- t*_j模型的頻譜殘差值，取最小殘差的λ 倍數值做為該組 fc,k - t*_j的Ω₀ (圖 3.12)。

該組 fc,k - t*j-Ω₀模型所計算出來的殘差值即為 Fkj，Θ 為所有 Fkj的總和平均，也代表模型點 X 的整體殘差值 (式 3-10)。經由 NA 的計算流程後，最小的 Θ’代表的模型X’就是最符合資料的模型。

Θ =

^∑ ^∑ ^F^𝑘𝑗

𝑛𝑗=1 𝑑𝑘=1

𝑛×𝑑 (式 3-10)

n 為測站數，d 為震源數。

圖 3.12、Ω₀定義說明，同一個震源-測站對 (λ 搜尋範圍相同)，兩個不同的 X 模型 點。 (a) f_c- t*組合： f_c = 5.71，t* = 0.062。(左圖) 黑線為資料頻譜，並計算同一 組 f_c- t*模型下不同 λ 倍數值的殘差值，並繪成右圖。紅線為 f_c- t*- Ω₀模型形成的 頻譜，灰實線為 f_c- t*-λ 模型形成的頻譜，灰虛線為 f_c- t*-0.1λ 模型形成的頻譜，

藍虛線為Ω0的搜尋區間界線。(右圖) 淺藍線為殘差值，紅點為最小殘差值 (F_kj)，

其所代表的λ 倍數值為 Ω0。(b)另一組 f_c- t*組合： f_c = 2.33， t* = 0.31 的 Ω₀搜尋過程。

(33)

第四章逆推結果與震源參數

本章節共分成三個部分，第一部分為振幅頻譜的擬合結果與 CEM-NA 逆推出 來的 fc與 t*範例。第二部分為利用 fc計算出僅考慮斷層錯動前後物理量變化的靜態震源參數。第三部分為考慮斷層錯動的過程中，物理量隨能量變化的動態震源參數。

4.1 CEM-NA 之逆推結果 4.1.1 資料擬合

本論文利用擬合模型的變異數縮減 (Variance Reduction, VR) 來判斷資料擬合的好壞 [(式 4-1) ~ (式 4-5)]。整體來說，模型擬合的越好，VR 值會越高。本論文總計 227 個群組，所有的資料 VR 均 > 70 % ； 99 % 的群組 VR > 80% ；有 68%

的群組 VR > 90% (圖 4.1)。

VR 計算，如果資料點權重 ω_𝑖 = 0 (式 3-9) 則不納入平均值與模型值計算，

令 xi為第 i 個觀測資料，yi為第 i 個模型點，i =1~n，n 為總點數。

𝑋_𝑖 = 𝑙𝑜𝑔₁₀(𝑥_𝑖) (式 4-1) 𝑌_𝑖 = 𝑙𝑜𝑔₁₀(𝑦_𝑖) (式 4-2) 資料變異數定義為

𝑉_{𝑑𝑎𝑡𝑎} = ∑^𝑛_𝑖=1(𝑋_𝑖 − X̅)² (式 4-3) X̅表示平均值，擬合模型變異數則定義為

𝑉_{𝑚𝑜𝑑𝑒𝑙} = ∑^𝑛_𝑖=1(𝑋_𝑖 − 𝑌_𝑖)² (式 4-4) 最後計算 Variance Reduction：

VR = [1 − (^𝑉^{𝑚𝑜𝑑𝑒𝑙}

𝑉_data )] × 100% (式 4-5)

(34)

圖 4.1、Variance Reduction 之統計圖。大部分 (99%) 的資料 VR 值在 80% 以上。

4.1.2 CEM-NA 的逆推範例地震

本節以群組組成測站數 4 與震源數 4 內的其中一個地震例子來呈現 CEM-NA 如何取得最佳解 (圖 4.2、圖 4.3 與圖 4.4)。

圖 4.2、CEM-NA 之逆推範例，範例群組的震央位置 (白點)，紅點為波形展示的地震位置，測站位置 (三角形)。範例地震規模 4，發震時間：2017.066.09.50.08，深度 109.2 km。其餘圖示與圖 1.2 相同。

(35)

圖 4.3、CEM-NA 之逆推範例，頻譜擬合。(左)波形。灰底部分為訊號時間窗，紅線圍成的區域為背景雜訊時間窗 (右)頻譜擬合之結果。黑線為訊號振幅頻譜，紅線為模型，灰虛線：雜訊頻譜。紅虛線：ω_𝑖 = 0的地方，不予列入計算。左下角 的數字如 (5.71 , 0.01) 代表 (f_c, t*)。16 Hz 為本論文模型擬合的最高頻。A, B, C, D，

分別對應圖 4.4 中的 A, B, C, D 全域殘差值分布。

(36)

圖 4.4、CEM-NA 之逆推範例，全域的殘差值分布。由最佳解 X 中僅改變其中一組 (fc, t*) 的值，並計算全部格點的殘差值。可以注意到模型解落在封閉區間形成的谷底，代表模型解的嚴謹性。A, B, C, D，分別對應圖 4.3 中的 A, B, C, D 頻譜擬合。

4.1.3 跨群組事件的 f_c

跨群組事件約束 CEM-NA 的結果。跨群組事件為同時屬於 2 個以上不同群組的震源 (圖 4.5)。如果 CEM-NA 逆推結果是嚴謹的話，理論上跨群組事件的模型參數應該要相似。我們利用兩個方式做歸納：一個是計算個別地震在不同群組的平均值與標準差，另一個計算整體跨群組事件的變異數。

(37)

在個別地震方面，82.4%的 fc 標準差< 2 Hz；有 58.33%的跨群組震源 fc 標準 差 < 1.5 Hz。而整體 f_c的變異數為 2.285 Hz。

整體變異數計算方式為式 4-6：令𝑓_𝑐,𝑖為單一震源位在第 i 個群組內的 fc，𝑓̅_𝑐為 該震源的 f_c平均，則整體變異數可寫成：

Variance = ^∑ ^∑ ^(𝑓^𝑐,𝑖^−𝑓^{̅ )}^𝑐

2 𝑛𝑖=1

𝑘𝑗=1

∑^𝑘_𝑗=1𝑛_𝑗 (式 4-6) 其中 k 代表跨群組的震源個數，本論文中有 108 個，n 代表該震源總共跨了多少個群組，∑^𝑘_𝑗=1𝑛_𝑗表示所有震源所跨的群組總個數。

圖 4.5、跨群組事件的 f_c。上圖為該震源的跨群組個數。下圖為 108 個震源的平均值與標準差。

4.2 靜態震源參數計算

本論文利用 fc與 Ω0求取靜態的震源參數。「靜態」是僅注意現象始末，不注意其發生過程的物理量，意即與動力學 (dynamic) 無關的物理參數。本節說明地震矩與靜態應力降。

4.2.1 地震矩 (Seismic moment, M0)

地震矩為指示地震釋放能量大小的物理量。可由剪切常數 μ、破裂面積 A 與

滑移量 D 求得 (式 4-7) [Kanamori, 1977]：

(38)

亦可由震矩張量解 (moment tensor solution) (式 4-8)求得：

𝑀₀ = ¹

√2(∑ 𝑀_𝑖𝑗 _𝑖𝑗²)¹² (式 4-8) 本論文利用 Brune-type 的公式 (式 2-6) 與 P 波頻譜擬合出來的 Ω₀，計算每個測站的結果 (式 4-9)，再將之平均作為此震源的地震矩 (式 4-10)：

𝑀_0,𝑖 = ^4𝜋𝜌^𝑖^𝑅^𝑖^𝛼³

𝑈𝜙𝜃 𝐶

Ω_0,𝑖 (式 4-9)

𝑀₀ = ¹

𝑛∑^𝑛_𝑖=1𝑀_0,𝑖 (式 4-10) 𝑈_𝜙𝜃

𝐶 是 P 波 radiation pattern 的球狀平均為 0.52 。 [Aki and Richard, 1980;

Abercrombie, 1995]， ρ 為震源附近的密度 (kg/m³)， R 為震源與測站之間的距離 (m)，α 是震源附近的 P 波速度 (m/s)。M0的單位: (牛頓-米, N-m)。震源附近的物 理參數使用阿拉斯加地區一維物理參數模型 [Jiyao et al., 2013] (圖 4.6)。

圖 4.6、阿拉斯加地區一維物理參數圖。 (a) Vp、Vs 隨深度的變化。 (b) 密度隨 深度的變化。虛線之間為本論文震源深度的範圍。 [資料取自 Jiyao et al., 2013;

YoungHee. et al., 2014]

(39)

另一方面，本論文利用阿拉斯加近震規模 (ML) 與震矩規模 (Mw) 的關係轉換式 (式 4-11 與式 4-12)，將 IRIS 發布的近震規模轉換為震矩規模 (圖 4.7)：

Alaska, depth < 40 km, Mw = ML (式 4-11) Alaska, depth > 40 km, M_w = M_L+0.03 (式 4-12)

圖 4.7、阿拉斯加地區近震規模 (ml) 與震矩規模 (mw) 轉換圖。[圖擷取自 Ruppert and Hansen, 2010]

並利用震矩規模轉換為地震矩 (式 4-13) [Hanks and Kanamori, 1979]：

M_w =²

3log₁₀(M₀) − 10.73 (式 4-13) M0的單位為達因-公分 (dyne-cm)。

藉由式 4-10 平均計算出來的 M₀比地震目錄發布的 M₀為大 (圖 4.8)，可能是式 4-11 與式 4-12 關係轉換式的誤差造成的。接下來的震源參數計算，我們將使用 P 波頻譜擬合出來的平均 M₀。

(40)

圖 4.8、阿拉斯加地區地震目錄 M0 (catalog M0) 與頻譜擬合 M0 (average M0) 的關係 圖。橫軸為地震目錄 M0，利用式 4-11 與式 4-12 轉換而來，縱軸為式 4-10 計算出 來的平均 M₀。

4.2.2 靜態應力降 (Stress drop, Δσ)

靜態應力降指的是在斷層錯動前後，斷層面上的剪切應力 (shear stress) 變化 值，以下簡稱為應力降。由 CEM 定出的 fc，我們可以計算 Brune-type 圓形斷層面的半徑 rc (式 4-14) [Madariaga, 1976]：

r_c= k^𝛽

f_c (式 4-14)

k為常數，破裂速度Vr = 0.9β時，k = 0.32。β是震源附近的S波波速 (圖4.6)。則應 力降可由地震矩與震源半徑求得 (式4-15) [Eshelby, 1957]：

∆σ = ⁷

16 M0

rc3 (式 4-15)

將式 4-14 帶入式 4-15 得到式 4-16：

f_c = 0.42𝛽[^∆σ

M0]^1/g (式 4-16) M0單位為牛頓-米(N-m)，Δσ 單位為 MPa，β 單位為 m/s。本論文使用逆推出來的 f_c與平均 M₀在對數軸上做線性迴歸，得到式 4-16 的 g 值為 4 (圖 4.9)。

(41)

圖 4.9、阿拉斯加地區平均 M₀與 f_c的關係圖。圖中的虛線都是利用式 4-16 並假設 g = 3 計算不同應力降的參考值，不同顏色為不同速度模型的計算結果 (藍虛線為 Vs，黃虛線為 Vp)。藍色實線為本論文計算出來的迴歸直線斜率為 1/4 (g = 4)。

f_c與 M₀在前人研究統計上呈現負相關的趨勢，g 值則有 g = 3 [Lio et al., 1986]

與 g = 4 [Hiramatsu et al., 2002] 2 種情況。不過 Hiramatsu et al. (2002) 指出，g = 4 的情形很可能是因為有限帶寬造成的假像。如果震源符合自我相似性規範，代表不同地震矩的地震有著相似的破裂行為，在規模較小 (ML < 1) 的地震 g = 3 依然會成立。本論文在後續的震源參數計算中，阿拉斯加的震源資料表現出自我相似性的特徵 (圖 4.11)，因此 g = 3 將被用在本論文的應力降計算。阿拉斯加地區的應力降中位數為 17.66 MPa，平均值與標準誤差為 18.34 ± 1.1 MPa (圖 4.10)。

(42)

4.3 動態震源參數計算

本節將計算動態的震源參數，考慮物理量在錯動過程中的變化。

4.3.1 地震體波輻射能 (Radiated Energy, E_R)

輻射能在本論文指的是地震體波 (P 與 S 波) 傳遞出來的總能量。輻射能可藉 由積分 P 波頻譜並回推 S 波能量的方式求得 (式 4-17) [Haskell, 1964; Boatwright, 1980; Boatwright and Fletcher, 1984; Venkataraman and Kanamori, 2004a, 2004b]：

E_R = ^8π

15ρα⁵(1 + q) ∫ 𝑓₀^∞ ²|Ṁ(f)|²d𝑓 (式 4-17) f 為頻率，ρ 為震源附近的密度 (kg/m³)，α 是 P 波速度 (m/s)，q 為 S 波與 P 波的 能量比值，如果兩者頻譜的 fc相同 (f_c^𝑃 = f_c^𝑠)，可寫成式 4-18，Ṁ為 Brune 模型的震源時間函數 (式 4-19)。

q =^E^R^𝑆

E_R^𝑃= 1.5(^α

β)²，if f_c^𝑃 = f_c^𝑠 (式 4-18) Ṁ(f) = ^M⁰

1+(^f

fc)² (式 4-19)

將式 4-18 與式 4-19 帶入式 4-17 得到式 4-20：

E_R= ^8π

15ρα⁵[1 + 1.5 (^α

β)²] ∫_0.001¹⁰⁰ 𝑓²| ^M⁰

1+(^f fc)²|

2

d𝑓 (式 4-20)

本論文輻射能計算頻率範圍 f 為 0.001~100 Hz。

輻射能大小會受到地震釋放能量大小的影響，直接比較各個震源的輻射能並不公平，所以本論文將輻射能與地震矩相除，引入 ER/M0這個參數 (scaled energy，

比例能量)，相當於計算傳遞出來的地震體波能量占整體能量的多少。阿拉斯加 ER/M0的平均值與標準誤差為 (3.05 ± 0.28) × 10^-5。本論文將比例能量與平均地震 矩作圖，並觀察不同 M0的比例能量變化 (圖 4.11)。我們發現比例能量基本上不隨平均地震矩變化，其暗示不同地震矩的地震有相似的破裂行為，呈現地震的自我相似性。

(43)

圖 4.11、E_R/M₀與 M₀的關係圖。M₀= 10^14.5與 10¹⁵ 將 M₀區分為 3 個區域。這 3 個區域的平均值與標準誤差分別為 2.44 ± 0.19、2.52 ± 0.14 與 4.99 ± 0.36，均需乘上 10^-5。E_R/M₀稍微與 M₀的變化有關，但整體來說還是平均值保持一致，其顯示地震理論的自我相似特徵。藍∆：區域平均值。