周成峰

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

0 引言

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

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

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

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

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

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

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

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

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

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

1 頻散計算及FCT校正

1.1 空間與時間頻散計算

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

(1)

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

物理量q可以反映數(shù)值速度與真實速度的關(guān)系。

為了方便公式推導(dǎo)和分析，引入無量綱參數(shù)ζ和H：

(2)

式中λ為波長；

Δt為時間步長；

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

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

(3)

其中：

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

則γx、γz為平面波傳播方向與坐標(biāo)軸x、z所成夾角，c1、c2為空間四階差分權(quán)系數(shù)。

當(dāng)x和z方向的空間步長Δx=Δz=Δ時，引入控制數(shù)值頻散的無量綱參數(shù)ζ=VP-solidΔt/Δ和網(wǎng)格尺寸無量綱參數(shù)H=Δ/λ，結(jié)合不考慮耗散的二維Biot各向同性介質(zhì)中三類波的速度表達(dá)式及式(3)可推得規(guī)則網(wǎng)格快縱波數(shù)值頻散qp-fast,慢縱波數(shù)值頻散qp-slow和橫波數(shù)值頻散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分別為固體基質(zhì)(固相)密度和孔隙流體(流相)密度；

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

φ為孔隙度；

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

γ是與流體運(yùn)動中骨架的微觀模型有關(guān)的因子，當(dāng)流體中的固體顆粒為球形時，γ=1/2。

1.2 FCT校正方法

FCT方法基本原理是假設(shè)所有的極值點都是由數(shù)值頻散引起的，對所有網(wǎng)格點進(jìn)行漫射校正處理，達(dá)到消除數(shù)值頻散的效果。再對非局部極值點進(jìn)行補(bǔ)償?shù)姆绰湫Ｕ谷魏尾恍枰涞牡胤降靡曰謴?fù)。實際上，漫射校正是一個線性的平滑過程，反漫射校正由于是有選擇的，所以是一個非線性的過程[6]。對于交錯網(wǎng)格有限差分速度-應(yīng)力波動方程，應(yīng)用FCT方法對差分?jǐn)?shù)值解進(jìn)行校正的過程中，可以只對速度進(jìn)行計算，也可以只對應(yīng)力進(jìn)行校正，因為在完成對速度或應(yīng)力的校正后，通過交錯網(wǎng)格有限差分速度-應(yīng)力波動方程對應(yīng)力或速度進(jìn)行了校正。顯然，兩種方法都需對4個參數(shù)進(jìn)行校正，運(yùn)算量是相同的，本文選擇對速度進(jìn)行校正。

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

2 頻散分析

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

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

2.1 空間頻散分析

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

表1 空間頻散計算參數(shù)

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

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

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

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

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

2.2 時間頻散分析

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

表2 時間頻散計算參數(shù)

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

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

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

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

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

3 FCT校正分析及探討

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

圖8 FCT壓制效果計算模型

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

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

(a)固相80ms波場快照

(b)固相檢波點波形圖

(c)流相80ms波場快照

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

4 結(jié)論

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

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

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