國立交通大學
土木工程學系
碩士論文
應用三維波傳模擬表面波震測於複合土體有效剪力
波速量測之探討
The study of effective shear wave velocity measurement
in composite soil by application of three-dimensional
simulation of surface wave
研究生：張友誠
指導教授：林志平 博士
應用三維波傳模擬表面波震測於
複合土體有效剪力波速量測之探討
The study of effective shear wave velocity
measurement in composite soil by application of
three-dimensional simulation of surface wave
研 究 生：張友誠 Student：Yu-Cheng Chang
指 導 教 授：林志平 博士 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 2013
Hsinchu, Taiwan, Republic of China 中華民國一○二年七月
應用三維波傳模擬表面波震測於
複合土體有效剪力波速量測之探討
研究生：張友誠 指導教授：林志平 博士 國立交通大學土木工程學系中文摘要
為了提昇軟弱地盤之強度與勁度，大地工程師們常用高壓噴射灌漿或 是擠壓砂樁等地盤改良工法以達到上述目的，由於改良後會在地中形成固 結物，使得土壤為高度非均質性。現有一課題值得探討：如何決定改良成 效？傳統檢核方法如鑽孔取樣、標準貫入試驗因取樣空間小、耗時且成本 高等限制，想要取得此地盤整體平均的改良率有相當的難度，相反的表面 波震測具有大範圍取樣，非破壞性及快速施作等優點，其反算得到的剪力 波速剖面實質上反映了複合土體之彈性行為，但表面波方法有一基本假設 是測線下方為水平一維層狀地層，使得非均質土壤在剪力波速量測上卻有 均質化的特性，為此文獻進行了二維數值模擬探討，評估地盤改良率與剪 力波速提昇率之關係，但二維模擬建立的地層模型無法反映實際地盤改良 場址，加上表面波震測會因側向地層變化造成訊號收錄有一定程度之影響， 進行三維數值模擬有其必要性。 本研究採用三維波譜元素法進行探討，研究成果顯示，表面波於改良 後場址之行為，若是簡化牆狀地層模型，其剪力波速之均質化會遵守一維 等效波速理論之下限；若於柱狀地層模型，剪力波速之均質化不再遵守一 維等效波速理論之下限，而是在等效模數理論之上下限之間，可根據模擬 資料提出均質化等效速度與改良率之關係。表面波的側向影響範圍約為波 長的 1.3 倍，在同樣範圍內，測線外異質土壤與測線內原土壤的速度比越大， 量測速度受到的影響越大；當速度比相同的情況下，原土壤速度越低，量 測速度受到的影響也越大。 關鍵字：表面波、表面波震測、波譜元素法、均質化、側向變化地層The study of effective shear wave velocity measurement
in composite soil by application of three-dimensional
simulation of surface wave
Student：Yu-Cheng Chang Advisor：Dr. Chih-Ping Lin Department of Civil Engineering
National Chiao Tung University
Abstract
Soft ground is often improved to increase strength and stiffness by methods such as jet grouting and stone column which result in heterogeneous ground with improved columns. Conventional methods used to assess such ground improvement are subjected to several limitations such as small sampling volume, time-consuming, and cost ineffectiveness. It’s difficult to assess the average property of the improved ground and the actual replacement ratio of ground improvement. The use of seismic surface wave method (i.e. multi-station analysis of surface wave, MASW) for such a purpose seems to be a good candidate. But the surface wave method is essentially a 1-D method assuming horizontally-layered medium. What MASW measures in the highly heterogeneous improved ground remains to be investigated. The feasibility of MASW in characterizing the heterogeneous ground with improved columns has been demonstrated and the relationship between replacement ratio and velocity increase has been studied by 2D model. This study aims to investigate the homogenization of shear wave velocity measured by MASW in real 3D condition. Testing results may also be affected by the survey line location relative to the improved columns. The lateral sampling space of the surface wave testing was also investigated by 3D modeling.
3D spectral element method was utilized for the numerical investigation. Numerical results show that the effective shear velocity can be derived based on lower bound of equivalent velocity model when the heterogeneous ground is a 2D wall-like structure. However, in real 3D column-embedded structure, the relation between the effective shear wave velocity and replacement ratio falls between upper bound and lower bound of the equivalent modulus model. A predictive model is derived based on the numerical simulations. The lateral sampling range is about 1.3 wavelength. In the same spatial range, the measured velocity is affected more by offline heterogeneous material for higher velocity contrast. At constant velocity contrast, lower velocity material along the survey line is affected more by offline heterogeneity.
Key words: surface wave, surface wave method, spectral element method, homogenization, lateral heterogeneous layer
致 謝
研究期間承蒙林志平教授無論在邏輯思考或是表達能力等方面指出明 顯的缺點與改進的方向，讓學生對於從未接觸過的領域有了突飛猛進的成 長；再來感謝俊宏與柏林兩位學長對於學生研究上的指導與協助，總是不 厭其煩地悉心解釋遭遇到的問題，得以有靈感嘗試找尋解決的辦法；此外 仍要謝謝口試委員們不吝指教與建議，使得本論文之研究能夠趨近完善； 最後要感謝就學期間遇到的諸位大地組教授們，不論是學業或是待人處事 等方面都獲益匪淺。 同時還要感謝研究室每一位成員，無論是分擔工作或是生活上彼此間 的照料，讓我有十足的信心可以面對困難，加上大地組同學們總是樂觀進 取，感染正值困惑的我，學習何謂永不放棄，才有今日的成就。 最後再獻上個人由衷的祝福，希望所有曾經幫助過我的人能夠事事盡 心，身體健康，千言萬語也感激不盡！目錄
中文摘要 ... I Abstract ... II 致 謝 ... III 目錄 ... IV 圖目錄 ... VI 表目錄 ... IX 第一章 緒論... 1 1.1 研究動機 ... 1 1.2 研究目的 ... 4 1.3 研究內容 ... 5 第二章 文獻回顧 ... 6 2.1 表面波震測 ... 6 2.1.1 表面波基本波傳原理 ... 6 2.1.2 表面波量測方法 ... 9 2.2 數值方法 ... 16 2.1.2 有限差分法 ... 16 2.2.2 波譜元素法 ... 20 2.3 表面波震測於地盤改良成效之評估 ... 28 2.3.1 現場實測案例 ... 28 2.3.2 理論與數值模型 ... 31 2.4 綜合評析 ... 40 第三章 研究方法 ... 41 3.1 三維震波模擬方法與選用軟體 ... 42 3.1.1 多頻道頻率波速轉換分析方法 ... 423.1.2 三維震波模擬軟體 SPECFEM3D 架構說明 ... 43 3.1.3 地層模型設計說明... 46 3.1.4 地層剪力波速反算軟體簡介 ... 48 3.2 表面波震測資料模擬參數設定 ... 49 3.2.1 表面波震測施測參數 ... 49 3.2.2 數值模擬參數 ... 49 第四章 結果與討論 ... 51 4.1 地層均質化傾向之行為模式 ... 53 4.2 地層側向變化之影響 ... 56 4.3 地盤改良成效檢測之適用性 ... 63 第五章 結論與建議 ... 68 5.1 結論... 68 5.2 建議... 68 參考文獻 ... 69
圖目錄
圖 1.1 (a) 表面波震測測線配置 (b) 改良前後剪力波速之變化 (c) 對應之 剪力波速提昇率 ... 3 圖 1.2 表面波震測側向取樣空間可能範圍(如圖中虛線)示意圖 ... 4 圖 2.1 表面波示意圖 (http://www.geo.mtu.edu/UPSeis/waves.html) ... 7 圖 2.2 表面波震測技術流程圖 ... 10 圖 2.3 柏松比與波速度之關係(Richart et al., 1970) ... 10 圖 2.4 SASW 現地施測示意圖 ... 11 圖 2.5 SASW 頻散曲線分析(相位差對頻率圖) ... 12 圖 2.6 SASW 頻散曲線分析結果 ... 12 圖 2.7 地層剪力波速剖面 ... 13 圖 2.8 MASW 現地施測示意圖 ... 13 圖 2.9 多頻道波譜法頻散曲線分析(相位-空間位置圖)... 15 圖 2.10 多頻道波譜法頻散曲線分析結果 ... 15 圖 2.11 多頻道波場轉換法頻散曲線分析結果 ... 16 圖 2.12 有限差分法網格系統示意圖 ... 17 圖 2.13 速度錯置網格架構，時間定義在 k+1/2 ... 19 圖 2.14 應力錯置網格架構，時間定義在 k+1 ... 19 圖 2.15 有限地球模型示意圖，震源 xs可置於模型內任意位置(Komatitsch and Tromp, 1999) ... 21圖 2.16 網格分割示意圖(Komatitsch and Tromp, 1999) ... 22
圖 2.17 六面體元素之幾何形狀示意圖，空心代表可以省略(Schuberth, 2005) ... 23
圖 2.18 一階與二階形狀函數示意圖 (a) 兩個節點 (b) 三個節點 (Schuberth, 2005) ... 23
圖 2.19 Lagrange 多項式於[-1, 1]內插，nL=4，有 nL+1=5 個控制點(Komatitsch et al., 2005)... 25 圖 2.20 二維波場離散化示意圖(Komatitsch et al., 2005) ... 26 圖 2.21 高壓噴射灌漿場址配置圖... 29 圖 2.22 表面波震測法高壓噴射灌漿場址：(a) 改良前後剪力波速剖面；(b) 改良前後剪力波速提升率 ... 29 圖 2.23 擠壓砂樁改良場址配置圖... 30 圖 2.24 擠壓砂樁改良場址剪力波速剖面 ... 30 圖 2.25 一維等值理論模型示意圖... 32 圖 2.26 一維等值理論模型下之剪力波速提升比與剪力波速比 Vm / Vc之關係 圖(改良率為 15%，ρm / ρc = 0.85)... 34 圖 2.27 一維等值理論模型下之剪力波速提升比與改良率之關係圖(剪力波 速比為 4.7，ρm / ρc = 0.85) ... 35 圖 2.28 數值模型參數與配置示意圖 ... 36 圖 2.29 數值模擬地層模型配置示意圖 ... 38 圖 2.30 表面波震測數值模擬分析結果：(a) 改良率 30%頻散曲線影像圖； (b) 改良率 30%頻散曲線影像圖；(c) 反算分析之剪力波速剖面；(d) 剪 力波速提升率 ... 39 圖 2.31 二維模擬下(A1~A6)地盤改良率與剪力波速提昇率之關係圖 ... 40 圖 3.1 研究方法流程圖 ... 41 圖 3.2 多頻道波場轉換法分析流程：(a) 時間域資料(t-x domain)；(b) 頻率 域資料(f-x domain)；(c) 空間域離散化至波數域；(d) 頻散曲線 ... 43 圖 3.3 Cubit 圖形化使用介面 ... 44 圖 3.4 SPECFEM3D 模擬流程圖 ... 45 圖 3.5 探討均質化行為地層模型簡圖 ... 47 圖 3.6 側向取樣空間範圍地層模型簡圖 ... 47
圖 3.7 含改良樁地層模型簡圖：(a) 俯視圖；(b) 側視圖... 48 圖 3.8 SurfSeis 使用者介面圖 ... 49 圖 3.9 高斯函數模擬震源時間函數示意圖 ... 50 圖 4.1 水平層狀地層模型與測線位置示意圖 ... 51 圖 4.2 頻散曲線分析結果 ... 52 圖 4.3 數值模擬地層模型與測線位置示意圖 ... 53 圖 4.4 表面波震測數值模擬分析結果：(a) 改良率 40% 時間域震波模擬資 料；(b) 改良率 40% 頻散曲線影像圖；(c) 反算分析之剪力波速剖面； (d) 對應剪力波速提升率 ... 55 圖 4.5 剪力波速提升率與改良率之關係圖 ... 56 圖 4.6 夾層地層模型與測線位置示意圖 ... 57 圖 4.7 (a) 基態相位速度對頻率的關係圖(剪力波速比為 4.7)；(b) 基態相位 速度對頻率的關係圖(剪力波速比為 2) ... 59 圖 4.8 波長與相位速度關係圖 ... 60 圖 4.9 基態相位速度對頻率的關係圖(剪力波速比為 4.7) ... 62 圖 4.10 測線至改良土壤之距離與表面波波長關係圖 ... 62 圖 4.11 柱狀地層模型示意圖 ... 64 圖 4.12 表面波震測數值模擬分析結果：(a) 改良率 25% 時間域震波模擬資 料；(b) 改良率 25% 頻散曲線影像圖；(c) 反算分析之剪力波速剖面； (d) 對應剪力波速提升率 ... 66 圖 4.13 剪力波速提升率與改良率之關係圖 ... 67
表目錄
表 2.1 參數數值表 ... 37 表 2.2 地層模型配置數值表 ... 37 表 3.1 探討均質化行為地層模型參數表 ... 46 表 3.2 表面波震測施測參數表 ... 49 表 3.3 數值模擬參數表 ... 50 表 4.1 土壤參數設定表 ... 51 表 4.2 土壤參數設定表 ... 54 表 4.3 材料參數設定表(剪力波速比固定為 4.7)... 57 表 4.4 材料參數設定表(剪力波速比固定為 2) ... 58 表 4.5 土壤參數設定表 ... 61 表 4.6 地層材料參數設定表 ... 64第一章 緒論
1.1 研究動機
近年地球物理探測方法快速興起，其快速、省時且經濟之優點可應用 於工程地質調查之地層構造探測、地層動態性質探測、地層震動特性探測、 地下埋設物探測及地下水調查等各方面，地球物理探測方法之選擇，須依 調查之目的而定，同時必須考慮探測現場之環境、地形、地質、土層性質 以及土地利用狀況，配合地質及鑽探資料解釋比對。因各種地球物理探測 方法各有其限制及解析能力，除要有經驗的探勘人員施測與解析之外，有 時更須數種地球物理探測方法同時併用，相輔相成以獲得更多資料作正確 之判釋，作為設計、施工之參考。 台灣地狹人稠、土地資源有限，又因地質多變，因此面臨承載力不足 與壓縮量過大之軟弱地盤的機會漸增，地盤改良技術成為大地工程師處理 問題的利器之一。在地盤改良方法中以增加土壤剪力強度以及承載力為目 標之改良方法如高壓噴射攪拌工法、機械攪拌工法、擠壓砂樁等工法，施 作後將會在地中形成改良柱體，該柱體的存在造成改良完成之土體具有高 度的不均質性。以往地盤改良的檢核，不外乎鑽孔取樣、標準貫入試驗及 CPT 等試驗，然而地盤改良後之土體為高度不均質的複合土體，影響改良 成效的因素眾多，且相當複雜，採用此些試驗方法所得之結果由於取樣空 間小，不具整體代表性，且易受局部變異而誤判全域狀況，因此採用可大 面積進行快速且具整體代表性之檢驗方式具有其必要性。 近年來表面波震測法技術研究日趨成熟，透過對其所得之剪力波速剖 面之解析，已成功應用於地層比對、鋪面厚度檢測、土壤改良成效檢測、 地下孔洞探測、液化潛能分析以及隧道襯砌與背填灌漿之品質管制等地工 與環境問題解決上。表面波震測法係利用人為方式在地表產生震波，於地 表傳遞時，因在不同波長下之影響深度不同，可將地層不同深度內之地層物理性質反應於波速上。表面波震測的施做是先於試驗場址佈設受波器後， 製造人為震源產生擾動，將此地表垂直位移訊號收錄後，對其進行訊號分 析而獲得表面波頻散曲線，再透過反算技巧對此頻散曲線進行反算分析便 可得到現地一維地層剪力波速剖面。故表面波施作上不需鑽孔破壞地表土 層，具有容易、快速又經濟的優點，且其最大特點在於所得為測線下方土 層之綜合表現，其取樣體積遠大於傳統之試驗方法，因此更能代表改良土 體整體之改良效果，具有應用地盤改良成效檢驗之潛力。 利用表面波震測進行地盤改良的檢核主要是針對具強度提升需求之地 盤改良，檢核的方法是在地層進行地盤改良前後於同一位置進行表面波震 測，以獲得地層地盤剪力波速剖面。以戴永政(2009)於宜蘭所進行之地盤改 良檢核為例，L1~L4 為配置之表面波震測測線(見圖 1.1a)，於改良前後得到 之剪力波速變化如圖 1.1b 所示，剪力波波速於改良前後提升率約 14~17%(圖 1.1c)，與設計改良樁體積替換率 14.43%相仿 (後續簡稱為改良率，其定義 為改良樁總截面積除以改良區域總面積) 戴永政認為兩者間具相關性，而提 出剪力波速提升率即可表示整體改良率之假說，此一假說後續則有陳逸倫 (2011)透過二維數值模擬嘗試以表面波波長與改良樁間距的角度檢視有無 規則可循。林志平等人(2012)持續配置不同的模擬參數，進一步發現，隨著 地盤改良率的提升，剪力波速提升率也隨之增加，且模擬的結果會遵循一 維等效波速理論行為之下限，但二維模擬在地層模型配置上，簡化改良樁 為牆狀，無法實際反應地盤改良場址，因此本研究將以三維波傳模擬角度 出發，針對二維設置的地層模型加以三維化，檢核是否仍依循一維等效波 速理論行為之下限，再進一步建立含有改良樁柱體之地層模型，期望可以 更加了解表面波於高度異質性土體下反應的有效剪力波速量測之行為模 式。
圖 1.1 (a) 表面波震測測線配置 (b) 改良前後剪力波速之變化 (c) 對應之 剪力波速提昇率
(b) (c)
此外，上述之地盤改良場址施做表面波震測發現另一現象：在測線規 劃上，L1 及 L4 測線有通過改良樁，而 L2 及 L3 並未通過，然而比較位於 樁上的 L1、L4 測線結果與佈設於樁間的 L2、L3 測線，如圖 1.1a 所示改良 後結果卻十分相近。此結果顯示表面波震測法在設計樁徑 1.2m、樁距 2.8m 之高壓噴射灌漿改良場址無法區別出樁體與樁間的差異，表面波施測結果 之最小波長約為 5 m，造成此一結果之可能原因在於表面波震測法在垂直測 線方向之影響範圍同時涵蓋樁體區域與樁間區域，使得量測結果為灌漿後 土體之整體綜合表現，因此本研究將嘗試界定表面波在地層中傳遞於側向 可能涵蓋之範圍(見圖 1.2)，提供初步的探討。 圖 1.2 表面波震測側向取樣空間可能範圍(如圖中虛線)示意圖
1.2 研究目的
本研究主要利用三維數值模擬，將二維數值模擬建立的地層模型三維 化，以期驗證表面波震測於複合土體之均質化行為；再針對表面波震測測 線在側向變化上剪力波速對比之地層中探討表面波震測受到側向影響之範 圍；最後重建地盤改良場址，檢視地盤改良率與剪力波速提昇率之關係， 提供表面波震測做為地盤改良成效檢測之依據。1.3 研究內容
本研究共分為五個章節做為陳述，於第一章 緒論中就引起此次研究的 動機做說明(1.1 研究動機)，並將此研究所欲達到的成果列出(1.2 研究目的)。 在第二章中，將由介紹表面波震測開始，了解表面波於地層中傳遞時有何 特性，並介紹地球物理探測常見的表面波震測方法(2.1 表面波震測)，接著 列舉常見的數值模擬方法，針對三維震波模擬軟體採用的數值方法做詳加 介紹，最後經由現場案例與二維數值模擬探討回顧，用以檢視表面波震測 於高度異質性土壤，量測得到的剪力波速可能的均質化行為模式，對照現 場改良率有何關係。 經由第二章的文獻回顧後，統整收集到的文獻資料，於第三章中說明 所使用之三維震波模擬軟體，以及模擬震測方法的施測及其分析的詳細過 程(3.1 三維震波模擬方法與選用軟體)，而後說明本研究中模擬所用的參數 (3.2 表面波震測資料模擬參數設定)。在第四章中以不同的模擬地層類型分 別討論表面波均質化的行為，及側向變化地層對於表面波震測測線量測震 波訊號所造成的影響。最後在第五章中提出本研究的結論，並以研究所得 的結果提出建議。第二章 文獻回顧
2.1 表面波震測
2.1.1 表面波基本波傳原理
在均質 (homogeneous) 、等向 (isotropic) 之半無限域 (halfspace) 中， 實體波 (body wave) 傳遞至地表面，由於邊界條件的限制，產生沿地表面 傳遞之振波，稱之為表面波 (surface wave) 。Rayleigh 於 1885 年首先進行 研究，沿彈性半無限空間自由表面傳遞之波，其擾動 (disturbance) 僅侷限 於自由邊界之有限範圍內，該波的影響隨深度的增加而迅速減低。因此此 彈性表面波日後即以雷利波 (Rayleigh wave) 命名之。表面波包括雷利波和 拉夫波 (Love wave) 兩種。雷利波在傳遞時，地層介質顆粒同時有與傳遞 方向垂直及平行的運動，亦即當其向+x 方向傳遞時，在地表面處，介質顆 粒會以原位置為中心，形成橢圓形的路徑，以逆時鐘方向震動。其所形成 的橢圓型平面只會在 x、z 平面上，而 y 方向上並無分量 (見圖 2.1a)，與同 在 x、z 平面上行進的壓力波 (compression wave，或稱 P 波) 及 x、z 平面 上的垂直向剪力波 (SV 波) 之分量，統稱為 P-SV 波。拉夫波 (圖 2.1b) 具 有水平向剪力波 (SH 波) 之分量，主要是因層狀土層因實體波在土層中發 生大量的折射及反射，產生強度集於介質表面的波形，質點運動方向為沿 著水平面與波傳行進方向垂直。 然而，於工程上使用之表面波震測是應用雷利波之特性做為分析標的， 為了便於解釋，後續之表面波專指雷利波。表面波傳遞時，其波傳影響範 圍大約侷限於一個波長之深度內，因此，各個不同波長的表面波所反映之 地層深度特性將有所不同，亦即表面波影響深度隨頻率之不同而異。當土 層剪力模數隨著深度變化，不同頻率所造成的表面波波速便有所不同，此 一特性稱之為頻散現象 (dispersion)，集合不同頻率所對應之表面波波速，
可得一曲線，稱之為頻散曲線 (dispersion curve)。表面波震測之基本原理即 藉由量測表面波於不同頻率之相位速度，獲得其頻散曲線後再進行剪力波 速度層構造之反算分析，近年來應用於非破壞檢測方面也相當廣泛，如土 層剖面及道路鋪面的探測、RC 結構強度的檢測及海洋地質構造的測量等。 戴永政(2009) 更提出應用表面波震測於地盤改良成效之評估，後續則有陳 逸倫(2011)、林志平等(2012) 嘗試利用數值模擬進行表面波震測法做為地盤 改良成效檢測工具之探討，後續小節將做進一步的介紹。 圖 2.1 表面波示意圖 (http://www.geo.mtu.edu/UPSeis/waves.html) 假設地層為均質、均向之線彈性體，若於地表面上某一方向 x，量測 地表隨時間 t 之垂直運動 u(x,t)，則就某一角頻率 ω (= 2πf ) 而言，雷利波 運動之通解可表示為： (x,t)=U0(ω)A(x)e-jkxejwt=U 0(w)A(x)e jw[t- x w/k] (2.1) 其中，U0為初始振幅之大小，與震源型式相關；A 為振幅隨空間改變因子， 與幾何阻尼相關；k 為空間頻率之大小，又稱波數(Wavenumber)。k 之倒數 λ = 2π/k 即為波長。式(2.1)中顯示波傳之主要行為，包含波傳之衰減與波動 (a) (b)
之速度，其中波傳速度(v)與材料之彈性模數有直接相關，且可定義為時間 頻率與空間頻率之比值，如下所示： v=ω k = 2πf k (2.2) (2.2)亦可表達成頻率(f)與波長(λ)之乘積： v = f · λ (2.3) 在均質的線彈性體中，(2.3)式之乘積為一定值，頻率較低之波長較長， 而頻率較高之波長較短，亦即波傳之速度為一常數，不隨頻率之改變而改 變。在均質線彈性體中，僅有單一模態之波傳，然而，一般層狀地層沿深 度方向為非均質，亦即其材料之彈性模數隨深度而變化。若考慮地層之非 均質性及材料阻尼，就某一角頻率 ω (= 2πf ) 而言，雷利波之通解可由振態 疊加求得： (x,t)= Sm(w,x)e-j(wt-kmx) m (2.4) 其中，m 為振態數；S 為振幅因子，包含震源、受波器、幾何阻尼及材料之 綜合影響。 由於表面波之波傳影響範圍大約侷限於一個波長之深度內，而不同頻 率震波具有不同之波長，當地層之速度隨深度而變化時，各振態之相位速 度變為頻率之函數，此速度隨頻率之變化稱為頻散現象。(2.4)可改寫成： (x,t)=U0(ω)A(α,x,ω)e-j(wt-ψ(x,ω)) (2.5) 其中 A 亦受到材料阻尼(α)之影響，且為頻率之函數，ψ為一複合相位函數， 描述相位隨著空間變化所產生之波傳現象。任一頻率簡諧波之特徵點 (例如 波峰或波谷) 可以下式表示： ωt-ψ(x,ω)=const (2.6) 因此，將(2.6)式對時間微分可以得到有效相位速度之定義：
v= ω ∂ψ(x,ω) ∂x (2.7) 由(2.7)式可得知，多重模態波傳之有效相位速度為區域性之物理量，亦即 不同位置具有不同之波傳速度。若地層為常態之速度剖面 (亦即，剪力模數 隨深度增加而增加)，則通常波傳由基態所控制，而有效相位速度趨近於基 態之速度 (林志平等，2002)。 2.1.2 表面波量測方法 表面波震測法主要可以區分為三大步驟：現場施作、頻散曲線分析以及 地層反算(如圖 2.2 所示)，先於欲試驗之場址安裝地音計(geophone)後，製 造人造震源產生擾動，將此地表垂直位移訊號收錄後，對其進行訊號分析 而獲得表面波頻散曲線，由於雷利波速度與剪力波速度相近(參考圖 2.3 Richart et al., 1970)，利用頻散曲線與地層波速變化之關係，再透過反算技 巧對此頻散曲線進行反算分析便可得到現地一維地層剪力波速剖面。由於 在現地施測時所採用接收器數目的不同，可將表面波震測法分為波譜分析 法 (Spectral Analysis of Surface Wave, 簡用 SASW) 以及多頻道表面波震測 法 (Multi-channel Analysis of Surface Wave Method, 簡用 MASW)，又由於 MASW 提供了較多空間上的訊息，使得其頻散曲線的分析上有不同的方法， 根據採用的方法不同，又可將 MASW 分為多頻道表面波譜法 (Multi-channel Spectral Analysis of Surface Wave Method, 簡用 MSASW) 以及多頻道波場 轉換法 (Multi-channel Wave Transform of Surface Wave Method, 簡用 MWTSW)。
圖 2.2 表面波震測技術流程圖 圖 2.3 柏松比與波速度之關係(Richart et al., 1970) 2.1.2.1 SASW 1. 現地施測 SASW 在現地的施測上(如圖 2.4 所示)僅採用兩個接收器，以一近站支 距(即震源至第一個接收器位置的距離，一般而言，近站支距 x0採用與兩個
接收器相隔間距 Δx 相同的距離，即 x0 =Δx。) 左敲、右敲，紀錄後改變接 收器間距，使得兩個接收器的中點位置固定在同一個地方，而後左敲、右 敲、紀錄，再繼續改變兩個接收器的間距，重複紀錄的動作。 圖 2.4 SASW 現地施測示意圖 2. 頻散曲線分析 透過現地施測，每一筆紀錄我們都將可收集到兩個接收器位置的地層 震動訊號，此訊號可表示為位置與時間的函數 u1(x1,t), u2(x2,t)。透過傅利葉 轉換(Fourier Transform)可得到此二函數在頻率域的函數 U1(x1,ω), U2(x2,ω) 以及各頻率在 x1, x2處的相位角。根據相位速度的計算式： c(ω)= ω Δψ(ω) Δx (2.8) 其中，c(ω)表示不同頻率時的雷利波波速；Δψ(ω)表示不同頻率時在 x1, x2 處相位角的相減；∆x 表示 x2-x1。將不同頻率所計算得到的雷利波相位速度 畫上，便可得到頻散曲線。在計算過程中，相位差∆ψ 可直接由 u1(x1,t), u2(x2,t)
互能頻譜(Cross Spectral Density)的相位角獲得，即
Δψ(ω)=Angel(CSD(u1(x1,t), u2(x2,t))) (2.9) 如圖 2.5(虛線)所示，由於頻譜分析所得之相位角差僅侷限在-π~π 之間，因
此在計算(2.8)式前須先將其摺開(unwrap)，得到如圖 2.5 (實線)的結果。由 圖 2.5 中亦可看見，在所得的頻譜分析結果中會有一段資料不佳段，不同的 接收器間距會有不同位置及寬度的資料不佳段，因此在施測時收錄多筆不 同間距的資料，分別計算完得到頻散曲線後取其算術平均數以為最後頻散 曲線結果(如圖 2.6 所示)。最後藉由反算技術，獲得以接收器終點為代表值 的地層剪力波速剖面(如圖 2.7 所示)。 圖 2.5 SASW 頻散曲線分析(相位差對頻率圖) 圖 2.6 SASW 頻散曲線分析結果
圖 2.7 地層剪力波速剖面
2.1.2.2 MASW
1. 現地施測
多頻道表面波震測法(Multi-channel Analysis of Surface Wave Method， 簡稱 MASW 法)，首先是由地球物理領域之學者提出(McMechan and Yedlin, 1981; Gabriels et al., 1987；Park et al., 1999)，其在現地施測時，採用多個在 同一直線上的接受器(如圖 2.8 所示)，只需敲擊一次，便可完成。在實際操 作上，一般採用 1~2 公尺之受波器間距，並設置 12 個以上之接收器接收震 源所發出之震波訊號。以 24 個受波器為例，在第一個受波器之線外取適當 近站支距，反覆在同一震源處施作，將其疊加以消除雜訊之影響，直至收 錄到清晰之表面波訊號為止。 圖 2.8 MASW 現地施測示意圖
2. 頻散曲線分析 多個接收器的收錄，提供了更多空間上的訊息，使得多頻道表面波震 測法在頻散曲線的分析上可利用二維訊號處理技術分析震測資料，求得訊 號品質優良之頻散關係曲線。根據不同的訊號處理技術，多頻道式表面波 量測法可再分為多頻道波譜分析法 (MSASW) 以及多頻道波場轉換法 (MWTSW) 兩種。 (1) 多頻道波譜分析法 多頻道波譜分析法與波譜分析法的頻散曲線分析類似，經由將多處接 收器位置之震測資料 ui(xi, t)透過傅利葉轉換後得到其在頻率域的函數 Ui(xi,ω)，以及各頻率在 xi處的相位角。如此一來，對於每一個特定頻率， 皆能獲得其在不同接收器位置的相位角(如圖 2.9 所示)，將頻譜分析所得之 相位角摺開後，由(2.8)式可知，由此斜率便可求得此特定頻率的相位波速。 在圖中亦可看見，線段中會有資料不佳的段落，依據不同的頻率資料不佳 的段落出現位置會不同，一般來說，低頻時會發生於距離震源較近的位置， 稱作近場效應(near field effect)，而高頻時會發生於距離震源較遠的位置， 稱作遠場效應(far field effect)。分析過程中應當要注意此兩個效應所可能造 成頻散曲線分析摺開時的錯誤(此種效應的影響在 SASW 的分析中不易發現， 使得 SASW 的分析有潛在的錯誤)，於計算相位速度相關之斜率∆ψ/∆x 時， 可避開資料不佳段。對不同頻率進行相同的動作，便可得到此地層表面波 之頻散曲線(如圖 2.10 所示)。最後藉由反算技術便可獲得以測線中點為代 表值的地層剪力波速剖面。
圖 2.9 多頻道波譜法頻散曲線分析(相位-空間位置圖)
圖 2.10 多頻道波譜法頻散曲線分析結果 (2) 多頻道波場轉換法
多頻道波場轉換法常見於地球物理領域之濾波處理，包括頻率波數轉 換法(Frequency-Wavenumber Transform, f-k Transform)及慢度頻率轉換法 (Slowness-Frequency Transform, p-w Transform)，其可用以區隔表面波與實 體波，亦可用來求取表面波之頻散曲線(McMechan and Yedlin, 1981； Gabriels et al., 1987；Park et al., 1999)。雖說有不同域的轉換，然而不同域 的轉換，實際的獨立變數僅有兩個，因此這些方法在數學上是相關連，只
是由於需求性的不同而使其表現在所需的物理量上(Lin and Chang, 2004)。 以頻率波速轉換法為例，將多處接收器位置之二維震測資料 u(x,t)透過在時 間上以及空間上的傅利葉轉換後得到其在頻率-空間頻率域的函數 U(k,ω)， 又c=ω k且f= ω 2π，對 U(k,ω)進行變數變換即可得到波速與頻率的函數 U(c,f)。 若以空間頻率的波譜大小為色階，可得到此波場轉換後的結果如圖 2.11 所 示。圖中白線便是頻散曲線位置所在，亦是色階最深色處所在。若仔細注 意，會發現在低頻位置的色階分佈較不集中，這是因為在離散的傅利葉轉 換中，會因為有限的空間位置而產生的遺漏(leakage)現象，越低頻其波長越 長，影響越顯著。待由此獲得頻散曲線後，便可透過反算技術獲得以測線 中點為代表值的地層剪力波速剖面。 圖 2.11 多頻道波場轉換法頻散曲線分析結果
2.2 數值方法
2.1.2 有限差分法 有限差分法簡單的說便是在離散的向度上，將控制方程式及邊界狀態 中的微分項以差分式替代，使得原先的微分方程式變成簡單的四則運算式。要將有限差分法應用至特定的控制方程式或微分方程式中時，總括來看有 三個主要步驟：先(1)建構離散的有限差分離散模型，亦即建構離散的網格， 並對微分方程式導入有限差分式。然後對導出的有限差分運算式進行(2)有 限差分離散模型的分析，確定離散模型的收斂性、穩定性及餘尾誤差級數。 最後使用計算系統協助(3)疊代計算(iteration)。 1. 網格系統(Grid) 使用有限差分法的第一步便是要先決定所採用的離散網格系統。所謂 離散網格系統的意思是說，在微分方程式中，其所定義的時間及空間皆是 連續，而為了將有限差分式導入，必須將連續的空間以及時間向度離散， 而以離散的點代替原本連續的物理量。這樣離散後的時間及空間向度，便 稱做為離散網格系統，或簡稱網格。 在傳統的有限差分法使用上，皆是採用均勻方型網格，亦即其將空間 離散為均勻分布的格點(如圖 2.12a 所示)。而後為因應不同的使用需求而有 均勻三角形網格(Uniform Trigular Grid)、非均勻網格(non-uniform grid)等網 格系統的使用(如圖 2.12b,c,d 所示)。
2. 標準有限差分式(Standard Finite Difference Approximation)： 當離散網格系統選定後，便可開始選用適當的有限差分式導入微分方 程式中。習慣上，對於一開始發展有限差方法時應用在均勻方型網格系統 的有限差分式，便稱其為標準有限差分式。假設有一函數 F(x)存在有連續 的一次微分式，則此微分式可以採用有限差分式(2.10)~(2.12)式進行導入： 前向差分式 '( ) ( ) ( ) O(h) h x F h x F x F (2.10) 後向差分式 '( ) ( ) ( - ) O(h) h h x F x F x F (2.11) 中央差分式 F'(x) F(xh)2hF(xh)O(h2) (2.12) 其中： h 為極小的數 O 為餘尾誤差 O(h)稱作一級餘尾誤差 O(h2)稱作二級餘尾誤差，級數越大表示誤差較小。 餘尾誤差用來表示有限差分式的誤差估計，誤差級數估計方式是使用泰勒 展開式： ) ( ! ... ) ( '' ! 2 ) ( ' ) ( ) ( ( ) 2 x F n h x F h x hF x F h x F n n (2.13) ) ( ! ... ) ( '' ! 2 ) ( ' ) ( ) ( ( ) 2 x F n h x F h x hF x F h x F n n 對 F(x + h)、F(x - h)等位於有限差分式右側的函數進行展開後，以單一式子 或以線性方式組成上述三式結果，而剩下未出現在差分式中的展開項便以 O 來表示。若其展開至一次微分項則寫做 O(h)；若其展開至二次微分項則寫 做 O(h2 )，以此類推。 3. 有限差分法震波模擬 在標準有限差分法的架構下所進行的震波模擬，容易因為高波松比
(Poison’s Ratio)地層材料的存在而發生即使已滿足穩定性條件，亦會有不穩 定以及嚴重的假頻散現象產生的問題。為解決這問題，目前在有限差分法 網格系統中受到大量採用的錯置網格(Staggered Grid)系統是由 Madariaga 在 1976 年提出，應用在地震時地層裂隙的動態模擬上。錯置網格與均勻方型 網格的主要差異是在於，在錯置網格中，各物理量的定義位置不在同一格 點上，且相距離散間距的一半(如圖 2.13、圖 2.14 所示)。
圖 2.13 速度錯置網格架構，時間定義在 k+1/2
2.2.2 波譜元素法
一般除了有限元素法為常見的數值模擬方法之外，Patera 於 1984 年結 合有限元素法 (Finite element method, FEM) 和擬頻譜法 (Pseudo-spectral method)，提出波譜元素法 (Spectral element method, SEM)，此技術基於運 動方程式之弱式(Weak formulation)，並結合有限元素法的彈性與擬頻譜法 的精度。網格主要以六面體元素組成，元素中波場的離散化透過高階拉格 朗日內插 (High-degree Lagrange interpolants)，積分的計算則經由 GLL 積分 法(Gauss-Lobatto-Legendre quadrature)，結合上述方法可以形成對角質量矩 陣，大大的簡化了運算。 1. 波譜元素法之數學式 在有限的地球模型中(體積 Ω，如圖 2.15)，吾人試圖尋求解決震波形成 的位移場，此模型體積之邊界包含自由應力表面∂Ω 與人造吸收邊界 Γ，與 所有邊界向外垂直的單位向量標記作n，座標採用位置向量 x=(x, y, z)，為了 方便表示採用指數型符號 xi, i=1, 2, 3, x1=x, x2=y,x3=z。+xi方向之單位向量標 記作xi，相對的偏微分標記作∂i。 運動方程式主導震波在彈性介質中的傳遞，方程式可透過強式(strong formulation)或弱式求解：強式直接處理微分形式之運動方程式並結合邊界 條件，有限差分法即為使用強式求解之例子；弱式則將運動方程式改寫成 積分項，如有限元素法以及本文提到的波譜元素法。 波傳方程式之微分項可以用下式來表示： ρ∂t2s=∇∙T+f (2.14) 其中 ρ 為密度；s 是位移；T 為應力張量，於彈性非等向性之介質，根據虎 克定律應力張量與位移梯度∇s 呈線性關係，可用下列式子表示： T=c:∇s (2.15) 震源則可以表示成點作用力 f，可用彎矩張量 M 寫成：
f=− M∙∇δ(x − xs)S(t) (2.16)
xs為點作用力的位置，δ(x − xs)為位於 xs之 Dirac delta distribution，S(t)為震
源時間函數(source time function)。
運動方程式(2.14)必須求解於自由應力邊界條件之表面∂Ω： T∙n=0 (2.17) n是垂直於自由表面之單位向量。 圖 2.15 有限地球模型示意圖，震源 xs可置於模型內任意位置(Komatitsch and Tromp, 1999) 然而震波能量傳遞至人造邊界 Γ 需作吸收處理，為此才可模擬一個半 無限域空間，所以經常採用平行軸方程式來降低波場邊緣震波的能量 (Clayton and Engquist, 1977)：
T∙n=ρ[vn(n∙∂ts)n+v1(t̂1∙∂ts)t̂1+v2(t̂2∙∂ts)t̂2] (2.18) t̂1, t̂2為相互垂直的單位向量，且與人造邊界 Γ 相切；vn為n方向傳遞之準壓 縮波波速；v1, v2分別為t̂1, t̂2方向偏振之準剪力波速度。此邊界條件能完美 地吸收單一方向垂直邊界之衝擊，但是對於擦過邊界傳遞的波卻無法有效 降低影響。 於 SEM 只考慮積分項表示的運動方程式，(2.14)須與一任意的向量 w 作內積，於 FEM 稱為測試函數，再對整個模型體積 Ω 作分部積分，最後加 入上述提到的邊界條件，得到：
ρw∙∂t2sd3x Ω =− ∇w:Td3x Ω +M:∇w(xs)S(t) + ρ[vn(n∙∂ts)n+v1(t̂1∙∂ts)t̂1+v2(t̂2∙∂ts)t̂2]∙ Γ wd2x (2.19)
可以發現 f 項已利用 Dirac delta distribution 之特性直接積分；而使用弱式還 有另一好處是分部積分式中於自由應力表面(2.17)之積分可以抵銷。 2. 網格設計 SEM 所要模擬的空間域 Ω 將被分為數量 ne相互不重疊之元素 Ωe， e=1,…,ne使得Ω = ⋃ Ωe ne e=1 (如圖 2.16)，同樣的吸收邊界 Γ 也被分為 nb數量 之元素 Γb，b=1,…,nb使得Γ = ⋃ Γb nb b=1 。
圖 2.16 網格分割示意圖(Komatitsch and Tromp, 1999)
一般來說，FEM 可以使用四面體或六面體作為體積元素，但 SEM 限制 只能使用六面體，同樣的定義吸收邊界其表面元素形狀為四邊形，接下來 的部分將會說明如何將這些元素參數化。 3. 映射函數(mapping function)：轉換物理座標至標準區間 經由座標轉換，每一六面體元素 Ωe之物理座標將轉換成參考座標 ξ = (ξ, η, ζ)，映射至標準區間 Λ = [-1,1]×[-1,1]×[-1,1]，定義為： x(ξ) = Na na a=1 (ξ) xa (2.20) na為每一元素中的節點數，三維元素可以用 8 個節點表示六面體，也可以 用 27 個節點表示有曲面的元素(如圖 2.17)
圖 2.17 六面體元素之幾何形狀示意圖，空心代表可以省略(Schuberth, 2005) Na為形狀函數，由一階或二階 Lagrange 多項式之乘積；若多項式的階次為 nL，將有 nL+1 個控制點，-1≤ξi≤1, i=1,…,nL定義如下： linL(ξ)= ξ− ξj ξi− ξj nL j=0 j≠i (2.21) 若是 nL=1 即可產生直線之形狀函數，nL=2 即可產生曲線的形狀函數(如圖 2.18) (a) (b) 圖 2.18 一階與二階形狀函數示意圖 (a) 兩個節點 (b) 三個節點 (Schuberth, 2005)
所以，一個節點數為 8 的三維元素，其形狀函數可以用 nL=1 之 Lagrange 多 項式表示成： N1(ξ, η, ζ)=l01(ξ)l01(η)l01(ζ) N2(ξ, η, ζ)=l11(ξ)l01(η)l01(ζ) ⋮ N8(ξ, η, ζ)=l11(ξ)l11(η)l11(ζ) (2.22) 然而，元素 Ωe採用式(2.20)作座標轉換，數學上必須考慮座標轉換對於積分 造成的影響，積分式中將會新增額外一項，即 Jacobian(Je) dx dy dz = Je dξ dη dζ (2.23) 為了計算 Jacobian 首先須對式(2.20)偏微分，得到 Jacobi-matrix(J)，再取行 列式值： J = ∂x(ξ) ∂ξ = ∂Na(ξ) ∂ξ xa na a=1 (2.24) Je = detJ = ∂(x, y, z) ∂(ξ, η, ζ) (2.25) Jacobi-matrix 於網格分割中扮演著相當重要的腳色，它的行列式值不能等於 零，才能確保座標轉換是單一且不可逆的，並且影響數值模擬的穩定。 同樣的，對於二維邊界元素之座標轉換如下： x(ξ, η) = Na na a=1 (ξ, η) xa (2.26) 垂直邊界 Γb之向量n： n = 1 Jb ∂x ∂ξ× ∂x ∂η (2.27) Jb表示 Jacobian Jb = ∂x ∂ξ× ∂x ∂η (2.28) dx dy = Jb dξ dη (2.29)
4. 表示元素格點之內插函數 (interpolation function) 在模型網格化與使用弱式表示運動方程式後，下一步便是定義基函數 來表示元素內未知的位移向量。因此解析式(2.19)，現在必須藉由離散函數， 使用適當的插值方法估計離散點的值。傳統的 FEM 使用低階多項式表示元 素之幾何與向量，而 SEM 也用低階多項式表示元素之幾何，但所有未知的 函數則採用高階 Lagrange 多項式來內插，典型的 SEM 常使用 4~10 階 Lagrange 多項式作內插，這是 SEM 與 FEM 兩者主要的差別。
再仔細觀察式(2.21)，多項式li的值在座標 ξj剛好是 1 或 0(圖 2.19)，而 在這些離散點以外多項式會有任何可能的值，因為只對這些離散點的值有 興趣，所以並不違反下列的數學推導式。 可以用 Kronecker-delta 來表示這個特性： linL(ξ j)=δij (2.30) 而定義於式(2.21)控制點 ξi, i=0…nL的選擇則是透過 Gauss-Lobatto-Legendre points，也是下列式子之根： (1 − ξ2)Pn'L(ξ)=0 (2.31)
Pn'L就是 nL階 Lagrange 多項式的一次微分(Canuto et al., 1988)。這樣的選擇
是基於 Lagrange 內插與 GLL 求積的結合可使質量矩陣成為對角矩陣，因此 一個二維元素將有(nL+1) 2 個格點(圖 2.20)，三維元素將有(nL+1) 3 個格點。 圖 2.19 Lagrange 多項式於[-1, 1]內插，nL=4，有 nL+1=5 個控制點(Komatitsch et al., 2005)
圖 2.20 二維波場離散化示意圖(Komatitsch et al., 2005) 對於函數 g 在體積元素Ωe的 nL階拉格朗日內插表示為： g(x(ξ, η, ζ))≈ g(x(ξi, ηj, ζk)) nL i,j,k=0 Lijk(ξ)= gijkli(ξ nL i,j,k=0 )lj(η)lk(ζ) (2.32) gijk=g(x(ξi, ηj, ζk)) 表示 GLL 點為 x(ξi, ηj, ζk) 時函數 g 之值；為了簡明註記 多項式階次 nL予以省略；L 即為三個座標方向上一維 Lagrange 多項式之乘 積，定義於空間中的張量基底li×lj×lk作為函數的內插。這個過程也被稱作 張量化，意思是每個座標方向上的內插是相互獨立的。 利用(2.32)，函數 g 之梯度可表示為： ∇g(x(ξ, η, ζ))≈ xl∂l 3 l=1 g(x(ξ, η, ζ)) = xl gijk 3 l=1 3 l=1 [li(ξ)lj(η)lk(ζ)∂lξ+li(ξ)lj(η)lk(ζ)∂lη+li(ξ)lj(η)lk(ζ)∂lζ] (2.33) 單引號表示微分，xl, l = 1, 2, 3 表示沿座標軸 x, y, z 方向增加之單位向量。 估計 GLL 點 x(ξi, ηj, ζk) 於函數 g 之梯度可簡化為：
∇g(x(ξ, η, ζ))≈ xl [ gijk i=0 3 l=1 li(ξi)∂lξ+ gijk j=0 lj(ηj)∂lη+ gijk k=0 lk(ζ)∂lζ] (2.34) 上式透過積分的連鎖律(chain rule)，∂ξ ∂x⁄ 必須透過(2.24)之反矩陣求得，因 此適當的網格化讓反矩陣存在且唯一是相當重要的。 5. 元素體積之積分函數 為了求解(2.19)之運動方程式，每一元素之積分需透過數值方法計算。 傳統的 FEM 使用高斯求積分(Gauss quadrature)，SEM 則是採用 GLL 積分
法，結合 GLL 點可形成對角質量矩陣。體積元素 Ωe之積分可估計為： g(x)dx= g(x(ξ, η, ζ))Je(ξ, η, ζ)dξ dη dζ 1 -1 1 -1 1 -1 ≈ ωiωjωk nL i,j,k=0 g(ξi, ηj, ζk)Je(ξi, ηj, ζk) (2.35)
ωi > 0, i=0…nL為 GLL 積分法之權重(Canuto et al., 1988)
6. 組結全域系統 (assemble global system)
個別處理每一元素之後，可以從圖 2.19 發現 GLL 積分所採用的控制點 必會包含-1 和 1，也就是說這些鄰近的元素不論是邊界或是角落都會共用這 些格點，因此在整體結構下必須能夠辨別這些共用點，因此首要步驟便是 定義全域(global)與地域(local)的映射關係，可透過 FEM 或分類演算法去搜 尋之，並指派獨一無二的全域座標。一旦映射定義完成，分別在元素間計 算的內力也需要於共用點上加總起來，也就是 FEM 所稱的系統組合，全域 系統下隨著時間演變的常微分方程式一般項如下： MÜ+CU̇+KU=F (2.36) 其中： U：位移向量
M：質量矩陣 C：吸收邊界矩陣 K：勁度矩陣 F：震源項
2.3 表面波震測於地盤改良成效之評估
2.3.1 現場實測案例 戴永政於 2009 年之研究指出，表面波震測應用於複合改良土體，期望 透過表面波震測法快速、經濟、取樣範圍大的優點提供剪力波速變化以利 複合改良土體成效檢核上的參考。並以高壓噴射灌漿改良場址(圖 2.21)與擠 壓砂樁改良場址(圖 2.23)為例，初步探討表面波震測的側向取樣空間於其進 行複合土體改良成效評估的影響，以及改良後波速變化所代表的意義。提 出的論點如下： 1. 表面波震測所得之剪力波速剖面為其影響範圍內波傳通過地層之波速 平均，由兩個場址不同改良樁徑、改良間距、改良率(IR%)做為測試， 其結果如圖 2.22 顯示表面波震測在高壓噴射灌漿改良場址改良樁徑 1.2m 樁距 2.8m 的情況下是無法區分改良樁與樁間的差異，因此可視為 整體改良區剪力波速的綜合表現，然而在改良樁徑 0.7m 樁距 6m 之擠 壓砂樁場址施做結果如圖 2.24 顯示表面波震測通過改良樁體與樁間波 速有所差異，表示表面波震測法於側向上具有一定程度之取樣空間，於 測線佈設時，應考量側向取樣空間之影響。 2. 由高壓噴射灌漿改良場址得到試驗結果，表面波震測得到之剪力波速改 良率在 15~17%，與設計改良率 14.43%相近，初步認為在高壓噴射灌漿 改良工法中，採用表面波震測所得之剪力波速提升率與改良率(IR%)有 關係，表面波震測反應改良面積與整體改良面積的比例關係。圖 2.21 高壓噴射灌漿場址配置圖 圖 2.22 表面波震測法高壓噴射灌漿場址：(a) 改良前後剪力波速剖面；(b) 改良前後剪力波速提升率 L1_before L1_after L2_before L2_after L3_before L3_after L4_before L4_after 150 175 200 225 250 0 2 4 6 8 10 12 14 16 18 20 剪力波速, m/s 深 度 , m L1 L2 L3 L4 0 2 4 6 8 10 12 14 16 18 20 深 度 , m 剪力波速提升率, % 0 10 20 30
(a)
(b)
圖 2.23 擠壓砂樁改良場址配置圖
圖 2.24 擠壓砂樁改良場址剪力波速剖面
基於現場施測結果，建議後續研究應於更多不同改良率之複合改良土 體場址進行試驗，以印證表面波震測所得剪力波速提升率與改良場址改良 率之關係。此外，透過理論、數值模擬以及現場案例的探討，進一步了解
表面波震測的有效側向取樣空間；以提供未來在地盤改良區域中佈設較佳 測線位置(避免無法反應整體改良效果)之用。 2.3.2 理論與數值模型 從上節可知，表面波震測在改良前後剪力波速剖面差異所反應之工程 性質並無法明確的量化，於高壓噴射灌漿改良場址案例中，觀察到由表面 波震測所得之改良後剪力波速提升率與設計改良率具高度相似，進而猜測 由表面波震測所得之剪力波速提升率可能是反應現場之實際改良率。 針對具高度側向變化之高壓噴射灌漿之地盤改良工法，分別由一維等 值理論模型以及數值模擬方法兩方面進行剪力波速提升率與現地改良率之 關係探討。 2.3.2.1 一維等值理論模型 表面波於地盤改良場址中傳遞時，其波同時經過改良段以及未改良段， 然而由於表面波震測法在進行反算分析時是採用一維水平地層之反算方法， 因此所得之剪力波速剖面對於具有高度側向變化之高壓噴射灌漿場址會有 一均質化的過程，此均質化過程較為複雜受其波長不同而有所不同，在此 利用一維等值理論模型（波速等值模型與模數等值模型）進行簡化，先行 觀察改良率與剪力波速提升率之關係，並做為後面數值模擬結果之比較基 準。 如圖 2.25 所示，假設改良樁彼此平行，且土體和水泥漿之間膠合良好， 則根據波速等值理論，其等值波速上下限可分別將之視為改良樁與土體之 並聯與串聯傳遞，因此等值波速 VequV之上下限可表示如下： 1 VequV = fc Vc + fm Vm (下限) (2.37) VequV = fcVc + fmVm (上限) (2.38) 其中，fc為改良樁佔複合物之寬度比（即改良率），fm為土體佔複合物之寬
度比，Vc為改良樁之剪力波速，Vm為土體之剪力波速。 根據模數等值理論，其等值波速 VequG可表示為式(2.39) VequG= Gequ ρequ (2.39) 其中，平均密度可表示如式(2.40) equ = fc c + fm m (2.40) 依據混和定則，Voigt(1889)將異質性複合材料視為一組並聯彈簧，此「定 應變狀態」可視為多相異質性複合材料模數之上限值，如式(2.41)；Reuss (1929) 將異質性複合材料視為一組串聯彈簧，亦即為定應力狀態，可得多 相複合材料之模數之下限值，如式(2.42)。 Gequ = fcGc + fmGm (2.41) 1 Gequ = fc Gc + fm Gm (2.42) 上述 Gc及 Gm分別是改良樁與土體的剪力模數，且根據彈性力學理論 Gc=ρcVc2，Gm=ρmVm2， c及 m分別是改良樁與土體的密度。 圖 2.25 一維等值理論模型示意圖 為可更清楚的表示出改良前後剪力波速提升率與改良率的關係，定義 剪力波速提升率 VIR 為
VIR=Vequ-Vm
Vm (2.43)
其中，Vequ可為波速等值速度或模數等值速度。
依循剪力波速提升率之關係，根據(2.37)式，波速等值速度之剪力波速提升
率 VIRV之下限(lower bound)可表示為：
VIRV, low = 1-Vm Vc 1 fc+ Vm Vc -1 (2.44)
而根據(2.38)式，波速等值速度之剪力波速提升率 VIRV之上限(upper bound)
可表示為： VIRV, up = fcVc Vm + fm - 1 (2.45) 此外，依據(2.39)、(2.40)及(2.42)式，模數等值速度之剪力波速提升率VIRG之 下限(lower bound)可表示為： VIRG, low= 1 fc ρc ρm-1 +1 [ fc ρm ρc Vm Vc 2 -1 +1] - 1 (2.46) 依據(2.39)~(2.41)式，模數等值速度之剪力波速提升率VIRG之上限(upper bound)可表示為： VIRG, up= ρc Vc Vm 2 +ρmfm fc ρc+ρmfm fc - 1 (2.47) 由(2.44)(2.46)式中便可清楚得知，剪力波速提升率主要受到改良率fc以 及改良樁與土體速度比 Vm / Vc所影響，雖尚無法確認表面波震測所得之剪 力波速為何種等值模式，但由此一維簡單模型可初步瞭解剪力波速提升率
與改良率之間具有關聯。為進一步瞭解此些影響參數與剪力波速提升率間 之行為，設定改良率為 15%，ρm / ρc = 0.85，則剪力波速提升率與改良率的 關係如圖 2.26 所示。 圖 2.26 一維等值理論模型下之剪力波速提升比與剪力波速比 Vm / Vc之關係 圖(改良率為 15%，ρm / ρc = 0.85) 由圖中可見，即使改良樁之剪力速度持續增大，其對複合土壤的剪力波 速提升效果是會有其限度，以模數等值速度下限來說，當改良樁與土體剪 力波速比達 4.7 倍時，剪力波速提升率便已達極限約 8.3%，而在波速等值 速度下限中其極限值較大，且在剪力波速比達 2 時便可使複合土體之剪力 波速提升至 8.3%，由此看前述之噴射灌漿案例，不論其剪力波速比為何， 其剪力波速提升率至少在 10%以上，已超越模數等值速度下限之極限 8.3%， 所以可初步推論該處之表面波震測結果較接近波速等值速度下限之行為。 然而，由此一角度觀察並不易瞭解剪力波速提升率與改良率之關係，因此， 再從固定剪力波速比的方向觀察式(2.44)(2.46)，一般而言，高壓噴射改良樁 所採用之水泥漿其剪力波速可達原土體之 4.7 倍以上，且當其剪力波速比達 4.7 倍時，波速等值速度下限已達其極限，可簡化問題，因此設定剪力波速
比為 4.7 進行觀察，如圖 2.27 所示，波速等值速度下限之結果在相同的剪 力波速比下，較低的改良率便可有較高的剪力波速提升率，且改良率在 30% 以下時，兩者間關係接近是 1:1。 圖 2.27 一維等值理論模型下之剪力波速提升比與改良率之關係圖(剪力波 速比為 4.7，ρm / ρc = 0.85) 根據初步推論，在前述之高壓噴射灌漿場址之表面波震測結果是依循波速 等值速度模型之下限，則以所得之剪力波速提升率約 15%對照圖 2.27，可 得其改良率約為 18.5%，是為合理結果。 2.3.2.2 二維數值模擬結果 由簡單的一維理論模型可初步推測剪力波速提升率與現地改良率有關， 但並無法確定表面波震測所得之剪力波速剖面在均質化的過程中是依循何 種模型，因此，林志平等(2012)進一步採用 2D 數值模擬之方式進行探討。 數值模擬方法是由 Levander 所發展之四階速度應力錯置網格有限差分 法(Forth-Order, stress-velocity, staggered grid finite difference method, 續將以 速度應力差分法簡稱之)。在邊界的定義上，主要需處理的是位於模型外圍 的一個水平自由表面(free surface)與三個輻射邊界(radiation condition)。在水
平自由邊界上採用 Levander 所提出的映像法(image method)，在三個幅射邊 界上採用 Cerjan et al.所提出之吸收邊界(absorbing boundary)。在數值模型的 配置上，雖在邊界處設有吸收邊界，然而為能更加有效的降低因人造邊界 所造成之回波影響，震源與最後一個地音計的區段(study zone)與吸收邊界 中都隔有一段距離(extra zone)，此段距離是以地表地層之剪力波速估算，儘 量使得在收錄時間內，由震緣處發出的波不會回傳至地音計處被收錄。此 外，在震源的施加上是以高斯函數一階微分做為點震源模擬函數，以垂直 向壓力施加，而網格間距(h)與步進間距(Δt)的取用需根據 Levander 所提出 之假頻散現象及穩定度的影響。數值模型之配置與各參數定義如圖 2.28 所 示，各參數於模擬中之數值如表 2.1 所列。 圖 2.28 數值模型參數與配置示意圖
表 2.1 參數數值表 參數 數值 Nx (cells) 2000 Ny(cells) 1000 Nt(steps) 8192 h(m) 0.2 Δt(ms) 0.0515
Source Wavelet g(t)2texp(t2)
α 2500 b(cells) 60 數值模擬之地層配置參考高壓噴射灌漿場址案例如圖 2.29 所示，在二 維的條件下，共進行六種不同現地改良率之情形如表 2.2，改良深度自地表 下 5 公尺至 25 公尺，當改良樁樁徑 1.2 公尺，淨空間距在改良率 15%時為 6.8 公尺；而在材料參數的設定上，原土壤之剪力波速為 170 m/s，壓縮波 速 330 m/s，密度為 1.7 g/cm3，改良樁之剪力波速為 800 m/s，壓縮波速 1920 m/s，密度為 2 g/cm3，此中，改良樁與原土壤之剪力波速比為 4.7。另方面， 在表面波震測資料的收錄模擬上，採用 10 公尺之近站支距，24 個接收器， 接收器間距為 2 公尺，組成共 46 公尺長之測線。 表 2.2 地層模型配置數值表 Model GIR(%) D(m) W(m) h(m) A1 15 1.2 6.8 0.2 A2 20 1.2 4.8 0.1 A3 30 1.2 2.8 0.2 A4 40 1.2 1.8 0.2 A5 50 1.2 1.2 0.2 A6 60 1.2 0.8 0.2
圖 2.29 數值模擬地層模型配置示意圖 以改良率 30%為例，數值模擬所得時間域資料，並進行頻散曲線分析 之結果如圖 2.30a、b 所示，圖中細實線為原土壤頻散曲線之理論值，白色 圓圈為分析所得之結果，從圖 2.30b 之結果可發現，頻散曲線在低於 20Hz 處的相位波速有明顯的提升，將此頻散曲線結果採用美國 Kansas Geological Survey 所發展之 Surfseis v1.8 進行反算分析，於反算參數中設定反算層數為 2 層，地表至地表下 5 公尺為一層，超過地表下 5 公尺為半無限域空間，反 算結果如圖 2.30c 所示，在圖 2.30c 中可見，未改良的深度內其剪力波速剖 面與原土壤相同皆為 170m/s，而有改良的深度，在改良率為 30%的狀況下， 其剪力波速提升到 210m/s(剪力波速提升率為 23.4%，如圖 2.30d 所示)。 將表 2.2 所有模擬結果繪入剪力波速提升率與改良率之關係圖中(如圖 2.31 所示，剪力波速比為 4.7)，由圖中各個資料點可以發現，隨著改良率之 增加剪力波速提昇率也隨之增加，模擬結果顯示：表面波震測在高度異質 性地層中的均質化行為較為貼近波速等值速度行為之下限，當改良樁的勁 度相較未改良土壤高出許多時，改良樁所佔的體積比例控制了有效剪力波 速的量測，而非改良樁本身的勁度。然而二維模擬最大的癥結在於簡化改 良樁為牆狀，無法真實反應實際地盤改良場址。因此表面波震測所量測得
到的剪力波速提昇率某種程度反映地盤改良率，此假說是否可成立，還有 待商榷。 圖 2.30 表面波震測數值模擬分析結果：(a) 改良率 30%頻散曲線影像圖； (b) 改良率 30%頻散曲線影像圖；(c) 反算分析之剪力波速剖面；(d) 剪力 波速提升率 Phase Velocity (m/s) F re q u e n c y ( H z ) 100 150 200 250 20 30 40 50 60 70 10 20 30 40 50 60 0 0.1 0.2 0.3 0.4 0.5 0.6 Offset (m) T im e ( s e c ) 150 170 190 210 230 250 0 1 2 3 4 5 6 7
Shear Wave Velocity (m/s)
D e p th ( m ) -0.1 0 0.1 0.2 0.3 0 1 2 3 4 5 6 7 Improving Rate D e p th ( m ) (a) (b) (c) (d)
IR=30%
IR=30%
圖 2.31 二維模擬下(A1~A6)地盤改良率與剪力波速提昇率之關係圖
2.4 綜合評析
綜合以上文獻回顧，得知表面波震測應用於大範圍非壞性檢測極具潛 力，二維數值模擬之結果似乎可以解釋表面波於高度異質性土壤均質化之 過程，並得到地盤改良率與對應剪力波速提昇率之關係可能依循一維等效 波速理論下限值，但實際地層之改良樁並不是簡化牆狀，而是地層中含有 許多柱體，因此需進一步作三維數值模擬重新檢視此均質化行為；另一方 面表面波震測測線施作時也會收到從測線兩側傳來的訊號，目前二維模擬 無法探討此一情形，唯有在三維空間才可針對側向變化進行模擬規畫，因 此本研究採用三維數值模擬工具。第三章 研究方法
如圖 3.1 所示，本研究採用數值模擬之方式，主要透過不同的地層模型 期望了解表面波在高度異質性土體所量測得到的剪力波速均質化之過程， 嘗試界定表面波震測於側向變化地層可能取樣的範圍，最後提出表面波均 值化可能依循的準則，與現地資料結果進行比對探討，使得表面波震測適 用性更廣，成為具有代表性的地盤改良檢測方式之一。 圖 3.1 研究方法流程圖3.1 三維震波模擬方法與選用軟體
3.1.1 多頻道頻率波速轉換分析方法
本研究主要是利用多頻道波場轉換法中頻率波數轉換法進行資料的分 析(Lin and Chang, 2004)。其在頻散曲線分析上(圖 3.2)，主要是將震測資料 由時間-空間域(t-x domain)以離散傅立葉轉換轉換至頻率-空間域(f-x
domain)，再經離散空間域傅立葉轉換(discrete-space Fourier Transform)將頻 率-空間域轉換至頻率-波數域(f-k domain, Prokis and Manolakis, 1992)，再透 過變數轉換而獲得頻率-波速域(f-v domain)的結果。在實際操作中，由於在 時間上的取樣點較為豐富，時間與頻率的轉換可以經由較為有效率的快速 傅立葉轉換來進行，而空間與波數的轉換上，依舊維持。由數學運算來看， 首先得到將多處接收器位置之二維震測資料 u(xn,tm)進行快速傅立葉轉換： U( i,xn)= 1 M u(tm,xn)exp(-j2πfntm) M-1 m=0 (3.1) 得到離散的頻率-空間域資料U( i,xn)。其中，M 是在時間軸上以∆t 為取樣時 間間距的離散數目，亦即在時間軸上所取的紀錄點數；n 表示接收器位置的 索引號； fi= i (M-1)∆t，表示頻率。而後再將U( i,xn)進行離散空間域傅立葉轉 換， Y( i,k)= U( i,xn)exp(-jkxn) N-1 n=0 (3.2) 得到在波數上連續的頻率-波數域資料Y( i,k)。其中，N 是接收器的數目。 又已知波速v=2πf k ，則欲得到頻率-波速域的結果可將(3.2)式改寫做 Y( i,v)= U( i,xn)exp(-j2πfi v xn) N-1 n=0 (3.3) 而於頻率-波速域中，單一頻率之頻譜波峰值即為其所對應之雷力波波速。
圖 3.2 多頻道波場轉換法分析流程：(a) 時間域資料(t-x domain)；(b) 頻率 域資料(f-x domain)；(c) 空間域離散化至波數域；(d) 頻散曲線
3.1.2 三維震波模擬軟體 SPECFEM3D 架構說明
本研究關於前處理建立模型的部分，軟體主要採用 Cubit 13.1，由聖地 亞國家實驗室(Sandia National Laboratories)於 1990 年初開發，開發者為 Ted Blacker 與 Michael Stephenson。Cubit 是一個高階的六面體網格快速生成工
(a) (b)
具，適用於 CFD 和 FEA 分析所需的結構化網格。其好處是它能夠從多種 CAD 軟體導入幾何體再加以劃分網格，並擁有直觀的圖形界面環境(如圖 3.3)，輔以大量的命令語言命令行界面。
圖 3.3 Cubit 圖形化使用介面
SPECFEM3D 這套模擬三維地震波之軟體最早是從 1998 年開始，由 Dimitri Komatitsch 和 Jeroen Tromp 兩位學者在美國哈佛大學和加州理工學 院發展出來，基於 1995 年至 1997 年由 Dimitri Komatitsch 和 Jean-Pierre Vilotte 在法國巴黎地球物理研究院早期之研究，根據 spectral element method 撰寫出 Fortran 90 程式碼，具有可移植性，且支援 Linux 系統平行運算。 軟體架構主要如圖 3.4 所示，Cubit 主要功能是建立模型、劃分網格、 給定材料參數並加上邊界條件；Scotch 可將模型劃分為較小的塊體，透過 Mass Passing interface(MPI)分配給 CPU 內各個核心去處理工作；
xgenerate_databases 將分析所需的格點與特定模型參數之資訊建立後，最後 再由 xspecfem3D 進行運算求解。
圖 3.4 SPECFEM3D 模擬流程圖 對於波譜元素模擬首要關鍵便是三維模型需有高質量的網格，主要受 到以下幾點的限制： (1)每一最小波長之格點數量(2)數值收斂穩定條件(3)數值計算代價與計算 機運算效能之平衡。 質量差的網格可能會產生諸多數值問題例如增加計算上的耗費，模擬缺乏 收斂性，或是不正確的結果。SEM 有一特性影響網格的建立是用於離散波 場之多項式階次 n，如果 n < 4 不精確性與標準 FEM 相似；當 n > 10 精確 性可以提昇但是數值模擬計算之代價也相對變得提高。n 之選擇與網格間距 Δh 也有關：每一最小波長的點數須等於或大於 5，也就是(n+1)λmin/Δh≧5。 網格設計另一限制條件為穩定的時間計畫，步進時間有一計算穩定之上限 值，定義庫倫穩定數：C =Δt(v/Δh)max，Δt 選擇的步進時間，(v/Δh)max壓縮 波速與網格間距之最大比值，此庫倫收斂穩定條件庫倫穩定數不應超過
![圖 2.19 Lagrange 多項式於[-1, 1]內插，n L =4，有 n L +1=5 個控制點(Komatitsch et al., 2005)](https://thumb-ap.123doks.com/thumbv2/9libinfo/8458807.183082/36.892.115.808.392.1130/圖2Lagrange多項式於內插nL=4有nL+1=5個控制點Komatitschetal2.webp)
