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

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

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

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

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

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

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

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

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

1.1 透明邊界條件

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

1.2 Clayton-Engquist邊界條件

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

1.3 PML邊界條件

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

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

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

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

因?yàn)樵诙S聲場中任意方向傳播的波都可以分解為沿x方向和沿z方向傳播的兩個(gè)波，所以，在匹配層中分別對這兩個(gè)波進(jìn)行衰減吸收，便可以達(dá)到對任意方向傳播的波進(jìn)行衰減的目的。PML邊界條件的策略就是在需要的邊界區(qū)域加入相應(yīng)的PML層，達(dá)到衰減模型人工邊界反射的目的。如圖1所示，中間部分為模型區(qū)域，而四周則是用來吸收模型邊界反射而鑲加上去的PML層。其中各個(gè)邊界區(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 模擬計(jì)算與結(jié)果

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表1 各記錄點(diǎn)反射波振幅的最大峰值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隨觀測點(diǎn)與模型邊界間距離的變化Fig.8 Variation of the peak ratios of CE/TBC and PML/TBC with the distance between the observation point and model boundary

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

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

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

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

進(jìn)行頻率試驗(yàn)的模型與前述試驗(yàn)的模型相同。在選取觀測點(diǎn)時(shí)，考慮到峰值比在D點(diǎn)的位置上已經(jīng)相對比較穩(wěn)定，所以選擇D點(diǎn)為源點(diǎn)/觀測點(diǎn)進(jìn)行頻率模擬實(shí)驗(yàn)。在5～150Hz的頻段范圍內(nèi)，以5Hz的間隔，做了31次頻率實(shí)驗(yàn)。用與3.1節(jié)的實(shí)驗(yàn)數(shù)據(jù)類似的處理方法，計(jì)算出各個(gè)頻率在D點(diǎn)的反射波峰值的CE/TBC與PML/TBC比值，繪成曲線（圖9）。從中可以看到CE/TBC在低頻時(shí)變化較大。但是當(dāng)頻率＞50Hz時(shí)，CE/TBC相對較穩(wěn)定，均在0.3左右；也就是說，在頻率＞50Hz的情況下，CE邊界條件的吸收效果是TBC的3倍左右；而PML/TBC呈線性上升關(guān)系，表明PML相對于TBC的吸收優(yōu)勢隨著頻率上升是逐漸變?nèi)醯?。但是，在我們?shí)驗(yàn)的150Hz有效范圍之內(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é)論是在模型邊界條件的垂直反射情況下得到的。為了了解不同反射角度情況下各邊界條件的吸收效果，做了一項(xiàng)模擬實(shí)驗(yàn)。由于需要覆蓋較寬的反射角度，對前面用過的模型進(jìn)行了擴(kuò)展，此模擬實(shí)驗(yàn)采用的模型大小為2km×2 km，模型介質(zhì)波速和模型離散間距不變。

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

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

圖10 反射角模擬實(shí)驗(yàn)觀測系統(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

在實(shí)驗(yàn)的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)勢顯著增強(qiáng)，在54°左右的時(shí)候，其吸收效果是TBC的27（＝1/0.036 7）倍左右。但是隨著反射角增大，其吸收效果優(yōu)勢開始顯著地變?nèi)酢Ｔ?5°反射情況下，CE的吸收效果相對于TBC的優(yōu)勢已經(jīng)很小了，此時(shí)的CE/TBC比值為0.67左右。也就是說，在75°反射情況下，CE的吸收效果是TBC的1.5倍左右，優(yōu)勢已經(jīng)大幅減低。

4 分析與結(jié)論

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

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

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

c.在實(shí)驗(yàn)的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)勢顯著增強(qiáng)，在54°左右的時(shí)候，其吸收效果是TBC的27倍左右。但是隨著反射角繼續(xù)增大，其吸收效果優(yōu)勢開始顯著地變?nèi)?。?5°反射情況下，CE的吸收效果相對于TBC的優(yōu)勢已經(jīng)很小了，此時(shí)的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］.西安石油大學(xué)學(xué)報(bào)：自然科學(xué)版，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è)大學(xué)學(xué)報(bào)，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ì)算機(jī)仿真初步研究［J］.計(jì)算機(jī)科學(xué)，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.