周成峰

(福建省建筑設計研究院有限公司 福建福州 350001)

0 引言

目前建立在經典彈性波理論基礎之上的反射波法只適合研究固相、液相或氣相等單相介質中地震波傳播規(guī)律，而巖石裂隙發(fā)育的富水區(qū)及大部分土是由固相和流相(包括氣相和液相)組成的雙相或多相介質[1]。與單相介質不同，地震波在雙相或多相介質中傳播時有其特殊的傳播規(guī)律及波動響應特征。Biot介質是一種典型的雙相介質，高階交錯網格波場數值模擬技術可有效地用來研究Biot介質波動響應特征。

Biot介質高階交錯網格波場數值模擬是用有限的離散值來代替連續(xù)介質，離散代替連續(xù)在頻率上總會存在計算誤差，頻率上的計算誤差會導致地震波具有不同的相速度和群速度，即頻散現象，其中數值頻散包括空間頻散和時間頻散。數值頻散嚴重干擾地震波場，極大影響數值模擬的精度和分辨率[2]。消除數值頻散的方法和原則主要有：

(1)在兼顧計算量的前提下盡可能減少空間網格間距；

(2)當空間網格間距一定時，在保證空間分辨率的情況下盡可能降低子波主頻；

(3)在保證穩(wěn)定性的前提下，減小時間步長；

(4)適當提高時間和空間差分階數，從而提高模擬精度，壓制數值頻散；

(5)對差分算子進行校正，補償離散化引起的誤差，但是利用差分算子校正進行壓制數值頻散會增加一定的計算量，在時間域遞推顯格式的差分算法有一定難度[3]；

(6)使用通量傳輸校正(FCT)方法平滑波場，壓制數值頻散[4]。

在使用FCT方法時會損失一定的能量和頻率，并需要更多的內存量，但FCT方法可以用較低階的差分方程就可達到較高精度的正演模擬，大大提高運算效率。同時，FCT方法可以放大時間和空間步長，從而抵消FCT帶來的計算量增加，因此FCT是很好的一種消除數值頻散方法[5]。

本文將在不考慮耗散的二維Biot各向同性介質中，運用Matlab編程對交錯網格差分格式具有空間四階精度和時間二階精度的數值頻散問題和FCT校正方法進行分析、探討。

1 頻散計算及FCT校正

1.1 空間與時間頻散計算

為了考察網格離散對波動傳播速度的影響，引入式(1)的無量綱物理量：

(1)

c為彈性波在真實連續(xù)介質中的傳播速度；

物理量q可以反映數值速度與真實速度的關系。

為了方便公式推導和分析，引入無量綱參數ζ和H：

(2)

式中λ為波長；

Δt為時間步長；

Δx為x方向的空間步長；

Δz為z方向的空間步長。

(3)

其中：

Ai=cisin(πHicosγi)+c2sin(3πHicosγi)，i=x,z。

則γx、γz為平面波傳播方向與坐標軸x、z所成夾角，c1、c2為空間四階差分權系數。

當x和z方向的空間步長Δx=Δz=Δ時，引入控制數值頻散的無量綱參數ζ=VP-solidΔt/Δ和網格尺寸無量綱參數H=Δ/λ，結合不考慮耗散的二維Biot各向同性介質中三類波的速度表達式及式(3)可推得規(guī)則網格快縱波數值頻散qp-fast,慢縱波數值頻散qp-slow和橫波數值頻散qs：

(4)

(5)

(6)

其中：

Ax=c1sin(πHcosγx)+c2sin(3πHcosγx)

Az=c1sin(πHcosγz)+c2sin(3πHcosγz)

式(4)～式(6)中，vshear-solid為骨架橫波速度；

vp-solid為骨架縱波速度；

vf為孔隙流體縱波速度；

ρs、ρf分別為固體基質(固相)密度和孔隙流體(流相)密度；

密度比ξ定義為ρf/ρs；

φ為孔隙度；

彎曲度τ是獨立于固體和流體密度的幾何度量；

γ是與流體運動中骨架的微觀模型有關的因子，當流體中的固體顆粒為球形時，γ=1/2。

1.2 FCT校正方法

FCT方法基本原理是假設所有的極值點都是由數值頻散引起的，對所有網格點進行漫射校正處理，達到消除數值頻散的效果。再對非局部極值點進行補償的反漫射校正，使任何不需要漫射的地方得以恢復。實際上，漫射校正是一個線性的平滑過程，反漫射校正由于是有選擇的，所以是一個非線性的過程[6]。對于交錯網格有限差分速度-應力波動方程，應用FCT方法對差分數值解進行校正的過程中，可以只對速度進行計算，也可以只對應力進行校正，因為在完成對速度或應力的校正后，通過交錯網格有限差分速度-應力波動方程對應力或速度進行了校正。顯然，兩種方法都需對4個參數進行校正，運算量是相同的，本文選擇對速度進行校正。

利用FCT方法進行校正，廣義地說主要有3步：交錯網格差分計算、通量漫射校正(漫射因子η1)和通量反漫射校正(反漫射因子η2)，具體執(zhí)行步驟詳見參考文獻[7]。

2 頻散分析

在確保數值計算穩(wěn)定的前提下，數值頻散及FCT校正方法分析的基本流程如圖1所示。

圖1 數值頻散及FCT校正方法分析基本流程

2.1 空間頻散分析

應用Matlab編程，計算式(4)～式(6)可得，空間差分網格離散引起的快縱波(快P波)、慢縱波(慢P波)和橫波的數值頻散曲線，如圖2～圖4所示。其中，橫坐標為單位波長內網格點數控制參數H，縱坐標為數值頻散qp-fast、qp-slow和qs。在二維平面內考慮到Ax、Az關于角度γx、γz取值的對稱性，考察3個不同的傳播方向：沿x軸(γx=0°，γz=90°)，沿xz平面三分線(γx=30°，γz=60°)，沿xz平面對角線(γx=45°，γz=45°)。相關參數取值如表1所示。

表1 空間頻散計算參數

其中，vp-solid、vshear-solid和vf的單位為m/s。

圖2 空間四階有限差分快P波數值頻散曲線

圖3 空間四階有限差分慢P波數值頻散曲線

圖4 空間四階有限差分橫波數值頻散曲線

分析圖2～圖4可知，無論是快縱波、慢縱波還是橫波，因空間網格離散，數值速度隨傳播方向不同大于或小于真實速度，即發(fā)生頻散；一般情況下，當00.2時，隨著單位波長內網格點數的減少而數值頻散現象越來越嚴重。

2.2 時間頻散分析

應用Matlab編程，計算式(4)～式(6)可得，時間差分網格離散引起的快縱波、慢縱波和橫波的數值頻散曲線，如圖5～圖8所示，其中橫坐標為時間步長Δt，縱坐標為數值頻散qp-fast、qp-slow和qs。在二維平面內，考慮到Ax、Az關于角度γx、γz取值的對稱性，考察3個不同的傳播方向：沿x軸(γx=0°，γz=90° )，沿xz平面三分線(γx=30°，γz=60° )，沿xz平面對角線(γx=45°，γz=45°)。相關參數取值如表2所示。

表2 時間頻散計算參數

其中，vp-solid、vshear-solid和的單位為m/s。

圖5 時間二階有限差分快P波數值頻散曲線

圖6 時間二階有限差分慢P波數值頻散曲線

圖7 時間二階有限差分橫波數值頻散曲線

分析圖5～圖7可知，無論是快縱波、慢縱波還是橫波，因時間離散，數值速度均大于真實速度，即發(fā)生頻散，但在不同傳播方向上頻散情況基本相同；一般情況下，當0<Δt≤0.2×10-3時，基本上不發(fā)生頻散，當Δt>0.2×10-3時，隨著時間步長的增加，數值頻散現象越來越嚴重。基此可知，空間步長和時間步長越小，數值頻散現象越弱，但相應地計算量會增加，計算效率降低，因此我們在選取空間步長和時間步長的值時，應在減弱數值頻散的同時盡可能地減少計算量。

3 FCT校正分析及探討

在如圖8所示的Biot雙相各向同性介質模型中(邊界長度單位為m)，基于FCT校正公式詳見文獻[7]，應用Matlab編程，對比分析有、無使用FCT法壓制數值頻散的效果圖，如圖9～圖12所示。

圖8 FCT壓制效果計算模型

其中，η1是漫射因子，它是常量或一個線性函數，其取值取決于有限差分階段頻散誤差的大小。在實際應用中，漫射因子的值可以通過簡單地質模型的數值模擬運算來確定，當漫射因子取值合理時，數值模擬的結果不隨漫射因子的變化而明顯變化。一般情況下，0.01≤η1≤0.05比較合適，模擬效果較好[8]。η2是反漫射因子，η2的取值與η1不同，這是因為振幅和分辨率的損失主要有兩方面的原因，一是傳統(tǒng)有限差分運算引起的頻散，二是人為加入的漫射。反漫射運算不僅要補償人為加入漫射引起的損失，也要補償傳統(tǒng)有限差分運算帶來的振幅損失。因此，η2的取值應比η1大10%～15%[7]。

應用FCT校正方法，取η1=0.001,η2=0.0011，得到的頻散壓制效果圖如圖9所示。

(a)固相80ms波場快照

(b)固相檢波點波形圖

(c)流相80ms波場快照

(d)流相檢波點波形圖圖9 頻散壓制效果圖一

取η1=0.01，η2=0.011，得到的頻散壓制效果圖如10所示。

(a)固相80ms波場快照

(b)固相檢波點波形圖

(c)流相80ms波場快照

(d)流相檢波點波形圖圖10 頻散壓制效果圖二

取η1=0.1，η2=0.11，得到的頻散壓制效果圖如11所示。

(a)固相80ms波場快照

(b)固相檢波點波形圖

(c)流相80ms波場快照

(d)流相檢波點波形圖圖11 頻散壓制效果圖三

不使用FCT校正，得到的壓制效果圖如12所示。

分析圖9～圖12可知：當不使用FCT方法或η1、η2值相對設置較小時，波場模擬中數值頻散現象較嚴重；當η1=0.01，η2=0.011時，波場模擬中數值頻散得到很好壓制；當η1、η2值相對設置較大時，漫射通量和反漫射通量校正過度，波場模擬中的有效波受到干擾，從而使得波場快照和檢波點波形顯得模糊。以上說明，可以用FCT法消除數值頻散，但η1、η2值應控制在合理范圍，取值過小則失去壓制數值頻散的效果；取值過大則真實波場會被平滑變得模糊不清，分辨率顯著下降。本文建議η1取[0.01，0.05]，η2比η1大10%～15%。

(a)固相80ms波場快照

(b)固相檢波點波形圖

(c)流相80ms波場快照

(d)流相檢波點波形圖圖12 頻散壓制效果圖四

4 結論

在無耗散的二維Biot雙相各向同性介質中，對具有空間四階精度與時間二階精度的交錯網格差分格式的頻散問題，及FCT校正方法，經Matlab軟件編程計算分析，本文的結論有：

(1)空間和時間的網格離散都會帶來數值頻散。對于空間網格離散，一般情況下，當00.2時，隨著單位波長內網格點數的減少而數值頻散現象越來越嚴重。對于時間網格離散，一般情況下，當0<Δt≤0.2×10-2時，基本上不發(fā)生數值頻散；當Δt>0.2×10-3時，隨著時間步長的增加而數值頻散現象越來越嚴重，但在不同傳播方向上頻散情況基本相同。

(2)空間步長和時間步長越小，數值頻散現象越弱，但相應的計算量會增加，計算效率降低。為了進一步提高計算效率，減弱數值頻散，根據能量守恒定律引進通量校正傳輸法(FCT)，對波場漫射通量和反漫射通量進行校正。經大量試算發(fā)現，當FCT方法的系數η1、η2取值過小時，會失去壓制數值頻散的效果；當η1、η2取值過大時，會把真實波場平滑變得模糊不清，分辨率顯著下降；一般情況下，η1取[0.01，0.05]，η2比η1大10%至15%，壓制數值頻散效果較好。