数据驱动的快速构造约束叠前地震多道反演方法

本发明针对油气地震勘探领域，涉及复杂地质构造区域和地震资料质量低的叠前地震多道反演，特别是高倾角构造区的弹性参数确定性预测。

背景技术：

叠前地震反演作为油气地球物理勘探领域获取地下介质弹性、物性、岩性和含流体性质的关键技术，是储层定量预测和含油气性检测的重要途径。基于地震振幅随偏移距或角度变化(amplitudevariationwithoffset/angle)信息，叠前地震反演通过地下介质弹性、物性参数与地层反射系数的定量关系，实现储层岩性圈闭刻画和含油气性展布研究。然而，叠前地震反演具有强不确定性和不稳定性，反演参数的空间连续性有待改善。基于此，将多道地震处理技术应用于叠前地震反演的研究越来越广泛。多道地震处理技术考虑地震数据相邻地震道的相关性，能够有效压制随机噪声，提高地震解释的稳定性。多道地震处理技术在地震勘探领域已经被广泛研究，例如：噪声压制、地震资料处理和地震资料插值等。

近年来，多道地震处理技术由于其稳定性强而发展迅速，尤其是多道地震反演技术。一维横向无约束多道地震反演(laterallyunconstrainedinversion(1d-lui))是一种经典的多道反演方法，与单道地震反演相比，能有效压制噪声。实际上，1d-lui仅求解反演参数的全局最优解，未考虑相邻地震道的内在联系，反演参数的连续性有待提高。与之相比，一维横向约束多道反演技术(laterallyconstrainedinversion(1d-lci))具有更好的横向连续性，该技术假设地下介质为水平层状分布，通常用一个全变差算子最为横向约束。1d-lci在地下结构简单、成层分布时，通常能取得好的应用效果，然而当地下结构复杂时，该方法容易出现地层边界过于平滑或边界模糊的现象，降低反演参数的分辨率，尤其在高倾角区域。基于此，基于地下构造或地震反射特征的横向约束方法相继被提出。目前，多道地震反演的研究趋向于发展更合理、更符合实际情况的横向约束或空间约束。然而，横向约束算子引入目标泛函通常采取降维策略，即将地震数据逐道排列成一列向量，这不可避免的会产生大型矩阵，严重影响了计算机存储和运算效率，制约其广泛实际应用。为克服这一问题，许多相应的快速算法相继发展，如blockcoordinatedencent(bcd)快速算法，该算法的思想是在每次迭代中逐道估计局部最优解直至收敛，这使得实现大规模的地震资料多道反演相对容易。然而，bcd算法在每次迭代中循环估计各个地震道，解决了计算机存储问题，但计算效率提高并不明显，并且逐道估计局部最优解来逼近全局最优解，不能保证反演参数为全局的最优解，且地震反演的精度仍需要进一步得到改善。

技术实现要素：

本发明要解决的就是复杂构造区域叠前地震反演缺乏稳定性和连续性的问题。

为实现上述目的，解决现有方法存在的关键问题，本发明采用下述技术方案：一种数据驱动的快速构造约束叠前地震多道反演方法，其包括如下步骤：

步骤1，估算地层倾角，并在旋转坐标系下分别构建平行于地层倾角和垂直于地层倾角的构造算子；

步骤2，基于地震数据的互相关系数描述地下介质的地震反射特征，并构建数据局部优化算子；

步骤3，构建数据驱动的快速构造约束叠前地震多道反演方法的目标泛函；

步骤4，构建数据驱动的快速构造约束叠前地震反演的优化算法。

优选的，所述步骤1包括：

将地层倾角定义为水平方向与地震数据变化梯度最小方向的夹角，则估算地层倾角的表达式为：

φ(x,z)＝tan^-1(rx(x,z)/rz(x,z))(1)，

式(1)中，φ为地层倾角，x、z分别为横向、纵向位置，rx(x,z)和rz(x,z)表示地震数据沿x、z方向的一阶偏导数；

利用地层倾角构建旋转算子的表达式为：

式(2)和(3)中，和为旋转坐标系下的旋转算子，为第i采样点估计的地层倾角；

设rx和rz分别为沿x和z方向的全变差算子，将rx和rz旋转，得到平行和垂直于地层倾角方向的构造算子的表达式为：

式(4)和(5)中，rparl和rperp分别为平行和垂直于地层倾角方向的构造算子。

优选的，所述步骤2包括：

根据互相关的定义，计算互相关系数的公式c表示为：

式(6)中，s(i,j)表示第j道的第i个采样点的地震记录，j'为第j道的相邻地震道，w为用于相关系数计算的时窗，k表示相邻地震道内用于相关性分析的采样点与分析点的上下漂移时间，k限定在[-33]内；

构建数据局部优化算子h的表达式为：

式(7)中，hi,j和ci,j分别为数据局部优化算子h和互相关系数c在第j道的第i个采样点值，c0为互相关系数的阈值。

优选的，所述步骤3包括：

叠前地震多道正演模型表示为：

式(8)中，s和m为二维地震记录和待反演参数，sⁱ为第i道地震数据，i＝1,...,n，每道地震数据包含的角度个数为h，g为正演算子，mⁱ(i＝1,2,…,n)表示由纵波阻抗、横波阻抗和密度的自然对数组成的第i道待反演参数，sⁱ，g和mⁱ具体表达式为：

式(9)、(10)中，c1＝1+tan²(θ)，c2＝-8γ²tan²(θ)，θi(i＝1,2,…,h)为不同角度的入射角，矩阵w表示由地震子波时移得到的子波核矩阵，w(θi)为不同角度地震子波核矩阵(i＝1,2,…,h)，ip、is和id分别为纵波阻抗、横波阻抗和密度的自然对数，l为一阶差分算子；

地层倾角为零的特殊情况时，地层水平成层分布，则有：

rparl＝rx,rperp＝rz(11)，

rparl和rperp分别为横向和纵向全变差算子，标准构造约束多道反演的目标泛函表示为：

式(12)中，分别为式(8)中矩阵s、m展开的对角矩阵，为式(10)中矩阵g的对角矩阵，λ＞0、α＞0分别为平行、垂直地层倾角的正则化参数；

基于全变差正则化的多道地震反演方法是基于地层水平成层分布，所述基于全变差正则化的多道地震反演方法对反演参数在纵向和横向上分别采用全变差算子进行约束，所述基于全变差正则化的多道地震反演方法的目标泛函为：

式(13)中，为佛罗贝尼乌斯范数，

在式(13)的基础上，根据构造算子的定义，采用阿达玛乘积算子将构造算子引入，构建快速构造约束叠前地震多道反演方法的目标泛函表达式为：

式(14)中，为阿达玛乘积算子，表示两个相同规模的矩阵对应元素相乘；qcos和qsin为与地震数据相同规模的旋转算子，其形式为：

式(15)、(16)中，为第i道第j个采样点的地层倾角；

通过阿达玛乘积算子将数据局部优化算子h引入式(14)，得到数据驱动的快速构造约束叠前地震多道反演方法的目标泛函：

式(17)中，算子h控制了每个采样点地震数据对反演的贡献，本发明进一步根据相邻地震道反演，通过项恢复地震道每个采样点处的反演参数，β＞0为地震反演横向约束的正则化算子。

优选的，所述步骤4包括：

将式(14)改写成一个约束极小化问题，表示为：

式(18)中，x1和x2为中间变量，表示反演参数沿纵向和横向的一阶偏导；

根据交替方向乘子算法，式(18)的拉格朗日形式为：

式(19)中，z1和z2分别为与x1和x2对应的对偶变量，采用交替方向乘子算法对式(19)进行求解，其中每个变量的更新迭代规则为：

根据式(20)至式(24)交替估计m、x1、x2、z1和z2这五个参数；

将式(17)改写成一个约束极小化问题，表示为：

式(25)中，为中间变量，根据交替方向乘子算法将式(25)进一步写成拉格朗日的形式：

式(26)中，z3为与对应的对偶变量，采用交替方向乘子算法对式(26)进行求解，其中每个变量的更新迭代规则为：

根据式(27)至式(33)，交替估计m、x1、x2、z1、z2和z3这七个参数，直至目标泛函收敛。

进一步验证反演方法的可行性，可行性验证包括两个部分：

第一步：数值实验测试方法可行性，设置理论弹性模型，基于常规数值模型建立合成地震记录，首先在无噪声干扰的情况下，分别采用快速构造约束叠前地震反演方法和数据驱动的快速构造约束叠前地震反演方法进行反演，测试反演结果与真实模型的差异，然后在合成地震记录上添加噪声进行反演，测试反演方法的稳定性；

第二步：实际地震资料处理测试本发明方法的有效性，选取地质构造复杂的二维地震剖面，同样分别采用快速构造约束叠前地震反演方法和数据驱动的构造约束叠前地震反演方法进行反演，通过与实际测井数据进行对比，能够提高反演参数的稳定性和可靠性。

本发明方法从地震反射系数近似方程出发，充分考虑地下介质局部连续性和地震资料本身信噪比对反演结果的影响，提高了地下介质各弹性参数反演的稳定性、连续性和可靠性，能实现地下介质各弹性参数稳定反演。

附图说明

附图是用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。为了更清楚地说明本发明的具体实施过程或技术方案，附图如下：

图1本发明公开的一种数据驱动的快速构造约束叠前地震多道反演方法的流程示意图。

图2估算地层倾角和旋转坐标系下构建构造算子示意图。

图3amarmousi2纵波阻抗真实模型示意图。

图3bmarmousi2横波阻抗真实模型示意图。

图4a受噪音干扰的合成小角度地震数据示意图。

图4b根据受噪音干扰的地震记录估算的地层倾角示意图。

图4c无噪音干扰的合成小角度地震记录示意图。

图4d根据无噪音干扰的地震记录估算的地层倾角示意图。

图5a基于常规构造约束反演的纵波阻抗示意图。

图5b基于常规构造约束反演的横波阻抗示意图。

图5c基于快速构造约束反演的纵波阻抗示意图。

图5d基于快速构造约束反演的横波阻抗示意图。

图5e基于数据驱动的快速构造约束反演的纵波阻抗示意图。

图5f基于数据驱动的快速构造约束反演的横波阻抗示意图。

图6a常规构造约束反演方法相对误差随迭代次数的变化示意图。

图6b常规构造约束反演、快速构造约束反演和数据驱动的快速构造约束反演方法相对误差对比示意图。

图7a实际资料基于常规构造约束反演的纵波阻抗示意图。

图7b实际资料基于常规构造约束反演的横波阻抗示意图。

图7c实际资料基于快速构造约束反演的纵波阻抗示意图。

图7d实际资料基于快速构造约束反演的横波阻抗示意图。

图7e实际资料基于数据驱动的快速构造约束反演的纵波阻抗示意图。

图7f实际资料基于数据驱动的快速构造约束反演的横波阻抗示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明充分考虑地下介质局部连续性和地震资料本身信噪比对反演结果的影响，研发了一种数据驱动的快速构造约束叠前地震多道反演技术(详见图1流程示意图)，提高了地下介质各弹性参数反演的稳定性、连续性和可靠性。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

本发明一种数据驱动的快速构造约束叠前地震多道反演方法，包括如下步骤：

步骤1，估算地层倾角，并在旋转坐标系下分别构建平行于地层倾角和垂直于地层倾角的构造算子；

步骤2，基于地震数据的互相关系数描述地下介质的地震反射特征，并构建数据局部优化算子；

步骤3，构建数据驱动的快速构造约束叠前地震多道反演方法的目标泛函；

步骤4，构建数据驱动的的快速构造约束叠前地震反演的优化算法。

步骤1中，将地层倾角定义为水平方向与地震数据变化梯度最小方向的夹角，则估算地层倾角的表达式为：

φ(x,z)＝tan^-1(rx(x,z)/rz(x,z))(1)，

式(1)中，φ为地层倾角，x、z分别为横向、纵向位置，rx(x,z)和rz(x,z)表示地震数据沿x、z方向的一阶偏导数。这种计算地层倾角的方法对信噪比高的地震数据通常能够准确的估计地层倾角，而当地震数据信噪比低(信噪比小于4)时，对地层构造描述精度有待改善，该种方法反正切函数会放大噪声的影响，对地震数据质量有一定的依赖性。通常采用平滑的方式，降低噪声影响，但这是以牺牲一些细节构造信息为代价。

本发明对地震多道反演添加平行和垂直地层倾角方向的构造约束，以保留反演参数的构造特征，这通常通过构造算子实现。估算地层倾角就是为了构建更合理的构造算子，在估算地层倾角之后，利用旋转坐标系(图2)可以构建构造算子。图中xr、zr坐标轴为旋转之后的坐标轴，旋转角度为地层倾角φ。

构造算子的构建可以从一阶全变差算子rx和rz(沿x和z方向的全变差算子)出发，利用地层倾角可以构建如下的旋转算子：

式(2)和(3)中，和为旋转坐标系下的旋转算子，为第i采样点估计的地层倾角。将rx和rz旋转，可以得到平行和垂直于地层倾角方向的构造算子，其具体形式为：

式(4)和(5)中，rparl和rperp分别为平行和垂直于地层倾角方向的构造算子。构造算子的精度直接关系到反演参数的稳定性和准确性。在估算地层倾角之前进行随机噪声压制、连续性处理等，可以有效地提高倾角估计的准确性。

步骤2中，基于地震数据互相关评价地震数据每个采样点的可靠性，利用相关系数控制每个采样点地震数据对反演参数的贡献，以降低噪音干扰严重、极度不连续或死道等采样点处地震数据对反演的贡献，提高地震反演的可靠性和连续性。根据互相关的定义，计算互相关系数的公式c表示为：

式(6)中，s(i,j)表示第j道的第i个采样点的地震记录，j′为第j道的相邻地震道，w为用于相关系数计算的时窗，k表示相邻地震道内用于相关性分析的采样点与分析点的上下漂移时间，根据地震数据在局部一般不会发生剧烈抖动，将k限定在[-33]内。本发明利用互相关系数控制每个采样点地震数据对反演的贡献，这一过程通过数据局部优化算子h实现，算子h的表达式为：

式(7)中，hi,j和ci,j分别为数据局部优化算子h和互相关系数c在第j道的第i个采样点值，c0为互相关系数的阈值。将互相关系数小于该值的采样点视为不可靠的点，控制该点地震数据对反演参数不做贡献，根据相邻地震道估计这些采样点的反演参数。

步骤3中，常规多道地震反演方法易产生大型矩阵，不利于计算机运算求解。本发明为避免大型矩阵的产生，提高反演效率，发展了一种引入构造算子的方法。

以反演纵波阻抗、横波阻抗和密度参数为例，叠前地震多道正演模型通常可以表示为：

式(8)中，s和m为二维地震记录和待反演参数，sⁱ为第i道地震数据i＝1,...,n，每道地震数据包含的角度个数为h，g为正演算子，mⁱ(i＝1,2,…,n)表示由纵波阻抗、横波阻抗和密度的第i道自然对数组成的待反演参数。sⁱ，g和mⁱ具体表达式为：

式(9)、(10)中，c1＝1+tan²(θ)，c2＝-8γ²tan²(θ)，θi(i＝1,2,…,h)为不同角度的入射角，矩阵w表示由地震子波时移得到的子波核矩阵，w(θi)为不同角度地震子波核矩阵(i＝1,2,…,h)，ip、is和id分别为纵、横波阻抗和密度的自然对数，l为一阶差分算子。

首先考虑地层倾角为零的特殊情况，构造算子为：

rparl＝rx,rperp＝rz(11)

即地层水平成层分布，构造算子分别为横向和纵向全变差算子。此时，标准构造约束多道反演的目标泛函可以表示为：

式(12)中，分别为式(8)中矩阵s、m展开的对角矩阵，为式(10)中矩阵g的对角矩阵，λ＞0、α＞0分别为平行、垂直地层倾角的正则化参数。

然而，基于全变差正则化(totalvariationregularization,简称tv正则化)的多道地震反演方法是基于地层水平成层分布，该方法对反演参数在纵向和横向上分别采用全变差算子进行约束，其反演的目标泛函为：

式(13)中，为佛罗贝尼乌斯范数(frobenius范数，简称f-范数)。可见，式(12)和式(13)是等价的，不同的是，式(13)具有极高的反演效率。本发明在式(13)的基础上，根据构造算子的定义，采用阿达玛乘积算子(即hadamard乘积算子)将构造算子引入目标泛函，其表达式为：

式(14)中，“”为hadamard乘积算子，表示两个相同规模的矩阵对应元素相乘。qcos和qsin为与地震数据相同规模的旋转算子，式(2)和式(3)中的和实际为矩阵qcos和qsin展开的对角阵，和矩阵规模远大于qcos和qsin，严重影响反演效率，qcos和qsin的形式为：

式(15)、(16)中，为第i道第j个采样点的地层倾角。

式(14)为本发明发展的快速构造约束叠前地震多道反演方法的目标泛函，该方法利用二维多道正演模型的运算优势，在全变差正则化的基础上采用hadamard乘积算子引入构造算子，避免了大型矩阵的产生。

通过hadamard乘积算子将数据局部优化算子h引入目标泛函(14)，可得到：

式(17)中，算子h控制了每个采样点地震数据对反演的贡献，本发明进一步根据相邻地震道反演，通过项恢复地震道每个采样点处的反演参数，β＞0为地震反演横向约束的正则化算子。式(17)即为数据驱动的快速构造约束叠前地震多道反演方法的目标泛函，该方法能有效降低地震数据质量对反演结果的影响。

步骤4中，直接求解式(17)具有一定困难，首先研究式(14)的求解策略。从形式上看，快速构造约束多道地震反演的目标泛函(14)与常规多道地震反演的目标泛函是不同的，它包括frobenius范数和hadamard乘积算子两种特殊运算符，这增加了目标泛函的求解难度。本发明从两种运算符的基本定义出发，基于交替方向乘子算法，即admm(thealternatingdirectionmethodofmultipliers,简称admm)算法，详细推导了其求解过程。

直接对式(14)进行求解是非常困难的，本发明首先将其改写成一个约束极小化问题，可以表示为：

式(18)中，x1和x2为中间变量，表示反演参数沿纵向和横向的一阶偏导。根据交替方向乘子算法，式(18)的拉格朗日形式为：

式(19)中，z1和z2分别为与x1和x2对应的对偶变量，交替方向乘子算法在对偶交替估计过程中使反问题的解逐步逼近最优解，具有更高的计算效率和精度。采用交替方向乘子算法对式(19)进行求解，其中每个变量的更新迭代规则为：

根据式(20)～(24)交替估计m、x1、x2、z1和z2这五个参数，通常就能够收敛。该算法不涉及大型矩阵，且在每次迭代中对所有地震道同时估计，避免了块坐标下降算法逐道估计局部极小值的不足，提高了反演效率和精度。

快速构造约束叠前地震多道反演算法的伪代码如下所示：

数据驱动的快速构造约束叠前地震多道反演的目标泛函同样包含frobenius范数和hadamard乘积算子两种特殊运算符，难以直接求解。与式(14)的求解过程一样，本发明基于交替方向乘子算法详细推导了求解过程。

首先将式(17)改写成一个约束极小化问题，可以表示为：

式(25)中，为中间变量，根据交替方向乘子算法将式(25)进一步写成拉格朗日的形式：

式(26)中，z3为与对应的对偶变量，采用交替方向乘子算法对式(26)进行求解，其中每个变量的更新迭代规则为：

该算法在快速构造约束多道反演方法的基础上降低了地震资料质量的影响。根据上述公式交替估计m、x1、x2、z1、z2和z3这七个参数，直至目标泛函收敛。

数据驱动的快速构造约束叠前地震多道反演算法的伪代码如下所示：

进一步验证反演方法的可行性，可行性验证包括两个步骤：

第一步：数值实验测试方法可行性，设置理论弹性模型，基于常规数值模型建立合成地震记录，首先在无噪声干扰的情况下，分别采用快速构造约束叠前地震反演方法和数据驱动的快速构造约束叠前地震反演方法进行反演，测试反演结果与真实模型的差异，然后在合成地震记录上添加噪声进行反演，测试反演方法的稳定性；

第二步：实际地震资料处理测试本发明方法的有效性，选取地质构造复杂的二维地震剖面，同样分别采用快速构造约束叠前地震反演方法和数据驱动的构造约束叠前地震反演方法进行反演，通过与实际测井数据进行对比，能够提高反演参数的稳定性和可靠性。

为进一步说明本发明的可行性和有效性，下面列举两个实施例：

实施例1：数值模型测试，详见图3、图4、图5和图6。

本发明选取部分marmousi-2纵、横波阻抗(图3a和图3b)和密度模型对该方法进行测试。用于测试的真实模型共包含851道，每道606个采样点，采样间隔为2ms。利用30hz雷克子波和真实阻抗模型合成地震记录，在合成地震记录中随机去除50％地震道的部分地震数据，并添加10％的高斯噪声以测试反演方法的稳定性。由于地震数据具有带限性质，反演结果缺乏低频成分，本文将真实阻抗模型进行0-5hz低通滤波得到低频阻抗模型，用来补充反演结果的低频成分。

本文基于叠后地震记录(图4c)估计地层倾角和构建构造算子。本发明提出的数据驱动的快速构造约束叠前地震多道反演方法需利用互相关来评价地震数据的合理性，首先对合成的角度道集进行相关性分析，并计算互相关系数。为方便成图，现有的常规构造约束反演方法可以简称为sci(structurallyconstrainedinversion)，本发明将提出的快速构造约束反演简称为fsci(faststructurallyconstrainedinversion)，其目标泛函表示为式(17)，本发明进一步提出了数据驱动的快速构造约束叠前地震多道反演简称为dfsci(data-drivenfaststructurallyconstrainedinversion)，其目标泛函表示为式(25)。将这三种方法的反演连续性、可靠性和效率分别进行了对比分析。

图5为三种方法纵、横波阻抗的反演结果。从整体上看，与真实阻抗模型相比，三种方法均与真实模型相近，表明了本发明反演测试工作是可靠的。另外，常规构造约束反演方法(sci)与快速构造约束反演方法(fsci)反演结果基本一样，构造特征明显，地层边界清晰，尤其在高倾角区域，这说明快速构造约束反演能够保留常规构造约束反演方法的优势。与前两种方法相比，数据驱动的快速构造约束反演方法(dfsci)具有更高的连续性和信噪比(图5e和图5f)，表明本发明提出的数据驱动的快速构造约束叠前地震多道反演方法能够保证反演参数的稳定性、可靠性和连续性。

图6对三种方法的反演效率进行了对比，从图6a中可看出，常规构造约束反演的收敛效率比较低，需要数十次甚至上百次迭代才能收敛，且从图6b可以看出，该方法迭代一次需要更长的时间，而本发明提出的两种方法收敛速度极快，收敛误差小，也说明了这两种方法具有高的计算精度和效率。

实施例2：实际资料处理，详见图7。

本发明在常规构造约束反演方法(sci)的基础上发展了快速构造约束反演方法(fsci)，为改善对地震资料质量的影响又进一步提出了数据驱动的快速构造约束反演方法(dfsci)。并且分别利用三种方法进行了叠前地震反演测试。图7a和7b为常规构造约束反演方法反演的纵横波阻抗，可以看出，该方法能够保证地层边界清晰，构造特征明显。图7c和7d为本文提出的快速构造约束反演方法反演的纵横波阻抗，与常规方法相比，该方法的构造特征和地层边界更为清晰，这是因为该方法为全局优化算法，具有更高的计算精度。图7e和7f在保证地层边界和构造特征的同时，进一步提高了反演参数的横向连续性和信噪比，这表明本文提出的两种方法具有更好的实际应用效果。