随着计算机技术的飞速发展, 越来越多的数值方法被应用于基于波动方程的地震波场数值模拟中, 主要包括有限差分法^[1]、有限元法^[2]、伪谱法^[3]等, 这些方法的应用提升了逆时偏移成像^[4-7]、全波形反演^[8-9]等技术解决实际问题的能力。高精度地震波场的模拟成为这些新技术的基础, 也是众多地球物理学家研究的重点^[10-11]。

实际情况下的地震波场是在半无限区域传播的, 而数值模拟只能对有限区域模拟计算, 这样在边界内侧和外侧形成了具有不同波阻抗的两种介质, 计算边界相当于一个反射面, 当波场能量传播到计算边界时会给计算区域带来虚假反射, 从而导致错误的结果, 因此需要人为设置边界介质将能量吸收, 避免能量反射。在过去几十年, 发展了海绵吸收^[12]、最优化条件^[13]、特征值分解法^[14]、球面精确吸收条件^[15]等吸收边界技术, 但是, 这些技术在应用中受到一定的限制, 如在地震波为掠射波情况下会产生虚假反射、所有角度入射会产生低频反射, 而且实施过程中存在大量数值计算。20世纪90年代, BERENGER^[16]为了对电磁波数据精确建模提出了一种新的边界条件, 称之为完全匹配层(perfectly matched layer, PML)吸收边界条件, 其基本思想是:在吸收边界区域匹配与计算区域相同波阻抗, 使入射波无反射的进入吸收层, 吸收层中是衰减介质, 使得入射波迅速衰减直至消除。理论上PML边界能够吸收任何入射角度和频率的入射波, 实践也证明PML吸收边界条件十分有效并被成功应用于声波^[17-20]和弹性波^[21]的研究。

如果在波场延拓的过程中采用基于波动方程时域有限差分进行求解, 那么其中的离散化处理会降低PML吸收边界的效果, 尤其在掠入射的情况下, 虚假反射更加明显^[22], 使得PML在对薄层介质、震源位置靠近网格边缘或者检波器处于远偏移距的模型的数值模拟效果并不理想。为了解决上述问题, KUZUOGLU等^[2]提出了一种复频移方法, 这种方法被RODEN等^[23]称为褶积完全匹配层(convolutional perfectly matched layer, CPML)吸收边界。

在前人的研究基础上, 本文对CPML吸收边界下的一阶速度-应力方程进行了详尽推导, 并将其与交错网格有限差分相结合, 进行波场数值模拟。对比分析波场模拟结果可知, CPML具有更好的吸收效果并且能够解决以往吸收边界在掠入射情况下吸收效果不理想的问题。最后用3个模型对本文方法进行了测试。

1 一阶速度-应力声波方程CPML的推导

当地震波以掠入射进入到经典PML吸收层时, 并没有传播至吸收层较深部位置, 而是沿着平行于吸收层传播, CPML的思想是在经典PML吸收边界的基础上, 引入更加广谱的滤波器将这部分能量衰减并且阻止其离开吸收层。由于交错网格在不增加计算量情况下能够压制数值频散, 为了便于将交错网格与CPML吸收边界相结合, 本文利用一阶速度应力声波方程进行推导。在直角坐标系中, 二维一阶速度-应力方程的时域形式^[11]如下:

$ \frac{{\partial v}}{{\partial t}} = \frac{1}{\rho }\frac{{\partial p}}{{\partial x}} $

(1)

$ \frac{{\partial w}}{{\partial t}} = \frac{1}{\rho }\frac{{\partial p}}{{\partial z}} $

(2)

$ \frac{{\partial p}}{{\partial t}} = k\left( {\frac{{\partial v}}{{\partial x}} + \frac{{\partial w}}{{\partial z}}} \right) $

(3)

式中:v和w分别为质点在x和z方向上介质震动速度; p代表声压; t为时间; ρ为介质密度; k为介质的体模量。

为了在波动方程中引入CPML吸收边界条件, 将声波方程进行傅里叶变换得:

$ - {\rm{i}}\omega \hat v = \frac{1}{\rho }\frac{{\partial \hat p}}{{\partial x}} $

(4)

$ - {\rm{i}}\omega \hat w = \frac{1}{\rho }\frac{{\partial \hat p}}{{\partial z}} $

(5)

$ - {\rm{i}}\omega \hat p = k\left( {\frac{{\partial \hat v}}{{\partial x}} + \frac{{\partial \hat w}}{{\partial z}}} \right) $

(6)

式中: $\hat{\upsilon },\hat{w},\hat{p}$ 分别是v, w, p的傅里叶变换; i为虚数单位; ω为频率。以x方向为例, 对(4)式进行频率域复坐标变换^[24], 经变换后的 ${\tilde{x}}$ 形式为:

$ \tilde x = {\kappa _x}x + \frac{1}{{{\alpha _x} + {\rm{i}}{\omega _0}}}\int\limits_0^x {{d_x}\left( s \right){\rm{d}}s} $

(7)

式中:d_x(s)是在x方向的衰减系数; 衰减程度依赖于传播距离; κ_x主要影响以掠射形式入射到衰减层地震波的衰减; α_x主要影响入射到衰减层的低频成分。如果κ_x=1, α_x=0, CPML边界条件就退化为经典PML边界条件。

将(4)式中的x替换为 ${\tilde{x}}$ , 则波动方程中涉及到对空间求导数的算子就转换为:

$ {\partial _{\tilde x}} = {\partial _x}\frac{{{\partial _x}}}{{{\partial _{\tilde x}}}} = \frac{1}{{{s_x}}}{\partial _x} $

(8)

式中:

$ {s_x} = {\kappa _x} + \frac{{{d_x}}}{{{\alpha _x} + {\rm{i}}\omega }} $

(9)

用 ${{{\tilde{s}}}_{x}}\left( t \right)$ 表示1/s_x的傅里叶反变换, 对方程(8)做时间维度的傅里叶反变换则有:

$ {\partial _{\tilde x}} = {{\tilde s}_x}\left( t \right) * {\partial _x} $

(10)

公式(9)有如下的形式:

$ \frac{1}{{{s_x}}} = \frac{1}{{{\kappa _x}}} - \frac{{{d_x}}}{{\kappa _x^2}}\frac{1}{{\left( {\frac{{{d_x}}}{{{\kappa _x}}} + {\alpha _x}} \right) + {\rm{i}}\omega }} $

(11)

因此,

$ {{\tilde s}_x}\left( t \right) = \frac{{\delta \left( t \right)}}{{{\kappa _x}}} - \frac{{{d_x}}}{{\kappa _x^2}}H\left( t \right){{\rm{e}}^{ - \left( {\frac{{{d_x}}}{{{\kappa _x}}} + {\alpha _x}} \right)t}} $

(12)

式中:δ是狄拉克函数; H是阶跃函数。令:

$ {\xi _x}\left( t \right) = - \frac{{{d_x}}}{{\kappa _x^2}}H\left( t \right){{\rm{e}}^{ - \left( {\frac{{{d_x}}}{{{\kappa _x}}} + {\alpha _x}} \right)t}} $

(13)

在时间域中公式(8)变为:

$ {\partial _{\tilde x}} = \frac{1}{{{\kappa _x}}}{\partial _x} + {\xi _x}\left( t \right) * {\partial _x} $

(14)

为了实现(14)式中褶积项的计算, 需要对时间进行离散化, 将时间长度离散成N个时间步长, 每个时间步长为Δt, 如果用ψ_xⁿ表示第n个时间步长褶积项的计算结果, 则有:

$ \psi _x^n = {\left[ {{\xi _x}\left( t \right) * {\partial _x}} \right]^n} = \int\limits_0^{n\Delta t} {{{\left( {{\partial _x}} \right)}^{n\Delta t - \tau }}{\xi _x}\left( \tau \right){\rm{d}}\tau } $

(15)

由于对时间的积分采用了交错网格, 所以对空间的导数可以定义在mΔt和(m+1)Δt中间的时刻上^[24], 因此(15)式变为:

$ \begin{array}{l} \psi _x^n = \sum\limits_{m = 0}^{n - 1} {\left\{ {{{\left( {{\partial _x}} \right)}^{\left[ {n - \left( {m + \frac{1}{2}} \right)} \right]\Delta t}}\int\limits_{m\Delta t}^{\left( {m + 1} \right)\Delta t} {{\xi _x}\left( \tau \right){\rm{d}}\tau } } \right\}} = \\ \sum\limits_{m = 0}^{n - 1} {\left\{ {{{\left( {{\partial _x}} \right)}^{\left[ {n - \left( {m + \frac{1}{2}} \right)} \right]\Delta t}}\left( { - \frac{{{d_x}}}{{\kappa _x^2}}} \right)\int\limits_{m\Delta t}^{\left( {m + 1} \right)\Delta t} {{{\rm{e}}^{ - \left( {\frac{{{d_x}}}{{{\kappa _x}}} + {\alpha _x}} \right)\tau }}{\rm{d}}\tau } } \right\}} = \\ \sum\limits_{m = 0}^{n - 1} {\left\{ {{{\left( {{\partial _x}} \right)}^{\left[ {n - \left( {m + \frac{1}{2}} \right)} \right]\Delta t}}\frac{{{d_x}}}{{{\kappa _x}\left( {{d_x} + {\kappa _x}{\alpha _x}} \right)}} \cdot } \right.} \\ \left. {\left[ {{{\rm{e}}^{ - \left( {\frac{{{d_x}}}{{{\kappa _x}}} + {\alpha _x}} \right)\Delta t}} - 1} \right]{{\rm{e}}^{ - \left( {\frac{{{d_x}}}{{{\kappa _x}}} + {\alpha _x}} \right)m\Delta t}}} \right\} = \\ \sum\limits_{m = 0}^{n - 1} {\left\{ {{a_x}{{\left( {{\partial _x}} \right)}^{\left[ {n - \left( {m + \frac{1}{2}} \right)} \right]\Delta t}}{{\rm{e}}^{ - \left( {\frac{{{d_x}}}{{{\kappa _x}}} + {\alpha _x}} \right)m\Delta t}}} \right\}} \end{array} $

(16)

式中:

$ {a_x} = \frac{{{d_x}}}{{{\kappa _x}\left( {{d_x} + {\kappa _x}{\alpha _x}} \right)}}\left( {{b_x} - 1} \right) $

$ {b_x} = {{\rm{e}}^{ - \left( {\frac{{{d_x}}}{{{\kappa _x}}} + {\alpha _x}} \right)\Delta t}} $

对(16)式进一步推导:

$ \begin{array}{l} \psi _x^n = {a_x}{\left( {{\partial _x}} \right)^{\left( {n - \frac{1}{2}} \right)\Delta t}} + \sum\limits_{m = 1}^{n - 1} {\left\{ {{a_x}{{\left( {{\partial _x}} \right)}^{\left[ {n - \left( {m + \frac{1}{2}} \right)} \right]\Delta t}} \cdot } \right.} \\ \;\;\;\left. {{{\rm{e}}^{ - \left( {\frac{{{d_x}}}{{{\kappa _x}}} + {\alpha _x}} \right)m\Delta t}}} \right\} = {a_x}{\left( {{\partial _x}} \right)^{\left( {n - \frac{1}{2}} \right)\Delta t}} + \\ \;\;\;{b_x}\sum\limits_{m = 1}^{n - 1} {\left\{ {{a_x}{{\left( {{\partial _x}} \right)}^{\left[ {n - \left( {m + \frac{1}{2}} \right)} \right]\Delta t}}{{\rm{e}}^{ - \left( {\frac{{{d_x}}}{{{\kappa _x}}} + {\alpha _x}} \right)\left( {m - 1} \right)\Delta t}}} \right\}} = \\ \;\;\;{a_x}{\left( {{\partial _x}} \right)^{\left( {n - \frac{1}{2}} \right)\Delta t}} + {b_x}\psi _x^{n - 1} \end{array} $

(17)

由(17)式可知, (14)式中的褶积项可以写成一个递推公式, 如果采用(16)式编程实现, 将会有繁重冗余的褶积计算, 增加了编程难度并且占用较大内存, 而采用(17)式编程实现, 只需每时刻更新记忆变量ψ_xⁿ, 这样就为编程实现带来了便利。从上述的推导可以看出, 要将CPML吸收边界条件引入到波动方程中只需将空间导数∂_x替换成∂_x/κ_x+ξ_x(t)*∂_x即可。将一阶速度-应力方程改写为具有CPML吸收边界条件的形式如下:

$ \frac{{\partial v}}{{\partial t}} = \frac{1}{\rho }\left( {\frac{1}{{{\kappa _x}}}\frac{{\partial p}}{{\partial x}} + {\psi _x}p} \right) $

(18)

$ \frac{{\partial w}}{{\partial t}} = \frac{1}{\rho }\left( {\frac{1}{{{\kappa _z}}}\frac{{\partial p}}{{\partial z}} + {\psi _z}p} \right) $

(19)

$ \frac{{\partial p}}{{\partial t}} = k\left( {\frac{1}{{{\kappa _x}}}\frac{{\partial v}}{{\partial x}} + {\psi _x}v + \frac{1}{{{\kappa _z}}}\frac{{\partial w}}{{\partial z}} + {\psi _z}w} \right) $

(20)

本文中α_x, d_x, κ_x的形式分别为:

$ {d_x} = - \frac{{3{v_{\max }}}}{{2L}}\ln R{\left( {\frac{x}{L}} \right)^2}\;\;\;\;0 < x < L $

(21)

$ {\kappa _x} = 1 + {\left( {\frac{x}{L}} \right)^2}\left( {\kappa _x^{\max } - 1} \right)\;\;\;0 < x < L $

(22)

$ {\alpha _x} = {\alpha ^{\max }}{\left( {1 + \frac{x}{L}} \right)^2}\;\;\;0 < x < L $

(23)

式中:x为CPML内的计算点到CPML吸收层内边界的距离; v_max为声波的最大速度; L是吸收层的厚度; R为理论反射系数; κ_x^max可以取1, 7, 20;α^max=πf₀, f₀为地震波的主频。

2 基于CPML的交错网格差分格式

求解带有CPML吸收边界的波动方程, 需要将其离散化, 交错网格较常规网格更能够有效压制频散^[25], 在不增加计算量的情况下可以提高数值模拟精度, 因此本文采用交错网格有限差分的方法对波场进行正演计算。波场及相关参数分布在相差半网格点的主体网格和交错网格点上, 以上定义能够满足递推要求。网格布局如图 1所示, 网格中介质波场及参数见表 1。

方程(18)至方程(20)的2N阶交错网格差分格式为:

$ \begin{array}{l} v_{i,j + 0.5}^{t + 0.5} = v_{i,j + 0.5}^{t - 0.5} + \frac{{\Delta t}}{{{\rho _{i,j + 0.5}}}}\left\{ {\frac{1}{{{\kappa _x}\Delta x}}\sum\limits_{n = 0}^{N - 1} {{a_n}\left( {p_{i,j + n + 1}^t - } \right.} } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\left. {p_{i,j - n}^t} \right) + \psi _x^tp} \right\} \end{array} $

(24)

$ \begin{array}{l} w_{i + 0.5,j}^{t + 0.5} = w_{i + 0.5,j}^{t - 0.5} + \frac{{\Delta t}}{{{\rho _{i + 0.5,j}}}}\left\{ {\frac{1}{{{\kappa _z}\Delta z}}\sum\limits_{n = 0}^{N - 1} {{a_n}\left( {p_{i + n + 1,j}^t - } \right.} } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\left. {p_{i - n,j}^t} \right) + \psi _z^tp} \right\} \end{array} $

(25)

$ \begin{array}{l} p_{i,j}^{t + 1} = p_{i,j}^t + {k_{i,j}}\Delta t\left\{ {\frac{1}{{{\kappa _x}\Delta x}}\sum\limits_{n = 0}^{N - 1} {{a_n}\left( {p_{i,j + n + 0.5}^{t + 0.5} - } \right.} } \right.\\ \;\;\;\;\;\;\;\;\;\;\left. {p_{i,j - n + 0.5}^{t + 0.5}} \right) + \psi _x^{t + 0.5}v + \frac{1}{{{\kappa _z}\Delta z}}\sum\limits_{n = 0}^{N - 1} {{a_n}\left( {p_{i + n + 0.5,j}^{t + 0.5} - } \right.} \\ \left. {\;\;\;\;\;\;\;\;\;\;\left. {p_{i - n + 0.5,j}^{t + 0.5}} \right) + \psi _z^{t + 0.5}w} \right\} \end{array} $

(26)

式中:a_n为差分系数。

3 数据测试 3.1 测试一

为了验证CPML边界条件的有效性, 进行了各向同性介质二维声波波场数值模拟。数值模型计算采用时间二阶、空间十阶交错网格有限差分法, 模型的计算区域为水平与垂直方向均为200网格点, 水平和垂直方向的空间网格大小均为10m, 介质速度为2500m/s, 密度1000kg/m³, 震源位于(100, 100)网格点处, 采用主频为20Hz的雷克子波, 时间采样间隔为1ms, 利用公式(21)至公式(23)对衰减系数进行设置, κ_x取1, 正演波场快照如图 2所示。图 2a表示计算区域外有20层CPML吸收层的波场快照, 图上基本上观察不到来自边界的反射, 计算区域的波场更加纯净, 图 2b表示计算区域外有20层经典PML吸收层的波场快照, 虽然经典PML层吸收了大部分的能量, 但是图上还是能看到来自边界较弱能量的反射, 这些反射能量对计算区域波场形成了噪声干扰, 不利于波场的精确模拟, 图 2c为计算区域外有50层经典PML吸收层的波场快照, 可看出, 当地震波传播到边界时, 没有出现明显的能量反射现象。对比分析可知, 在相同衰减层数的情况下, CPML吸收边界的吸收效果明显优于经典PML吸收边界; 当增加吸收层数时, 经典PML吸收边界也能够达到理想效果。因此, 为实现相同的吸收效果, CPML吸收边界条件所需吸收层数更少, 从而降低了计算所需的内存空间。

3.2 测试二

逆时偏移成像过程中, 吸收边界条件的选取直接关系到波场模拟的精度, 最终影响到成像剖面的刻画。为了说明在掠入射情况下, CPML吸收边界较经典PML具有更好的吸收效果, 特设计均匀介质模型如图 3所示。计算区域为水平方向440网格点, 垂直方向240网格点, 震源s位于(220, 25)网格点处, 在(60, 30)网格点布置检波器1, 在(220, 25)网格点布置检波器2, 在计算区域外均采用30层的吸收层, 其它参数同测试一, 正演模拟波场快照如图 4所示。

从图 4a可以看到, 地震波与吸收层入射角较小时, 经典PML吸收边界对入射地震波的吸收效果较好, 但是对掠入射地震波的吸收能力尚有不足, 有来自边界的反射现象, 这对有效信号形成噪声干扰, 不利于波场的精确模拟; 而图 4b采用CPML吸收边界, 不但具有更高的吸收能力, 更重要的是解决了地震波与吸收层形成掠入射情况下的吸收问题, 计算区域波场较为纯净。这是由于经典PML吸收边界条件的衰减系数依赖于入射波的方向, 当入射地震波与吸收层入射角较小时吸收系数有较大值, 当入射角较大直至发生掠入射的情况, 吸收系数有极小值, 因此能看到图 4a箭头位置的边界反射现象, 而在CPML吸收边界条件中引入的是复坐标变换, 一方面可以将入射波校正到垂直入射的情况; 另一方面可以削弱低频反射现象, 从而提高了正演模拟精度。

为了更加清晰地观察吸收边界条件对波场正演模拟精度的影响, 采集图 3所示模型中检波器1和检波器2的单道信号, 并与理论解进行比较。本文所述的理论解为在波场正演模拟中, 对模拟区域进行扩展, 从而保证地震波场不含任何计算边界影响时所得到的用于参考的信号记录。图 5为分别采用CPML和经典PML吸收边界条件采集的检波器1单道信号, 由于检波器1距激发点较远, 且位置接近上边界, 而地震波入射到上边界属于掠入射。图中红色曲线为理论解, 黑色虚线为检波器记录的单道信号, 由图 5a可见, 采用CPML吸收边界条件的单道记录与理论单道记录吻合较好, 虽然属于掠入射情况, 但是CPML边界能对掠射波进行很好的吸收; 由图 5b可见, 经典PML吸收边界条件对掠射波吸收的效果并不理想, 正演波场出现误差。图 6为分别采用CMPL和经典PML吸收边界采集的检波器2单道记录, 两种单道记录与理论记录相吻合, 说明在地震波与吸收层成小角度入射情况下, 两种边界条件吸收效果均较理想。

3.3 测试三

为考察不同吸收边界对地震记录的影响, 设计水平层状介质模型, 第1层厚度620m, 速度为2300m/s; 第2层厚度460m, 速度为2600m/s; 第3层厚度420m, 速度为2800m/s。对该3层水平层状介质进行正演模拟, 计算区域为水平方向300网格点, 垂直方向150网格点, 震源位于(150, 2)网格点, 地震记录时长1.3s, 两种吸收边界的层数均为20层, 其它参数均相同, 地震记录如图 7所示。从图 7a箭头处可以明显观察到来自边界的反射, 而图 7b中只有经过2个界面反射的双曲同相轴。对比可知, 利用CPML吸收边界条件可以得到更加纯净的地震记录。

4 结论与讨论

本文详细推导了CPML吸收边界条件下的一阶速度-应力声波方程, 并将CPML技术与交错网格有限差分相结合, 对均匀介质进行了数值模拟。对比经典PML和CPML吸收边界条件吸收效果可知, 在相同吸收层数的情况下, CPML吸收边界条件吸收效果明显优于经典PML吸收边界条件, 并且CPML吸收边界能够解决低频反射和掠入射的问题, 从而提高地震波场正演模拟精度。若将CPML吸收边界应用于逆时偏移成像等技术中, 可以提升波场传播精度, 而且较少的吸收层就能够达到理想效果, 节约了计算存储。理论上, CPML吸收边界也可以应用到基于qP波动方程的各向异性介质正演模拟中, 相对于PML吸收边界, 它可以对qP波和qS波在远偏移距的情况下形成更好的吸收效果。