波動方程正演模擬邊界條件的比較分析

付小波，韓 超，原健龍，余嘉順，2

（1.成都理工大學 地球物理學院，成都610059；2.新西蘭皇家地質(zhì)與核科學研究所，惠靈頓）

邊界條件是地震波數(shù)值模擬方法技術中的一項重要內(nèi)容。許多專家學者從不同角度提出多種構造邊界條件的方法。1977年，Clayton與Engquist［1］根據(jù)旁軸近似理論（Claerbout［2］；Claerbout與Johnson［3］），提出利用一系列不同近似精度的單程波動方程來吸收模型邊界的反射能量，稱作Clayton-Engquist（下文采用簡略記號CE來表達）邊界條件。1978年，Reynolds通過對波動方程的分解得到了透明邊界條件［4］（下文采用簡略記號TBC來表達）。1994年，Berenger針對電磁波傳播模擬，提出了一種高效的吸收邊界條件［5］，即完全匹配層（Perfectly Matched Layer，PML）。Berenger結(jié)合實際理論模型證明了該方法可以很好地吸收來自各個方向、各種頻率的電磁波。1995年，F(xiàn)rank，Hastings等將PML的算法思想引入彈性波模擬中，建立了彈性波場模擬的PML吸收邊界條件［6］。

此后，人們針對這些邊界條件的質(zhì)量與效率做了一些比較性研究。2006年，李景葉、陳小宏等分別采用CE邊界條件和PML邊界條件做一階的速度－應力彈性波方程數(shù)值模擬，得到不同的模擬圖像［7］。經(jīng)過觀察比較，得出PML邊界的效果優(yōu)于CE邊界的定性認識。同年，邢麗也分別用CE邊界條件和PML邊界條件進行二維聲波方程數(shù)值模擬，以圖像顯示觀察比較的方式，說明了PML邊界條件明顯好于CE邊界條件［8］。

上述基于圖件觀察比較得出的定性結(jié)論對一般性的勘探成像模擬是具有指導意義的；但是對于弱能量波場現(xiàn)象的模擬來說，就顯得不夠精細了。比如在模擬粗糙界面的弱散射現(xiàn)象時，其散射能量與邊界條件殘存反射能量在強度上可能是同級的，如果邊界條件選擇不當，邊界反射造成的干擾甚至會將目標弱散射完全淹沒。在此情況下，必須對邊界條件的吸收效果有定量的精確把握，才能保證目標現(xiàn)象模擬的可行性。因此，對上述幾種邊界條件的吸收效果作高精度比較分析，得出精確的定量結(jié)論，就顯得至關重要。本文對TBC邊界條件、CE邊界條件和PML邊界條件的吸收效果進行了精細的定量模擬比較研究。

1 聲波模擬方程與邊界條件

為了對3種吸收邊界條件的數(shù)值計算效果做比較，我們考慮一個均勻的二維聲場介質(zhì)空間Ω＝｛（x，z）｜0≤x≤X，0≤z≤Z｝以及時間區(qū)間τ＝｛t｜0≤t≤T｝，在此空間上聲波的傳播速度為v。在此空間上的聲波傳播過程由如下聲波方程［9］決定

我們用δ＞0對Ω做空間等距分割，用Δ＞0對τ做時間等距分割，并對（1）式中的時間二階導數(shù)和空間二階導數(shù)分別采用二階中心差分和四階中心差分來離散。整理后得到如下關于波場傳播的離散遞推關系

上式中表示k時刻聲波場在坐標位置為（i，j）的那個質(zhì)點上的振幅。其中G為空間導數(shù)的四階差分

1.1 透明邊界條件

Reynolds根據(jù)標量聲波方程（1）得到的透明邊界條件為

1.2 Clayton-Engquist邊界條件

Clayton和Engquist提出的一系列不同精度的吸收邊界條件中，如下二階近似邊界條件［1］應用得最為普遍

1.3 PML邊界條件

二維聲波方程的PML控制方程［10］為

其中ux，uz，u′x，u′z是中間變量；d（x），d（z）分別是x方向和z方向的衰減系數(shù)，其函數(shù)表達式［11］為

式中：R為反射系數(shù)；L為邊界層厚度；w為計算網(wǎng)格點到內(nèi)邊界的距離。

圖1 PML邊界區(qū)域分布示意圖Fig.1 Schematic diagram for PML boundary

因為在二維聲場中任意方向傳播的波都可以分解為沿x方向和沿z方向傳播的兩個波，所以，在匹配層中分別對這兩個波進行衰減吸收，便可以達到對任意方向傳播的波進行衰減的目的。PML邊界條件的策略就是在需要的邊界區(qū)域加入相應的PML層，達到衰減模型人工邊界反射的目的。如圖1所示，中間部分為模型區(qū)域，而四周則是用來吸收模型邊界反射而鑲加上去的PML層。其中各個邊界區(qū)域具有如下性質(zhì)：在區(qū)域1中，d（x）≠0，d（z）≠0；在區(qū)域2中，d（x）＝0，d（z）≠0；在區(qū)域3中，d（x）≠0，d（z）＝0。

2 模擬計算與結(jié)果

我們考察一個點聲源所產(chǎn)生的波的邊界反射情況，模擬計算在一個600m×600m的區(qū)域內(nèi)進行。模型介質(zhì)的波傳播速度為1.5km/s，為理想水介質(zhì)的聲波速度。點聲源位于模型的中央。我們用1m的間隔對模型進行離散化，采用0.000 1s的時間間隔進行模擬計算。聲源信號子波函數(shù)為

其中：震源位置坐標xs＝300m，zs＝300m；信號子波主頻f＝100Hz；信號延遲因子t0＝9/（4f）；信號衰減因子

圖2為分別采用TBC邊界條件、CE邊界條件和PML邊界條件進行模擬得到的在0.25s時刻的波場快照，從中能夠看到3種吸收邊界條件對模型邊界反射都有明顯的吸收效果。通過對圖件的觀察可以得到初步定性認識：CE邊界條件與PML邊界條件的吸收效果明顯好于TBC邊界條件。

盡管從圖2中采用PML邊界進行模擬得到的波場快照圖幾乎看不見任何邊界反射的痕跡，但這并不表明邊界反射不存在。事實上，如果對波場快照的能量進行適當增益之后再繪圖，PML邊界條件的反射波就顯示出來了（圖3）。

那么，PML邊界反射波的能量與其他2種邊界條件的反射波的能量的定量差別有多大呢？我們將通過下面的比較模擬例子來尋找此問題的答案。

3 邊界條件吸收效果的模擬比較

為了更好地了解不同邊界條件的吸收效果差異，我們對模型邊界的垂直反射的吸收效果，不同頻率成分的吸收效果、以及不同反射角的吸收效果進行比較分析。

圖2 3種不同邊界條件情況下在0.25s時刻的波場快照Fig.2 Snapshots of the wave field at 0.25sunder three boundary conditions

圖3 PML邊界條件波場快照提高增益后的顯示結(jié)果Fig.3 Snapshot of the wave field with gaining process using PML under boundary conditions

3.1 垂直反射的吸收效果

為了獲得3種邊界條件的吸收效果，本文采用同一個模型，分別鑲上不同的吸收邊界條件后進行模擬。所采用的模型大?。?.2km×1.2 km，模型介質(zhì)波速：1.5km/s。為了取得一個比較穩(wěn)定的結(jié)果，我們在每種邊界條件下均進行了5次模擬試驗，各次試驗的震源與接收點位置重疊，分別位于圖4中的A、B、C、D與E點的位置上。其位置坐標分別為（10，600）、（20，600）、（40，600）、（80，600）、（160，600）。所用的信號源子波函數(shù)由（18）式給定，子波主頻為100Hz。模擬過程是在空間離散間隔為δ＝1m的離散網(wǎng)格上用前述的有限差分方法進行的，模擬的時間采樣間隔Δ＝0.000 1s，模擬時長為0.3s。

我們的目的是對模型左邊界的反射波能量進行分析。考慮到試驗記錄中不但含有邊界反射波，而且還含有直達波，為此需要將直達波從記錄上消除，用純凈的反射波來進行比較分析，可以通過對模擬記錄與僅包含直達波的記錄相減來完成。

圖4 垂直反射模擬實驗觀測系統(tǒng)示意圖Fig.4 Simulation observation system of vertical reflection

為了得到一個僅包含直達波的記錄，將模型擴大為10km×10km來進行一次模擬計算。在此情況下，由于模型邊界距離源點和接收點非常遙遠，模型邊界反射尚未來得及傳回到接收點的位置上，模擬過程（0.3s）就已經(jīng)結(jié)束了。因此，這一模擬過程得到的記錄只含有直達波，沒有反射波。

圖5和圖6具體地展示了上述方法。圖5的Ⅰ、Ⅱ、Ⅲ、Ⅳ分別為在擴展模型鑲透明邊界、CE邊界以及PML邊界情況下在D點激發(fā)和接收的記錄。圖6的3條時間記錄曲線（Ⅱ’、Ⅲ’、Ⅳ’）分別為圖5中Ⅱ、Ⅲ、Ⅳ減去Ⅰ的結(jié)果，為所需要的用于分析邊界條件吸收效果的純凈反射波。

讀取純凈反射波的波峰值，列于表1之中，可對3種邊界條件的吸收差異作對比分析。從表1的數(shù)據(jù)可以看出，采用不同的邊界條件，模型邊界反射的強弱存在顯著的差異。其中CE邊界條件下的反射波峰值（表中第3列）小于TBC邊界條件（表中第2列），而PML邊界條件下的反射波峰值（第4列）又小于CE邊界條件下的反射波峰值。這種情況在圖7中體現(xiàn)也非常明顯。而且這種相對差異在A、B、C、D與E震源/接收情況體現(xiàn)得非常一致。這在圖8的關于CE/TBC和PML/TBC的反射波峰值比值的圖示中體現(xiàn)得非常直觀。

圖5 經(jīng)過直達波峰值歸一化處理的全程波形記錄Fig.5 Full waveform records normalised with the peak amplitude of the direct wave

圖6 經(jīng)過直達波峰值歸一化處理的反射波波形記錄Fig.6 Reflection waveform records normalised with the peak amplitude of the direct wave

表1 各記錄點反射波振幅的最大峰值Table 1 Peak values of the reflected wave amplitudes recorded at each point

圖7 3種邊界條件反射波峰值的比較Fig.7 Comparison of the reflection peak values under three boundary conditions

圖8 振幅峰值比CE/TBC及PML/TBC隨觀測點與模型邊界間距離的變化Fig.8 Variation of the peak ratios of CE/TBC and PML/TBC with the distance between the observation point and model boundary

當源點/接收點位置不同，上述峰值的比值存在差異。幾何擴散會對反射波峰隨著位置產(chǎn)生變化；但是，CE/TBC與PML/TBC比值已經(jīng)消除了幾何擴散的影響，因此CE/TBC與PML/TBC隨著位置的變化必然另有原因。我們認為，這種相對位置造成的比值差異，可能是由于不同邊界條件的頻率響應差異引起的。當一個邊界條件的響應頻段較低的時候，由于對應的主頻較長，就可對距離模型邊界較遠的模型內(nèi)部波場產(chǎn)生影響。反之，如果一個邊界條件的響應頻段較高，其對距離模型邊界較遠的模型內(nèi)部波場產(chǎn)生的影響就較弱。因此，CE/TBC與PML/TBC的值隨著位置的變化而變化，反映了3種邊界條件頻率響應的差異。

通過比較試驗得到如下認識：（1）在垂直于模型邊界反射的情況下，CE邊界條件的吸收效果優(yōu)于TBC的吸收效果。在距離邊界條件較遠的區(qū)域，兩者的相對比值有逐漸進入穩(wěn)定狀態(tài)的趨勢（圖8），據(jù)此，其穩(wěn)定狀態(tài)比值應該不大于0.287 5（表1）。也就是說，在進入穩(wěn)定狀態(tài)之后，CE邊界條件的吸收效果至少比TBC邊界條件好3.5（＝1/0.287 5）倍。（2）PML邊界條件的吸收效果又進一步顯著地優(yōu)于CE邊界條件的吸收效果；并且在兩者逐漸進入相對穩(wěn)定的狀態(tài)下（圖8），PML邊界條件的吸收效果至少比TBC邊界條件好16.5（＝1/0.060 4）倍。

3.2 不同頻率成分的吸收效果

前面關于PML和CE邊界條件優(yōu)于TBC邊界條件的結(jié)論是在源信號主頻為100Hz的情況下得到的。這一結(jié)論在其他頻段是否依然成立是一個需要研究的問題。為此，對于3種邊界條件對不同頻率波場成分的吸收效果，采用不同頻率的源信號進行了一系列的模擬分析。

進行頻率試驗的模型與前述試驗的模型相同。在選取觀測點時，考慮到峰值比在D點的位置上已經(jīng)相對比較穩(wěn)定，所以選擇D點為源點/觀測點進行頻率模擬實驗。在5～150Hz的頻段范圍內(nèi)，以5Hz的間隔，做了31次頻率實驗。用與3.1節(jié)的實驗數(shù)據(jù)類似的處理方法，計算出各個頻率在D點的反射波峰值的CE/TBC與PML/TBC比值，繪成曲線（圖9）。從中可以看到CE/TBC在低頻時變化較大。但是當頻率＞50Hz時，CE/TBC相對較穩(wěn)定，均在0.3左右；也就是說，在頻率＞50Hz的情況下，CE邊界條件的吸收效果是TBC的3倍左右；而PML/TBC呈線性上升關系，表明PML相對于TBC的吸收優(yōu)勢隨著頻率上升是逐漸變?nèi)醯?。但是，在我們實驗?50Hz有效范圍之內(nèi)，PML的吸收優(yōu)勢依然十分顯著，即便是在150Hz的相對高頻上，其吸收效果依然比TBC好10倍（PML/TBC＜0.1）以上。

圖9 反射波峰值振幅比CE/TBC與PML/TBC隨頻率的變化曲線Fig.9 Peak reflected amplitude ratios CE/TBC and PML/TBC versus frequency change

3.3 不同反射角的吸收效果比較

上述結(jié)論是在模型邊界條件的垂直反射情況下得到的。為了了解不同反射角度情況下各邊界條件的吸收效果，做了一項模擬實驗。由于需要覆蓋較寬的反射角度，對前面用過的模型進行了擴展，此模擬實驗采用的模型大小為2km×2 km，模型介質(zhì)波速和模型離散間距不變。

此項模擬的觀測系統(tǒng)如圖10所示。源信號子波頻率為100Hz，源點位置為S（80m，1 000 m），與模型邊界的距離與原先的D點相同。在與模型左邊界平行，由S點向上布置了一個700 m長、間隔為1m的觀測接收排列。S點觀測的反射角度為0°，而排列末端最后一個觀測點接收到的反射波的反射角為

所以這個觀測系統(tǒng)排列覆蓋了［0°，77°］的反射角度范圍。用與3.1節(jié)的實驗數(shù)據(jù)類似的處理方法，計算出各個反射角情況下反射波峰值的CE/TBC與PML/TBC比值曲線（圖11）。

圖10 反射角模擬實驗觀測系統(tǒng)示意圖Fig.10 Schematic diagram showing the simulation experiment observation system at the different reflection angles

圖11 反射波峰值振幅比CE/TBC與PML/TBC隨反射角的變化曲線Fig.11 Peak reflection amplitude ratios CE/TBC and PML/TBC versus reflection angles

在實驗的75°反射角范圍之內(nèi)，PML相對于TBC的吸收效果優(yōu)勢比較穩(wěn)定，受反射角的影響不大。而CE邊界條件相對于TBC的吸收效果優(yōu)勢則隨著反射角的變化存在較大的差異。在30°以內(nèi)的小反射角范圍內(nèi)，其吸收優(yōu)勢穩(wěn)定地保持在3.5倍左右。隨著反射角變大到50多度，其吸收優(yōu)勢顯著增強，在54°左右的時候，其吸收效果是TBC的27（＝1/0.036 7）倍左右。但是隨著反射角增大，其吸收效果優(yōu)勢開始顯著地變?nèi)?。?5°反射情況下，CE的吸收效果相對于TBC的優(yōu)勢已經(jīng)很小了，此時的CE/TBC比值為0.67左右。也就是說，在75°反射情況下，CE的吸收效果是TBC的1.5倍左右，優(yōu)勢已經(jīng)大幅減低。

4 分析與結(jié)論

本文針對透明邊界條件、Clayton-Engquist邊界以及PML邊界條件的吸收效果差異進行了模擬比較分析，得出如下定量結(jié)論。

a.以對模型邊界垂直反射的吸收效果來衡量，CE邊界條件對100Hz的波場成分的吸收效果是TBC吸收效果的3.5倍；而PML吸收效果則是TBC吸收效果的16.5倍。在3種邊界條件中，PML顯著地優(yōu)于其他兩種邊界條件。

b.CE邊界條件相對于TBC邊界條件的吸收優(yōu)勢在低頻時不明顯，而且變化較大；但是在頻率達到50Hz后CE的吸收效果是TBC的3.5倍左右。而PML相對于TBC的吸收優(yōu)勢隨著頻率上升是逐漸變?nèi)醯?；但是，在實驗?50Hz有效范圍之內(nèi)，PML的吸收優(yōu)勢依然十分顯著，即便是在150Hz的相對高頻上，其吸收效果依然比TBC好10倍以上。

c.在實驗的75°反射角范圍之內(nèi)，PML相對于TBC的吸收效果優(yōu)勢比較穩(wěn)定，受反射角變化的影響不大。而CE邊界條件相對于TBC的吸收效果在30°以內(nèi)的小反射角范圍內(nèi)，其吸收優(yōu)勢穩(wěn)定地保持在3.5倍左右。隨著反射角變大到50多度，其吸收優(yōu)勢顯著增強，在54°左右的時候，其吸收效果是TBC的27倍左右。但是隨著反射角繼續(xù)增大，其吸收效果優(yōu)勢開始顯著地變?nèi)?。?5°反射情況下，CE的吸收效果相對于TBC的優(yōu)勢已經(jīng)很小了，此時的CE/TBC為0.67左右。也就是說，在75°反射情況下，CE的吸收效果是TBC的1.5倍左右，優(yōu)勢已經(jīng)大幅減低。如果二者的比值隨角度的變化發(fā)展趨勢繼續(xù)保持下去，CE對TBC的優(yōu)勢在＞75°的大角度情況下出現(xiàn)逆轉(zhuǎn)是完全可能的。

總的來說，在3種邊界條件中，PML顯著地優(yōu)于其他2種邊界條件。相對于其他2種邊界條件，PML的吸收效果在150Hz頻率范圍內(nèi)和75°反射角范圍內(nèi)始終保持穩(wěn)定的相對優(yōu)勢。

［1］Robert Clayton，Bjorn Engquist.Absorbing boundary conditions for acoustic and elastic wave equation［J］.Bulletin of the Seismological Society of America，1977，66（11）：1529－1540.

［2］Claerbout J F.Coarse grid calculations of waves in inhomogeneous media with application to delineation of complicated seismic structure［J］.Geophysics，1970，35：407－418.

［3］Claerbout J F，Johnson A G.Extrapolation of time dependent wave forms along their path of propagation［J］.Geophysics Journal，1971，26：285－295.

［4］Reynolds A C.Boundary conditions for the numerical solution of wave propagation problems［J］.Geophysics，1978，43（6）：1099－1110.

［5］Jean-Pierre Berenger.A perfectly matched layer for the absorption of electromagnetic waves［J］.Journal of Computational Physics，1994，114：185－200.

［6］Hastings F D，Schneider J B，Broschat S L.Application of the perfectly matched layer（PML）absorbing boundary condition to elastic wave propagation［J］.Journal Acoustical Society of America，1996，100（5）：3061－3068.

［7］李景葉，陳小宏.TI介質(zhì)地震波場數(shù)值模擬邊界條件處理［J］.西安石油大學學報：自然科學版，2006，21（4）：20－23.Li J Y，Chen X H.Treatment of the boundary conditions in the numerical simulation of the seismic wave field in transversely isotropic（TI）medium［J］.Journal of Xi’an Shiyou University（Natural Science Edition），2006，21（4）：20－23.（In Chinese）

［8］邢麗.地震波數(shù)值模擬中的吸收邊界條件［J］.上海第二工業(yè)大學學報，2006，23（4）：272－278.Xing L.Absorbing boundary conditions for the numerical simulation of acoustic waves［J］.Journal of Shanghai Second Polytechnic University，2006，23（4）：272－278.（In Chinese）

［9］田力，余嘉順，唐紅.封閉空間聲場的計算機仿真初步研究［J］.計算機科學，2003，9（31）：222－223.Tian L，Yu J S，Tang H.Preliminary study on computer simulation sound wave-field in an enclosed space［J］.Computer Science，2003，9（31）：222－223.（In Chinese）

［10］王守東.聲波方程完全匹配層吸收邊界條件［J］.石油地球物理勘探，2003，18（1）：31－34.Wang S D.Absorbing boundary condition for acoustic wave equation by perfectly matched layer［J］.Oil Geophysical Prospecting，2003，38（1）：31－34.（In Chinese）

［11］Francis Collino，Chrysoula Tsogka.Application of the PML absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media［J］.Geophysics，2001，66（1）：294－307.