一种基于改进共轭梯度法的弹性波全波形反演速度分层建模方法

2023-07-08王立德王小卫张志强王建乐王德英

岩性油气藏 2023年4期

王立德，王小卫，周辉，吴杰，张志强，王建乐，王德英，冯刚

（1.中国石油勘探开发研究院西北分院，兰州 730020；2.中国石油大学（北京），北京昌平 102249；3.中国石油大庆油田有限责任公司采气分公司，黑龙江大庆 163000）

0 引言

随着油气资源的开发利用，勘探开发成本低、类型相对简单的构造型油气藏大多已被发现和开发，而复杂构造区往往是一些成藏质量较差的区域，地震成像复杂，勘探难度大，因此高精度地震成像需求越来越迫切。弹性波全波形反演（EFWI）是利用多分量地震记录的多种波场信息进行岩性参数反演，相比其他的建模方法具有更高的精度［1］。以全波形反演为代表的地震反演成像技术在地震地质建模综合研究中也扮演着越来越重要的角色［2］，但是EFWI 存在很强的非线性，初始模型的建立、低频数据的缺失等问题均可能导致反演陷入局部极小值，从而制约了EFWI 在实际资料（尤其是陆上实际资料）中的推广应用。

全波形反演（FWI）需要直接求取目标函数对模型参数的导数从而进行速度更新，这将带来巨大的计算量和存储量。为克服大计算量问题，Lailly［3］提出了基于震源波场和残差波场的梯度求取方法，其具体做法是对2 种波场进行一次互相关。该方法对于每个单独激发的震源只要正演2 次就可进行1次更新，能够减少全波形反演的计算量，由于该方法是基于线性理论的，对模型的大尺度背景恢复起来很困难。Tarantola［4］基于时间域纵波方程将非线性原理应用到最小二乘反演中，并推导出了利用伴随方程法求取介质参数更新方向的公式。Tarantola［5］基于弹性介质中的波动方程形式提出并实现了利用弹性波进行全波形反演的思路，并将FWI 理论推广到了EFWI。任浩然等［6］利用多种地震信息进行建模，但多重因素的叠加加大了反演的非线性并降低了实用性。为降低EFWI 的非线性程度，许多学者做了针对性的研究，最常用的是自地震数据低频成分到高频成分逐渐匹配的多尺度反演策略，该策略可以在时间域［7］、频率域［8］和拉普拉斯域［9］实现，但后2种域的计算过程复杂，对计算资源的硬件性能要求更高。在时间域可以通过对地震数据的滤波实现多尺度反演，由于地震数据本身就是时空域的，对地震数据进行时间和空间域相关处理更方便。从小炮检距数据到大炮检距数据的反演，通过限制参与反演数据的炮检距来实现。大炮检距可以有效重建低波数信息，小炮检距则对高波数成分具有很好的恢复作用［10-12］。另外，比较常用的降低EFWI非线性程度的方法是利用初至波进行近地表速度反演，因为EFWI 的非线性大部分由反射波引起的，当初始速度精度不够时，反射点很难归位［13］，因此，避免反射波参与反演可有效地降低反演的非线性。尽管很多学者针对EFWI的非线性问题进行了大量研究，并在理论模型上得到了验证，但对于实际资料（尤其是陆上资料）而言，EFWI 的速度建模方法依然难以见到效果，这是因为陆上资料更容易受震源能量、激发响应、噪音污染、近地表传播方式等的影响［14］，从而加剧EFWI 的非线性程度。王华忠等［15］指出EFWI的非线性取决于介质模型的复杂性和描述地震波物理传播过程的正算子的复杂性，并在分析了陆上和海上地震数据特点的基础上，提出了全波形反演走向实用化的策略，即利用特征波场反演而不是全波场，同时增加相位信息，尽量充分考虑初始模型的先验信息。

EFWI 的非线性很大部分来自于深层反射波，这将直接影响纵、横波速度场的反演精度。为提高速度反演的精度，提出一种分层反演策略：①利用时窗选取包含浅层速度信息的较早时间到达的地震数据，用小偏移距早至波数据，自低频到高频多尺度重建浅层速度场，以减少中深层反射波及噪声带来的强非线性影响，提高浅层速度的反演质量。早至波与初至波一样，并不是特指某一类型的波，而是包含了很多波的信息。小偏移距早至波主要包括浅层直达波、反射波及小折射波等。筛选的具体标准是根据反演浅层界面的深度以及相应深度的初始速度场，估算出地震记录时间，这样选择的偏移距范围内的地震记录时间能够约束地震波穿透深度，从而获得包含浅层速度场的波场记录。在EFWI 中，浅层速度的反演精度将直接决定深层反演的好坏。在这个过程中用小尺度交错网格有限差分模拟浅层波场。②固定好浅层速度，用全波场、全偏移距数据对中深层速度进行从低频到高频的多尺度反演，用大尺度交错网格有限差分模拟中深层波场，在减少计算量的同时保证波场模拟的精度。③最后使用伴随波场法求取纵、横波速度梯度，引入深度学习优化器RMSprop［16］的梯度处理方式，提高反演的稳定性和收敛速率，进一步提高深层速度场的反演精度。将该分层反演策略应用到理论模型和陆上实际资料，以期有效提高浅层和中深层速度场的反演精度。

1 方法原理

1.1 时间域EFWI 基本原理

全波形反演的基本思想是用模拟数据不断地逼近观测获得的实际数据，使它们的误差最小化，此时拟合结果对应的模型参数最接近真实的模型参数。基于模拟和观测记录，全波形反演目标泛函往往可以被定义成多种不同的误差形式，不同形式的目标泛函采用不同的参考标准匹配模拟数据和观测数据。理论上如果目标泛函能够达到一定范围内的极小值，此时的模型参数就会更加接近实际情况。此次研究采用的EFWI 目标函数表达式如下：

式中：uobs(x,t)为观测记录；u(m;x,t)为正演模拟记录，在分层反演方法中，模拟记录按照分层时窗大小进行调整；dΩ 表示对检波点接收区域各道的积分求和；dt表示对记录时间的积分求和；x为地震道；m为模型参数，本文中指浅层及中深层的纵、横波速度；T表示地震记录的时间，s；Ω 为检波点接收区域。

直接求取目标函数关于速度的梯度非常困难，比较常用的方法是通过伴随波场求取的伴随方法。Zhang 等［17］从Lagrange 乘子法原理出发，通过定义褶积型目标函数，得出了基于弹性波波动方程的伴随方程，并给出了相应模型参数的梯度公式，本文采用此种方法求取模型参数的梯度。二维弹性波波动方程通过伴随方程法得到的纵、横波速度场的梯度表达式如下：

式中：ρ为密度，g/cm3，此次研究不涉及密度反演，因而为常数；vP，vS分别为纵、横波速度，m/s；λ为伴随波场值；v为震源波场值。

1.2 基于RMSProp 的共轭梯度优化方法

1.2.1 基本理论

EFWI 主要利用非线性最优化方法最小化观测记录和理论记录之间的误差，以重建地下介质的物性参数。为了最小化目标函数，通常使用梯度类优化算法来解决非线性问题。牛顿法具有二阶收敛的优势，相较于一阶收敛的最速下降法和介于一阶与二阶收敛的共轭梯度法，它的收敛更快，但是需要直接求取目标函数的海森矩阵，难度较大，且需要巨大的计算量和巨大的内存资源成本，因此该方法在EFWI 中很难推广应用。共轭梯度法（CG）的收敛速率介于一阶和二阶方法之间，可求解无约束二次规划问题，由于其表达式简单，实现起来容易，在大规模数据的无约束优化问题中具有独特的优势［18］。其主要逻辑是找到一组共轭正交的方向向量，对每个方向上的模型参数依次进行迭代更新，当每个方向上的误差达到或接近0 时，迭代结束。该方法的收敛速度大于最速下降法等一阶的优化方法。它的模型参数更新方式可以用下式表示：

式中：gn为公式（2）和（3）所计算的纵波和横波速度的梯度场；n为迭代次数；dn为通过共轭梯度法计算的更新方向，将其代入下式进行模型参数的更新

式中：α为更新步长。

RMSProp 来源于深度学习，本文将这种优化理论应用到共轭梯度法中，相比于传统共轭梯度方法，新方法可使EFWI 收敛更快、反演精度更高。优化后需对纵、横波速度场每个位置的梯度gn做如下处理：

式中：gn,i为纵、横波速度梯度gn的第i个元素；为对gn,i移动加权平均的结果；η为常数，一般取值为0.99；为处理后的纵、横波速度梯度元素，所有形成了新的梯度场。将代入到公式（4）和（5）便可进行纵、横波速度的迭代更新。相对于传统的共轭梯度法，优化后的方法具有如下优势：①可使当前参与更新的梯度包含先前若干次迭代时的梯度信息，可保证整个反演过程的稳定性并提高收敛效率。②公式（7）的作用相当于实现了搜索步长的自适应调整，当反演梯度gn,i非常大时，计算得到的也很大，将变小，代入公式（4）计算的dn,i相应变小，按照公式（5）的步长α进行速度更新时更新量亦相应减小，可避免发生因更新过量而越过目标泛函极小点且在极小点附近摆动的现象，从而提高了反演精度。反之，梯度gn,i非常小时，通过增大梯度可以加速目标函数的下降，使收敛更快。③在EFWI 中，深层的梯度能量很弱，导致深层区域纵、横波速度场的重建效果不理想。本文提出的公式（7）的处理方式隐含浅层和中深层速度场梯度能量的均衡化处理，使深层的梯度能量相对提升，提高了深层速度的反演质量。

1.2.2 数值算例

为验证基于RMSProp 改进的CG 方法对EFWI效果的提升，选择了二维Overthrust 模型进行反演测试。在反演参数相同的情况下，分别用传统的CG 法和改进后的CG 法重建该模型（图1）。

图1 Overthrust 模型的真实速度和初始速度Fig.1 True and initial velocity of Overthrust model

该二维Overthrust 模型共计800×170 个网格点，交错网格有限差分空间步长为8 m；使用雷克子波生成震源，模型设计50 炮，每炮800 个检波器；为降低非线性影响，选择多尺度反演，共设置4 个频带，各频带的主频依次是4 Hz，7 Hz，11 Hz 和15 Hz，每个频带迭代20 次。对比2 种方法对应的反演结果可以看出，无论是纵波速度还是横波速度，改进后的CG 方法重构的速度场精度更高，尤其在深层，构造更加清晰，速度值更接近真实值，说明本文方法能够有效平衡浅层和深层的梯度能量，从而改善深层反演的结果（图2）。

图2 不同方法反演的Overthrust 模型速度场对比Fig.2 Comparison of velocity field of Overthrust model inversed by different methods

为了更加直观地对比2 种方法的反演精度，分别求取了2 种方法在整个反演过程中的速度误差值，该误差是基于真实速度模型的相对误差。在4个频带共80 次的迭代过程中，改进CG 反演的速度误差曲线均位于传统方法误差曲线下方，误差曲线下降更快，说明优化后反演的速度场精度更高，收敛更快，通过更少的迭代次数即可达到传统CG 反演的精度（图3）。

图3 不同频带反演速度的相对误差Fig.3 Relative errors of inversion velocity in different frequency bands

1.3 分层反演策略

1.3.1 基本理论

EFWI 强非线性是由资料采集等客观因素引起的地震数据与反演参数之间的强不确定性，因此通过筛选不同数据子集，分步骤、多尺度的反演是降低EFWI 非线性程度的有效策略（图4）。基于这一规律的时间域EFWI 分层反演策略，可降低反演的非线性程度，提高反演质量。具体反演过程如下：

图4 速度场分层反演流程Fig.4 Flowsheet of layered inversion of velocity field

（1）利用反演浅层界面的深度以及相应深度的初始速度，估算出相应地震记录时间，用该时间窗口选取较早到达的地震数据，同时选定偏移距范围。用包含丰富浅层速度信息的小偏移距直达波、小折射波［19］、浅层反射波等早至波反演浅层速度，规避主要由中深层反射波带来的强非线性干扰。

（2）通过从低频到高频逐渐匹配的多尺度反演手段，重建浅层速度。反演时采用小网格剖分的交错网格有限差分方法进行波场模拟，以保证浅层速度波场模拟的精度，提高浅层反演质量。

（3）固定反演重建的浅层速度场，约束中深层反演，降低反演的不确定性，利用全采样时间的全偏移距数据反演中深层速度。

（4）由低频到高频多尺度反演中深层速度，由于中深层速度较大，采用大网格的交错网格有限差分模拟波场，利用含有中深层速度场信息的地震记录段反演中深层速度场。该方法在减少计算量的同时可获得高精度的中深层反演结果。

分层反演时速度场梯度需要按照以下公式进行处理

式中：d1表示浅层速度场的深度，m。梯度gn经过公式（9）定义的分段处理函数p(z) 处理后，浅层梯度值为0，不再更新，而深层梯度的能量得到提升。

1.3.2 数值算例

为验证分层反演策略的有效性，利用二维Marmousi 模型进行测试（图5）。速度更新时，对二维Marmousi 模型的初始速度场分别进行常规EFWI和分层反演策略下的EFWI。整个速度场的地震记录采样时长为4 s，采样间隔为1 ms，共40 炮激发，用0～250 ms 的时间窗筛选地震记录以确定偏移距范围，进行多尺度反演的雷克子波主频依次为4 Hz，8 Hz，12 Hz，18 Hz 和25 Hz。浅层深度为300 m，真实纵、横波速度场深度为0～300 m 时为低匀速层，有限差分网格为5 m×5 m。浅层反演结束后固定浅层速度，用整个4 s 的全偏移距地震记录进行多频带反演中深层速度场，有限差分网格为10 m×10 m。对比基于2 种反演策略得到的纵、横波速度场，可见浅层低匀速层反演结果在反演参数完全相同的情况下，无论是纵波速度还是横波速度，分层反演的低速层均更干净，更接近真实速度场，而常规反演在低匀速层存在干扰；同时，在深度为0～1.5 km，距离为0～6 km 的区域内，分层反演得到的纵、横波速度场均更加清晰，局部构造的分辨率更高，弱反射层得到了较好的恢复（图6）。

图5 二维Marmousi 模型的真实速度和初始速度Fig.5 True and initial velocity of 2D Marmousi model

2 实际资料应用

为了进一步测试新方法对实际资料的有效性和适应性，对四川某三维区块的二维测线进行EFWI速度分层建模测试。该三维区块为正交观测系统采集，最小炮检距为28.28 m，覆盖次数为11，道距为40 m，子波主频为25 Hz，三分量地震记录的采样时长为5 s，采样间隔为1 ms。X分量地震记录信噪比低，很难看到有效信息，Y分量信噪比更低，Z分量信噪比相对较高（图7）。通过层析反演得到该区块的初始速度，待反演的二维速度场横向距离为5 km，纵向深度为7 km。

分层反演时，给定纵、横波初始速度场（图8），采用0～1 s 时窗反演深度为0～2 500 m 的速度场。反演该深度时有限差分的空间采样间隔为15 m，反演更深的速度场时空间采样间隔为25 m。共选取300 炮近偏移距数据参与浅层速度的反演，反演设置了5个频带，采用主频依次为6 Hz，10 Hz，15 Hz，20 Hz 和25 Hz 的雷克子波进行波场模拟。从浅层速度场的反演结果（图9a，9b）来看，浅层速度场得到了有效更新。由于浅层速度已经被重建，与待反演的深层速度在交界位置产生了不耦合层，这将在下一步的反演中产生虚假反射影响深层速度建模，因此在反演深层速度场时，需要对该不耦合层进行光滑处理以降低虚假反射的影响。在平滑不耦合层后，固定浅层速度场，迭代更新深层速度，由于深度2 000 m 附近的速度场被重新平滑，因此反演深层速度时需要从这一深度开始反演。深层速度场反演时采用了300 炮的全偏移距全记录时长的数据，除空间采样间隔和所用地震数据范围外，所有的反演参数与浅层一致。从全速度场最终的反演结果（图9c，9d）可以看出，浅、中、深层速度场均被更新，可见水平形态的构造趋势。

图8 四川某三维区块的初始速度场Fig.8 Initial velocity field of a 3D block in Sichuan

图9 四川某三维区块分层反演后的速度场Fig.9 Velocity field after layered inversion of a 3 D block in Sichuan

高斯束弹性波偏移共成像点（CIP）道集的拉平程度可用于验证速度场的精度，分别利用初始速度场和反演重建速度场进行偏移并抽取相同位置的共成像点道集进行对比，同相轴拉平程度高的对应速度场越精确。对比四川某工区实际资料的纵、横波共成像点道集可以看出，相比初始速度场，使用分层反演后的纵、横波速度偏移得到的共成像点道集同相轴拉平程度均更高、同相轴更清晰，能量亦更为集中，说明优化后速度场的精度更高（图10）。