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

1中國石油大學(xué)（北京）地球物理與信息工程學(xué)院；油氣資源與探測國家重點實驗室；CNPC物探重點實驗室，北京 102249

2中海油研究總院，北京 100027

3Tulsa大學(xué)地球科學(xué)學(xué)院，Oklahoma，74104USA

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

1 引言

地震波數(shù)值正演是研究地震波在地下介質(zhì)傳播規(guī)律的重要手段（金勝汶等，1994；劉禮農(nóng)等，2004；謝桂生等，2005）.在使用數(shù)值方法模擬無限或半無限大介質(zhì)中的地震波傳播會因計算機(jī)內(nèi)存產(chǎn)生截斷邊界，從而產(chǎn)生虛假數(shù)值反射，吸收邊界條件技術(shù)可以用于消除截斷邊界效應(yīng)產(chǎn)生的虛假反射.經(jīng)過三十多年的發(fā)展，在地震波數(shù)值領(lǐng)域中主要產(chǎn)生了兩大類吸收邊界技術(shù)：海綿邊界條件（Cerjan etal.，1985；Sochacki etal.，1987）和旁軸近似邊界條件（Clayton and Engquist，1977；Stacey，1988；Higdon，1991；Quarteroni etal.，1988），但是這些吸收邊界技術(shù)在很多情況下表現(xiàn)不佳，比如對于大角度入射波和低頻能量不能很好的吸收.Bérenger（1994）在Maxwell電磁波方程中提出了完全匹配層（Perfectly Matched Layer，PML），用于消除時域有限差分法電磁波模擬中的虛假反射.PML對于較大范圍入射角和頻率成分的波具有高效的吸收作用，理論上對于入射波具有零反射系數(shù)，目前已經(jīng)成為常用的吸收邊界技術(shù).之后，PML很快被證明在地震波場數(shù)值模擬同樣高效（Chew and Liu，1996；Hasting etal.，1996）.

后來，PML從數(shù)學(xué)上證明可以同復(fù)拉伸坐標(biāo)的概念聯(lián)系起來（Chew and Weedon，1994）.但Kuzuoglu和Mittra（1996）指出，使用傳統(tǒng)復(fù)拉伸坐標(biāo)算子的PML對于隱矢波和極低頻波成分吸收效能會出現(xiàn)下降，并且在大角度特別是切入射情況下（如震源靠近PML吸收邊界的情況下）會產(chǎn)生較大程度的虛假反射，因此為解決該問題，Kuzuoglu和Mittra（1996）在電磁波模擬中引入了復(fù)頻移拉伸算子用于提高PML的吸收效能.在傳統(tǒng)PML算法中，為在時間域中求解復(fù)拉伸坐標(biāo)下的波動方程，需要對波場變量進(jìn)行非物理分裂，但經(jīng)典的分裂場形式的PML（Chew and Weedon，1994）不能直接使用復(fù)頻移拉伸算子.在PML實現(xiàn)形式上，一些非分裂場格式被提出（Cheng and Zhao，2011）.Roden和Gedney（2000）在Maxwell方程模擬中提出了卷積PML來實施復(fù)頻移拉伸算子，在PML方程中出現(xiàn)卷積項從而得到非分裂場的形式.之后，使用復(fù)頻移拉伸算子的卷積PML也被引入地震波模擬中（Drossaert and Giannopoulos，2007；Komatitsch and Martin，2007）.

Li和Dai（2008）在電磁波模擬中提出了匹配Z變換PML（Matched Z-transform PML，MZT-PML），同卷積PML是基于Fourier變換理論導(dǎo)出不同，MZT-PML是基于匹配Z變換理論得到非分裂場PML格式，也避免了卷積PML中的片斷常數(shù)近似（Roden and Gedney，2000），其推導(dǎo)過程也相對簡單（Li and Dai，2008）.具有非分裂場格式的MZTPML技術(shù)可以在地震波方程模擬中直接應(yīng)用復(fù)頻移拉伸算子（Shi etal.，2012），本文推導(dǎo)了適用于速度-壓力形式黏聲波方程的MZT-PML吸收邊界控制方程，并對相應(yīng)的數(shù)值正演結(jié)果進(jìn)行了吸收效果分析.

2 標(biāo)準(zhǔn)線性固體黏聲波方程

在實際介質(zhì)中，由于地層的吸收衰減作用引起的速度耗散，地震波頻率會隨著傳播距離的增大而降低，因此近年來地震勘探研究越對黏彈性介質(zhì)的關(guān)注不斷提高（Carcione etal.，1988；Liu and Chang，1996）.常用的黏彈性介質(zhì)模型有Kelvin-Voigt模型、Maxwell模型和標(biāo)準(zhǔn)線性固體（Standard Linear Solid，SLS）模型，但Maxwell模型不能描述彈性蠕變，Kelvin-Voigt不能描述應(yīng)力弛豫.SLS模型（Carcione，2007）可以同時描述上述兩種現(xiàn)象，因此相對符合實際黏彈性介質(zhì)中的地震波傳播規(guī)律.本文主要目的在于說明PML吸收邊界的應(yīng)用，因此為簡便起見使用基于單SLS模型的黏聲波方程.在二維Cartesian直角坐標(biāo)系中，速度-壓力格式的

SLS模型黏聲波方程如下（Carcione，2007）：

其中，ω0為震源的中心角頻率.由（1）式可見，當(dāng)弛豫時間常數(shù)τε＝τσ時，（1）式即退化為聲波方程.

3 匹配Z變換完全匹配層

將（1）式變換到頻率域，有

其中，i為虛數(shù)單位，ω為角頻率，根據(jù)復(fù)拉伸坐標(biāo)變換（Chew and Weedon，1994）的概念，將空間偏導(dǎo)數(shù)?η（η＝x，z）替換為復(fù)拉伸坐標(biāo)?~η＝S-1η?η.其中Sη為拉伸坐標(biāo)算子，對PML的吸收特性有重要影響.傳統(tǒng)拉伸坐標(biāo)算子為

其中，實系數(shù)dη≥0，使波場振幅在PML內(nèi)部呈指數(shù)級別衰減.使用傳統(tǒng)拉伸算子（4）式的PML存在大角度入射波的吸收效能不足的問題，這是因為傳統(tǒng)拉伸算子只對行波具有吸收作用，而對于沿求解區(qū)域和PML區(qū)域界面平行方向傳播的隱矢波不起作用，另外對于極低頻成分也存在奇異性因而也不能較好地消除極低頻反射（Chew and Liu，1996）.為解決該問題，復(fù)頻移拉伸坐標(biāo)算子被提出：

式（5）中通過引入輔助衰減系數(shù)κη≥1和頻移系數(shù)αη≥0，將復(fù)平面的極點從從實軸移動到負(fù)虛半平面上，以提高對隱矢波和低頻波的吸收.易見（5）式比（4）式更具有一般性，令κη＝1和αη＝0即為傳統(tǒng)拉伸算子（4）式.在復(fù)拉伸坐標(biāo)下，（3）式改寫為因為拉伸坐標(biāo)算子Sη為頻率ω的函數(shù)，因此將（6）式變換回時間域會產(chǎn)生卷積項（Roden and Gedney，2000）.考慮到頻率域的乘法運算和時間域的卷積運算都等價于Z域的乘法運算，因此將（6）式變換到Z域可以使推導(dǎo)過程相對簡單.

MZT-PML的主要思想在于將拉伸坐標(biāo)算子（5）式的倒數(shù)分裂為兩部分的乘積，且令iω僅存在于其中一項中：

利用對應(yīng)關(guān)系iω?s將（6）式變換到S域，有

其中，算符ZT［］·表示匹配Z變換運算，Δt表示時間采樣間隔.為簡短起見，這里只以式（6a）的MZTPML格式為例進(jìn)行說明；另外這里假設(shè)PML吸收邊界由η＝0（η＝x，z）起始.將（6a）由頻率域變換到Z域有：

其中，（1-z-1）／Δt是iω的Z變換，考慮映射關(guān)系iω??t和?t?（1-z-1）／Δt，將（9）式代入（10）式可得

將（12）轉(zhuǎn)化為

式（15）即為關(guān)于聲壓場p連同輔助變量Ψvx，x和Ψvz，z的有限差分時間更新格式.Ψvx，x和Ψvz，z需要在更新速度場v和壓力場p前計算，注意到式（15c）中同時含有不同離散時刻（n和n+1時間步）的Ψ，因此Ψn應(yīng)當(dāng)在更新前暫存于一個臨時變量中.（6）式中關(guān)于速度場vx、vz和記憶變量R的MZT-PML格式可以同理導(dǎo)出.

4 數(shù)值正演

在數(shù)值模擬時，空間上采用八階交錯網(wǎng)格有限差分格式（Virieux，1986），并采用比Taylor展開系數(shù)精度更高的Holberg差分系數(shù)（Holberg，1987）來計算所有空間導(dǎo)數(shù).Holberg系數(shù)通過在給定階數(shù)截斷計算得到，降低了群速度產(chǎn)生的誤差從而整體上減小了數(shù)值頻散.例如，x方向空間偏導(dǎo)數(shù)在i+1／2網(wǎng)格點處的差分離散格式為

其中，八階精度Holberg差分系數(shù)C1＝1.231666，C2＝-0.1041182，C3＝0.02063707，C4＝-0.003570998.模型計算區(qū)域大小為5250m×1500m的長方形區(qū)域，由701×201個網(wǎng)格點構(gòu)成，有限差分元胞邊長為Δx＝Δz＝7.5m.介質(zhì)密度ρ＝2000kg·m-3，縱波速度cp＝3000m／s，品質(zhì)因子Q＝30.時間上，采用二階顯式差分格式，因此時間更新步長Δt由需滿足Courant穩(wěn)定性條件：

取Δt＝1.2ms，模擬持續(xù)時長4.8s，共計4000時間步.采用點聲壓震源激發(fā)，震源子波為一導(dǎo)數(shù)高斯函數(shù)其中震源中心頻率f0＝20Hz，子波時延t0＝1.2／f0＝0.06s，激發(fā)點位于（517.5m，180m）處；這里設(shè)置2個檢波點，檢波點1和檢波點2分別位于（742.5m，1380m）和（5017.5m，142.5m）處.計算區(qū)域四周由PML吸收邊界包圍，厚度為75m，占10個有限差分元胞.故震源靠近頂部吸收邊界，檢波點1和2分別靠近底部和頂部吸收邊界（圖1）.一般在PML區(qū)域內(nèi)，要將復(fù)頻移算子（5）式中各參數(shù)設(shè)置為沿PML區(qū)域和求解區(qū)域的交界面外法線方向漸進(jìn)變化：

其中，η為PML內(nèi)一點到與交界面的距離，L為PML的厚度.各衰減參數(shù)取值如下（Li and Matar，2010；Festa and Vilotte，2005）：d0＝-（Pd+1）cplog10（Rc）／2L，κ0＝-（Pκ+1）blog10（Rc）／2L，α0＝πf0；其中Rc＝10-5，Pd＝Pκ＝2，Pα＝1，b＝10Δh，Rc為理論反射系數(shù)，Δh為空間元胞邊長.這里也給出使用復(fù)拉伸坐標(biāo)算子（3）式的傳統(tǒng)分裂場PML吸收邊界條件下的黏聲波方程模擬（Liu and Tao，1997）結(jié)果，用以分析對比.另外，為檢驗吸收邊界效果需要同參考解對比，本文參考解是在不含吸收邊界的足夠大區(qū)域（模擬時間內(nèi)邊界反射不會到達(dá)檢波器）使用相同的數(shù)值離散格式計算得到，震源和檢波器相對位置不變.

圖1 計算區(qū)域示意圖圖中黑色內(nèi)框線表示PML與求解區(qū)域的交界，圓圈表示震源位置，三角表示檢波點位置Fig.1 Diagrammatic figure of the computational model The inner black lines denote the interface between PML and solution domain，the circle denotes the source location and the triangles denote receiver locations

圖2 傳統(tǒng)PML邊界條件下（a）0.48s，（b）0.96s和（c）1.68s的黏聲波聲壓場p的波場快照（振幅已加10倍增益）Fig.2 Snapshots of the pressure field pof viscoacoustic wave at（a）0.48s，（b）0.96sand（c）1.68s in the classical PML condition，and the amplitudes has been scaled by a factor of 10

圖3 MZT-PML邊界條件下（a）0.48s，（b）0.96s和（c）1.68s的黏聲波聲壓場p的波場快照（振幅已加10倍增益）Fig.3 Snapshots of the pressure field pof viscoacoustic wave at（a）0.48s，（b）0.96sand（c）1.68s in the MZT-PML condition，and the amplitudes has been scaled by a factor of 10

該算例中，由于震源靠近頂部吸收邊界，形成了大角度入射PML的情況.從波場快照圖2可見，傳統(tǒng)PML對大角度入射波吸收不完全，產(chǎn)生了低頻虛假反射，且虛假反射能量隨著入射角度和傳播距離的增大而增強(qiáng).從1.68s的波場可見，在遠(yuǎn)偏移距低頻反射尤為明顯.而圖3中MZT-PML邊界由于使用復(fù)頻移拉伸算子，對大角度入射波吸收很好，沒有出現(xiàn)圖1中的低頻反射；在圖3中1.68s波場中，即使在遠(yuǎn)偏移距也沒有可見的虛假反射.

從圖4a可見，在檢波點1處無論是傳統(tǒng)PML還是MZT-PML的計算結(jié)果都和參考解吻合良好，說明在小角度入射情況下，兩種吸收邊界都可以高效吸收外行波場能量.而從圖4b中檢波點2的地震記錄可見，PML情況下的地震記錄上產(chǎn)生了明顯的虛假反射抖動，與參考解吻合不佳，這是因為檢波器2記錄到了大角度入射波產(chǎn)生的低頻虛假反射，這些虛假反射在遠(yuǎn)偏移距更強(qiáng).而在檢波器2處，MZT-PML情況下的記錄與參考解吻合良好，說明MZT-PML相對于傳統(tǒng)PML，對切入射波的吸收有明顯改善.這里也給出完全彈性介質(zhì)情況下（令弛豫時間常數(shù)τε＝τσ）聲波方程在MZT-PML邊界條件下的模擬結(jié)果：圖5中，對比聲波方程和黏聲波方程的地震記錄，由于介質(zhì)的黏滯性，導(dǎo)致聲波振幅衰減強(qiáng)烈（如圖5b）.但MZT-PML邊界條件下，無論是聲波方程還是黏聲波方程的計算結(jié)果都與各自參考解吻合良好.

另外，長時間數(shù)值正演中的穩(wěn)定性也是吸收邊界效能的一個重要指標(biāo)（Komatitsch and Martin，2007；Chen and Bording，2010；Chen，2011），本文將波場傳播時間延長至120s（共105時間步），并記錄求解區(qū)域Ω中（不包括PML吸收邊界）的波場總能量E為

圖4 傳統(tǒng)PML和MZT-PML邊界條件下得到的（a）檢波點1和（b）檢波點2的黏聲波方程地震記錄Fig.4 Seismograms calculated in classical PML and MZT-PML conditions at（a）receiver 1and（b）receiver 2

圖5 MZT-PML邊界條件下得到的（a）檢波點1和（b）檢波點2的聲波、黏聲波方程地震記錄Fig.5 Seismograms of acoustic and viscoacoustic wave calculated in MZT-PML condition at（a）receiver 1and（b）receiver 2

圖6 吸收邊界條件下黏聲波場的總能量衰減（a）0～6s（b）0～120s，縱坐標(biāo)軸為半對數(shù)坐標(biāo)顯示Fig.6 Total energy decay of the viscoacoustic wave field on a semi-logarithmic scale for y-axis in the figure（a）0～6s（b）0～120s

從圖6中展示了彈性和黏彈性介質(zhì)中兩種吸收邊界條件下的波場總能量衰減情況.由圖6a可見，當(dāng)數(shù)值正演開始時（約0.1s左右），震源將能量注入介質(zhì)使得總能量急劇上升，之后由于黏彈性介質(zhì)本身的衰減，使得黏聲波情況下的總量明顯低于彈性聲波.由于PML吸收邊界的作用，也導(dǎo)致兩種介質(zhì)中的波場能量開始逐漸下降.在約1.7s時能量出現(xiàn)一個陡降，因為這時直達(dá)波完全離開求解區(qū)域，理論上講該時刻之后求解區(qū)域中的能量全部是虛假能量，這時的能量衰減情況可以較為準(zhǔn)確地反映吸收邊界的效能，從圖6a可見MZT-PML情況下波場剩余能量低于PML情況（紅藍(lán)色曲線之間，以及黑色實虛線之間的對比）.由圖6b可見，在數(shù)值正演后期（6s之后），能量仍舊持續(xù)下降，說明MZTPML在彈性聲波、黏聲波方程模擬中具有長時間數(shù)值穩(wěn)定性.在黏聲波介質(zhì)中，MZT-PML情況下總能量持續(xù)下降至10-30J數(shù)量級，低于PML情況下的10-13J數(shù)量級.

5 討論和結(jié)論

本文基于匹配Z變換的思想建立了用于黏聲波方程正演的非分裂MZT-PML，這種非分裂PML技術(shù)可以直接使用復(fù)頻移拉伸算子.在有限差分黏聲波方程的數(shù)值正演實驗中，同傳統(tǒng)分裂場PML邊界條件的應(yīng)用效果進(jìn)行了比較.數(shù)值實驗結(jié)果說明，同傳統(tǒng)PML相比，當(dāng)震源靠近吸收邊界時，MZT-PML可以明顯改善對大角度入射波的吸收.另外，MZT-PML在完全彈性介質(zhì)和黏彈性介質(zhì)情況下均表現(xiàn)良好.從105時間步的能量衰減計算可見，MZT-PML邊界條件下波場總能量在數(shù)值正演后期總體上呈持續(xù)下降趨勢，說明MZT-PML在黏聲波方程模擬中具有長時間穩(wěn)定性.該方法具有比傳統(tǒng)PML更好的吸收衰減效能.

Bérenger J P.1994.A perfectly matched layer for the absorption of electromagnetic waves.J.Comput.Phys.，114（2）：185-200.

Carcione J M，Kosloff D，Kosloff R.1988.Wave propagation simulation in a linear viscoelastic medium.Geophys.J.Int.，95（3）：597-611.

Carcione J M.2007.Wave Fields in Real Media：Theory and Numerical Simulation of Wave Propagation in Anisotropic，Anelastic，Porous and Electromagnetic Media.2nd edn.Amsterdam，The Netherlands：Elsevier Science.

Cerjan C，Kosloff D，Kosloff R，etal.1985.A nonreflecting boundary condition for discrete acoustic and elastic wave equations.Geophysics，50（4）：705-708.

Chen J Y，Bording R P.2010.Application of the nearly perfectly matched layer to the propagation of low-frequency acoustic waves.J.Geophys.Eng.，7（3）：277-283.

Chen J Y.2011.Application of the nearly perfectly matched layer for seismic wave propagation in 2Dhomogeneous isotropic media.Geophysical Prospecting，59（4）：662-672.

Cheng J Y，Zhao J G.2011.Application of the nearly perfectly matched layer to seismic-wave propagation modeling in elastic anisotropic media.B.Seismol.Soc.Am.，101（6）：2866-2871.

Chew W C，Liu Q H.1996.Perfectly matched layers for elastodynamics：A new absorbing boundary condition.J.Comput.Acoust.，4（4），341-359.

Chew W C，Weedon W H.1994.A 3Dperfectly matched medium from modified Maxwell′s equations with stretched coordinates.Microw.Opt.Tech.Lett.，7（13）：599-604.

Clayton R，Engquist B.1977.Absorbing boundary conditions for acoustic and elastic wave equations.B.Seismol.Soc.Am.，67（6）：1529-1540.

Drossaert F H，Giannopoulos A.2007.Complex frequency shifted convolution PML for FDTD modelling of elastic waves.Wave Motion，44（7-8）：593-604.

Engquist B，Majda A.1977.Absorbing boundary conditions for the numerical simulation of waves.Mathematics of Computation，31（139）：629-651.

Festa G，Vilotte J P.2005.The Newmark scheme as velocity–stress time-staggering：an efficient PML implementation for spectral element simulations of elastodynamics.Geophys.J.Int.，161（3）：789-812.

Hasting F D，Schneider L B，Broschat S L.1996.Application of the perfectly matched layer（PML）absorbing boundary condition to elastic wave propagation.J.Acoust.Soc.Am.，100（5）：3061-3069.

Higdon R L.1991.Absorbing boundary conditions for elastic waves.Geophysics，56（2）：231-241.

Holberg O.1987.Computational aspects of the choice of operator and sampling interval for numerical differentiation in large-scale simulation of wave phenomena.Geophysical Prospecting，35（6）：629-655.

Jin S W，Chen B Y，Ma Z T.1994.Three-dimensional wave equation forward modeling by the finite-difference methods.Acta Geophysica Sinica（in Chinese），37（6）：804-810.

Komatitsch D，Martin R.2007.An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation.Geophysics，72（5）：SM155-SM167.

Kuzuoglu M，Mittra R.1996.Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers.IEEE Microwave Guided Wave Letters，6（12）：447-449.

Li J X，Dai J F.2008.Modified Z-transform-based FDTD algorithm for the anisotropic perfectly matched layer.International Journal of Numerical Modelling：Electronic Networks，Devices and Fields，21（5）：279-286.

Li Y F，Matar O B.2010.Convolutional perfectly matched layer for elastic second-order wave equation.J.Acoust.Soc.Am.，127（3）：1318-1327.

Liu L N，Cui F L，Zhang J F.2004.Seismic modeling with one-way wave equation in 3-D complex structures.Chinese J.Geophys.（in Chinese），47（3）：514-520.

Liu Q H，Chang C.1996.Compressional head waves in attenuative formations：forward modeling and inversion.Geophysics，61（6）：1908-1920.

Liu Q H，Tao J P.1997.The perfectly matched layer for acoustic waves in absorptive media.J.Acoust.Soc.Am.，102（4）：2072-2082.

Quarteroni A，Tagliani A，Zampieri E.1988.Generalized Galerkin approximations of elastic waves with absorbing boundary conditions.Comput.Method.Appl.M.，163（1-4）：323-341.

Roden J A，Gedney S D.2000.Convolution PML（CPML）：An efficient FDTD implementation of the CFS-PML for arbitrary media.Microw.Opt.Tech.Lett.，27（5）：334-339.

Shi RQ，Wang SX，Zhao JG.2012.An unsplit complex-frequencyshifted PML based on matched Z-transform for FDTD modelling of seismic wave equations.J.Geophys.Eng.，9（2）：218-229.Sochacki J，Kubichek R，George J，etal.1987.Absorbing boundary conditions and surface waves.Geophysics，52（1）：60-71.

Stacey R.1988.Improved transparent boundary formulations for the elastic-wave equation.B.Seismol.Soc.Am.，78（6）：2089-2097.

Virieux J.1986.P-SV wave propagation in heterogeneous media：velocity-stress finite-difference method.Geophysics，51（4）：889-901.

Xie G S，Liu H，Li Y M，etal.2005.Seismic modeling by reflection／transmission operator and one-way wave equation under the condition of fluctuating reflectors.Chinese J.Geophys.（in Chinese），48（5）：1172-1178.

附中文參考文獻(xiàn)

金勝汶，陳必遠(yuǎn)，馬在田.1994.三維波動方程有限差分正演方法.地球物理學(xué)報，37（6）：804-810.

劉禮農(nóng)，崔風(fēng)林，張劍鋒.2004.三維復(fù)雜構(gòu)造中地震波模擬的單程波方法.地球物理學(xué)報，47（3）：514-520.

謝桂生，劉洪，李幼銘等.2005.界面起伏條件下反射／透射算子+單程波方程的地震波模擬方法.地球物理學(xué)報，48（5）：1172-1178.