基于各向异性马尔科夫随机域的叠前地震参数反演方法与流程

本发明涉及一种基于各向异性马尔科夫随机域的叠前地震参数反演方法，属于地震资料反演的技术领域。

背景技术：

地震反演是定量解释地震资料和预测地层弹性参数的重要手段，在油气资源的勘探和开发过程中应用广泛。根据地震反演所采用的正演模型，可以分为基于波动理论的波动方程反演和基于robinson褶积模型反演两大类。在实际工作中主要使用基于褶积模型的反演方法，该类方法易于实现，并能够满足精度要求。叠前地震反演是一种基于褶积模型的反演，它基于振幅随炮检距变化理论，直接利用信息丰富的叠前地震道集数据，可以同步获得纵波速度、横波速度和密度等多种弹性参数，在弹性参数反演和储层流体识别中发挥着不可替代的作用。

但叠前地震反演尚存诸多问题，如关于计算叠前道集反射系数，现有的文献及商业化软件多数使用精确左普利兹方程的线性近似式，但这些近似式需要一定的假设条件，会造成计算误差。此外，叠前地震的多参数反演是一个高度非线性问题，需要解决反演病态问题、地震带限问题。有关反演病态问题可以通过正则化方法得到解决，尤其是基于马尔科夫随机域的边界保护正则化，而带限问题可以使用测井数据的约束反演。在叠前地震反演方面，马尔科夫随机域技术没有得到足够的重视，现有研究中主要使用一些各向同性的马尔科夫随机域(或者各向同性的平滑因子)，这忽略了马尔科夫随机域的各向异性特点，尤其是对于层状和复杂介质。在层状和复杂介质中，水平、垂直和对角方向上模型梯度值的先验作用差异很大。因此，对传统的马尔科夫随机域技术需要进行改进和完善，以适应对于层状或复杂地层的叠前地震反演，在解决叠前多参数反演病态问题的同时，获取符合真实地质构造的弹性参数模型，以更好的指导复杂地质构造下油气资源的勘探、开发和利用。

技术实现要素：

本发明所要解决的技术问题在于克服现有技术的不足，提供一种基于各向异性马尔科夫随机域的叠前地震参数反演方法，解决现有方法忽略了马尔科夫随机域的各向异性特点，尤其是对于层状和复杂介质，存在叠前多参数反演的非线性和不稳定问题。

本发明具体采用以下技术方案解决上述技术问题：

基于各向异性马尔科夫随机域的叠前地震参数反演方法，包括以下步骤：

步骤一、建立叠前地震参数的目标函数，所述目标函数包含数据项和先验约束项；

步骤二、利用佐普利兹方程计算地层的纵波反射系数，得到的纵波反射系数与零相位ricker子波褶积后生成合成地震记录；根据测量所得地震记录与合成地震记录之间误差的二范数获得目标函数的数据项；

步骤三、利用各向异性扩散法中的边界终止函数计算各向异性马尔科夫随机域权系数，以得到数据点所在不同方向的权系数值；

步骤四、基于各向异性马尔科夫随机域将所得数据点所在不同方向的权系数值和纵波速度、横波速度先和密度先验约束项结合，获得目标函数的先验约束项；

步骤五、提取待反演区域的测井数据，对测井数据进行反演参数的对数线性拟合，得到拟合系数及拟合误差，用以控制目标函数极小寻优过程的稳定性；

步骤六、根据所得数据项和先验约束项代入所建立的叠前地震参数目标函数，及确定目标函数，并基于快速模拟退火算法和拟合系数及拟合误差进行目标函数的极小寻优；

步骤七、完成对目标函数的迭代寻优，输出反演结果。

进一步地，作为本发明的一种优选技术方案，所述步骤一建立叠前地震参数的目标函数为：

j(vp,vs,ρ)＝j1(vp,vs,ρ)+j2(vp,vs,ρ)

其中，j1为数据项，j2为先验约束项，vp、vs和ρ分别为待反演地层的纵波速度、横波速度和密度参数。

进一步地，作为本发明的一种优选技术方案，所述步骤二获得目标函数的数据项为：

其中，d为测量所得地震记录；w为零相位ricker子波；r为纵波反射系数；θ为入射角；t为地震记录的采样时间。

进一步地，作为本发明的一种优选技术方案，所述步骤三计算各向异性马尔科夫随机域权系数α采用公式:

其中，g(▽v)为边界终止函数；▽v为梯度值；q为扩散系数，为预设定的常数。

进一步地，作为本发明的一种优选技术方案，所述步骤四获得目标函数的先验约束项为：

其中，c为数据点邻域内点集；αc为各向异性马尔科夫随机域权系数；d为相邻两个数据点的一阶差分；φ(t)为势函数；vp、vs和ρ分别为待反演地层的纵波速度、横波速度和密度参数，δp、δs和δd分别为与之对应的正则化参数。

进一步地，作为本发明的一种优选技术方案，所述步骤五对测井数据进行反演参数的对数线性拟合采用公式：

lnvs＝k·lnvp+kc+δls

lnρ＝m·lnvp+mc+δld

其中，vp为待反演的纵波速度；vs为待反演的横波速度；ρ为待反演的密度参数；δls和δld为拟合误差；kc和mc为拟合系数。

本发明采用上述技术方案，能产生如下技术效果：

本发明是一种基于各向异性马尔科夫随机域的叠前地震多参数的反演方法，该方法通过各向异性扩散法计算各向异性马尔科夫权系数，用以校正层状或复杂地层对反演目标函数先验约束项的影响，并基于精确佐普里兹方程进行正演计算，并采用测井数据的对数线性拟合关系来控制目标函数全局优化反演的稳定性。分别对合成地震数据及实际地震资料进行叠前反演，并取得准确的反演结果。该方法有效解决叠前多参数反演的非线性和不稳定问题，可以有效的保护反演地层参数的横向连续性和纵向分层性，特别适用于复杂地质条件下的弹性参数反演和储层参数预测，对油气资源的勘探和开发具有重要指导和现实意义。

本发明的方法，可以有效解决叠前多参数反演的非线性和不稳定问题，对于层状或复杂地层，各向异性马尔科夫权系数可以校正地层的各向异性对反演结果的影响，反演的模型参数可以准确反映地层的层状特征，并对断层和边界起到很好的保护效果。

附图说明

图1为本发明方法的流程图。

图2为各向异性马尔科夫随机域权系数的示意图。

图3(a)为合成地震数据测试的纵波速度理论模型；图3(b)为合成地震数据测试的横波速度理论模型；图3(c)为合成地震数据测试的密度速度理论模型。

图4(a)为采用纵波速度模型计算的水平方向各向异性马尔科夫权系数；图4(b)为采用纵波速度模型计算的垂直方向各向异性马尔科夫权系数。

图5(a)为合成地震数据测试的纵波速度反演结果；图5(b)为合成地震数据测试的横波速度反演结果；图5(c)为合成地震数据测试的密度速度反演结果。

图6(a)实际地震资料的纵波速度反演结果；图6(b)实际地震资料的横波速度反演结果。图6(c)实际地震资料的密度反演结果。

具体实施方式

下面结合说明书附图对本发明的实施方式进行描述。

本发明提出的基于各向异性马尔科夫随机域的叠前地震多参数的反演方法，该方法如图1所示，具体包括以下步骤：

步骤一、设置叠前地震参数的目标函数。本反演方法的目标函数包含数据项和先验约束项：

j(vp,vs,ρ)＝j1(vp,vs,ρ)+j2(vp,vs,ρ)(1)

其中，j1为数据项，j2为先验约束项，vp、vs和ρ分别为待反演地层的纵波速度、横波速度和密度参数。

步骤二、使用佐普利兹方程计算地层的纵波反射系数，得到的纵波反射系数与零相位ricker子波褶积后生成合成地震记录。目标函数的数据项为测量所得地震记录与合成地震记录之间误差的二范数，具体表现形式如下：

其中，d为测量所得地震记录；w为震源子波，这里采用零相位ricker子波近似代替；r为地层的纵波反射系数；θ为入射角；t为地震记录的采样时间。

步骤三、计算各向异性马尔科夫随机域权系数。本发明使用各向异性扩散法中的边界终止函数g(▽v)计算层状或复杂地层中的各向异性马尔科夫随机域权系数α，其具体表达形式如下：

其中，▽表示梯度算子，计算所需要的梯度值▽v采用地层的纵波速度参数。q为扩散系数，为预设定的常数。

本发明采用一阶马尔科夫邻域计算各向异性马尔科夫随机域权系数，一阶马尔科夫邻域对应的有8个不同方向的▽v，具体计算形式如下：

其中，下标(i,j)表示纵波速度参数的坐标；下标s和n表示2个垂直方向，e和w表示2个水平方向，wn、ws、en和es表示4个对角方向。

因此，可得到数据点周围8个不同方向的权系数值，如图2所示，即αn，αs，αe，αw，αwn，αws，αen，αns。

步骤四、计算先验约束项，先验约束项由纵波速度、横波速度先和密度先验约束项3部分组成：

j2(vp,vs,ρ)＝j2p(vp)+j2s(vs)+j2d(ρ)(5)

针对层状地层的各向异性，基于各向异性马尔科夫随机域的先验约束项为：

其中，c为数据点邻域内点集，这里使用一阶马尔科夫领域，对于某个数据点，对应与之相邻的8个方向的点集；d表示相邻两个数据点的一阶差分；δp、δs和δd分别为与之对应的正则化参数。此外，先验约束项是通过施加正则化势函数φ实现的，这里采用具有边界保护性质的势函数：

各向异性马尔科夫随机域权系数αc由步骤三得到，用以校正层状地层的各向异性。对层状地层，通常水平方向的权系数大于垂直方向的权系数，确保横向梯度值较小，以保护地层的横向连续性；而层间较大梯度值得以保留，以体现地层的纵向分层性。

步骤五：提取待反演区域的测井数据，对测井数据分别进行反演参数vp和vs及vp和ρ的对数线性拟合，得到对应的拟合系数kc和mc及拟合误差δls和δld，用以控制目标函数极小寻优过程的稳定性，其具体形式如下：

反演参数的对数线性拟合关系式控制步骤六中快速模拟退火过程的稳定性。

步骤六：根据所得数据项和先验约束项代入所建立的叠前地震参数目标函数，确定目标函数及并基于快速模拟退火算法进行目标函数的极小寻优，实现最优化。通过步骤二、步骤三和步骤四建立的目标函数的具体形式为：

使用基于全局优化算法的快速模拟退火算法进行目标函数的极小寻优，快速模拟退火算法结合了三参数的同步反演。此外，对于横波速度的反演，实际反演的参数是δls，最终vs通过由步骤五建立的对数线性拟合关系式得到。

三参数扰动量的具体表现形式为：

lnvs^(m+1)＝k·lnvp^(m+1)+kc+δls^(m+1)(11)

密度参数扰动量需要满足附加条件：

m·lnvp^(m+1)+mc-δld≤lnρ^(m+1)≤m·lnvp^(m+1)+mc+δld(12)

式中vp^(m)、δls^(m)和ρ^(m)是当前模型参数值；vp^(m+1)、δls^(m+1)和ρ^(m+1)是扰动后的参数值；t^(m)是当前的温度值；[vpmin,vpmax]、[δlsmin,δlsmax]和[ρmin,ρmax]是3个参数取值的范围；ξ为分布在[0，1]的随机数；sign(·)为符号函数；k和kc是测井数据中vp和vs的拟合系数；m和mc是测井数据中vp和ρ的拟合系数；δls和δld分别是对应的拟合误差。

步骤七、完成对目标函数的迭代寻优，输出反演结果。

为了验证本发明的方法能够有效解决叠前多参数反演的非线性和不稳定问题，对于层状或复杂地层，各向异性马尔科夫权系数可以校正地层的各向异性对反演结果的影响，本发明列举验证例进行验证说明。

图2为各向异性马尔科夫随机域权系数示意图。以待反演的纵波速度模型的数据点vi,j为例，各权系数αn，αs，αe，αw，αwn，αws，αen，αns对应梯度值的方位如虚线所示，不同灰度的背景表明地层速度的横向连续性和纵向分层性。垂直方向上的权系数值较小，而水平方向上的权系数较大，以校正地层的各项异性对先验约束项的影响。以两个应用实例对本发明进行说明。

实例一，合成地震数据测试：

图3(a)、图3(b)、图3(c)为合成地震数据测试的理论模型，其分别为纵波速度、横波速度和密度三组参数模型。

使用精确佐普里兹方程求解该理论模型的纵波反射系数，与主频为45hz的零相位理论ricker子波进行褶积，获得观测地震记录。观测地震记录有61个角度道集，每个道集有5个角度道(角度间隔为3°，角度覆盖范围为0-15°)。

按公式(2)建立目标函数的数据项。使用精确佐普里兹方程对初始模型进行正演计算，获得合成地震记录。

按公式(3)和(4)计算各向异性马尔科夫随机域的权系数，使用纵波速度模型作为梯度值，扩散系数q取850。除模型边界外，数据点处均可得8个对应不同方向的权系数。

图4(a)和图4(b)分别为水平方向和垂直方向的各向异性马尔科夫随机域权系数值。可以看到，在地层分界面处，垂直方向上的权系数值较小，而水平方向上的权系数较大。

按公式(5)、(6)和(7)建立目标函数的先验约束项。

将15和45道作为虚拟井，并对测井数据按公式(8)进行vp和vs及vp和ρ的对数线性拟合，得到对应的拟合系数和拟合误差。

在本测试中，λ1和λ2的初值为0.3和0.6，vp、δld和密度参数对应的δ的初值分别为250.0、20.0和0.18。

快速模拟退火的初始温度为0.05，终止温度为0.00001，温度衰减系数为0.9，模型参数扰动按公式(10)、(11)和(12)。

对所有数据点依次进行迭代寻优，降温后重复以上步骤，直到达到终止温度，输出最终的反演结果。图5(a)、图5(b)、图5(c)为合成地震数据测试的反演结果，其分别为纵波速度反演结果，横波速度反演结果，密度速度反演结果；可以看到反演结果与理论模型非常吻合，特别是地层的层状特征得到很好保护。

实例二，实际地震资料反演：

实际地震资料中为海上地震资料的一条二维任意线。该任意线共有1981个角度道集，每个道集有15个角度道，角度范围为3–45°，角度间隔为3°，道间距为12.5m，时间采样率为2ms。使用工区内几个主要地层的参数值建立vp和密度的初始模型，δls的初始模型为0。

采用工区内主要地层的层速度计算各向异性马尔科夫权系数。其他参数设置与实例一相同。

图6(a)、图6(b)、图6(c)为实际地震资料的反演结果，其分别为纵波速度反演结果，横波速度反演结果，密度速度反演结果。可以看到，反演结果较好的体现地层的层状特征，目的层的连续性很好。

综上，本发明的方法有效解决叠前多参数反演的非线性和不稳定问题，可以有效的保护反演地层参数的横向连续性和纵向分层性，特别适用于复杂地质条件下的弹性参数反演和储层参数预测。

上面结合附图对本发明的实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下做出各种变化。