非分裂完全匹配层边界存储时间域全波形反演

成景旺 毛宁波* 吕晓春 常锁亮 严 皓 仲 华

(①油气资源与勘探技术教育部重点实验室(长江大学)，湖北武汉 430100； ②华北水利水电大学资源与环境学院，河南郑州 450011； ③太原理工大学矿业工程学院，山西太原 030024； ④中海石油(中国)天津分公司，天津 300452)

1 引言

全波形反演集地震子波估计、初始模型建立、正演模拟、反演于一体，是一套完整的理论体系，已被证明是一种建立高精度速度模型的有效方法[1]，可在时间域、频率域、Laplace域实现[2]。Tarantola[3]利用伴随状态法高效率地求取梯度，实现了二维时间域全波形反演；Pratt[4]将全波形反演理论推广到频率域；Shin等[5]针对地震数据带宽有限、初始模型获取困难等问题，提出利用阻尼波场零频分量反演低频模型作为频率域波形反演的初始模型，即Laplace域全波形反演。频率域全波形反演由于其固有的多尺度特征使其理论研究和实际应用得到快速发展[6-10]。然而频域率反演最大的问题就是大型稀疏方程组的存储以及求解过程中巨大的内存需求，大多借用MUMPS线性方程组求解软件包进行求解。但对于大规模尤其是三维情况下，频率域正演对内存的超大要求，一般的计算机集群仍无法满足要求，因此近年来时间域全波形反演研究成为热点[11，12]，尤其是时间域正演联合频率域反演的混合算法[13-15]。而制约时间域全波形反演的关键问题之一就是正演波场的存储或重建。由于全波形反演的一次梯度求取过程与逆时偏移过程相同，故可将逆时偏移的边界存储策略应用于全波形反演。国内外学者在解决逆时偏移巨大内存需求上已做了大量的研究，Symes[16]提出了采用设置检查点的方法以降低逆时偏移的存储量。Dussaud等[17]指出检查点技术虽然对存储量的需求最小，但是却明显地增加了计算量。Clapp[18]提出了边界存储策略，该策略要求存储边界网格层内所有时刻的波场和整个空间最后时刻的波场，反传时作为边界条件和初值条件，只需额外正演一次即可重建正演波场；Clapp[19]又提出了随机边界方法，只保存最后一个时刻所有空间点的波场值作为反传时的初值，存储量进一步减少但是在偏移剖面中会引入噪声。在所有这些策略中，Clapp的边界储存法具有存储量小、计算量相对较少且适用于任何边界条件等优点[20]，因此该边界存储法被广泛研究并被用于逆时偏移。

不同的边界条件的吸收效果以及存储量都会不同。完全匹配层(PML)边界条件被认为是最好的吸收边界条件，实现方法主要有全局分裂完全匹配层(SPML)边界条件、局部SPML边界和非分裂完全匹配层(NPML)边界条件。其中全局SPML对边界区域和计算区域使用同样形式的PML波动方程，编程简单，但是由于波场分成垂直边界和水平边界方向两部分(二维为2个方向，三维为3个方向)，因此在整个计算过程中所需存储量是常规波动方程的2倍甚至3倍；局部SPML边界条件只有在PML吸收层内采用PML波动方程，而在计算区域采用常规波动方程，因此只有在吸收层内需要额外增加存储量。但局部SPML在计算区域和边界区域使用不同形式的波动方程，因此实现起来较繁琐，除了计算区域外，需要考虑模型边界、边角的问题，编程较困难。目前局部SPML边界条件已被广泛用于基于边界储存的波场重建[20，21]。但是传统的SPML边界条件存在一定的缺陷，对大角度入射产生的掠射波吸收效果不佳，为此引入了复频移(CFS)伸展函数改进坐标变换，该边界条件称为复频移(CFS-PML)吸收边界条件[22-24]。采用复频移伸展函数后，PML不易采用传统的分裂形式实现，而采用不分裂卷积算法时需要进行大量卷积计算。为此，Komatitsch等[25]采用了不分裂递推卷积方法实现了CFS-NPML吸收边界，并进行了弹性波的数值模拟；Drossaret等[26]提出了基于递归积分的非分裂CFS-NPML边界条件，通过引入辅助变量实现递归积分，同样避免了卷积计算。这两种NPML边界条件虽实现方法不同，但最终计算使用的离散公式形式相同。综合考虑，CFS-NPML不需要对波场进行分裂，实现简单，边界处所需变量个数少，因此更加适合于基于边界存储的波场重建。本文分析了CFS-NPML条件的实现原理和吸收效果，提出基于CFS-NPML边界条件的有效边界存储策略进行正演波场重建，实现时间域全波形反演，并通过理论模型论证其可行性。

2 基于边界存储的波场重建方法

2.1 CFS-NPML边界条件

采用时间域一阶应力—速度方程进行全波形反演。在二维情况下，弹性介质一阶应力—速度波动方程可表示为

(1)

(2)

式中sp为伸展函数。CFS-NPML与SPML边界相比，主要的区别就是引入了复频移伸展(CFS)函数。CFS函数定义为

(3)

其中κp(p=x，z)和αp(p=x，z)为复频移伸展函数中的两个参数，满足κ≥1和α≥0。当κ=1且α=0时就变为常规的PML边界条件。参数κ主要用于吸收广角入射时内边界产生的瞬逝波，而参数α主要影响对波的低频成分的吸收。σp为伸展坐标系下的衰减因子。本文的参数计算公式[24]为

(4)

式中：L为PML边界网格层数；l表示与内边界的网格点距离；Pd、Pκ、Pα为多项式衰减因子的系数，取值范围一般为[1，4]，通常情况下取2；κmax通常介于1～20；αmax=πf0，f0是震源的主频；cp是PML内纵波的传播速度；Rc为边界的吸收系数，一般取0.0001。

引入两个辅助变量，经详细推导可得到时间域CFS-NMPL边界条件的一阶应力速度方程[22]为

(5)

式中Ω和Ψ为推导过程中引入的辅助变量，其对应的控制方程为

(6)

(7)

将式(6)和式(7)写成统一形式的一阶微分方程

(8)

则辅助变量可通过下式迭代求解[24]

(9)

式中： Δt为时间步长；n为时间采样序号。比较式(5)与式(1)可以看出，通过引入辅助变量，将波动方程右端分成了正常项和衰减项两部分。在计算区域内，式(5)演变为式(1)，可以按照正常波动方程进行计算。而在边界吸收层内可先计算正常项，然后减去由辅助变量表示的衰减项即可。因此CFS-NPML边界条件不需要对波场进行分裂，只需在边界吸收层内额外增加几个辅助变量的存储量即可。对于SPML边界条件，需要对波场分量分为x方向和z方向两部分(三维情况下为三部分)，如vx=vxx+vxz。对全局SMPL、局部SPML以及CFS-NPML在边界吸收层内需要申请的变量个数进行对比，如表1所示。其中全局SPML在二维情况下需要申请10个变量，三维情况下需要申请27个。但是由于全局SPML是整个模型保存变量，因此所需内存最大。局部SPML在二维情况下边界内需申请15个变量，三维情况下需申请36个变量；CFS-NMPL在二维情况下边界内需保存13个变量，三维情况下需申请27个变量。所以CFS-NPML边界条件不仅实现方便，而且在吸收层内所需的变量个数最少。

表1 不同PML吸收边界条件边界层内所需变量个数统计

注：i=1代表x方向；i=2代表z方向；i=3代表y方向

为了对比CFS-CPML边界与局部SPML边界吸收效果，建立一个网格数为300×300的均匀介质理论模型，纵波速度为3300m/s，横波速度为1905m/s，空间步长为5m，时间步长为0.5ms。震源采用主频为20Hz的Ricker子波，位于模型中心(750m，750m)，选择较少的边界网格点数进行边界吸收效果对比。当边界吸收网格点数为10时(图1)，局部SPML边界条件出现了较弱的边界反射，而CFS-NPML的边界吸收依然较好，即使能量放大到10倍也看不见边界反射，能量放大到100倍时可以看到微弱的反射，明显好于局部SPML边界条件。当边界吸收网格点数为5时(图2)，局部SPML和CFS-NPML均有边界反射，但CFS-NPML的边界反射能量较弱，即使能量放大10倍也要优于边界网格点数为10时的局部SPML条件。在(750m，300m)处取单道波形进行对比(图3)，未做边界吸收处理的反射波与直达波能量相当，其他四种方式的边界条件均对边界反射(0.35s后)有吸收作用。将该道记录中的边界反射局部放大，可以清楚地看到CFS-NMPL边界条件的反射能量均小于SPML边界条件。设E为未作边界处理的边界反射波最大振幅，E′为吸收边界处理后的边界反射波最大振幅，将比值(E-E′)/E作为分析边界反射好坏的标准。采用局部SPML(边界网格数为5)的边界反射吸收率为88%，采用局部SPML(边界网格为10)的边界反射吸收率为92.9%； 采用CFS-PML(边界网格数为5)的边界反射吸收率为98.1%，采用CFS-PML(边界网格数为10)的边界反射吸收率为99.9%。因此为了进一步减少边界存储内存，在尽量减少边界网格点数的前提下，CFS-NMPL比常规PML边界条件具有更好吸收效果。除此之外CFS-NPML边界条件不需要对波场进行分裂，计算效率高，编程难度小，是时间域全波形反演波场重建最有效的正演边界处理方法。

图1 边界网格点数为10时SPML(上)和CFS-NPML(下)边界条件吸收效果对比

2.2 基于边界存储的波场重建实现原理

基于边界储存的波场重建就是要利用最后一个时刻(二阶时间离散)的波场值通过反传计算出前面任意时刻的波场值。无论是哪种方式的PML边界条件，其在时间方向上都是不可逆的，因此需要存储边界区域每个时刻的全部波场值。反向传播作为边界条件替换边界处的波场值，此时只需计算内部区域，不需再做边界处理就可以完全计算出前面任意时刻的波场值。

基于边界存储的波场重建可按如下步骤进行：

图2 边界网格点数为5时SPML和CFS-NPML边界条件吸收效果对比

图3 不同边界条件(750m，300m)处的地震记录

(7)若n>1，则返回步骤(2)；n=1结束，完成波场重建。

在上述波场重建过程中，应特别注意应力分量和速度分量的更新迭代次序一定要与正向传播时的相反，否则随着时间的迭代，在计算区域与边界区域交界处会产生新的绕射源，导致重建波场中存在干扰波。为了验证上述CFS-NPML边界存储波场重建的有效性，应用理论模型进行测试。将Marmousi模型上面加入一海水层，并等比例缩小。最终模型网格数为300×116，网格间距为10m，时间步长为0.5ms，震源采用位于(1500m，120m)处20Hz的Ricker子波，为了尽量减少边界保存所需内存，取PML吸收层网格点数为10。图4为正向传播和逆向重建的200ms、450ms和800ms波场快照。图4左边一列为正向传播波场，中间一列为速度分量和应力分量的更新次序与正向传播一样时得到的重建波场，最右边一列速度分量与应力分量更新次序与正向传播相反的重建波场。可以看出，利用上述重建步骤得到的重建波场(右列)与正向传播波场完全一致，说明利用该边界存储策略完全能够重建正演波场。而当速度分量和应力分量更新次序与正向传播一致时，可以看出得到的重建波场大致与正向传播波场相同。800ms时刻的重建波场与正向传播波场一致，但随着波场由边界向里边传播，在计算区域与边界交界处出现了新的干扰波场，并随着正常波场一起向模型内部传播。从200ms和450ms的波场快照图中可明显看出存在干扰波，这说明在波场重建过程中速度分量和应力分量的更新次序对重建波场的重要性。

图4 200(上)、450(中)、800ms(下)正向传播波场与重建波场快照对比

3 全波形反演模型试算

为验证本文方法的正确性，采用Marmousi模型进行试算。基于一阶应力—速度弹性波波动方程，采用预条件共轭梯度法进行FWI[27]。为了减少反演时间，将Marmousi模型等比例缩小，并在模型上面加入一海水层，最终模型网格数为300×116。反演初始模型(图5)为真实模型的二维高斯函数光滑结果，光滑过程中纵横向的相关长度均取100m。空间网格间距为10m，时间步长为0.5ms，记录长度为2s。震源采用主频为20Hz的Ricker子波，CFS-NMPL边界网格数取为10。

震源点和接收点深度均设置为140m。一共激发19炮，炮间距为150m，每一炮都有固定的275道接收，道间距为10m。由于时间域FWI是全频带数据反演，反演容易陷入局部极小值，因此采用滤波器进行多尺度反演[28]，从低到高依次给定截止频率实现不同频带的反演。本文采用巴特沃斯低通滤波器进行滤波，给定的低通截止频率分别为5Hz、9Hz、16.36Hz和29.2Hz。反演过程中每个频带设置的最低迭代次数为20次。 正演过程采用基于区域分解的并行计算[29]，本次数值测试一共采用10个进程进行计算。为了减弱多参数弹性波FWI中各个参数之间的相互耦合性，反演过程采用杨积忠等[30]提出的多参数反演策略。最终四个频带的纵、横波速度和密度反演结果如图6～图8所示，可以看出随着反演截止频率的增加，模型的细节构造越来越明显。为了进一步更加精确地对比反演结果，提取横向750m、1250m和2000m处的速度进行对比(图9～图11)，可以看出，反演得到的纵、横波速度、密度与真实相值吻合。

图5 Marmousi模型(左)及反演初始模型(右)

表2给出了常规保存波场FWI和本文方法所需的内存存储量和梯度计算所用时间。根据纵、横波速度和密度的反演梯度，需要保存vx、vz、τxz、τxx和τzz共5个变量，其中常规FWI需保存所有时刻所有网格点的波场值，而本文方法只需保存边界处(边界网格点数仅为10)的波场值即可(计算中波场采用单精度浮点数)。

表2 两种反演方法内存和梯度计算时间对比

从表2可以看出，基于非分裂PML边界存储的FWI在内存需求上明显减少，仅仅是常规FWI内存需求的23.9%。由于常规保存波场FWI在一次梯度求解中相当于进行两次正演模拟，而基于边界存储法在一次梯度求取中相当于要进行三次正演模拟，因此基于边界存储的全波形反演必然会增加计算时间。但从表2中可以看出，由于采用并行计算技术，所有的梯度计算时间总和并没有显著增加，其计算时间比常规FWI增加了16%。本文提出的基于边界存储的FWI既可有效减少内存需求，又能保证计算效率。

图6 不同截止频率反演的纵波速度

图7 不同截止频率反演的横波速度

图8 不同截止频率反演的密度

图9 不同横向位置处的纵波速度真实模型(黑色)、初始模型(红色)和反演结果(蓝色)对比

图11 不同横向位置处的密度真实模型(黑色)、初始模型(红色)和反演结果(蓝色)对比

4 结论

本文将逆时偏移中的边界存储重建波场技术引入时间域全波形反演，并采用CFS-NPML边界条件代替传统PML边界条件解决了时间域全波形反演的巨大存储量问题。CFS-NPML边界条件不需要对波场进行分裂，在边界吸收层内所需的变量个数最少，且在边界网格数较少的情况下，CFS-NPML边界条件的吸收效果要优于传统PML边界条件。因此基于CFS-NPML边界条件实现波场重建，可最大化地减少边界存储内存需求，是时间域全波形反演波场重建的有效方法。同时，应用并行计算技术，可保证迭代过程中的梯度计算时间没有显著增加。

本文提出的基于边界存储的全波形反演实现方便简单，既有效减少了内存需求，又保证了计算效率。该边界条件还可用于黏弹性或各向异性等复杂介质的正演波场重建，也可直接推广到三维。基于CFS-NPML边界存储的复杂介质三维时间域全波形反演是进一步的研究方向。