本发明属于地震勘探技术领域，尤其涉及一种基于杨氏模量、泊松比的两项式地震叠前反演方法。
背景技术：
当前全球地震勘探技术发展的几个主要趋势包括：从常规油气向常规与非常规油气并重局面发展、从简单地表和浅水区向复杂地表和深水区发展、从叠后地震反演向叠前弹性反演发展。页岩油气的开采需压裂页岩生储层，因此岩区的易压裂程度成为评价优质页岩油气储层非常重要的指标。岩石物理研究表明，杨氏模量和泊松比能够较好的表征岩体的造缝难易程度，标识页岩气的“甜点”区；由于地震测线长度与储层埋深的比值通常很小，故通常很难获得超深储层的大反射倾角叠前地震数据。叠前avo反演可以充分利用叠前地震资料包含的地质、岩性和流体信息，获取丰富的地下介质弹性参数。avo：(amplitudevariationwithoffset)利用地震反射波振幅随偏移距的变化特征来探讨反射系数响应随入射角的变化,进而确定反射界面上覆、下伏介质的岩性特征及物性参数的技术。主要用于石油天然气勘探。叠前avo反演的核心是反射系数方程。通常反射系数方程由三项构成，其中之一为密度项。密度项通常不作为直接用来预测油气储层的指标，原因是具体工区的密度变化范围不大，而各种岩性的理论密度范围交集较广，且其对油气储层特征的直接相关性差，单由密度项无法进行岩性、流体的判断。在缺乏大角度地震资料或者地震资料数据信噪比较低的情况下，大角度反演的结果品质较差，提取参数的过程中，参数矩阵的稳定性较差，造成最终提取的弹性参数误差较大。另外，由于密度项的参数权值在角度有限情况下较小，反演的难度更加大。以ypd三项反射系数方程为基础，构建的ypd叠前地震反演方法，同时反演了杨氏模量、泊松比与密度。该方法为目前较新的杨氏模量、泊松比反演方法(2014年建立)。具体做法如下：假设待反演杨氏模量、泊松比及密度反射系数服从柯西分布，假设似然函数服从高斯分布，同时在反演目标函数中加入初始模型约束，并通过初始模型建立各道去相关矩阵，消除待反演参数间的互相关性，在贝叶斯反演框架下进行反演。该方法在叠前地震数据远角道集缺失的情况下，大角度入射角的引入会使参数矩阵的稳定性降低，造成最终提取的弹性参数存在较大误差。
综上所述，现有技术存在的问题是：针对深层、超深层页岩气储层的地震勘探，由于远角道集数据的缺失，使得利用avo叠前反演较难获取地下介质准确的杨氏模量、泊松比等岩石物理弹性参数。因此，提高参数矩阵稳定性，降低计算带来的累积误差成为亟需解决的问题。
技术实现要素：
针对现有技术存在的问题，本发明提供了一种基于杨氏模量、泊松比的两项式地震叠前反演方法。
本发明是这样实现的，一种基于杨氏模量、泊松比的两项式地震叠前反演方法，所述基于杨氏模量、泊松比的两项式地震叠前反演方法包括以下步骤：
步骤一，叠前地震资料预处理；
步骤二，构造层位解释；
步骤三，测井资料外推与模型构建；
步骤四，多角度子波提取；
步骤五，最终利用模拟退火反演方法获得杨氏模量和泊松比数据体。
进一步，所述基于杨氏模量、泊松比的两项式地震叠前反演方法具体包括以下步骤：
(1)基于aki-richards近似方程推导得到yp两项式反射系数近似方程，作为深层杨氏模量、泊松比叠前反演方法的理论基础；
aki-richards近似方程表述为：
其中，vp、vs分别为纵、横波速度，ρ为密度，θ为入射角。
(2)采用幂指数拟合的方式，引入参数可变的关系式消除密度项，由密度与纵波速度的微商关系：
dρ＝flvpl-1dvp(2)
可得
纵横波模量反射系数与纵横波速度、密度反射系数的关系如下所示：
将(3)带入(4)(5)得
令
(3)整理杨氏模量、泊松比关于纵横波速度的关系，并融合以上各式，得
将(9)、(10)、(11)带入aki-richard方程(1)，得到新的两项式近似方程
其中，k为横、纵波速度比的平方，l为幂指数。
进一步，所述基于杨氏模量、泊松比的两项式地震叠前反演方法在进行杨氏模量与泊松比反演计算以前，需从先验信息中求得参数l、k；在实施操作过程中，其中k为测井资料中横、纵波速度比的平方，l则通过密度与纵波速度的幂指数拟合计算得到。
待反演参数的计算可简化为求解：
其中，ce和cσ与角度有关；re和rσ分别是杨氏模量、泊松比的反射系数时间序列
反射系数与某一角度子波卷积形成一个固定角度的地震合成记录。将式(13)等号两边同时与子波卷积可推导为式
spp(t，θi)是入射角为θi的角道集，w(t，θi)是一个子波核矩阵。
卷积运算之后，可得
若设d＝gm，spp表示包含向量d的叠前角度道集，we和wσ构成小波核矩阵g，re和rσ是包含m的反射率项。
反演的目标函数为
f(v)＝||s-d||→min(21)
式中，v为待反演弹性参数构成的参数向量；s(t)＝w*r为初始低频模型在各角度下的合成地震记录，w为一组不同角度地震子波，r为计算得到的低频地震模型的反射系数；d为实际部分叠加道集地震数据；
构造f(v)的jacobian矩阵j，通过模拟退火算法对j进行求解，得到对v的修正量δv0，将v0＝(v+δv0)作为新的反演参数向量带入式(21)，反复迭代求解f(v)的极小值，直至求得使f(v)达到允许误差的δvn，此时认为vn＝(vn-1+δvn)为最优解；对m的多次迭代最终值即为反演得到的最终杨氏模量、泊松比数据体。
本发明的优点及积极效果为：通过发展两项avo反射系数方程，将密度项以合理的方式与纵波速度建立数学关系，消除密度项，可以在缺少大角度叠前地震资料的情况下进行叠前直接反演，得到目标弹性参数，这样可以减小传统方法带来的累积误差；实际数据来自cgghampson-russell公司开发软件的avo模块的测试数据。对地震资料进行预处理后抽取角道集数据，利用hrs软件依据测井资料对角道集数据提取五组不同角度的子波，并依据所建立的不同子波核矩阵与yp近似方程计算的反射系数合成出不同角度域地震剖面，带入反演流程进行迭代计算；得到最终反演结果的杨氏模量、泊松比反演剖面，反演结果与实际情况较为吻合，显示了本发明较好的应用潜力。
附图说明
图1是本发明实施例提供的基于杨氏模量、泊松比的两项式地震叠前反演方法流程图。
图2是本发明实施例提供的yp近似方程地震叠前反演流程图。
图3是本发明实施例提供的正波阻抗界面反射系数对比示意图；
图中：(a)不同方程反射系数对比；(b)不同方程反射系数差值对比。
图4是本发明实施例提供的负波阻抗界面反射系数对比示意图；
图中：(a)不同方程反射系数对比；(b)不同方程反射系数差值对比。
图5是本发明实施例提供的层状介质模型反演结果示意图；
图中：(a)角度域合成地震记录；(b)杨氏模量、泊松比反演曲线对比。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
下面结合附图对本发明的应用原理作详细的描述。
如图1所示，本发明实施例提供的基于杨氏模量、泊松比的两项式地震叠前反演方法包括以下步骤：
s101：叠前地震资料预处理；
s102：地质层位解释；
s103：测井资料外推与模型构建；
s104：多角度子波提取；
s105：最终利用模拟退火反演方法获得杨氏模量和泊松比数据体。
本发明实施例提供的基于杨氏模量、泊松比的两项式地震叠前反演方法具体包括以下步骤：
首先，基于aki-richards近似方程推导得到yp两项式反射系数近似方程，作为深层杨氏模量、泊松比叠前反演方法的理论基础。
aki-richards近似方程表述为：
其中，vp、vs分别为纵、横波速度，ρ为密度，θ为入射角。
由于可认为密度项是一个与泊松比等信息相关的非独立参数，且相对于其他参数，其变化范围不大，因此可以与纵波速度建立关系。gardner针对砂泥岩建立了密度与纵波速度的数学联系，而后gardner公式成为了将岩体密度拟合纵波速度的经典公式。采用幂指数拟合的方式，引入参数可变的关系式消除密度项，以此保证对实际工区拟合的准确性。
由密度与纵波速度的微商关系
dρ＝flvpl-1dvp(2)
可得
纵横波模量反射系数与纵横波速度、密度反射系数的关系如下所示：
将(3)带入(4)(5)得
令
整理杨氏模量、泊松比关于纵横波速度的关系，并融合以上各式，得
将(9)、(10)、(11)带入aki-richard方程(1)，得到新的两项式近似方程
其中，k为横、纵波速度比的平方，l为幂指数。可将上式称为yp近似方程，该式建立了杨氏模量反射率、泊松比反射率与反射系数的直接线性关系。
本发明开展基于yp近似方程的两项avo反演，其主要流程包含叠前地震资料预处理、地质层位解释、测井资料外推与模型构建、多角度子波提取等环节，最终利用模拟退火反演方法获得杨氏模量和泊松比数据体。主要反演流程如图所示：在进行杨氏模量与泊松比反演计算以前，需从先验信息中求得参数l、k。在实施操作过程中，为确保参数拟合结果准确可靠，针对目的层段可采用测井资料的纵、横波速度与密度确定l、k，其中k为测井资料中横、纵波速度比的平方，l则通过密度与纵波速度的幂指数拟合计算得到。
待反演参数的计算可简化为求解
其中，ce和cσ与角度有关；re和rσ分别是杨氏模量、泊松比的反射系数时间序列
反射系数与某一角度子波卷积形成一个固定角度的地震合成记录。将式(13)等号两边同时与子波卷积可推导为式
spp(t，θi)是入射角为θi的角道集，w(t，θi)是一个子波核矩阵。
卷积运算之后，可得
若设d＝gm，spp表示包含向量d的叠前角度道集，we和wσ构成小波核矩阵g，re和rσ是包含m的反射率项。
反演的目标函数为
f(v)＝||s-d||→min(21)
式中，v为待反演弹性参数构成的参数向量；s(t)＝w*r为初始低频模型在各角度下的合成地震记录，w为一组不同角度地震子波，r为计算得到的低频地震模型的反射系数；d为实际部分叠加道集地震数据；
构造f(v)的jacobian矩阵j，通过模拟退火算法对j进行求解，得到对v的修正量δv0，将v0＝(v+δv0)作为新的反演参数向量带入式(21)，反复迭代求解f(v)的极小值，直至求得使f(v)达到允许误差的δvn，此时认为vn＝(vn-1+δvn)为最优解；对m的多次迭代最终值即为反演得到的最终杨氏模量、泊松比数据体。下面结合具体应用实施例对本发明的应用效果作详细的描述。
位于四川的元坝气田是迄今为止中国发现的埋藏最深的海相气田，选择元坝气田实测井资料构造地质模型对精确zoeppritz方程、aki-richards近似方程与yp近似方程的计算精度进行实验对比分析，模型参数如表1所示。
表1含水饱和度(2014)地震模型
波速与密度的乘积被称为波阻抗，正波阻抗界面的入射介质波阻抗小于投射介质波阻抗，负波阻抗界面则相反。在本模型中，上覆地层与气层的分界面为负波阻抗界面，气层与下伏地层的分界面为正波阻抗界面。分别用精确的zoeppritz方程、aki-richard近似方程、yp近似方程计算上述模型不同界面处的反射系数及近似方程与精确方程的残差。图3和图4分别为正、负波阻抗界面对比分析，图3a、图4a为分别用精确的zoeppritz方程、aki-richard近似方程、yp近似方程计算得到的反射系数随入射角的变化，图3b、图4b为aki-richard近似方程、yp近似方程计算得到的反射系数与精确zoeppritz方程计算得到的反射系数的差值随入射角的变化。
由图3a、图4a可见，基于yp近似方程计算的反射系数与精确zoeppritz方程有一定误差，但精度值稳定优于aki-richard近似方程的计算结果。由图3b、图4b可见，aki-richard近似方程对zoeppritz方程的残差在入射角为50°时开始急剧增大。相对而言，yp近似方程的误差在入射角为60°以内时相当稳定而微小，证明了yp近似方程具有的实际可行性。
为验证yp两项avo反演的可行性与稳定性，本发明采用多种模型进行实验。首先设置一组具有一定地层厚度、共16层的多层水平层状介质模型，赋予其纵、横波速度与密度。子波采用主频为40hz的雷克子波，采用精确zoeppritz方程进行正演得到角度域叠前角道集，采用本发明所述yp近似方程进行杨氏模量、泊松比的直接反演，结果如图5所示。
图5(a)为正演所得角度域叠前角道集，图5(b)蓝色曲线为地层参数计算所得的杨氏模量、泊松比曲线，红色曲线为反演所得杨氏模量、泊松比曲线。由反演结果可见，利用yp近似反演能够获取与真实值基本吻合的杨氏模量和泊松比。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。