趙建國，史瑞其，陳競一，趙維俊，王宏斌，潘建國

1.中國石油大學(xué)（北京）地球物理與信息工程學(xué)院/油氣資源與探測國家重點(diǎn)實(shí)驗(yàn)室/CNPC物探重點(diǎn)實(shí)驗(yàn)室，北京 102249

2.中國海洋石油研究總院，北京 100027

3.塔爾薩大學(xué)地球科學(xué)學(xué)院，美國 塔爾薩 74104

4.沈陽地質(zhì)礦產(chǎn)研究所，沈陽 110034

5.中石油勘探開發(fā)研究院西北分院，蘭州 730000

0 引言

在地震波數(shù)值模擬中，為了在有限計(jì)算區(qū)域模擬無限大區(qū)域中地震波的傳播，需要對計(jì)算區(qū)域截斷邊界給出吸收邊界條件以消除邊界處產(chǎn)生的虛假數(shù)值反射。目前，廣泛使用的主要有2種構(gòu)造吸收邊界條件的方式：一種是通過波動方程因式分解獲得單行波方程的旁軸近似邊界條件，如Clayton－Engquist吸收邊界條件［1］、廖氏吸收邊界條件［2］、Higdon吸收邊界條件［3］等；另一種是在邊界上引入吸收介質(zhì)，使得地震波在無反射地進(jìn)入吸收介質(zhì)后被衰減掉，如完全匹配層（pefectly matched layer，PML）［4］。PML作為一種高效的吸收邊界，在地震波正演中被廣泛使用。PML最初采用對波場變量進(jìn)行非物理場分裂的方式來實(shí)現(xiàn)，使得計(jì)算內(nèi)域與PML吸收邊界區(qū)域的數(shù)值計(jì)算格式完全不同，實(shí)現(xiàn)方式較為繁瑣。后來，復(fù)拉伸坐標(biāo)［5］概念的引入使得PML可以通過非分裂場方式實(shí)現(xiàn)。

在使用非分裂場方式實(shí)現(xiàn)PML的算法中，需要引入輔助變量并在時間域迭代中對其進(jìn)行更新。通過輔助變量更新方式的不同，可以劃分出一系列不同的PML算法，如通過遞歸積分［6］和遞歸卷積［7］算法。以上算法使得輔助變量的時間計(jì)算格式局限于時間二階精度。但是，在波動方程數(shù)值模擬中，為了獲得更高的模擬精度，常常需要使用高階時間遞推格式；因此，上述非分裂PML在這種情況下的應(yīng)用受到限制。針對該問題，作者介紹了一種通過輔助微分方程求解輔助變量的 ADE－PML（auxiliary differential equation PML）［8］算法，將其用于高階時間精度（四階Runge－Kutta）的聲波方程數(shù)值模擬中。通過數(shù)值模擬，考察ADE－PML在高階精度時間遞推格式聲波方程模擬中的應(yīng)用效果并進(jìn)行理論分析。

1 方法原理

在二維直角坐標(biāo)系下，聲波方程的一階壓力－速度形式在時間域如下：

其中：p為聲壓；體積模量K＝ρc2，ρ為質(zhì)量密度，c為聲波在介質(zhì)中的傳播速度；vx和vz分別為沿x和z方向的速度場分量。注意，這里vx和vz是指介質(zhì)中質(zhì)點(diǎn)振動速度，c為波在介質(zhì)中傳播的群速度。

本節(jié)將介紹一種使用輔助微分方程在時間域更新輔助變量的非分裂完全匹配層算法，這種方法可以適應(yīng)高階精度的時間遞推格式。本文以四階Runge－Kutta時間遞推格式為例說明ADE－PML在高階時間遞推格式中的應(yīng)用。為構(gòu)造PML，需要對空間偏導(dǎo)數(shù)做復(fù)拉伸坐標(biāo)變換。為簡單起見，首先以x方向?yàn)槔挥懻搙正半坐標(biāo)軸的形式，并假設(shè)PML區(qū)域從x＝0開始。然后將x坐標(biāo)替換為復(fù)拉伸坐標(biāo)：

其中：sx是復(fù)拉伸函數(shù)，決定著PML的吸收效果。則關(guān)于x方向的偏導(dǎo)可以改寫為

z方向同理。因此，只需將關(guān)于x，z方向的偏導(dǎo)?x，?z替換為復(fù)拉伸坐標(biāo)下的??x和??z后即可。由于PML需在頻率域?qū)С?，以式?－1）的頻率域方程為例：

其中：i為虛數(shù)單位；ω為角頻率。當(dāng)確定復(fù)拉伸函數(shù)sx，sz的具體形式后，即可將式（4）變換回時間域。傳統(tǒng)的復(fù)拉伸函數(shù)為

其中：dx≥0是使地震波振幅在PML內(nèi)呈指數(shù)級下降的衰減系數(shù)。由于式（5）只對行波有效，對平行于PML和計(jì)算內(nèi)域交界傳播的隱失波和極低頻波無效［9－10］，因此，只能增加 PML厚度或?qū)⒄鹪捶胖糜谶h(yuǎn)離PML的位置，用以增加內(nèi)存和計(jì)算時間。為解決這個問題，Kuzuoglu等［11］將復(fù)拉伸函數(shù)的概念擴(kuò)展，提出了復(fù)頻移拉伸函數(shù)：

其中：κx和αx為輔助衰減系數(shù)；κx≥1，作用是將傾斜入射波矯正為垂直入射波；αx≥0使復(fù)平面的極點(diǎn)從實(shí)軸移動到虛負(fù)半平面上，作用是消除切入射波引起的低頻波和隱失波。為避免數(shù)值反射，各衰減系數(shù)需要漸進(jìn)變化。當(dāng)取κx＝1和αx＝0時，式（6）即退化為傳統(tǒng)復(fù)拉伸函數(shù)（5）。筆者主要采用指數(shù)函數(shù)［12－14］來定義 PML 的衰減系數(shù)分布：

其中：L是PML厚度；pd，pα和pκ為常數(shù)，通常取pd＝2，pα＝1，pκ＝2；dmax，αmax和κmax的最優(yōu)值依賴于計(jì)算模型參數(shù)或所選數(shù)值模擬方法。根據(jù)文獻(xiàn)［12］，衰減系數(shù)dmax通過關(guān)于垂直入射平面波的目標(biāo)反射系數(shù)R計(jì)算：

其中：cmax為計(jì)算模型的最大聲波速度。根據(jù)文獻(xiàn)［13］和［14］，輔助衰減系數(shù)αmax和κmax由以下關(guān)系確定：

其中：cmin為計(jì)算模型的最小聲波速度；W0為當(dāng)前所用空間離散格式每個波長所需的最小采樣點(diǎn)數(shù)；Δh為最大網(wǎng)格間距；f0為震源中心頻率。

ADE－PML的主要思想為將1/sx分解為2項(xiàng)之和，并使得iω只存在于其中一項(xiàng)的分母中：

則式（4）轉(zhuǎn)化為

其中，為與空間偏導(dǎo)數(shù)?xvx相關(guān)的輔助變量，定義為

將式（12）改寫為

將式（11）代入式（5），由頻率域變換回時間域后，得到非分裂形式的PML方程（同理處理項(xiàng)）：

可見要得到時域PML方程，只需將空間偏導(dǎo)項(xiàng)?xA替換為?xA/κx＋ψx即可。將式（13）變換回時間域，得到關(guān)于輔助變量的微分方程：

同理可得關(guān)于的輔助微分方程。ADE－PML的特點(diǎn)在于通過直接在時間域離散微分方程式（15）來求解輔助變量，因此適用于任意時間精度的遞推格式的數(shù)值模擬。下面以四階Runge－Kutta（以下簡稱RK4）格式為例說明ADE－PML在高階時間精度聲波方程數(shù)值模擬中的應(yīng)用。

在時間域直接對式（15）進(jìn)行RK4時間離散，輔助變量滿足其中：Δt為離散時間步長；下標(biāo)i對應(yīng)每個RK4循環(huán)中第i步迭代；當(dāng)θ取0，1/2和1時，分別對應(yīng)求取輔助變量的顯式、半隱式和隱式格式。對式（16）化簡可得

其中：

下面給出求解一階速度－壓力聲波方程式（1）的RK4時間遞推算法。設(shè)w＝（vx，vz，p），則由n時刻的wn更新n＋1時刻的wn＋1的過程可表達(dá)為

其中：算子L（·）表示式（1）中等式的右端項(xiàng)；RK4系數(shù)α2＝1/2，α3＝1/2，α4＝1；β1＝1/6，β2＝2/6，β3＝2/6，β4＝1/6。

2 數(shù)值模擬

為驗(yàn)證ADE－PML的吸收效果，使用四階精度旋轉(zhuǎn)交錯網(wǎng)格［15－16］進(jìn)行空間離散。模型計(jì)算區(qū)域?yàn)橐? 200m×1 600m的長方形區(qū)域，由651×201個網(wǎng)格點(diǎn)構(gòu)成，有限差分元胞邊長為Δx＝Δz＝8 m。介質(zhì)密度ρ＝2 000kg/m3，聲波速度cp＝3 000 m/s。時間上采用RK4高階時間更新格式，這里取Δt＝1.2ms，模擬持續(xù)時長3.6s，共計(jì)3 000時間步。采用點(diǎn)聲壓震源激發(fā)，震源子波為一導(dǎo)數(shù)高斯函數(shù)，其中a＝（πf0）2，震源中心頻率f0＝15Hz，子波時延t0＝1.2/f0＝0.06s，激發(fā)點(diǎn)位于（176m，496m）處；設(shè)置2個檢波點(diǎn)，檢波點(diǎn)1和檢波點(diǎn)2分別位于（1 456m，496 m）和（176m，4 896m）處；計(jì)算區(qū)域四周由PML吸收邊界包圍，厚度為80m，占10個有限差分元胞。震源和檢波點(diǎn)2靠近左側(cè)PML，檢波點(diǎn)1靠近右側(cè)PML。

為檢驗(yàn)復(fù)頻移拉伸算子的效果，分別給出使用傳統(tǒng)拉伸算子（5）的ADE－PML（這里稱為非復(fù)頻移ADE－PML）和復(fù)頻移拉伸算子（6）的 ADE－PML條件下的模擬計(jì)算結(jié)果。

該算例中，由于震源靠近左側(cè)吸收邊界，形成了高角度入射PML的情況。從波場快照圖1a可見，非復(fù)頻移ADE－PML對高角度入射波衰減不完全，入射波在左側(cè)PML內(nèi)形成損耗波，這些波產(chǎn)生了低頻虛假反射，且虛假反射能量隨著入射角度和傳播距離的增大而增強(qiáng)：從1.68s的波場可見，在遠(yuǎn)偏移距低頻虛假反射明顯影響了求解區(qū)域的波場。而圖1b中復(fù)頻移ADE－PML邊界由于使用復(fù)頻移拉伸算子，對高角度入射波吸收很好，在左側(cè)PML內(nèi)沒有產(chǎn)生損耗波和低頻反射，在圖1b中1.68s波場中，即使在遠(yuǎn)偏移距也沒有可見的虛假反射。

在地震記錄圖2a中，隨著偏移距的增大，高角度入射產(chǎn)生的虛假反射振幅逐漸變強(qiáng)，在直達(dá)波同相軸之后出現(xiàn)了一些“拖尾”現(xiàn)象，這是由于波場快照圖1中在左側(cè)吸收邊界附近產(chǎn)生的低頻虛假反射所致，從圖1中1.68s的波場可見，在遠(yuǎn)偏移距這些低頻反射尤為強(qiáng)烈。而圖2中復(fù)頻移ADE－PML由于使用了復(fù)頻移拉伸算子，對高角度入射波吸收很好，沒有出現(xiàn)圖1中的低頻虛假反射；在圖2中1.68s波場中，即使在遠(yuǎn)偏移距也沒有產(chǎn)生可見的虛假反射，在地震記錄圖3b中，同非復(fù)頻移ADEPML相比（圖3a），復(fù)頻移ADE－PML條件下僅可見直達(dá)波同相軸而沒有產(chǎn)生圖3a中的“拖尾”現(xiàn)象。

從檢波點(diǎn)1地震記錄（圖3a）可見，在檢波點(diǎn)1處，非復(fù)頻移ADE－PML和復(fù)頻移ADE－PML的計(jì)算結(jié)果都與精確解（即解析解）吻合良好，說明在垂直或小角度入射情況下，使用非復(fù)頻移和復(fù)頻移拉伸算子都可以高效吸收外行波場能量。因此，復(fù)頻移拉伸算子對垂直入射波的吸收無明顯改善。

而從檢波點(diǎn)2的地震記錄（圖3b）中可見：非復(fù)頻移ADE－PML條件下的地震記錄上產(chǎn)生了明顯的虛假反射抖動，與參考解吻合不佳，這是因?yàn)闄z波器2記錄到了高角度入射波產(chǎn)生的低頻虛假反射，這些虛假能量隨著偏移距增大而增強(qiáng)；而在檢波器2處，復(fù)頻移ADE－PML情況下的記錄與參考解吻合良好。這說明，復(fù)頻移拉伸算子對切入射波的吸收有明顯改善。

圖1 非復(fù)頻移（a）和復(fù)頻移（b）ADE－PML邊界條件下的聲壓場p的波場快照Fig.1 Snapshots of the pressure field pof acoustic wave in the non－CFS（a）and CFS（b）ADE－PML condition

圖2 非復(fù)頻移（a）和復(fù)頻移（b）ADE－PML條件下在直線x＝xs上接收到的地震記錄Fig.2 Seismograms recorded at line x＝xsin the non－CFS（a）and CFS（b）ADE－PML condition

圖3 非復(fù)頻移和復(fù)頻移ADE－PML邊界條件下得到的檢波點(diǎn)1（a）和檢波點(diǎn)2（b）的聲波方程地震記錄同精確解的對比Fig.3 Seismograms calculated in non－CFS ADE－PML and CFS ADE－PML conditions at receiver 1（a）and receiver 2（b）

圖4 顯式、隱式和半隱式更新ADE－PML輔助變量條件下在檢波器1（a）和檢波器2（b）處得到的地震記錄Fig.4 Seismograms calculated using the explicit，implicit and semi－implicit update scheme for the ADE－PML auxiliary variables at receiver 1（a）and receiver 2（b）

圖5 顯式、隱式和半隱式更新ADE－PML輔助變量條件下在檢波器1（a）和檢波器2（b）處得到的測試解同精確解的差異Fig.5 Differences between the test and reference solution calculated using the explicit，implicit and semi－implicit update scheme for the ADE－PML auxiliary variables at receiver 1（a）and receiver 2（b）

圖6 波場總能量衰減Fig.6 Total energy decay of the wavefield

如前文式（15）中所示，對于ADE－PML輔助變量的計(jì)算可使用顯式、隱式、和半隱式更新格式，這里分別對3種輔助變量的更新格式下的ADE－PML吸收效果進(jìn)行了檢驗(yàn)。如前述復(fù)頻移拉伸算子可以改善吸收效果，因此在復(fù)頻移ADE－PML邊界條件下計(jì)算。從圖4中可見，3種更新格式得到的地震記錄均無虛假反射振動，因?yàn)閺?fù)頻移拉伸算子的作用，在檢波點(diǎn)2處也與精確解吻合良好。圖5中繪出了3種更新格式同精確解的差異曲線，用于細(xì)微比較三者的差異。從圖中差異曲線的擺動幅度可見，由于數(shù)值頻散效應(yīng)在遠(yuǎn)偏移矩的積累，圖5b的差異曲線比圖5a中的擺動大，無論是垂直入射的檢波點(diǎn)1還是高角度入射的檢波點(diǎn)2處，顯式格式比隱式和半隱式格式同精確解的吻合程度好。

圖6為3種更新格式在長時間數(shù)值模擬中（這里進(jìn)行了2×105時間步，共240s的模擬計(jì)算）的波場總能量變化情況，從波場能量的變化情況可檢驗(yàn) ADE－PML吸收邊界的穩(wěn)定性［7，13］。求解區(qū)域Ω中（不包括吸收邊界區(qū)域）的總能量E通過下式計(jì)算：

從圖6a可見：當(dāng)數(shù)值模擬開始時（約0.1s左右），震源將能量注入介質(zhì)使得總能量急劇上升；之后由于吸收邊界作用，波場能量開始逐漸下降。在約1.7 s時能量出現(xiàn)一個陡降，因?yàn)檫@時聲波完全離開求解區(qū)域，理論上講該時刻之后求解區(qū)域中的能量為虛假能量。由圖6b可見，在數(shù)值模擬后期（6s之后），波場總能量持續(xù)下降至10－20J數(shù)量級，說明顯式、隱式和半隱式記憶變量更新格式都具有長時間穩(wěn)定性，并且仔細(xì)觀察可發(fā)現(xiàn)3種更新格式的能量衰減速度存在差異，隱式情況下能量衰減最快。

3 結(jié)論和討論

數(shù)值模擬表明，使用復(fù)頻移拉伸坐標(biāo)算子的ADE－PML比使用非復(fù)頻移拉伸算子（即傳統(tǒng)拉伸坐標(biāo)算子）的ADE－PML更能有效消除高角度入射情況下的虛假反射。對ADE－PML輔助變量的3種更新格式的檢驗(yàn)說明，3種格式在吸收效果上不存在明顯差異，但顯式條件下可以獲得相對較好的吸收效果；長時間的能量衰減計(jì)算表明，3種更新格式條件下的ADE－PML都穩(wěn)定至2×105時間步，說明ADE－PML在數(shù)值模擬中具有長時間穩(wěn)定性。由于ADE－PML技術(shù)具有非分裂場格式，降低了編程難度，且由于適應(yīng)高階精度時間格式，可以進(jìn)一步推廣到波動方程高精度偏移成像中。

（Reference）：

［1］Clayton R，Engquist B.Absorbing Boundary Conditions for Acoustic and Elastic Wave Equations［J］.Bulletin of the Seismological Society of America，1977，67（6）：1529－1540.

［2］Liao Z P，Wong H L，Yang B P，et al.A Transmitting Boundary for Transient Wave Analysis［J］.Scientia Sinica，1984，27（10）：1063－1076.

［3］Higdon R.Absorbing Boundary Conditions for Difference Approximations to the Multidimensional Wave Equation［J］.Mathematics of Computation，1986，47（176）：437－459.

［4］Bérenger J P.A Perfectly Matched Layer for the Absorption of Electromagnetic Waves［J］.Journal of Computational Physics，1994，114（2）：185－200.

［5］Chew W C，Weedon W H.A 3－D Perfectly Matched Medium from Modified Maxwell’s Equations with Stretched Coordinates［J］.Microwave and Optical Technology Letters，1994，7（13）：599－604.

［6］Drossaert F H，Giannopoulos A.A Nonsplit Complex Frequency Shifted PML Based on Recursive Integration for FDTD Modeling of Elastic Waves［J］.Geophysics，2007，72（2）：9－17.

［7］Komatitsch D，Martin R.An Unsplit Convolutional Perfectly Matched Layer Improved at Grazing Incidence for the Seismic Wave Equation［J］.Geophysics，2007，72（5）：155－167.

［8］Ramadan O.Auxiliary Differential Equation Formulation：An Efficient Implementation of the Perfectly Matched Layer［J］.IEEE Microwave and Wireless Components Letters，2003，13（2）：69－71.

［9］Moerloose J D，Stuchly M A.Behavior of Berenger’s ABC for Evanescent Waves［J］.IEEE Microwave and Wireless Components Letters，1995，5（10）：344－346.

［10］Moerloose J D，Stuchly M A.Reflection Analysis of PML ABC’s for Low－Frequency Applications［J］.IEEE Microwave and Guided Wave Letters，1996，6（4）：177－179.

［11］Kuzuoglu M，Mittra R.Frequency Dependence of the Constitutive Parameters of Causal Perfectly Matched Anisotropic Absorbers［J］.IEEE Microwave and Guided Wave Letters，1996，6（12）：447－449.

［12］Collino F，Tsogka C.Application of the PML Absorbing Layer Model to the Linear Elastodynamic Problem in Anisotropic Heterogeneous Media［J］.Geophysics，2001，66（1）：294－307.

［13］Festa G，Vilotte J P.The Newmark Scheme as Velocity－Stress Time－Staggering：An Efficient PML Implementation for Spectral－Element Simulations of Elastodynamics［J］.Geophysical Journal International，2005，161（3）：789－812.

［14］Zhang W，Shen Y.Unsplit Complex Frequency－Shifted PML Implementation Using Auxiliary Differential Equations for Seismic Wave Modeling［J］.Geophysics，2010，75（4）：T141－T154.

［15］周輝，謝春臨，王尚旭，等.復(fù)雜地表地震資料疊前深度偏移成像［J］.吉林大學(xué)學(xué)報：地球科學(xué)版，2012，42（1）：262－268.Zhou Hui，Xie Chunlin，Wang Shangxu，et al.Prestack Depth Migration for Geological Structures with Complicated Surface［J］.Journal of Jilin University：Earth Science Edition，2012，42（1）：262－268.

［16］孟慶生，樊玉清，張珂，等.高階有限差分法管波傳播數(shù)值模擬［J］.吉林大學(xué)學(xué)報：地球科學(xué)版，2011，41（1）：292－298.Meng Qingsheng，F(xiàn)an Yuqing，Zhang Ke，et al.Tube Wave Propagation Numerical Simulation Based on High Order Finite－Difference Method［J］.Journal of Jilin University：Earth Science Edition，2011，41（1）：292－298.