Pseudo-section概念於表面波震測應用之數值模擬探討

(1)

國立交通大學

土木工程學系碩士班

碩士論文

Pseudo-Section 概念於表面波震測應用

概念於表面波震測應用

之數值模擬

之數值模擬探討

探討

A Study of Pseudo-Section Concept

in Surface-Wave Testing, by Numerical Simulation

研究生：林俊宏

指導教授：林志平博士

(2)

Pseudo-section 概念於

概念於

表面波震測應用之

表面波震測應用之數值模擬探討

數值模擬探討

Study of Pseudo-section Concept in

Surface-Wave Testing, by Numerical Simulation

研究生：林俊宏 Student：Chun-Hung Lin 指導教授：林志平博士 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 2005

Hsinchu, Taiwan, Republic of China 中華民國九十四年七月

(3)

國

立

交

通

大

學

國

立

交

通

大

學

國

立

交

通

大

學

博碩士論文全文電子

博碩士論文全文電子檔著作權授權書

檔著作權授權書

本授權書所授權之學位論文，為本人於國立交通大學土木工程所大地工程組，93 學年度第二學期取得碩士學位之論文。論文題目：Pseudo-section 概念於表面波震測應用之數值模擬探討指導教授：林志平博士 ■ 同意 □不同意本人茲將本著作，以非專屬、無償授權國立交通大學與台灣聯合大學系統圖書館：基於推動讀者間「資源共享、互惠合作」之理念，與回饋社會與學術研究之目的，國立交通大學及台灣聯合大學系統圖書館得不限地域、時間與次數，以紙本、光碟或數位化等各種方法收錄、重製與利用；於著作權法合理使用範圍內，讀者得進行線上檢索、閱覽、下載或列印。論文全文上載網路公開之範圍及時間：本校及台灣聯合大學系統區域網路 ■ 中華民國 94 年 8 月 29 日公開校外網際網路 ■ 中華民國 94 年 8 月 29 日公開授授授授權權權權人人：人人：：：林俊宏林俊宏林俊宏林俊宏親筆簽名親筆簽名親筆簽名親筆簽名：：：：________________________________________________________________________________________ 中華民國中華民國中華民國中華民國九十四九十四九十四九十四年年年年八八八八月月月月二十九二十九二十九二十九日日日日

(4)

Pseudo-section 概念於表面波震測應用之

概念於表面波震測應用之

數值模擬探討

研究生：林俊宏指導教授：林志平博士國立交通大學土木工程學系

中文摘要

目前表面波震測法藉由對施測參數(近站支距、展距以及接收器間距) 的適度配合，已可在施測參數與施測目標的互制情況下獲得較佳結果，然而對於互制的行為卻依舊沒有得到解決。本研究將針對最近所提出的 Pseudo-section 概念解決方案以數值模擬方式評估此概念應用於表面波震測的可行性。本研究採用四級速度-應力有限差分法作為數值模擬工具，以不同的側向變化地層為例，對表面波在側向變化地層的行為，以及 pseudo-section 概念在側向變化地層中的表現進行模擬及討論。根據研究結果顯示，表面波的頻散曲線在遭遇側向變化地層時會有假頻散以及調整帶的現象產生，此現象間接的造成以水平層狀地層為假設的反算結果錯誤。此外，以 pseudo-section 概念應用於表面波震測中施測，數值模擬結果顯示在無反向行進的反射波訊號存在或此反向行進的反射波與前行的表面波訊號有足夠的距離而無干擾的情況下，可得到與傳統測線所得頻散曲線差不多或更為良好的結果。對於有反向行進的反射波的情形，應加以發展訊號分析方法將其分離出，以使 pseudo-section 概念於表面波的應用能有更加的表現，即使無法將其分離出，也可採用 pseudo-section 測線作為側向變化地層存在的檢驗手段。關鍵字：表面波、表面波震測、有限差分法、側向變化地層

(5)

A Study of Pseudo-Section Concept in

Surface-Wave Testing, by Numerical Simulation

Student：Chun-Hung Lin Advisor：Dr. Chih-Ping Lin Institute of Civil Engineering

National Chiao Tung University

Abstract

Presently, by adjusting the survey line parameters, such as near offset, receiver spacing, and offset range, surface wave testing can provide reasonable results even though tradeoffs are involved when determining the field configuration parameters. However, we still have problems with the tradeoffs. We use numerical simulation to study the application of pseudo-section concept, which is introduced to provide a solution of the tradeoffs, in surface wave testing.

Forth-order velocity-stress finite difference method is used to simulate the surface wave testing data based on conventional and pseudo-section survey line in earth models with different lateral heterogeneity. According to the results, false dispersion and interference zone exist when surface wave propagating in the earth with lateral heterogeneity. Furthermore, comparing the results based on conventional survey line and pseudo-section concept with no backward reflecting waves or well separation of the backward reflecting and forward propagating waves, the performance of pseudo-section survey line is the same or even better. When pseudo-section survey is interfered by backward reflecting

(6)

wave, it is necessarily to develop a new signal processing to filter out the backward reflecting wave to achieve better performance. In its current form, pseudo-section survey is still a good way to detect the existence of lateral heterogeneity.

Key words ： surface wave, surface wave testing, finite difference, lateral heterogeneity

(7)

致

致謝

謝

首先須感謝林志平老師於論文撰稿期間在架構以及研究辨証上的協助，由於老師的用心指導得以補足我個人在思考上的不足以及缺失，使得無經驗的我亦可有此成果，並獲得評審的肯定。此外，感謝美國聖路易斯大學地球與大氣科學系 Herrmann 教授以及中央大學地球科學系王乾盈教授在地震波波傳原理上所提供的協助，使得不曾接觸地震波波傳的我得以解決相關問題；同時亦要感謝中央大學地球科學系陳浩維教授在最後於數值模擬上的經驗分享以及程式勘誤，而使本研究得以順利完成。最後感謝在我身圍週遭的家人以及朋友，因為有你們時時的鼓勵與打氣才讓我有不斷前進的動力。再次感謝在我撰文期間所有給予我協助的人，因為你們熱心的支援，而造就了本論文，感謝化與再多也不嫌，向各位致上最敬意。

(8)

中文摘要中文摘要中文摘要中文摘要...I Abstract ... II 致致致致謝謝謝 ...IV謝目錄目錄目錄目錄 ... V 表目錄表目錄表目錄表目錄 ...VIII 圖目錄圖目錄圖目錄圖目錄 ...IX 第一章第一章第一章第一章緒論緒論緒論緒論... 1 1.1 研究動機... 1 1.2 研究目的... 3 1.3 研究內容... 4 第二章第二章第二章第二章文獻回顧文獻回顧文獻回顧文獻回顧... 5 2.1 震波種類與其基本理論... 5 2.1.1 壓縮波與剪力波... 7 2.1.2 表面波及其頻散曲線... 12 2.2 表面波量測方法及其施測問題... 30 2.2.1 表面波量測方法... 30 2.2.1.1 SASW... 30 2.2.1.2 MASW ... 33 2.2.2 多頻道表面波震測之施測問題詳析... 36 2.3 Pseudo-Section 概念 ... 38 2.4 P-SV 震波模擬方法 ... 40 2.4.1 解析方法... 40

(9)

2.4.1.1 射線理論... 42 2.4.1.2 空間頻率/慢度積分法... 44 2.4.1.3 震態疊加法... 45 2.4.2 數值方法... 47 2.4.2.1 有限差分法... 48 2.4.2.2 有限元素法... 56 2.4.2.3 擬頻譜法... 63 第三章第三章第三章第三章研究方法與數值模擬參數研究方法與數值模擬參數研究方法與數值模擬參數研究方法與數值模擬參數... 70 3.1 研究方法... 70 3.1.1 速度-應力差分法運算式 ... 72 3.1.2 傳統測線與 pseudo-section 測線施測... 82 3.1.3 多頻道頻率波速轉換法... 85 3.1.4 地層模型設計說明... 88 3.2 表面波震測資料模擬參數設定... 91 3.2.1 表面波震測施測參數... 91 3.2.2 數值模擬參數... 95 第四章第四章第四章第四章結果與討論結果與討論結果與討論... 97結果與討論 4.1 A 類地層 ... 97 4.2 B 類地層 ... 99 4.3 C 類地層 ... 103 4.3.1 地層側向變化之影響... 103 4.3.2 pseudo-section 概念可行性探討... 115 第五章第五章第五章第五章結論與建議結論與建議結論與建議結論與建議... 122 5.1 結論 ... 122

(10)

5.2 建議 ... 123 參考文獻參考文獻參考文獻參考文獻... 125 附錄一附錄一附錄一附錄一雷利波傳遞矩陣計算式雷利波傳遞矩陣計算式雷利波傳遞矩陣計算式雷利波傳遞矩陣計算式... 140 附錄二附錄二附錄二附錄二四級有限差分法運算元推導四級有限差分法運算元推導四級有限差分法運算元推導四級有限差分法運算元推導... 142 附錄三附錄三附錄三附錄三有限差分法自由表面運算式有限差分法自由表面運算式有限差分法自由表面運算式有限差分法自由表面運算式... 144

(11)

表目錄

表 2.1 卜松比與半無限域地層中雷利波波速對照表*... 16 表 3.1 地層模型參數... 89 表 3.2 表面波震測施測參數表... 92 表 3.3 數值模擬參數表... 96

(12)

圖目錄

圖 1.1 表面波震測法步驟示意圖... 1 圖 1.2 多頻道表面波震測法現地施測參數圖... 2 圖 1.3 pseudo-section 概念與傳統測線頻散曲線結果比較 ... 3 圖 2.1 實體波示意圖 (Bolt, 1976) ... 6 圖 2.2 表面波示意圖 (Bolt, 1976) ... 7 圖 2.3 應力定義示意圖... 8 圖 2.4 卜松比與半無限域雷利波速關係圖... 16 圖 2.5 水平層狀地層模型... 18 圖 2.6 多震態雷利波頻散曲線－Thomson-Haskell 法... 28 圖 2.7 多震態雷利波頻散曲線－Delta Matrix 法 ... 29 圖 2.8 SASW 現地施測示意圖... 30 圖 2.12 MASW 現地施測示意圖 ... 33 圖 2.13 多頻道波譜法頻散曲線分析(相位-空間位置圖)... 34 圖 2.14 多頻道波譜法頻散曲線分析結果... 34 圖 2.15 多頻道波場轉換法頻散曲線分析結果... 35 圖 2.16 空間映頻問題示意圖... 36 圖 2.17 空間映頻於多頻道波場轉換法影響示意圖... 37 圖 2.18 資料遺漏於多頻道波場轉換法影響示意圖... 37 圖 2.19 pseudo-section 概念示意圖... 39 圖 2.20 廣義射線理論示意圖... 43 圖 2.21 複數空間頻率平面積分路徑示意圖... 47 圖 2.22 有限差分法網格系統示意圖... 49

(13)

圖 2.23 速度錯置網格架構，時間定義在 k+1/2 ... 54 圖 2.24 應力錯置網格架構，時間定義在 k+1 ... 55 圖 2.25 二維模型離散化示意圖... 57 圖 2.26 二維空間問題常用元素... 57 圖 2.27 節點自由度示意圖... 58 圖 2.28 元素與內插方式示意圖... 59 圖 2.29 一維空間問題內插權重係數示意圖... 60 圖 2.30 擬頻譜法網格示意圖... 68 圖 3.1 研究方法流程圖... 70 圖 3.2 自由表面邊界設定示意圖... 78 圖 3.3 吸收邊界示意圖... 79 圖 3.4 衝擊式震源α_s參數的影響 ... 80 圖 3.5 程式計算流程圖... 81 圖 3.6 MASW 傳統測線施測示意圖 ... 82 圖 3.7 MASW pseudo-section 測線示意圖 ... 83 圖 3.8 震測資料矩陣格式... 83

圖 3.9 MASW pseudo-section 測線靜態誤差示意圖(Lin et al., in press)... 85

圖 3.10 震測資料空間-頻率域矩陣示意圖 ... 85

圖 3.11 U(f_i,x_n)實部於頻率-空間上之反應 (Lin et al., in press)... 87

圖 3.12 最佳展距選取範圍 (Lin et al., in press) ... 88

圖 3.13 多頻道頻率波速轉換法頻散曲線分析的流程圖... 88

圖 3.14 A 類地層模型示意圖... 89

圖 3.15 Ｂ類地層模型示意圖... 90

(14)

圖 4. 1 A-1, A-2 傳統測線與 pseudo-section 測線震測資料 ... 97

圖 4. 2 A-1, A-2 傳統測線與 pseudo-section 測線施測頻散曲線比較圖 . 98 圖 4. 3 B-1, B-2, B-3 傳統測線與 pseudo-section 測線施測頻散曲線比較圖 ... 100 圖 4. 4 B-1, B-2, B-3 傳統測線與 pseudo-section 測線施測頻散曲線比較圖 ... 101 圖 4. 5 B-2E 表面波陣測頻散曲線 ... 102 圖 4. 6 C1-1 傳統測線震測資料圖... 104 圖 4. 7 C1-1 表面波震測頻散曲線圖... 104 圖 4. 8 C1-2 傳統測線測資料圖... 105 圖 4. 9 C1-2 表面波震測頻散曲線圖... 106 圖 4. 10 C2-1 傳統測線震測資料圖... 107 圖 4. 11 C2-1 表面波震測頻散曲線... 107 圖 4. 12 C2-2 傳統測線震測資料圖... 109 圖 4. 13 C2-2 表面波震測頻散曲線... 109 圖 4. 14 C3 傳統測線震測資料圖 ... 110 圖 4. 15 C3 表面波震測頻散曲線 ... 111 圖 4. 16 C1-1 與 C3 傳統測線施測頻散曲線比較... 112 圖 4. 17 C4 傳統測線施測震測資料圖 ... 114 圖 4. 18 C4 表面波震測頻散曲線 ... 114 圖 4. 19 C1-1 pseudo-section 測線震測資料圖... 115 圖 4. 20 C1-1 擬震測資料的反射干擾... 116 圖 4. 21 C1-2 pseudo-section 測線震測資料圖... 117 圖 4. 22 C1-2E 表面波震測頻散曲線 ... 117

(15)

圖 4. 23 C2-1 pseudo-section 測線震測資料圖... 118

圖 4. 24 C2-1 擬震測資料表面波反射訊號圖... 119

圖 4. 25 C2-2 pseudo-section 測線震測資料圖... 119

圖 4. 26 C3 pseudo-section 測線震測資料圖 ... 120

(16)

第一章

第一章緒論

緒論

1.1 研究動機

研究動機

經過數十年的發展，表面波震測法除了其非侵入性及非破壞性的優點使其在難以鑽探的地層可以快速地進行工址調查外，還因其具有分析自動化以及結果影像化的潛力而獲得親睞。目前其已成功的應用在工址調查、液化潛能分析、地層分層、地盤改良效果以及道路鋪面厚度檢測等領域上。表面波震測法概略的說，主要有三個步驟(如圖 1.1 所示)：先至 1) 現地施測，而後對施測所得的地層震動訊號，進行 2) 頻散曲線分析，最後透過 3) 地層反算，而得到地層的剪力波速剖面。圖圖圖圖 1.11.11.1 1.1 表面波震測法步驟示意圖表面波震測法步驟示意圖表面波震測法步驟示意圖表面波震測法步驟示意圖根據在現地施測時所使用的接收器數目，表面波震測法可分為表面波譜法(Spectral Analysis of Surface Wave, SASW, Heisey et al., 1982；Nazarian te al., 1984；Stoke te al., 1994)以及多頻道表面波震測法(McMechan and

(17)

Yedlin, 1981; Gabriels et al.,1978；Park et al., 1999)。多頻道表面波震測法採用在一直線上的多個接收器收錄資料，在施測上有三個施測參數需要考慮：近站支距、接收器間距以及測線展距(如圖 1.2 所示)。然而就震測目標 (探測深度深，解析度好)以及訊號分析的需求，此三個參數存在著互制的衝突性，而致使此法在施作上有一定的困難性。陳逸龍(2004)提出了對於施測參數上的建議，降低這些互制的衝突對結果的影響，並以實際試驗測試 pseudo-section 概念(Lin et al., in press)對於徹底解決這互制問題的可行性。

圖圖圖圖 1.1.1.1.2222 多多多多頻道表面波震測法現地施測參數圖頻道表面波震測法現地施測參數圖頻道表面波震測法現地施測參數圖頻道表面波震測法現地施測參數圖 Pseudo-Section 的概念是固定接收器位置，藉由多次改變震源擊發處，而後將這些資料透過縫合(seaming)的技巧，將不同震源處的資料整合，模擬出一極長展距的震測資料(Lin et al., in press)，如此一來，在施測上便可以有限的接收器，採用較小的接收器間距而又能獲得夠大的測線展距，解決傳統多頻道表面波震測法施測參數互制衝突的現象。根據其測試的結果(見圖 1.3)，以 pseudo-section 概念施測的震測資料使用較少的接收器數量，但其所得到的頻散曲線結果，雖有預期外的跳動產生，仍能與使用較多接收器一次接收的傳統測線大致吻合，增加了我們對 pseudo-section 概念在實際施作上可行性的信心。

(18)

圖圖圖

圖 1.1.1.1.3333 pseudopseudopseudopseudo----sectionsectionsection 概念與傳統測線頻散曲線結果比較section概念與傳統測線頻散曲線結果比較概念與傳統測線頻散曲線結果比較概念與傳統測線頻散曲線結果比較

在此測試中所出現預期外的跳動情形，初步檢討認為是由於地層具有的側向變化所造成的影響，然而，由於對現地情況的無法掌控，以及其他可能存在因素，如震源的不穩定、雜訊干擾或施作上的疏忽等，使得欲以實際試驗確認此跳動情形的原因有所困難。因此，為先排除其他可能因素的干擾，欲採用模擬的方式以人工製造在具有側向變化的土層中的震測資料，以對於 pseudo-section 概念在可控制地層條件下進行測試。

1.2 研究目的

研究目的

本研究將採用數值模擬的方式，以已知的地層狀況製造傳統多頻道表面波震測資料以及採用 pseudo-section 概念所收錄的表面波震測資料，測試在純粹只有側向變化地層的影響時，表面波的行為差異以及 pseudo-section 概念的使用情形，並藉由觀察具有側向變化的地層對於震測資料的影響，檢討 pseudo-section 概念應用上的可能問題，而後統合觀察到的問題，提出目前應用 pseudo-section 概念的限制及可能解決方案，以期解決或降低施測

(19)

參數互制而造成震測目標無法同時達到要求以及不佳的離散訊號分析結果的問題。

1.3 研究內容

研究內容

本研究以四個章目做為陳述，於第一章緒論中就引起此次研究的動機做說明(1.1 研究動機)，並將此研究所欲達到的成果列出(1.2 研究目的)。在第二章中，將由基本的波傳理論開始，了解在地層中傳遞的不同震波種類及其特性(2.1 震波種類與其基本理論)，而後介紹由表面波特性發展出的不同表面波震測法及其施測上的問題(2.2 表面波量測方法及其施測問題)以及 pseudo-seciton 概念是如何嘗試解決這些問題(2.3 pseudo-section 概念)。製造人造震測資料可由解析方法以及數值方法下手，在 2.4 P-SV 震波模擬一節中，將對一些常見的模擬方式做簡單介紹，藉以瞭解在模擬震波的領域上有哪些技術可供使用。經由第二章的文獻回顧後，統整收集到的文獻資料，於第三章中說明將要使用何種震波模擬方法，以及模擬震測方法的施測及其分析的詳細過程(3.1 研究方法)，而後說明本研究中模擬所用的參數(3.2 表面波震測模擬參數設定)。在第四章中以不同的模擬地層類型分別討論表面波的行為及 pseudo-section 測線所得到的結果。最後在第五章中提出本研究的結論，並以研究所得的結果提出建議。

(20)

第二章

第二章文獻回顧

文獻回顧

2.1 震波種類與其基本理論

震波種類與其基本理論

在地層中傳遞的震波大致可分為兩類，一類稱作實體波(Body Wave)，而另一類稱作表面波(Surface Wave)。實體波指的是在物質內部傳遞的波，由其不同的傳遞特性可分為 P 波及 S 波(如圖 2.1 所示)。P 波又可稱為壓縮波，P 是指最初(Primary)的意思，表示在震波紀錄上最先到達。其在地層中傳遞時，地層介質顆粒的運動方向與傳遞的行進方向平行，亦即若 P 波在 x 方向上傳遞，則介質顆粒會沿著ｘ方向，以原來位置為中心作左右的震盪，而使得介質顆粒有疏密的情形產生(如圖 2.1a 所示，介質顆粒在空間上相鄰最密或最疏處的距離是一個波長的距離)。S 波又可稱為剪力波，S 是指第二(Secondary)的意思，表示在震波紀錄上較 P 波慢到達。其在地層中傳遞時，地層介質顆粒的運動方向與傳遞方向相互垂直，亦即若 S 波在 x 方向上傳遞，則介質顆粒會沿著ｚ方向，以原來位置為中心做上下的震盪，而使得介質顆粒的排列有高下之分如波浪狀(如圖 2.1b 所示，介質顆粒在空間上相鄰等高處的距離是一個波長的距離)。表面波與實體波不同，表面波是存在於具有自由表面存在的介質中，其沿著自由表面行進，且會隨著深度而有能量衰減的現象。表面波依其傳遞特性可分為雷利波(Rayleigh Wave)及勒爾波(Love Wave)(如圖 2.2 所示)。雷利波在傳遞時，地層介質顆粒同時有與傳遞方向垂直及平行的運動，亦即當其向+ｘ方向傳遞時，在地表面處，介質顆粒會以原位置為中心，形成橢圓形的路徑，以逆時鐘方向震動。其所形成的橢圓型平面只會在 x、z 平面上，而 y 方向上並無分量(如圖 2.2a 所示，介質顆粒在空間中相鄰最密或

(21)

圖圖圖

圖 2.2.2.2.1111 實體波示意圖實體波示意圖實體波示意圖實體波示意圖 (Bolt, 1976) (Bolt, 1976) (Bolt, 1976) (Bolt, 1976)

最舒的距離為一個波長)，與同在 x、z 平面上行進的 P 波及 x、z 平面上的 S 波分量(以 SV 表示)，統稱為 P-SV 波。勒爾波與表面波不同，雖亦是在地表面處傳遞，然而其傳遞方向與地層介質顆粒的運動方向互相垂直，且只有在 y 方向上有分量，亦即當其在 x 方向傳遞時，地層介質顆粒會以原位置為中心，沿著水平 y 方向做前後的運動(如圖 2.2b 所示，介質顆粒於空間中在ｙ方向位移相同的最短距離為一個波長)。其單純只有 y 方向的分量，是由地底介面反射回地表的 S 波在 y 方向的分量組成，可以以 SH 波名之。在沒有介面存在的半無限域空間中，地表處不會有勒爾波的存在。在接下來的兩小節中，我們將以量化的方式，以平面波(plane wave)為例來討論實體波與表面波的性質以及這兩種波與拉瑪彈性材料常數(Lame’s constant)λ,µ的關係。

(22)

圖圖圖

圖 2222....2 2 2 表面波示意圖2 表面波示意圖表面波示意圖表面波示意圖 (Bolt, 1976) (Bolt, 1976) (Bolt, 1976) (Bolt, 1976) 2.1.1 壓縮波與剪力波壓縮波與剪力波壓縮波與剪力波壓縮波與剪力波地層在受到擾動後的運動情形是受到運動方程式(equation of motion)以及虎克定律(Hooke’s Law)控制。在三維的空間中，以卡氏座標系為座標(如圖 2.3 所示)，運動方程式可寫做： z y x f t u xx yx _zx x x ∂ ∂ + ∂ ∂ + ∂ ∂ + = ∂ ∂ σ σ σ ρ 2₂ z y x f t u _xy _yy _zy y y ∂ ∂ + ∂ ∂ + ∂ ∂ + = ∂ ∂ σ σ σ ρ ₂ 2 z y x f t u xz yz _zz z z ∂ ∂ + ∂ ∂ + ∂ ∂ + = ∂ ∂ σ σ σ ρ 2₂ (2. 1) 其中，ρ表示地層密度；u_x,u_y,u_z分別表示在x,y,z方向上的位移量；t表示時間；f_x, f_y, f_z分別表示在x,y,z方向上的實體力(Body force)；σ_ij表示作用在法

(23)

線方向為 i 方向平面上 j 方向的應力，且σ_ij =σ_ji。例如σ_xy便是表示作用在 x 平面上 y 方向的力(如圖 2.3 所示)且σ_xy =σ_yx。假設在沒有震源影響或是實體力的作用下，((2. 1)式可改寫為： z y x t ux xx yx _zx ∂ ∂ + ∂ ∂ + ∂ ∂ = ∂ ∂ σ σ σ ρ 2₂ (2. 2a) z y x t uy xy yy zy ∂ ∂ + ∂ ∂ + ∂ ∂ = ∂ ∂ σ σ σ ρ ₂ 2 (2.2b) z y x t uz xz yz _zz ∂ ∂ + ∂ ∂ + ∂ ∂ = ∂ ∂ σ σ σ ρ 2₂ (2.2c) 圖圖圖圖 2.2.2.3333 應力定義示意圖2. 應力定義示意圖應力定義示意圖應力定義示意圖而對於均值均向的彈性材料來說，其應力應變之組合律可以兩個獨立的拉瑪彈性常數λ,µ (Lame’s constants)透過虎克定律表示之，虎克定律的通用式可寫做： ij ij kk ij λε δ µε σ = +2 _{(2. 3)} 其中，δ_ij為 Kronet delta，表示其在i= j時為 1，i≠ j時為 0；λ,µ表示是地層材料的拉瑪彈性常數；ε_ij表示 i 方向上的線段在 j 方向上的應變量，且依

(24)

位移與應變的關係可得 _      ∂ ∂ + ∂ ∂ = i u j ui j ij 2 1 ε 。例如： x u z u y u x ux y _z _x xx _∂ ∂ +       ∂ ∂ + ∂ ∂ + ∂ ∂ =λ µ σ 2 ) ( x u y ux y xy _∂ ∂ + ∂ ∂ =µ σ       ∂ ∂ + ∂ ∂ = x z z x xz u u u u µ σ (2. 4) 將((2. 4)式代入((2. 2a)式中可得：

(

)

_      ∂ ∂ + ∂ ∂ + ∂ ∂ +       ∂ ∂ + ∂ ∂ + ∂ ∂ ∂ ∂ + = ∂ ∂ 2 2 2 2 2 2 2 2 z u y u x u z u y u x u x t ux λ µ x y _z µ _x _x _x ρ _{(2. 5a)} 同理，由((2. 2b),(2.2c)式可得：

(

)

_       ∂ ∂ + ∂ ∂ + ∂ ∂ +       ∂ ∂ + ∂ ∂ + ∂ ∂ ∂ ∂ + = ∂ ∂ 2 2 2 2 2 2 2 2 z u y u x u z u y u x u y t u _y _y _y z y x y µ µ λ ρ (2. 5b)

(

)

_      ∂ ∂ + ∂ ∂ + ∂ ∂ +       ∂ ∂ + ∂ ∂ + ∂ ∂ ∂ ∂ + = ∂ ∂ 2 2 2 2 2 2 2 2 z u y u x u z u y u x u z t uz λ µ x y _z µ _z _z _z ρ _{(2. 5c)} 將((2. 5)式以向量型式表示，可表示為：ｕｕｕｕｕｕｕｕｕｕｕｕ.. 2 ) ( ) 2 ( + ∇ ∇⋅ + ∇ = λ µ µ ρ 　　 (2. 6) 便得到在三維空間中，均質均向線彈性材料的運動方程向量式。一般為了後續計算上的便利會採用向量關係式： ) ( -) ( 2 _ｕｕｕｕｕｕｕｕｕｕｕｕ=∇ ∇⋅ ∇×∇× ∇ 　　 (2. 7) 將其帶入((2. 6)式：

)

(

-)

(

)

2 (

ｕ

..

=

λ

+

µ

∇

⋅

µ

∇

×

∇

×

ρ

(2. 8)

(25)

根據 Helmholtz 理論(Lay and Wallace, 1995)，任何向量場ｕ皆可以純量

潛能φ (scalar potential)與向量潛能

ψ

(vector potential)表示為：

ψ ｕｕｕｕ=∇φ +∇× (2. 9) 其中，

0 =

×

∇

φ

，

∇ ψ

⋅

=

0

(2. 10) 將((2. 9)式帶入((2. 8)式整理，便可獲得 0 ] [ ] ) 2 [( .. 2 .. 2 ₋ ₊_∇_× _∇ ₋ ₌ ∇ + ∇ λ µ φ ρφ µ ψ ρψ (2. 11) 令 ρ µ λ α = +2 ， ρ µ β = (2. 12) 則((2. 11)式在符合下列的情況時， 0 1 .. 2 2 ₋ ₌ ∇ φ α φ ， 1 0 .. 2 2 ₋ ₌ ∇ ψ β ψ _{(2. 13)} 可獲得其解。((2. 13)式即為波動方程式，其中，α、β分別為此彈性體的壓縮波(Ｐ波)與剪力波(Ｓ波)波速。由((2. 13)式可知，對於地層受到擾動後的位移情形，可以將其分為 P 波組以及 S 波組討論，而且地層的實體波波速與頻率無關，只與地層材料的拉瑪彈性常數有關。對於一般地層材料而言，λ與µ的值相差不大，由((2. 12)式可知，α >β ，亦即 P 波的波速較 S 波快，也因此在震測資料上會先看見 P 波造成的位移，而後才有 S 波的位移。接著利用變數分離法(separation of variables)對((2. 13)式波動方程式求解。假設

)

(

)

(

)

(

)

(

)

,

(

x

y

z

t

=

X

x

Y

y

Z

z

T

t

φ

(2. 14) 帶入(2. 13 式可得平面波在任何方向行進時的波動純量潛能φ通解為

(26)

)] ( exp[ ) , , , (x y z t A i ωt k_α_xx k_α_y y k_α_zz φ = ± ± ± ± _{(2. 15)} 其中， 2 2 2 2 2 α ω α α αx + k y + k z = k (2. 16) 定義了於特定頻率ω 下，以 P 波波速α 在卡氏坐標系中朝(k_α_x,k_α_y,k_α_z)方向傳遞的平面波表面，而 A 為常數，表示其振幅；同理可得波動向量潛能

ψ

的通解： )] ( exp[ ) , , , (x y z t = B ±i ωt ± kβxx ± kβyy ± kβzz ψ _{(2. 17)} 其中， 2 2 2 2 2 β ω β β βx + k y + k z = k (2. 18) 定義了於特定頻率ω 下，以 S 波波速β 在卡氏坐標系中朝(kβx,kβy,kβz)方向傳遞的平面波表面，Ｂ為常數，表示其振幅。由地層擾動所造成的位移u ，基於 P 波與 S 波在行進時並不會受到互相的干擾的假設，我們可以以Ｐ波與Ｓ波震動的線性組合表示之。所以根據 Helmholtz 理論，((2. 9)式可寫做：

ｕ

=

ｕ

p

+

ｕ

s =∇

φ

+∇×ψ (2. 19) 其中，u ,pｕｕｕｕｓｓｓｓ分別表示由 P 波及 S 波所造成的位移。又 ∧ ∧ ∧ ∂ ∂ + ∂ ∂ + ∂ ∂ ∇ = z z y y φ φ φ φ ｘｘ＝　ｐｐｐｐｕｕｕｕ _{(2. 20)} ∧ ∧ ∧       ∂ ∂ ∂ ∂ +       ∂ ∂ ∂ ∂ +       ∂ ∂ ∂ ∂ × ∇ = z y -x y x -z x z -y ＝ ψz ψy ψx ψz ψy ψx ψ ｕｕｕｕ_ｓ_ｓ_ｓ_ｓ _{(2. 21)} 其中，ψx,ψy,ψz分別表示ψ 在x,y,z上的分量，符合((2. 17)式通式解的形

(27)

式；∧x,∧y,z∧分別表示x,y,z方向的單位向量。為了簡化問題，假設入射波為一 0 , _y = y k k_α _β 的平面波，亦即此波在 y 方向上並無變化， →0 ∂ ∂ y i ψ _{。則((2. 21} 式可改寫為： ∧ ∧ ∧       ∂ ∂ +       ∂ ∂ − ∂ ∂ +       ∂ ∂ × ∇ z z y z x z ψy ψx ψ ψy x z ＝ ψ _{(2. 22)} 將((2. 20)及((2. 22)式代入((2. 19)式，得 ∧ ∧ ∧       ∂ ∂ + ∂ ∂ +       ∂ ∂ − ∂ ∂ +       ∂ ∂ ∂ ∂ z x y z x z φ ψy ψx ψ φ ψy ｚｘ＝ x z ｕｕｕｕ _{(2. 23)} 由((2. 23)式的結果，我們可以很清楚的看見，在 y 方向的位移只與 S 波有關，而在 x，y 方向上的位移則同時牽涉到 P 波與 S 波的反應，但其 S 波牽涉的與 y 方向牽涉的 S 波不同，所以我們可以很成功的，以適當的座標方式使 P 波的在一方向上，而後便可將擾動後地層的位移分解成 x-z 平面上，以及 y 方向上來討論。y 方向的波動位移以 SH 波表示，x-z 平面上的位移以 P-SV 波表示。 2.1.2 表面波及其頻散曲線表面波及其頻散曲線表面波及其頻散曲線表面波及其頻散曲線表面波於半無限域地層表面波於半無限域地層表面波於半無限域地層表面波於半無限域地層：表面波的存在需有地表自由表面，其特徵是傳遞能量的大小會隨著深度的增加而遞減。在此我們先將探討在半無限域地層中，表面波存在的情形。由上一小節的結果可知，各個方向上的位移u_x,u_y,u_z 可表示為： z x u_x y ∂ ∂ − ∂ ∂ = φ ψ (2. 24a) x z u x z y _∂ ∂ − ∂ ∂ = ψ ψ (2. 24b)

(28)

x z u_z y ∂ ∂ + ∂ ∂ = φ ψ (2. 24c) 然而，由於u_y可被獨立分開來，其位移僅牽涉到 S 波的向量潛能，是否以向量潛能的方式計算並不影響其通式的形式，因此將其保留為u_y，而不需特以向量潛能表示。由得到實體波通式的形式經驗，並已知其在 z 方向上有能量遞減的特性，我們對於一個在+ x 方向上具有頻率ω，以波速 c 傳遞的波，可假設其φ,ψ_y,u_y分別為：

( )

z

[

ik

(

x ct

)

]

f − = exp φ (2. 25a)

( )

z

[

ik

(

x ct

)

]

g y = exp − ψ _{(2. 25b)}

( )

z

[

ik

(

x ct

)

]

h u_y = exp − (2. 25c) 其中，f(z)、g(z)、h(z)用以表示其隨著地層深度變化而有能量大小變化的現象； c k =ω 表示空間頻率。將式((2. 25a)，((2. 25b)式代入((2. 13)式，可得 0 ) ( ) ( 2 2 '' ₊ ₌ z f r k z f (2. 26a) 0 ) ( ) ( 2 2 '' ₊ ₌ z f s k z g (2. 26b) 其中， 2 / 1 2 2 1      − =

α

c r (2. 27a) 2 / 1 2 2 1      − =

β

c s (2. 27b) 而對於u_y，若以假設其向量潛能的方式代入((2. 13)式中，亦可發現 0 ) ( ) ( 2 2 '' ₊ ₌ z h s k z h (2. 26c) 將((2. 26)式求解，可得 ) exp( ) exp( ' ) (z A ikrz A ikrz f = + − (2. 28a)

(29)

) exp( ) exp( ' ) (z B iksz B iksz g = + − (2. 28b) ) exp( ) exp( ' ) (z C iksz C iksz h = + − (2. 28c) 其中，

A

,

A

'

,

B

,

B

'

,

C

,

C

'

為任意常數，但由於表面波有隨著深度而能量傳遞遞減的行為，所以

A

'

=

B

'

=

C

'

=

0

；且r 為虛數，s為實數，亦即表示 c > >β α 。將((2. 28)式代回((2. 25)式，整理後得到：

(

)

[

ikrz

ik

x

ct

]

A

−

+

−

= exp

φ

(2. 29a)

(

)

[

iksz

ik

x

ct

]

B

y

= exp

−

+

−

ψ

(2. 29b)

(

)

[

iksz

ik

x

ct

]

C

u

_y

= exp

−

+

−

(2. 29c) 為了要將Ａ、Ｂ、Ｃ的常數值找出來，我們需要給予邊界條件。在自由表面處，在ｚ方向上的應力會為零，亦即

σ

_xz

=

σ

_yz

=

σ

_zz

=

0

。採用((2. 3)式及位移與應變的關係，我們可獲得邊界條件： 0 =       ∂ ∂ + ∂ ∂ = z u x u_z _x zx

µ

σ

(2. 30a) 0 =       ∂ ∂ + ∂ ∂ = y u z u_y _z zy

µ

σ

(2. 30b)

(

2

)

_=0      ∂ ∂ + ∂ ∂ + ∂ ∂ + = y u x u z u_z _x y zz

λ

µ

λ

σ

(2. 30c) 將((2. 24)式代入((2. 30)式，可得 0 2 ₂ 2 2 2 2 = ∂ ∂ − ∂ ∂ + ∂ ∂ ∂ z x x z y y ψ ψ φ (2. 31a) 0 = ∂ ∂ z u_y (2. 31b)

(

2

)

2 0 2 2 2 2 = ∂ ∂ ∂ + ∂ ∂ + ∂ ∂ ∂ + z x x z z y ψ µ φ λ φ µ λ (2. 31c)

(30)

先將((2. 29c)式代入((2. 31b)可得 C = 0。由此便可知，勒爾波在半無限域地層中是不存在的，其需要地層有介面的存在時，才會成形。再將((2. 29a)， ((2. 29b)式分別代入((2. 31a)，((2. 31c)式中，將可獲得一組聯立方程式：

(

1

)

0 2rA− −s2 B= (2. 32a)

(

)

[

α2 r2 +1 −2β2

]

A−2β2sB =0 (2. 32b) 則除了 A = B = 0 的這組解之外，欲使此聯立方程式有解，必須

(

)

(

)

0 2 2 1 1 2 2 2 2 2 2 = − − + − − s r s r β β α (2. 33) 解出此行列式，得到

(

)

[

α2 r2 +1 −2β2

]

(

1−s2

)

−4rsβ2 =0 (2. 34) 將((2. 27)式代入((2. 34)式，整理後可得 2 / 1 2 2 2 / 1 2 2 2 2 2 1 1 4 2 _      −       − =       − β α β c c c (2. 35) ((2. 34)式便是在 1887 年首先由 Rayleigh 所提出的表面波波速解方程式，在此式中，並沒有頻率的出現，表示了在半無限域地層中，雷利波波速 c 與頻率無關，亦即沒有頻散現象的發生。為了更加了解此一方程式的特性，令 2 2 , _      =       = α β β ξ c q (2. 36) 代入((2. 35)式，且將等式兩邊平方後，去除ξ =0的解，可得：

(

3 2

)

16

(

1

)

0 8 8 2 3 ₋ ₊ ₋ ₊ ₋ ₌ q qξ ξ ξ (2. 37) 在此式中，當半無限域地層的彈性參數確定後， q 便是常數，ξ為未知數， ((2. 37)式為一元三次方程式，可解得三個解。由於雷利波的波速 c 必符合 c > >β α 的要求，因此只有符合ξ <1的解，才是雷利波波速的解。

(31)

透過((2. 36)式的變數變換，可知在半無限域空間中，雷利波的波速與地層的卜松比(Poisson’s Ratio)有關，如圖 2.4 及表 2.1 所示，q 與卜松比有一對一的關係。在可能的卜松比範圍(0 ~ 0.5)中，ξ在 0.7640 至 0.9128 之間變動，而由((2. 36)式，可得到雷利波與 S 波波速的關係。在這要特別提出， 9128 . 0 = ξ 是對應於 0.5 的卜松比，亦即 S 波波速為 0 的狀況，其雷利波波圖圖圖圖 2.42.42.4 2.4 卜松比與半無限域雷利波速關係圖卜松比與半無限域雷利波速關係圖卜松比與半無限域雷利波速關係圖卜松比與半無限域雷利波速關係圖表表表表 2222....1 1 1 卜松比與半1 卜松比與半卜松比與半卜松比與半無限域地層中雷利波波速對照表無限域地層中雷利波波速對照表無限域地層中雷利波波速對照表無限域地層中雷利波波速對照表**** 卜松比

ν

**

q

ξ

_{雷利波波速}

c

_R 0 ₂1 0.7640 0.8741β 8 1 7 3 0.8059 0.8977β 4 1 3 1 0.8453 0.9194β 2 1 0 0.9128 0 * Udias, 1999 **

(

q

)

q − − = 1 2 2 1 ν 速亦為 0，由此亦可知，雷利波是不存在於純彈性液體的半無限域空間中

(32)

(Udias, 1999)。雷利雷利雷利雷利波於波於波於水平層狀波於水平層狀水平層狀水平層狀地層地層地層地層：：：：本研究中採用的表面波震測法，其所收集分析的波實際是指在 P-SV 波中的雷利波。雷利波很重要的一個特性便是其因地層隨身度變化而產生頻散(dispersion)現象，亦即雷利波波速會因頻率的不同而有所不同，然而由半無限域地層的雷利波波速解可知，在半無限空間中，雷利波的波速固定。為表達雷利波的頻散現象，接下來將討論在均質均向的層狀水平地層(後續將以水平層狀地層表示之)中，雷利波波速與頻率的關係。假設一水平層狀地層(如圖 2.5 所示)，地層的參數只與 z 有關。在接下來的推導中，將只討論 P-SV 波的部份，因此在 y 方向的位移分量uy將令其為 0 則一個在+x 方向上傳遞的平面波可表示為：

(

)

[

ik x ct

]

c z k U u_x = ( , , )exp − _{(2. 38)}

(

)

[

ik x ct

]

c z k iW u_z = ( , , )exp − (2. 39) 其中，U(k,z,c),W(k,z,c)分別表示在 x 及 z 方向位移量的振幅大小，其具有隨著 z 增加而遞減，隨著不同空間頻率 k 以及雷利波波速 c 而改變的特性，為求便易，後續將以U(z),W(z)表示；此外，

c,

k

相乘便是頻率ω。由於在層與層之間，在ｚ面上的應力有連續的現象，因此令

(

)

[

ik x ct

]

c z k X xz = ( , , )exp − σ (2. 40)

(

)

[

ik x ct

]

c z k iZ zz = ( , , )exp − σ (2. 41) 其中，X(k,z,c),Z(k,z,c)分別表示在 z 平面上，x 及 z 方向應力的大小，其具有隨著 z 增加而遞減，隨著不同空間頻率 k 以及雷利波波速 c 而改變的特性，為求便易，後續將以X(z),Z(z)表示。

(33)

圖圖圖圖 2.2.2.5555 水平層狀地層模型2. 水平層狀地層模型水平層狀地層模型水平層狀地層模型採用((2. 3)式及位移與應變的關係，我們可獲得應力與位移的關係為：       ∂ ∂ + ∂ ∂ = z u x u z z x zx µ( ) σ (2. 42a)

(

)

      ∂ ∂ + ∂ ∂ + = x u z u z z z x zz λ µ λ σ ( ) 2 ( ) (2. 42b) 將((2. 38), ((2. 39)式代入((2. 42)式中，

(

)

[

]

[

(

)

]

     ₋ ₋ ₊ ₋ = ik x ct dz z dU ct x ik k kW z zx exp ) ( exp ) ( ) ( µ σ (2. 43)

(

)

exp

[

ik

(

x-ct

)

]

i kU(z)exp

[

ik

(

x-ct

)

]

dz dW(z) ) ( 2 ) ( i λ µ λ σ_zz = z + z + (2. 44) 再將((2. 40), ((2. 41)式代入((2. 43), ((2. 44)式中，整理後可得 ) ( ) ( 1 ) ( z kW z X dz z dU + = µ (2. 45)

(34)

) ( 2 ) ( 2 1 ) ( z U k z Z z z dW µ λ λ µ λ+ − + = _{(2. 46)} 又，前述所假設的u ,x uz必需符合運動方程式，由((2. 5a), ((2. 5c)式可知其須符合：

(

)

_      ∂ ∂ + ∂ ∂ +       ∂ ∂ + ∂ ∂ ∂ ∂ + = ∂ ∂ 2 2 2 2 2 2 ) ( ) ( ) ( ) ( z u x u z z u x u x z z t u z x λ µ x z µ x x ρ _{(2. 47)}

(

)

_      ∂ ∂ + ∂ ∂ +       ∂ ∂ + ∂ ∂ ∂ ∂ + = ∂ ∂ 2 2 2 2 2 2 ) ( ) ( ) ( z u x u z u x u z z z t u z z λ µ x z µ z z ρ _{(2. 48)} 將((2. 15), ((2. 16)式代入，

(

)

[

−

]

= −ρ(z)U(z)c2k2expik x ct

(

)

[

ik

(

x ct

)

]

dz dW z kU z z k  −      ₊ + − λ( ) µ( ) ( ) exp

(

)

[

ik x ct

]

dz z U d z U k z _ −      − − ( ) ( ) ₂( ) exp 2 2 µ (2. 49)

(

)

[

−

]

= −iρ(z)W(z)c2k2expik x ct

(

)

[

ik

(

x ct

)

]

dz z W d dz z dU k z z i _ −      + +2 ( ) ( ) ( ) exp ) ( ₂ 2 µ λ

(

)

[

ik x ct

]

dz W d z W k z i _ −      + + ( ) ( ) ₂ exp 2 2 µ (2. 50) 再將((2. 45), ((2. 46)式代入((2. 49), ((2. 50)式，整理過後可得

[

]

dz z dW z k dz z dX z U z z c z k2 −ρ( ) 2 +λ( )+2µ( ) ( )= ( )− λ( ) ( ) (2. 51)

(35)

[

]

dz z dZ dz z dU z k z W z c z k2 −ρ( ) 2+µ( ) ( )= µ( ) ( )+ ( ) (2. 52) 最後將((2. 46)式代入((2. 51)式且將((2. 45)式代入((2. 52)式，整理後可得：

(

)

) ( ) ( 2 ) ( ) ( k U(z) ) ( 2 ) ( ) ( ) ( ) ( 4 ) ( ) ( 2 2 z Z z z z z z z z z c z k dz z dX µ λ λ µ λ µ λ µ ρ + +       + + + − = _{(2. 53)} ) ( ) ( ) ( ) ( 2 2 z kX z W k c z dz z dZ − − = ρ (2. 54) 由((2. 45), ((2. 46), ((2. 53), ((2. 54)式，可獲得運動-應力向量(motion-stress

vector, Aki and Richard, 2002) [U(z),W(z),X(z),Z(z)]的一階微分方程式

                                − − + − + + − =             ) ( ) ( ) ( ) ( 0 ) ( 0 ) ( 2 ) ( ) ( 0 0 ) ( ) ( ) ( 2 ) ( 1 0 0 ) ( 2 ) ( ) ( 0 ) ( 1 0 ) ( ) ( ) ( ) ( 2 2 2 2 2 z Z z X z W z U k k c z z z z k k c z k z z z z z k z z k z Z z X z W z U dz d ρ µ λ λ ρ ζ µ λ µ λ λ µ (2. 55) 其中，

[

]

) ( 2 ) ( ) ( ) ( ) ( 4 ) ( z z z z z z µ λ µ λ µ ζ + + = _{(2. 56)} 傳遞矩陣傳遞矩陣傳遞矩陣傳遞矩陣(Propagator Matrix)：欲解出((2. 55)式的方法有許多，包括了數值積分法(numerical integration method, Takeuchi and Saito, 1972)以及傳遞矩陣法(Propagation matix method, Thomson, 1950；Haskell, 1953；Gilbert and Backus, 1966)。目前較為使用的是傳遞矩陣法，其能符合我們將地層假設為水平層狀的情形。利用矩陣式求解的方法是由 Thomson(1950)先行提出，經由 Haskell(1953)對其中錯誤改正，最後在 1966 年由 Gilbert and Backus 將傳遞矩陣的方法引入，而得知 Thomson-Haskell 的矩陣求解方式是傳遞矩陣法的一個特例情況。以下將介紹傳遞矩陣的求解方式。

(36)

在求解前，為了後續參數使用的方便，採用((2. 12)式將((2. 55)式改為：                                   − −       −       − − =             ) ( ) ( ) ( ) ( 0 ) ( 0 ) ( ) ( 2 1 0 0 ) ( ) ( ) ( ) ( 1 0 0 ) ( ) ( 2 1 0 ) ( ) ( 1 0 ) ( ) ( ) ( ) ( 2 2 2 2 2 2 2 2 2 2 z Z z X z W z U k k c z z z k z z c k z z z z k z z k z Z z X z W z U dz d ρ α β ϖ ρ α ρ α β β ρ (2. 57) 其中， 1 ) ( ) ( 1 ) ( 4 ) ( ₂ 2 2 2 −       − = z z c z z α β β ϖ _{(2. 58)} 令

[

]

T Z(z) X(z) W(z) U(z) ) (z = b ，則((2. 57)式可改寫為 ) ( ) ( ) ( z z A dz z d b b = (2. 59) 其中，                     − −       −       − − = 0 ) ( 0 ) ( ) ( 2 1 0 0 ) ( ) ( ) ( ) ( 1 0 0 ) ( ) ( 2 1 0 ) ( ) ( 1 0 ) ( 2 2 2 2 2 2 2 2 2 k k c z z z k z h z c k z z z z k z z k z A ρ α β ρ α ρ α β β ρ (2. 60) 定義傳遞矩陣為

∫

+

∫

+ + = z z z 1 2 2 z z 1 1 1 0 0 1 0 0 ... d )d A( ) A( )d A( I ) , ( ζ ζ ζ ζ ζ ζ ζ z z P _{(2. 61)} 其中，Ｉ為 4x4 的單位矩陣；zref表示為深度參考點；ζi為表示深度的不同變數。由此定義，可以清楚的發現，其符合((2. 59)式： ) , ( ) ( ) , ( ref ref z z P z A dz z z dP = (2. 62) 假設 z 是在第 i 層地層中zi ~ zi+1變化的深度變數，則Ａ將與深度無關，

(37)

成為常數方陣，記為Ai。由此假設，此層中的傳遞矩陣可簡化為

(

)

(z-z ) A A ... 2 1 A I ) ,

(z z_ref = + z−z_ref _i + _ref 2 _i _i + P =exp

[

(

z-zref

)

Ai

]

(2. 63) 根據 Slyvester 的展開式(Hildebrand, 1952)，若一方陣 A 具有 n 個不同值的 eigenvaluesλk，則此方陣的函數P可表示為：

(

)

(

)

∑

_∏

∏

= ≠ ≠ = n 1 k k r r k k r r -I -A ) ( ) A ( λ λ λ λ_k P P (2. 64) 對Ai求其 eigenvaluesλk，即是求

0 I

-A

_i

λ

=

(2. 65) 的解，可得 1 ₂ 2 i k c k α λ =± − 及 1 ₂ 2 i c k β − ± ，代入(2. 64 式，便可獲得每一層的傳遞矩陣(見附錄一)。由傳遞矩陣的定義((2. 61)式可知： I ) , (z_ref z_ref = P _{(2. 66)} 亦即對於在任何深度 z，以zref為深度參考點的傳遞矩陣可寫做 ) , ( ) , ( ) ,

(z z_ref P z z_ref P z_ref z_ref

P = _{(2. 67)} 假設b(z0)為((2. 59)式以z0為深度參考的解，則 ) ( ) , ( ) (z P z z₀ b z₀ b = _{(2. 68)} 為((2. 59)式的另一個解。在此我們將可得到傳遞矩陣的兩個特性：其一，使用((2. 68)式的關係，則在z0至zn間被ｎ個不連續界面分開的水平層狀地層模型(如圖 2.5 所示)，其解可寫做 b(z_n)=P(z_n,z₀)b(z₀) =P(z_n,z_n₋₁)P(z_n₋₁,z₀)b(z₀) (2. 69)

(38)

=P(z_n,z_n₋₁)P(z_n₋₁,z_n₋₂)...P(z₁,z₀)b(z₀) 亦即，由自由表面開始至半無限域空間的全域傳遞矩陣可以以每層地層所得到的區域傳遞矩陣相乘獲得。其二，傳遞矩陣具有 ) , ( ) , ( ₁ ₀ ₀ ₁ 1 z z P z z P− = (2. 70) 的特性。由((2. 69)式可知，n=2 時，其解可寫做 ) ( ) , ( ) (z₂ P z₂ z₁ b z₁ b = = P(z₂,z₁)P(z₁,z₀)b(z₀) (2. 71) 若令z2 =z0，則 ) , ( ) , ( I= P z₀ z₁ P z₁ z₀ _{(2. 72)} 便得到如((2. 70)式的第二特性。為了對傳遞矩陣法的物理意義更加顯著，在這再由 Thomson-Haskell 方法切入做討論。前已提及 Thomson-Haskell 方法是傳遞矩陣中的ㄧ個特殊解，假設

(

)

[

k ref

]

kexp z-z v (z)= λ b (2. 73) 為((2. 59)式的解，其中vk為對應於 A 的 eigenvalueλk的 eigenvector。若令一個矩陣Ｂ，其欄位的位置是((2. 73)式形式的不同解，則 4x4 的常數矩陣 A 就有四個λk及vk，此矩陣Ｂ也就同為 4x4 的方陣。((2. 59)式的通解便為((2. 73 式形式的解線性疊加，則其解 b 可寫做 a b(z)=B(z-z_ref) (2. 74) 其中，a 為 4x1 的向量，表示不同解疊加時的權重係數。已知Ｂ為((2. 73) 式形式的解的收集，將((2. 65)式所得到的解λk代入((2. 73)式，則我們可以將其與 z 相關的部分與非 z 相關的部分分離出，得到 ) z -QE(z ) z -B(z _ref = _ref (2. 75) 其中，

(39)

            − − − = s s t t t t r r r r s s µ µ µ µ µ µ µ µ 2 2 2 2 1 1 1 1 Q _{(2. 76)}

[

]

(

)

[

]

(

)

[

]

(

)

[

]

             − − − − − − = z z ks z z ks z z kr z z kr ref ref ref ref exp 0 0 0 0 exp 0 0 0 0 exp 0 0 0 0 ) ( exp ) z -E(z _ref _{(2. 77)} 2 2 c -2 t β = (2. 78) _ 2 2 2 2 _ 2 2 2 2 s i 1 -c i c -1 s i 1 c i c -1 r ) (c ) c ( = = = − = > > < < β β α α β α β α r c c (2. 79)

藉由波的特性給予((2. 74)式邊界條件(詳細推導請見 Aki and Richard, 2002)，可以解得 T PA     = PA↓ ↑ SA↓ SA↑ a (2. 80) 其中， _PA↓ PA↑ SA↓ SA↑ _分別表示 P 波向下傳遞、向上傳遞、S 波向下傳

遞及向上傳遞的量度 (amplitude) ，可稱其為量度向量 (amplitude vector,

Buchen and Ben-Hador, 1996)。由上述的式子可知，Ｑ與ａａａａ僅與地層參數有

關，亦即在同一地層中，不論在地層中的任何位置，Ｑ與ａａａ皆相同，因此ａ在水平層狀的地層中，第ｉ層的下部介面zi處會滿足， i a b(z_i)=B(z_i -z_i_-₁) (2. 81) 而在其上部介面zi−1處會滿足

(40)

i a b(z_i_-₁)=B(z_i_-₁-z_i_-₁) (2. 82) 將((2. 76)式代入((2. 82)式中可得 i i 1 -i

)

Q

(z

a

b

=

(2. 83) 整理後可得

)

(z

Q

_i-1 _i_-₁ i

b

a

=

(2. 84) 代入((2. 81)式 ) (z )Q (-d E Q ) ( i-1 -1 i i i i b b zi = (2. 85) 其中，d_i =z_i −z_i₋₁，即該層厚度。又由((2. 69)式知 ) (z ) z , (z ) ( _i _i_-₁ b _i_-₁ b z_i =P (2. 86) 因而得到傳遞矩陣的定義便是 -1 i i i i 1) Q E (-d )Q , (zi zi− = P (2. 87) 這所表示的傳遞矩陣是指，由zi往zi−1方向傳遞的傳遞矩陣，而若對於由zi−1 往zi方向傳遞，則由((2. 70)式的關係可得 -1 i i i i 1, ) Q E (d )Q (z_i₋ z_i = P (2. 88) 從上述式子，已可了解，傳遞矩陣在連結同一地層中不同深度的運動-應力向量是將各地層中的不變量 Q 以 E 矩陣依著其由參考深度開始的傳遞距離及方向做相位的變化後，而重新得到的一個矩陣量。也就是，傳遞矩陣最後是將每個地層的地層資訊以及波在傳遞時的相位變化總和起來，而得到整個地層的訊息。下段將代入邊界條件進而說明這一概念。雷利波的頻散曲線雷利波的頻散曲線雷利波的頻散曲線雷利波的頻散曲線以水平層狀地層的模型來說，P-SV 波須符合三個邊界條件：1) 在自由表面處其在 z 面上的應力為零。2) 層與層的介面處，其運動-應力向量連續。3) 在最底層層面以下，波不會有向上傳遞的量。亦即， 0 (0) ) 0 ( ₁ 1 = Z = X (2. 89a) (0) ) (d_i _i ₁ i = b + b (2. 88b)

(41)

0 ) 0 ( SA (0) PAn 1 = +1 = ↑ + ↑ n (2. 88c) 由((2. 74)式知 i i i

(0)

Q

a

b

=

(2. 90) i i i i i i

(d

)

Q

E

(-d

)

a

b

=

(2. 91) 將((2. 91)式重新排列可得

)

(d

(d)Q

E

_i _i _i _i i

b

a

=

(2. 92) 代入((2. 90)式

)

d

(

)Q

(d

E

Q

(0)

_i _i _i -1_i _i _i i

b

=

(2. 93) 將邊界條件((2. 89)式代入

0)

(

)Q

(d

E

Q

(0)

_i _i _i _i-1 _i ₁ i

=

b

+

b

(2. 94) 由此可發現，層與層之間的關係式已經得到，又由((2. 88)式，已可清楚的了解傳遞矩陣扮演的角色(後續使用上，((2. 88)式中的傳遞矩陣將記為P_i)。利用這樣傳遞的關係，最後可得到 (0) ... (0) ₁ ₂ ₃ _n ₁ 1 = b + b PPP P_n (2. 95) 代入((2. 90)式， 1 n 1 n 3 2 1 1(0)= ... Q +a + b PPP P_n (2. 96) 此式便是 Thomson-Haskell 方法所得到的關係式。令b1(0)=Tan+1，並將邊界條件((2. 89a), ((2. 89c)式代入，                           =             + ↓ + ↓ 0 SA 0 PA T T T T T T T T T T T T T T T T 0 0 (0) W (0) U 1 n 1 n 44 43 42 41 34 33 32 31 24 23 22 21 14 13 12 11 1 1 (2. 97)

(42)

則在PA 1 ,SAn+1 ↓ + ↓ n 皆不為零的情況下，欲求((2. 97)式的 eigenvalue 需由 0 T T T T 43 41 33 31 ₌ (2. 98) 獲得，滿足此一式的ｃ與ｋ便是雷利波頻散曲線的位置。為使此式更便於表示頻散曲線方程式，令P=P1P2P3...Pn，則((2. 98)式可寫做 0 PV U' k) D(c, = = (2. 99) 其中，       = 1 0 0 0 0 1 0 0 U' (2. 100)             = = + + + + + + + + + + + 1 n 1 n 1 n 1 n 1 n 1 n 1 n 1 n 1 n 1 n 13 1 n s 2 t t r 2 1 r s 1 e Q V

µ

(2. 101) 如此一來，便得到隱性解表示的雷利波頻散曲線，亦可知道，在水平層狀的地層中，雷利波的波速是與頻率(ω =ck)相關的。在求解上，雖說以數值計算設備進行試誤法找出符合((2. 99)式的 c,k 值，便可獲得其頻散曲線，然而，若僅是使用((2. 99)式的計算方式，當頻率(或空間頻率)變大，很容易會有計算不穩定的情況發生(如圖2.6 所示，黑白的交界處為頻散曲線的位置，對應於其所使用的地層參數，約在頻率1Hz 處便發生不穩定的現象)。造成其不穩定的原因，主要是因為Ｅ矩陣中有指數項的存在，當以數值計算設備進行計算時，因為數值計算設備在計算上有位數的限制，指數項的相乘將使得其值產生計算上的誤差，而造成錯誤。

(43)

Knopoff, 1964；Dunkin, 1965；Watson, 1970；Schwab, 1970；Schwab and Knopoff, 1972；Kennett, 1974；Abo-Zena, 1979；Kennett and Kerry, 1979)，

目前對於水平層狀地層的計算，在使用上以Delta Matrix法( Dunkin, 1965)

以及 Schwab and Knopoff ( Knopoff, 1964；Schwab, 1970；Schwab and Knopoff, 1972)法較為熱門。這部份的討論已超出本研究範圍，不予贅述，

可逕行參考 Schwab and Knopff (1972)、Dunkin(1965)或是 Buchen and

Ben-Hador(1996)的文章。在此須再一提，由Kennett 團隊所發展的RT Matrix

法( Kennett, 1974；Kennett and Kerry, 1979)，將傳遞矩陣與反射傳遞係數的關係做出連結，說明了表面波是由許多不同的反射或折射波到達地面後被鎖在地表的反應。在計算水平層狀的地層時，以RT Matrix法做計算，其計算速度上的表現較不出色，然而，當目前開始考慮地層的黏滯性以及非均向性的水平層狀地層時，其在時間上的優異表現將會凸顯( Buchen and Ben-Hardor, 1996)。圖圖圖

圖 2222....6666 多多多多震態雷利波頻散曲線震態雷利波頻散曲線震態雷利波頻散曲線震態雷利波頻散曲線－－Thomson－－ThomsonThomson----HaskellThomson HaskellHaskellHaskell 法法法法

如圖2.7所示，雷利波在同一頻率下有可能會以多種波速在行進，而使

(44)

模態(mode)，一般稱最下邊的這條頻散曲線震動模態為基態(fundmental mode)，而其餘的為高次態(higher mode)或是過頻(overtone)。如前所述，傳遞矩陣會攜帶不同地層參數的訊息在其中，因此，不同的地層參數，亦會有不同的雷利波頻散曲線存在。目前，藉由此一特性，以及其與剪力波的高相關性，我們可以透過反算的方式得到其相對應的地層剪力波數，而使得表面波震測法得以使用。只是，目前的反算只考慮基態的影響，而忽略高次態的貢獻，這部份擁有極好的發展空間，可使得反算結果更為可信。圖圖圖

圖 2222....7777 多震態雷利波頻散曲線多震態雷利波頻散曲線多震態雷利波頻散曲線多震態雷利波頻散曲線－－－Delta Matrix－Delta MatrixDelta Matrix 法Delta Matrix法法法

雷利波的頻散現象上有一值得提及之處。以純彈性的地層假設而言，雷利波的傳遞並不受到震源能量分布的影響，震源對雷利波的影響僅在於，所產生的震源是否存在有涵蓋目標頻散曲線中所需頻率的能量，因此，雖然上述的推導過程中，僅假設了平面波的波傳形式，但，即使是點震源亦會有相同的結果(Wang, 1981)。在這並不深入推導，可逕行參考列於 2.4 節中的參考文獻，或是 Wang(1981)及 Kennett(1985)。但應注意，在實際的地層中，地層材料並非是全彈性，其會因實際狀況而具有黏滯性或非均向性的行為，而與全彈性的假設結果有所差異。

(45)

2.2 表面波量測方法及其施測問題

表面波量測方法及其施測問題

2.2.1 表面波量測方法表面波量測方法表面波量測方法表面波量測方法

如前所述，表面波震測法主要可分為三個步驟：先至1) 現地施測，而

後對施測所得的地層震動訊號，進行 2) 頻散曲線分析，最後透過 3) 地層

反算，而得到地層的剪力波速剖面。由於在現地施測時所採用接收器數目

的不同，可將表面波震測法分為波譜分析法(Spectal Analysis of Surface Wave,

簡用 SASW)以及多頻道表面波震測法(Multi-channel Analysis of Surface Wave Method, 簡用 MASW)，又由於MASW 提供了較多空間上的訊息，使

得其頻散曲線的分析上有不同的方法，根據採用的方法不同，又可將MASW

分為多頻道表面波譜法(Multi-channel Spectral Analysis of Surface Wave

Method)以及多頻道波場轉換法(Multi-channel Wave Transform of Surface Wave Method)。 2.2.1.1 SASW 現地施測現地施測現地施測現地施測：SASW 在現地的施測上(如圖2.8所示)僅採用兩個接收器，以一圖圖圖

圖 2.82.82.8 2.8 SASWSASWSASWSASW 現地施測示意圖現地施測示意圖現地施測示意圖現地施測示意圖

(46)

與兩個接收器相隔間距∆x相同的距離，即x₀ =∆x。)，左敲、右敲、紀錄後改變接收器間距，使得兩個接收器的中點位置固定在同一個地方，而後左敲、右敲、紀錄，再繼續改變兩個接收器的間距，重複紀錄的動作。頻散曲線分析頻散曲線分析頻散曲線分析頻散曲線分析：透過現地施測，每一筆紀錄我們都將可收集到兩個接收器位置的地層震動訊號，此訊號可表示為位置與時間的函數u₁(x₁,t) ,u₂(x₂,t)。透過傅利葉轉換(Fourier Transform)可得到此二函數在頻率域的函數 ) , ( , ) , ( ₁ ₂ ₂ 1 x ω U x ω U 以及各頻率在x₁ x , ₂處的相位角。根據相位速度的計算式 x ) ( ) c( ∆ ∆ =

_φ

ω

_ω

ω

(2. 102) 其中，c(ω)表示不同頻率時的雷利波波速；∆φ(ω)表示不同頻率時在x₁ x , ₂處相位角的相減；∆x表示x₂ −x₁。將不同頻率所計算得到的雷利波相位速度畫上，便可得到頻散曲線。在計算過程中，相位差∆φ可直接由u₁(x₁,t) ,u₂(x₂,t)

互能頻譜(Cross Spectral Density)的相位角獲得，即

))) , ( ), , ( ( ( ) ( = Angel CSD u₁ x₁ t u₂ x₂ t ∆

φ

ω

(2. 103) 如圖 2.9(虛線)所示，由於頻譜分析所得之相位角差僅侷限在−π ~π之間，因此在計算((2. 102)式前須先將其摺開(unwrap)，得到如圖 2.9 (實線)的結果。由圖2.9中亦可看見，在所得的頻譜分析結果中會有一段資料不佳段，不同的接收器間距會有不同位置及寬度的資料不佳段，因此在施測時收錄多筆不同間距的資料，分別計算完得到頻散曲線後取其算術平均數以為最後頻散曲線結果(如圖 2.10 所示)。最後藉由反算技術，獲得以接收器終點為代表值的地層剪力波速剖面(如圖2.11所示)。

(47)

圖2.9 SASW 頻散曲線分析(相位差-頻率圖)

圖2.10 SASW頻散曲線分析結果

(48)

2.2.1.2 MASW 現地施測

現地施測現地施測

現地施測：多頻道表面波震測法(Multi-channel Analysis of Surface Wave

Method，簡稱MASW法)，首先是由地球物理領域之學者提出(McMechan and Yedlin, 1981; Gabriels et al.,1978；Park et al., 1999)，其在現地施測時，採用

多個在同一直線上的接受器(如圖 2.12 所示)，只需敲擊一次，便可完成。在實際操作上，一般採用1~2公尺之受波器間距，並設置12個以上之接收器接收震源所發出之震波訊號。以24個受波器為例，在第一個受波器之線外取適當近站支距，反覆在同一震源處施作，將其疊加以消除雜訊之影響，直至收錄到清晰之表面波訊號為止。圖圖圖

圖 2.2.2.122.1212 12 MMASWMMASWASWASW 現地施測現地施測現地施測現地施測示意圖示意圖示意圖示意圖頻散曲線分析頻散曲線分析頻散曲線分析頻散曲線分析：多個接收器的收錄，提供了更多空間上的訊息，使得多頻道表面波震測法在頻散曲線的分析上可利用二維訊號處理技術分析震測資料，求得訊號品質優良之頻散關係曲線。根據不同的訊號處理技術，多頻道式表面波量測法可再分為多頻道波譜分析以及多頻道波場轉換法兩種。 1.多頻道波譜分析法多頻道波譜分析法多頻道波譜分析法多頻道波譜分析法：多頻道波譜分析法與波譜分析法的頻散曲線分析類似，經由將多處接收器位置之震測資料u_i(x_i,t)透過傅利葉轉換後得到其在頻率域的函數U_i(x_i,ω)，以及各頻率在x_i處的相位角。如此一來，對於每一個特定頻率，我們皆有其在不同接收器位置的相位角(如圖 2.13 所示)，將頻譜分析所得之相位角摺開後，由((2. 102)式可知，由此斜率便可求得此特定

(49)

頻率的相位波速。在圖中亦可看見，線段中會有資料不佳的段落，依據不同的頻率資料不佳的段落出現位置會不同，一般來說，低頻時會發生於距

離震源較近的位置，稱作近場效應(near field effect)，而高頻時會發生於距

離震源較遠的位置，稱作遠場效應(far field effect) 。分析過程中應當要注意

此兩個效應所可能造成頻散曲線分析摺開時的錯誤(此種效應的影響在 SASW的分析中不易發現，使得SASW的分析有潛在的錯誤)，於計算相位速度相關之斜率∆φ /∆x時，可避開資料不佳段。對不同頻率進行相同的動作，便可得到此地層表面波之頻散曲線(如圖 2.14 所示)。最後藉由反算技術便可獲得以測線中點為代表值的地層剪力波速剖面。圖圖圖圖 2.12.12.13333 2.1 多頻道波譜法多頻道波譜法多頻道波譜法多頻道波譜法頻散曲線分析頻散曲線分析((((相位頻散曲線分析頻散曲線分析相位相位相位----空間位置圖空間位置圖空間位置圖空間位置圖)))) 圖圖圖圖 2.12.12.12.14444 多頻道波譜法多頻道波譜法多頻道波譜法多頻道波譜法頻頻頻頻散曲線分析結果散曲線分析結果散曲線分析結果散曲線分析結果

(50)

2.多頻道波場轉換法多頻道波場轉換法多頻道波場轉換法：多頻道波場轉換法常見於地球物理領域之濾波處理，多頻道波場轉換法

包括頻率波數轉換法（Frequency-Wavenumber Transform, f-k Transform）及

慢度頻率轉換法（Slowness- Frequency Transform, p-w Transform），其可用

以區隔表面波與實體波，亦可用來求取表面波之頻散曲線(McMechan and

Yedlin，1981；Gabriels et al.，1978；Park et al.，1998)。雖說有不同域的轉換，然而不同域的轉換，實際的獨立變數僅有兩個，因此這些方法在數學上是相關連，只是由於需求性的不同而使其表現在所需的物理量上(Lin and Chang, 2004)。以頻率波速轉換法為例，將多處接收器位置之二維震測資料 ) , ( u x t 透過在時間上以及空間上的傅利葉轉換後得到其在頻率-空間頻率域的函數U(k,ω)，又 k c=ω 且 π ω 2 = f ，對U(k,ω)進行變數變換即可得到波速與頻率的函數U(c, f)。若以空間頻率的波譜大小為色階，可得到此波場轉換後的結果如圖 2.15 所示。圖中白線便是頻散曲線位置所在，亦是色階最深色處所在。若仔細注意，會發現在低頻位置的色階分佈較不集中，這是因為在離散的傅利葉轉換中，會因為有限的空間位置而產生的遺漏(leakage)現象，越低頻其波長越長，影響越顯著。待由此獲得頻散曲線後，便可透過反算技術獲得以測線中點為代表值的地層剪力波速剖面。圖圖圖圖 2.2.2.11115555 多2. 多多頻道波場轉換法頻散曲多頻道波場轉換法頻散曲頻道波場轉換法頻散曲線分析結果頻道波場轉換法頻散曲線分析結果線分析結果線分析結果

Pseudo-section概念於表面波震測應用之數值模擬探討

國立交通大學

國立交通大學

國立交通大學

國立交通大學

土木工程學系碩士班

土木工程學系碩士班

土木工程學系碩士班

土木工程學系碩士班

碩士論文

碩士論文

碩士論文

碩士論文