一种基于同伦延拓的叠前AVO三参数反演方法和系统与流程

本发明涉及一种基于同伦延拓的叠前avo三参数反演方法和系统，属于地震勘探技术领域。

背景技术：

avo反演是一种估计地下介质弹性参数的叠前反演技术。根据zoeppritz方程，叠前道集中同一反射点的反射振幅和入射角及反射界面上下介质的弹性参数性质有关，avo反演正是根据这一关系，利用反演方法由叠前道集计算地下介质的弹性参数。zoeppritz方程是一个非常复杂的公式，大多数avo反演方法都是在zoeppritz方程的各种近似公式的基础上进行的。考虑到avo反演的不稳定性，许多学者通过降维的方法来提高反演过程的稳定性，ursenbach等人对双参数avo反演的各种方法进行了总结，证明了各种双参数反演方法具有等价的信息量，反演结果利用适当的公式可以进行互相转换。考虑到密度信息在流体预测中的重要作用，许多学者开展了三参数反演方法研究。buland、印兴耀、陈建江等人利用贝叶斯理论研究了三参数avo反演方法，用cauchy分布描述先验模型参数的分布，利用参数协方差矩阵对反演过程进行约束，利用先验地质信息作为avo反演的约束条件，提高了反演过程的稳定性。考虑到avo反演过程的局部极小问题，misra等人研究了avo反演的全局优化方法，他们在快速模拟退火方法的基础上提出了一种混合全局优化方法，通过一种保边界的光滑预处理方法将先验信息引入到反演过程中，提高反演过程的稳定性。kuzma等人研究了基于支持向量机的avo反演方法，利用已知的模型数据和观测数据对支持向量机进行训练，得到近似的反演过程关系式，该方法具有计算速度快的优点。hennenfent等人将曲波变换和小波变换引入到avo反演过程中，通过曲波变换和小波变换提高反演过程的稳定性。

虽然目前已经有了一些有效的三参数avo波形反演方法，但是现有技术中的反演方法并没有很好的解决avo三参数反演中局部极小的问题，在反演过程中的具体表现是，当给定的初始模型与准确结果差异较小时，反演结果准确性较高，但是当给定的初始模型与准确结果相差较大时，往往无法得到正确的反演结果，这就使得反演过程中初始模型对反演结果非常重要。而在实际反演时，初始模型往往与准确结果差异较大，不能满足要求。

技术实现要素：

针对上述问题，本发明的目的是提供一种基于同伦延拓的叠前avo三参数反演方法和系统，其将同伦延拓算法引入到叠前avo三参数反演过程中，拓宽了反演的收敛范围，降低了计算结果对初始模型的依赖性，提高了反演结果的可靠性。

为实现上述目的，本发明采取以下技术方案：一种基于同伦延拓的叠前avo三参数反演方法，包括以下步骤：s1通过地震数据获得角度域共成像点道集，并根据角度域共成像点道集获得初始参数模型；s2读取角度域共成像点道集中的一个共成像点道集s^*，将共成像点道集s^*代入初始参数模型，获得初始参数向量p⁰；s3采用同伦延拓算法获得参数向量的迭代公式，使得地震记录与角度域共成像点道集中的地震记录的差最小；s4将初始参数向量p⁰代入迭代公式，获得最终的参数向量pⁿ，n为迭代次数。

进一步，角道集地震记录公式为：

s＝wq(p)

其中，w为由地震子波形成的矩阵，q(p)为反射系数向量。

进一步，步骤s3中地震记录与角度域共成像点道集中的地震记录的差最小，通过求解以下泛函o(p)的极小值获得：

其中，s(p)是角度域共成像点道集中的地震记录，s^*是已知的叠前地震记录，mp是参数向量的低频分量，p^l是其它资料得到的低频分量，α和β分别是约束因子和正则因子，m和r均为二阶常数矩阵。

进一步，泛函o(p)的极小值通过使泛函o(p)对p的导数等于零，并通过同伦映射h(p,t)对其求解，从而获得最终的迭代公式；其中，同伦映射的公式为：

其中，t为时间，t为转置矩阵。

进一步，迭代公式为：

其中，m和r为常数矩阵，p为参数向量，p^k为第k次迭代时对应的参数向量，α为约束因子，β为正则因子。

进一步，的计算方法为：将初始参数向量p⁰带入角道集地震记录公式获得初始地震记录s(p⁰)，对初始地震记录s(p⁰)求导，得到的计算方法为：将第k次迭代时对应的参数向量p^k带入角道集地震记录公式获得第k次迭代时的地震记录s(p^k)，对第k次迭代时的地震记录s(p^k)求导，得到

进一步，参数向量p包括和ρ^j，和ρ^j分别为j时刻的纵波速度、横波速度和介质密度，j＝1,2,……,nt-1，nt是记录时间。

进一步，步骤s4结束后，以最终的参数向量pⁿ为新的迭代初始值，利用gauss-newton法对迭代方程的解进行精细化。

本发明还公开了一种基于同伦延拓的叠前avo三参数反演系统，包括：初始参数模型建立模块，用于通过地震数据获得角度域共成像点道集，并根据角度域共成像点道集获得初始参数模型；初始参数向量获取模块，用于读取角度域共成像点道集中的一个共成像点道集s^*，将共成像点道集s^*代入初始参数模型，获得初始参数向量p0；迭代公式建立模块，用于采用同伦延拓算法获得参数向量的迭代公式，使得地震记录与角度域共成像点道集中的地震记录的差最小；输出模块，用于将初始参数向量p0代入迭代公式，获得最终的参数向量pⁿ，n为迭代次数。

进一步，迭代公式为：

其中，m和r为常数矩阵，p为参数向量，p^k为第k次迭代时对应的参数向量，α为约束因子，β为正则因子。

本发明由于采取以上技术方案，其具有以下优点：1、本发明利用精确的zoeppritz方程描述叠前地震记录的生成过程，具有较高的精度，能够很好的描述地震波在地下介质中的传播过程。2、本发明利用地震资料中的叠前角度道集数据开展反演工作，能充分考虑多层介质情况下的地震波的耦合问题，反演结果具有更高的精度。3、本发明可以根据叠前道集数据同时反演出三个弹性参数(纵波速度、横波速度及密度)，三个弹性参数比两个弹性参数在储层预测中具有优势，尤其是密度参数，往往能够更好地预测储层和油气的分布情况。4、本发明在avo反演过程中，引入了同伦延拓算法，扩大了反演的收敛范围，提高了反演结果的可靠性。

附图说明

图1是本发明的基于同伦延拓的叠前avo三参数反演方法的结构示意图；

图2是本发明一实施例中单道模型的实测参数向量，图2(a)、(b)、(c)分别为该单道模型的实测纵波速度、横波速度和密度的折线图；

图3是本发明一实施例中单道模型通过反演获得的参数向量，图3(a)、(b)、(c)分别为该单道模型反演获得的纵波速度、横波速度和密度的曲线图；

图4是本发明一实施例中二维模型的纵波速度，图4(a)是实测的纵波速度，图4(b)是通过反演获得的纵波速度；

图5是本发明一实施例中二维模型的横波速度，图5(a)是实测的横波速度，图5(b)是通过反演获得的横波速度；

图6是本发明一实施例中二维模型的密度，图6(a)是实测的横波速度，图6(b)是通过反演获得的横波速度。

具体实施方式

为了使本领域技术人员更好的理解本发明的技术方向，通过具体实施例对本发明进行详细的描绘。然而应当理解，具体实施方式的提供仅为了更好地理解本发明，它们不应该理解成对本发明的限制。在本发明的描述中，需要理解的是，所用到的术语仅仅是用于描述的目的，而不能理解为指示或暗示相对重要性。

实施例一

本实施例涉及一种基于同伦延拓的叠前avo三参数反演方法，如图1所示，包括以下步骤：s1通过地震数据获得角度域共成像点道集，并根据角度域共成像点道集获得初始参数模型；s2读取角度域共成像点道集中的一个共成像点道集s^*，将共成像点道集s^*代入初始参数模型，获得初始参数向量p⁰；s3采用同伦延拓算法获得参数向量的迭代公式，使得地震记录与角度域共成像点道集中的地震记录的差最小；s4将初始参数向量p⁰代入迭代公式，获得最终的参数向量pⁿ，n为迭代次数，在本实施例中迭代次数优选为30次。

叠前avo三参数反演方法的理论基础是zoeppritz方程。在各向同性弹性介质中，当一个平面纵波入射到两种介质的分界面上，就要产生反射波和透射波。在界面上，根据应力连续性和位移连续性，并引入反射系数、透射系数，就可以得出相应波的位移振幅方程，即zoeppritz方程。假设上层介质的纵波速度、横波速度和密度分别为vp1、vs1和ρ1，第二层介质的纵波速度、横波速度和密度分别为vp2、vs2和ρ2，入射纵波、反射横波、透射纵波、透射横波的角度分别为θp1、θs1、θp2、θs2，则纵波反射系数qpp、横波反射系数qps、纵波透射系数tpp、横波透射系数tps满足zoeppritz方程：

求解zoeppritz方程可得到纵波、横波的反射系数和透射系数的表达式，其中，纵波的反射系数qpp为：

其中，

射线参数u满足snell定律：而

纵波的反射系数qpp公式给出了单个界面的反射系数的计算方法，对于多层介质，每一个反射界面都用上述公式计算反射系数，从而得到叠前角道集反射系数q(θ,t)，将其与地震子波w(t)褶积即可得出叠前角道集地震记录s(θ,t)的计算公式：

s(θ,t)＝w(t)*q(θ,t)

下面讨论叠前角道集地震记录s(θ,t)的计算公式的离散形式，首先假设，角度θ离散为θi,i＝1,2,…,nθ，时间t离散为tj,j＝1,2,…,nt，nθ是离散角度的个数，nt是记录时间。j时刻的离散的纵波速度、横波速度和密度分别为ρ^j,j＝1,2,…,nt-1。反射系数q(θ,t)离散为：

qi,j＝q(θi,tj),i＝1,2,…,nθ,j＝1,2,…,nt-2

其中，qi,j是ρ^j、ρ^j+1的非线性函数，将所有qi,j排列成向量q，向量q即为反射系数向量，将ρ^j排成向量p，则q是p的非线性函数：

q＝q(p)

将地震记录s(θ,t)离散为：

si,j＝s(θi,tj),i＝1,2,…,nθ,j＝1,2,…,nt-2

再将所有地震记录si,j排列成向量s，向量s即为角道集地震记录向量，将其与地震子波形成的矩阵w褶积即可得出角道集地震记录向量的公式为：

s＝wq(p)

利用上述公式可以由参数向量p计算出角道集地震记录向量，实现叠前avo三参数正演模拟。

叠前avo三参数反演模拟与上述过程相反，是根据已知的叠前地震记录s^*，通过反演的方法求解参数向量p，使得计算出的地震记录s(p)与已知的叠前地震记录s^*的差最小。即地震记录s(p)与已知的叠前地震记录s^*可以通过下式表示：

考虑到avo三参数反演的不稳定性，需要在上述公式中加入低频约束项和正则化项，则上述公式被修正为：

这里m、r是常数矩阵。

其中，mp是参数向量的低频分量，p^l是其它资料得到的低频分量，α和β分别是约束因子和正则因子。

下面基于同伦延拓算法计算上述经过修正的方程的极小值，若上述方程在p点达到极小值，则必须满足其导数为零，即：

将上述方程代入得：

由于是关于p的非线性函数，所以上式为一个非线性方程组，该非线性方程组通过同伦法进行求解。取初值p⁰为其它资料得到的低频分量p^l，并假设mp^l＝p^l，则构造同伦映射，即可将上式写成：

可以验证p⁰为h(p,0)＝0的解，而h(p,1)＝0就是经过同伦映射前的方程，利用同伦的方法对上式进行求解。首先将t在[0,1]上离散，即可以写成：

其中，n为迭代次数。p^k为第k次迭代时对应的参数向量。为了避免二阶导数的计算，令

则经过同伦映射后获得的方程可以写为：

将上式两边关于t求导，得：

通过求解上述微分方程得到参数向量的迭代公式：

其中，m和r为常数矩阵，p为参数向量，p^k为第k次迭代时对应的参数向量，α为约束因子，β为正则因子。

进一步，在步骤s4结束后，以最终的参数向量pⁿ为新的迭代初始值，利用gauss-newton法对迭代方程的解进行精细化。

实施例二

基于相同的发明构思，本实施例公开了一种基于同伦延拓的叠前avo三参数反演系统，包括：

初始参数模型建立模块，用于通过地震数据获得角度域共成像点道集，并根据角度域共成像点道集获得初始参数模型；

初始参数向量获取模块，用于读取角度域共成像点道集中的一个共成像点道集s^*，将共成像点道集s^*代入初始参数模型，获得初始参数向量p⁰；

迭代公式建立模块，用于采用同伦延拓算法获得参数向量的迭代公式，使得地震记录与角度域共成像点道集中的地震记录的差最小；输出模块，用于将初始参数向量p⁰代入迭代公式，获得最终的参数向量pⁿ，n为迭代次数。

其中，迭代公式为：

在上式中，m和r为常数矩阵，p为参数向量，p^k为第k次迭代时对应的参数向量，α为约束因子，β为正则因子。

实施例三

为了验证本发明中的反演方法是否可以有效的对参数向量进行反演，本实施例引入了一个单道模型。如图2所示，图2(a)、(b)、(c)分别为该单道模型的实测纵波速度、横波速度和密度的折线图，图3(a)、(b)、(c)中折线是与图2相同的该单道模型的实测纵波速度、横波速度和密度的折线图，图3(a)、(b)、(c)中黑色折线是实际测得的实测纵波速度、横波速度和密度图。灰色曲线是采用本发明中反演方法获得的实测纵波速度、横波速度和密度的曲线图，可以看出图3(a)、(b)、(c)中曲线与折线基本重合，只有比较少的部分存在较小的偏差，说明对于一个单道的模型，本发明中方法能够较为准确的反演出参数向量，反演出的参数向量基本与实测值相符。

实施例四

实施例三中以证明本发明中反演方法在一个单道的模型中能够比较准确的反演出参数向量，但其结构较为简单，故本实施例引入一个第二维模型，图4(a)、图5(a)、图6(a)分别为该二维模型中实际测量的地下介质的纵波速度、横波速度及密度，作为对比，图4(b)、图5(b)、图6(b)分别为根据本发明中反演方法获得的地下介质的纵波速度、横波速度及密度。通过对比可知两组图片形状和颜色都基本相同，可见对于二维模型，本发明中方法也能够较为准确的反演出参数向量。

最后应当说明的是：以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详细的说明，所属领域的普通技术人员应当理解：依然可以对本发明的具体实施方式进行修改或者等同替换，而未脱离本发明精神和范围的任何修改或者等同替换，其均应涵盖在本发明的权利要求保护范围之内。上述内容仅为本申请的具体实施方式，但本申请的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本申请揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本申请的保护范围之内。因此，本申请的保护范围应以权利要求的保护范围。