一种基于精确Zoeppritz方程的时移地震AVA差异反演方法与流程

本发明涉及油气田、页岩气地震勘探和储层参数预测方法，并且更具体地，涉及一种基于精确zoeppritz方程的时移地震ava差异反演方法。

背景技术：

时移地震油藏监测技术能够测量由于油气开采而引起的地震响应的变化，然后借助地震反演技术对这些变化进行量化，得到储层弹性参数的变化，最后根据岩石物理关系得到储层含油气饱和度、压力以及温度等物性参数的变化，为开发中后期剩余油分布以及储层精细描述提供重要的地震信息，因此该技术被广泛应用于储层表征、监测和管理中。

对于时移地震反演，sarkar等在2003年(sarkars,gouveiawp,johnstondh.2003.ontheinversionoftime-lapseseismicdata.73thannualinternationalmeeting,seg,expandedabstract,1489-1492)通过对比分析指出，不同年份时移地震数据间的反演有必要进行一些耦合，否则将引起最终模型的假象，从而导致不正确的解释结果，差异反演能够提供必要的耦合条件。所以如果基础测量和监测测量数据间的可重复性是合适的，差异反演是最合适的时移地震反演方法。为了进一步提高反演的精度，theune等在2010年(theune,u.,jensasi.o.,andeidsvik,j.2010,analysisofpriormodelsforablockyinversionofseismicavadata:geophysics,vol.75,no.3,25-35)在贝叶斯理论差异反演的基础上，研究了针对时移地震ava(振幅随采集平面的方位角的变化，amplitudevariationwithangle)数据的块反演方法，实际数据处理结果显示该方法能更好地反映时移弹性参数的变化。王守东等在2012年(王守东,王波.时移地震资料贝叶斯avo波形反演.地球物理学报，2012,55卷7期,2422-2431)提出了基于贝叶斯理论的avo波形差异反演方法，但是在方法中正演算子仍然采用的近似公式，导致算子的精度有很大的损失。此外，在先验约束中仅依靠改进的cauchy约束，很难反映差异反演中边界的突变情况，所以该方法在差异反演的精度以及准确性方面有一定的缺陷，需要在先验约束中加入块约束来提高反演精度。基于zoeppritz方程近似公式的反演方法在很多方面限制了叠前反演精度的提高，所以利用精确zoeppritz方程来进行求解是很重要的。此外在拥有良好的多波数据的情况下，联合反演能够进一步增强反演过程的稳定性，提升反演结果的精度。

以上介绍的方法中多是运用分别反演再做差得到差异数据，或者仅是非差异数据的反演，未能实现直接对差异数据进行反演来直接获得三参数的差异数值。并且这些方法中大多是基于zoeppritz方程的近似公式来进行正演算子的构建，这在反演精度上会有一定的降低。

技术实现要素：

针对以上问题，本发明提供了一种基于精确zoeppritz方程的时移地震ava差异反演方法，利用精确zoeppritz方程得到准确度高的正演算子，并运用贝叶斯理论，在先验模型中加入有效突出地层边界效应的块约束，并运用多波联合反演，由此得到高精度、高分辨率的差异反演结果，避免了传统时移地震反演方法分别对开发前、后地震进行反演时引入的人为误差，为直接差异反演提供了技术支持，并最终获得弹性参数的差异数值。

本发明提供的基于精确zoeppritz方程的时移地震ava差异反演方法，包括以下步骤：步骤1，基于该精确zoeppritz方程建立差异数据正演模型；步骤2，基于贝叶斯理论对该差异数据正演模型进行时移地震ava多波联合差异反演，获取油藏弹性参数的变化。

在一个优选的实施方式中，差异数据正演模型以下式表示：

其中，δd为时移地震差异数据，δm为油藏弹性参数的改变量，m1为基础数据对应的弹性参数，l为基于zoeppritz方程的正演算子在基础数据对应的弹性参数m1处的一阶偏导数，e是与差异数据直接相关的误差项。

在一个优选的实施方式中，步骤2可以包括如下子步骤：步骤21，根据步骤1中的该差异数据正演模型确定似然函数；步骤22，构建先验模型，其中，该先验模型包括高斯分布项和服从微分拉普拉斯分布的块约束项；步骤23，根据步骤21中确定的该似然函数以及步骤22中确定的该先验模型，确定反演目标函数；步骤24，对该反演目标函数求解，获得该油藏弹性参数的变化。

在一个优选的实施方式中，步骤21中，该似然函数的负对数以下式表示：

其中，f(δd|δm)为该似然函数的负对数，cd为噪声的协方差矩阵，const为常数项。

在一个优选的实施方式中，步骤22中，该先验模型的负对数以下式表示：

f(δm)＝const+f1(δm)+f2(δm)

其中，f(δm)、f1(δm)和f2(δm)分别表示该先验模型的负对数、高斯分布项的负对数以及块分布项的负对数，cδm为包含三个差异数据相关性的协方差矩阵，μ为差异数据的均值向量，n为模型参数的长度，d为一阶微分算子，kl为尺度因子，l＝1、2或3。

在一个优选的实施方式中，针对多波地震数据，该反演目标函数以下式表示：

α＝σδpp/σδps

β＝σδpp

其中，j(δm)为该反演目标函数，α控制ps波数据的比重，β控制先验信息的比重，σδpp为pp波差异地震数据的噪声方差，σδps为ps波差异地震数据的噪声方差。

在一个优选的实施方式中，在步骤24中，采用迭代重加权最小二乘法进行迭代求解，以获得油藏弹性参数的变化。

在一个优选的实施方式中，该油藏弹性参数包括纵波速度、横波速度以及密度。

本发明基于精确zoeppritz方程构建了时移地震差异数据的正演算子，并由此构建反演目标函数，很好的避免了近似公式带来的计算误差。此外本发明直接利用差异数据进行差异反演，摆脱了分别反演做差带来的计算量的问题，有效的提高了运算效率，有利于实际数据的处理；另外，本发明不仅能够稳定合理的从差异地震数据中准确获取纵波速度、横波速度两个弹性参数变化信息，而且同时能够准确估计密度信息的变化量，验证了新方法的可行性和有效性，为开发中后期储层精细描述和井位部署提供重要的数据。

上述技术特征可以各种适合的方式组合或由等效的技术特征来替代，只要能够达到本发明的目的。

附图说明

在下文中将基于实施例并参考附图来对本发明进行更详细的描述。其中：

图1显示了根据本发明实施例的基于精确zoeppritz方程的时移地震ava差异反演方法的示意性流程图；

图2显示了根据本发明实施例的基于贝叶斯理论进行反演的示意性流程；

图3显示了根据本发明实施例的单井流体替换前后的测井曲线对比示意图；

图4显示了根据本发明实施例的油藏弹性参数的真实差异数值以及初始模型图；

图5显示了根据本发明实施例的采用精确zoeppritz方程进行正演得到的合成叠前角道集；

图6(a)显示了采用现有技术中分别反演然后做差获得弹性参数变化量；

图6(b)显示了根据本发明实施例的反演方法获得的弹性参数变化量；

图6(c)显示了采用aki-richards近似公式的反演方法对差异数据直接反演得到的弹性参数的变化量；

图7显示了根据本发明实施例的含噪时移地震差异反演结果示意图；

图8显示了根据本发明实施例的基于精确zoeppritz方程的时移地震差异数据联合反演结果图。

在附图中，相同的部件使用相同的附图标记。附图并未按照实际的比例。

具体实施方式

下面将结合附图对本发明作进一步说明。

图1显示了本发明提供的基于精确zoeppritz方程的时移地震ava差异反演方法100，如图1所示，该方法100包括以下步骤：

s110，基于所述精确zoeppritz方程建立差异数据正演模型；

s120，基于贝叶斯理论对所述差异数据正演模型进行时移地震ava多波联合差异反演，获取油藏弹性参数的变化。

具体地，针对时移地震待反演参数分布不连续的特征，引入更加合理的先验模型，本发明在假设弹性参数服从高斯分布的同时加入服从微分拉普拉斯分布的垂向块约束项，避免了传统时移地震反演方法分别对开发前、后地震进行反演时引入的人为误差，提升了反演结果的精度，突出了储层弹性参数的变化，为直接差异反演提供了技术支持，并最终获得弹性参数的差异数值。在反演的过程中利用多波联合反演，加入了ps波的信息，使得反演的结果会更加准确。

在步骤s110中，精确zoeppritz方程描述的平面波入射到两个半无限均匀各向同性弹性介质的分界面时的反射和透射情况(zoeppritz,1919；akiandrichards,1980)。其表达式如下：

在时移地震反演中，对于匹配良好的时移地震数据，地震响应的变化通常限定在特定的地层中，而在非油藏部分，基础数据和监测数据应该是相同的。差异数据是将监测数据和基础数据进行相减，这样做可以使得不同时间采集的两次地震数据获得一些必要的耦合，将相减得到的差异数据直接进行反演，进而获取由开采引起的油藏弹性参数的变化。

时移地震基础数据和监测数据的正演过程分别如下：

d1＝g(m1)+n1(2)

d2＝g(m2)+n2(3)

其中，g为zoeppritz方程正演算子。d1为基础数据，d2为监测数据。m1为基础数据对应的弹性参数，m2为监测数据对应的弹性参数，n1和n2分别为两次数据采集时的噪声数据。

将等式(3)在基础数据对应的弹性参数m1处进行泰勒展开，

将等式(4)和等式(2)相减得，

令则有

δd＝lδm+e(6)

其中，δd是差异地震数据，δm是弹性参数的改变量，l(m1)是基于zoeppritz方程的正演算子关于纵横波速度以及密度三个弹性参数在基础数据对应的弹性参数m1处的一阶偏导数，e是与差异数据直接相关的误差项，也可以理解为噪声。

在s120中，基于贝叶斯理论对所述差异数据正演模型进行时移地震ava多波联合差异反演，获取油藏弹性参数的变化。具体地，该步骤可以包括如下的子步骤：

s121，根据步骤s110中的所述差异数据正演模型确定似然函数；

s122，构建先验模型，其中，先验模型包括高斯分布项和服从微分拉普拉斯分布的块约束项；

s123，根据步骤s121中确定的所述似然函数以及步骤s122中确定的所述先验模型，确定反演目标函数；

s124，对所述反演目标函数求解，获得所述油藏弹性参数的变化。

首先，在确定似然函数时，假设等式(6)中的误差项e服从零均值高斯分布，则似然函数可表示为，

上式中，cd是噪声的协方差矩阵，nd是观测数据的长度。对上式取负对数，

时移地震反演通常只是对目的层段弹性参数的变化量进行反演，因此可以认为这一小段的弹性参数是服从正态分布的。由于时移地震差异数据存在分块现象，因此考虑在先验模型中再加入块约束项，并根据theuneuetal(2010)的研究，假设其服从微分拉普拉斯分布。综上所述，先验模型由两部分组成：(1)包含先验低频趋势和不同模型参数之间内部协方差的高斯分布项；(2)具有长尾巴分布特征的块约束项。

假设先验模型服从高斯分布，则有：

对上式取负对数，

其中，cδm为包含三个差异数据相关性的协方差矩阵，n为模型参数的长度，μ为差异数据的均值向量(三个弹性参数需要分别求取)。为了获得剧烈的弹性参数变化，研究在先验模型中加入垂向的块约束，则先验模型的负对数定义为，

f(δm)＝const+f1(δm)+f2(δm)(11)

其中f2(δm)代表块约束项，其表达式如下：

其中，d是一阶微分算子，kl,l＝1,2,3是尺度因子，三个弹性参数之间可能会不一样。

基于贝叶斯理论，约去常数项，可以得到如下所示的反演目标函数：

在实际处理过程中，通常假设观测数据的噪声是不相关的，则噪声的协方差矩阵会简化为一个对角矩阵，即cd＝σn²i，其中σn²表示噪声的方差，i是nd×nd的单位矩阵，nd是观测数据的长度。对于多波地震数据，可以将上述目标函数进行拓展，并假设pp波差异地震数据的噪声方差为σδpp，ps波差异地震数据的噪声方差为σδps，那么上述目标函数可以表示为，

其中α＝σδpp/σδps控制ps波数据的比重，β＝σδpp控制先验信息的比重。

将上述目标函数关于δm求导并令导数为零，

整理得，

其中，

由于上述等式右端项元素a也包含δm，因此采用迭代重加权最小二乘算法进行迭代求解。由于正演算子l(m1)中包含基于精确zoeppritz方程的正演算子关于纵横波速度以及密度三个弹性参数在基础数据对应的储层弹性参数m1处的一阶偏导数，因此需要对已匹配好的基础数据进行基于精确zoeppritz方程的叠前反演来获取正演算子矩阵l。这一项可以借助zhouetal(2016)提出的基于精确zoeppritz方程的非线性avo三参数反演方法得到。

然而，本领域技术人员应当理解的是，基于zoeppritz方程的叠前ava反演算法相比于基于近似公式的反演算法计算量会稍大一些，按照上面所述，需要在一次叠前反演的基础上再进行一次时移地震差异反演，如果处理的实际工区数据较大可能会增大计算量。通过计算对比可以发现，正演算子l在取值m1和m0(利用基础数据单独反演后的反演结果)时计算得到的矩阵元素差异较小，因此为了节省计算量也可以由来近似代替

为了检验该反演算法的可行性和稳定性，申请人在勘探工作中进行了实践。图3中的实线是流体替换前即开发前的测井曲线，虚线是利用gassmann方程对实际油藏区域进行流体替换后的测井曲线，其中，横坐标表示的是波速，纵坐标表示的是时间。图4中实线是将图3所示的开发前后的测井曲线相减得到的真实弹性参数变化量，虚线是对真实变化量进行多次平滑得到的初始模型，其中，横坐标表示的是波速，纵坐标表示的是时间。采用精确zoeppritz方程进行正演得到合成叠前角道集(以4°为间隔，4°到40°之间的10个角度)，如图5所示，其中，横坐标表示的是角度，纵坐标表示的是时间。正演过程中采用的是主频为30hz的雷克子波，采样间隔为2ms。为了检验反演算法的抗噪性，在合成的基础数据和监测数据上分别加入信噪比(均方根振幅比)为2的随机噪声，并利用该算法进行反演。

图6a是利用基于精确zoeppritz的反演方法对合成的基础数据和监测数据分别反演之后相减得到的弹性参数变化量，其中，横坐标表示的是波速，纵坐标表示的是时间(下同)，图6b是利用新方法反演得到的弹性参数变化量，从图6a和图6b的对比可以看出，差异数据的反演结果明显优于分别反演的结果，尤其是横波速度和密度。为了进一步对比分析新方法的优越性，本文在此给出了常规时移地震反演的结果，图6c即是利用基于aki-richards近似公式的反演方法对差异数据直接反演得到的弹性参数的变化量。由于近似公式计算精度有限，且存在诸多假设条件，从图6b和图6c的对比可以看出，新方法的反演结果的精度明显高于常规反演结果的精度。图7是在基础数据和监测数据上分别加入信噪比为2的随机噪声后利用研究方法进行反演得到的反演结果，其中，横坐标表示的是波速，纵坐标表示的是时间，从反演结果可以看出，该方法具有很好的抗噪性。另外，如果存在匹配良好的多波地震数据，利用该方法还可以进行多波联合反演，图8即为基于合成数据的纵横波联合反演结果，其中，横坐标表示的是波速，纵坐标表示的是时间，从图8和图6b的对比可以看出，联合反演能够进一步增强反演过程的稳定性，提升反演结果的精度。

本发明基于精确zoeppritz方程构建了时移地震差异数据的正演算子，并由此构建反演目标函数，很好的避免了近似公式带来的计算误差。此外本发明直接利用差异数据进行差异反演，摆脱了分别反演做差带来的计算量的问题，有效的提高了运算效率，有利于实际数据的处理。

在本发明的描述中，需要理解的是，术语“上”、“下”、“底”、“顶”、“前”、“后”、“内”、“外”、“左”、“右”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。

虽然在本文中参照了特定的实施方式来描述本发明，但是应该理解的是，这些实施例仅仅是本发明的原理和应用的示例。因此应该理解的是，可以对示例性的实施例进行许多修改，并且可以设计出其他的布置，只要不偏离所附权利要求所限定的本发明的精神和范围。应该理解的是，可以通过不同于原始权利要求所描述的方式来结合不同的从属权利要求和本文中所述的特征。还可以理解的是，结合单独实施例所描述的特征可以使用在其他所述实施例中。