地震声波数值模拟中人工边界条件的差别与组合

2022-03-10韩复兴王若雯孙章庆高正辉

吉林大学学报（地球科学版） 2022年1期

韩复兴，王若雯，孙章庆，高正辉

吉林大学地球探测科学与技术学院，长春 130026

0 引言

在地震声波数值模拟计算过程中,选取的有限计算区域会产生边界反射，干扰正常波场计算，因此需要引入人工边界条件降低边界反射的影响。人工边界条件按方法原理可以分为吸收类边界条件和衰减类边界条件两大类。吸收类边界条件通过在边界处使用单程波动方程模拟地震波的传播，使入射波只向计算区域外传播而不会产生边界反射；基于傍轴近似原理的Clayton-Engquist(CE)边界条件是吸收类边界条件中最经典的一种边界条件。衰减类边界条件的方法原理是在计算区域外引入一层衰减区域，入射波在衰减层内传播时，入射波能量逐渐衰减从而不会产生边界反射；其中完全匹配层(PML)边界条件应用效果最优，近年来得到广泛应用[1]。

近年来许多相关专家学者都对CE边界条件和PML边界条件的应用效果进行了比较，包括定性比较和定量比较。定性比较主要通过观察波场快照或合成地震记录的边界反射细节来比较边界条件的吸收效果。如：王守东[2]导出了完全匹配层法的声波方程表达形式，并通过增强地震记录振幅的显示方式得出PML边界条件吸收效果优于CE边界条件的结论；邢丽[3]通过数值算例分析认为，CE边界条件在程序方面易于实现，但1阶CE边界条件精度低，2阶CE边界条件有时又会出现不稳定的现象，PML边界条件虽然程序实现较为繁琐，但计算过程稳定，选取适当的参数时，几乎可以吸收全部的边界反射；王永刚等[4]分别采用井间和地面地质模型进行数值模拟，结果表明PML边界条件相对于旁轴近似法具有优越性；裴俊勇等[5]利用CE边界条件和PML边界条件进行有限差分正演模拟，模拟图像观察结果说明PML边界条件吸收边界效果较好。

定量比较通过量化边界反射数值进行，主要方法为根据反射波峰值高低比较边界条件吸收效果强弱。如，付小波等[6]对CE边界条件和PML边界条件进行了精细的定量比较性研究，发现CE边界条件相对PML边界条件受入射波角度限制更大，当入射角增大时，CE边界条件吸收效果大大降低，PML边界条件则不受入射波角度的限制。

本文通过对PML边界条件和CE边界条件进行详细的分析研究，考虑到其各自的优缺点，将CE边界条件和PML边界条件组成新的组合边界，并以数值模拟结果验证算法的有效性，以期达到在减少衰减带厚度的同时提高计算效率的目的。

1 PML边界条件两种差分形式对比

PML衰减介质中的波动方程可以看作是常规波动方程的推广，入射波传播到介质中振幅呈指数衰减。当衰减因子和衰减带厚度选择适当时，理论上入射波不会再反射回计算区域[7]。

PML的实现形式有分裂式(SPML)和非分裂式(NPML)两种[8-9]。SPML是原始的PML形式，其是在分裂原始波场分量的过程中加入PML衰减因子得出基于PML边界的波动方程。在计算区域和衰减区域都使用基于PML边界的波动方程，将计算区域的PML衰减因子设置为0，在衰减区域设置衰减因子，这种方式称为全局SPML。全局SPML虽然数据存储量增大，但相对而言编程实现简单[10]，也是本文考虑和使用的PML形式。在计算区域使用加入PML衰减因子前的普通波动方程、在衰减区域使用基于PML边界波动方程的SPML形式则为局部SPML。这种形式节省了大量的数据存储量，但要考虑不同区域的PML波动方程，编程难度也有所增加。NPML不需要分裂波场，有一定的优势，但在计算过程中要进行复杂的卷积计算，计算效率低[11-12]。

本文将讨论在地震声波有限差分数值模拟过程中实现全局SPML的非解耦与解耦差分形式的基本原理和实现效果。

1.1 全局SPML非解耦差分形式

均匀介质中二维声波波动方程为

(1)

式中：u为声波波场；(x,z)为空间坐标；v为声波波速；t为时间。孙林洁等[13]推导出基于PML边界的三维声波波动方程，根据其推导过程可以得到基于PML边界的二维声波波动方程[14]：

(2)

式中：d(x,z)为任意点(x,z)的衰减因子。对式(2)进行时间2阶、空间2m阶的有限差分。设有限差分空间采样步长为Δx、Δz，时间采样步长为Δt，x=iΔx、z=jΔz，t=kΔt，可以得到

(3)

1.2 全局SPML解耦差分形式

(4)

(5)

式中：ux、uz分别为x、z方向的波场值；dx、dz分别为x、z方向的衰减因子。

将式(4)(5)作向前差分，可以得到有限差分格式：

(6)

(7)

(8)

(9)

(10)

1.3 解耦与非解耦PML差分形式数值模拟对比分析

设定Δx=Δz=6 m，Δt=0.8 ms，总采样时间0.8 s，这里取衰减带厚度180 m。主频为30 Hz，雷克子波为激发震源。建立大小为1800 m×1800 m均匀介质模型(图1a)、层状介质模型(图1b)和复杂介质模型(图1c)进行解耦与非解耦两种差分形式的数值模拟，震源坐标分别为(900 m, 900 m)、(900 m, 240 m)、(900 m, 240 m)。

图1 均匀介质(a)、层状介质(b)、复杂介质(c)模型

均匀介质模型震源位于模型中心，层状介质和复杂介质模型震源位于模型上层，为更明确地观察波场的传播状态，均匀介质模型波场快照取0.56 s时刻，层状模型和复杂模型取0.72 s时刻。图2、图3、图4分别为均匀介质模型、层状介质模型和复杂介质模型使用时间2阶、空间2阶或4阶有限差分格式，非解耦及解耦PML差分形式进行数值模拟得到的波场快照。

由图2可以看出，当有限差分阶数为2阶时，非解耦与解耦PML差分形式数值模拟皆出现了少许的频散现象(图2a、b)，使用解耦差分形式时的频散现象更弱；当空间阶数升级为4阶时，两种PML差分形式都得到了较好的数值模拟结果(图2c、d)。

a.非解耦，时间2阶、空间2阶；b.解耦，时间2阶、空间2阶；c.非解耦，时间2阶、空间4阶；d.解耦，时间2阶、空间4阶。

从图3、图4可以看出，当地质模型变得更加复杂时，使用精度低的差分格式出现的频散会更加严重，解耦算法相对非解耦算法而言数值模拟效果要好一些(图3a、b，图4a、b)；空间阶数提高，解耦与非解耦算法都可以改善频散现象(图3c、d，图4c、d)。

a.非解耦，时间2阶、空间2阶；b.解耦，时间2阶、空间2阶；c.非解耦，时间2阶、空间4阶；d.解耦，时间2阶、空间4阶。

在现实生产试验中，我们通常使用空间8阶或更高空间阶数的差分格式进行计算，因此选择计算和编程更加方便的非解耦PML差分形式进行计算更具优势。

表1为使用3种模型(图1a、b、c)时，解耦与非解耦PML差分形式不同差分精度数值模拟的计算时间比较。可以看出，相比解耦PML差分形式，非解耦PML差分形式所用计算时间更少。这是因为解耦PML差分形式比非解耦PML差分形式增加了两项中间量的差分运算(对Ax、Az的差分计算)。同时根据计算原理可以得知，非解耦PML差分形式在计算上较为简单，在使用高精度差分格式时，计算和编程都比较容易，计算效率高、易于实现；而解耦PML差分形式在计算和编程上都较为繁琐。

表1 PML差分形式计算时间比较

2 CE边界条件与PML边界条件组合

2005年，吴国忱等[15]提出了一种将吸收类边界条件和衰减类边界条件进行组合的组合边界思想，即在衰减类边界条件衰减区域的最外层使用吸收类边界条件对边界反射进行再次吸收，数值模拟结果显示这种组合方式有较好的吸收效果及较少的计算量。2009年，杜启振等[16]运用这一组合思想将改进的扩边衰减边界与特征分析法组合，组合边界减少了衰减边界中衰减带的厚度，数值模拟结果显示这种组合方式有较好的吸收效果。本文采用这种组合思想，使用更具优势的非解耦PML差分形式，先使用PML衰减边界衰减到达边界的入射波，再在衰减边界外区域设置占用内存更小、方法更简单的2阶CE吸收边界吸收未衰减完全的入射波。因为使用2阶CE边界条件进行二次吸收，在PML衰减带内的入射波可以不衰减为0，所以组合边界条件使用较小的衰减带厚度。

分别使用均匀介质模型(图1a)、层状介质模型(图1b)以及复杂介质模型(图1c)验证所提出组合边界条件的有效性。采用时间2阶、空间12阶的高阶有限差分进行声波波场数值模拟。为说明组合边界的衰减效果，分别给出衰减带厚度为90 m时，以及增加衰减带厚度为180 m后的波场模拟快照。图5为均匀介质模型组合边界条件和PML边界条件在0.56 s时刻的波场模拟快照。图6、图7分别为层状介质模型和复杂介质模型组合边界条件和PML边界条件在0.72 s时刻的波场模拟快照。可以看出，在使用较小的衰减带厚度时，组合边界条件可以较好地吸收边界反射(图5a、图6a、图7a)。对比图5a与图5b、图6a与图6b、图7a与图7b可以看出，在采用相同衰减带厚度时，组合边界条件可以取得较好的吸收效果；对比图5a与图5c、图6a与图6c、图7a与图7c可以看出，组合边界条件在采用较少的衰减带厚度时也能达到很好的吸收效果。

a.组合边界，衰减带厚度90 m；b.PML边界，衰减带厚度90 m；c.PML边界，衰减带厚度180 m。

组合边界条件和PML边界条件吸收边界反射所用时间如表2所示。从表2中可以看出：3种模型组合边界所用时间与同样衰减带厚度下PML边界所用时间相差不大，这表明组合边界中所使用的2阶CE边界条件所占内存很小；而组合边界所用时间小于增加衰减带厚度后PML边界所用时间，说明组合边界条件的应用在保证吸收效果的同时可以有效提高计算效率，验证了这种组合边界的有效性。