CN103245970A - 叠前地震宽角度反演方法 - Google Patents

Abstract

本发明提供一种叠前地震宽角度反演方法，该叠前地震宽角度反演方法包括：输入至少三个角叠加地震数据体、至少三个对应子波以及纵波速度、横波速度以及密度初始模型；建立叠前反演目标函数；将该叠前反演目标函数泰勒级数展开；通过Zoeppritz方程矩阵形式对参量求导，得到反射系数对该纵波速度、该横波速度和该密度的偏导数；根据实际资料，通过矩阵分块减少计算量；求取反演参数摄动量；以及通过求取的该反演参数摄动量控制该反演目标函数的收敛性，获得叠前反演参数。该叠前地震宽角度反演方法解决了现有技术中难以满足Zoeppritz方程成立的近似条件和误差大的问题，具有适用范围广，误差小，反演精度高的优点。

Description

叠前地震宽角度反演方法

技术领域

本发明涉及地球物理勘探技术，特别是涉及到一种叠前地震宽角度反演方法。

背景技术

叠前地震反演可提供多种有效的地层参数数据成果，是复杂油气藏描述的一项重要技术。

Zoeppritz方程描述了波在上下介质传播的反射、透射关系，是叠前地震反演的理论基础。

由于求取Zoeppritz方程精确解的计算量大，且反演系数矩阵是“病态”的原因，实际应用过程中会遇到计算不稳定，目标函数不收敛等难题。

目前的叠前反演方法几乎都是基于Zoeppritz方程的近似式实现的，仅适应于弱反射介质界面、中小角度（或小偏移距）的反射问题，但实际情况往往难以满足公式成立的近似条件。

在小角度范围内只要两层介质的地层参数变化稍大，这种近似在反复迭代计算过程中就会产生一定的计算误差，而且Zoeppritz方程近似式在大角度计算时，误差比较大，难以提高大角度地震资料的叠前反演精度，无法满足勘探开发高精度研究的需求。

为此我们发明了一种新的叠前地震宽角度反演方法，解决了以上技术问题。

发明内容

 本发明的目的是提供一种适用范围广，误差小，反演精度高的叠前地震宽角度反演方法。

    

本发明的目的可通过如下技术措施来实现：该叠前地震宽角度反演方法包括：（1）输入至少三个角叠加地震数据体、至少三个对应子波以及纵波速度、横波速度以及密度初始模型；（2））建立叠前反演目标函数，

                                                 

式中

表示该纵波速度，

表示该横波速度，

为该密度，D为实际角道集记录，

式为地震模型响应，

为使用Zoeppritz方程计算的纵波反射系数，W为地震子波；（3）将该叠前反演目标函数泰勒级数展开；（4）通过Zoeppritz方程矩阵形式对参量求导，得到反射系数对该纵波速度、该横波速度和该密度的偏导数；（5）根据实际资料，通过矩阵分块减少计算量；（6）利用贝叶斯理论，加入岩石物理先验信息，并求取反演参数摄动量；以及（7）通过求取的该反演参数摄动量控制该反演目标函数的收敛性，获得叠前反演参数。

本发明的目的还可通过如下技术措施来实现：

在步骤（4）中，将Zoeppritz方程写为矩阵形式如下：

       

                                      

令

来代替上下两地层岩石物理参数，将该矩阵

两边分别对m中的各分量求偏导数得到该反射系数对该上下地层岩石物理参数的偏导数Jacobi矩阵： 

   

                       

计算

，求解该偏导数Jacobi矩阵得到该反射系数对该纵波速度、该横波速度和该密度的导数。

在步骤（5）中，将偏导数Jacobi矩阵进行划分，将（4×1）维的划分为（2×1）维，将（4×4）维的划分为（2×2）维。

在步骤（3）中，将合成记录进行泰勒展开，并省去二阶高阶项，求取纵横波速度变化量的一阶导数，并令其为0，建立该合成记录与地震记录残差、该合成记录与各参数偏导之间的关系式。

在步骤（3）中，有三个角度的角道集叠加数据，分别为小角度D1、中角度D2和大角度D3，则该关系式的矩阵形式为： 

 

                        

式中

=

，

=

，

=

式中，

和

分别为该纵波速度、该横波速度和该密度的变化量，

。

 在步骤（6）中，建立实际地震记录与模型记录之间的关系得到似然函数，利用贝叶斯公式建立先验概率分布函数，对该先验概率分布函数两边取对数，对模型参数扰动量取导数，并令导数为0，得到参数扰动量表达式，进行迭代运算计算出三个弹性参数的摄动量

。

 在步骤（6）中，建立该先验概率分布函数所用公式如下：

 

 其中，

为模型参数摄动量的方差，

 得到该摄动量的表达式：

该摄动量随迭代次数的变化而变化。

在步骤（7）中，计算参数扰动量与模型相加得到实际地层的弹性参数，对该弹性参数进行参数范围约束、纵横波速度关系约束、纵波速度与密度关系约束，获得叠前反演弹性参数。

在步骤（7）中，在进行该参数范围约束时，三个反演参数的最大最小值用于反演三个参数的变化范围约束，在有井区时，将研究区有井资料统计得到三个弹性参数的最小值的80%和最大值的120%作为反演迭代过程中参数的约束范围，在无井资料或井资料较少时，将沉积岩中三个弹性参数的变化范围作为约束。

 在步骤（7）中，在进行该纵横波速度约束时，实际工区岩石物理分析拟合得到的纵横波速度的关系作为约束条件，对该横波速度进行约束，其中约束公式采用下式：

式中Vp为该纵波速度，Vs为该横波速度，a和b由拟合得到。

在步骤（7）中，在进行该纵波速度与密度关系约束时，采用井资料纵波速度与密度公式拟合，将拟合得到的该纵波速度与该密度的关系作为约束条件对该密度进行约束，其中拟合公式为：

其中ρ为该密度，Vp为该纵波速度，a，b通过井资料拟合得到。

本发明中的叠前地震宽角度反演方法，应用了Zoeppritz方程精确解实现了地震宽角度反演，建立了利用Zoeppritz方程建立对纵、横波速度、介质密度偏导矩阵计算方法，导出Zoeppritz方程与储层有效弹性参数有关的Jacobi偏导数矩阵方程，实现基于Zoeppritz方程反射系数梯度矩阵的精确计算，利用反射系数梯度矩阵的精确解实现叠前地层参数反演的思路。

在实施过程中，引入贝叶斯的方法提高反演的稳定性，并通过优化Jacobi矩阵的结构和引入扰动变量来提高运算效率，确保目标函数的收敛性。

 最终，根据以上内容，建立一套完整的叠前地震宽角度精确解反演方法，从而实现对地层参数更加精确的描述。

本发明中的叠前地震宽角度反演方法，应用Zoeppritz方程推导反射系数对各参数的偏导数Jacobi矩阵，避免了由于近似公式带来的条件限制和近似误差的影响，使反演能够适用于宽入射角度道集，从而更加有效利用地震资料中的大角度信息。

理论研究和实际应用表明，基于Zoeppritz方程的叠前宽角度反演方法能够提高叠前弹性参数的反演精度，尤其是横波速度和密度。

 附图说明

图1为本发明的一实施例中叠前地震宽角度反演方法的流程图；

图2为本发明的一实施例中叠前宽角度纵波速度反演剖面的示意图；

图3为本发明的一实施例中叠前宽角度横波速度剖面的示意图；

图4为本发明的一实施例中叠前宽角度密度剖面的示意图；

图5为本发明的一实施例中A砂组叠前反演纵横波速度比平面图；

图6为本发明的一实施例中加噪音反演结果的示意图；

图7为本发明的一实施例中Zoeppritz方程精确解与Zoeppritz近似方程的解与模型的误差分析对比图。

具体实施方式

为使本发明的上述和其他目的、特征和优点能更明显易懂，下文特举出较佳实施例，并配合所附图式，作详细说明如下。

如图1所示，图1为本发明的叠前地震宽角度反演方法的流程图。

在步骤101，输入三个角叠加地震数据体、三个对应子波以及纵波速度、横波速度以及密度初始模型，角叠加地震数据体角度范围大小视具体情况而定。

在一实施例中，输入的角叠加地震数据体和对应子波的个数为至少三个。

流程进入到步骤102。

在步骤102，建立叠前反演目标函数。

 根据广义线性反演的思路，基于褶积模型和Zoeppritz方程建立如下目标函数：

 式中

表示纵波速度，

表示横波速度，为密度，D为实际角道集记录。

式为地震模型响应。

为使用Zoeppritz方程计算的纵波反射系数，W为地震子波。

 流程进入到步骤103。

在步骤103，将目标函数泰勒级数展开。

将合成记录进行泰勒展开，并省去二阶高阶项，求取纵横波速度变化量的一阶导数，并令其为0，建立起合成记录与地震记录残差与合成记录与各参数偏导之间的关系式。

步骤103中有三个角度的角道集叠加数据，分别为小角度D1、中角度D2和大角度D3，则该关系式的矩阵形式为： 

                        （2）

式中=

，

=

，

=

式中

，

和

分别为纵波速度、横波速度和密度的变化量。

流程进入到步骤104。

在步骤104中，通过Zoeppritz方程矩阵形式对参量求导，得到反射系数对纵波速度、横波速度、密度的导数。

由于Zoeppritz方程形式复杂，通常将其写为矩阵形式如下：

       

                                      （3）

令

来代替上下两层介质的岩石物理参数，将（3）两边分别对m中的各分量求偏导数得： 

       

                       （4）

方程（4）即反射系数对上下地层岩石物理参数的偏导数Jacobi矩阵。

根据实际资料，只要计算出：

，求解方程（4）即可以得到反射系数对各参数的偏导数。

流程进入到步骤105。

在步骤105中，偏导矩阵分块优化。

步骤104中由于方程组的数量十分庞大，每次迭代都需要求解这些方程组，必然带来很大的计算量，而且方程组的稳定性也给求解方程带来很大的麻烦，对此分析矩阵方程，对矩阵进行必要的简化，减少求解方程的计算量，提高反演的稳定性。

由于参数间存在着不关联性，在Zoeppritz方程偏导数方程系数矩阵中存在着大量的0值，如果根据这些0值得分布特点巧妙地避开大量的0+x或0*x等无效运算，对Zoeppritz偏导方程系数矩阵进行划分，将（4×1）维的划分为（2×1）维的，将（4×4）维的划分为（2×2）维的，可以节省大量的运算时间。

流程进入到步骤106。

在步骤106，利用贝叶斯理论，加入岩石物理先验信息，并求取反演参数摄动量。

 建立实际地震记录与模型记录之间的关系得到似然函数，利用贝叶斯公式建立先验概率分布函数，对概率分布函数两边取对数，对模型参数扰动量取导数，并令导数为0，可得参数扰动量表达式，进行迭代运算可以完整的计算出三个弹性参数的摄动量

。   

步骤106建立先验概率分布函数所用公式如下：

其中，

为模型参数摄动量的方差，

步骤106中得到摄动量的表达式：

 

该摄动量随迭代次数的变化而变化，可以有效提高反演运算效率和收敛效果，在一定程度上能够克服因地震资料噪音带来的反演结果不稳定。

 流程进入到步骤107。

在步骤107，通过求取的该反演参数摄动量控制目标函数的收敛性，最后获得叠前反演参数。

计算参数扰动量与模型相加得到实际地层的弹性参数，对弹性参数进行参数范围约束、纵横波速度关系约束、纵波速度与密度关系约束，最终获得叠前反演弹性参数。

步骤107中反演参数范围约束：三个反演参数的最大最小值用于反演三个参数的变化范围约束。

在有井区，将研究区有井资料统计得到三个弹性参数的最小值的80%和最大值的120%作为反演迭代过程中参数的约束范围。

若无井资料或井资料较少，可根据前人经验，将沉积岩中三个弹性参数的变化范围作为约束。

步骤107中纵横波速度约束：

实际工区岩石物理分析拟合得到的纵横波速度的关系作为约束条件，对横波速度进行约束。

其中约束公式采用下式：

 

式中Vp为纵波速度，Vs为横波速度，a和b由拟合得到

步骤107中纵波速度与密度关系约束：采用井资料纵波速度与密度公式拟合，将拟合得到的纵波速度与密度的关系作为约束条件对密度进行约束。

其中拟合公式为：

其中ρ为密度，Vp为纵波速度，a，b通过井资料拟合得到。

流程结束。

从图1中的流程图可以看出，经过zeoppritz方程精确解建立的反射系数与地层两侧纵波速度、横波速度、密度的关系，避免了采用近似式带来的条件限制和近似误差的影响，精度高。

并且可以更加有效利用地震资料中的大角度信息，得到叠前弹性参数，更加准确。

本发明所述的方法，可得到更高精度的叠前弹性参数，对于进行储层分布规律研究尤其是含流体预测，更加可靠。

从图2到图4可以看出，反演得到的纵波速度、横波速度以及密度的分辨率较高，能够较好的识别薄储层，而且能够很好的刻画地质沉积特征。

图中层位断层明显，断点清晰，能够很好的刻画砂体特征。

纵向上基本可以分辨出五米以上的储层，且与测井吻合较好，能够反映岩性变化，横向上砂体展布清楚，变化自然。

如图5为本发明的一实施例中A砂组叠前反演纵横波速度比平面图，图中左上角的低值区为含油气有利区域，低值区域很好地反应了河流相的储层展布特征。

油井分布在低值区，水井分布在高值区，预测结果与钻井动态资料吻合较好。

图6为本发明的一实施例中加噪音反演结果的示意图，左图为纵波速度反演结果及精度对比图，中图为横波速度反演结果及精度对比图，右图为密度反演结果及精度对比图，其中模型为方波形曲线，另外两条曲线为反演结果，其中粗线为信噪比5:1时的反演结果，细线为信噪比1:1时的反演结果。

从图6可以看出，在信噪比较高时（5:1），反演结果与理论模型吻合率较高，即使是薄层的相对误差变大，但对由岩性变化引起的参数相对增减的表征依然准确；当信噪比降低至1:1时，反演结果出现抖动，但是总体变化规律与模型仍然保持一致。

证明本方法具有较好的抗噪性。

图7为本发明的一实施例中Zoeppritz方程精确解与Zoeppritz近似方程的解与模型的误差分析对比图。

在图7中，左图为纵波速度反演结果及精度对比图，中图为横波速度反演结果及精度对比图，右图为密度反演结果及精度对比图，其中模型为方波形曲线，另外两条曲线为反演结果，粗线为常规反演结果，细线为本方法反演结果。

从对比结果看，采用本方法求得的纵波速度、横波速度和密度比近似解的结果有10%~15%的提高。 

Claims (10)

1.叠前地震宽角度反演方法, 其特征在于，该叠前地震宽角度反演方法包括： 

（1）输入至少三个角叠加地震数据体、至少三个对应子波以及纵波速度、横波速度以及密度初始模型；

（2）建立叠前反演目标函数，

        

式中表示该纵波速度，

表示该横波速度，为该密度，D为实际角道集记录，

式为地震模型响应，

为使用Zoeppritz方程计算的纵波反射系数，W为地震子波；

（3）将该叠前反演目标函数泰勒级数展开；

（4）通过Zoeppritz方程矩阵形式对参量求导，得到反射系数对该纵波速度、该横波速度和该密度的偏导数；

（5）根据实际资料，通过矩阵分块减少计算量；

（6）利用贝叶斯理论，加入岩石物理先验信息，并求取反演参数摄动量；以及

（7）通过求取的该反演参数摄动量控制该反演目标函数的收敛性，获得叠前反演参数。

2.根据权利要求1所述的叠前地震宽角度反演方法，其特征在于，在步骤（4）中，将Zoeppritz方程写为矩阵形式如下：

       

                                      

令

来代替上下两地层岩石物理参数，将该矩阵两边分别对m中的各分量求偏导数得到该反射系数对该上下地层岩石物理参数的偏导数Jacobi矩阵： 

                              

计算

，求解该偏导数Jacobi矩阵得到该反射系数对该纵波速度、该横波速度和该密度的导数。

3.根据权利要求2所述的叠前地震宽角度反演方法，其特征在于，在步骤（5）中，将偏导数Jacobi矩阵进行划分，将（4×1）维的划分为（2×1）维，将（4×4）维的划分为（2×2）维。

4.根据权利要求1所述的叠前地震宽角度反演方法，其特征在于，在步骤（3）中，将合成记录进行泰勒展开，并省去二阶高阶项，求取纵横波速度变化量的一阶导数，并令其为0，建立该合成记录与地震记录残差、该合成记录与各参数偏导之间的关系式。

5.根据权利要求4所述的叠前地震宽角度反演方法，其特征在于，在步骤（3）中，有三个角度的角道集叠加数据，分别为小角度D1、中角度D2和大角度D3，则该关系式的矩阵形式为： 

                        

式中

=

，

=

，

=

式中

，

和

分别为该纵波速度、该横波速度和该密度的变化量，

。

6.根据权利要求1所述的叠前地震宽角度反演方法，其特征在于，在步骤（6）中，建立实际地震记录与模型记录之间的关系得到似然函数，利用贝叶斯公式建立先验概率分布函数，对该先验概率分布函数两边取对数，对模型参数扰动量取导数，并令导数为0，得到参数扰动量表达式，进行迭代运算计算出三个弹性参数的摄动量

和

。

7.根据权利要求6所述的叠前地震宽角度反演方法，其特征在于，在步骤（6）中，建立该先验概率分布函数所用公式如下：

    

其中，

为模型参数摄动量的方差，

 得到该摄动量的表达式：

 

该摄动量随迭代次数的变化而变化。

8. 根据权利要求1所述的叠前地震宽角度反演方法，其特征在于，在步骤（7）中，计算参数扰动量与模型相加得到实际地层的弹性参数，对该弹性参数进行参数范围约束、纵横波速度关系约束、纵波速度与密度关系约束，获得叠前反演弹性参数。

9.根据权利要求8所述的叠前地震宽角度反演方法，其特征在于，在步骤（7）中，在进行该参数范围约束时，三个反演参数的最大最小值用于反演三个参数的变化范围约束，在有井区时，将研究区有井资料统计得到三个弹性参数的最小值的80%和最大值的120%作为反演迭代过程中参数的约束范围，在无井资料或井资料较少时，将沉积岩中三个弹性参数的变化范围作为约束。

10.根据权利要求8所述的叠前地震宽角度反演方法，其特征在于，在步骤（7）中，在进行该纵横波速度约束时，实际工区岩石物理分析拟合得到的纵横波速度的关系作为约束条件，对该横波速度进行约束，其中约束公式采用下式：

式中Vp为该纵波速度，Vs为该横波速度，a和b由拟合得到。

Images