一种起伏地表的地下构造成像方法与流程

本发明涉及一种起伏地表的地下构造成像方法，属于油气物探工程领域。

背景技术：

随着我国勘探开发的深入，勘探重点区域逐渐转移到西部探区。西部探区的勘探开发存在如下主要特点：(1)地表起伏剧烈，近地表构造复杂；(2)背景构造复杂，例如逆掩推覆和盐下成像等；(3)储层埋藏较深，信号较弱。因此，西部探区勘探的重点问题是如何实现起伏地表条件下的高效、高精度保幅成像。

当前的常规反演成像方法通过最速下降法、共轭梯度法等迭代算法求解，虽然可以得到高质量的成像结果，但计算量过于庞大。如果对观测数据进行随机采样，在每次更新模型时只计算部分数据，可以减小计算量。然而，随机采样策略导致反演成像时的误差函数和梯度都具有随机性。由于对观测数据随机采样，当地下构造复杂或地表起伏剧烈时，相邻两次迭代中偏移数据间相关性较差，随机梯度下降法的收敛曲线常出现剧烈抖动，在进行地下结构成像时成像质量较差，收敛速度较慢，计算效率低。

技术实现要素：

本发明的目的在于，提供一种起伏地表的地下构造成像方法，它可以解决当前技术中存在的问题，提高起伏地表条件下成像时的计算效率，同时不影响其成像质量。

为解决上述技术问题，本发明采用如下的技术方案：一种起伏地表的地下构造成像方法，包括以下步骤：

s1，输入观测数据，偏移速度场，同时预先设置反射系数模型的初始值；

s2，对观测数据进行基于相关系数的采样或重新采样，得到偏移数据；

s3，根据当前的反射系数模型和偏移数据获得数据残差；

s4，根据数据残差，确定随机梯度和/或随机共轭梯度，基于随机梯度和/或随机共轭梯度进行反射系数模型迭代，获得最终的反射系数模型。

与现有技术相比，用本发明的方法在起伏地表条件下成像时，在减小偏移数据量的同时能有效地减小偏移数据随迭代次数的变化，进而减弱更换偏移数据对梯度方向共轭性的破坏及收敛曲线的抖动现象，提高了迭代时的收敛稳定性及其成像质量。此外，本发明将最小二乘反演思想中引入起伏地表逆时偏移，能够压制成像噪音、均衡成像振幅、提高成像分辨率，能够对近地表构造及深层构造高精度成像。

附图说明

图1为本发明的一个实施例的流程框图；

图2为本发明的一个实施例的偏移速度场；

图3为本发明的一个实施例的起伏地表高程信息；

图4为本发明的一个实施例的观测地震数据；

图5为本发明的一个实施例的对观测数据进行基于相关系数炮数据采样后得到的偏移数据；

图6为本发明的一个实施例的起伏地表逆时偏移结果；

图7为基于随机梯度下降法的起伏地表最小二乘逆时偏移结果；

图8为本发明的一个实施例的起伏地表最小二乘逆时偏移结果；

图9为本发明一个实施例与其他两种成像方法结果的单道振幅对比；

图10为本发明一个实施例与其他两种成像方法结果的归一化数据残差收敛曲线；

图11为本发明一个实施例与其他两种成像方法的计算时间对比；

图12是本发明的一个实施例流程示意图。

下面结合附图和具体实施方式对本发明作进一步的说明。

具体实施方式

本发明的实施例1：一种起伏地表的地下构造成像方法，如图1所示，包括以下步骤：

s1，输入观测数据，偏移速度场，同时预先设置反射系数模型的初始值；这里的观测数据包括但不限于起伏地表高程、观测系统、野外观测记录(即炮数据)和可以预设的偏移参数等，其中偏移参数包括偏移速度场横向采样点nx及纵向采样点nz、空间采样间隔dx、时间采样间隔dt、时间采样点数nt、主频f0、采样后得到的偏移数据中的炮数ns及模型更新方向控制因子o和迭代终止的阈值。通常的，反射系数模型的初始值被设置为0，偏移速度场是对观测数据进行常规的速度分析后得到。

s2，对观测数据进行基于相关系数的采样或重新采样，得到偏移数据；这里的相关系数计算方法为：

其中，dj(x,z),为上次迭代时偏移数据中的单炮数据及其均值，di(x,z),为剩余数据中的单炮数据及其均值。在公式(1)中，首先对di(x,t),dj(x,t)进行去均值处理，减弱非零均值噪声对相关系数的影响，因此该公式具有一定的抗噪性。然后，分别计算di(x,t)与上次迭代时偏移数据中每个单炮数据的相关系数。最后，对所有的相关系数求和得到di(x,t)与偏移数据的相关系数。由式(1)-(3)可知，相关系数可以反映某一单炮数据与上次迭代时偏移数据的相关性，以相关系数作为参考进而发展出基于相关系数的如下采样策略：

在首次采样时，就是在第1次迭代时，按照固定的间隔采样ns个沿偏移距方向均匀分布的观测数据作为偏移数据，其中采样间隔大于炮间隔(采样间隔等于炮间隔时采样结果为整个数据体；只有采样间隔大于炮间隔时采样结果才是部分炮数据)。当迭代次数大于1并且满足重采样条件时，计算剩余数据中每炮数据与上一次迭代时偏移数据的相关系数，选取相关系数最大的ns个单炮数据作为本次迭代的偏移数据；当剩余数据中单炮数据的数量不足ns时，不足的部分从已计算过的数据中随机选取。当不满足重采样条件时，直接将上一次迭代时偏移数据作为本次迭代的偏移数据。需要说明的是，这里的剩余数据是指观测数据中去掉以前所有迭代时所采样的偏移数据后，剩余的数据。重新采样条件为：当迭代次数k除以模型更新方向控制因子o的余数为1时，对观测数据进行重采样，更换偏移数据；否则，无需对观测数据进行重采样，直接将上次迭代时的偏移数据作为本次迭代的偏移数据。

在随机炮数据采样中，偏移数据随迭代次数随机变化，当相邻两次迭代中偏移数据相关性较差时，其收敛曲线常产生剧烈抖动，收敛速度较慢，计算效率低。而基于相关系数的炮数据采样策略能够保证相邻两次迭代中偏移数据的差异最小，也就是该采样策略不受地下构造或地表条件的影响，能有效地减小偏移数据随迭代次数的变化，进而减弱更换偏移数据对梯度方向共轭性的破坏及收敛曲线的抖动现象，提高了迭代时的收敛稳定性，进而提高了计算效率。同时，在基于相关系数的炮数据采样策略中加入了重采样判定条件，不满足重采样条件时，偏移数据不变，这样就保证了梯度方向的共轭性完全成立，改善迭代时的收敛效率和成像质量。

s3，根据当前的反射系数模型和偏移数据获得数据残差；具体的说，对于反射系数模型，通过逆时反偏移模拟获得模型数据p(x)，逆时反偏移模拟模型数据所用的计算公式为：

p(x)＝ω²∫m(x')f(ω)g0(x'；xs)g0(x；x')dx'(4)

其中，p(x)代表逆时反偏移得到的模型数据，m(x')代表反射系数模型，f(ω)为震源，ω为频率，g0(x'；xs),g0(x；x')为背景速度对应的格林函数，xs为炮点位置，x,x'分别为地下任意一点位置。

式(4)可通过两步正演模拟求解：

其中，p0(x')为在偏移速度v0(x)下求得的背景波场，将其与反射系数模型的乘积作为震源再次激发，即可求得逆时反偏移模拟的模型数据。根据式(5)和式(6)进行波场延拓时，每经过一次延拓都需要对波场施加自由边界条件：

pπ(x,z)＝0(7)

其中，π代表起伏地表的位置，即将起伏地表之上的波场赋值为0。

基于上述说明，公式(4)可以用矩阵算子形式简化表示为：

d'＝lm(8)

其中，m为反射系数模型的矩阵形式，d'为模型数据的矩阵形式，l为逆时反偏移算子，即起伏地表最小二乘逆时偏移中的正算子。也就是说通过对反射系数模型做逆时反偏移运算即可得到模型数据p(x)，在本发明方法中，在每次迭代中只需要计算偏移数据对应的模型数据，表示为：

d's＝lsm(9)

其中，d's,ls分别为偏移数据对应的模型数据及逆时反偏移算子，得到模型数据后，与当前的偏移数据相减，获得数据残差和分别为第k次迭代时的偏移数据及其对应的逆时反偏移算子，m^(k-1)为第k-1次迭代计算得到的反射系数模型。

s4，根据数据残差，确定模型更新梯度，基于模型更新梯度进行反射系数模型迭代，获得最终的反射系数模型，这里的梯度包括了随机梯度和/或随机共轭，具体的说，通过数据残差计算随机梯度的计算公式为：

其中，为的共轭转置，即逆时偏移算子。

随机共轭梯度修正因子的计算公式为：

则随机共轭梯度的计算公式为：

然后，每一步的迭代中，都要基于判定条件来选择随机梯度或随机共轭梯度，作为模型更新梯度进行反射系数模型迭代，获得最终的反射系数模型，具体包括以下步骤：

s41，判断迭代次数k除以模型更新方向控制因子o的余数是否为1，若是，则选择随机梯度为模型更新梯度，作为反射系数模型迭代的方向；若否，则选择随机共轭梯度为模型更新梯度，作为反射系数模型迭代的方向。

s42，当所述的模型更新梯度小于阈值时，进入步骤s44，这里的阈值是可以根据经验和实际精度的要求调整设置的，当所述的模型更新梯度大于阈值时，根据模型更新梯度计算更新步长，更新步长的计算公式为：

其中，为模型更新梯度，包括随机梯度和/或随机共轭梯度，具体迭代中选用哪个参数在步骤s41中确定。

s43，基于更新步长，迭代更新反射系数模型，反射系数模型的更新公式为：

并返回步骤s2，

s44，判定此时获得的反射系数模型为最终的反射系数模型，由最终的反射系数模型获得地下构造成像结果。采用本发明的方法，通过引入基于相关系数的采样策略，避免了随机采样，在起伏地表的条件下成像时，在减小偏移数据量的同时能有效地避免偏移数据随迭代次数产生巨大变化，进而减弱更换偏移数据对梯度方向共轭性的破坏及收敛曲线的抖动现象，提高了迭代时的计算效率。此外，本发明将最小二乘反演思想中引入起伏地表逆时偏移，能够压制成像噪音、均衡成像振幅、提高成像分辨率，对近地表构造及深层构造高精度成像；迭代过程中利用混合梯度，在一部分迭代中使用了随机共轭梯度求解，相对于随机梯度收敛更加稳定，成像效果更好。

实施例2

为了更具体地说明本发明的方法，以加拿大逆掩断层模型为例说明本发明的方法。具体流程如图12所示。输入观测数据，包括偏移速度场(如图2所示)、起伏地表高程(如图3所示)、观测系统、野外观测数据(如图4所示)、迭代终止的阈值及预设偏移参数，并设置初始反射系数模型的值为0，观测系统分布为：在起伏地表以30米间隔均匀分布278炮，每炮都是556个检波点接收，检波点间隔为15米，偏移参数如下：偏移速度场横向采样点为556及纵向采样点为250，空间采样间隔为15米，时间采样间隔为0.5毫秒，时间采样点数为4000个，主频为25赫兹，采样后偏移数据的炮数为40，模型更新方向控制因子为5；对观测数据进行基于相关系数的炮数据采样，采样后的观测数据如图5所示；对于当前反射系数模型，通过逆时反偏移模拟模型数据，与偏移数据相减，计算数据残差；根据数据残差计算随机梯度下降方向；判断是否需要修正随机梯度方向，若满足判断条件，将随机梯度修正为随机共轭梯度，作为模型更新方向，否则，将随机梯度作为模型更新方向；判断模型更新方向是否满足阈值条件，若不满足阈值条件，根据模型更新方向计算更新步长并更新反射系数模型，再次对观测数据进行基于相关系数的炮数据采样、逆时反偏移、计算模型更新方向，直至模型更新方向满足阈值条件，若满足阈值条件，输出最终的反射系数模型，即最终的成像结果(如图8所示)。

从成像结果上看，采用本发明的方法，相对于起伏地表逆时偏移结果(如图6所示)，压制了成像噪音，补偿了深部的成像能量，近地表构造成像及中深部断层更加清晰，提高了成像分辨率；相对于基于随机梯度下降法的起伏地表最小二乘逆时偏移成像结果(如图7所示)，成像噪音更弱，成像振幅更加均衡，对近地表构造的成像更加清晰。

图9为从基于共轭梯度法的起伏地表最小二乘逆时偏移成像结果、基于随机梯度下降法的起伏地表最小二乘逆时偏移成像结果和本发明的成像结果中distance＝1000米处抽取的单道振幅曲线，本发明的成像结果振幅与基于共轭梯度法的起伏地表最小二乘逆时偏移成像结果振幅相近，相对于基于随机梯度下降法的起伏地表最小二乘逆时偏移成像结果振幅具有更好的保幅性。

图10为三种成像方法的归一化数据残差收敛曲线，本发明的数据残差与基于共轭梯度法的起伏地表最小二乘逆时偏移方法的数据残差都能够稳定收敛，而基于随机梯度下降法的起伏地表最小二乘逆时偏移成像方法的数据残差具有明显波动，并且收敛的速度慢于本发明方法。

图11为几种方法的运行时间对比，成像测试所用cpu型号为intel(r)xeon(r)cpue5-2650v2@2.60ghz。本发明成像方法的计算时间远远小于基于共轭梯度法的起伏地表最小二乘逆时偏移方法的计算时间，但二者成像质量相当，相对于基于随机梯度下降法的起伏地表最小二乘逆时偏移方法计算时间相近，但成像质量较高。因此，本发明方法相对与现有技术方法，同时兼顾较高的计算效率和成像质量，对西部探区的勘探具有重要意义。