田 坤 李振春 黄建平 李 娜 孔 雪 刘玉金
（中国石油大学（华东）地球科学与技术学院，山东青岛 266555）
基金项目：国家973课题（编号2011CB202402），石油大学创新基金（编号27R1001046A）及（Y090104）
作者简介：田坤，男，在读博士研究生，现从事地震波正演模拟研究。Email：tiankunwudi＠163.com。
摘 要：利用波动方程进行数值模拟时，由于介质的计算范围肯定是有限的，这样就会人为地造成 计算边界，不可避免地产生边界反射，因此需要使用吸收边界条件使入射到边界的能量被吸收掉，减小 边界反射所带来的影响。PML吸收层技术已经被证明是非常有效的边界吸收技术，对体波和面波的吸收 都具有非常好的效果，已经被广泛应用于弹性波的数值模拟中。但是在大角度入射即掠射的情况下传统 的PML技术还是存在一定的问题，掠射情况下衰减不够，反射系数比较大，离散后会产生比较严重的 假反射，降低吸收效果，对真实波场产生比较严重的影响，这种影响随入射角、偏移距的增加而增大，而且会对实际应用中的成像、反演等其他处理产生不利影响。而掠射是普遍存在的，比如薄片区域、震 源位置接近研究区域边缘、大偏移距接收等等。本文基于一阶速度-应力方程，提出了一种卷积完全匹 配层（CPML）技术，并采用交错网格有限差分方法对其进行了实现，结果表明CPML技术有效改善了 掠射情况下的吸收效果，尤其是对于散射场等弱能量波场，相对来说效果更是明显。并且在实现过程中 不用分裂变量，应用更加方便简单，易于编程实现，计算卷积时采用递推的形式，不会增加计算量，存 储量也没有太大的变化。
关键词：PML；吸收边界条件；弹性波；数值模拟；卷积完全匹配层技术
An Improved Perfectly Matched Layer Absorbing Boundary Condition for the Carbonate Karst Reservoir
Tian Kun，Li Zhenchun，Huang Jianping，Li Na，Kong Xue，Liu Yujin
（School of Geoscience，China University of Petroleum（East China），Qingdao 266555，China）
Abstract：When we simulate from wave equation，it's can cause computing bounds artificially and the reflections from artificial boundaries inevitably，because the calculation range is limited so it need absorbing boundary conditions to make the energy that incidenting on the bound is absorbed and minimal the impact caused by boundary reflections.The perfectly matched layer（PML）absorbing boundary condition has proven to be a very efficient method for the numerical forward modeling of elastic wave equation to absorb both body waves and surface waves.But at grazing incidence the classical discrete PML method still exist some problems，that is the attenuation is not enough，the reflection coefficient is big to some extent，these canproduce severe spurious reflection after dispersion and reduce the effect to absorb，then cause seriously influences to the real field and these influences are growing with increasing incident angle and offset，besides，these can negatively affect the imaging and inversion and other treatment in practical application.But grazing incidence is natural，for instance in the case of very thin mesh slices，as well as sources located close to the edge of the mesh，or in the case of receivers located at very large offset.In this paper，we present an improved PML absorbing boundary condition at grazing incidence which is called convolution PML（CPML）based on the first-order velocity-stress partial differential wave equation.And it is implemented by a staggered grids finite-difference method.The results show that the CPML technique canobtain a better result in the case of grazing incidence，especially for scattering field and such like this weak reflection wave field.And its implementation is convenient and easy to be used and programmed because of the un-split variable.The cost of the calculation is not increase because of the recursion when convolution is computed.The memory storage is also similar to that of the classical PML
Key words：PML；absorbing boundary condition；elastic wave；Forward modeling；CPML
引言
由于计算能力的快速发展，在过去几十年里对地震波在复杂介质中传播的数值模拟方法的研究也越 来越广泛和深入。其中应用最广泛的是有限差分法［1～3］，而在有限区域的有限差分正演模拟中，无论是 求解麦克斯韦电磁波方程组，还是求解弹性波动力学方程组，Bérenger（1994）提出的完全匹配层（perfectly matched layer，PML）吸收边界条件是最有效的吸收边界条件［4］。完全匹配层中的波动方程可 以看做是常规的波动方程的推广，波在其中传播时相位改变而振幅随指数衰减；PML和主研究区域弹 性参数相同而衰减系数不同，波阻抗完全匹配，理论上不会发生反射。Chew and Weedon（1994）引入 复数伸展坐标系对PML吸收边界条件进行公式化［5］；Rappapport（1995）证明了PML介质等价于在吸 收边界区域引入各向异性介质［6］；大量研究［5，7］表明，PML吸收边界条件比指数衰减吸收边界条件［8，9］、 廖氏吸收边界条件［10］、Higdon吸收边界条件［11］和旁轴近似吸收边界条件［12］等具有更优越的吸收性能。Chen et al.（1997）和Wang and Oristaglio（2000）将PML吸收边界条件成功应用到电磁波方程组的求 解中［7，13］。近年来，PML吸收边界条件也被应用到声波和弹性波的有限差分正演模拟中［13～16］。另外，Teixeir and Chew（1999）和Chew and Liu（1996）在柱坐标系和球坐标系中实现了PML吸收边界 条件［17，18］。
但是传统的PML也存在一定的缺陷，离散后的PML层反射系数不严格为0，尤其在大入射角 的情况下更为严重，在这种情况下，相当一部分能量以反射波的形式被送回到主研究区域，而大 入射角的情况是普遍存在的，比如薄片区域、震源位置接近研究区域边缘、大偏移距接收等等。为了克服这个问题，一种修改复坐标变换以改进离散后大入射角吸收效果的PML方法被提了出来 并被应用到麦克斯韦方程组的求解中（Roden and Gedney，2000）并被命名为卷积完全匹配层（CPML）［19］。本文采用这种方法来求解弹性波方程，结果表明CPML技术有效改善了掠射的吸收 效果，并且在实现过程中不用分裂变量，应用更加方便简单，计算卷积时采用递推的形式，不会 增加计算量，存储量也没有太大的变化。
1 CPML基本原理
PML技术本质上是将波动方程在PML层内进行复坐标变换，对于变换后的坐标，方程及其解的形 式是不变的，但是对于原坐标解是衰减的。传统的PML中频率域坐标变换为（以x为例）：
国际非常规油气勘探开发（青岛）大会论文集
在CPML中对（3）式进行了扩展，使其形式更一般化：
国际非常规油气勘探开发（青岛）大会论文集
其中αx≥0，κx≥1，可以看出传统PML是CPML在αx＝0，κx＝1情况下的特例。
将坐标变换关系变换回时域，用 表示1/sx的傅里叶反变换，可以得到
国际非常规油气勘探开发（青岛）大会论文集
根据（4）式可得
国际非常规油气勘探开发（青岛）大会论文集
δ与1和e-atH（t）与1/（a＋iω）是傅里叶变换对，所以
国际非常规油气勘探开发（青岛）大会论文集
这样CPML中坐标变换关系就分为两部分，其中第一部分容易计算，只要计算空间偏导除以系数κ 就可以了，最主要的是计算第二部分的卷积。在离散的交错网格下将n时刻的卷积写成
国际非常规油气勘探开发（青岛）大会论文集
交错网格的情况下可以写为
国际非常规油气勘探开发（青岛）大会论文集
其中
国际非常规油气勘探开发（青岛）大会论文集
式中：
国际非常规油气勘探开发（青岛）大会论文集
因为（12）式是简单的指数形式，所以可以把（11）式写成递推的形式：
国际非常规油气勘探开发（青岛）大会论文集
这样在主研究区域进行正常计算，在PML层内通过下面的关系进行坐标变换后的计算就可以了：
国际非常规油气勘探开发（青岛）大会论文集
其中Ψx可以通过（13）式的递推得到。
这就是CPML的基本原理，可以很容易地推广到y，z方向。角上的情形与传统的PML类似，所有 方向同时考虑就可以了。
2 正演模拟计算与对比
为了验证CPML的有效性与优越性，下面对各向同性介质进行了试算，并与传统的PML的计算结 果进行了对比。图1是层状介质模型，采样点为301×301，网格间距5m×5m，上层纵横波速度分别为 3000m/s，1600m/s，下层纵横波速度分别为3300m/s，1900m/s，密度均为2800kg/m3，震源坐标（150，5），道间距5m，采样率0.5ms，采样时间1.5s，图2和图3分别是主频为20Hz的z分量炮记录 和波场快照对比，可以清楚地看到PML在上界面产生了比较强的反射，而且对炮记录有所影响，而 CPML的效果要更好一些。
图4和图5分别是主频为30Hz的z分量炮记录和波场快照对比，也可以看出CPML比PML的 吸收效果要更好一些。图6是3个单道的波形对比图，单道1、2、3坐标分别是（201，5），（251，5），（301，5），图6左侧是0～1.2s的波形对比，右侧是0.3～1.2s的波形对比，也可以 看出CPML要比PML效果好，而且随着偏移距增大，入射角增加，对比越明显，传统的PML中的 反射影响越严重。
图1 层状介质模型
图2 主频为20Hz的炮记录对比
图3 主频为20Hz的0.3s波场快照（上）及局部放大图（下）对比
图4 主频为30Hz的炮记录对比
图7是含散射体介质模型，在（225，15）处有一散射体，散射体纵横波速度分别为2700m/s，1300m/s，其他参数与前相同。图8和图9分别是主频为30Hz的z分量炮记录和波场快照对比，可以看 出PML的虚假反射对散射场有相当的影响，这会对实际应用中的成像、反演等其他处理产生不利影响。
图5 主频为30Hz的0.3s波场快照（上）及局部放大图（下）对比
3 结论
本文基于不分裂的卷积方法在交错网格有限差分情形下对传统的PML技术进行了一般化扩展，推 导了CPML的时域坐标变换关系的递推公式，并通过模型试算进行了对比分析，得到以下几点结论： (1)传统的PML在离散后反射系数不为0，尤其在大入射角、大偏移距的情况下会产生很强的假反射，降低吸收效果，对真实波场产生比较严重地影响，这种影响随入射角、偏移距的增加而增大；(2)改进 后的CPML技术能够增强掠射情况的波场衰减，降低反射系数，有效改善匹配层的吸收效果；(3)对于 含散射体的复杂介质模型，传统的PML也会产生假反射，并更加严重地影响散射场等弱能量波场，而 且会对实际应用中的偏移、反演等处理产生不利影响，而CPML技术同样可以改善吸收效果，降低不利 影响；(4)CPML技术采用不分裂变量的卷积方法进行计算，并通过递推来计算卷积，使其存储量和计算 量与传统PML相比都没有太大变化，不会增加计算成本。
本文的CPML技术能够有效改善掠射的吸收效果，而且不增加计算成本，这对很多如薄片区域、震 源接近研究区域边缘、大偏移距接收等实际情况有很好的应用前景。对于偏移成像等处理中的吸收边界 条件，这也是一种比较好的选择。但是CPML没有改变传统的PML的基本思想，所以它还是存在一些 问题，比如对于一些各向异性介质有固有的不稳定性、离散后匹配层反射系数不严格为零等。本文在匹 配层外围采用的是Dirichlet边界条件，后续研究中可以将其替换为旁轴吸收边界条件以进一步提高离散 后的吸收效果，另外对于二阶位移波动方程的CPML吸收边界条件也有待进一步研究。CPML技术以及 后续相关方法的进一步研究，将有利于西部碳酸盐岩探区复杂近地表速度模型下的正演模拟，尤其是对 碳酸盐岩探区近地表散射波机理认识的研究，可为将来碳酸盐岩探区勘探开发服务。
图6 不同单道的CPML和PML波形对比图
图7 含散射体介质模型
图8 含散射体模型的z分量炮记录对比
图9 含散射体模型的z分量0.3s波场快照（上）及局部放大图（下）对比
参考文献
［1］Alterman Z and Karal F C.Propagation of elastic waves in layered media by finite difference methods.Bulletin of the Seismological Society of America，1968，58：367-398.
［2］Madariaga，R..Dynamics of an expanding circular fault.Bulletin of the Seismological Society of America，1976，65： 163-182.
［3］Virieux，J.P-SV wave propagation in heterogeneous media：velocity-stress finite-difference method.Geophysics，1986，51： 889-901.
［4］Bérenger，J.P .A perfectly matched layer for the absorption of electromagnetic waves.Journal of Computational Physics，1994，114：185-200.
［5］Chew W.C.，and W.H.Weedon.A 3-D perfectly matched medium from modified Maxwell's equations with stretched coordinates.Microwave and Optical Technology Letters，1994，7：599-604.
［6］Rappaport，C.M..Perfectly matched absorbing boundary conditions based on anisotropic lossy mapping of space.IEEE Microwave Guided Wave Lett.，1995，5：90-92.
［7］Chen，Y.H.，Chew，W.C.，Oristaglio，M.L..Application of perfectly matched layers to the transient modeling of subsurface EM problems.Geophysics，1997，62：1730-1736.
［8］Marfurt，K.J..Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations.Geophysics，1984，49：533-549.
［9］Shin，C..Sponge boundary condition for frequency domain modeling.Geophysics，1995，60：1870-1874.
［10］Liao，Z.P.，Wong，H.L.，Yang，B.P.，et al..A transmitting boundary for transient wave analysis.Scientia Sinica，1984，27（10）：1063-1076.
［11］Higdon，R.L.Absorbing boundary condition for elastic waves.Geophysics，1991，56：231-241.
［12］Engquist，B.，Majda，A..Absorbing boundary conditions for the numerical simulation of waves.Mathematics of Computation，1977，31（9）：629-651.
［13］Wang，T.，Oristaglio，M.L.3D simulation of GPR surveys over pipes in dispersive soils.Geophysics，2000，65：1560- 1568.
［14］Liu，Q.H.，Tao，J.P..The perfectly matched layer for acoustic waves in absorptive media.Acoust.Soc.Am.，1997，102：2072-2082.
［15］Zeng，Y.，He，J.，Liu，Q.H..The application of the perfectly matched layer in numerical modeling of wave propagation in poroelastic media.Geophysics，2001，66：1258-1266.
［16］Imhof，M.G..Calculating the seismic effect of 3D underground structures and topography with the finite-difference method.72nd Internat.Mtg.，Soc.Expl.Geophys.，Expanded Abstracts，2002：1939-1942.
［17］Teixeira，F.L.，Chew，W.C..On causality and dynamic stability of perfectly matched layers for FDTD simulations.IEEE Transactions on Microwave Theory and Techniques，1999，47：775-785.
［18］Chew，W.C.，Liu，Q.H..Perfectly matched layers for elastodynamics：A new absorbing boundary condition.Journal of Computational Acoustics，1996，4：341-359.
［19］Roden，J.A.，and S.D.Gedney.Convolution PML（CPML）：An efficient FDTD implementation of the CFS-PML for arbitrary media.Microwave and Optical Technology Letters，2000，27：334-339.
内容全部来源于网络收集，如有侵权，请联系网站删除：QQ：24596024