杜 鑫
(中国石化石油工程地球物理有限公司华东分公司，江苏南京 211112)
全波形反演一词最早是由Pan等人在1986年提出的，全波形反演的核心思想是基于上世纪80年代初Tarantola[1]提出的利用正传波场与反传波场的广义最小二乘拟合而形成的。1988年，Tarantola[2]在已有的全波形反演基础上，给出了完整的时间域全波形反演方法的理论框架。随后，Bunks[3]等引入多尺度全波形反演思想，使反演结果更加稳定。目前，全波形反演已经从时间域发展到频率域并取得良好的效果[4-7]。
全波形反演目前应用较多的求解方法是收敛速度较快的梯度类局部优化迭代反演方法,主要包括最速下降法、共轭梯度法、拟牛顿法等[8-9]。Nocedal[10]提出了拟牛顿法(BFGS)，对收敛速度加速的同时能提高深层地层的分辨率，Pratt[11]等提出利用Hessian矩阵预处理梯度值以改善迭代的效果，张广智[12]等借鉴Beck[13]等的思想，对共轭梯度法的迭代算法进行改造，提出快速共轭梯度法并将这种方法应用到频率域全波形反演中。
采用林穗华[14-17]等提出的方法，通过增加谱因子的方式对共轭梯度法进行改造，结合Hestenes和Stiefe提出的算法与Polak、Ribieren和Polyakn提出的算法优点形成全新的修正谱共轭梯度法，并用这种方法来改善声波全波形反演效果。Marmousi模型的测试结果表明，该方法能在改善收敛效果的基础上加快迭代收敛速度。
全波形反演既可以在时间域实现也可以在频率域或拉普拉斯域实现，具有在复杂地质条件下精确刻画构造细节及岩性变化的潜力。
1.1 时间域全波形反演基本框架
全波形反演是一种高精度速度建模方法，全波形反演包含初始模型的建立、地震波场正演数值模拟、目标函数的建立、优化方法的选择等。全波形反演基本步骤为：①给定初始模型，并对其做正演模拟，设定迭代停止的精度要求及最大迭代次数；
②利用正演模拟合成的地震记录与实际观测的地震记录建立目标函数；
③判定目标函数反演误差值是否满足终止条件或迭代次数达到最大设定值，满足条件则输出该模型并作为最终反演结果；
若不满足则转至④；
④利用优化算法对目标函数的最小二乘结果进行迭代求解，利用梯度类局部优化算法进行求解，更新速度场；
⑤利用模型更新方程和速度模型，转至①。
1.2 时间域声波方程正演计算
时间域各向同性介质声波方程为：
(1)
在使用基于交错网格的有限差分算法进行数值模拟时，需要将声波方程转化为一阶声波方程：
(2)
式中：v为纵波速度，m/s；
vx,vz分别为x和y方向速度分量，m/s；
P为地层压力分量，Pa；
κ=ρv2为体积模量在空间的分布，Pa；
t为传播时间，s；
x为平行于地表方向的距离，m；
z为垂直与于地表的方向距离，m。
在有限差分地震波数值模拟过程中，需要引入具有吸收衰减属性的边界进行数值截断，通常采用完全匹配层(PML)吸收边界条件。其基本思想是将地震波场分裂为两组:一组垂直于波场传播方向，另一组平行于波场传播方向。对到达PML边界的波场，将平行界面法向传播的平面波进行迅速衰减，而对于其他方向传播的波并不进行衰减。在进行数值模拟时，衰减因子采用Collino推导的衰减公式如下：
(3)
式中：vmax为最大纵波速度，m/s；
x为离PML内层边界计算点的距离，m；
δ为PML匹配层厚度，m；
R为介质的反射系数。
1.3 时间域全波形反演
类似于其他反演问题，全波形反演的最小二乘意义下的目标函数可写为：
(4)
式中：dsyn为正演模拟得到的波场，dobs为实际地震数据的炮集波场。
通过迭代的方法更新初始速度模型，即：
vk+1=vk+1+αkdk
(5)
式中：α为迭代步长(可以通过线性搜索、抛物线插值等最优化方法取得)，k为迭代次数，dk为搜索方向。
全波形反演多次迭代极小化实际观测数据与合成地震数据的残差即目标函数，是实现反演速度场v的过程。目标函数的极小化以初始速度v0为前进方向进行搜索，通过多次迭代使目标函数在v0附近收敛到局部极小值。
共轭梯度法的搜索方向表达式的构成为当前次目标函数的梯度与前一次搜索方向的线性组合，具体公式如下：
(6)
式中：βk为共轭梯度修正因子；
gk为梯度方向。
为了提高共轭梯度法的计算效率，2001年Birgin和Martinez[18]结合谱梯度方法和CS，提出了SCG，其搜索方向如下：
(7)
式中：θk为谱系数。
当时θk=1时，式(7)转换为式(6)，实现谱共轭梯度法向对应的共轭梯度法的转换，因此共轭梯度法的优点自然地出现在谱共轭梯度法上。谱共轭梯度法适合求复杂的无约束优化问题，能够得到更好的数值计算结果。
林穗华提出的谱共轭梯度法如下：
(8)
利用谱梯度法进行全波形反演，需要计算梯度方向谱因子，还需要确定迭代步长ak。常用的步长求取方法有线搜索方法和抛物线搜索方法。本文采用的是改进弱Wolfe非精确线搜索方法。其具体条件为：
(9)
(10)
利用线搜索计算步长时，在给定搜索区间[0,αmax]求取能够使式(9)和(10)成立的最大步长αk。其中式(9)的目的是要求函数下降值大于沿切线下降。式(10)可以防止步长选取过小，确保函数值有足够的下降量。
共轭梯度算子取值为谱PRP算法以及谱HS算法的最小值，该方法结合了PRP共轭梯度法，能较好地克服小步长时收敛缓慢以及HS共轭梯度法不依赖于步长的优点，都满足共扼关系，具有充分的下降性。具体流程如下：①给定初始速度模型v1,终止条件ε>0，k=1，令初始方向d1=-g1；
②判别gk与ε的关系，决定是否终止计算；
③利用线搜索方法计算初始步长αk；
④由式(8)计算θk，βk,再由式(7)计算搜索方向dk+1(④对算法的效果影响最大)；
⑤令k=k+1,计算vk+1。转回②。
为了测试谱共轭梯度法全波形反演效果，利用Marmousi模型进行测试。不同算法选用相同的反演参数，雷克子波主频为5 Hz，空间采样间隔为dx=dz=25 m，正演时间采样为2 ms，单道记录长度8 s，迭代40次。本次模型实验采用PML吸收边界条件，其横向与纵向的网格数均设置为50。采用地表放炮，共41炮，炮间距为250 m，共481个检波器，检波点道距为25 m。计算平台主要配置为：处理器i5 3570，8 GB内存，500 GB硬盘存储，Ubuntu18操作系统。图1是Marmousi模型和初始速度模型。算法目前只支持二维地震数据以及速度模型输入。
图2是Marmousi模型两种方法迭代反演的速度模型对比。由图1和图2可以看出，两种方法都能得到基本符合目标模型的反演结果。谱共轭梯度法(用时56 h)用时略少于共轭梯度法(用时60 h)，深层反演结果的准确性及分辨率更高，更接近实际速度模型(图2圈住部分)，验证了本文所提方法的优势。
图1 Marmousi模型(上)与初始速度模型(下)对比
图2 谱共轭梯度法反演结果(上)与共轭梯度法反演结果(下)对比
在进行全波形反演迭代时，初始搜索方向均为负梯度方向。将每次迭代误差值除以初始误差值得到归一化误差，通过对比CG与SCG方法的归一化误差曲线(图3)，可见SCG方法在22次迭代时的误差与传统CG方法在35次迭代时的误差相近，SCG方法收敛速度更快且收敛过程相对稳定。
图3 两种方法反演误差曲线对比
相比于CG方法，SCG方法的额外计算量为谱因子θk的计算，其计算公式中的共轭因子βk已知，而gkTdk-1及范数求解在共轭因子计算的过程中可以得到。因此，谱因子计算过程中的额外计算量主要为数乘运算，与整个全波形反演计算过程中的矩阵乘法相比完全可以忽略。
(1)谱共轭梯度法具有共轭梯度法的全部优点，得益于谱系数的引入，利用了更多的目标函数信息，相对于共轭梯度法有更好的数值表现。
(2)Marmousi模型的反演计算结果在提高速度模型反演结果准确性的基础上能够加快迭代收敛速度。
猜你喜欢共轭步长梯度一个带重启步的改进PRP型谱共轭梯度法数学物理学报(2022年1期)2022-03-16基于Armijo搜索步长的BFGS与DFP拟牛顿法的比较研究成都信息工程大学学报(2021年5期)2021-12-30一个改进的WYL型三项共轭梯度法数学物理学报(2021年6期)2021-12-21基于随机森林回归的智能手机用步长估计模型中国惯性技术学报(2020年2期)2020-07-24巧用共轭妙解题河北理科教学研究(2020年1期)2020-07-24一种自适应Dai-Liao共轭梯度法应用数学(2020年2期)2020-06-24基于Armijo搜索步长的几种共轭梯度法的分析对比成都信息工程大学学报(2019年2期)2019-08-28一个具梯度项的p-Laplace 方程弱解的存在性华东师范大学学报（自然科学版）(2019年3期)2019-06-24基于动态步长的无人机三维实时航迹规划北京航空航天大学学报(2016年12期)2016-02-27