研究方法

如圖

3.1

所示，本研究將對震源產生的地層表面擾動情形進行數值模

圖 圖 圖

圖 3.13.13.13.1 研究方法流程圖研究方法流程圖研究方法流程圖研究方法流程圖

擬。以已知的地層情況，製造震源穩定且無雜訊等可能干擾因素存在的震 測資料，而後觀察在不同地層況狀下以傳統測線施測以及以

pseudo-section

概念施測的震測資料，並以多頻道波譜轉換法對其進行訊號分析，以求得 雷利波的頻散曲線。在已知的地層情況中，均質均向半無限域地層以及水 平層狀地層中的頻散曲線具有解析解

(

詳見

2.1.2

節

)

，因此，我們將以解析 解的結果與數值法結果做一個比較，討論施測參數對頻散曲線的影響，並 以此處的結果作為經驗值，避免對具有側向變化地層的數值模擬結果做過 度的解讀。而對於兩種不同施測方式結果的震測資料 ，將以時間

/

空間域的 表面擾動情形以及在頻率

/

波數域上的頻散曲線差異為討論重心，避免牽涉

需反算過程的震測地層波速剖面結果。本研究中的數值模擬是以Ｍ

athworks

公司所發展的

Matlab 6.5.1

為程式撰寫引擎。

震測模擬震測模擬

震測模擬震測模擬：在本研究中，由於欲討論具有側向變化的土層，因此震測模擬 將採用由

Levander(1988)

在

Madariaga(1976)

、

Virieux(1984;1985)

使用的速度

應力

(velocity-stress)

有限差分法架構下所發展的四級錯置網格有限差分法

(fourth-order staggered-grid finite-difference method)

作為運算式

(

以下將以速 度應力差分法表示

)

。

速度應力差分法其運算式與傳統解決波傳方式不同，其不直接對二次 微分的微分方程進行差分動作，而是對兩組一次微分方程進行差分；在網 格上採取錯置網格

(staggered-grid)

，讓不同的物理量不在同一位置上，而可 使下一時間中，空間上一點之物理量位於這一時間點空間物理量的中心，

美好的符合差分式所需

(

其細節將在

3.1.1

節中加以說明

)

。在時間上採用中 央差分式

(central difference approximation)

，但由於其錯置網格的定義，使其 為顯性解法

(

亦即其下一時間的結果可由這一時間決定

)

。此方法除具有容易 程式化、高記算效率以及由於有限差分法在計算上的區域性

(loacal)

而不需 龐大的記憶體等一般有限差分法所具有的優點外，還具有

(1)

僅需處理一次 微分方程的單純性

(2)

較不受地層物質卜松比

(Poission’s ratio)

而影響計算的 容許性

(3)

較不因格網而產生嚴重人造頻散現象的正確性

(4)

可依情況使用應 力式震源或速度式震源的選擇性；

(5)

穩定並正確定義自由表面邊界狀態的 容易性；

(6)

可依需求再提高有限差分運算式精度級數

(order)

的擴充性；

(7)

可因應不同計算特性的程式而採用純量式、向量式或平行電腦處理的多元 性（（（（

Levander, 1988;Yao and Margrave, 2000

））））等優點。

表面波震測法 表面波震測法 表面波震測法

表面波震測法：所欲模擬的表面波震測法是採用多頻道波場轉換法。多頻 道波場轉換法與僅使用兩個接受器的波譜分析法比起，

1)

不需多次變換間

距而可有較好的頻散曲線頻寬；

2)

可避開進行摺開的動作而無須擔心其造

) 複雜的二級海格波利克系統

(second-order hyperbolic system)

，因此，在速度 應力差分法中並不直接對此兩組控制方程進行差分運算元應用，反之，放

錯置網格架構 錯置網格架構 錯置網格架構

錯置網格架構：：：在有限差分法的使用上需先將問題討論的時間、空間進行： 離散。如前所述，速度應力差分法採用錯置的網格將不同的物理量定義在 不同的離散點上，使物理量定義處能自然符合運算式的需求

(

為使差分運算 式中物理量表示易為閱讀，且讓離散後的物理量與理論上的物理量有所區 分，將使用U,V,Σ ,,Τ Ξ 分別代表 v_x,v_z,τ_xx,τ_zz,τ_xz)。在時間上，速度

(U

、

V)

與 應力

(

Σ、Ｔ、Ξ

)

並不定義在同一時間點上，而是將其錯開相差，亦即將速

度定義在

k+1/2

時間點上，而將應力定義在

k+1

時間點上

(

如圖

3.2

及圖

3.3

所示

)

。對空間的離散上，速度

U

與

V

並不在同一點：

U

是定義在格點，即 點

i,j

上，而

V

則是平移錯開至四個相鄰格點的中央，即點

i+1/2,j+1/2

處

(

如 圖

3.2

所示

)

；應力Σ、Ｔ定義在

i+1/2,j

的位置，而應力Ξ定義在

i,j+1/2

位 置上

(

如圖

3.3

所示

)

。除了會隨著時間而有所變化的速度與應力兩個物理量 外，材料的性質

(

λ,µ,ρ

)

是不變的，因而不需在時間上做離散的動作，然而 在空間上的離散卻是避免不了。此錯置網格架構將λ定義在

i+1/2,j

的位 置；將μ同時定義在

i+1/2,j

及

i,j+1/2

處；亦將ρ同時定義在

i,j

及

i+1/2,j+1/2

兩個位置上

(

如圖

2.23

及圖

2.24

所示

)

。這樣的網格架構將可使得下一時間 物理量能自然在這一時間物理量的中央位置，而使四級有限差分式可以應 用。詳細的運算關係將在下一段做介紹。

四級錯置網格有限差分法運算式 四級錯置網格有限差分法運算式 四級錯置網格有限差分法運算式

四級錯置網格有限差分法運算式：：：：在將有限差分法導入

(3. 2

式時，需要空 間上以及時間上的有限差分運算元，亦即將微分式以差分式替換的轉換方 式。四級錯置網格有限差分法在時間的一次微分上是採用標準的二級中央 差分式

(second –order central difference)

，其運算元如下：

) 1 ( _{+1/2 −1/2}

∆ −

= ^{k k}

k

t a a

a t

D

(3. 3)

其中，D_t⁺為對時間一次微分的差分運算元，用以逼近控制方程中的



穩定性與假頻散問題的控制

制。根據

Levander(1988)

對本研究中速度應力差分法所提出的建議是：

穩定性控制 的一個水平自由表面

(free surface)

與三個輻射邊界

(radiation condition)

。 在錯置網格的架構下要正確且穩定的定義自由表面在垂直面上沒有應 力的狀態是很容易的((((

Levander, 1988

))))。在自由表面的位置，其邊界條件是

0

0 )

( ₀

0 =

∂ +∂

∂

= ∂ ₌

= z

x z z

xz x

u z µ u τ

為了讓離散後的地層模型能符合

(3. 14

式，我們需在自由表面處及其上方位 置增加格點，對格點在自由表面上下的應力值予以限制，使其以自由表面 為中心，而成反對稱的值((((

Levander, 1988,

又稱作虛擬影像法

(image method),

如圖

3.2

所示))))。亦即，若假設自由表面處是網格中

z

方向上代號數為

0( j = 0 )

，則

0 =0

T ，T₋₁ =−T₁

2 / 1 2 /

1 =−Ξ

Ξ ，Ξ_−3/2 =−Ξ_3/2

(3. 14)

此外，為解決Ｕ、Ｗ在自由表面處上方的值的問題，採用

Robertsson(1996)

所建議的，將自由表面以上的Ｕ、

V

皆設為零

(

處理自由表面的詳細計算式 請見附錄三

)

。

圖圖

圖圖 3.23.23.23.2 自由自由自由自由表面邊界設定示意圖表面邊界設定示意圖表面邊界設定示意圖表面邊界設定示意圖

而至於位於左側、右側以及底部的邊界，在輻射的情況下，不會有反 射面存在，且沒有能量會在這裡擾動，但若只是將邊界的速度值及應力值 設為零，會因為有限的模型大小而無法真實模擬無線遠的邊界，而造成虛

擬的反射面存在。為了達到無限遠的邊界定義，有許多不同種類的吸收邊 界

(absorbing boundary condition, Cerjan et al., 1985

；

Sochacki et al., 1987

；

Chew and Liu, 1996

；

Festa and Nielsen, 2003)

或是傳導邊界

(transmitting boundary condition, Clayton and Engquist, 1977

；

Reynolds, 1978

；

Emmerman and Stephen, 1983

；

Higdon, 1990,1991,1992,1994

；

Zhu, 1999)

已被發展用以 消除在此三個邊界的反射情形。

在此，我們將採用較易為使用的吸收邊界((((

Cerjan et al., 1985

))))，沿著三 個邊界建立起三個寬帶，使得波在此寬帶內傳遞會有遞減的情形

(

如圖

3.3

所示

)

。將這樣的概念應用在此的方法是建立起一個權重矩陣

W(x,z)

，規範

圖 圖 圖

圖 3.3.3.3333 3. 吸收邊界示意圖吸收邊界示意圖吸收邊界示意圖吸收邊界示意圖

著計算出的速度值應以多少的比例減少，越接近寬帶外緣，其值越接近零；

越接近寬帶內緣，其值越接近

1

。其數學式為：

W(x,z) =

[

( )

]

²

exp− a b−x

, for x< b or x > Lx - b

[

^{( )}

]

²

exp− a b−z

, for z > Lz - b 1,

其他位置

(3. 15)

其中，

a

是可控制遞減速度的係數，數字越大遞減越快接近零

(

一般盡量不 使其小於

0.92)

；

b

是控制遞減帶寬的係數

(Cerjan et al.(1985)

建議至少需

20)

；

Lx

、

Lz

分別是模型中

x

方向和

z

方向的大小。

震源設定 震源設定 震源設定

震源設定：：：：震源的形式有許多不同種類，例如有爆炸型、剪力型、水平或 垂直型的點震源等等，在震源的設定上會以不同的應力或速度值又或是震 源所分佈在格點上的多寡或是影響範圍作為不同震源類型的區隔((((

Aboudi, 1971

))))。常見用於數值模擬的衝擊式震源

(impulse)

函數有三種：

高斯函數

(Gaussian function)

) ) ( exp(

)

(t t t₀ ²

g = −α_s −

(3. 16)

高斯函數一階微分 ^g(t)=−2α_s^(t−t₀^)exp(−α_s^(t−t₀⁾²⁾

(3. 17)

高斯函數二階微分 g(t)=(4α_s²(t−t₀)²exp(−α_s(t−t₀)²)

(3. 18)

其中，α_s是控制衝擊式震源到達最大位移所需的時間，其值越大，越快到 達

(

如圖

3.4

所示

)

，亦即表示其所包含的高頻波能量將較α_s小的能量增加；

t0是指震源開始作用的時間，在此之前皆為零。

圖 圖圖

圖 3.3.3.3.4444 衝擊式震源衝擊式震源衝擊式震源衝擊式震源α_s參數的影響參數的影響參數的影響參數的影響

αs

=2500

αs

=50

) ) ( exp(

) ( 2 )

(t t t₀ t t₀ ²

g =−α_s − −α_s −

對於點震源來說，可依目的的不同而有不同置入模型的方式，大致上 有壓力式震源

(pressure source)

、剪力波式震源

(shear wave source)

以及正向應 力

/

速度式震源

(normal stress/velocity source)

((((

Yao and Margrave, 2000

))))，不同 震源方式所擾動的物理量值如下：

壓力式 Σ_i,_j =g⁽t^),Τ_i,_j =g⁽t⁾

(3. 19a)

剪力波式 ^U_i,j ^=g(t),U_i,j−₁ ^=−g(t)， V_i,j = g(t),V_i+1,j =g(t)

(3. 19b)

正向應力/速度式 T_i,_j = g⁽t⁾或^V_i,_j =g⁽t⁾

(3. 19c)

本次研究中採用高斯函數一階微分為正向應力的衝擊式震源函數，函數的 控制參數取為α_s =2500, t₀ =0。

程式計算流程 程式計算流程 程式計算流程

程式計算流程 程式的計算流程非常的簡單，大抵上是由六個動作完成：

(1)

定義地層參數及程式控制參數；

(2)

檢驗離散後的穩定性及假頻散現象的控 制；

(3)

激發震源；

(4)

以速度應力差分法運算式計算震波傳遞；

(5)

套用邊界 定義；

(6)

將結果記錄下輸出

(

如圖

3.5

所示

)

。

圖圖圖

圖 3.53.53.53.5

程式計算流程圖程式計算流程圖 程式計算流程圖程式計算流程圖

3.1.2 傳統測線與傳統測線與傳統測線與傳統測線與 pseudo-section 測線測線測線施測測線施測施測施測

如圖

3.6

所示，由

2.2.1.2

節處已知多頻道表面波震測法的施測參數有 三個：近站支距x₀、接收器間距∆x以及測線長度L。已知測線長度與接收 器數目N 有關：L=(N −1)∆x，為了比較上的便利，在後續模擬中的參數說 明將以接收器數目取代測線長度。又由於在具有側向變化的地層上，我們 關心地層側向變化狀況與所得到的模擬結果關係，因此增加震源位置x_s的 表示。則，一個以傳統測線施測的模擬中，表達其施測參數的會有四個：

震源位置x_s、近站支距x₀、接收器間距∆x以及接收器數目N。

圖圖圖

圖 3.3.3.3.6666 MASWMASWMASWMASW 傳傳傳統測線施測示意圖傳統測線施測示意圖統測線施測示意圖 統測線施測示意圖

對於

pseudo-section

測線而言

(

如圖

3.7

所示

)

，其需要改變震源位置，

以得到其震測結果。就

pseudo-section

概念而言，不同震源位置並無一定要 求，只需其拓展的展距必須與原展距有所重疊，在此，本研究採用在擴展 的展距中，只有一個接收器重疊，亦即，每一次的震源位置會改變一個展 距長。例如第一次的震源位置為x_c，展距長L=(N−1)∆x，則下一次震源的

以得到其震測結果。就

pseudo-section

概念而言，不同震源位置並無一定要 求，只需其拓展的展距必須與原展距有所重疊，在此，本研究採用在擴展 的展距中，只有一個接收器重疊，亦即，每一次的震源位置會改變一個展 距長。例如第一次的震源位置為x_c，展距長L=(N−1)∆x，則下一次震源的

在文檔中 Pseudo-section概念於表面波震測應用之數值模擬探討 (頁 85-0)