P波震源加速度時間函數之研究－以集集地震為例

全文

(1)國立台灣師範大學地球科學系碩士班. 碩士論文. P 波震源加速度時間函數之研究－以集集地震為例. 研究生：高毓琦指導教授：陳光榮博士中華民國 107 年 2 月.

(2) 摘要本研究選取強地動觀測計畫(Taiwan Strong Motion Instrumental Program, TSMIP)的集集地震事件序列之垂直分量、南北分量、東西分量等三分量合成之總量記錄，對於收錄到集集主震的四個分區：台中分區(TCU)、花蓮分區(HWA)、嘉南分區(CHY)與台東分區(TTN)進行分析。使用複數解載(CMD)與三維 Qp 值模型解算出 4~10 Hz 內各測站之震源加速度函數；先做同一分區的場址效應平均化，再做各個分區間比較。以同樣的方式再做集集地震序列之一餘震各分區之平均震源加速度函數，同樣進行各分區間比較，並進行主震與餘震事件的比較。研究結果有以下發現：集集主震與餘震的震源加速度函數，去除異常紀錄之台中測站，各分區的震源加速度時間函數，都有共同的狀態。首先從 0 至 4 秒之間僅有微量的水平振幅，約在 0.5~1 gal 之間，4 秒以後開始有所增幅，至第 7 秒多達到峰值，然後振幅逐漸減低，之後因不同分區而有不同的減低幅度，直至原來前 3~4 秒時的振幅。50 個測站中有 16 個測站其 S-P 波走時差在 10.24 秒以內的花蓮分區、與 35 個測站中有 5 個測站其 S-P 波走時差在 10.24 秒以內的嘉南分區，震源加速度時間函數最大振幅皆達到 3 gal。全部測站 S-P 走時差皆在 10.24 秒以上之台東分區，震源加速度時間函數最大約 1 gal。 I.

(3) 集集餘震的震源加速度函數，因所有分區走時差在 6 秒以上，且起始與結束部分波動幅度差不多大，這代表了所有分區的震源加速度函數只有 P 波到達，或是有部分 S 波滲入，並未像集集主震台中分區之震源加速度時間函數，受到表面波滲入影響，在時間末端的加速度時間函數有極大的振幅。將集集餘震所有分區的震源加速度函數進行比較，會發現此次餘震事件中，台中分區與花蓮分區的測站走時差皆在 10.24 秒內，花蓮分區走時差中位數約在 7.711 秒、台中分區走時差中位數約 7.716 秒，也表示震央距中位數兩者相差不大，但台中分區的震源加速度函數最大振幅卻仍達到±15 gal，花蓮分區的震源加速度函數最大振幅卻僅有±5 gal，由於分區裏所有測站之 S-P 波走時差皆在 6 秒以上，可以排除是表面波到達所影響，造成兩個分區震源加速度函數差異甚大原因，推測台中分區與同一事件其他分區相比，場址效應明顯大於同一事件其他分區。而餘震事件中摒除較大場址效應之台中分區不計，其震源加速度時間函數之變化，在 0~5 秒以內，花蓮分區有 3~4 gal、台東分區有 2~3 gal、嘉南分區有 0.5~1.5 gal 的起伏，5 秒過後開始出現反向振幅，在 7 秒時達到峰值，直到 10.24 秒結束為止；花蓮分區最大振幅可達到 5gal、台東分區最大振幅可達到 3gal、嘉南分區最大振幅僅有 2gal；由以上結果發現，震央距越遠，震源加速度函數之最大振幅越小。 II.

(4) 關鍵字：震源函數、複數解載、路徑效應、Q 值模型、集集地震、自回歸－赤池訊息準則 Keyword: source function, complex demodulation, path effect, Q model, Chi-Chi earthquake, Autoregressive - Akaike Information criterion. III.

(5) 致謝本篇論文的誕生，一路上經過了許多的波折，遇到許多的人事來來去去，有消失不見的人事物，也有新出現的人事物。經歷過這些大風大浪，我才知道學校這個環境的珍貴，能夠讀書真的是何其有幸，可以和實驗室裡的成員們共享資訊、推心置腹而不會有所顧忌，是在任何一個職場上永遠都找不到的價值。撰寫這篇論文，花費了我四年的碩班就讀時間，與職場工作的三年時間。在此我最想感謝的是我的指導教授陳光榮，他雖然有著豐富的地球物理學識，卻從來不會強迫學生們吸收所有的知識，等著有興趣的學生自由探索與追尋，也因此我能夠自由自在地想像、撰寫程式並付諸執行計算，實驗室有一半的程式是我整合前人之心血結晶，研發出較簡化使用步驟的新程式。感謝口委們對於我論文所提供的建議，為此也讓我來回奔波多次，才有這樣的成果。感謝就學期間建勳、仁聖與貞伶，雖然只陪我前兩年，但有了他們的論文當基礎，才有這篇論文的產生；還有外系研究助理昱達，謝謝你維持實驗室內計劃案的執行；感謝一路相挺的正宇學長、丞皞學長、中途退出的興邦學長，三位學長在各自之學校皆有任職，卻還是抽空回來攻讀博班，並擔負起指導後進的責任；中途加入的立程學弟、若涵學妹，你們在其他領域的表現應用在研究上，真的讓我覺得後生可畏；從大學時代到碩士 IV.

(6) 班，一路看著我就讀師大地科的師母，謝謝您一路上對我的照顧，還在我在氣象局發生事故時出手相救；另外還有好大、花紋很特別的 Dino，每次總是靜靜地在一旁聽著我們跟老師的討論，甚至四處走動討東西吃，但是只要摸到你的毛，心情就會放鬆下來，也為了要來看你，增加了我跟指導教授討論的次數。貴賓演講時間所接待的呂副主任佩玲，開啟了我前往地震測報中心新視界的大門；進入地震測報中心之後認識了許多主管、同事，還有在地球物理學術界知名之人物。也要感謝中央氣象局地震測報中心的 GDMS 資料庫，供應資料給與我們研究。進入地震中心這三年來，我因為無法適應頻出狀況，有大家的包容和協助，我才能在身體病痛不斷的情形下，撐著寫完這篇論文。最後，我想感謝的人，是一路上陪伴我的男友。我在就讀研究所之前，因家庭因素與家人發生糾紛，賭氣之下離家出走獨自在外居住，從教師實習期間到就讀研究所，僅靠著實習津貼與教學研究助理費用生活，他資助我這段期間的所有生活費用，讓我能專心於研究上；有時也要陪同我就醫，包容我生病期間發作的情緒。同時也在這段期間，我認識了男友的親戚家人們，他們接受我了的情況並偶然帶來物品給我，一路上不離不棄。一路走下來經歷的風風雨雨之後，想對你說：有你真好。 V.

(7) 目次摘要 ....................................................................................................................... I 致謝 ..................................................................................................................... IV 目次 ..................................................................................................................... VI 圖目錄 ............................................................................................................... VIII 表目錄 ................................................................................................................. XI 第一章. 前言................................................................................................... 1. 1.1 研究動機 ............................................................................................................ 1 1.2 文獻回顧 ............................................................................................................ 5 1.3 研究目的 ............................................................................................................ 6 1.4 章節內容 ............................................................................................................ 7 第二章. 研究方法 ........................................................................................... 8. 2.1 P 波振幅取得 ..................................................................................................... 8 2.2 AR-AIC 震波到時判斷........................................................................................ 8 2.3 複數解載 .......................................................................................................... 12 2.4 消除路徑效應之「振幅解算」 ...................................................................... 15 2.5 由模數(MODULE)重建震源加速度函數 ........................................................... 19 2.6 震源平均函數 .................................................................................................. 23 第三章. 資料處理 ......................................................................................... 25. 3.1 資料來源 .......................................................................................................... 25 3.2 資料選取與統計 .............................................................................................. 27 3.3 處理程序 .......................................................................................................... 32 VI.

(8) 第四章. 結果與討論 ..................................................................................... 36. 4.1 集集地震四分區比較與討論 .......................................................................... 36 4.2 與其他集集地震序列的比較 .......................................................................... 51 第五章. 結論................................................................................................. 66. 參考文獻 ........................................................................................................ 68 附錄 A 一維 AR-AIC 之證明......................................................................... 72 附錄 B 本研究之 QP 模型 ............................................................................ 77 附錄 C 集集地震與另一餘震之資料表......................................................... 86 C.1 集集地震之資料表 .......................................................................................... 86 C.2 集集餘震(09220014)之資料表 ....................................................................... 96. VII.

(9) 圖目錄圖 1.1、地震記錄受到三種效應所影響：震源、路徑、場址效應(工程科學研究室， 2017) ................................................................................................................ 2 圖 1.2、常見簡單震源之時間函數 ............................................................................. 3 圖 1.3、常見複雜震源之時間函數 ............................................................................. 4 圖 2.1、判斷 P 波「到達時刻」，此時間點並非開始偏離平衡點之時刻；而是第一個出露於背景值之時刻。 ............................................................................ 9 圖 2.2、使用前田直樹之 AR-AIC 判斷到達時刻，每一個區域極小值(LOCAL MINIMA) 皆代表一個相位到達而引起之變化；若擷取範圍適當，同樣方式也可以判斷 S 波或表面波到達時刻 .............................................................................. 12 圖 2.3、本圖由上到下三部分分別顯示頻率域振幅之處理，包含截取中央頻率，擷取兩倍振幅，以及中央頻率之平移，完成複數解載的過程。 .............. 15 圖 2.4、計算 P 波在各區塊走時示意圖（陳，1993） ........................................... 18 圖 2.5、Q 模型示意圖；A、B、C、D 表示測站；數字表示區塊編號；星號表示震源位置（陳，1993）...................................................................................... 18 圖 2.6、將每個中央頻率的時間函數 AI(T)乘以 EXP(ΠFI T/Q)，所獲得之各頻率震源加速度時間函數，頻率越高，解算前後的倍數放的越大。 ...................... 22 圖 2.7、各個單一中央頻率波形函數從 4HZ、4.5HZ…至 10HZ 為止，疊加上去所解算出來的震源時間加速度函數。 .................................................................. 22 圖 2.8、同一分區內各測站的震源加速度時間函數，進行相對極大值對齊；再對每一條數值進行平均，大於平均值兩個標準差之外的紀錄挑去，進行 1~2 次加總平均。 .................................................................................................. 24 圖 3.1、微地動、P 波與表面波所涵蓋的主要頻率範圍 ........................................ 27 圖 3.2、收錄不良之 TCU101 三分量與總向量波形圖，初達 P 波波相不明顯 .... 28 圖 3.3、收錄不良的 CHY062 三分量與總向量波形圖，僅收錄到 S 波之後記錄 29 VIII.

(10) 圖 3.4、收錄不良的 HWA053 三分量與總向量波形圖，單一分量無記錄........... 29 圖 3.5、收錄不良的 TTN009 三分量與總向量波形圖，無 P 波到達前時間段 .... 30 圖 3.6、收錄良好 TCU067 三分量與總向量波形圖，有明顯 P 波到達相位 ........ 30 圖 3.7、地震資料處理流程圖 ................................................................................... 33 圖 3.8、經過相對極大值對齊、一致化上下動之後，所獲得之台中分區平均震源加速度時間函數圖。 ...................................................................................... 34 圖 3.9(續圖 3.8)、去除離群加速度時間函數之後所獲得之台中分區平均震源加速度時間函數圖 .................................................................................................. 35 圖 4.1、台中分區(TCU)各測站震源加速度函數，經過相對極大值對齊、平均挑選過後的平均分區震源加速度時間函數圖 ...................................................... 37 圖 4.2、台中分區(TCU)與集集地震震央分布圖 ...................................................... 38 圖 4.3、嘉南分區(CHY)各測站震源加速度函數，經過相對極大值對齊、平均挑選過後的平均分區震源加速度時間函數圖 ...................................................... 40 圖 4.4、嘉南分區(CHY)測站與集集地震震央分布圖 .............................................. 41 圖 4.5、花蓮分區(HWA)各測站震源加速度函數，經過相對極大值對齊、平均挑選過後的平均分區震源加速度時間函數圖 .................................................. 43 圖 4.6、花蓮測站(HWA)與集集地震震央分布圖 .................................................... 44 圖 4.7、台東分區(TTN)各測站震源加速度函數，經過相對極大值對齊、平均挑選過後的平均分區震源加速度時間函數圖 ...................................................... 46 圖 4.8、台東分區(TTN)測站與集集地震震央關係圖 .............................................. 47 圖 4.9、將台中分區、嘉南分區、花蓮分區與台東分區進行震源平均函數比對；紅色為花蓮分區平均震源函數，藍色為嘉南分區平均震源函數，綠色為台東分區平均震源函數，桃色為台中分區平均震源函數。 .......................... 50 圖 4.10、將嘉南分區、花蓮分區與台東分區之震源加速度函數放大比較 ......... 50 圖 4.11、集集餘震花蓮分區(HWA)各測站震源加速度函數，經過相對極大值對齊、 IX.

(11) 平均挑選過後的平均分區震源加速度時間函數圖 ...................................... 52 圖 4.12、花蓮分區測站與兩事件分布圖 ................................................................. 53 圖 4.13、集集餘震台中分區(TCU)各測站震源加速度函數，經過相對極大值對齊、平均挑選過後的平均分區震源加速度時間函數圖 ...................................... 55 圖 4.14、台中分區測站與兩事件分布圖 ................................................................. 56 圖 4.15、集集餘震台東分區(TTN)各測站震源加速度函數，經過相對極大值對齊、平均挑選過後的平均分區震源加速度時間函數圖 ...................................... 58 圖 4.16、台東分區測站與兩事件分布圖 ................................................................. 59 圖 4.17、集集餘震嘉南分區(CHY)各測站震源加速度函數，經過相對極大值對齊、平均挑選過後的平均分區震源加速度時間函數圖 ...................................... 61 圖 4.18、嘉南分區測站與兩事件分布圖 ................................................................. 62 圖 4.19、09220014 花蓮分區、台中分區、台東分區與嘉南分區之平均震源加速度時間函數之比較 .......................................................................................... 63 圖 4.20、將集集餘震之花蓮分區、台東分區、嘉南分區之震源加速度函數放大進行比較 .......................................................................................................... 64. X.

(12) 表目錄表 3.1、中央氣象局各種強震儀的種類及特性比較 ............................................... 26. XI.

(13) 第一章. 前言. 1.1 研究動機台灣位於環太平洋地震帶上，受到歐亞板塊與菲律賓海板塊互相碰撞的影響，1991 年至 2010 年間的每年平均無感地震為 10271 次，有感地震為 658 次（中央氣象局，2014）。預防地震所造成之破壞、減低震災的損失，就一直是政府的施政重點。也因此，了解何處容易釋放較大能量及何處容易發生地震是相當重要的研究方向。地震防災規畫的核心問題架構在準確的強地動值（包括最大加速度值、加速度反應譜與強地動延時等參數）。最終目的在於預估可能發生具有破壞性地震的位置與規模，模擬各地的強地動特性，可以提供參考數據予地震工程與建築設計訂定標準之所需。地震災害的程度，與當地所受到的地動加速度大小相關，而影響地震波傳至地面時地動大小的因素有三個，包括震源特性、路徑傳播效應及場址效應（圖 1.1）。當一個地震發生時，因破裂特性不同而造成不同規模大小地震，是為震源效應；由震源產生之震波在地層中傳遞時會隨傳播之路徑而衰減，是為路徑傳播效應；當震波到達測站局部場址時，因局部地層結構之不同，而對當地測站所產生的地震紀錄振幅大小、時間長短時，需考慮局部之場址效應。所有這些因子加上水文、地質、山崩潛勢等資料庫之建立，由細部的模型建立及模擬， 1.

(14) 即可達到評估區域地震災害潛勢之目的。而不管是由震源特性來估計場址地面振動，或是由地表振動紀錄來解算震源特性，皆必須以路徑效應當作修正的要素。台灣歷年來利用實體波及表面波走時來推求地球內部速度構造之研究已有很好的進展(Rau & Wu, 1995 ; Ma et al., 1996)。如果已知場址效應及路徑效應，則由震源規模大小及位置等已知條件將不難推估地表地動情況。進而假設該地區可能發生之最大規模及位置之地震，即可推估其最大地表地動之空間分佈狀況。而逆推三維衰減構造之方法也已發展成熟（Chen et al., 1996; Chen, 1998）。近年得利於寬頻帶，以及高動態範圍的地震儀及強震儀在國內外大量設置，許多不同大小和遠近的地震都可完整的記錄下來，提供了進行振幅研究的寶貴素材。. 圖 1.1、地震記錄受到三種效應所影響：震源、路徑、場址效應(工程科學研究室， 2017). 在一般情況下，震源會同時產生疏密波(pressure wave)與剪力波 (shear wave)之時間函數，因為本研究是由最先到達測站之震波解算出震源加速度時間函數(Source Accleration Time Function)，故在此僅 2.

(15) 討論疏密波之震源加速度函數。僅有單一破裂過程，產生是簡單的震源(Simple source)時，其震源時間函數通常具有如圖 1.2 所示之情況，然而當震源甚為複雜時，如圖 1.3 所示，代表一震源可能含有三次的震源事件(Source events)。因此如果能將震源時間函數求出，對於震源破裂過程（包括破裂位移、破裂速度、加速度時間函數）將可迎刃而解。. 圖 1.2、常見簡單震源之時間函數. 3.

(16) 圖 1.3、常見複雜震源之時間函數. 由已知震源函數，加上路徑效應及場址效應，進而預估出地表之地動場大小之研究已由陳(2009)完成。然真正震源函數，尤其是震源時間函數，常隨著應力型態及地體構造之不同而差異甚大。台灣地區之孕震帶非常複雜，震源時間函數之型態可能及多樣性，甚難掌握。因此，在地面地動值之預估上，將造成不可預期之誤差。因此對於地震震源時間函數之推算及分析，將顯得舉足輕重。本研究將使用三維速度模型(Roecker et al., 1987)計算震源走時及路徑，利用三維衰減構造因子 Qp（陳，1993）計算需消除之路徑效應；利用複數解載(CMD)計算測站地面振動內各種頻率之絕對能量。由 Qp 與 CMD 這兩種方式，復原出震源加速度時間函數。並將同一分區內的加速度函數進行峰值對齊，以單一分區內多站平均法，平均 4.

(17) 化分區內場址效應與儀器響應，進而求得地表剩餘能量比，藉此求得各測站平均加速度函數。. 1.2 文獻回顧陳(1993)曾經對台灣地區三維區塊模型進行 P 波之 Q 值模型建立，當時假設在 4~25 Hz 震波頻率內，其 Qp 值在此範圍內將維持一定之常數值，Qp 值不會隨著頻率變化而變動，因此頻率最低限為 4 Hz，而 Chopra S. (2011)等人在印度 Gujarat 地區研究 Q 值時，將 Qp 研究頻率範圍限制在 10 Hz 以內，對應的條件是震矩規模在 2.0 以上，規模越大，會伴隨著偏向低頻的 P 波與 S 波；另外以赤澤隆士(Azakawa, 2004a, 2004b)提出自動判別 P 波到達時刻之方式，首先針對一段時間序列進行 5-8 Hz 濾波，而後再搭配波形包絡函數(envelop function)與波形函數進行短時窗平均與長時窗平均比值(STA/LTA)比值和赤池訊息準則(AIC)計算，以判斷出 P 波到時。由這三篇文獻，本論文探討 P 波震源函數時，其頻率範圍將限制在 4-10 Hz 之間。另外劉(2009)曾經利用複數解載法，探討竹子湖強震站與周邊鄰近測站場址效應所造成的特定頻率放大現象；王(2011)也曾利用相同方法探討不同地震事件反應在地動加速度譜值異常的仁愛國小測站。此兩篇文獻皆利用了複數解載(CMD)針對不同頻率地震波，在不同時刻之振幅變化特性分析。 5.

(18) 1.3 研究目的本研究主要使用中央氣象局強地動觀測計畫 (Taiwan Strong Motion Instrumental Program, TSMIP)之中，GMT 時間 1999 年 9 月 20 日、震矩規模(Mw)7.6 之集集地震；採用已知全台灣速度構造與三維衰減構造（陳，1993；Chen et al, 1996; Chen, 1998），挑選最靠近集集地震之台中(TCU)分區測站，加上距離稍遠但測站方位角較為集中之嘉南(CHY)分區測站，還有橫跨中央山脈地帶之花蓮分區(HWA)與跨越海岸山脈之台東分區(TTN)，由這四個分區強震站所有收到地震記錄的測站之中，進行複數解載。複數解載的方式是將一整段時間域的紀錄進行複立葉變換至頻率域，在頻率域記錄中，對於進行逆複立葉變換的取樣頻率在 4-10 Hz 範圍內、每隔 0.5 Hz(4 Hz、4.5 Hz、5 Hz…10 Hz)、頻寬為中央頻率左右各 0.5 Hz、經過帶狀濾波後進行逆變換獲得單一頻率的時間域記錄，即加速度時間函數。而後針對各個測站於 P 波到達後 10.24 秒內進行複數解載(CMD)，找出中央頻率進行逆變換所獲得的單一頻率之加速度時間函數；以取樣率 200~250 點／秒的儀器來說，10.24 秒紀錄點數為 2 的 n 次方之整數倍。並針對震波到達測站後 10.24 秒內，疊加各頻率的時間加速度函數，進而求得 4~10 Hz 內單一測站之震源加速度時間函數。對於同一分區的測站進行峰值疊加，並計算各測站震源加速度函數之平均值，即獲得單一分區之加速 6.

(19) 度時間函數。製作集集地震與餘震之一的各分區震源加速度函數，同一事件進行各分區比較，再對兩個地震事件進行震源加速度時間函數比較。由兩次不同震源效應的地震事件，來證明計算震源加速度時間函數的可行性。. 1.4 章節內容本論文內容共分為五章，第一章為研究動機、目的與內容簡介；第二章研究方法介紹震源加速度函數的解算方式，從 AIC 計算法、使用逆推之 Qp 值到振幅解算過程；第三章資料處理介紹，從資料來源、資料篩選方式、到組合前一章所敘述之計算方式來進行解算；第四章為本論文對於台中分區、嘉南分區、花蓮分區與台東分區之研究結果，並加入另一事件之同樣四個分區來做比較；第五章為本研究所獲得之結論。. 7.

(20) 第二章研究方法 2.1 P 波振幅取得由中央氣象局所提供之強地動觀測計畫所提供的強地動測站紀錄為時間、垂直向分量、南北向分量、與東西向分量之加速度振幅紀錄。本研究為了求得與傳遞方向平行之 P 波質點運動方向(Partical motion)，首先將三個分量計算向量和，求得振幅之純量；再由垂直向之正負，決定振幅之正負值，進而還原成 P 波質點運動振幅之總量。. 2.2 AR-AIC 震波到時判斷判斷一個波相之「到達時刻」或「結束時刻」如圖 2.1 所示，並非地震記錄開始偏離平衡點之時刻，而是記錄開始從背景值冒出之時刻；相對地，「結束時刻」之判斷即是最後一個時間點沒入背景值之時刻。一般來說藉由人工判斷到達時刻仍有數秒內之誤差，此時數值化的 AR-AIC 算法可彌補此問題，提升判斷的精準度。. 8.

(21) 圖 2.1、判斷 P 波「到達時刻」，此時間點並非開始偏離平衡點之時刻；而是第一個出露於背景值之時刻。（資料來源：Yoshida, et al., 1994）. 為了判斷計算震源函數的範圍，我們在此使用赤池資訊準則 (Akaike Information Criterion)，其方式是由赤池弘次所發展出的統計學計算法(Akaike, 1979)，其方式是針對某一段時間序列，假設可被分成背景值部分與訊號部分（包含波相到達時刻）等兩個區間。則 AIC 值計算方式如下：. . . ˆ , ˆ | x  2  estimated parameter AIC  2  log L K , M ,  1 2  K  M log  12,max  N  M  K log  22,max 2.  N  2M log 2  1  2 mi  1 i 1. (2.1). 在方程式(2.1)中，N 為時間序列總點數；K 為將一段時間序列分成兩部分計算的分隔點(driving point)；m1 與 m2 為背景值與訊號值之 9.

(22) 自迴歸階數(Auto-regressive Order)，有時我們會將兩者的階數統一為 M 個；而 σ1,max2 與 σ2,max2，是有著理想白噪訊特性的背景值殘差項與訊號值殘差項之變異數，當 log[L(K, M, Θ1, Θ2| x)]獲得最大概似估計值 (maximum approximate likelihood value)時，則 σ12 與 σ22 皆會得到最大值 σ1,max2 與 σ2,max2。由上式結果中，每個時刻對應計算出的 AIC 值，彼此之間有著函數關係時，任一時刻 AIC 之區域極小值(local minima)，即為任一波相開始或結束之時刻。在本論文中我們僅介紹判別 P 波到時方式，因此在 AIC 函數上，最接近 P 波到達時刻之相對極小值，即為 P 波到達時刻。事實上只要用以計算的時間點數擷取範圍恰當，同樣方式也可延伸到判別其他波相之到達時間。TSMIP 在正常運作下蒐集地震紀錄的方式，是感應到加速度振幅達到 3.92 gal 臨界值時，觸發紀錄運作，將此段地震連同臨界值出現前 20 秒之前置時間(pre-event)，到達小於臨界值 15 秒後置時間(post-event)一併做紀錄，由此種特性，判斷 P 波到達測站時應該在開始紀錄後 20~25 秒內，因此可以用來判斷較為精確的 P 波到達時間。但由於原有的 AIC 值計算公式涉及到須計算自回歸階數、對應之自回歸係數、以及理想白噪訊之變異數、計算方式較為繁複。因此前田直樹(Maeda, 1985)改進上面計算 AIC 值的公式，修改新的計算方式如下： 10.

(23) AIC  K logvar1, K   N  K 1logvarK  1, N   K  log  12,max  N  K  1log  22,max. (2.2). 方程式與(2.1)與(2.2)兩者，計算 AIC 的方式各有著相同與相異的特性：共同的部分是，當背景值與訊號值之分隔點正好在波相到達時，皆會使計算變異數 σ1,max2 與 σ2,max2 的部分獲得最大值，而後再對於似然估計值(approximate likelihood value)取對數化數值時，使 AIC 達到最小值；不同的部分是，Maeda 的方法可以直接針對原始之地震訊號分段進行計算，而不需設定自回歸階數 mj，也不需計算自回歸係數 aj、白噪訊變異數 σ12 與 σ22 …等數值。針對震源函數之分析，將 AIC 值的計算引入其中，則相對於人工判讀大量地震事件之震波到時，更能大幅提高判讀 P 波到時之精準程度。一維 AR-AIC 之詳細證明過程如附錄 A。. 11.

(24) 圖 2.2、使用前田直樹之 AR-AIC 判斷到達時刻，每一個區域極小值(local minima) 皆代表一個相位到達而引起之變化；若擷取範圍適當，同樣方式也可以判斷 S 波或表面波到達時刻. 2.3 複數解載複數解載為時序資料分析方式之一，藉由擷取特定頻寬範圍之頻率值，以分析某一特定頻率分量之振幅以及相位隨時間變化之函數值，整個複數解載過程所使用到的處理工具主要有：快速傅立葉變換(FFT)、帶狀窗型濾波(bandpass filter)、能夠避免混疊(alising)效應的餘弦鐘形濾波(cosine bell filter)、頻率平移、以及傅立葉逆變換(Inverse FFT)。複數解載程序，應用於時間序列 X(t)（指每隔一定時間間隔所記錄之離散且週期性之數值），由兩項簡易數學運算所定義。首先，在頻譜上的頻帶，以一個複數指數函數 Exp （而ω是位移帶的中央頻率），乘以時間序列的變數位移到零的頻率。一個新的序列是透過每個頻帶 12.

(25) 所產生的：. X s  ' , t   X t exp  i' t . (2.3). 其次，使用一組權重序列 ak (k=-m~+m)經低通濾波所產生位移的頻率序列，將此頻譜進行傅立葉逆變換，轉換到時間域中，此序列正是帶有解載頻率之時間序列函數 Xd(ω’, t)：. X d  ' , t  . k  m.  a X ', t  kt . k  m. k. s. (2.4). 解載函數可以用更簡便的形式表示 X d  ' , t   X d  ' , t  exp  i ' , t . (2.5). 用 Xd(ω’, t)中之模數(Modulus)和相位(phase)來表示每一項之結果。讓我們假設一個含有週期性的時間序列，其可在傅立葉頻譜上的頻率 ω0 產生峰值：. X t   At cos0t   . (2.6). 此簡單的複數解載的中心在頻率 ω’=ω+δω 成為.  At  i 0t t  i (0t  ) i (0t  ) X s  ' , t    e e e  2 . .  At  i t   i20  t    e e  2  . . .  (2.7). 這些包含了頻率-δω 和 – (2ω0 +δω) 位移的頻率序列再經由低通濾波，如期地把較高頻的(-2ω0 + δω) 13.

(26) 完全去除，剩下.  At  i t   X d    t    e  2  . (2.8). 在此我們假設濾波器有單一反應，並且在頻率-δω 上沒有相位的平移。如果我們知道資料含有頻率 ω0 的週期性變化，那麼中央頻率自然而然選擇 ω’=ω0 作為複數解載所在的中央頻率，接著頻帶平移到零（頻帶右翼移至頻譜左側；頻帶左翼移至頻譜右側），頻帶隨即反方向至 ω0 處。在這個例子中，δω=0：.  At  i X d 0 , t    e  2 . (2.9). 由於：.  At  i  At  e  R  iI   2   2 cos   i sin    At   At   cos   i   2  sin   2   . (2.10). 因此： At  . 2R 2I  cos  sin . (2.11). 如下圖 2.3 的步驟，說明將時間序列進行 FFT 獲得頻率序列之後，. 14.

(27) 進行 CMD 處理，將各頻率的時間函數對齊，以等高線之二維圖形得到三維函數之時頻圖。圖 2.3 中，左半邊為頻率域的振幅分布情形，右半邊為進行變換所對應的時間域振幅。其中黑線部分為複數指數函數，乘以位移至 0 Hz 之前的時間域序列；紅線為經過低通濾波，將中央頻率位移至 0 Hz 之時間域序列。. 圖 2.3、本圖由上到下三部分分別顯示頻率域振幅之處理，包含截取中央頻率，擷取兩倍振幅，以及中央頻率之平移，完成複數解載的過程。. 2.4 消除路徑效應之「振幅解算」簡諧運動(Simple Harmonic Motion)表示的是震波每經一個週期並沒有振幅的增加或減少，但自然界中每個週期振動的能量與振幅，受到阻尼或外加力量影響，皆會有不同程度的減少或增加。僅受到阻 15.

(28) 尼影響時，每經過一個振動週期，則會衰減一定的能量；每個衰減週期所減少的能量，與前一個週期過後所剩餘之能量呈現正比關係，如下(2.12)所示： . dE 1  E dt 2f. (2.12). . dE 1   kE dt 2f. (2.13). 寫成完整的微分形式. 由於物質之滯彈性(anelasticity)，能量衰減後剩餘值 E’與衰減前量值 E 之比值(E’=kE)介於 0 至 1 之間(k: 0~1)，為一固定比值，若將比值寫成分母形式，則倒數形式會使衰減因子係數 Q 介於∞至 1 之間。故微分形式表示法可改寫如下： . dE 1 1   E dt 2f Q. (2.14). 使用指數函數積分，計算能量衰減函數： dE 2f   dt E Q 1.  E dE   . 2f 2ft dt  ln E    const. Q Q.  E t   k  e. . 2ft Q.  E0 e. . 2ft Q. (2.15). 對於一個距離測站為 r 的某地震之加速度譜，可定義為(Hough et 16.

(29) al., 1988) Ar, f   2f  S  f Gr, f e ft* 2. (2.16). 方程式(2.16)中的 G(r,f)為幾何展開(Geometric spreading)，S(f)為震源之加速度譜，其中之 t*，可用如下之數學式表示： t dt 1 dr   0 QV 0 Q. t*  . t. (2.17). 若已知震波經過多個區塊路徑，則可改寫方程式(2.17)之積分形式，將 t*改以不連續之加總函數表示：. t* . n t t t1 t 2   ...  n   k Q1 Q2 Qn k 1 Qk. (2.18). V 為震波速度，Q 為衰減因子，表示為一個週期性波動或振動在每個週期的能量損耗率。已知能量衰減前後的倍數關係，且振幅之平方值與能量成正比. A2  E  A  E ，則衰減前後之振幅比值，可從衰減前後之能量比值求得  ai E0  E   e Ai E0. 2f i t Q. e. . f i t Q.  e fit*. (2.19). 如果已知某地震之震源位置，利用波線追跡法，找出 P 波波線路徑，計算 P 波線在各區塊內走時(如圖 2.4 所示)，並該波線所經 Q (簡化之 Q 模型如圖 2.5)代入(2.19)計算出單一頻率 fi 之 17. ai 值。 Ai.

(30) 圖 2.4、計算 P 波在各區塊走時示意圖（陳，1993）. 圖 2.5、Q 模型示意圖；A、B、C、D 表示測站；數字表示區塊編號；星號表示震源位置（陳，1993）. (2.16)式中，因為計算理論地動峰值所使用的觀測值為同一地震，故震源效應項 S(f)可消去；幾何展開項 G(r,f)在小範圍內也可消去。本研究消除路徑效應的方式，是採用三維衰減構造 Qp 模型（陳，1993） 18.

(31) 及速度構造(Roecker et al., 1987)在 4-25Hz 振幅範圍內，對於 P 波到達時的速度區塊模型進行 t*值的計算；本論文內採用 4-10Hz 範圍，進一步求出各區塊所對應之 Q 值。利用三維衰減構造模型所求出的 Qp 值；同時配合初始之資料庫格式所定義之 Qp>0，以四周區塊 Qp 內插法計算出同一深度層中其他缺少 Qp 數值區塊之 Qp 值，使用已經求得內插法之三維衰減構造模型，配合速度構造在台灣地區的區塊分割方式，計算已知經緯度、深度之震源，並配合已建立參數的測站（包含名稱、經緯度、海拔），來計算出本次地震傳輸到各地測站之 tp* 值。將此 tp*值代入簡化過的方程式(2.16)中，求出各頻率在各時間點的地動加速度值。計算出的 tp*值配合前一節之 AIC 到時判斷法，與下一節要探討的波形函數重建來配合，以計算出震源加速度函數。Qp 之三維衰減構造模型列於附錄 B。. 2.5 由模數(Module)重建震源加速度函數對於一段長度為時間之序列 X(t)，實際進行 FFT 計算時，通常都會搭配一段與 X(t)相同點數，但振幅為 0.0 之實數序列，一同放入程式之實部與虛部進行計算，最常見的做法是把時間序列 X(t)值放入 FFT 之實部、振幅值 0.0 序列放入虛部，進行因式分解、排序轉換等 FFT 之計算流程。 19.

(32) 假設震源時間函數為 S(t)，且考慮進行 CMD 之後分成實部與虛部，各對應餘弦函數與正弦函數之模數(Modulus)，則復原成震源時間函數時，其 S(t)如下所示 S t   Areal,1 t cos1 t   Areal, 2 t cos2 t   ...  Areal,n t cosn t  n.   Areal,i t cosi t . (2.20). i 1.   Aimage,i t   2 2 1   Ai t   Areal,i t   Aimage,i t  i t   tan    Areal,i t    其中 Areal,i 與 Aimage,i 各為解算至震源、第 i 個頻率所對應之實部振幅與虛部振幅；而φ為相位角。其中 Ai 為第 i 個 CMD 頻率所對應之震源振幅向量和，現在令測站之地震紀錄為 R(t)，則測站之震源時間函數如下所示： Rt   areal,1 t cos1 t   areal, 2 t cos2 t   ...  areal,n t cosn t  n.   areal,i t cosi t . (2.21). i 1.   aimage,i t   2 2 1  ai t   areal,i t   aimage,i t  i t   tan    areal,i t   . 其中 areal,i 與 aimage,i 各為在單一測站、第 i 個頻率所對應之實部振幅與虛部振幅；而φ為相位角。其中 ai 為第 i 個 CMD 頻率所對應之測站振幅向量和。 20.

(33) 前一節所探討在單一頻率 fi 下，震源振幅與測站振幅之關係式，如(2.22)所示。利用複數解載(Complex demodulation) 可以將位於測站 a1(t)、a2(t)……an(t)之絕對值求出，將其值代入(2.23)式可以求得位於震源時的 A1(t)、A2(t)……An(t)，即. ai t   Ai t   e. . f i t Q. (2.22). m f i t  t  f t  Ai t   ai t   e Q  ai t   exp  i   ai t    f i  k  Q   k 1 Qk.   . (2.23). 若 P 波波線經過 m 個區塊(block), 其 Qp 值分別為 Q1, Q2, …,Qm ，而在該對應之區塊內之走時為 t1, t2, …, tm, 。對應特定範圍之中央頻率 f1, f2, …fn，則.  n  m tj S (t )   Ai (t )  ai t   exp   f i   i 1  i 1  j 1 Q j n.    . 即可獲得震源加速度時間函數為 S(t)（公式(2.20)、(2.24)）. 21. (2.24).

(34) 圖 2.6、將每個中央頻率的時間函數 ai(t)乘以 Exp(πfi t/Q)，所獲得之各頻率震源加速度時間函數，頻率越高，解算前後的倍數放的越大。. 圖 2.7、各個單一中央頻率波形函數從 4Hz、4.5Hz…至 10Hz 為止，疊加上去所解算出來的震源時間加速度函數。 22.

(35) 2.6 震源平均函數本資料來源為中央氣象局的強地動觀測計畫而來，由所收錄到的地震事件進行處理，每個震源事件傳播到各測站後，必定會根據頻率大小，產生不同相位、不同振幅的時頻函數。同一個地震事件紀錄裡，經過 CMD 分割解算過後，儀器響應影響了各頻率與相位的函數值，故各頻率之時間函數其振幅大小、與相位正負，可能會與原來之時間函數出現上下相反的情形（梁，2016）；在截取 4~10 Hz 進行解算並疊加之後，也必定會對震源加速度函數之初動記錄有所影響。以台灣地區七大分區為界線，為了使場址效應與儀器效應在同一個分區內平均化，並減少初始滑動抵消掉震源平均函數的計算值。在此我們一律將震源加速度函數改為相對極大值之曲線函數，並對於同一分區內測站紀錄之波峰進行對齊，平均化同一個分區內所有測站之場址與儀器效應，計算各分區之平均加速度函數。. 23.

(36) 圖 2.8、同一分區內各測站的震源加速度時間函數，進行相對極大值對齊；再對每一條數值進行平均，大於平均值兩個標準差之外的紀錄挑去，進行 1~2 次加總平均。. 24.

(37) 第三章資料處理 3.1 資料來源為了深入研究地震，中央氣象局自民國 91 年 7 月起，配合六年國建，實施強地動觀測計畫，將全台灣劃分成七大區域，分別為台北 (TAP)、台中(TCU)、嘉南(CHY)、高屏(KAU)、宜蘭(ILA)、花蓮(HWA)與台東(TTN)地區，並設置了七百多部自由場強震儀觀測站（中央氣象局地震測報中心，2009）。本研究取用自由場強地動觀測網，所用的紀錄儀器型號包括 A800、A900、A900A、IDS3602、IDS3602A、K2、ETNA、CV-574C、CV-575C 以及 SMART24A 等十種強震儀；除 IDS 系列取樣率為 250 點／秒以外，其餘取樣率皆為 200 點／秒；除 A800 最大紀錄為正負 1G 範圍，其餘地震儀皆可記錄至最大正負 2G 範圍之訊號。強震儀以觸發模式來判定地震是否發生，當強震儀感應到最大振幅的千分之二時（±3.92 gal，相當於震度二級）即觸發儀器將此次事件紀錄下來，直至波動振幅下降至臨界值（約±3.92 gal）才會停止；在訊號觸發儀器紀錄時，觸發前 20 秒之前置時間(pre-event)與下降至臨界值後之 15 秒之後置時間 (post-event)紀錄也一併寫入記憶體中，包含大部分事件之主相位歷時，各測站之地震紀錄，全長可達 90 秒甚至 120 秒；而每個測站配合強震儀所設置之 6 MB 靜態記憶體，在正常運作下，最多可紀錄高達 58 25.

(38) 筆地震訊號。TSMIP 所使用之強震儀，其詳細特性及列表見下表 3.1 所示。表 3.1、中央氣象局各種強震儀的種類及特性比較加速度儀特性製造商. Teledyne Geotech. 型號. Full Scalea. 解析度. 記憶體. 記錄時間. 取樣率. 動態範圍. 即時輸出. A800. ±1G. 12-bit. 1 MB. 14 MIN. 200/s. 72 db. No. A900. ±2G. 16-bit. 6 MB. 200/s. 96 db. No. A900A. ±2G. 16-bit. 6.5 MB. 90 MIN. 200/s. 96 db. Yes. SMART24A. ±2G. 24-bit. 插卡. 插卡. 200/s. 120db 以上. Yes. IDS-3602. ±2G. 16-bit. 2 MB. 23 MIN. 250/s. 96 db. No. IDS-3602A. ±2G. 16. 4 MB. 90 MIN. 250/s. 96 db. Yes. SSA-16. ±2G. 24. 2 MB. 28 MIN. 200/s. 96 db. No. K2. ±2G. 18. 插卡. 200/s. 114db. Yes. ETNA. ±2G. 16. 2 MB 插卡. 插卡. 200/s. 96 db. Yes. CV-574C. ±2G. 16. 插卡. 插卡. 200/s. 96 db. Yes. CV-575C. ±2G. 24. 插卡. 插卡. 200/s. 120db 以上. Yes. Terra Tech. Kinemetric. 日本東京測震株式會社. 87 MIN. a. 最大加速度範圍(Full Scale)一欄裡，G 為重力加速度之單位，1G=980gal；進入西元 2000 年後，IDS 系列之地震儀被逐步更換成取樣率 200 點／秒的其他廠商所製作之地震儀。（資料來源：李錫堤、邱宏智、林柏伸、鄭錦桐(2006)。）. 本研究之地震事件，從交通部中央氣象局 TSMIP 繪出三分量波形挑選。有的測站紀錄未紀錄到前置時間(pre-event)區段，因此無完整的 P 波區段；無 P 波到達前後之地震紀錄，將會被當作收錄不良的 P 26.

(39) 波波形，捨棄不用。而後在所有地震事件中，挑選取樣率為 200 點／秒者，大部分的 Anderson 系列（A800、A900、A900A）、K2、ETNA、CV 系列皆為，僅部分測站在較早期事件使用取樣率 250 點／秒之 IDS 系列（IDS3602、 IDS3602A）。對於取樣率所造成的點數差異，以 SigmaPlot 軟體內附的函數「內插法」(interpolate)功能，將所有點數均調整成取樣率 200 點／秒之地震紀錄。. 3.2 資料選取與統計微地動、P 波與表面波所涵蓋之頻率範圍，依照規模大小之不同而有所不同，如圖 3.1 所示：. 圖 3.1、微地動、P 波與表面波所涵蓋的主要頻率範圍（資料來源：Waves, Seismometers, and Seismograms, 2017）. 測站紀錄資料處理的方式分為以下步驟： 27.

(40) 1.. 去除三分量紀錄異常的測站，對所有測站原始三分量紀錄，進行總量之計算，並以垂直向來決定 P 波質點運動之正負值，求出 P 波質點運動方向之加速度函數。測站紀錄不良的例子如圖 3.2 至圖 3.5 所示，收錄良好的波形圖如圖 3.6 所示：. 圖 3.2、收錄不良之 TCU101 三分量與總向量波形圖，初達 P 波波相不明顯. 28.

(41) 圖 3.3、收錄不良的 CHY062 三分量與總向量波形圖，僅收錄到 S 波之後記錄. 圖 3.4、收錄不良的 HWA053 三分量與總向量波形圖，單一分量無記錄. 29.

(42) 圖 3.5、收錄不良的 TTN009 三分量與總向量波形圖，無 P 波到達前時間段. 圖 3.6、收錄良好 TCU067 三分量與總向量波形圖，有明顯 P 波到達相位. 30.

(43) 2.. 針對有 P 波相位到達前後之加速度函數之紀錄，進行所有測站等頻率、等頻寬之複數解載；本研究從中央頻率 ω0=4Hz 開始，每隔 0.5Hz 取其中央頻率 ω0，則中央頻率 ω0= 4.0、4.5、5.0……、 10.0Hz，取樣至 ω0=10Hz 為止；頻寬為 2Δω=0.5Hz，獲得各個中央頻率函數隨時間的變化值。. 3.. 對單一測站之複數解載，調整至 P 波到達測站時刻為 0.0 秒，在 10.24 秒內，代入陳(1993)對於台灣地區推算之三維 Qp 值模型，進行 P 波振幅解算，由測站一定範圍內各頻率之時間函數，求得各測站在各頻率之時間波形函數，即在震源時的 CMD，換句話說即是求得各頻率時間函數，在震源時所對應之振幅 A(f,t)。. 4.. 將單一測站所有頻率進行疊加，可求得單一測站 P 波到達後 10.24 秒內之震源加速度函數 A(t)。將所有測站疊加後之波峰值進行對齊，計算平均化震源加速度函數。為了平均化同一分區場址與儀器效應之影響；需去除與平均值距離大於兩個個標準差之離群值(outlier，異常值、極端值、偏離值)之異常振幅去除，再次將剩下的振幅進行平均，重新計算平均值，以求得更穩定之平均值；經多次去除處理，可獲得各個分區所對應之穩定振幅平均 (robust mean)A(f,t)。. 5.. 針對 4-10Hz 範圍內之 Qp 值，當作不隨頻率變動之常數進行解算 31.

(44) （張，2004），以去除路徑效應；由於地震儀對於不同頻率之地震波振幅及相位之影響有所不同，還原成震源函數時要將振幅大小以及相位上下相反也考慮進去（梁，2016）；將頻率範圍內的加速度函數疊加後，從到達測站當作時間 0.0 秒解算時，其加速度波形會隨機地往上或往下振動，並在函數時間範圍內產生相對極值。為了將儀器效應、場址效應與剩餘路徑效應在同一分區內平均化，解算震源加速度時間函數之後，將每個測站解算結果之第一個峰值相互疊合，將此峰值當作相對極值，繪出上凸或下凹之曲線函數，作為計算震源加速度平均函數的方式。 6.. 將所有測站以分區為一組，計算平均震源加速度函數過後，計算其他有收到紀錄的分區，將同樣一個地震事件震源加速度函數會出，進行比較。而後再選規模大於六以上，較接近原集集地震震央之地震事件，處理第二個事件進比較。. 3.3 處理程序如下之流程圖 3.7 所示，另外我們列舉集集地震之台中分區為例，示範分區資料處理過程（圖 3.7、圖 3.8、圖 3.9），其他分區的資料處理結果，詳見附錄 C：. 32.

(45) 集集地震事件. No. 波形圖判斷是否. 廢棄不用. 有完整 P 波. Yes 進行 TSMIP 紀錄前. 以 P 波原始向量. 20~25 秒 P 波原始. 計算 AIC 值，判斷. 波形分析. P 波到達時刻. 使用 FFT 找出頻譜；並合併 IFFT、. 獲得 P 波到達各測. cosine bell filter 找. 站之各時刻. 出 CMD. 對 CMD 時間範圍從到達時刻至 10.24 秒、頻率範圍 4~10Hz 內進行解算. 將所有去除路徑效應之加速度函數，進行尖峰值對齊，平均化場址與儀器效應. 針對自由場七大分區最靠近震源四大分區 (TCU、CHY、HWA、TTN) 進行兩震源事件比較. 33. 圖 3.7、地震資料處理流程圖.

(46) 圖 3.8、經過相對極大值對齊、一致化上下動之後，所獲得之台中分區平均震源加速度時間函數圖。. 34.

(47) 圖 3.9(續圖 3.8)、去除離群加速度時間函數之後所獲得之台中分區平均震源加速度時間函數圖. 35.

(48) 第四章結果與討論 4.1 集集地震四分區比較與討論各分區平均震源加速度函數，詳細分析結果如以下對各分區之說明所示。 1.. 台中分區：台中分區是所有分區裏，大部份測站最靠近震央的分區。其中最接近震央的 TCU065（26.13 公里），S 波與 P 波到時差約 3.716 秒；較遠離震央的 TCU064（58.48 公里），S 波與 P 波到時差約 7.959 秒。所有測站並無任何一個測站震央距在 75 公里以上（相當於 S-P 波走時差在 10.24 秒以上），所有測站之 S-P 波走時差皆在 10.24 秒以內。代入三維 Qp 值模型後，各測站震央距、走時差、與分區內各測站震源函數經挑選過後，平均加速度函數波形圖，如圖 4.1 所示；台中測站與震央分布圖，如圖 4.2 所示。. 36.

(49) 圖 4.1、台中分區(TCU)各測站震源加速度函數，經過相對極大值對齊、平均挑選過後的平均分區震源加速度時間函數圖. 37.

(50) 圖 4.2、台中分區(TCU)與集集地震震央分布圖 38.

(51) 2.. 嘉南分區：嘉南分區的測站群，距離震央最近的測站為 CHY088 其震央距有 68.48 公里，S-P 波走時差約 9.673 秒；距離震央最遠的測站為 CHY066，其震央距有 120.39 公里，S-P 波走時差 15.678 秒，大約有前五個測站的 S-P 波走時差在 10.24 秒以內，也就是震央距不到 75 公里。代入三維 Qp 值模型後，嘉南分區各測站之震央距、走時差與解算過後震源加速度函數挑選過後計算而得的平均加速度函數，如圖 4.3 所示；嘉南分區測站分布與震央之關係圖，如圖 4.4 所示。. 39.

(52) 圖 4.3、嘉南分區(CHY)各測站震源加速度函數，經過相對極大值對齊、平均挑選過後的平均分區震源加速度時間函數圖. 40.

(53) 圖 4.4、嘉南分區(CHY)測站與集集地震震央分布圖 41.

(54) 3.. 花蓮分區：花蓮分區的測站，幾乎都沿著海岸線的方向架設。已知全部收錄到事件的 50 個測站中，前 16 個之 S-P 波之走時差在 10.24 秒以內；距離震央最近的測站是較南邊的 HWA032，震央距約 62.78 公里，S-P 波走時差約 8.76 秒，距離震央最遠的測站是北邊的 HWA045，震央距約為 106.83 公里，S-P 波走時差大約 13.549 秒，S-P 波走時差不到 10.24 秒的測站，大部分都集中在花蓮縱谷的中南段。代入 Qp 三維模型後，各測站震央距、走時差、與分區內各測站震源函數經挑選過後，平均加速度函數波形圖，如圖 4.5 所示；花蓮分區測站與震央分布圖，如圖 4.6 所示。. 42.

(55) 圖 4.5、花蓮分區(HWA)各測站震源加速度函數，經過相對極大值對齊、平均挑選過後的平均分區震源加速度時間函數圖 43.

(56) 圖 4.6、花蓮測站(HWA)與集集地震震央分布圖 44.

(57) 4.. 台東分區：台東分區所解算出來的的震源加速度函數，是四個分區內振幅最小的。是因為在台東分區計算時，所有的測站 S-P 波走時差皆在 10.24 秒以上。所有測站因震央距大於 75 公里以上，也就是 S-P 波到時差在 10.24 秒以上，所解算的震源函數幾乎都只有 P 波，沒有 S 波或表面波滲入之影響，是造成震源加速度時間函數振幅偏小的原因。代入 Qp 三維模型後，各測站震央距、走時差、與分區內各測站震源函數經挑選過後，平均加速度函數波形圖，如圖 4.7 所示；台東分區測站與震央分布，如圖 4.8 所示。. 45.

(58) 圖 4.7、台東分區(TTN)各測站震源加速度函數，經過相對極大值對齊、平均挑選過後的平均分區震源加速度時間函數圖. 46.

(59) 圖 4.8、台東分區(TTN)測站與集集地震震央關係圖. 47.

(60) 由以上四個分區之特性，我們用前一章所描述之資料處理方法解算集集地震之震源函數，使用 CMD 解算 10.24 秒內各測站隨時間變化之頻率函數，再代入三維 Qp 值模型，從各測站回推在震源時隨時間變化之時間頻率函數，最後將所有頻率之時間函數疊加起來，獲得震源頻率加速度函數。接著再進行峰值對齊，挑去與平均值誤差太大的震源加速度函數，獲得各分區平均化之震源加速度函數，並進行分區間之比較（圖 4.9）。在台中分區(TCU)震源函數部分，即使經過波形篩選，將場址效應與儀器響應在分區內平均化，前六秒的資料仍然集中於±10gal 以內；而且之後的震源加速度函數變化幅度非常大，在 10.24 秒時震源加速度函數的數據就已達到將近 30gal 左右。最有可能之因素，是因為台中分區測站群距離震央較近，大部分測站的震央距不到 50 公里、S-P 波走時差在 3~8 秒範圍內，緊跟在 S 波之後的表面波也有可能在 10.24 秒範圍之內一起到達。我們使用三維模型 Qp 值解算震源函數，再用峰值對齊進行波形挑選截取平均值時，會連帶地將較後面的 S 波、甚至是更大振幅之表面波也一起納入 10.24 秒內的範圍進行解算，造成台中分區的測站群有著振幅極大之震源加速度函數。嘉南分區 37 個測站之中有 5 個，其 S-P 波到時差在 10.24 秒以內；而花蓮分區 50 個測站中有 16 個，其 S-P 波到時差在 10.24 秒以內。比起台中測站大部分測站的走時差要遠上許多，即使 S 波滲入震源加速度之內被解算，但是並不會將表面波也一起納入解算 10.24 秒的範圍。造成嘉南分區震源加速度函數如此變化，最主要的原因是最靠近震央測站 S-P 波的走時差有 9.673 秒，這樣的結果雖然仍會有 S 波之滲入進行解算震源平均加速度函數，但表面波並沒有大量地滲入 48.

(61) 10.24 秒範圍內，造成震源加速度函數內一起被解算。故兩個分區整段震源加速度函數的數值，大約在±4 gal 的範圍內進行振動，兩個震源加速度函數之數值極為相似，即使將整段函數放大到上下限接近± 4 gal 範圍內，仍可看的出來，兩分區的加速度誤差僅有不到±1 gal 之差異。只在 2~7 秒的範圍，花蓮分區內的振幅稍為比嘉南分區測站之振幅大一點點，之後將近 10 秒，嘉南分比起花蓮分區也只有大一點點的振幅差距，不到± 1gal 之變化。與台中分區比較，並沒有像台中分區那樣有±10 gal 到±20 gal 以上的加速度函數變化。花蓮分區震央距中位數約為 77.05 公里，走時差中位數約 10.648 秒；嘉南分區震央距中位數 96.2 公里，走時差中位數 12.639 秒，嘉南分區卻與花蓮分區的震源加速度函數振幅相差不遠，有可能因挑波或其他特定因素使得兩者振幅差異不大，具體成因仍有待進一步研究。台東分區因為所有的測站震央距皆在 75 公里以上，意味著此分區之 S-P 波走時差皆在 10.24 秒以上，所有的解算結果皆只有 P 波之震源加速度函數，在到達測站後約前五秒以內，與嘉南分區和花蓮分區相似，大約有 2 gal 的振幅值；但在四、五秒過後，台東分區的震源函數因為無 S 波與其他表面波的干擾，解算過後的振幅值大約為±1gal，比起嘉南或花蓮分區的震源加速度函數振幅更為平緩。為了進一步分析震源加速度函數細部，除台中分區之外，針對其餘三個分區的加速度函數細節部分，另外繪出去除台中測站曲線圖形，將剩下的三個分區進行放大比較（圖 4.10），會發現嘉南與花蓮地區震源加速度函數振幅非常近似，僅有波峰位置稍有不同，兩者之間的相關係數大約是 0.98；台東分區之平均震源加速度函數，因為僅有 P 波被解算回震源，故在震源加速度函數振幅數值最小也最平緩。 49.

(62) 圖 4.9、將台中分區、嘉南分區、花蓮分區與台東分區進行震源平均函數比對；紅色為花蓮分區平均震源函數，藍色為嘉南分區平均震源函數，綠色為台東分區平均震源函數，桃色為台中分區平均震源函數。. 圖 4.10、將嘉南分區、花蓮分區與台東分區之震源加速度函數放大比較. 50.

(63) 4.2 與其他集集地震序列的比較為了證明此方法的可行性，我們使用集集餘震之一、大於規模 6 以上地震餘震事件當做對照組，再次進行震源加速度函數之計算。本次餘震與集集大地震相比較，發現此餘震的位置有以下特性： 1.. 花蓮分區距離此餘震位置最近，花蓮分區最接近餘震震央的測站為 HWA032，震央距 39.38 公里，S-P 波走時差 5.855 秒；離餘震最遠的測站為 HWA045，震央距 88.57 公里，S-P 波走時差 11.003 秒。全部測站之中，只有最遠的測站 HWA045，S-P 波走時差在 10.24 秒以上。代入 Qp 三維模型後，各測站震央距、走時差，如圖 4.11 所示；各測站與震央之分布圖如圖 4.12 所示。. 51.

(64) 圖 4.11、集集餘震花蓮分區(HWA)各測站震源加速度函數，經過相對極大值對齊、平均挑選過後的平均分區震源加速度時間函數圖 52.

(65) 圖 4.12、花蓮分區測站與兩事件分布圖 53.

(66) 2.. 台中分區與花蓮分區的測站分布差不多遠，距離餘震震央最近的測站為 TCU067，震央距 44.31 公里，S-P 波走時差 6.212 秒；距離餘震震央最遠的測站為 TCU064，震央距 72.64 公里，S-P 波走時差 9.503 秒。全部有觸發的測站，其走時差皆在 10.24 秒內。代入 Qp 三維模型後，各測站震央距、走時差，如圖 4.13 所示；各測站與震央之分布圖如圖 4.14 所示。. 54.

(67) 圖 4.13、集集餘震台中分區(TCU)各測站震源加速度函數，經過相對極大值對齊、平均挑選過後的平均分區震源加速度時間函數圖. 55.

(68) 圖 4.14、台中分區測站與兩事件分布圖 56.

(69) 3.. 台東分區是其次遠的測站，距離震央最近的測站為 TTN031，震央距 66.92 公里，S-P 波走時差約 8.45 秒，距離震央最遠的測站為 TTN050，震央距約 127.61 公里，S-P 波走時差 15.451 秒；26 個接收到資料的測站，有 7 個測站 S-P 波走時差在 10.24 秒以內。代入 Qp 三維模型後，各測站震央距、走時差，如圖 4.15 所示，測站分布與震央分布圖，如圖 4.16 所示。. 57.

(70) 圖 4.15、集集餘震台東分區(TTN)各測站震源加速度函數，經過相對極大值對齊、平均挑選過後的平均分區震源加速度時間函數圖. 58.

(71) 圖 4.16、台東分區測站與兩事件分布圖 59.

(72) 4.. 嘉南分區是這四個分區中距離餘震震央最遠的分區，最接近震央的測站為 CHY052，震央距 81.64 公里，S-P 波走時差 10.664 秒，距離震央最遠的是 CHY069，震央距 129.29 公里，S-P 走時差 16.077 秒；全部測站 S-P 波走時皆在 10.24 秒以上。代入 Qp 三維模型後，各測站震央距、走時差，如圖 4.17 所示；嘉南分區測站分布與震央關係圖，如圖 4.18 所示。. 60.

(73) 圖 4.17、集集餘震嘉南分區(CHY)各測站震源加速度函數，經過相對極大值對齊、平均挑選過後的平均分區震源加速度時間函數圖 61.

(74) 圖 4.18、嘉南分區測站與兩事件分布圖 62.

(75) 由以上四個分區之特性，我們另外使用同樣的方式解算集集餘震之震源加速度時間函數，與解算集集主震時所用的相同步驟，獲得各分區場址平均化之震源加速度時間函數。解算結果如圖 4.19 所示、去除台中分區，將剩下的三個分區放大繪圖所得如圖 4.20。. 圖 4.19、09220014 花蓮分區、台中分區、台東分區與嘉南分區之平均震源加速度時間函數之比較. 63.

(76) 圖 4.20、將集集餘震之花蓮分區、台東分區、嘉南分區之震源加速度函數放大進行比較. 本次地震規模約 6.2。若將四個分區拿來做比較，可以看出 S 波滲入測站的數量多寡，會影響解算時震源加速度時間函數時振幅相對之大小。但在圖 4.19 中可以發現，雖然台中分區相對於其他測站比較，所有測站 S-P 波走時差皆在 10.24 秒內，與花蓮分區測站群之距離差不多，但是解算出來的震源加速度時間函數，卻與其他分區之震源加速度時間函數大相逕庭；花蓮分區測站只有一個測站，其 S-P 波走時差在 10.24 秒以上，卻反而比較貼近台東分區與嘉南分區的震源加速度時間函數。造成台中分區震源加速度時間函數這樣的結果，各 64.

(77) 測站所使用來記錄之儀器型號，在各地為隨機分布，並無任何一特定型號之地震儀特別集中於某一分區內，故儀器效應可以忽略；最有可能的原因，是因為整個台中分區相對於其他分區比較，本身就是一個範圍較大的場址效應區，造成在函數振幅上過大的數值。此外，本餘震研究結果再次證明，越多震央距較遠的測站，所解算出來的震源加速度函數值會越小。另外，此次的集集餘震之震源加速度函數與集集主震之震源加速度函數相比，規模 7.6 之集集主震從花蓮、嘉南、台東分區所解算的震源加速度函數，從前 4 秒的 1gal，到 5 秒開始後不同分區峰值有 1~3gal 之變化；規模 6.2 之集集餘震之震源加速度函數中，花蓮分區震央距中位數 57.24 公里、台東分區震央距中位數 91.42 公里、嘉南分區震央距中位數 105.03 公里，在一開始就因不同分區震央距就已有不同振幅大小，之後再出現反向振幅；花蓮分區約 4~5gal、台東分區約 3~4gal、嘉南分區約 1~1.5gal 振幅變動。. 65.

(78) 第五章. 結論. 整理本研究所獲得之結論如下： 1.. 集集主震的震源加速度函數，前 3~4 秒內約在 0~1gal 之間有近乎水平，不到 1gal 以內的波動；4 秒之後開始有劇烈的曲線波動，有部分測站 S-P 走時差在 10.24 秒以內的花蓮與嘉南分區，振幅變動達到 3gal.之後，在 7 秒時達到最大峰值，之後振幅逐漸減弱，至 10.24 秒為止恢復成前 3~4 秒內的振幅值約 1gal.；嘉南分區震央距中位數大於花蓮分區震央距之中位數，在此次事件卻表現出幾近相同的振幅，表示嘉南分區之沖積層場址效應在此次主震有顯現出來；而所有測站 S-P 波走時差皆在 10.24 秒以上、震央距約在 76 公里以上的台東分區，約在 4~5 秒內開始變動，直到 7 秒達到最大峰值，但其峰值大約只有 1~2gal.左右，之後振幅同樣逐漸減弱，直至恢復成前 3~4 秒間的振幅值，因無 S 波或其他波介入之影響，台東分區之平均震源加速度函數振幅最小。台中分區測站因距離震央過近，在加上有隨後的表面波跟來，其震源加速度函數在前 4 秒內有微幅約 5gal 的擾動，5 秒後開始出現劇烈振幅變化，至 7~8 秒之間達到最大峰值，高峰值可到達 10gal，之後再逐漸減弱，卻因為有表面波的滲入，再形成另一個反方向的高峰值，大約在±10~20gal.之間波動。. 2.. 分析兩天後所發生的集集地震之餘震事件(09220014)，與主震集集事件做比較；集集餘震事件之震源加速度函數，在 0~5 秒時，呈現水平至微幅抬升的波動。含有 S-P 波走時差 10.24 秒以內的花蓮與台東分區，較近的花蓮分區可達 3~4gal、次之台東分區也 66.

(79) 有 2~3gal，全部測站走時差皆在 10.24 秒以上的嘉南分區，幾乎沒有波動；5 秒開始之後，開始產生另一波動，直至 7 秒達到最大峰值，其中花蓮分區達到 5gal、台東分區達到 2.5gal、嘉南分區只有 1~2gal 左右的峰值，7 秒之後開始減弱，至 10.24 秒時回復至前 3~4 秒內振幅。在此一提，花蓮分區與台中分區的測站 S-P 波走時差，無論是由平均或是由中位數來看，同樣有 8~9 秒走時差，而且由頭尾兩邊的振幅大小相近來看，兩個分區都只受到 S 波之滲入，並無表面波之干擾，但台中分區的震源加速度函數之振幅，與花蓮分區的相比仍然呈現大幅度的波動，峰值可達 12~13gal.，兩個分區的震源加速度函數呈現如此的差異，推測台中分區與其他三個分區相比，是場址效應較明顯分區，造成振幅高達±10gal 以上量級。兩個地震事件皆證明了震央距是主要影響解算過後震源加速度函數結果的最大因素，而且集集主震對於台中分區，還有表面波的滲入，造成解算過後較大的震源加速度函數。由於此次餘震在嘉南分區之沖積層效應並未像集集主震事件一樣，顯現在較小振幅上，未來若改由岩盤站或沖積層站來做台灣地區測站之分區分界，是否能看得出場址效應對於震源加速度函數的影響，有待進一步研究。. 67.

(80) 參考文獻 Akaike, H. (1971).Autoregressive model fitting for control, Annals of the Institute of Statistical. Mathematics, 23, 163-180. Akaike, H. (1979). A bayesian extension of the minimum AIC procedure of autoregressive model fitting, Biometrika, 66, 237-242. Allen, R.M., Kanamori, H.(2003) The Potential for Earthquake Early Warning in Southern California, Science, 300, 786-789. Azakawa, T. (2004a). Automatic Detection Procedure of Onset Time of Pand S-Phase in Seismic Records. Architectural Institude of Japan, 21246, 491-492. [Text in Japanese] Azakawa, T. (2004b). A Technique for Automatic Detection of Onset Time of P- and S-Phase in Strong Motion Records. 13th World Conference on Earthquake Engineering, Symposium held at Vancouver, B. C., Canada. Banks, R. J. (1975). Complex demodulation of geomagnetic data and the estimation of transfer function. Geophys. J. R. astr. Soc., 43, 87-101. Brigham, E. O. (1974). The fast fourier transform (2nd ed.). Englewood Cliffs, NJ: Prentice-Hall International. Brigham, E. O. （1985）。快速傅立葉變換（The Fast Fourier Transform）（黎文明譯）。台南市：復漢出版社。（原作 1974 年出版） Chen, K. J. (1998). S-Wave Attenuation Structure in the Taiwan Area and its Correlation to Seismicity. Terrestrial, Atmospheric and Oceanic Sciences, 9(1), 97-118. Chen, K. J., Yeh, Y. H., & Shyu, C. T. (1996). Qp Structure in the Taiwan Area and Its Correlation to Seismicity. Terrestrial, Atmospheric and Oceanic Sciences, 7(4), 409-429. Chen, K.J., Yeh, Y. H., & Tsai, Y. B. (1985). A study on complex demodulation of the geomagnetic total intensity data of Taiwan. Bull. Inst. Earth Sci., Academia Sinica, 5, 125-136. 68.

(81) Chen, K. J., Ho, Y. R., Chiu B., Wang, J.S., Sun, R., Lin, C. H. (2007). Correlation Between Time Change in Modulus of Short-Period Geomagnetic Variation and Seismicity in Taiwan. Terrestrial, Atmospheric and Oceanic Sciences, 18(3), 577-591. Chopra,. S., Kumar, D., & Rastogi, B. K. (2011). Attenuation of high frequency P and S waves in the Gujarat region, India. Pure and applied geophysics, 168(5), 797-813.. Geo SIG Presentation Seismic Signals and Sensors, Retrieved August 07, 2014, from http://www.geosig.com/files/GS_Presentation_Seismic_Signals_an d_Sensors.pps Hough, S. E., Anderson, J. G., Brune, J., Vernon III. F., Berger, J., Fletcher, J., Haar, L., Hanks, T., & Baker L. (1988). Attenuation near Anza, California. Bulletin of the Seismological Society of America, 78(2), 672-691. Kanamori, H. (1967). Spectrum of Short-period core phases in relation to the attenuation in the mantle. Journal of Geophysical Research, 72(8), 2181-2186. Maeda, N., (1985). A method for reading and checking phase times in auto-processing system of seismic wave data. Zisin Ser. 2, 38(3), 365-379. [Text in Japanese] Ma, K. F., Wang, J. H., Zhao, D. (1996). Three-Dimensional Seismic Velocity Structure of the Crust and Uppermost Mantle beneath Taiwan. Journal of Physics of the Earth, 44(2), 85-105. Rau, R. J.,Wu. Y. M. (1995). Tomographic imaging of lithospheric structures under Taiwan. Earth and Planetary Science Letters, 133(3-4), 517-532. Roecker, S. W., Yeh, Y. H., & Tsai, Y. B. (1987). Three-dimensional P and S wave velocity structures beneath Taiwan: deep structure beneath an arc-continent collision. Journal of Geophysical Research, 92(B10), 10547-10570. 69.

(82) Sleeman, R., Torild van Eck (1999). Robust automatic P-phase picking: an on-line implementation in the analysis of broadband seismogram recordings. Physics of the Earth and Planetary Interiors, 113, 265-275. WPSMap - Purchase-Build a Seismometer, Retrived Novenber 23, from http://www.wpsmap.com/menu/seismometer.php Waves, Seismometers, and Seismograms, Retrived August 07, 2017, from http://eqseis.geosc.psu.edu/~cammon/HTML/Classes/IntroQuakes/ Notes/seismometers.html Yoshida, M., Chiba, H., Hagiwara, H., Ogata, Y., and Mizoue, M. (1994). Summary of Real Time Processing Methods for Micro-earthquake Automatic Processing System (1981-1993). Bull. Earthq. Res. Inst., 69, 121-138. [Text in Japanese] 工程地震學研究室(2017)。工程地震學研究室。國立中正大學地球與環境科學系。上網日期：2017 年 9 月，檢自： http://www.eq.ccu.edu.tw/lab/lab520/introduction.html 中央氣象局地球物理資料管理系統。自由場強地動觀測網。上網日期： 2014 年 9 月，檢自： http://gdms.cwb.gov.tw/gdms-freefield.php 中央氣象局地震測報中心(2009)。中央氣象局地震測報中心 20 週年專刊. 中央氣象局地震測報中心. Retrieved from http://scweb.cwb.gov.tw/Twenty.aspx?ItemId=8&loc=tw 王仁聖(2011)。以時頻分析探討地動加速度譜之異常。國立臺灣師範大學地球科學系碩士論文。地震活動彙整，中央氣象局。上網日期：2014 年 9 月，檢自 http：//www.cwb.gov.tw/V7/earthquake/rtd_eq.htm 林孝維(2004)。利用傅氏振幅譜比法分析全台灣強震站的場址。國立中央大學地球物理研究所碩士論文。何英銳(2003)。地磁複變解調與地震活動之相關性：以台灣崙坪全磁力資料為例。國立臺灣師範大學地球科學研究所碩士論文。李錫堤、邱宏智、林柏伸、鄭錦桐(2006)。強震資料庫建置與維護及 70.

(83) 使用者平台開發(一)。國立中央大學財團法人中興工程顧問社技術報告。梁文宗(2016)。氣象局新一代地球物理管理資料系統規劃。交通部中央氣象局地震類委託研究計畫期末報告。陳光榮(1993)。台灣地區 Qp 值之空間分佈及其特性。國立台灣大學海洋研究所博士論文。陳光榮(2010)。台灣北部地區中大型地震震波時頻即時分析系統。交通部中央氣象局 101 年度地震類委託研究計劃。張議仁(2004)。大屯火山區三維 Vp、Vp/Vs 及 Qp 構造反演。國立中正大學地震研究所碩士論文。劉貞伶(2009)。以複變解調法分析竹子湖強震站最大地動加速度的異常放大現象。國立臺灣師範大學地球科學系碩士論文。. 71.