以頻域及時頻分析輔助孔內震測走時分析自動化

全文

(1)國立交通大學土木工程學系. 碩士論文. 以頻域及時頻分析輔助孔內震測走時分析自動化 Automation of Borehole Seismic Travel-time by Spectral and Time-Frequency Analysis. 研究生：胡成泓指導教授：林志平博士. 中華民國一百零一年七月.

(2) 以頻域及時頻分析輔助孔內震測走時分析自動化 Automation of Borehole Seismic Travel-time by Spectral and Time-Frequency Analysis 研究生：胡成泓指導教授：林志平. Student：Cheng-Hung Hu 博士. Advisor：Dr. Chih-Ping Lin. 國立交通大學土木工程學系碩士論文. A Thesis Submitted to Department of Civil Engineering College of Engineering National Chiao Tung University In Partial Fulfillment of the Requirements for the Degree of Master in Civil Engineering. July 2012 Hsinchu, Taiwan, Republic of China 中華民國一百零一年七月.

(3) 以頻域及時頻分析輔助孔內震測走時分析自動化學生：胡成泓. 指導教授：林志平博士國立交通大學土木工程學系. 中文摘要現有的孔內震測方法多採用以人工手動挑取初達波的走時分析方法，然而手動挑波結果往往因分析人員的經驗、訊號品質優劣而有差異，當資料量多繁雜時，相當費時費工；而當震波訊號的訊雜比較低時，判讀上的困難將使其難保正確性；此外，手動挑取之初達時間往往只在時間域上判釋，並無法得知其所得速度之頻率範圍。雖然陸續一直有學者進行初達波自動化選取之相關研究，但所提出之方法基本上皆是基於時間域之分析，且常常是有條件限制之適用性，對於訊號品質較差的資料易有判讀錯誤的情形。故本研究從訊號分析的角度，利用頻率域及時間－頻率域之訊號分析方法，透過有效的濾波，一方面由不同場域觀察震波訊號之特性，另一方面藉以訂定出半自動挑波法，並透過現地資料之應用探討不同孔內震測技術之適用性。本研究在頻率域下採用頻譜分析法，透過簡化之均質均向土層假設，利用各頻率下之相位角差斜率以獲得相位速度；而在時間－頻率域上採用多重濾波法及小波轉換繪製時頻譜，並提出在時頻譜上利用振幅起伏起始點之自動挑取方法獲得頻帶波相平均初達時間，進而計算該頻帶下之初達波速。經現地資料測試，頻率域之簡化頻譜分析法較適用於震源穩定且有兩個接受器之孔內懸盪式震測法，此分析可同時提供剪力波速並可有額外之頻率訊息，可提供工程師據以進行進一步之動態分析，雖然此方法具有自動化分析之高度潛力，但其對於訊雜比低之訊號結果誤差較大，需進一步進行分析探討；而採用時頻分析法所訂定出之半自動挑法，無論採用多重 I.

(4) 濾波法或小波轉換法於現地懸盪式震測及下孔震測之資料驗證，皆有良好表現，整體波速趨勢與手動挑波結果一致，且相對穩定，並具有良好的抗雜訊能力，然而此時頻譜半自動挑波方法受限於時頻譜分析技術之時頻域間的解析度限制，針對不同震源類型的孔內震測試驗須調整分析參數達到最佳化，以便半自動挑波法能正確判釋波相初達時間，本研究根據現地試驗資料針對孔內懸盪式震測法以及下孔式震測法分別提出建議之使用參數。綜合成果顯示，頻譜分析法較適合應用於震源穩定且高頻的懸盪式震測法，惟對訊雜比低之訊號需後續進一步探討；而時頻分析法對於各項孔內震測皆可有良好應用成效，且有良好之抗雜訊能力，可廣為推廣應用。. 關鍵字：孔內震測方法、走時分析、初達波、自動挑波法、時頻分析. II.

(5) Automation of Borehole Seismic Travel-time by Spectral and Time-Frequency Analysis. Student：Cheng-Hung Hu. Advisor：Dr. Chih-Ping Lin. Department of Civil Engineering National Chiao Tung University. Abstract The data reduction of existing borehole seismic methods are currently based on travel-time analysis in which first-arrival times of the seismogram are manually picked. However, manual picking is not efficient and depends heavily on the data quality and analyst’s experience. Besides the difficulty of picking first arrivals in the time domain, no information is obtained regarding the effective frequency of the analysis. Several automatic picking methods were proposed to increase the efficiency of data processing. But most of them are based on time-domain analysis and work well only on some specific conditions. Data reduction methods based on the frequency domain and time-frequency domain were proposed in the study. The aim of the study was to observe the behavior of the borehole seismic waves from different domains (i.e. frequency and time-frequency domain) and propose a semi-automatic procedure for more objective first-arrival time picking. Field data from PS-Log testing and down-hole testing were used to evaluate the feasibility of the proposed procedure, in particular for the shear wave velocity. In the frequency domain, the phase velocity can be calculated from the phase shift of two receiving stations in each frequency. To avoid the difficulty of unwrapping phase angles in lack of low frequency components, the formation around the borehole is assumed homogeneous and isotropic such that the phase velocity can be calculated from the slope of the phase difference versus frequency at certain frequency range. A better understanding of the seismic III.

(6) signal can be achieved by time-frequency analysis, such as multi-filter analysis or wavelet transform. In the time-frequency domain, a semi-automatic procedure was proposed to determine the first arrival times from the envelope of the time-frequency spectrum at specific frequency range. Evaluation using field data shows that the spectral analysis method is suitable for PS-Log testing with stable source and two simultaneous recording channels. The method is easy to be automated and yields Vs and the corresponding effective frequency. However, more study should be conducted to overcome. occasional. failure. in. low. S/N. data.. The. semi-automatic. time-frequency spectrum analysis, either based on multi-filter analysis or Morlet wavelet transformation, works well both in PS-Log testing and down-hole testing. The picked first arrival times from the proposed method agree well with the manually-picked results in high quality data. The procedure is stable and yields reasonable results in data with low S/N. But it should be noted that the method inherit the tradeoff between frequency and time resolution from the time-domain analysis. The parameters used in the time-frequency analysis should be optimized for different source types. The parameters for PS-Log testing and down-hole testing were suggested respectively in the study and verified by the field data.. Keywords : borehole seismic method, travel time analysis, first-arrival time, automatic first-arrival picking, time-frequency analysis. IV.

(7) 致謝本文得以順利完成，首要感謝指導教授林志平博士於研究所兩年來的悉心教誨，給予學生在論文研究以及專業知識上的諸多指導，使學生得以突破研究瓶頸並獲得學識上的精進。此外，於求學態度與待人處世方面亦給予學生許多啟發，尤其對於邏輯思考及簡報表達的能力，更是大大提升。很慶幸能夠進入林老師門下，在此致上最誠摯的謝意。求學期間，承蒙廖志中教授、潘以文教授、單信瑜教授、方永壽教授、黃安斌教授於學識上的指導，在此致上萬分謝意。而口試期間，良受馮正一教授、余騰鐸教授、葛宇甯教授、張宗盛博士給予本論文許多寶貴的意見與建議，使得本論文更加完善，在此表示由衷感謝。研究期間承蒙俊宏、志忠、柏林、瑋晉等學長兩年的照顧與協助。尤其感謝俊宏學長總是幫我解決很多研究上的問題，指導我的簡報與論文。感謝王凱、士輝、逸倫學長帶我度過充實歡樂的碩一生涯。感謝同門志昇兩年來的共患難，以及其他同甘共苦的大地組同學們。也感謝友誠、昶笙、瀅郁等學弟妹的幫忙，讓我的碩二生涯能專心致力於論文研究。最後，感謝家人們對我的支持與栽培，這二十幾年來如果沒有您們包容與教誨便沒有今天的我，謝謝您們。此外特別感謝五年來一直陪在我身邊的女友郁珊，妳的陪伴與鼓勵是我一直以來努力的泉源。僅以此文獻給我最親愛的家人與女友郁珊。. V.

(8) 目錄中文摘要 ................................................................................................................ I Abstract .............................................................................................................. III 致謝 ...................................................................................................................... V 目錄 ..................................................................................................................... VI 圖目錄 ................................................................................................................. IX 表目錄 ............................................................................................................... XV 第一章前言 ....................................................................................................... 1 1.1 研究背景與動機 .................................................................................... 1 1.2 研究目的 ................................................................................................ 5 1.3 研究流程與論文架構 ............................................................................ 6 第二章文獻回顧 ............................................................................................... 8 2.1 孔內震測方法 ........................................................................................ 8 2.1.1 下孔式震測法 ............................................................................. 9 2.1.2 跨孔式震測法 ........................................................................... 16 2.1.3 懸盪式震測法 ........................................................................... 19 2.2 走時分析方法與評析 .......................................................................... 22 2.3 自動挑波方法與評析 .......................................................................... 25 2.3.1 STA/LTA 方法 ........................................................................... 25 2.3.2 AIC 方法 ................................................................................... 27 2.3.3 類神經網路 ............................................................................... 29 2.3.4 能量比值法 ............................................................................... 32 2.4 訊號分析方法 ...................................................................................... 35 2.4.1 時間域濾波分析 ....................................................................... 35 2.4.2 頻率域濾波分析 ....................................................................... 38 VI.

(9) 2.4.3 時頻譜濾波分析 ....................................................................... 40 第三章研究方法 ............................................................................................. 52 3.1 研究方法流程說明 .............................................................................. 53 3.2 速度分析方法 ...................................................................................... 54 3.2.1 時間域速度分析法 ................................................................... 54 3.2.2 頻率域相位速度分析法........................................................... 55 3.2.3 時間－頻率域速度分析法....................................................... 60 3.3 振幅起伏起始點挑取自動化分析法 .................................................. 62 3.4 試驗資料 .............................................................................................. 64 3.4.1 半自動化挑波程序訂定測試資料........................................... 64 3.4.2 半自動化挑波程序驗證用資料............................................... 69 第四章試驗結果與討論 ................................................................................. 72 4.1 半自動化挑波流程訂定 ...................................................................... 72 4.1.1 頻率域分析法 ........................................................................... 72 4.1.2 時間－頻率域分析法 ............................................................... 79 4.1.3 時窗函數前處理影響 ............................................................. 101 4.1.4 分析方法比較 ......................................................................... 104 4.1.5 半自動化挑波分析法流程..................................................... 107 4.2 半自動化挑波法驗證 ........................................................................ 107 4.2.1 懸盪式震測資料 ..................................................................... 107 4.2.2 下孔式震測資料 ..................................................................... 119 4.2.3 綜合比較 ................................................................................. 125 第五章結論與建議 ....................................................................................... 126 5.1 結論 .................................................................................................... 126 5.2 建議 .................................................................................................... 126 VII.

(10) 參考文獻 ......................................................................................................... 128 附錄(一) .......................................................................................................... 132. VIII.

(11) 圖目錄圖 1.1 結合 SPT 與上孔震測，以 SPT 的敲擊作為震源，在孔旁設置受波器。 ................................................................................................................... 2 圖 1.2 下孔震測 Test-01 資料及走時軌跡 ......................................................... 3 圖 1.3 下孔震測 Test-02 資料及走時軌跡 ......................................................... 4 圖 1.4 Test-01、Test-02 走時軌跡比較 .............................................................. 4 圖 1.5 研究流程圖 ............................................................................................... 6 圖 2.1 孔內震測法示意圖(Seismic Shear Wave Studies) .................................. 9 圖 2.2 非破壞性地表震測法示意圖(Seismic Shear Wave Studies) .................. 9 圖 2.3 下孔式震測法配置圖(D7400-08, 2008) ................................................ 10 圖 2.4 直接法 (direct method)示意圖(Kim et al., 2004) ................................. 13 圖 2.5 間距法 (interval method)示意圖(Kim et al., 2004) .............................. 13 圖 2.6 修正間距法 (modified-interval method) 示意圖(Kim et al., 2004) .... 14 圖 2.7 逆推法 (inversion method) (Kim et al., 2004) ...................................... 14 圖 2.8 下孔震測利用左敲、右敲輔助凸顯剪力波的對稱性(D7400-08, 2008) ................................................................................................................. 15 圖 2.9 跨孔式震測法配置圖(D4428/D4428M-07, 2004) ................................ 17 圖 2.10 利用跨孔式震測法施作斷層剖面速度影像方法(Angioni et al., 2003) ................................................................................................................. 18 圖 2.11 懸盪式波速量測系統配置圖(OYO corporation) ................................ 20 圖 2.12 震源原理示意圖(OYO suspension PS logger Manual) ....................... 20 圖 2.13 懸盪式震測法訊號波形(國家地震工程研究中心) ............................ 21 圖 2.14 不同選點方式的正規化波速比較(Stokoe and Hoar, 1978) ............... 23 圖 2.15 P 波中 SV 波的近場效應(Sánchez-Salinero, 1987) .......................... 24 圖 2.16 SV 波中 P 波的近場效應(Sánchez-Salinero, 1987) .......................... 24 IX.

(12) 圖 2.17 STA/LTA 方法示意圖(Han et al., 2010) .............................................. 26 圖 2.18 AIC 方法判讀兩種波相到時示意圖(Zhang et al., 2003) ................... 27 圖 2.19 長、中、短三個時窗判斷器，並以一並聯判斷器取其交集判斷(Zhao and Takano, 1999) ................................................................................... 30 圖 2.20 長、中、短三個時窗判斷器示意圖(Zhao and Takano, 1999) ......... 31 圖 2.21 固定時窗(罗进华 et al., 2008) ............................................................ 33 圖 2.22 滑動時窗示意圖 (徐钰 et al., 2012) .................................................. 33 圖 2.23 修正能量比值法(MER)示意圖 ........................................................... 34 圖 2.24 闡述非週期訊號取樣時造成遺漏的原因(CRYSTAL instruments)... 36 圖 2.25 由於資料截斷造成頻譜能量遺漏(测试技术实验指导) .................... 36 圖 2.26 SASW 頻譜分析(相位差-頻率圖)(林俊宏, 2005) .............................. 39 圖 2.27 多重濾波法案例。(賴雅娟, 2009) ...................................................... 43 圖 2.28 傅立葉傳換、短時傅立葉轉換、小波轉換解析度比較(Chun-Lin, 2010) ................................................................................................................. 45 圖 2.29 小波函數的尺度變化(李立愷, 2005) .................................................. 47 圖 2.30 小波函數的平移變化(李立愷, 2005) .................................................. 47 圖 2.31 各種母小波函數；(a)Haar (b)Mexican hat (c)Meyer (d)Morlet (e)Dauvechies(db6) ............................................................. 49 圖 2.32 Morlet 小波的實部(實線)、虛部(虛線)及其傅立葉轉換頻譜 ......... 50 圖 3.1 本研究之試驗階段流程 ......................................................................... 53 圖 3.2 時間域訊號波形及手動挑取初波達時間點......................................... 55 圖 3.3 兩時間訊號之相位差(Joh, 1996) .......................................................... 56 圖 3.4 懸盪式震測法及表面波譜法之設置比較............................................. 56 圖 3.5 頻譜分析於懸盪式震測資料之應用步驟............................................. 59 圖 3.6 訊號經時頻分析方法轉換成時頻譜 ..................................................... 60 圖 3.7 訊號各頻率下之波形及其包絡線 ......................................................... 61 X.

(13) 圖 3.8 振幅起伏起始點挑取自動化方法 ......................................................... 63 圖 3.9 受波器與震源 ......................................................................................... 65 圖 3.10 控制監測裝置 ....................................................................................... 65 圖 3.11 絞盤 ....................................................................................................... 65 圖 3.12 #008 訊號波形及初波達時間點 .......................................................... 66 圖 3.13 #009 訊號波形及初波達時間點 .......................................................... 67 圖 3.14 #010 訊號波形及初波達時間點 .......................................................... 67 圖 3.15 #082 訊號波形及初波達時間點 .......................................................... 68 圖 3.16 #083 訊號波形及初波達時間點 .......................................................... 68 圖 3.17 #084 訊號波形及初波達時間點 .......................................................... 69 圖 3.18 交通大學博愛校區下孔震測鑽孔位置............................................... 70 圖 3.19 下孔震測施測設備示意圖(鄒和翰, 2005) .......................................... 71 圖 3.20 震測儀與孔內受波器 ........................................................................... 71 圖 4.1 #008 於頻譜分析法所得能量分佈與相位角差頻譜............................ 73 圖 4.2 #009 於頻譜分析法所得能量分佈與相位角差頻譜............................ 74 圖 4.3 #010 於頻譜分析法所得能量分佈與相位角差頻譜............................ 75 圖 4.4 #082 於頻譜分析法所得能量分佈與相位角差頻譜............................ 76 圖 4.5 #083 於頻譜分析法所得能量分佈與相位角差頻譜............................ 77 圖 4.6 #084 於頻譜分析法所得能量分佈與相位角差頻譜............................ 78 圖 4.7 多重濾波法參數設置比較 ..................................................................... 80 圖 4.8 多重濾波法於測試資料#008 之應用(受波器 H1 資料) ...................... 81 圖 4.9 多重濾波法於測試資料#008 之應用(受波器 H2 資料) ...................... 82 圖 4.10 多重濾波法於測試資料#008 之應用結果.......................................... 82 圖 4.11 多重濾波法於測試資料#009 之應用(受波器 H1 資料) .................... 83 圖 4.12 多重濾波法於測試資料#009 之應用(受波器 H2 資料) .................... 83 XI.

(14) 圖 4.13 多重濾波法於測試資料#009 之應用結果.......................................... 84 圖 4.14 多重濾波法於測試資料#010 之應用(受波器 H1 資料) .................... 84 圖 4.15 多重濾波法於測試資料#010 之應用(受波器 H2 資料) .................... 85 圖 4.16 多重濾波法於測試資料#010 之應用結果.......................................... 85 圖 4.17 多重濾波法於測試資料#082 之應用(受波器 H1 資料) .................... 86 圖 4.18 多重濾波法於測試資料#082 之應用(受波器 H2 資料) .................... 86 圖 4.19 多重濾波法於測試資料#082 之應用結果.......................................... 87 圖 4.20 多重濾波法於測試資料#083 之應用(受波器 H1 資料) .................... 87 圖 4.21 多重濾波法於測試資料#083 之應用(受波器 H2 資料) .................... 88 圖 4.22 多重濾波法於測試資料#083 之應用結果.......................................... 88 圖 4.23 多重濾波法於測試資料#084 之應用(受波器 H1 資料) .................... 89 圖 4.24 多重濾波法於測試資料#084 之應用(受波器 H2 資料) .................... 89 圖 4.25 多重濾波法於測試資料#084 之應用結果.......................................... 90 圖 4.26 小波轉換法於測試資料#008 之應用(受波器 H1 資料) .................... 92 圖 4.27 小波轉換法於測試資料#008 之應用(受波器 H2 資料) .................... 92 圖 4.28 小波轉換法於測試資料#008 之應用結果.......................................... 93 圖 4.29 小波轉換法於測試資料#009 之應用(受波器 H1 資料) .................... 93 圖 4.30 小波轉換法於測試資料#009 之應用(受波器 H2 資料) .................... 94 圖 4.31 小波轉換法於測試資料#009 之應用結果.......................................... 94 圖 4.32 小波轉換法於測試資料#010 之應用(受波器 H1 資料) .................... 95 圖 4.33 小波轉換法於測試資料#010 之應用(受波器 H2 資料) .................... 95 圖 4.34 小波轉換法於測試資料#010 之應用結果.......................................... 96 圖 4.35 小波轉換法於測試資料#082 之應用(受波器 H1 資料) .................... 96 圖 4.36 小波轉換法於測試資料#082 之應用(受波器 H2 資料) .................... 97 圖 4.37 小波轉換法於測試資料#082 之應用結果.......................................... 97 XII.

(15) 圖 4.38 小波轉換法於測試資料#083 之應用(受波器 H1 資料) .................... 98 圖 4.39 小波轉換法於測試資料#083 之應用(受波器 H2 資料) .................... 98 圖 4.40 小波轉換法於測試資料#083 之應用結果.......................................... 99 圖 4.41 小波轉換法於測試資料#084 之應用(受波器 H1 資料) .................... 99 圖 4.42 小波轉換法於測試資料#084 之應用(受波器 H2 資料) .................. 100 圖 4.43 小波轉換法於測試資料#084 之應用結果........................................ 100 圖 4.44 Tukey Window 在時間域及頻率域之波形 ....................................... 102 圖 4.45 套用時窗函數至訊號做前處理 ......................................................... 103 圖 4.46 時窗函數前處理影響於頻譜分析法之比較..................................... 103 圖 4.47 時窗函數前處理影響於時頻分析法之比較..................................... 104 圖 4.48 訊雜比高之測試資料於各分析法波速比較圖................................. 105 圖 4.49 訊雜比低之測試資料於各分析法波速比較表................................. 106 圖 4.50 半自動挑波法以多重濾波法作時頻譜，應用於懸盪式震測資料之走時曲線比較 ........................................................................................... 108 圖 4.51 (a)半自動挑波法以多重濾波法作時頻譜，應用於懸盪式震測資料之剪力波速剖面比較; (b)頻帶波速剖面對應之頻率範圍 .................... 109 圖 4.52 半自動挑波法以小波轉換作時頻譜，應用於懸盪式震測資料之走時曲線比較 ............................................................................................... 110 圖 4.53 (a)半自動挑波法以小波轉換作時頻譜，應用於懸盪式震測資料之剪力波速剖面比較; (b)頻帶波速剖面對應之頻率範圍 ........................ 111 圖 4.54 懸盪式震測資料#069 波形及手動挑取剪力波初達時間 ............... 112 圖 4.55 懸盪式震測資料#069 套用時窗函數及其時頻譜 ........................... 113 圖 4.56 懸盪式震測資料#069 半自動挑取頻帶下剪力波初達時間 ........... 113 圖 4.57 懸盪式震測資料#070 波形及手動挑取剪力波初達時間 ............... 114 圖 4.58 懸盪式震測資料#070 套用時窗函數及其時頻譜 ........................... 114 XIII.

(16) 圖 4.59 懸盪式震測資料#070 半自動挑取頻帶下剪力波初達時間 ........... 115 圖 4.60 懸盪式震測資料#071 波形及手動挑取剪力波初達時間 ............... 115 圖 4.61 懸盪式震測資料#071 套用時窗函數及其時頻譜 ........................... 116 圖 4.62 懸盪式震測資料#071 半自動挑取頻帶下剪力波初達時間 ........... 116 圖 4.63 懸盪式震測資料#151 波形及手動挑取剪力波初達時間 ............... 117 圖 4.64 懸盪式震測資料#151 套用時窗函數及其時頻譜 ........................... 118 圖 4.65 懸盪式震測資料#151 半自動挑取頻帶下剪力波初達時間 ........... 118 圖 4.66 下孔式震測資料於手動挑波法及半自動挑波法之走時曲線比較(以多重濾波法作時頻譜) ......................................................................... 119 圖 4.67 下孔式震測資料於手動挑波法及半自動挑波法之剪力波速剖面(以多重濾波法作時頻譜) ......................................................................... 120 圖 4.68 下孔式震測驗證資料於半自動挑波法之頻帶波速剖面對應之頻率範圍(以多重濾波法作時頻譜) ................................................................ 120 圖 4.69 下孔式震測資料於手動挑波法及半自動挑波法之走時曲線比較(以小波轉換作時頻譜) ............................................................................. 121 圖 4.70 下孔式震測資料於手動挑波法及半自動挑波法之剪力波速剖面(以小波轉換作時頻譜) ............................................................................. 121 圖 4.71 下孔式震測驗證資料於半自動挑波法之頻帶波速剖面對應之頻率範圍(以小波轉換作時頻譜) .................................................................... 122 圖 4.72 小波參數 WinLen 於下孔式震測資料之時頻譜比較...................... 123 圖 4.73 小波參數 WinLen 於下孔式震測資料之半自動挑波比較.............. 124 圖 4.74 因小波參數 WinLen 異同造成走時曲線偏移.................................. 124. XIV.

(17) 表目錄表 2.1 小波轉換三大類比較 ............................................................................. 46 表 3.1 半自動化挑波程序訂定測試資料於時間域初達波分析剪力波速資料 ................................................................................................................. 69 表 4.1 測試資料於頻率域頻譜分析相位波速資料......................................... 79 表 4.2 測試資料於多重濾波法分析頻帶初達波速資料................................. 90 表 4.3 測試資料於小波轉換法分析頻帶初達波速資料............................... 101 表 4.4 訊雜比高之測試資料於各分析法波速比較表................................... 105 表 4.5 訊雜比低之測試資料於各分析法波速比較表................................... 106 表 4.6 波速特異資料之手動挑波與半自動挑波比較表............................... 117. XV.

(18) 第一章前言 1.1 研究背景與動機土層動態參數在地震工程中扮演極重要的角色，其用途包括地震危害度分析、土壤液化潛能分析、邊坡與土壩之動態穩定分析、及土壤結構動態之互制行為之研究、震動基礎設計分析等有關土壤動態行為方面的研究ㄧ直是被受矚目，固此求得正確土層動態參數，在大地工程動態分析與設計評估上為很重要的ㄧ環。為了求得土層剪力波速，現地震波量測方式大致可分為孔內震測及非破壞性地表震測兩類。孔內震測(Borehole Seismics)是利用機械方法或技術入侵材料之中，探測材料的各種特性，其方法包括上孔法、下孔法、跨孔法及懸盪式波速量測法等。非破壞性地表震測(Non- Destructive Surface Seismics)是利用物理、機械方法(或技術)，在不破壞材料的情況下，檢測材料是否產生缺陷，或是探測材料物理、機械性質等檢測技術，其方法包括折射震測法、反射震測法、表面波震測法等。另一方面，標準貫入試驗 N 值(Standard Penetration Test N Value, SPT-N)在大地工程之地質鑽探中，係用於研判地層之軟弱或緊密程度，可說是應用最廣、資料最豐富且最為經濟之調查手段之一。自(Terzaghi and Peck, 1948) 提出 N 值與砂性土壤之相對密度、容許承載力；N 值與黏土強度、容許承載力等間之關係以來，已累積相當多的資料可供大地工程設計之參考。(林士誠, 1999) 既然土層剪力波速及標準貫入試驗 N 值對於大地工程都是如此重要的一個環節，若能結合兩種試驗方法，不論在經濟效益、時間考量或是工程的需求上都有相當程度的提升。而在土層剪力波速之獲得上，由於孔內震測法在理論及計算上遠較無鑽孔探測法簡單，且可與標準貫入試驗同時施作，如圖 1.1 所示，以上孔震測法的配置方式，可以以標準貫入試驗之夯垂敲擊作為震源，配合標準貫入試驗在現地預定深度的鑽孔進行檢測，同 1.

(19) 時獲取各深度的 SPT-N 值及剪力波速。. 圖 1.1 結合 SPT 與上孔震測，以 SPT 的敲擊作為震源，在孔旁設置受波器。. 然而現有的孔內震測法多採用以人工手動挑取初達波的走時分析方法，相當費時費工、當震波訊號的訊雜比較低時，也有判讀的困難，故若想結合標準貫入試驗及孔內震測，解決上述問題才是關鍵。舉例來說，圖 1.2、圖 1.3 為在交通大學博愛校區一鑽孔所施測之下孔震測資料，並在同一鑽孔施測兩次，命名為 Test-01、Test-02，很明顯可以看出 Test-02 資料淺層部分由於高頻雜訊影響大，降低訊號品質，難以判斷出剪力波相之初達時間，尤其以圖左之 H1 方向分量最為明顯。圖 1.4 為兩次試驗之走時軌跡比較圖，明顯顯示出淺層軌跡之偏差確實較大。同理，若在其他試驗測得如 Test-02 資料資料品質較差的狀況，又沒有品質良好之對照組，便難以用人工手動挑波方法判斷出一致的初達時間。 2.

(20) 故發展自動化挑波技術對於地震工程及震測方法等相關研究中是十分重要的議題，雖然陸續一直有學者 (Allen, 1978; Coppens, 1985; Han et al., 2010; Leonard, 2000; Zhang et al., 2003)在進行相關研究，但所提出的方法常常是有條件的適用性限制，對於訊雜比較低的資料也易有分析錯誤的情形。除此之外，波速常有受阻尼衰變的現象，若只是單純針對時間域挑取初達波，對於挑取的波形只能知道其時間域的訊息，並無法其所表示之頻率為何，物理意義含量較少。所以本研究希望從訊號分析的角度，透過時間域、頻率域、及時頻分析進行有效的濾波，一方面觀察震波訊號在這些分析狀況下的情形，一方面探討其作為自動化挑波方法之可行性。. Left Hit +Right Hit (H2). 2. 2. 4. 4. 6. 6. 8. 8. Depth (m). Depth (m). Left Hit +Right Hit (H1). 10. 10. 12. 12. 14. 14. 16. 16. 18. 18. 0. 0.01. 0.02. 0.03. 0.04. 0.05. 0.06. 0.07. 0.08. 0.09. 0.1. 0. Time (s). 0.01. 0.02. 0.03. 0.04. 0.05 Time (s). 圖 1.2 下孔震測 Test-01 資料及走時軌跡. 3. 0.06. 0.07. 0.08. 0.09. 0.1.

(21) Left Hit +Right Hit (H2). 2. 2. 4. 4. 6. 6. 8. 8. Depth (m). Depth (m). Left Hit +Right Hit (H1). 10. 10. 12. 12. 14. 14. 16. 16. 18. 18. 0. 0.01. 0.02. 0.03. 0.04. 0.05. 0.06. 0.07. 0.08. 0.09. 0.1. 0. 0.01. 0.02. 0.03. 0.04. Time (s). 0.05. 0.06. 0.07. 0.08. 0.09. 0.1. Time (s). 圖 1.3 下孔震測 Test-02 資料及走時軌跡 H1. H2 2 Test-01. Test-01. Test-02. Test-02. 4. 4. 6. 6. 8. 8. Depth (m). Depth (m). 2. 10. 10. 12. 12. 14. 14. 16. 16. 18 0.015. 0.02. 0.025. 0.03. 0.035. 0.04. 18 0.015. 0.045. Time (s). 0.02. 0.025. 0.03 Time (s). 圖 1.4 Test-01、Test-02 走時軌跡比較 4. 0.035. 0.04.

(22) 1.2 研究目的 1. 透過時間域、頻率域與時頻域觀察孔內震測波形之行為，並據以訂定自動化/半自動化波速分析流程。 2. 比較自動化/半自動化波速分析流程於不同孔內震測資料的適用性差異。. 5.

(23) 1.3 研究流程與論文架構研究動機與目的. 準備階段. 測試階段. 文獻回顧. . . 震測方法. 自動化方法. 訊號分析. S/N高之資料 Time Domain. Frequency Domain. 手動挑波. Window. Window. Window. SA. MFT. WT. Vph. Vmf. Vwt. Vs. Time-Frequency Domain. 整合比較建立規範套用S/N低之資料. . 應用階段. . 評量階段. . NG. OK. PS-logger. Downhole Test 自動化分析. 研究結果分析結論與建議. 圖 1.5 研究流程圖. 6.

(24) 圖 1.5 為本研究之流程，分為五個階段：研究動機與目的、準備階段、試驗階段、應用階段及評量階段。這五大階段分別對應於本論文一至五章，論文架構如下：第一章為本研究之背景、動機與目的，還有簡單的研究方法與流程介紹。第二章即是準備階段，為相關研究的文獻回顧部分： (1) 介紹各種初達波分析式震測方法之原理及分析概念。 (2) 回顧現有初達波自動化挑取方法及其優劣性。 (3) 回顧時域、頻域及時—頻域常用的訊號分析方法。第三章闡述本研究之方法，分為三大部分—分別說明以時域、頻域及時頻域定義計算波速的方法及步驟，並說明半自動化挑取程序，以及試驗規劃說明。第四章則為試驗階段，包含了圖 1.5 的測試階段及應用階段： (1) 測試階段分別挑選出三組訊雜比較高及訊雜比較低之懸盪式震測資料作為初步測試，並從第三章的各個分析方法中規劃出半自動化挑波程序。 (2) 應用階段則是為了驗證測試階段所規劃出的半自動化挑波程序，而將其分別應用在懸盪式及下孔式震測資料。第五章為最後的評量階段，就本研究之試驗提出適用性評估及結果討論，並提出自動化分析之可行性及相關建議。. 7.

(25) 第二章文獻回顧本章是研究相關的文獻回顧，內容包含從現有的初達波分析式震測方法介紹，及目前採用人工手動挑波的走時分析方式，到各學者針對這種初達波分析提出的自動化方法，以及本研究在第三章希望利用其來輔助探討的各種訊號方式。 2.1 孔內震測方法現今量測地層剪力波速的方法主要以利用震波之性質及原理進行試驗為主，故通稱為震測法；並可被分為破壞性之孔內震測法及非破壞性之地表震測法。孔內震測法根據震源與受波器相對位置的不同，可區分為跨孔式震測法(crosshole seismic)、上孔式震測(uphole seismic)、下孔式震測法 (downhole seismic)以及懸盪式震測法(suspension PS logging)等，如圖 2.1。有鑽孔之震測法單價較高，但可以量測到比較精確之剪力波速度剖面。而非破壞性方法則如折射波震測法(refraction seismic)、反射波震測法 (reflection seismic)、表面波震測法(surface seismic)等，如圖 2.2。然而上述的震測方法，目前除了表面波震測法採用頻率域分析以外，其餘皆採用初達波分析，藉由手動挑選震波的初達時間 (first arrival time) 來進行波速的計算。本研究主要針對鑽孔內震測法作探討，故以下分別對下孔式震測法、跨孔式震測法、以及懸盪式震測法做簡單的概念性介紹。. 8.

(26) (a) 跨孔式震測法. (b) 下孔式震測法. (c) 懸盪式震測法. 圖 2.1 孔內震測法示意圖(Seismic Shear Wave Studies). (a) 折射波震測法. (b) 反射波震測法. (c) 表面波震測法. 圖 2.2 非破壞性地表震測法示意圖(Seismic Shear Wave Studies) 2.1.1 下孔式震測法下孔式震測法所需之設備包括震源、孔內受波器、震測儀，將孔內受波器置於不同深度位置並固定於孔壁，由震源在鑽孔上方激發一震波，並啟動震測儀紀錄受波器之震波訊號，由各深度受波器之震波訊號可分析震波由震源孔傳達至受波器孔之走時，據以決定地層之波速變化，如圖 2.3 所示。除了震源至受波器之直達波，若地層存在阻抗不連續面，受波器亦可記錄其反射波，由反射訊號可決定阻抗不連續面之位置。震源位置通常距離孔口約 3m，以降低管波之影響。. 9.

(27) 圖 2.3 下孔式震測法配置圖(D7400-08, 2008) 至於波傳距離與波傳時間的決定上，由於分析方法的不同，在時間與距離的取捨上也會有所差異，以下是下孔式震測法常用的幾種分析方法(Kim et al., 2004)； (1) 直接法 (direct method) 以直接法計算剪力波速為目前應用最普遍的方式，傳統上乃先將波傳時間(t)與波傳距離(R)的關係轉換成波傳時間(tc)與土層深度(D)的關係(如圖 2.4 所示)關係式如下： 𝑡 𝑅 ∆𝐷 𝑉𝑑 = ∆𝑡𝑐. (2-1). 𝑡𝑐 = 𝐷. (2-2) 10.

(28) D：受波器深度； R：震源與受波器距離； tc：直接傳遞之波傳時間 t：量得之波傳時間 𝑉𝑑 為各深度區間平均波速由深度與波傳時間關係之走時曲線圖，依其變化取線性迴歸得一相連折線，求各直線段斜率之倒數，即為該相對土層之剪力波速值，如(2-2) 式。 (2) 間距法 (interval method) 剪力波速值係以鑽孔中兩組受波器之波傳距離差與波傳時間相除而得(如圖 2.5 所示) ，關係式如下： 𝑉=. 𝑅2 − 𝑅1 𝑡2 −𝑡1. (2-3). 𝑅1 ：震源與上受波器之距離； 𝑅2 ：震源與下受波器之距離； 𝑡1 、 𝑡2 ：震源至上、下受波器之波傳時間由於是兩組受波器接收同一震源產生之震波，所以精度上較其他兩法要高，但土層變化時，將受波器間距與兩受波器間連繫桿件型式之影響。 (3) 擬間距法 (pseudo-interval method) 關係式同間距法，但僅有一組受波器，R1、R2 與𝑡1 、 𝑡2 分別為不同次敲擊時之深度與波傳時間，以此法可得較多剖面之剪力波速值，但相較於間距法由於前後兩次敲擊能量未必相同，因此以此法得土層剪力波速顯得較具爭議。(Patel, 1981) (4) 修正間距法 (modified-interval method) 施作方法同間距法，但假設現地為水平層狀地層，故震波在各分層會有不同的波傳速度(圖 2.6)，關係式如下： 𝑖−1. 𝑇𝑖,𝑢. 𝐿𝑖(𝑖−1),𝑢 𝐿𝑖𝑗,𝑢 𝐿𝑖1,𝑢 𝐿𝑖2,𝑢 =∑ = + + ⋯+ 𝑉𝑗 𝑉1 𝑉2 𝑉𝑖−1 𝑗=1. 11. (2-4).

(29) 𝑖. 𝐿𝑖𝑗,𝑙 𝐿𝑖1,𝑙 𝐿𝑖2,𝑙 𝐿𝑖𝑖,𝑙 𝑇𝑖,𝑙 = ∑ = + +⋯+ 𝑉𝑗 𝑉1 𝑉2 𝑉𝑖. (2-5). 𝑇𝑖,𝑙 = 𝐷𝑇𝑖 + 𝑇𝑖,𝑢. (2-6). 𝑅𝑖,𝑙 × 𝑍𝑗 𝐷𝑖,𝑙. (2-7). 𝑗=1. 𝐿𝑖𝑗,𝑙 = 𝑉𝑖 =. 𝐿𝑖𝑖,𝑙. (2-8). 𝑇𝑖,𝑙 − ∑𝑖−1 𝑗=1. 𝐿𝑖𝑗, 𝑉𝑗. 此法雖然比間距法多考慮了波傳路徑上有不同的分層速度，但卻假設為直線路徑。 (5) 逆推法 (inversion method) 逆推法是建立在修正間距法的基礎上，除了考慮到分層速度不同以外，更考慮到震波在水平層狀地層的傳遞路徑會遵守司乃爾定律 (Snell’s law) ，在不同介質下產生折射(圖 2.7)，關係式如下： sin 𝜃𝑖𝑗 sin 𝜃𝑖1 sin 𝜃𝑖2 sin 𝜃𝑖𝑖 = =⋯= =⋯= 𝑉1 𝑉2 𝑉𝑗 𝑉𝑖. (2-9). 𝑍1 tan 𝜃𝑖1 + ⋯ + 𝑍𝑗 tan 𝜃𝑖𝑗 + ⋯ +𝑍𝑖 tan 𝜃𝑖𝑖 = 𝑆. (2-10). 𝐿𝑖𝑗 =. 𝑍𝑗 cos 𝜃𝑖𝑗. (2-11). 只有最上層的波傳採用直線路徑，以下的地層則採用司乃爾定律折射，並利用迭代的方法對各層波速進行反算，直到反算前後速度誤差小於 0.01%。. 12.

(30) 圖 2.4 直接法 (direct method)示意圖(Kim et al., 2004). 圖 2.5 間距法 (interval method)示意圖(Kim et al., 2004). 13.

(31) 圖 2.6 修正間距法 (modified-interval method) 示意圖(Kim et al., 2004). 圖 2.7 逆推法 (inversion method) (Kim et al., 2004). 14.

(32) 而下孔震測在施作時，通常會製造一左一右方向的水平敲擊，對受波器來說即是受到正反向的水平震波，由於剪力波的特性會使得正反向震波振幅對稱，故藉由這樣的敲擊方法幫助區別出剪力波初達時間，如圖 2.8。. 圖 2.8 下孔震測利用左敲、右敲輔助凸顯剪力波的對稱性(D7400-08, 2008). 15.

(33) 2.1.2 跨孔式震測法跨孔式震測法乃地球物理探勘於土層內部構造的應用之一，震源以及受波器皆配置於地表以下。且至少需要兩個垂直地表面的鑽孔，以三個鑽孔為最佳，在其中一個鑽孔設置震源，而剩下的鑽孔則設置受波器。根據 ASTM D4428/D4428M-07 (2004)之規範，受波器以及震源需保持平行，亦即兩者必須維持等深度(如圖 2.9 所示)，以接受震源所發出的壓力波和垂直向剪力波訊號，同時依據相鄰兩鑽孔之間距及波傳走時，計算兩鑽孔間土層之平均波傳速度。並藉由改變 0.5 或 1 公尺的垂直間距以探測鑽孔深度範圍內之土層之波速剖面。跨孔之水平間距不宜過長，以避免受波器收到之初達波為折射波而非剪力波，例如在夾於兩高速地層中之低速層，便會因折射現象造成波傳以最短路徑通過高速層，而高估了低速層之速度。 Roblee et al. (1994)曾指出跨孔法相對其他震測法有以下優勢： (1) 只要震源和接收器的距離固定，量測的精度和深度無關， (2) 每次量測的結果不會互相干擾，可視為獨立的， (3) 可確認測孔間的土壤是未受擾動的， (4) 各種類型的波都會被量測到， (5) 可以變化接收路徑的角度，進行跨孔震測法。除此之外，跨孔式震測法的一個特點是可以藉由分別移動震源及受波器的位置，進行孔間的斷層剖面速度影像化，見圖 2.10。. 16.

(34) 圖 2.9 跨孔式震測法配置圖(D4428/D4428M-07, 2004). 17.

(35) 圖 2.10 利用跨孔式震測法施作斷層剖面速度影像方法(Angioni et al., 2003). 18.

(36) 2.1.3 懸盪式震測法上述之孔內震測方法基本上其震源與受波器為各自獨立的單元。而日本 OYO 公司所研發的懸盪式波速量測系統(suspension PS logger)，將震源與受波器結合為一體。震源與受波器間距離固定而不因深度不同而改變，且保持緊密相連，可作連續性之量測。因震源與受波器距離相近，震源能量不會因距離過大而消散不辨識；且不易受地層變化影響，深層量測時所得之訊號解析度較傳統量測方式清晰，準確性佳。(蔡璧嬬, 2007) 懸盪式波速量測系統之完整架構圖如圖 2.11 所示，包括控制監控裝置 (suspension PS logger)、絞盤(winch)以及經由訊號傳輸電纜所連接、懸吊於固定支架上之探測管(suspension probe)。探測管由上而下可再分為放大器 (head reducer)、上⁄下受波器(upper and lower receiver)、濾波軟管(filter tube)、震源(source)、震源驅動器(source driver)、及呆重(weight)等單元。此探測管這種一體化的創新設計，可直接懸吊於孔內水中量測 P 波和 S 波速度，而不需將受波器緊靠於孔壁。兩組受波器的間距是 1 米，其可同時量測 P 波和 S 波的傳遞時間。又此種懸吊式受波器的構造幾乎和水的密度相等，故即使在水中亦可有效量測 S 波。而震源是使用間接振動型來產生 S 波，利用水當介質。它的原理如圖 2.12 所示是利用水平螺線線圈及擊槌在孔內液體中產生一壓力波，然後經過孔壁轉換成實體波採放射狀的傳播到岩體。當此實體波傳遞到每一個感應器所相對應之孔壁時，在孔內液體中先轉換成壓力波，再透過受波器接收訊號。. 19.

(37) 圖 2.11 懸盪式波速量測系統配置圖(OYO corporation). 圖 2.12 震源原理示意圖(OYO suspension PS logger Manual) 20.

(38) 由於上⁄下受波器與震源距離不同，在不同時間接受到震波訊號，此訊號經由探測管上端的放大器放大後回傳於控制監控裝置中，記錄的訊號經由走時分析，判讀上⁄下受波器收錄震波之時間間隔，即可計算上⁄下受波器間所在之土層的波速度資料。由於壓力波與剪力波組成之頻率範園不同，每一組受波器由壓電式地震檢波器(hydrophone)以及兩個磁電式地震檢波器(geophone)組成，其中壓電式地震檢波器收錄壓力波，另兩組磁電式地震檢波器則以反向配置，其目地在收錄剪力波並檢核其正確性，同時亦做為剪力波初達位置之研判使用。將訊號依時間序列可繪成如圖 2.13 之訊號波形，橫座標為時間軸 (單位為毫秒 ms)。下接收器(H2、/H2)距離震源較近，較早接收到訊號，上接收器(H1、/H1)接收到訊號時間較慢；只要判讀出受波器記錄之初達時間，上下接收器之初達時間差即為震波經過上下接收器間土層之時間間隔，進而計算波速。. 圖 2.13 懸盪式震測法訊號波形(國家地震工程研究中心) 懸盪式波速量測系統之基本理念為利用遠近受波器間之時間差以求取該段距離之地層波速，若其中一接收器之訊號太差而無法反應初達波位置，可藉由震源到另一接收器之初達時間做為依據，協助判定波速。且一般會 21.

(39) 藉由深度序列資料加以輔助判讀初達波。所得的結果也可利用上述方法進行交叉驗證以定出真實波速度值。 2.2 走時分析方法與評析以上所提的震測方法皆屬於走時分析的方式，也是所謂初達波式分析方法。需先經人為判定其初達波到時，才能進一步計算剪力波速度，而初達波是一種特殊類型波，具有起跳時間最早、能量較強的特性；在震波資料處理上，精確且有效率地挑取初達時間對於其後續處理有著重要的意義。然而，當鑽孔條件不良、鑽井井壁與套管間存在縫隙、受到特定管波干擾、或是任何可能影響剪力波訊號的因素發生，若擷取到剪力波訊號雜亂或受到干擾而不易判別，無法精準的定出初達波到時，則影響剪力波速度結果的正確性。(蔡璧嬬, 2007) 因此，初達時間的判釋結果影響計算出的波傳速度的正確性。為獲得準確的初達波位置，有下列幾個要點： (1) 判別初達波到達位置需觀察同一受波器之正向剪力波及反向剪力波訊號是否確實反向對稱。 (2) 判別初達波到達位置需確認其受波頻率屬地層材料之剪力波無誤。 (3) 近地表地層處波形常受到管波(tube wave)影響，在真實初達波到達之前已有管波序列存在，造成該處之波速分析易產生誤判。此時可借助較深處之波速分析結果加以推算淺層處之初達波位置，再自估算之初達波之後找尋對應波序。 (4) 上⁄下受波器之剪力波資料有主控頻率衰變現象，並且兩受波器之衰變情形彼此相異，故以上⁄下受波器對應之波序決定波速時應儘量找尋初達位置附近者。當波形訊號內容符合下列條件時，可被認定為品質良好之資料，所得之剪力波速度也較可信賴： 22.

(40) (1) 剪力波形良好，無雜訊干擾。 (2) 同一受波器之正向剪力波及反向剪力波訊號確實反向對稱。 (3) 初達波明確可辨。目前的走時分析大多由人工手動挑波為主，雖說有上述的挑波概念可以參考，但也承如上述所言，實際資料常受到雜訊干擾，使得挑波不易，相較之下，挑波的經驗成為關鍵，每個人所挑的初達時間也不盡相同，常見問題如下： Stokoe 等學者(Hoar and Stokoe, 1981; Stokoe and Hoar, 1978; Stokoe and Woods, 1972)利用跨孔式震測法量測剪力波速，其文獻顯示初達時間分析方式的差異會產生不同的剪力波速，如圖 2.14。在相同的間距下使用不同的選點方式計算剪力波速，將波速正規化之結果比較，顯示選點方式的差異就會造成求得的波速不同。(曾志瑋, 2009). 圖 2.14 不同選點方式的正規化波速比較(Stokoe and Hoar, 1978) 而 Sánchez-Salinero (1987)的研究中也發現，在二維的模型中，在 d/λ <2（d：震源至接收處的距離；λ：波長）、無阻尼、柏松比為 0.25 的條件 23.

(41) 下，會產生近場效應，如圖 2.15，即為同平面縱向運動的 P 波受到同平面橫向運動 SV 波干擾的情形；圖 2.16 則是同平面橫向運動 SV 波受到同平面縱向運動 P 波干擾的情形。. 圖 2.15 P 波中 SV 波的近場效應(Sánchez-Salinero, 1987). 圖 2.16 SV 波中 P 波的近場效應(Sánchez-Salinero, 1987) 24.

(42) 以上兩種情形皆是手動挑波方法常見的問題，使得挑波準則不一。除此之外，通常震波資料分析工作量大，以三向度的受波器來說，上百筆至上千筆的資料皆有可能，如此，初達波的挑取即變成一項費時、費工、十分繁雜的工作。因此也一直有學者想發展出一套自動化的挑波方法，不論在地震預警、或是震測分析上都有學者陸續提出新方法。 2.3 自動挑波方法與評析如上小節所述，不論是在地震警示、或是震測分析上都有學者陸續提出新的自動挑初達波方法，從最早開始也最普遍應用的 STA/LTA 方法、利用統計模型準則的 AIC 方法、或用類神經網路訓練判讀波相到時、以及近年來發展出的能量比值法及其修正方法。但是目前為止並沒有一種單一的方法可以辨識所有的震相，而且每種方法都有一定的侷限性和其應用範圍，必須用綜合分析的方式在才能在一定程度上彌補其缺點，利用多種方法聯合辨識，綜合其辨識結果，才能得出各個震相比較精確的到時。 2.3.1 STA/LTA 方法 STA/LTA 方法是目前地震預警系統中常用的 P 波判斷方法之一，由於地震訊號和背景雜訊在震幅特徵和頻率組成有很大不同，地震訊號特徵大多有明顯初達波、高頻且衰減快；而背景雜訊則可由低振幅、低頻的特徵來區分，故藉由定義 STA/LTA 的特徵函數(CF)，來表示所提取的訊號與背景雜訊的特徵不同來自動判讀波相。其中 STA 即為「短時窗訊號平均值」，可衡量地震訊號的變化，用來監測地震，變化較快；而 LTA 為「長時窗訊號平均值」，用來衡量背景雜訊，變化較緩慢。當地震波到達時，STA/LTA 比值會遽變，當超過設定的觸發門檻值(THR)時，則判識為 P 波到達。其中，特徵函數的選取直接影響了波相挑取的精度，選取方法有很多種，包括 Ambuter and Solomon (1974)、Anderson (1978)、McEvilly and Majer (1982) 等利用地震紀錄的振幅絕對值|𝑥𝑖 |作為特徵函數；Swindell and Snell (1977) 25.

(43) 則利用振幅的平方值𝑥 2 作為特徵函數；Earle and Shearer (1994)利用 E(𝑡) = √𝑥(𝑡)2 + 𝑥̅ (𝑡)2 作為特徵函數，其中𝑥(𝑡)是地震紀錄，𝑥̅ (𝑡)是其 Hilbert 轉換；而 Allen (1978)則是利用E(𝑡) = 𝑥(𝑡)2 + 𝑥 ′ (𝑡)2 + 𝐶2 作為特徵函數，其中𝑥(𝑡)為地震訊號，𝑥 ′ (𝑡)則為其一階微分。藉由這些特徵函數可以在時間序列的移動計算 STA 和 LTA，並根據場址的背景雜訊來調整對雜訊的敏感度。 STA/LTA 具體算法如下： 𝑖−𝑛𝑠. 1 STA(𝑖) = ∑ 𝐶𝐹𝑗 𝑛𝑠. (2-12). 𝑗=𝑖. 𝑖−𝑛𝑙. 1 LTA(𝑖) = ∑ 𝐶𝐹𝑗 𝑛𝑙. (2-13). 𝑗=𝑖. STA/LTA(𝑖) =. 𝑆𝑇𝐴(𝑖) 𝐿𝑇𝐴(𝑖). (2-14). 上式中，STA(𝑖)和LTA(𝑖)分別代表在𝑖時刻的短時窗訊號平均值和長時窗訊號平均值；CF 為在 i 時刻的特徵函數取值；𝑛𝑠及𝑛𝑙分別為短時窗和長時窗內的資料點數，如圖 2.17。. 圖 2.17 STA/LTA 方法示意圖(Han et al., 2010) 此法的優點是計算方便、費時短，適合即時處理；然而也存在一些缺點，當訊雜比較低、或者初達波不明顯時判讀效果不好。例如，如果在 P 波初達前存在一些脈衝雜訊及訊號異常跳動引起的振動，STA 值將變大後又迅 26.

(44) 速變小，因此 STA/LTA 比值將超過觸發門檻值，又突然降下，因而造成誤判。除此之外，短時窗長度相對於長時窗長度越小，STA/LTA 判讀敏感度越高，設置門檻值(THR)應該大一些；短時窗長度相對於長時窗長度越大，則敏感度下降，設置門檻值(THR)則應小一些。因此如何選取 STA、LTA、 THR 無疑對判讀初達波精度有重要影響。 2.3.2 AIC 方法 AIC (Akaike’s Information Criterion)為赤池資訊準則的縮寫，是一種衡量統計模型擬合優良性的標準，由日本統計學家赤池弘次(Akaike, H.)創立與發展。赤池資訊準則建立在熵(Entropy)的概念基礎上，可以權衡所估計模型的複雜度和此模型擬合數據的優良性。(Akaike, 1974) 而 AIC 在震相辨識的應用上，即是求解背景雜訊及地震訊號最佳劃分點的過程，此點與 AIC 曲線極小值點相對應，故 AIC 曲線的極小值即為震相到時點，如圖 2.18。. 圖 2.18 AIC 方法判讀兩種波相到時示意圖(Zhang et al., 2003) 27.

(45) 常用的 AIC 方法有：AR-AIC 方法、基於類神經網路的 AIC 方法、及基於小波轉換的 AIC 方法等。AIC 一般計算式如下： AIC = 2k − 2 ln 𝐿. (2-15). k 是模型參數的數量，L 是模型的最大概似函數。 AR-AIC 方法(Leonard, 2000; Leonard and Kennett, 1999; Sleeman and van Eck, 1999)是先假設可以把地震紀錄分成兩個平穩過程，在每個部分分別建立自回歸模型，這兩個 AR 模型是完全不同的，即一個只包含背景雜訊而另一個只包含地震訊號。當 AIC 取最小值時，AR 模型取最優的階數，此時 AR 模型與地震紀錄這個時間序列最相符，也就是雜訊與地震訊號的最優劃分點，即為波相到時點，計算式如下： AIC(k) = −(𝑘 − 1) log 𝑉𝐹2 − (𝑛 − 𝑘 + 1) log 𝑉𝐵2 + 2(𝑚𝐹 + 𝑚𝐵 ). (2-16). 𝑉𝐹2 為雜訊模型的變異數；𝑉𝐵2 為訊號模型的變異數；𝑚𝐹 為雜訊 AR 模型係數的個數，𝑚𝐵 為訊號 AR 模型係數的個數。以上做法是假設地震紀錄中的背景雜訊和地震訊號是平穩的，而實際是非平穩的，所以對於低訊雜比且初達波較隱沒的地震訊號，所估計的初達時間誤差較大。此外兩端的 AR 模型係數還必須用試誤法(try and error)才能決定，相當麻煩。而 Maeda (1985)提出了不同於 AR-AIC 算法，其做法是直接從地震記錄圖中計算 AIC 值，而不取 AR 係數。其計算式為： AIC(k) = k log(𝑉𝑎𝑟𝑖𝑎𝑛𝑐𝑒(𝑥[1, 𝑘])) + (𝑛 − 𝑘 + 1) log(𝑉𝑎𝑟𝑖𝑎𝑛𝑐𝑒(𝑥[𝑘 + 1, 𝑛])). (2-17). 然後在選定的時窗內對 k 逐點搜索，AIC 極小值的點即為 P 波到時點。但是對於訊雜比低時且初達波隱沒的訊號也是表現不佳，除此之外，AIC 方法如何能挑選到合理的初達時間點受到時窗的選取影響很大。 Zhang (2001)則採用 Maeda (1985)的計算公式，分別結合了多層感知器類神經網路(Multilayer Perceptron Neural Network)，成為 MLP-AIC 方法、以及 28.

(46) 結合小波轉換 Zhang et al. (2003)的 W-AIC 方法來進行波相辨識，藉此改進受時窗選取影響的問題。但是這些方法都只能針對 P 波震進行辨識，還不能用於 S 波震相辨識。 2.3.3 類神經網路然而不論是何種地震研究，精確的判讀出地震波的到時永遠是最基礎、也是最耗時的工作之一。傳統上利用電腦判讀地震波到時的方法有許多種，其中最廣為流傳的是所謂 STA/LTA 運算程序(Allen, 1978; Baer and Kradolfer, 1987)，即 2.3.1 所述。此運算程序根據地震儀的特徵週期來設計其長、短時窗的時間長度，因此用於傳統窄頻地震儀（如長週期或短週期地震儀）所得之地震記錄波形時，可以得到非常好的結果。然而若應用於寬頻地震波形資料時，卻有誤判的可能(Zhao and Takano, 1999)。由於寬頻地震波形記錄的頻率範圍較寬，因此傳統的 STA/LTA 運算程序無法找出一個適當的長、短時窗的時間長度，藉以判斷地震波的到時。因此，發展一個適當的運算程序，能夠從寬頻地震波形資料中讀取正確的到時，成為一項必要的工作。為了達到上述目的，日本東京大學的研究人員發展了一套類神經網路的寬頻地震資料之波相到時判讀程序(Zhao and Takano, 1999)。整個判斷程序的核心是由倒傳遞類神經網路（Backpropagation Neural Networks, BPNNs）所組成，而這種類似的波相到時判讀器在過去也曾應用於短週期地震資料 (Dai and MacBeth, 1995; Murat and Rudman, 1992; Wang and Teng, 1995)，並且相當成功。但用於寬頻地震資料方面，則必須加以修改。而這個程序是對同一組輸入信號，同時設置長、中、短三個時間長度不同的波相到時判斷器，希望同時利用短時窗能準確決定波相到達時間之特性及長時窗能準確判斷是否有波相到達的優點。並同時利用數組判斷標準來加以檢視，再將數個判斷結果以一個並聯判斷器連接，在最少誤判的情況下，藉以決定是否有地震波相到達，及其正確的到達時間，見圖 2.19、 29.

(47) 圖 2.20。然而該研究亦指出，此種類神經網路的判讀程序其成功關鍵在於適當的學習過程，特別是針對資料的特性，應有不同數目的判斷標準及相關的參數，而這些數目及參數的設定則必須經由實驗得到，並無一定的規定，故對於不同地區的寬頻地震資料會有不同的參數結果。. 圖 2.19 長、中、短三個時窗判斷器，並以一並聯判斷器取其交集判斷 (Zhao and Takano, 1999). 30.

(48) 圖 2.20 長、中、短三個時窗判斷器示意圖(Zhao and Takano, 1999) 而高弘 (2001)也嘗試利用此法，希望發展一套適用於台灣地區寬頻地震資料的類神經網路波相到時判讀程序，並將這套程序實際應用於現有的 BATS(Broadbadnd Array in Taiwan for Seismology)寬頻地震資料，驗證其未來全面使用的可行性。但究其結果卻不盡理想，因此提出以下結論： (1) 類神經網路技術確實能達到判別地震波相的目的，但是利用單一時窗的判讀器則會產生許多誤判的情形。這個情況可以長、中、短三個時窗同時判讀後取其交集而獲得顯著改進。 (2) 類神經網路波相到時判讀器必須經過相當程度的訓練，而訓練資料的準備及其是否具有足夠的資訊提供類神經網路學習，以便能於實際應用時達到預期的表現，則是這個技術是否成功的關鍵之一，必須特別加以注意。 (3) 將受過相當訓練的類神經網路波相到時判讀器應用於 BATS 所紀錄到的寬頻地震波形，可以得到滿意的結果。但因類神經網路本身的計算較 31.

(49) 為複雜，不利於用來處理即時資料，因此這部分的應用技術仍有待突破。 (高弘, 2001) 2.3.4 能量比值法在地震紀錄上，初達波是一個很特殊類型的波，具有跳起時間最早，能量強的特性；而在初達時間點之前的地震有效訊號為零，存在的只是雜訊，而在它之後是非常重要的地震訊號。因此初達時間點前後的能量特徵有非常大的差異，故 Coppens (1985)提出了能量比值法(Energy Ratio)，其定義為一個週期內的訊號能量與總時窗能量的比值，即： 𝑟. R(𝜏) =. ∫𝑟−𝐿 𝑥 2 (𝑡)𝑑𝑡. (2-18). 𝑟. ∫0 𝑥 2 (𝑡)𝑑𝑡. 式中，R(𝜏)為能量比值函數，𝑥(𝑡)為記錄訊號，L 為視週期的長度。離散後的公式為： 𝐴𝑖 =. 2 ∑𝑇𝑡=𝑇 𝑥 2 (𝑡) 0. (2-19). 0 ∑𝑇𝑡=𝑇 𝑥 2 (𝑡) 1. 式中，對資料𝑥(𝑡)沿時間方向取兩個時窗：設𝑇1 為前一時窗起點，𝑇0 為前一時窗終點，亦即為後一時窗起點，𝑇2 為後一時窗終點。依照上述公式，傳統的計算上，時窗從頭到尾是被固定劃分成多個能量搜索窗，窗口大小由使用者自己定義，各個時窗長度相同，通過比較某一時窗與其前面時窗能量大小的比值，來求取地震紀錄的初達時間。但這種固定時窗的方法，若時窗長度選取不當或首個時窗的起始點處於特殊位置及雜訊較大時，都會對準確度造成影響，如圖 2.21。故有多位學者對此提出了改進方法，其中最常用的滑動時窗能量比值法 (左国平 et al., 2004)，滑動時窗是為了克服固定時窗的缺點提出的，其主要原理是，時窗並不是固定不變的，而是從訊號時間序列的頭到尾，如圖 2.22 所示。∆x為每次移動的擷取間隔，並根據實際情況進行調整。時窗每移動 32.

(50) 一次，就計算一次前後能量大小的比值，然後從這些能量比值中找出最大值，再從對應的後時窗找出振幅最大的資料點，所對應的時間即為初達時間。但滑動時窗法也存在著缺點，例如時窗大小的選擇直接影響了挑取時間的準確性，不同時窗長度挑取的結果往往差別很大，這位正確選取時窗長度增加了困難度。因此张伟 et al. (2009)提出了變換時窗長度統計法來克服以上兩種方法的缺點。它利用改變不同長度的時窗來計算一系列的初達時間，用統計的方法選擇初達最集中的那個時間來做真正的初達時間。這樣的方法相較之下有效地提升了挑取精度，使時窗長度設置不當等人為因素引起的誤差變小了。. 取樣點數圖 2.21 固定時窗(罗进华 et al., 2008). 時窗長度. 滑動方向. 取樣點數圖 2.22 滑動時窗示意圖 (徐钰 et al., 2012) 33.

(51) 以上方法主要針對能量比值法的擷取時窗造成的問題提出修正，另一方面 Wong et al. (2009)則針對 Coppens (1985)所提出的式子做修正，採用了修正的能量比值法(Modified Energy Ratio)，簡稱 MER。以能量比值法的概念來說，假設𝑥𝑖 代表時間序列𝑖 = 1,2 … . N上對應的資料點，時窗內的資料點數為 L，則在時間序列𝑖時，其前後時窗內資料的能量比值為: er =. 2 ∑𝑖+𝐿 𝑗=𝑖 𝑥𝑗. ， ( Energy Ratio ). (2-20). ∑𝑖𝑗=𝑖−𝐿 𝑥𝑗 2. 為了讓能量比值法對於訊雜比低的資料能夠更凸顯，進行了以下修正項： mer𝑖 = [𝑒𝑟 ∗ 𝑎𝑏𝑠(𝑥𝑖 )]3 ， ( Modified Energy Ratio ). (2-21). MER 方法示意如圖 2.23：. 圖 2.23 修正能量比值法(MER)示意圖 Han et al. (2010)亦利用此法跟 2.3.1 所提到的 STA/LTA 方法做比較，顯示在訊雜比低時( SNR = 3.5 )，MER 能判讀出較精確的初達時間，但當訊雜比更低時( SNR = 1.5 )，則必須依靠中值濾波器、套用時窗、或是重複挑取的方式輔助，以獲得較精確的初達時間。. 34.

(52) 2.4 訊號分析方法由上一小節我們可以得知，其實現有的自動化分析方法並不能適用於所有的情況，也沒有一個單一的方法可以完整解決，甚至訊雜比較低時也會有判斷錯誤的情形。為此，本研究希望透過現有的訊號分析的方法，進一步詳觀震波訊號的特性，分別從時間域濾波的濾波、頻率域的濾波方法、以及時頻的濾波分析方法來探討挑取初達波的概念，推演到自動化挑波的可行性。故以下大致分成三大類進行回顧。. 2.4.1 時間域濾波分析窗函數(Window Function)是頻譜分析中一個重要的部分。在訊號處理中，窗函數是一種除了在給定區間之外取值均為 0 的實函數。例如：在給定區間內為常數而在區間外為 0 的窗函數即稱矩形窗。任何函數與窗函數的積仍為窗函數，所以相乘的結果就像透過窗口「看」其他函數一樣。用於修正訊號的非週期性及因訊號遺漏造成的不準確。訊號處理的主要數學工具是傅立葉轉換(Fourier Transform)，它是研究整個時間域和頻率域的關係。然而，當運用電腦實現訊號處理時，不可能對無限長的訊號進行測量和運算，而是取其有限的時間片段進行分析，故當輸入波形是採樣信號而非連續信號時，傅立葉分析通常對訊號擷取並用離散傅立葉變換。做法是從訊號中截取一個時間片段，這個區段就叫時窗 (Time Window)，將每段的序列當成一個固定的訊號，並進行週期延拓處理，得到虛擬的無限長的信號，然後就可以對信號進行傅立葉轉換、相關分析等數學處理。但由於擷取出來的時窗不易為訊號週期的整數倍，故在時窗間的接點處會有不連續的情況發生，造成遺漏(Leakage)效應。如圖 2.24 所示：. 35.

(53) 圖 2.24 闡述非週期訊號取樣時造成遺漏的原因(CRYSTAL instruments) 圖中可以看到時域上明顯沒有按週期截斷，才造成頻率域上的能量遺漏。無限長的訊號被截斷後，其頻譜發生改變，原本集中在 f0 處的能量被分散到兩個較寬的頻帶中去了，這種現象稱為頻譜能量遺漏(Proakis and Manolakis, 1992)。. 圖 2.25 由於資料截斷造成頻譜能量遺漏(测试技术实验指导) 為了克服時窗的不連續，傳統方法中是對訊號乘上一個窗函數(Window 36.

(54) Fuction)，使時窗的接點可以更加平滑且有連續性。而訊號截斷的作法也就相當把訊號乘上一個矩形窗 w(t)，如圖 2.25。遺漏和窗函數頻譜的兩側旁瓣有關，若兩側半的高度趨於零，而使能量相對集中在主瓣，就可以較接近真實的頻譜，因此可在時間域採用不同的窗函數來截斷訊號。而加窗在時域上採用點乘，在頻率域則為摺積。摺積可被看成是一個平滑的過程。此平滑過程可視為是由一組具有特定函數形狀的濾波器，因此，原始訊號中在某一頻率點上的能量會結合濾波器的形狀表現出來，進而減少遺漏，基於這個原理，故通常直接在時域加窗。但由於加窗計算中衰減了原始訊號的部分能量，也因此對於最後的結果還須加上修正係數. 37.

(55) 2.4.2 頻率域濾波分析頻譜是指一個時域的信號在頻率域下表示方式，可以針對信號進行傅立葉轉換而得，所得的結果會是以分別以振幅及相位為縱軸，頻率為橫軸的兩張圖，不過有時也會省略相位的資訊，只有不同頻率下對應振幅的資料。有時也以「振幅頻譜」表示振幅隨頻率變化的情形，「相位頻譜」表示相位隨頻率變化的情形。簡單來說，頻譜可以表示一個訊號是由哪些頻率的弦波所組成，也可以看出各頻率弦波的大小及相位等資訊。故在頻率域下，我們採用頻譜分析(Spectral Analysis)，是一種將複雜訊號分解為較簡單訊號的技術。許多物理訊號均可以表示為許多不同頻率簡單訊號的和，而找出一個訊號在不同頻率下的資訊（可能是振幅、功率、強度或相位等）的作法就是頻譜分析。頻譜分析可以對整個訊號進行。不過有時也會將訊號分割成幾段，再針對各段的訊號進行頻譜分析。週期函數最適合只考慮一個週期的訊號來進行頻譜分析；而傅立葉轉換中有許多分析非週期函數時需要的數學工具。一個函數的傅立葉轉換包括了原始訊號中的所有資訊，只是表示的型式不同，因此可以用反傅立葉變換重組原始的訊號。若要完整的重組原始訊號，需要有每個頻率下的振幅及其相位，這些資訊可以用二維向量、複數、或是極座標下的大小及角度來表示。在訊號處理中常常考慮振幅的平方，也就是功率，所得的就是功率譜密度(power spectral density, PSD)。實際上，大部份的儀器及軟體都用快速傅立葉變換(FFT)來產生頻譜的訊號。快速傅立葉變換是一種針對取樣訊號計算離散傅里葉變換的數學工具，可以近似傅立葉變換的結果。而表面波震測的分析方式也採用了頻譜分析法，亦即表面波譜法 (Spectral Analysis of Surface Wave, SASW)，其施測方式主要由一衝擊式震源及兩個受波器所組成(Nazarian et al., 1983; Stokoe et al., 1994)，假設受波器 38.

(56) 訊號可表示為位置與時間的函數𝑢1 (𝑥1 , 𝑡)，𝑢2 (𝑥2 , 𝑡)。透過傅立葉轉換可得到此二函數在頻率域的函數𝑈1 (𝑥1 , 𝜔)，𝑈2 (𝑥2 , 𝜔)以及各頻率在𝑥1 ，𝑥2 處的相位角。根據相位速度的計算式： V𝑝ℎ (𝜔) =. 𝜔 ∆𝜙(𝜔) ∆𝑥. (2-22). 其中V𝑝ℎ (𝜔)代表不同頻率時的雷利波波速；∆𝜙(𝜔)表示不同頻率時在𝑥1 ，𝑥2 處的相位角差值；∆𝑥表示𝑥2 − 𝑥1 。計算過程上兩受波器各頻率之相位角差(∆𝜙(𝜔) = 𝜙2 (𝜔) − 𝜙1 (𝜔))可由互能頻譜(Cross Spectral Density, CSD)之相位角獲得，即 ∆𝜙(𝜔) = Angle (𝐶𝑆𝐷 (𝑢1 (𝑥1 , 𝑡)，𝑢2 (𝑥2 , 𝑡))). (2-23). 如圖 2.26 虛線所示，兩受波器之相位角差，一般隨頻率增加而增加，但由於頻譜分析所得之相位角差僅侷限於−π~π之間，因此在計算(2-22)式之前必須先將相位角差∆𝜙(𝜔)摺開(unwrap)，以得各頻率真正的相位角差，如圖 2.26 實線的結果。. 圖 2.26 SASW 頻譜分析(相位差-頻率圖)(林俊宏, 2005). 39.

(57) 表面波譜法對於表面波應用於地質調查之推廣扮演了舉足輕重的地位，但以現有的技術看來，其施測方式與分析方式並非最佳。然而其利用頻譜分析的方式可以自動將訊號轉至頻率域，進而計算出頻散曲線，這點對於本論文的目標是一致的，概念上也有很大幫助及啟發，尤其對於同為單一震源、兩個受波器的配置，更有異曲同工之妙；對於這部分的應用在第三章會有更詳細的介紹。 2.4.3 時頻譜濾波分析頻率為時間訊號的重要特徵，傳統用傅立葉轉換進行頻譜分析，來了解一段時間內頻譜的分佈，但有時我們更有興趣的是頻率隨時間變化的情形，分析各種不同頻率隨時間變化的情形稱為時頻分析。時頻分析相較於頻譜分析多了頻率對時間的解析，利用時頻分析可以看出頻率隨著時間變化的情形，甚至在每個時間中可以看出訊號在一定頻率附近的能量聚集狀況。其主要研究非穩態訊號或時變訊號，為當今訊號處理研究極為重視的一個領域，特別是在二十世紀的 80 年代以後有著突破性的發展，各種時頻聯合分析方法得到了廣泛的研究及運用，逐漸形成了一套獨特的理論體系。故在本小節回顧了應用於表面波震測頻散現象的多重濾波法(Multiple Filter Technique)、及目前廣為應用的小波轉換方法(Wavelet Transform)，藉此作為第三章研究的應用。. 40.