胡建林 宋維琪 張建坤 邢文軍 徐文會(huì)

(①中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島 266555； ②中國(guó)石油冀東油田公司勘探開(kāi)發(fā)研究院，河北唐山 063004)

1 引言

復(fù)雜地下介質(zhì)中，地震波的傳播過(guò)程繁冗，難以得到解析解，因此，一般是通過(guò)正演模擬探究地下地震波的傳播。在地震波正演模擬中，利用波動(dòng)方程的正演模擬比用運(yùn)動(dòng)學(xué)的射線追蹤法可獲得更豐富的動(dòng)力學(xué)信息，因此地震波場(chǎng)的數(shù)值模擬是地震波場(chǎng)傳播研究中的重要手段之一[1-8]。

對(duì)聲波方程進(jìn)行網(wǎng)格離散時(shí)，采用交錯(cuò)網(wǎng)格的差分格式，可較好地壓制數(shù)值頻散，得到更高精度的正演模擬結(jié)果[2]。交錯(cuò)網(wǎng)格是指把速度及壓力P分別定義于兩套不同網(wǎng)格上的網(wǎng)格系統(tǒng)，該方法能更好地處理一階波動(dòng)方程中速度和壓力的耦合關(guān)系。

在數(shù)值模擬時(shí)，邊界問(wèn)題就是找到一種合適的邊界條件，使人為產(chǎn)生的計(jì)算區(qū)域邊界在實(shí)際效果上是透明的，即不產(chǎn)生邊界反射[3]。吸收邊界條件正是一種產(chǎn)生反射非常小甚至不產(chǎn)生反射的人工邊界條件。Smith[9]將Dirichlet邊界條件與Neumann邊界條件結(jié)合起來(lái)，由于兩種邊界條件產(chǎn)生的反射符號(hào)相反，兩者相加抵消反射，但其計(jì)算量較大； Cerjan等[10]、Burns[11]提出基于擴(kuò)大計(jì)算區(qū)域的吸收邊界條件，該方法在擴(kuò)大后的邊界區(qū)使用阻尼機(jī)制減少反射波，以較大的計(jì)算量換取邊界的吸收； Kosloff等[12]提出黏性阻尼方法以減少邊界反射，但需進(jìn)行波數(shù)域的變換。Clayton等[13]在網(wǎng)格邊界附近用單向波動(dòng)方程模擬聲波和彈性波的傳播，該方法對(duì)于入射角小于45°的入射波很有效； Liao等[14]提出邊界外推方法，用來(lái)處理聲波和彈性波的邊界吸收問(wèn)題； 基于頻率域內(nèi)聲波和彈性波的吸收邊界條件，Randall[15]提出分別處理脹縮波和旋轉(zhuǎn)波分量的方法。Keys[16]和Higdon[17]在Clayton等[13]研究的基礎(chǔ)上提出一種改進(jìn)的吸收邊界條件，成功消除了任意角入射波的邊界反射，并具有良好穩(wěn)定性。Bécache等[18]研究了時(shí)間和空間上的m階Higdon邊界體條件，發(fā)現(xiàn)當(dāng)階數(shù)m較大時(shí)，Higdon邊界將變得異常復(fù)雜。

Berenger[19]在使用時(shí)域有限差分法時(shí)，為截?cái)喽S網(wǎng)格空間，提出PML吸收邊界概念； Chew等[20]將PML引入彈性波模擬，采用擴(kuò)展坐標(biāo)系的PML方法； Rodent等[21]進(jìn)一步提出CPML(convolutional complex frequency-shifted PML)邊界條件，該邊界在非線性、各向異性介質(zhì)時(shí)效果更好； Gedney等[22]提出ADE-PML(the PML with Auxiliary Differential Equations)邊界條件，該方法可將差分?jǐn)U展到更高時(shí)間階數(shù)； Liu等[23，24]提出一種混合吸收邊界，采用一種過(guò)渡帶方式從內(nèi)部區(qū)域逐漸過(guò)渡到邊界，進(jìn)一步削弱了邊界反射； Gao等[25]回顧了大部分典型的聲波方程邊界條件，并通過(guò)數(shù)值實(shí)驗(yàn)對(duì)比了各種方法。

至于將PML邊界與Higdon邊界的聯(lián)合使用，尚未找到相關(guān)文獻(xiàn)。本文將PML邊界與Higdon邊界構(gòu)成聯(lián)合吸收邊界同時(shí)使用，比單獨(dú)使用PML邊界減小了邊界厚度，也比單獨(dú)使用Higdon邊界具有更好的吸收效果。

2 三維聲波方程交錯(cuò)網(wǎng)格有限差分正演模擬

交錯(cuò)網(wǎng)格有限差分法與常規(guī)網(wǎng)格有限差分法相比，計(jì)算效率更高，且能進(jìn)一步提高數(shù)值模擬的運(yùn)算精度，有效壓制數(shù)值頻散。三維聲波方程交錯(cuò)網(wǎng)格有限差分正演模擬主要由傳播方程、震源、邊界條件三部分組成。

2.1 傳播方程

在交錯(cuò)網(wǎng)格有限差分正演中，使用一階壓力—速度方程。在非均勻各向同性介質(zhì)中，三維聲波方程中的一階壓力—速度方程為

(1)

(2)

式中：P為壓力；vx、vy、vz代表x、y、z方向的速度分量；V為縱波速度(聲波速度)；ρ為介質(zhì)密度。

將式(1)、式(2)有限差分離散，即可得到一階聲波方程的時(shí)域迭代計(jì)算。

2.2 震源

震源子波函數(shù)有雷克子波、高斯子波、正弦衰減子波等。本文選用雷克子波

f(t)=A[1-2(πf0t)2]exp[-(πf0t)2]

(3)

式中f0為子波主頻。

加入震源后，在震源處，式(1)變?yōu)?/p>

(4)

式中f(t)為震源子波函數(shù)，其他參數(shù)同式(1)。

2.3 邊界條件

在數(shù)值模擬中，只能在有限的區(qū)域范圍內(nèi)進(jìn)行模擬，如果不加任何處理，即在邊界處壓力值和速度值默認(rèn)設(shè)置為0，那么地震波到達(dá)人工邊界處時(shí)，將會(huì)產(chǎn)生非常強(qiáng)的反射，嚴(yán)重干擾內(nèi)部波場(chǎng)。為了減少人工邊界產(chǎn)生的反射波，需要再添加吸收邊界。

人工吸收邊界主要有兩種吸收方式：匹配吸收邊界條件(最外層單層匹配吸收)和衰減吸收邊界層(有一定厚度的多層衰減)。文中聯(lián)合吸收邊界中的Higdon邊界屬于前者，而PML吸收邊界則是后者。這兩種邊界的具體討論見(jiàn)下文。

除了上面的三個(gè)交錯(cuò)網(wǎng)格有限差分的重要組成部分，有限差分離散時(shí)還需對(duì)離散格式的穩(wěn)定性和頻散進(jìn)行分析。由于這部分并不是本文討論的重點(diǎn)，所以略過(guò)。

3 PML吸收邊界

PML吸收邊界通過(guò)在垂直邊界方向引入阻尼因子，選取適當(dāng)?shù)淖枘釁?shù)，使向邊界外傳播的波在匹配層中不斷衰減，從而使反射波最大程度上得到削弱，實(shí)現(xiàn)消除邊界反射的目的。該方法對(duì)一階壓力—速度方程時(shí)間域進(jìn)行分裂，得到含有吸收衰減因子的吸收邊界方程組

(5)

P=Px+Py+Pz

(6)

(7)

式中：dx、dy、dz分別為x、y、z方向的衰減因子，計(jì)算如下[26]

(8)

(9)

式中：L為PML邊界總厚度，lx、ly、lz為PML邊界內(nèi)的點(diǎn)到最外層點(diǎn)的距離；d0為吸收系數(shù)的初始值；R=1.0×10-6為理論反射系數(shù)[27]；β為控制PML吸收邊界的吸收效果參數(shù)，常常設(shè)置為匹配層內(nèi)的最大縱波速度。

PML吸收邊界為有一定厚度的多層邊界衰減，以沿-x軸方向傳播的vx和Px分量為例，吸收邊界的模式如圖1所示。

圖1 PML邊界層沿-x軸方向傳播的vx和Px分量吸收示意圖線的粗細(xì)代表能量大小

由圖1可見(jiàn)，要使PML吸收層達(dá)到很好的吸收效果，可將PML層加厚。這種簡(jiǎn)單增加層厚度就能實(shí)現(xiàn)好的吸收效果是PML邊界的優(yōu)點(diǎn)，同時(shí)也會(huì)出現(xiàn)計(jì)算量問(wèn)題。由于厚度的增加，使三維模型明顯胖了一圈，導(dǎo)致計(jì)算量更為龐大。

4 Higdon吸收邊界條件

Higdon[17，28]提出了具有下面形式的吸收邊界條件(以x軸方向上的邊界為例)

(10)

式中：m為吸收邊界條件的階數(shù)；c是入射波傳播速度；αj是吸收入射角，第j個(gè)算子能完全吸收入射角為±αj的波。式中的反射系數(shù)為

(11)

式中θ為入射角。

m=1時(shí)，為一階Higdon吸收邊界,式(11)變?yōu)?/p>

(12)

m=2時(shí)，為二階Higdon吸收邊界，就有

(13)

在差分網(wǎng)格上離散后，吸收邊界差分算子為

(14)

式中： 系數(shù)βj=cosαj，為正數(shù)； Δt、Δx分別為時(shí)間步長(zhǎng)和空間步長(zhǎng)；I為單位算子，滿足

(15)

Higdon吸收邊界為最外層邊界匹配吸收，以沿-x軸方向傳播的vx和Px分量為例，吸收邊界的模式如圖2。

Higdon邊界僅有單層厚度，不會(huì)明顯增加計(jì)算量，且能消除任意角度入射波的邊界反射，具有良好的穩(wěn)定性。但要想得到更好的吸收效果，必須使用高階的Higdon吸收邊界形式(m>2)，而高階時(shí)離散化后過(guò)于復(fù)雜，這是Higdon邊界的一大缺陷。

圖2 Higdon邊界沿-x軸方向傳播的vx和Px分量吸收示意圖

5 聯(lián)合吸收邊界

5.1 聯(lián)合吸收邊界方法

由上文可知，PML吸收效果是通過(guò)一定厚度邊界中的垂直分量衰減來(lái)吸收向邊界外傳播的波，而Higdon邊界是根據(jù)偏微分算子消除一定角度的入射波。根據(jù)PML吸收邊界和Higdon邊界的吸收原理，結(jié)合兩種吸收邊界的優(yōu)劣，研究了新的聯(lián)合吸收邊界，使向邊界傳播的波先經(jīng)過(guò)PML邊界進(jìn)行衰減，然后到達(dá)最外層Higdon邊界時(shí)，再對(duì)殘余的波進(jìn)行匹配吸收，最大限度地削弱邊界反射。以x軸方向的左邊垂直分量的邊界為例，吸收邊界的模式如圖3所示。

圖3 聯(lián)合吸收邊界上沿-x軸方向傳播的vx和Px分量吸收示意圖線的粗細(xì)代表能量大小

由于在PML層最外層聯(lián)合使用了Higdon吸收邊界，因此在相同吸收效果下，可減小PML層厚度，從而減少計(jì)算量。

5.2 聯(lián)合吸收邊界方程推導(dǎo)

針對(duì)研究的聯(lián)合吸收邊界方法，對(duì)新的吸收邊界方程進(jìn)行推導(dǎo)。為了最大程度簡(jiǎn)化推導(dǎo)，先考慮在PML吸收邊界內(nèi)沿x軸方向傳播的平面波，此時(shí)PML邊界內(nèi)的一階聲波方程(式(5)～式(7))簡(jiǎn)化為

(16)

設(shè)vx0、Px0分別為分量vx、Px的幅值，則PML介質(zhì)中傳播的波場(chǎng)可表示為[29]

(17)

將式(17)代入式(16)，可解得

(18)

將式(18)代入式(17)可得PML介質(zhì)中傳播的波場(chǎng)

(19)

將式(19)寫為更一般的平面波形式，可得PML邊界中沿著坐標(biāo)軸傳播的平面波方程，以速度V沿x軸向x>0空間運(yùn)動(dòng)的平面波表示為

(20)

由上式可推導(dǎo)得到x=0邊界上的Higdon形式的吸收邊界條件算子

(21)

將上面的算子應(yīng)用于沿x軸傳播的平面波時(shí)，反射波場(chǎng)的值該為零。因此在x=0邊界上，PML吸收層內(nèi)，m階的Higdon形式的吸收邊界條件可寫為

(22)

同理，可對(duì)波場(chǎng)區(qū)域的其他邊界上的PML層內(nèi)的Higdon形式吸收邊界進(jìn)行類似的推導(dǎo)，得到對(duì)應(yīng)的Higdon邊界條件。

5.3 聯(lián)合吸收邊界算法及其計(jì)算步驟

將推導(dǎo)得到的聯(lián)合吸收邊界方程，加入到有限差分離散方程中，通過(guò)差分算法的研究，形成了聯(lián)合吸收邊界的交錯(cuò)網(wǎng)格有限差分算法，算法的波場(chǎng)更新計(jì)算如下：

一、根據(jù)三維聲波方程中的一階壓力—速度方程(式(5))對(duì)邊界內(nèi)部與x，y，z有關(guān)的壓力分量單獨(dú)進(jìn)行波場(chǎng)更新；

二、對(duì)PML邊界內(nèi)的壓力分量采用含衰減因子的PML吸收邊界方程組(式(5)、式(8))進(jìn)行邊界內(nèi)壓力分量波場(chǎng)衰減更新；

三、對(duì)最外層Higdon邊界進(jìn)行壓力分量波場(chǎng)更新；

四、使用P=Px+Py+Pz(式(6))對(duì)各方向的壓力分量進(jìn)行求和；

五、根據(jù)三維聲波方程中的一階壓力—速度方程(式(7))對(duì)邊界內(nèi)部與x，y，z有關(guān)的速度分量單獨(dú)進(jìn)行波場(chǎng)更新；

六、對(duì)PML邊界內(nèi)的壓力分量采用含衰減因子的PML吸收邊界方程組(式(7)、式(8))進(jìn)行邊界內(nèi)速度分量波場(chǎng)衰減更新；

七、對(duì)最外層的Higdon邊界采用Higdon邊界吸收算子(式(22))進(jìn)行速度分量波場(chǎng)更新；

根據(jù)上述的波場(chǎng)更新步驟，反復(fù)進(jìn)行迭代，即可得到邊界吸收良好的聯(lián)合吸收邊界正演模擬。

6 數(shù)值分析

設(shè)計(jì)均勻各向同性的層狀介質(zhì)模型，介質(zhì)層數(shù)為7層，模型尺寸為500m×600m×700m，模型網(wǎng)格尺寸為5m×5m×5m，網(wǎng)格數(shù)為Nx×Ny×Nz=101×121×141。震源使用40Hz的雷克子波，設(shè)置位置為(250m，300m，0m)，接收線過(guò)震源，且與y方向平行，同樣間隔5m，共101個(gè)檢波器。層狀模型和炮檢分布如圖4，各層參數(shù)分布如表1所示。

使用高階交錯(cuò)網(wǎng)格有限差分算法，對(duì)層狀介質(zhì)進(jìn)行正演計(jì)算，當(dāng)PML厚度足夠時(shí)，能夠較好地吸收邊界的反射。當(dāng)PML邊界厚度為15個(gè)點(diǎn)，觀察270ms時(shí)壓力P的波場(chǎng)快照和480ms長(zhǎng)度的炮集記錄，可見(jiàn)有明顯的邊界反射(圖5b、圖6b)；在PML邊界外，再用一階Higdon邊界進(jìn)行吸收，可見(jiàn)邊界反射得到明顯減弱(圖5c、圖6c)；在PML邊界外，再用二階Higdon邊界進(jìn)行吸收，可見(jiàn)邊界反射得到進(jìn)一步減弱(圖5d、圖6d)，此時(shí)吸收效果與30個(gè)點(diǎn)厚度的PML邊界(圖5a、圖6a)相近。通過(guò)結(jié)合Higdon邊界，彌補(bǔ)因邊界厚度不足導(dǎo)致的吸收不足，體現(xiàn)出聯(lián)合吸收邊界的優(yōu)越性。

從表2正演模擬耗時(shí)統(tǒng)計(jì)可發(fā)現(xiàn)，與單純的PML邊界相比，達(dá)到相同的吸收效果時(shí)，聯(lián)合吸收邊界大大減小了運(yùn)算量，PML邊界(15點(diǎn)厚度)聯(lián)合二階Higdon邊界與PML邊界(30點(diǎn)厚度)相比，同樣吸收干凈，但耗時(shí)約為后者的62%，大大減少了計(jì)算量。

圖4 層狀模型和炮檢分布示意圖

層號(hào)層網(wǎng)格深度m速度m/s密度g/cm310～15025002.252150～25030002.263250～35035002.284350～45030002.275450～55036002.296550～65038002.307650～70036002.30

圖5 270ms時(shí)接收線剖面上P分量波場(chǎng)快照

(a)PML邊界(30個(gè)點(diǎn)邊界厚度)； (b)PML邊界(15個(gè)點(diǎn)邊界厚度)； (c)PML邊界(15個(gè)點(diǎn)邊界厚度)+一階Higdon邊界； (d)PML邊界(15個(gè)點(diǎn)邊界厚度)+二階Higdon邊界

圖6 P分量炮集記錄

(a)PML邊界(30個(gè)點(diǎn)邊界厚度)； (b)PML邊界(15個(gè)點(diǎn)邊界厚度)； (c)PML邊界(15個(gè)點(diǎn)邊界厚度)+一階Higdon邊界； (d)PML邊界(15個(gè)點(diǎn)邊界厚度)+二階Higdon邊界

7 結(jié)論

本文將PML邊界與Higdon邊界結(jié)合起來(lái)，得到聯(lián)合吸收邊界，對(duì)各向同性多層介質(zhì)的聲波波場(chǎng)模擬，并且通過(guò)對(duì)比分析，可發(fā)現(xiàn)聯(lián)合吸收邊界具有如下優(yōu)勢(shì)：

(1)聯(lián)合吸收邊界相比單純的PML吸收邊界，可使用較小的PML邊界厚度，使模型大小明顯減小，較明顯地減少了計(jì)算量；

(2)與Higdon邊界相比，聯(lián)合吸收邊界結(jié)合了PML的吸收特性，能夠更有效地吸收邊界反射，獲得更好的模擬結(jié)果。