2.4 P-SV 震波模擬方法
2.4.1 解析方法
2.4.1.1 射線理論
0 , , (t r z=
u 可單純表示為時間以及地表位置與震源水平距離
r
的函數。對於 地層位移函數u(t,r,z)，可透過雙轉換(Double Transform)
轉換至頻率-
慢度域) , , (w p z
U 。我們可將時間
-
空間域上的位移u(t,r,z)表示為：∫
∫
∞−
∞
∞
−
⋅
⋅
=
0
0( ) ( , , )
) exp(
) , ,
(t r z dw i t dp J pr U p z
u ω ω ω
(2. 104)
其中，
U
為頻率-
慢度域中的地層位移函數；ω為頻率；p為慢度；J0是一 般使用的貝索函數(Bessel Function)
。對於層狀地層的位移情形的解析解來 說，處理的技術差異主要就在於如何解出U
的控制方程式解及如何處理不 同域的轉換。從式((2. 104)
中可知，時間-
空間域中的地層位移情形可以透過 對頻率-
慢度域中的位移函數積分獲得。因此，對於函數U
在慢度/
頻率的複 數平面上性質的了解是很重要的，這些特性將牽涉欲獲得時間-
空間域地層 位移情形時應如何劃分所需的積分路徑，及在積分過程中需做哪些的假 設，而這些不同的劃分及假設將會影響這些理論方法所能適用的範圍(Clarke, 1989)
。2.4.1.1 射線理論射線理論射線理論射線理論
在震波模擬上，有一派學者採用射線理論作為基礎，認為地層的反應 是由於波的射線所造成，其最基本假設便是此波的頻率極大。一般常見的 有廣義射線理論、全波理論以及ＷＫＢＪ理論。
廣義射線理論 廣義射線理論 廣義射線理論
廣義射線理論：完整的頻率
-
慢度域地層反應U
可視為是在地層中以不同能量管道
(Route)
傳遞的射線路徑(Raypath)
所造成的地層反應總和，而這些地層反應是由在地層中傳遞的
P
波、S
波及其在介面上的折射及反射所造成的結果。廣義射線理論是藉由清楚分離出不同射線路徑的到達波
(Isolated
Arrival)
，然後逐一將這些波對同一位置的影響作加總而得到模擬的地層反應。如圖
2.20
所示，假設有二組以不同慢度p
前進的射線，因為不同的p
而會有不同的路徑產生。對於A
點來說，造成其位移的除了路徑1
的反射 波外尚有路徑2
的折射波，將這兩種波的影響皆考慮便是廣義射線理論的 基本精神。圖 圖 圖
圖 2.2.2.2.202020 20 廣廣廣廣義射線理論示意圖義射線理論示意圖義射線理論示意圖義射線理論示意圖
1968
年，Helmberger
藉由Cagniard
所發展而由de Hoop(1960)
修正的理論，採用
Laplace
轉換將時間-
空間域的位移函數轉換至頻率-
慢度域，並對慢度
p
進行變數變換，使得原本由零積至無限大的積分變成在複數平面上 有限的封閉曲限，亦即Cagniard
積分路徑的積分。如此一來，在頻率及慢 度上的積分將能互相抵消而只剩下時間上的積分。在這技術中最困難的要屬將
Cagniard
積分路徑找出，因為這積分路徑是以非顯性的方程式呈現，並且會因為不同的射線路徑及不同的計算位置而不同，常需耗費許多時間
(Clarke, 1989)
。此方法對於高頻及計算相位角上有極好的結果，但對於有較多 水 平 地 層 的 狀 況 或 是 須 計 算 遠 距 離 處 的 位 移 便 不 適 合
(Wang and Herrmann, 1980)
。全波理論 全波理論 全波理論
全波理論：雖然廣義射線理論已經納入了折射波與反射波的影響，然而在 實際上，地層的震動情形並不僅止於折射波與反射波的影響，其尚有波前
或是沿著不連續介面前進的波等非地層幾何不連續
(nongeometic effect)
所造 成的影響。為能將這些影響包括於理論中，全波理論直接對慢度域及頻率 域積分，以近似法轉換射線理論(Ray Theory)
中頻率-
慢度域上的地層反應函 數後，將原本是實數慢度的積分分解為在複數慢度平面上的積分，並以劇 降近似(steepest descent approximation)
技巧解決鞍點(saddle point)
上的積分 困難，維持頻率上的積分不變(Kennett, 1983)
，而得到最後時間-空間域上的 地層反應。這樣的方式使得其可處理鞍點處的震波影響，進而增加了模擬 震波的準確性。WKBJ 理論理論理論理論：
WKBJ
理論是由Chapman
在1978
年時提出在有低速帶夾層 或在caustic
處(
請參考Aki and Richards, 2002)
時，射線理論會無法計算或得 到 錯 誤 的地 層反應 的 修 正方 法。此 理 論 名稱 的由來 是Chapman
使 用Wentzel-Kramers-Brillouin-Jeffreys
四人所提出的理論，得到頻率-
慢度域中U
的近似值，先對頻率做積分，再在複數慢度平面上延著實部軸進行慢度的 積分。此方法值得一提的是它可自然的對地層反應進行平滑化的動作，而使得
caustic
的存在不成問題。然而與全波理論比起，此理論較為無法考慮頻率相關的問題，但可有較為快速的計算
(Muller, 1985)
。 2.4.1.2 空間頻率空間頻率空間頻率空間頻率/慢度積分法慢度積分法慢度積分法慢度積分法概念上，這一類的方法是將進行傅立葉轉換後得到的頻率
-
慢度域上的U
以積分的方式得到時間-
空間域的 u。這樣的作法似乎與上述方法相同，然而其差別在於，上述三個方法的理論基礎都是建立在射線理論
(Ray
Theory)
上，亦即考慮波場中頻率極大時的情形，僅考慮較為高頻的實體波P
波與S
波的影響。但此類方式的理論基礎不將P
波與S
波的影響抽離出 來看待，不僅考量頻率趨近於無限大時的地層反應，亦考慮了其他頻率的 狀態而得到完整的震波模擬地層反應。這類方法中較為人所熟悉並常被使用的是反射係數法
(Reflectivity Method, Fuchs, 1968
；Fuchs and Muller, 1971
；Kennett and Kerry, 1979
；Kennett, 1983)
。反射係數法的主要精神是將頻率-
慢度域中的U
分解成向上行進及向下行進兩部份，藉由使用傳遞矩陣
(Propagator Matrix)
的概念，使 其符合層面間應力須連續的邊界定義。重複進行此動作直至將每層的反射 訊息獲得，而後使用快速傅立葉轉換及數值漢克爾轉換(Numerical Hankel
Transform)
得到前述的雙轉換結果。一開始所發展的反射係數法
(Fuchs, 1968
；Fuchs and Muller, 1971)
只侷 限於無自由表面存在且震源在反射係數收集區的另一側的地層模型，而在 慢度上的積分也只能是慢度對應在入射波的角度為實數時的範圍(Clark, 1989)
。經過Kennett and Kelly (1979)
及Kennett (1983)
對此方法的改良，其 將所要計算的反射係數分成震源以上以及震源以下兩部份進行計算，而後 分開積分加總，使得自由表面存在並且任意震源位置的地層模型震波模擬 得以實現。此外，在這基礎上只需將在慢度的積分範圍擴大至基態表面波 的慢度處，便可獲得完整的震波模擬(Clark, 1989)
。2.4.1.3 震態疊加法震態疊加法震態疊加法震態疊加法
震態疊加法的要義是將地層反應視為所有存在的表面波震態反應的總 和
(Takeuchi and Saito, 1972)
， 藉 由 將 在 頻 率-
空 間 頻 率 域(
空 間 頻 率(wavenumber)
與頻率及慢度有關，頻率-
空間頻率域是頻率-
慢度域的一體兩面，僅是變數的變換
)
中的singular
位置(
即表面波頻散曲線的位置)
找出，我 們便可使用殘數定理(Residue Theorem)
將頻率-
空間頻率域中的空間頻率積 分改以殘數數列(Residue Series)
表示。當有足夠數目的數列加總在一起，我 們便可得到完整的地層反應(Harvey, 1981)
。震態疊加法在起初的發展中，採用
Tomson and Haskell
方法求得pole
的位置，然而此方法在採用數值計算時，在高頻區會有精準度的問題，這 使得早期的震態疊加法僅適用於取得頻率小於
1Hz
的震波反應。對於Tomson and Haskell
方法在高頻會有精準度的問題，現已有許多處理方式可解決
(
請參閱2.1.2
節所列參考資料)
，最常為使用的是Knopoff’s Method (Knopoff, 1964
；Randall, 1967
；Schwab, 1970
；Schwab and Knopoff, 1972)
及Delta Matrices (Pestel and Leckie, 1963
；Thrower, 1965
；Dunkin, 1965
；Watson, 1970)
。如此一來即使至10000Hz
部份的pole
位置亦能準確的求出(Schwab et
al., 1984
)，使得震態疊加法模擬出的地層反應更加準確。到此，震態疊加法已將在複數空間頻率平面上積分時
pole
對地層反應 的影響考慮，然而，如圖2.21
，此圖是表示欲取得在頻率/
空間頻率域中完 整的地層的總反應時在空間頻率複數平面上的積分情形，其不僅要考慮複 數空間頻率上pole
的影響，亦應要同時考慮branch cut
、branch point
的部 分，若不考慮此二處對地層反應的影響，將會漏失掉反應初期Ｐ波或是Ｓ 波的震動，而無法獲得全部的地層反應。Harvey
於1981
年提出了震態鎖近 似法(Locked Mode Approxiation)
，避開在branch cut
的積分，而改以震態疊 加的方式處理。其作法是在極深的位置置入一極高速地層，使得原本在Reimann
面上的遺漏震態(Leaking Mode, Gilbert, 1964
；Aki and Richard, 2002)
能在空間頻率的實數軸上出現，再以震態疊加的方式計算其對地層的影 響，如此便可獲得完整的地層反應模擬(Harvey, 1981)
。在使用此方法時有 些許問題需要注意，例如需小心放置高速地層的深度，不讓其介面反射回 地表的波干擾了原本地層的反應；此外，在尋找由於使用震態鎖近似法而 鎖住的遺漏震態需要耗費極大的時間，使其極其相鄰的pole
能被精準的找出
(Wang, 1981)
，才可有較好的結果。圖 圖 圖
圖 2.2.2.2.22221111
複數空間頻率平面積分路徑示意圖複數空間頻率平面積分路徑示意圖複數空間頻率平面積分路徑示意圖複數空間頻率平面積分路徑示意圖 2.4.2 數值方法數值方法數值方法數值方法
震測資料的數值模擬目前已被廣泛的用來協助現地震測資料的解讀、
提供震測試驗訊號分析和研究試驗參數設定所需的人造震測資料，以及對 震波在非均質、非均向地層傳遞行為的瞭解。經過二、三十年的演進，目 前 已 有 許 多 不 同 的 數 值 方 法 可 供 選 用 。 除 了 廣 為 人 知 的 有 限 差 分 法
（（
（（
Alterman and Karal, 1968
；Boore, 1972
；Kelly et al., 1976
；Virieux, 1986
；Levander, 1988
；Takeuchi and Geller, 2000
；Zhang, 2004
等）））、有限元素法）（
（
（
（
Lysmer and Drake, 1972
；Schlue, 1979
；Chen, 198
4；Kay and Krebes, 1999
）））） 外，尚有擬頻譜法（pseudo-spectral method
）（（（（Gazdag, 1973
；Kosloff and Baysal, 1982
；Johnson, 1984
；Reshef et al., 1988a, 1988b
；Huang and Yeh, 1991
等）））、） 混合數值法（hybrid method
）（（（（Shtivelman, 1985
；Kummer, Behle & Dorau, 1987
；Ven den Burg, 1984
；Emmerich, 1992
；Moczo et al., 1997
等））））等技術。除此之外，許多合於上述數值方法概念但有不同操作上的應用變形，使得
數值模擬的成果越加豐碩。在下面三個小節中，作者將對有限差分法、有 限元素法及擬頻譜法在震波模擬使用上作回顧並簡單介紹各方法的基本概 念。
2.4.2.1 有限差分法有限差分法有限差分法有限差分法
有限差分法簡單的說便是在離散的向度上，將控制方程式及邊界狀態 中的微分項以差分式替代
(Boore, 1972)
，使得原先的微分方程式變成簡單的 四則運算式(Mocoz, Kristek and Halada, 2004)
。當我們要將有限差分法應用 至特定的控制方程式或微分方程式中時，總括來看有三個主要步驟：先(1)
建構離散的有限差分離散模型，亦即建構離散的網格，並對微分方程式導 入有限差分式。然後對導出的有限差分運算式進行(2)
有限差分離散模型的有限差分法簡單的說便是在離散的向度上，將控制方程式及邊界狀態 中的微分項以差分式替代