彈性波方程正演混合吸收邊界的改進(jìn)*

崔麗蘋，張會(huì)星

(中國(guó)海洋大學(xué)海底科學(xué)與探測(cè)技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，海洋地球科學(xué)學(xué)院，山東 青島 266100)

波動(dòng)方程正演模擬是研究地震波在地下介質(zhì)中的傳播規(guī)律、波場(chǎng)耦合關(guān)系以及空間變化特征的有效工具[1-2]?；诓▌?dòng)方程模擬地震波在地下介質(zhì)中的傳播，由于計(jì)算空間有限而必須引入人工截?cái)噙吔?，人工邊界若不做特殊處理，其?huì)在邊界處產(chǎn)生反射導(dǎo)致在中心波場(chǎng)內(nèi)產(chǎn)生虛假波場(chǎng)，降低波場(chǎng)質(zhì)量。因此人工邊界通常要進(jìn)行特殊處理以消除人工邊界的反射，人工邊界的處理一直是地震波正演模擬中的重要研究課題。

波動(dòng)方程數(shù)值模擬中主要是應(yīng)用吸收邊界方法消除人工邊界反射，目前常用的吸收邊界條件有三種[3]：①單程波邊界條件，包括旁軸近似法[4]、多向透射邊界[5]等。這類方法計(jì)算簡(jiǎn)單、運(yùn)算時(shí)間短，但其效果依賴于波的入射角度。當(dāng)入射角較小時(shí)吸收效果較好，入射角較大時(shí)，吸收效果變差；②衰減型吸收邊界[6]。這類方法擴(kuò)充計(jì)算區(qū)域，在擴(kuò)充區(qū)內(nèi)設(shè)計(jì)阻尼因子對(duì)地震波的振幅進(jìn)行衰減，實(shí)現(xiàn)對(duì)地震波的吸收。但這種方法一方面阻尼因子比較難確定，吸收效果不佳，另一方面擴(kuò)充計(jì)算區(qū)域增加了計(jì)算量；③完全匹配層(PML)吸收邊界[7]。這種方法在計(jì)算區(qū)域外增加匹配層，在匹配層使用具有衰減效應(yīng)的波動(dòng)方程來(lái)吸收邊界反射。PML邊界吸收效果強(qiáng)于前兩種，是目前最常用的方法，同時(shí)計(jì)算成本也最高。

常規(guī)單程波吸收邊界條件通過(guò)求解雙程波方程得到內(nèi)部區(qū)域的波場(chǎng)，求解單程波方程得到邊界波場(chǎng)。由于方程的差異而在內(nèi)區(qū)和邊界區(qū)域產(chǎn)生波場(chǎng)突變，進(jìn)而產(chǎn)生邊界反射。為了解決這一問題，劉洋等[8]提出了混合吸收邊界條件，即在內(nèi)部區(qū)域和邊界之間加入一個(gè)過(guò)渡區(qū)，在過(guò)渡區(qū)內(nèi)對(duì)單程波和雙程波的波場(chǎng)進(jìn)行線性加權(quán)疊加，從而實(shí)現(xiàn)波場(chǎng)的平滑，并將其應(yīng)用到二階彈性波位移—應(yīng)力方程的正演模擬[9]。任志明等[3]在此基礎(chǔ)上提出了基于一階彈性波速度—應(yīng)力方程的混合吸收邊界條件實(shí)現(xiàn)方法。研究發(fā)現(xiàn)，在過(guò)渡區(qū)內(nèi)對(duì)單程波和雙程波波場(chǎng)線性加權(quán)疊加，其吸收效果與過(guò)渡帶的層數(shù)密切相關(guān)，當(dāng)過(guò)渡區(qū)層數(shù)較少時(shí)邊界反射不能被完全吸收，與其它吸收邊界條件相比混合吸收邊界的吸收效果沒有存在明顯優(yōu)勢(shì)[3,8]。

本文基于一階彈性波速度-應(yīng)力方程的混合吸收邊界條件存在的問題，從一階彈性波速度-應(yīng)力方程出發(fā)，對(duì)混合吸收邊界條件進(jìn)行改進(jìn)，提出一種過(guò)渡區(qū)指數(shù)型加權(quán)疊加方法。數(shù)值模擬結(jié)果表明：與常規(guī)過(guò)渡區(qū)線性加權(quán)疊加方法相比，本文方法吸收效果更好；同時(shí)與常規(guī)PML方法相比，本文方法在吸收效果、計(jì)算效率及存儲(chǔ)方面都有優(yōu)勢(shì)。

1 彈性波一階速度-應(yīng)力方程及高階交錯(cuò)網(wǎng)格差分格式

二維各向同性介質(zhì)中，彈性波一階速度-應(yīng)力方程[10]如式(1)所示：

(1)

其中:vx和vz分別代表質(zhì)點(diǎn)在x方向和z方向的速度分量；τxx、τzz為正應(yīng)力；τxz為剪切應(yīng)力；ρ為介質(zhì)密度；λ和μ是拉梅系數(shù)；t為時(shí)間；x和z代表空間位置。

在圖1所示的交錯(cuò)網(wǎng)格空間內(nèi)，根據(jù)一階彈性波方程交錯(cuò)網(wǎng)格有限差分算法[11-12]，得到時(shí)間2階、空間2N階的差分格式為：

圖1 交錯(cuò)網(wǎng)格示意圖Fig.1 Staggered-grid sketch map

(2)

(3)

當(dāng)交錯(cuò)網(wǎng)格差分格式的精度為時(shí)間2階空間2N階時(shí)，穩(wěn)定性條件[12]滿足:

(4)

2 混合吸收邊界條件

2.1 一階Higdon單程波方程及其差分格式

以左邊界的水平分量為例，Higdon給出的一階Higdon單程波差分格式[5]為：

(5)

(6)

2.2 混合吸收邊界條件計(jì)算方法

混合吸收邊界的計(jì)算思路為：求解內(nèi)部區(qū)域和過(guò)渡區(qū)域的雙程波方程，求解過(guò)渡區(qū)和邊界區(qū)的單程波方程，并對(duì)過(guò)渡區(qū)得到的單程波場(chǎng)和雙程波場(chǎng)進(jìn)行加權(quán)疊加，得到無(wú)邊界反射的波場(chǎng)。在圖2所示的計(jì)算空間內(nèi)，具體計(jì)算步驟如下所示：

圖2 混合吸收邊界示意圖Fig.2 Sketch map of hybrid absorbing boundary

(2) 在過(guò)渡區(qū)Ⅱ和邊界區(qū)域Ⅲ，根據(jù)公式(5)和(6)給出

(3) 在過(guò)渡區(qū)Ⅱ進(jìn)行單程波場(chǎng)和雙程波場(chǎng)的加權(quán)疊加，得到過(guò)渡區(qū)最終波場(chǎng)值。以左邊界為例，計(jì)算公式如(7)式所示：

(7)

常規(guī)線性加權(quán)疊加參數(shù)

wAi=(i-0.5)/N(i=1,2,LN)

wBi=(i-1)/N(i=2,3,LN)。

(8)

本文提出一種新的指數(shù)加權(quán)疊加參數(shù)如下所示：

(9)

圖3是疊加參數(shù)曲線的示意圖，其中左圖表示指數(shù)型疊加參數(shù)曲線，右圖以wBi的疊加參數(shù)為例，展示了線性疊加參數(shù)和指數(shù)型疊加參數(shù)曲線。指數(shù)型加權(quán)疊加參數(shù)的設(shè)置原則為：在內(nèi)部區(qū)域Ⅰ和過(guò)渡區(qū)Ⅱ的交界處，即i=N時(shí)使得雙程波的權(quán)重趨近于1，單程波的權(quán)重趨近于0；在過(guò)渡區(qū)Ⅱ和邊界區(qū)域Ⅲ的交界處，即i值最小時(shí)，雙程波權(quán)重和單程波權(quán)重都趨近于1/2。與線性加權(quán)疊加參數(shù)相比，指數(shù)型加權(quán)疊加參數(shù)在內(nèi)部區(qū)域Ⅰ和過(guò)渡區(qū)Ⅱ的交界處相對(duì)平緩，更能有效保留雙程波的信息，避免出現(xiàn)波場(chǎng)值突變出現(xiàn)的邊界反射。從圖上可以看出，在內(nèi)部區(qū)域Ⅰ和過(guò)渡區(qū)Ⅱ的交界處，即i=N時(shí)，相較線性疊加參數(shù)，指數(shù)疊加參數(shù)的變化更為平緩。

圖3 疊加參數(shù)曲線示意圖Fig.3 Sketch map of weighted parameters

3 數(shù)值模型試算

3.1 均勻介質(zhì)模型

采用均勻介質(zhì)模型進(jìn)行試驗(yàn)，模型大小水平方向?yàn)? 000 m，垂直方向?yàn)? 000 m，網(wǎng)格步長(zhǎng)為Δx=Δz=5 m，采樣間隔Δt=0.5 ms，震源采用主頻為30Hz的雷克子波，在模型中間激發(fā)，過(guò)渡區(qū)層數(shù)為10層，采用時(shí)間2階，空間12階差分格式對(duì)其進(jìn)行數(shù)值模擬，得到過(guò)渡區(qū)為使用線性加權(quán)參數(shù)和指數(shù)型加權(quán)參數(shù)的數(shù)值模擬結(jié)果。

比較圖4和圖5的波場(chǎng)快照，可以看出過(guò)渡區(qū)單程波和雙程波為線性疊加參數(shù)時(shí)，在邊界處仍有反射波不能被完全壓制。改進(jìn)指數(shù)型疊加參數(shù)能有效壓制邊界反射，避免邊界反射對(duì)波場(chǎng)的影響。

圖4 線性加權(quán)參數(shù)混合吸收邊界波場(chǎng)Fig.4 Wave field of linear weighted hybrid absorbing boundary

圖5 指數(shù)型加權(quán)參數(shù)混合吸收邊界波場(chǎng)Fig.5 Wave field of exponential weighted hybrid absorbing boundary condition

3.2 復(fù)雜模型試算

本文設(shè)置復(fù)雜模型進(jìn)行數(shù)值模擬來(lái)檢驗(yàn)算法的有效性和實(shí)用性。模型大小水平方向?yàn)? 000 m，垂直方向2 000 m，速度模型如圖所示。數(shù)值模擬采用網(wǎng)格步長(zhǎng)為Δx=Δz=5 m，采樣間隔Δt=0.5 ms，震源采用主頻為30 Hz的雷克子波，炮點(diǎn)位置為(1 000 m，5 m)，過(guò)渡區(qū)層數(shù)為30層，采用時(shí)間2階，空間12階差分格式對(duì)其進(jìn)行數(shù)值模擬，得到過(guò)渡區(qū)為使用線性加權(quán)參數(shù)和指數(shù)型加權(quán)參數(shù)地震記錄如圖所示。

對(duì)比圖7和圖8得到的過(guò)渡區(qū)分別為使用線性加權(quán)參數(shù)和指數(shù)型加權(quán)參數(shù)的數(shù)值模擬地震記錄，可以看出當(dāng)使用線性加權(quán)參數(shù)時(shí)，混合吸收邊界條件不能完全吸收掉因邊界產(chǎn)生的反射波，在地震記錄上可以看到邊界反射的存在(圖7箭頭指向處)，指數(shù)型加權(quán)參數(shù)得到的地震記錄則對(duì)邊界反射壓制效果較好，能夠有效壓制邊界偽反射。

圖6 復(fù)雜模型示意圖Fig.6 Sketch map of complex model

圖7 線性加權(quán)參數(shù)混合吸收邊界地震記錄Fig.7 Seismic record of linear weighted hybrid absorbing boundary

圖8 指數(shù)型加權(quán)參數(shù)混合吸收邊界地震記錄Fig.8 Seismic record of exponential weighted hybrid absorbing boundary condition

4 結(jié)語(yǔ)

本文從彈性波速度-應(yīng)力方程出發(fā)，對(duì)一階混合吸收邊界條件進(jìn)行改進(jìn)，提出了一種過(guò)渡區(qū)指數(shù)型加權(quán)疊加方法。模型數(shù)據(jù)檢驗(yàn)說(shuō)明，相對(duì)于常規(guī)過(guò)渡區(qū)線性加權(quán)疊加方法，本文方法對(duì)邊界反射吸收效果較好，能夠有效壓制虛假波場(chǎng)，改善波場(chǎng)質(zhì)量。