交錯網格有限差分正演模擬的聯合吸收邊界

胡建林 宋維琪 張建坤 邢文軍 徐文會

(①中國石油大學(華東)地球科學與技術學院，山東青島 266555； ②中國石油冀東油田公司勘探開發研究院，河北唐山 063004)

1 引言

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

對聲波方程進行網格離散時，采用交錯網格的差分格式，可較好地壓制數值頻散，得到更高精度的正演模擬結果[2]。交錯網格是指把速度及壓力P分別定義于兩套不同網格上的網格系統，該方法能更好地處理一階波動方程中速度和壓力的耦合關系。

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

Berenger[19]在使用時域有限差分法時，為截斷二維網格空間，提出PML吸收邊界概念； Chew等[20]將PML引入彈性波模擬，采用擴展坐標系的PML方法； Rodent等[21]進一步提出CPML(convolutional complex frequency-shifted PML)邊界條件，該邊界在非線性、各向異性介質時效果更好； Gedney等[22]提出ADE-PML(the PML with Auxiliary Differential Equations)邊界條件，該方法可將差分擴展到更高時間階數； Liu等[23，24]提出一種混合吸收邊界，采用一種過渡帶方式從內部區域逐漸過渡到邊界，進一步削弱了邊界反射； Gao等[25]回顧了大部分典型的聲波方程邊界條件，并通過數值實驗對比了各種方法。

至于將PML邊界與Higdon邊界的聯合使用，尚未找到相關文獻。本文將PML邊界與Higdon邊界構成聯合吸收邊界同時使用，比單獨使用PML邊界減小了邊界厚度，也比單獨使用Higdon邊界具有更好的吸收效果。

2 三維聲波方程交錯網格有限差分正演模擬

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

2.1 傳播方程

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

(1)

(2)

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

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

2.2 震源

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

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

(3)

式中f0為子波主頻。

加入震源后，在震源處，式(1)變為

(4)

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

2.3 邊界條件

在數值模擬中，只能在有限的區域范圍內進行模擬，如果不加任何處理，即在邊界處壓力值和速度值默認設置為0，那么地震波到達人工邊界處時，將會產生非常強的反射，嚴重干擾內部波場。為了減少人工邊界產生的反射波，需要再添加吸收邊界。

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

除了上面的三個交錯網格有限差分的重要組成部分，有限差分離散時還需對離散格式的穩定性和頻散進行分析。由于這部分并不是本文討論的重點，所以略過。

3 PML吸收邊界

PML吸收邊界通過在垂直邊界方向引入阻尼因子，選取適當的阻尼參數，使向邊界外傳播的波在匹配層中不斷衰減，從而使反射波最大程度上得到削弱，實現消除邊界反射的目的。該方法對一階壓力—速度方程時間域進行分裂，得到含有吸收衰減因子的吸收邊界方程組

(5)

P=Px+Py+Pz

(6)

(7)

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

(8)

(9)

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

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

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

由圖1可見，要使PML吸收層達到很好的吸收效果，可將PML層加厚。這種簡單增加層厚度就能實現好的吸收效果是PML邊界的優點，同時也會出現計算量問題。由于厚度的增加，使三維模型明顯胖了一圈，導致計算量更為龐大。

4 Higdon吸收邊界條件

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

(10)

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

(11)

式中θ為入射角。

m=1時，為一階Higdon吸收邊界,式(11)變為

(12)

m=2時，為二階Higdon吸收邊界，就有

(13)

在差分網格上離散后，吸收邊界差分算子為

(14)

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

(15)

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

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

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

5 聯合吸收邊界

5.1 聯合吸收邊界方法

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

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

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

5.2 聯合吸收邊界方程推導

針對研究的聯合吸收邊界方法，對新的吸收邊界方程進行推導。為了最大程度簡化推導，先考慮在PML吸收邊界內沿x軸方向傳播的平面波，此時PML邊界內的一階聲波方程(式(5)～式(7))簡化為

(16)

設vx0、Px0分別為分量vx、Px的幅值，則PML介質中傳播的波場可表示為[29]

(17)

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

(18)

將式(18)代入式(17)可得PML介質中傳播的波場

(19)

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

(20)

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

(21)

將上面的算子應用于沿x軸傳播的平面波時，反射波場的值該為零。因此在x=0邊界上，PML吸收層內，m階的Higdon形式的吸收邊界條件可寫為

(22)

同理，可對波場區域的其他邊界上的PML層內的Higdon形式吸收邊界進行類似的推導，得到對應的Higdon邊界條件。

5.3 聯合吸收邊界算法及其計算步驟

將推導得到的聯合吸收邊界方程，加入到有限差分離散方程中，通過差分算法的研究，形成了聯合吸收邊界的交錯網格有限差分算法，算法的波場更新計算如下：

一、根據三維聲波方程中的一階壓力—速度方程(式(5))對邊界內部與x，y，z有關的壓力分量單獨進行波場更新；

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

三、對最外層Higdon邊界進行壓力分量波場更新；

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

五、根據三維聲波方程中的一階壓力—速度方程(式(7))對邊界內部與x，y，z有關的速度分量單獨進行波場更新；

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

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

根據上述的波場更新步驟，反復進行迭代，即可得到邊界吸收良好的聯合吸收邊界正演模擬。

6 數值分析

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

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

從表2正演模擬耗時統計可發現，與單純的PML邊界相比，達到相同的吸收效果時，聯合吸收邊界大大減小了運算量，PML邊界(15點厚度)聯合二階Higdon邊界與PML邊界(30點厚度)相比，同樣吸收干凈，但耗時約為后者的62%，大大減少了計算量。

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

層號層網格深度m速度m/s密度g/cm310～15025002.252150～25030002.263250～35035002.284350～45030002.275450～55036002.296550～65038002.307650～70036002.30

圖5 270ms時接收線剖面上P分量波場快照

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

圖6 P分量炮集記錄

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

7 結論

本文將PML邊界與Higdon邊界結合起來，得到聯合吸收邊界，對各向同性多層介質的聲波波場模擬，并且通過對比分析，可發現聯合吸收邊界具有如下優勢：

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

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