一种三维全波形反演能量加权梯度预处理方法与流程

本发明涉及地球物理勘查技术领域，主要用于提高三维全波形反演纵波速度场的精度。

背景技术：

速度是描述地下介质情况的重要参数，地球物理勘探的重点在于如何恢复地下介质所有尺度的速度信息。速度分析的方法有很多，但是始终达不到理想的效果。射线层析、波动方程层析等层析成像的方法可恢复低波数的速度信息，偏移的方法可提供高波数的反射率信息，然而这些方法都不能同时恢复所有波数的速度参数。全波形反演利用叠前地震数据的全部信息，可提供地下介质的高分辨率成像。常规全波形反演是非线性梯度类最优化方法，以观测数据与模拟数据的最小二乘目标函数值最小为标准更新地下速度模型。三维全波形反演主要存在以下难点：(1)巨大的计算量问题：三维全波形反演基于三维波动方程正演模拟的迭代反演方法，在一次迭代过程中，需要计算多炮的波动方程正演模拟，误差波场反向传播，这两部分每个时间步都需要进行差分运算，伴随着巨大的计算量。(2)海量存储的问题：基于伴随状态梯度算子的全波形反演方法要求存储整个正演波场，三维情况下的存储量是不可容忍的，可以说大部分的时间消耗在数据的I/O读写上。计算量大可以通过多种并行算法进行加速，而I/O读写速度只能靠提高硬件性能来改善。(3)梯度算子的求取效率与精度：常规全波形反演梯度算子是由震源正传波场的时间二阶偏导数与残差波场反传波场互相关所得，震源正传波场包含着球面波传播过程中几何扩散的能量损失，震源波场与残差反传波场互相关就导致了梯度深、浅层能量更加不均衡，会引起反演速度深层精度不足。(4)海森矩阵的求取：对梯度预处理同样是全波形反演中的一个重要环节，可有效的提高反演的精度。然而海森矩阵的求取伴随着巨大的存储量和计算量，如何避开海森矩阵的求取或是对海森矩阵进一步的近似都是需要考虑的问题。如何在提高效率、降低耗费且保证精度是三维全波形反演的首要问题。因此需要发展一种三维全波形反演能量加权梯度预处理方法。

技术实现要素：

本发明的目的正是为了解决现有技术存在的问题，提供一种高效高精度的三维全波形反演能量加权梯度预处理方法。

本发明的技术方案为：

一种三维全波形反演能量加权梯度预处理方法，具体步骤如下：

(1)三维高精度有限差分正演模拟

第一，根据初始速度的道头信息确定三维正演模拟的观测系统，主要确定炮检波点的位置；第二，根据初始速度的最大最小值求取满足有限差分数值模拟差分稳定性和频散关系的三维正演模拟参数；第三，引入完全匹配层边界条件，确定其需要的参数，用于消除正演模拟的边界反射干扰；第四，用时间二阶空间十阶有限差分方法进行三维正演模拟，存储边界波场用于波场重建，并计算存储到达每个速度网格点振幅最大值。

(2)源波场重建及误差波场逆时传播求取梯度

正演模拟得到的模拟炮集与实际观测炮集对应做差求取残差波场。读取步骤(1)中存储的边界波场作为边界条件用时间二阶空间十阶有限差分方法进行源波场重建。源波场重建的同时进行残差波场逆时传播。根据式(1)的伴随状态法进行梯度计算，将源波场重建的波场与残差波场逆时传播得到的波场对应时刻进行零延迟互相关得到梯度算子。

其中E为目标函数，表示梯度，m为模型参数，v表示各网格点的速度值，x_s表示震源点的位置，x表示各网格点的位置，t表示每个时间步，T表示最大时间步，为源波场关于时间的二阶偏导数，q为以残差波场为震源的逆时反传波场。

(3)能量加权梯度预处理

取存储的波场到达每个速度网格点能量的最大值，即进而求取每点能量的最大值为波场到达每个速度网格点能量的最大值即初至波的能量值，表征的是波传播球面波几何扩散的过程，用此对梯度进行预处理，得到能量加权梯度算子。

其中，表示能量加权梯度算子。

(4)求取合适步长迭代更新速度

首先给一个试探步长，再用Armijo条件的一维线搜索方法求取合适的步长作用与能量加权梯度算子对速度进行迭代更新。

本发明的技术效果体现在：

常规全波形反演梯度算子是由震源正传波场的时间二阶偏导数与残差波场反传波场互相关所得，震源正传波场包含着球面波传播过程中几何扩散的能量损失，震源波场与残差反传波场互相关就导致了梯度深、浅层能量更加不均衡，会引起反演速度深层精度不足。利用海森矩阵对梯度进行预处理可有效的均衡梯度的能量，但是三维情况下的海森矩阵构建困难，计算量与存储量巨大，因此我们只提取海森矩阵中球面波传播几何扩散的信息对梯度进行预处理达到均衡梯度能量的效果，提高速度深层的反演精度。

该方法基于GPU集群搭建多级异构并行算法解决三维全波形反演的运算效率问题；利用源波场重建策略解决三维情况下的海量存储问题；在保证运算效率的情况下提高梯度算子构建精度。

附图说明

图1为截取的部分二维Marmousi模型并重采样作为真实速度

图2为对图1的真实速度平滑作为初始速度

图3为不同全波形反演算法迭代60次速度的对比，其中图a为常规算法反演速度；图b为本文算法反演速度

图4为常规全波形算法迭代60次与本文算法迭代60次误差下降曲线的对比

图5为不同全波形反演算法迭代60次单道速度的对比，其中图a为常规算法反演速度；图b为本文算法反演速度

图6为全波形反演能量加权梯度预处理方法流程图

图7为三维SEG/EAGE推覆体速度模型，其中图a为真实速度；图b为初始速度

图8为三维SEG/EAGE推覆体模型不同迭代次数全波形反演速度，其中图a为第1次迭代的三维速度体与三个不同方向的剖面；图b为第10次迭代的三维速度体与三个不同方向的剖面；图c为第40次迭代的三维速度体与三个不同方向的剖面

图9为三维SEG/EAGE推覆体模型不同迭代次数的误差下降曲线

图10为三维SEG/EAGE推覆体模型最终反演速度与真实速度以及初始速度的对比，其中图a为真实速度；图b为初始速度；图c为反演速度

具体实施方式

下面结合附图说明本发明的具体实施方式：

通过模型测试说明具体的技术方案：

第一步：三维高精度有限差分正演模拟

第一，根据初始速度的道头信息确定三维正演模拟的观测系统，主要确定炮检波点的位置；第二，根据初始速度的最大最小值求取满足有限差分数值模拟差分稳定性和频散关系的三维正演模拟参数；第三，引入完全匹配层边界条件，确定边界厚度以及衰减吸收系数，波场能量到达边界处逐渐衰减为零，从而消除三维正演模拟的边界反射干扰；第四，用时间二阶空间十阶有限差分方法对速度模型进行三维正演模拟，存储边界波场用于波场重建，并计算波场到达每个速度网格点振幅最大值。

第二步：源波场重建及误差波场逆时传播求取梯度

第一步中三维正演模拟得到的模拟炮集与实际观测炮集对应做差求取残差波场。读取第一步中存储的边界波场作为边界条件，用时间二阶空间十阶有限差分方法进行源波场重建。源波场重建的同时进行残差波场逆时传播。根据下式(1)的伴随状态法进行梯度计算，将源波场重建的波场与残差波场逆时传播得到的波场对应时刻进行零延迟互相关得到梯度算子。

式中E为目标函数，表示梯度，m为模型参数，v表示各网格点的速度值，x_s表示震源点的位置，x表示各网格点的位置，t表示每个时间步，T表示最大时间步，为源波场关于时间的二阶偏导数，q为以残差波场为震源的逆时反传波场。

第三步：能量加权梯度预处理

读取存储的源波场到达每个速度网格点振幅的最大值，即进而求取每点能量的最大值为(前一个公式是最大振幅值，后一个公式是最大能量值)。源波场到达每个速度网格点能量的最大值即初至波的能量值，表征的是波传播球面波几何扩散的过程，用此对梯度进行预处理，得到能量加权梯度算子，

其中，表示能量加权梯度算子。

第四步：求取合适步长迭代更新速度

每次迭代，首先给一个试探步长，再用一维线搜索方法求取新的步长，此步长满足Armijo条件使误差下降，然后将求得的步长作用于能量加权梯度算子对速度进行迭代更新。图3到图5展示的是二维Marmousi模型反演测试的效果图。从图3可以看出图b中的速度成像浅层与深层的能量较图a更加均衡。而从图4误差曲线下降图的也可以看出文中基于能量加权梯度算子的全波形反演算法误差收敛更快，且收敛的误差值更小。从图5为单道速度反演结果的对比可以看出基于能量加权梯度算子的全波形反演算法深部速度更新更加准确。

图8到图10展示的是按照图6流程进行三维SEG/EAGE推覆体模型反演测试的效果图。图8展示的是不同迭代次数三维SEG/EAGE推覆体模型全波形反演的速度结果。其中图a为第1次迭代的三维速度体与三个不同方向的剖面；图b为第10次迭代的三维速度体与三个不同方向的剖面；图c为第40次迭代的三维速度体与三个不同方向的剖面。三个不同方向分别为横纵测线方向8000m处的剖面与深度3000m处的切片。从图7中可以看出随着迭代次数的增加，构造细节成像越来越清晰，特别对于3000m处目的层切片的河道处的细节恢复明显。图9是迭代过程中的误差下降曲线，误差做过归一化。从图中可以看出，17次之前误差下降迅速，17次之后逐渐趋向收敛。图10是三维SEG/EAGE推覆体模型全波形反演迭代第40次的最终反演速度。图a是初始速度，图b是真实速度，图c是全波形反演速度。从图中可以看出速度深浅层成像效果均衡，并可以看出全波形反演在构造细节上优势明显，可以对类似于河道等的构造细节精细刻画。