CN104597490A - 基于精确Zoeppritz方程的多波AVO储层弹性参数反演方法 - Google Patents

Abstract

本发明涉及基于精确Zoeppritz方程的多波AVO储层弹性参数反演方法，所述方法包括a.确定振幅缩放因子；b.统计模型参数的先验信息，并计算包含三参数相关性的协方差矩阵；c.建立时间域的初始弹性参数模型；d.基于所述弹性参数模型和精确Zoeppritz方程正演模拟角度域的多波地震叠前道集，由正演模拟记录和实际记录计算PP波与PS波反演残差；e.构建反演目标函数，对所述反演目标函数求解，得到模型参数扰动量的求解表达式；f.根据所述反演残差和所述扰动量求解表达式计算模型参数的扰动量；g.重复迭代步骤d、e、f获得最优模型参数。本方法有效提高现有叠前AVO反演方法的计算精度与稳定性，满足地震叠前反演识别油气储层，特别是页岩气储层表征的要求。

Description

基于精确Zoeppritz方程的多波AVO储层弹性参数反演方法

技术领域

本发明涉及油气田、页岩气地震勘探和储层参数预测方法，特别是一种基于精确Zoeppritz方程的多波AVO储层弹性参数反演方法。

背景技术

利用有限的测井数据，对储层进行空间上的连续描述是非常困难的。基于井数据插值、外推或模拟得到储层的3D空间模型精度取决于储层周围井的数量和空间分布。当储层周围井的数量有限时，基于测井的储层描述横向分辨率低，难于实现井间储层的有效表征。地震数据具有较好的空间连续性，在储层表征中起着非常重要的作用。高精度反演方法使反演结果具有重要价值，可以提高地震数据解释的可信度。利用反演结果可以更好地实现储层空间描述，很大程度地提高储层表征的价值。叠后地震反演对规模相对较大，砂泥岩阻抗差异比较大的储层表征精度高。面对空间展布范围小，具有复杂结构的储层，具有先天的缺陷的叠后反演显得力不从心。由于AVO(AmplitudeVersus Offset，振幅随偏移距的变化)角度道集保留了叠前不同入射角的反射信息，对岩性和流体的区分度高，因此叠前AVO反演可以得到纵波阻抗、横波阻抗和密度等弹性参数信息，而联合三种地震弹性参数可以更好地区分储层的岩性和流体，最终达到通过AVO反演进行储层表征的目的。由于精确Zoeppritz方程复杂，目前AVO反演主要基于Zoeppritz方程近似公式利用叠前道集上振幅随偏移距变化的信息来提取弹性参数，但近似公式不仅不能直接对纵波速度、横波速度、密度三个弹性参数直接反演，而且假设条件太多，尤其是在入射角(或偏移距)较大时，存在较大的计算误差，在计算精度上不能满足实际油藏储层精细表征要求。而全波形反演虽然能够利用全波场的信息来预测模型参数，但其计算量巨大，在反演尺度和计算效率上不能满足实际油藏储层精细表征要求，尤其对于实际三维大偏移距多波多分量地震数据。

综上所述，目前基于叠前地震数据的储层弹性参数反演方法研究存在以下问题：基于Zoeppritz方程近似公式的AVO反演无法提取大偏移距信息，适应性较差；油气藏储层特别是页岩气储层最常见的结构是泥砂岩薄互层结构，纵向弹性参数差异较大，极大地降低了Zoeppritz近似公式的求解精度，基于Zoeppritz方程近似表达式的AVO反演方法不再适用于这类储层；目前常规的AVO反演方法都不直接对纵波速度、横波速度、密度等弹性参数进行反演，而是间接地求解得到，降低了反演的精度；传统的叠前AVO反演方法一般只考虑PP波数据，PP波数据主要包含纵波速度信息，但是对模型速度和密度不敏感，单独利用PP地震AVO数据反演横波速度信息和密度具有较大的不确定性；如何有效获得地下介质弹性参数模型信息作为先验信息，以降低反演的不确定性，提高反演的精度，有待进一步研究。

现有技术中，如公开号CN104237936A的专利文献公布了一种油气检测的频变反演方法，具体包括：(1)输入叠前地震记录；(2)根据步骤(1)输入的叠前地震记录生成叠前角道集；(3)对步骤(2)得到的叠前角道集进行频谱分解获得分频角道集记录；(4)利用步骤(3)的分频角道集记录以及分频角道集记录对应的频率进行频变AVO反演，获得纵波频散梯度；(5)利用步骤(4)得到的纵波频散梯度预测油气储层。其所述算法不直接对纵波速度、横波速度、密度等弹性参数进行反演，而是间接地求解得到，降低了反演的精度；只考虑PP波数据，PP波数据主要包含纵波速度信息，但是对模型速度和密度不敏感。

发明内容

针对背景技术中出现的问题，本发明提出一种基于精确Zoeppritz方程的多波AVO储层弹性参数反演方法，包括以下步骤：

步骤a.确定振幅缩放因子；

步骤b.统计模型参数的先验信息，并计算包含三参数相关性的协方差矩阵；

步骤c.建立时间域的初始弹性参数模型；

步骤d.基于所述弹性参数模型和精确Zoeppritz方程正演模拟角度域的多波地震叠前道集，由正演模拟记录和实际记录计算PP波与PS波反演残差；

步骤e.构建反演目标函数，对所述反演目标函数求解，得到模型参数扰动量的求解表达式；

步骤f.根据所述反演残差和所述扰动量求解表达式计算模型参数的扰动量；

步骤g.重复迭代步骤d、e、f获得最优模型参数。

优选的是，所述步骤a确定振幅缩放因子包括：

基于地震数据提取角度依赖的子波；

基于测井数据和Zoeppritz正演模拟地震角道集并结合实际井旁地震数据确定振幅缩放因子。

在上述任一方案中优选的是，所述步骤b包括：

采用三变量修正Cauchy分布函数作为先验分布函数；

基于工区内所有测井数据统计得到各模型参数；

求取所述各模型参数的自相关系数和互相关系数，构建三参数的协方差矩阵，形成符合该工区的模型参数先验分布函数。

在上述任一方案中优选的是，所述先验信息包括：纵波速度、横波速度和密度。

在上述任一方案中优选的是，所述步骤c建立初始弹性参数模型包括：

利用地震构造解释资料，基于沉积模式建立地质模型；

将测井资料按构造模式进行插值和外推，得到每条测线的初始弹性参数模型。

在上述任一方案中优选的是，所述弹性参数模型包括：纵波速度模型、横波速度模型和密度模型。

在上述任一方案中优选的是，所述步骤d包括：

以时间域的初始模型作为输入，利用精确Zoeppritz方程计算每个角度的PP波与PS波的反射系数向量；

将子波与反射系数进行褶积，获得角度域的PP波与PS波角道集；

将模拟的角道集数据与实际角度域地震道集作差，获得多波联合反演残差。

在上述任一方案中优选的是，所述步骤e包括：

根据Bayesian原理，综合反演似然函数和先验分布函数得到后验概率分布函数；

根据后验概率确定反演目标函数；

借助于泰勒展开对目标函数进行简化，得到关于模型参数的扰动量函数；

将目标函数对扰动量求导并令导数等于零，得到扰动量的迭代求解公式。

在上述任一方案中优选的是，所述对反演目标函数的求解方法包括：高斯-牛顿法。

在上述任一方案中优选的是，所述步骤f中计算模型参数的扰动量的方法包括：迭代加权最小二乘法。

在上述任一方案中优选的是，所述步骤g中迭代次数通过反演残差控制。

在上述任一方案中优选的是，当地震数据信噪比较高时，使用统一的振幅缩放因子；当信噪比低时，分近、中、远偏移距分别计算振幅缩放因子。

在上述任一方案中优选的是，所述步骤c中建立弹性模型包括：

利用散点插值对各个层位的数据进行插值，完成地质层位建模；

根据地质层位进行弹性参数横向插值，即将测井信息进行横向插值，计算得到地下每个点上的弹性参数值，完成初始弹性参数建模的任务。

在上述任一方案中优选的是，基于精确Zoeppritz方程的正演过程表示如下：

$\left\{\begin{array}{c}{G}_{\mathrm{pp}}\left(m\right)=W*r_{\mathrm{pp}}\left(m\right)\\ G_{\mathrm{ps}}\left(m\right)=W*r_{\mathrm{ps}}\end{array}\right..$

在上述任一方案中优选的是，所述后验概率函数表示为：

$P\left(m\right|d)=\frac{P\left(d\right|m)P\left(m\right)}{P\left(d\right)}.$

本技术方案所述的地震叠前AVO反演基于精确Zoeppritz方程，相对于Zoeppritz近似方程而言，精确Zoeppritz方程克服了因入射角(或偏移距)过大而带来计算误差的问题，从而使得能够利用大偏移距信息进行储层参数预测，具有更好的适应性，满足泥砂岩薄互层结构，纵向弹性参数差异大的储层条件，如页岩气储层，有效提高叠前反演精度。本技术方案中联合PP波与PS波进行三参数反演，由于PS波地震数据包含更多的横波速度和密度信息，该方案能有效提高反演的稳定性，同时可以直接对纵波速度、横波速度、密度三个弹性参数进行反演，很好的消除了间接求解对反演精度的影响。本方案通过三变量修正柯西先验分布引入模型参数之间的相关性，在降低反演不确定性的同时，很好的保护了弱反射信息，提高了反演的精度，引入的模型参数的先验分布中包含大量横波速度和密度信息，因此在反演时不仅可以直接反演横波速和度密度，而且还具有很高的精度，即使在只有纵波地震资料的情况下也能很好的反演横波速度和密度。

附图说明

图1是按照本发明的基于精确Zoeppritz方程的多波AVO反演方法流程图。

图2是按照本发明一实施例输入的叠前AVO(Amplitude Versus Offset，振幅随偏移距的变化)PP波角道集地震数据示意图。

图3是按照本发明一实施例输入的叠前AVO(Amplitude Versus Offset，振幅随偏移距的变化)PS波角道集地震数据示意图。

图4是按照本发明一实施例输入的含噪叠前AVO(Amplitude Versus Offset，振幅随偏移距的变化)PP波角道集地震数据。

图5是按照本发明一实施例输入的含噪叠前AVO(Amplitude Versus Offset，振幅随偏移距的变化)PS波角道集地震数据。

图6是按照本发明一实施例无噪声情况下采用基于精确Zoeppritz方程的多波AVO反演方法得到的纵波速度(左)、横波速度(中)和密度(右)模型示意图。

图7是按照本发明一实施例信噪比为4:1的情况下采用基于精确Zoeppritz方程的多波AVO反演方法得到的纵波速度(左)、横波速度(中)和密度(右)模型示意图。

具体实施方式

下面参照附图结合示例性的实施例对本发明进行详细描述。

如图1所示，是按照本发明的基于精确Zoeppritz方程的多波AVO反演方法流程图，具体包括：

101、基于地震数据提取角度依赖的子波；基于测井数据和精确Zoeppritz方程正演模拟地震角道集并结合实际井旁地震数据确定振幅缩放因子：首先，基于实际的地震叠前道集和测井数据采取统计方法提取若干个依赖于入射角度的地震子波；以测井数据为输入模型利用精确Zoeppritz方程正演模拟角度域的PP波与PS波道集，与实际井旁角度域地震道集对比，计算振幅缩放因子，并应用于所提取的地震子波。

102、基于工区内所有测井数据统计模型参数的先验信息，包括纵波速度、横波速度和密度，并计算包含三参数相关性的协方差矩阵：假设模型参数服从三变量修正柯西分布，通过测井数据分析获得需要的模型参数，求取各参数的自相关系数和互相关系数，构建三参数相关的协方差矩阵，形成符合该工区的模型参数先验分布函数。

103、利用地震构造解释资料和测井数据，基于沉积模式建立时间域的初始弹性参数模型：利用地震构造解释资料，基于沉积模式建立地质模型，并将测井资料，按构造模式进行插值和外推，得到每条测线的初始弹性参数模型，包括纵波速度模型、横波速度模型和密度模型。

104、基于时间域的初始弹性参数模型和精确Zoeppritz方程正演模拟角度域的多波地震叠前道集，由正演模拟记录和实际记录直接计算PP波与PS波反演残差：以时间域的初始模型为输入，直接利用精确Zoeppritz方程计算每一个角度的PP波与PS波反射系数向量，将子波与反射系数进行褶积，获得角度域的PP波与PS波角道集，并与实际角度域地震道集作差，获得多波联合反演残差。

105、基于Bayesian原理、先验信息和精确Zoeppritz方程正演算子构建最大后验概率意义下的反演目标函数，利用高斯-牛顿法对反演目标函数进行求解，得到模型参数扰动量的求解表达式：利用以上所述先验模型对反演结果进行约束，并假设地震数据噪声服从Gauss(高斯)分布，先验信息服从修正柯西分布，则反演似然函数和先验概率分布分别满足Gauss分布和修正Cauchy(柯西)分布。根据Bayesian原理，综合反演似然函数和先验分布函数得到后验概率分布函数。根据后验概率确定反演的目标函数，由于模型参数不好直接求解，因此考虑借助于泰勒展开对目标函数进行简化，变成关于模型参数扰动量的函数，然后将目标函数对扰动量进行求导并令导数等于零，得到扰动量的迭代求解公式。其中泰勒展开需要求取精确Zoeppritz方程正演算子关于模型参数的一阶偏导数，该导数可以通过解析或数值方法求得。

106、根据模型参数扰动量的求解表达式公式和反演残差，采用迭代加权最小二乘算法计算模型参数的扰动量，重复迭代以上步骤104、105和106，并通过反演残差控制最大迭代次数，获得最优模型参数，包括纵波速度、横波速度和密度：上述目标函数是非线性方程，本发明实施例采用高斯-牛顿法迭代求解使目标函数取极小值的最优模型参数。利用以上所述加上扰动量后的参数模型为输入，重复迭代以上所述步骤104、105和106，并通过反演残差和最大迭代次数确定最优的弹性参数模型作为地震三参数叠前AVO反演最终结果，包括纵波速度、横波速度和密度模型。

本实施例具体采取以下工作步骤来实现上述技术方案：基于地震数据提取角度依赖的子波，通过测井数据正演模拟和井旁道地震数据确定振幅缩放因子→基于工区内所有测井数据统计模型参数的先验信息→利用地震构造解释资料和测井数据，基于沉积模式建立初始弹性参数模型→利用精确Zoeppritz方程正演模拟角度域的多波地震叠前道集并计算反演残差→基于广义线性反演思想将目标函数改写成关于模型参数扰动量的函数，然后对扰动量求导并令导数等于零得到扰动量的迭代求解公式→采用迭代加权最小二乘算法对模型参数扰动量进行求解→根据高斯-牛顿迭代公式更新模型参数→重复以上工作步骤直至反演残差达到要求或达到最大迭代次数→输出最终模型参数，包括纵波速度、横波速度、密度模型，用于油气藏储层预测。技术方案与工作步骤详细叙述如下：

(1)本方案假设反演之前地震子波已知，因而，需要基于实际的地震叠前道集和测井数据采取统计方法提取子波，子波在传播过程中受地层的影响会发生波形或频率变化，提取依赖于入射角度的地震子波能有效提高振幅匹配程度。

实际的地震振幅往往是相对值，采用精确Zoeppritz方程正演模拟的地震数据振幅与实际振幅存在一定数值差异。以测井数据为输入模型利用精确Zoeppritz方程正演模拟角度域的PP波与PS波道集，与实际井旁角度域地震道集对比，计算振幅缩放因子，并应用于所提取的地震子波，达到模拟记录与实际记录的振幅匹配。当地震数据信噪比较高时，为角道集的每一道使用统一的振幅缩放因子，以保证振幅随偏移距的变化关系；当信噪比低时，可分近、中、远偏移距分别计算振幅缩放因子，保证模拟记录与实际记录的最佳匹配，减少噪声对反演过程的影响。

(2)建立模型参数的先验信息。为了降低地震反演的不确定性，提高反演过程的稳定性，需要从其它途径获得地下介质地震弹性参数模型的信息作为先验信息。本发明采用三变量修正Cauchy分布函数作为先验分布函数，基于工区内所有测井数据统计得到的各模型参数，求取各参数的自相关系数和互相关系数，构建三参数的协方差矩阵，形成符合该工区的模型参数先验分布函数。在后续的反演目标函数中三变量修正柯西分布函数相应的规则化表达式为：

$R\left(m\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{m^T{\mathrm{\Φ}}^{i}m}{1+m^T{\mathrm{\Φ}}^{i}m}---\left(1\right)$

其中，m＝[vp,vs,den]^T为模型参数，上式中，Φⁱ＝(Dⁱ)^TΨ^-1Dⁱ，Ψ是3×3的协方差矩阵，它包含三个AVO参数之间统计相关性，D是3×3N的矩阵，具体形式如下所示：

$\left[D_{\mathrm{nl}}^i\right]=\left\{\begin{array}{cc}1,& \mathrm{if\; n}=1\mathrm{and\; l}=i\\ 1,& \mathrm{if\; n}=2\mathrm{and\; l}=i+N\\ 1,& \mathrm{if\; n}=3\mathrm{and\; l}=i+2N\\ 0,& \mathrm{other}\end{array}\right.---\left(2\right)$

(3)利用地震构造解释资料，基于沉积模式建立地质模型，并将测井资料，按构造模式进行插值和外推，得到每条测线的初始弹性参数模型，包括纵波速度模型、横波速度模型和密度模型。

建立弹性参数模型主要利用三维空间插值方法，其的技术流程是首先利用散点插值的方法对各个层位的数据进行插值，完成地质层位建模，然后根据地质层位进行弹性参数横向插值，即将测井信息进行横向插值，计算得到地下每个点上的弹性参数值，完成初始弹性参数建模的任务。

(4)基于时间域的初始弹性参数模型和精确Zoeppritz方程正演模拟角度域的多波地震叠前道集，具体包括：以时间域的初始模型为输入，直接利用精确Zoeppritz方程计算每一个角度的PP波与PS波反射系数向量，将子波与反射系数进行褶积，获得角度域的PP波与PS波角道集。将模拟的角道集数据并与实际角度域地震道集作差(如图2-5所示，是本实施例输入的实际叠前AVO角道集地震数据，图2为不含噪的PP波，图3为不含噪的PS波，图4为含噪的PP波，图5为含噪的PS波，图中纵轴表示时间，单位为秒，横轴表示入射角度)，获得多波联合反演残差。基于精确Zoeppritz方程的正演过程可表示如下：

$\left\{\begin{array}{c}{G}_{\mathrm{pp}}\left(m\right)=W*r_{\mathrm{pp}}\left(m\right)\\ G_{\mathrm{ps}}\left(m\right)=W*r_{\mathrm{ps}}\end{array}\right.---\left(3\right)$

(5)基于Bayesian原理、先验信息和精确Zoeppritz方程正演算子构建最大后验概率意义下的反演目标函数，利用高斯-牛顿法对反演目标函数进行求解，得到模型参数扰动量的求解表达式：利用以上所述先验模型对反演结果进行约束，并假设地震数据噪声服从Gauss(高斯)分布，而先验信息服从修正Cauchy(柯西)分布，则反演似然函数和先验概率分布分别满足Gauss分布和修正柯西分布。根据Bayesian原理，综合反演似然函数和先验分布函数得到后验概率分布函数。根据后验概率确定反演的目标函数，由于模型参数不好直接求解，因此考虑借助于泰勒展开对目标函数进行简化，变成关于模型参数扰动量的函数，然后将目标函数对扰动量进行求导并令导数等于零，得到扰动量的迭代求解公式。其中泰勒展开需要求取精确Zoeppritz方程正演算子关于模型参数的一阶导数，该导数可以通过解析或数值方法求得。设模型参数为m^T＝(m₁,m₂,…,m_n)^T，观测地震数据为d^T＝(d₁,d₂,…,d_n)^T，由贝叶斯理论可以得到，在已知叠前地震数据的情况下，反演地下介质弹性参数的问题可以归结为求解一个后验概率函数，

$P\left(m\right|d)=\frac{P\left(d\right|m)P\left(m\right)}{P\left(d\right)}---\left(4\right)$

其中P(d)＝∫p(d|m)p(m)dm为归一化因子，可以看作是常数，P(d|m)为似然函数，P(m)为先验概率分布。假设似然函数和先验分布分别满足高斯分布和修正柯西分布，

$\begin{array}{c}P\left(m\right)==\frac{1}{{\mathrm{\π}}^{\left(2N\right)}{\left|\mathrm{\ψ}\right|}^{N/2}}\exp (-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{m^T{\mathrm{\Φ}}^{i}m}{1+m^T{\mathrm{\Φ}}^{i}m})\\ P\left(d\right|m)={\left({\left(2\mathrm{\π}\right)}^n\right|C_D\left|\right)}^{1/2}\exp (-\frac{1}{2}{(d-G\left(m\right))}^T{C}_D^{-1}(d-G\left(m\right)))\end{array}---\left(5\right)$

则多波地震数据后验概率分布可以表示为，

$\begin{array}{c}P\left(m\right|d)\∝P\left(d\right|m)P\left(m\right)\\ \∝{\left({\left(2\mathrm{\π}\right)}^{\mathrm{Npp}}\right|C_{\mathrm{Dpp}}\left|\right)}^{1/2}\exp (-\frac{1}{2}{(d_{\mathrm{pp}}-G_{\mathrm{pp}}\left(m\right))}^T{C}_{\mathrm{Dpp}}^{-1}(d_{\mathrm{pp}}-G_{\mathrm{pp}}\left(m\right)))\\ *{\left({\left(2\mathrm{\π}\right)}^{\mathrm{Nps}}\right|C_{\mathrm{Dps}}\left|\right)}^{1/2}\exp (-\frac{1}{2}{(d_{\mathrm{ps}}-G_{\mathrm{ps}}\left(m\right))}^T)C_{\mathrm{Dps}}^{-1}(d_{\mathrm{ps}}-G_{\mathrm{ps}}\left(m\right))\\ *\exp (-R\left(m\right))\end{array}---\left(6\right)$

一般可通过求解上式的最大值作为最优解，等价于求解下式目标函数的最小值所对应的解，

$\begin{array}{c}S\left(m\right)=\frac{1}{2}{(d_{\mathrm{PP}}-G_{\mathrm{PP}}\left(m\right))}^T{C}_{\mathrm{Dpp}}^{-1}(D_{\mathrm{PP}}-G_{\mathrm{PP}}\left(m\right))\\ +\frac{1}{2}{(d_{\mathrm{PS}}-G_{\mathrm{PS}}\left(m\right))}^T{C}_{\mathrm{Dps}}^{-1}(d_{\mathrm{PS}}-G_{\mathrm{PS}}\left(m\right))+R\left(m\right)\end{array}---\left(7\right)$

可通过高斯-牛顿方法迭代求解该目标函数，为了简化求解过程，减少计算量，考虑将目标函数进行简化改写。将上式中的精确Zoeppritz方程正演算子在初始模型参数m⁰(m¹＝m⁰+Δm)处进行泰勒展开，如下：

$G\left(m^1\right)=G_{\mathrm{pp}}\left(m^0\right)+\frac{\∂G_{\mathrm{pp}}\left(m^0\right)}{\∂m}\×\mathrm{\Δm}+G_{\mathrm{ps}}\left(m_0\right)+\frac{\∂G_{\mathrm{ps}}\left(m^0\right)}{\∂m}\×\mathrm{\Δm}+\·\·\·\left(8\right)$

对上式取一阶近似，即利用线性关系近似表达非线性关系。将上式带入式(7)得：

$\begin{array}{c}S\left(m\right)=\frac{1}{2}{(d_{\mathrm{pp}}-G_{\mathrm{pp}}\left(m^0\right)-\frac{\∂G_{\mathrm{pp}}\left(m^0\right)}{\∂m}\×\mathrm{\Δm})}^T{C}_{\mathrm{Dpp}}^{-1}(d_{\mathrm{pp}}-G_{\mathrm{pp}}\left(m^0\right)-\frac{\∂G_{\mathrm{pp}}\left(m^0\right)}{\∂m}\×\mathrm{\Δm})\\ +\frac{1}{2}{(d_{\mathrm{ps}}-g_{\mathrm{ps}}\left(m^0\right)-\frac{{\∂G}_{\mathrm{ps}}\left(m^0\right)}{\∂m}\×\mathrm{\Δm})}^T{C}_{\mathrm{Dps}}^{-1}(d_{\mathrm{ps}}-G_{\mathrm{ps}}\left(m^0\right)-\frac{\∂G_{\mathrm{ps}}\left(m^0\right)}{\∂}\×\mathrm{\Δm})+R\left(m\right)\end{array}---\left(9\right)$

对于上式中的R(m)进行如下改写：

$R\left(m\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{m^T{\mathrm{\Φ}}^{i}m}{1+m^T{\mathrm{\Φ}}^{i}m}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{{(m^0+\mathrm{\Δm})}^T{\mathrm{\Φ}}^i(m^0+\mathrm{\Δm})}{1+{(m^0+\mathrm{\Δm})}^T{\mathrm{\Φ}}^i(m^0+\mathrm{\Δm})}---\left(10\right)$

目标函数变为如下形式：

$\begin{array}{c}S\left(m\right)=\frac{1}{2}{(d_{\mathrm{pp}}-G_{\mathrm{pp}}\left(m^0\right)-\frac{{\∂G}_{\mathrm{pp}}\left(m^0\right)}{\∂m}\×\mathrm{\Δm})}^T{C}_{\mathrm{Dpp}}^{-1}(d_{\mathrm{pp}}-G_{\mathrm{pp}}\left(m^0\right)-\frac{\∂G_{\mathrm{pp}}\left(m^0\right)}{\∂m}\×\mathrm{\Δm})\\ +\frac{1}{2}{(d_{\mathrm{ps}}-G_{\mathrm{ps}}\left(m^0\right)-\frac{{\∂G}_{\mathrm{ps}}\left(m^0\right)}{\∂m}\×\mathrm{\Δm})}^T{C}_{\mathrm{Dps}}^{-1}(d_{\mathrm{ps}}-G_{\mathrm{ps}}\left(m^0\right)-\frac{{\∂G}_{\mathrm{ps}}\left(m^0\right)}{\∂m}\×\mathrm{\Δm})\\ +\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{{(m^0+\mathrm{\Δm})}^T{\mathrm{\Φ}}^i(m^0+\mathrm{\Δm})}{1+{(m^0+\mathrm{\Δm})}^T{\mathrm{\Φ}}^i(m^0+\mathrm{\Δm})}\end{array}---\left(11\right)$

若假设地震数据中随机噪声之间相互独立，并服从高斯分布，PP波地震数据的噪音方差为σ_pp，PS波地震数据的噪音方差为σ_ps，则目标函数中的协方差矩阵和可以简化对角线为常数的对角阵。上式可简化为：

$\begin{array}{c}S\left(m\right)=\frac{1}{2}{(d_{\mathrm{pp}}-G_{\mathrm{pp}}\left(m^0\right)-\frac{{\∂G}_{\mathrm{pp}}\left(m^0\right)}{\∂m}\×\mathrm{\Δm})}^T(d_{\mathrm{pp}}-G_{\mathrm{pp}}\left(m^0\right)-\frac{\∂G_{\mathrm{pp}}\left(m^0\right)}{\∂m}\×\mathrm{\Δm})\\ +\frac{\mathrm{\α}}{2}{(d_{\mathrm{ps}}-G_{\mathrm{ps}}\left(m^0\right)-\frac{{\∂G}_{\mathrm{ps}}\left(m^0\right)}{\∂m}\×\mathrm{\Δm})}^T(d_{\mathrm{ps}}-G_{\mathrm{ps}}\left(m^0\right)-\frac{{\∂G}_{\mathrm{ps}}\left(m^0\right)}{\∂m}\×\mathrm{\Δm})\\ +\mathrm{\β}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{{(m^0+\mathrm{\Δm})}^T{\mathrm{\Φ}}^i(m^0+\mathrm{\Δm})}{1+{(m^0+\mathrm{\Δm})}^T{\mathrm{\Φ}}^i(m^0+\mathrm{\Δm})}\end{array}---\left(12\right)$

其中α＝σ_pp/σ_ps控制纵横波数据的使用比重，β＝σ_pp控制先验信息的权重。将上面得到的目标函数对弹性参数扰动量Δm进行求导得：

$\begin{array}{c}\frac{\∂S\left(m\right)}{\∂\mathrm{\Δm}}={\left(\frac{{\∂G}_{\mathrm{pp}}}{\∂m}\right)}^T\left(\frac{{\∂G}_{\mathrm{pp}}}{\∂m}\right)\×\mathrm{\Δm}-{\left(\frac{{\∂G}_{\mathrm{pp}}}{\∂m}\right)}^T(d_{\mathrm{pp}}-G_{\mathrm{pp}}\left(m^0\right))\\ +\mathrm{\α}{\left(\frac{{\∂G}_{\mathrm{ps}}}{\∂m}\right)}^T\left(\frac{{\∂G}_{\mathrm{ps}}}{\∂m}\right)\×\mathrm{\Δm}-\mathrm{\α}{\left(\frac{{\∂G}_{\mathrm{ps}}}{\∂m}\right)}^T(d_{\mathrm{ps}}-G_{\mathrm{ps}}\left(m^0\right))\\ +\mathrm{\β}(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{2{\mathrm{\Φ}}^i{m}^0}{1+{(m^0+\mathrm{\Δm})}^T{\mathrm{\Φ}}^i{(m^0+\mathrm{\Δm})}^T}+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{2{\mathrm{\Φ}}^i\mathrm{\Δm}}{1+{(m^0+\mathrm{\Δm})}^T{\mathrm{\Φ}}^i(m^0+\mathrm{\Δm})})\end{array}---\left(13\right)$

令等于零，可以模型参数扰动量的表达式：

$\mathrm{\Δm}=H_1^{-1}{\mathrm{\γ}}_1---\left(14\right)$

其中， $H_1={\left(\frac{{\∂G}_{\mathrm{pp}}\left(m^0\right)}{\∂m}\right)}^T\left(\frac{{\∂G}_{\mathrm{pp}}\left(m^0\right)}{\∂m}\right)+\mathrm{\α}{\left(\frac{{\∂G}_{\mathrm{pp}}\left(m^0\right)}{\∂m}\right)}^T\left(\frac{{\∂G}_{\mathrm{ps}}\left(m^0\right)}{\∂m}\right)+\mathrm{\β}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{2{\mathrm{\Φ}}^i}{1+{(m^0+\mathrm{\Δm})}^T{\mathrm{\Φ}}^i(m^0+\mathrm{\Δm})},$ ${\mathrm{\γ}}_1={\left(\frac{{\∂G}_{\mathrm{pp}}\left(m^0\right)}{\∂m}\right)}^T(d_{\mathrm{pp}}-G_{\mathrm{pp}}\left(m^0\right))+\mathrm{\α}{\left(\frac{{\∂G}_{\mathrm{pp}}\left(m^0\right)}{\∂m}\right)}^T(d_{\mathrm{ps}}-G_{\mathrm{ps}}\left(m^0\right))-\mathrm{\β}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{2{\mathrm{\Φ}}^i{m}^0}{1+{(m^0+\mathrm{\Δm})}^T{\mathrm{\Φ}}^i(m^0+\mathrm{\Δm})}.$ 由于等式两边都含有Δm，因此采用迭代加权最小平方算法进行求解，初始扰动量Δm⁰取为零。由于该方法收敛速度快，并且该方法属于多层迭代，为了减少计算量，一般在此计算扰动量时只迭代一次。则扰动量的求解表达式最终可以简化为：

Δm_k＝H^-1γ(15)

其中， $H={\left(\frac{{\∂G}_{\mathrm{pp}}\left(m^k\right)}{\∂m}\right)}^T\left(\frac{{\∂G}_{\mathrm{pp}}\left(m^k\right)}{\∂m}\right)+\mathrm{\α}{\left(\frac{{\∂G}_{\mathrm{ps}}\left(m^k\right)}{\∂m}\right)}^T\left(\frac{{\∂G}_{\mathrm{ps}}\left(m^k\right)}{\∂m}\right)+\mathrm{\β}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{2{\mathrm{\Φ}}^i}{1+{\left(m^k\right)}^T{\mathrm{\Φ}}^i\left(m^k\right)},$ $\mathrm{\γ}={\left(\frac{{\∂G}_{\mathrm{pp}}\left(m^k\right)}{\∂m}\right)}^T(d_{\mathrm{pp}}-G_{\mathrm{pp}}\left(m^k\right))+\mathrm{\α}{\left(\frac{{\∂G}_{\mathrm{ps}}\left(m^k\right)}{\∂m}\right)}^T(d_{\mathrm{ps}}-G_{\mathrm{ps}}\left(m^k\right))-\mathrm{\β}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{2{\mathrm{\Φ}}^i{m}^k}{1+{\left(m^k\right)}^T{\mathrm{\Φ}}^i\left(m^k\right)}.$ 其中正演算子关于模型参数的一阶导数与可通过解析的方法求得，

$\begin{array}{c}\frac{{\∂G}_{\mathrm{pp}}}{\∂m}=W*\frac{{\mathrm{\δr}}_{\mathrm{pp}}}{\mathrm{\δm}}\\ \frac{{\∂G}_{\mathrm{ps}}}{\∂m}=W*\frac{{\mathrm{\δr}}_{\mathrm{ps}}}{\mathrm{\δm}}\end{array}---\left(16\right)$

对单一界面，纵、横波反射系数r_pp和r_ps可以由下列方程计算得到：

MS＝N   (17)

其中，S为反射、透射系数矩阵，形式如下：

$S=\left[\begin{array}{cccc}{r}_{\mathrm{pp}}^D& r_{\mathrm{sp}}^D& r_{\mathrm{pp}}^U& r_{\mathrm{sp}}^U\\ r_{\mathrm{ps}}^D& r_{\mathrm{ss}}^D& r_{\mathrm{ps}}^U& r_{\mathrm{ss}}^U\\ t_{\mathrm{pp}}^D& t_{\mathrm{sp}}^D& t_{\mathrm{pp}}^U& t_{\mathrm{sp}}^U\\ t_{\mathrm{ps}}^D& t_{\mathrm{ss}}^D& t_{\mathrm{ps}}^U& t_{\mathrm{ss}}^U\end{array}\right]---\left(18\right)$

上式中，r和t分别表示弹性波入射到分界面时的反射、透射系数，D和U表示入射波为下行波、上行波。下标的第一个字母代表入射波的波型，第二个字母代表反射波或透射波的波型。本文是基于精确Zoeppritz方程的贝叶斯非线性反演，假设叠前地震数据已经过透射损失校正、几何扩散效应补偿以及衰减补偿。对于PP波和PS波联合AVO非线性反演，只需计算PP波反射系数和PS转换波反射系数即可。

M表示为：

$M=\left[\begin{array}{cccc}-{\mathrm{\α}}_{1}p& -\cos j_1& {\mathrm{\α}}_{2}p& \cos j_2\\ \cos i_1& -{\mathrm{\β}}_{1}p& \cos i_2& {-\mathrm{\β}}_{2}p\\ 2{\mathrm{\ρ}}_1{\mathrm{\β}}_1^{2}p\cos i_1& {\mathrm{\ρ}}_1{\mathrm{\β}}_1(1-2{\mathrm{\β}}_1^2{p}^2)& 2{\mathrm{\ρ}}_2{\mathrm{\β}}_2^{2}p\cos i_2& {\mathrm{\ρ}}_2{\mathrm{\β}}_2(1-2{\mathrm{\β}}_2^2{p}^2)\\ -{\mathrm{\ρ}}_1{\mathrm{\α}}_1(1-2{\mathrm{\β}}_1^2{p}^2)& 2{\mathrm{\ρ}}_1{\mathrm{\β}}_1^{2}p\cos j_1& {\mathrm{\ρ}}_2{\mathrm{\α}}_2(1-2{\mathrm{\β}}_2^2{p}^2)& -2{\mathrm{\ρ}}_2{\mathrm{\β}}_2^{2}p\cos j_2\end{array}\right]---\left(19\right)$

N表示为：

$N=\left[\begin{array}{cccc}-{\mathrm{\α}}_{1}p& \cos j_1& -{\mathrm{\α}}_{2}p& -\cos j_2\\ \cos i_1& -{\mathrm{\β}}_{1}p& \cos i_2& {-\mathrm{\β}}_{2}p\\ 2{\mathrm{\ρ}}_1{\mathrm{\β}}_1^{2}p\cos i_1& {\mathrm{\ρ}}_1{\mathrm{\β}}_1(1-2{\mathrm{\β}}_1^2{p}^2)& 2{\mathrm{\ρ}}_2{\mathrm{\β}}_2^{2}p\cos i_2& {\mathrm{\ρ}}_2{\mathrm{\β}}_2(1-2{\mathrm{\β}}_2^2{p}^2)\\ {\mathrm{\ρ}}_1{\mathrm{\α}}_1(1-2{\mathrm{\β}}_1^2{p}^2)& -2{\mathrm{\ρ}}_1{\mathrm{\β}}_1^{2}p\cos j_1& -{\mathrm{\ρ}}_2{\mathrm{\α}}_2(1-2{\mathrm{\β}}_2^2{p}^2)& 2{\mathrm{\ρ}}_2{\mathrm{\β}}_2^{2}p\cos j_2\end{array}\right]---\left(20\right)$

假设c^j＝[c₁ c₂ c₃ c₄ c₅ c₆]＝[α₁ α₂ β₁ β₂ ρ₁ ρ₂]为第j分个界面两侧地层的弹性参数，c₁-c₆为分界面上侧地层的纵波速度、横波速度、密度以及分界面下侧地层的纵波速度、横波速度、密度。公式(16)两边对求导得：

$\frac{\∂M}{\∂c_i^j}S+M\frac{\∂S}{\∂c_i^j}=\frac{\∂N}{\∂c_i^j}---\left(21\right)$

整理，写成散射系数关于模型参数的一阶偏导数公式：

$\frac{\∂S}{\∂c_i^j}=M^{-1}(\frac{\∂N}{\∂c_i^j}=\frac{\∂M}{\∂c_i^j}{M}^{-1}N)---\left(22\right)$

根据式(20)，我们只需要让等式右边对应行和列相乘就能得到第j个界面的纵波反射系数和转换波反射系数关于界面两侧的纵波速度、横波速度、密度的一阶偏导数。最终，我们可以得到模型参数的更新迭代公式：

m^k+1＝m^k+Δm_k   (23)

(6)根据模型参数扰动量的求解表达式公式和反演残差，采用迭代加权最小二乘算法计算模型参数的扰动量，重复迭代以上步骤4、5和6，并通过反演残差控制最大迭代次数，获得最优模型参数，包括纵波速度、横波速度和密度(如图6所示，是本实施例基于精确Zoeppritz方程的多波AVO反演得到的纵波速度(左)、横波速度(中)和密度(右)模型，图中纵轴表示时间，单位秒，横轴自左至右表示纵波速度(单位：千米/秒)、横波速度(单位：千米/秒)和密度(单位：克/立方厘米))。基于精确Zoeppritz方程的的多波AVO反演方法不仅对纵波速度、横波速度有精确的预测结果，由于引入了包含密度信息的先验分布，并采用PS横波进行联合反演，其对密度模型也预测准确。如图7所示为加入信噪比为4时的随机噪声的反演结果，三变量修正柯西先验约束不仅对保持反演过程稳定起到了关键作用，而且为反演带来了稀疏解并同时很好的保护了弱反射信息。

本发明由于采取以上技术方案，其具有以下优点：1、本技术方案所述的地震叠前AVO反演基于精确Zoeppritz方程，相对于常规AVO反演而言，很好的克服了大偏移距带来的计算误差，从而使得我们能够同时利用多种信息进行储层参数预测；2、本技术方案是基于精确Zoeppritz方程的AVO反演方法，具有更好的适应性，满足泥砂岩薄互层结构，纵向弹性参数差异大的储层条件，如页岩气储层，有效提高叠前反演精度；3、本技术方案联合PP波与PS波进行三参数反演，由于PS波地震数据包含更多的横波速度和密度信息，该方案能有效提高反演的稳定性；4、本技术方案可以直接对纵波速度、横波速度、密度三个弹性参数进行反演，很好的消除了间接求解对反演精度的影响；5、本技术方案通过三变量修正柯西先验分布引入模型参数之间的相关性，在降低反演不确定性的同时，很好的保护了弱反射信息，提高了反演的精度。6、本技术方案引入的模型参数的先验分布中包含大量横波速度和密度信息，因此在反演时不仅可以直接反演横波速和度密度，而且还具有很高的精度，即使在只有纵波地震资料的情况下也能很好的反演横波速度和密度。

为了更好地理解本发明，以上结合具体实施例对本发明作了详细说明。但是，显然可对本发明进行不同的变型和改型而不超出权利要求限定的本发明更宽的精神和范围。因此，以上实施例具有示例性而没有限制的含义。

Claims (10)

1.基于精确Zoeppritz方程的多波AVO储层弹性参数反演方法，包括以下步骤：

步骤a. 确定振幅缩放因子；

步骤b. 统计模型参数的先验信息，并计算包含三参数相关性的协方差矩阵；

步骤c. 建立时间域的初始弹性参数模型；

步骤d. 基于所述弹性参数模型和精确Zoeppritz方程正演模拟角度域的多波地震叠前道集，由正演模拟记录和实际记录计算PP波与PS波反演残差；

步骤e. 构建反演目标函数，对所述反演目标函数求解，得到模型参数扰动量的求解表达式；

步骤f. 根据所述反演残差和所述扰动量求解表达式计算模型参数的扰动量；

步骤g. 重复迭代步骤d、e、f获得最优模型参数。

2.根据权利要求1所述的基于精确Zoeppritz方程的多波AVO储层弹性参数反演方法，其特征在于，所述步骤a确定振幅缩放因子包括：

基于地震数据提取角度依赖的子波；

基于测井数据和Zoeppritz正演模拟地震角道集并结合实际井旁地震数据确定振幅缩放因子。

3.根据权利要求1所述的基于精确Zoeppritz方程的多波AVO储层弹性参数反演方法，其特征在于，所述步骤b包括：

采用三变量修正Cauchy分布函数作为先验分布函数；

基于工区内所有测井数据统计得到各模型参数；

求取所述各模型参数的自相关系数和互相关系数，构建三参数的协方差矩阵，形成符合该工区的模型参数先验分布函数。

4.根据权利要求1所述的基于精确Zoeppritz方程的多波AVO储层弹性参数反演方法，其特征在于，所述先验信息包括：纵波速度、横波速度和密度。

5.根据权利要求1所述的基于精确Zoeppritz方程的多波AVO储层弹性参数反演方法，其特征在于，所述步骤c建立初始弹性参数模型包括：

利用地震构造解释资料，基于沉积模式建立地质模型；

将测井资料按构造模式进行插值和外推，得到每条测线的初始弹性参数模型。

6.根据权利要求1所述的基于精确Zoeppritz方程的多波AVO储层弹性参数反演方法，其特征在于，所述弹性参数模型包括：纵波速度模型、横波速度模型和密度模型。

7.根据权利要求1所述的基于精确Zoeppritz方程的多波AVO储层弹性参数反演方法，其特征在于，所述步骤d包括：

以时间域的初始模型作为输入，利用精确Zoeppritz方程计算每个角度的PP波与PS波的反射系数向量；

将子波与反射系数进行褶积，获得角度域的PP波与PS波角道集；

将模拟的角道集数据与实际角度域地震道集作差，获得多波联合反演残差。

8.根据权利要求1所述的基于精确Zoeppritz方程的多波AVO储层弹性参数反演方法，其特征在于，所述步骤e包括：

根据Bayesian原理，综合反演似然函数和先验分布函数得到后验概率分布函数；

根据后验概率确定反演目标函数；

借助于泰勒展开对目标函数进行简化，得到关于模型参数的扰动量函数；

将目标函数对扰动量求导并令导数等于零，得到扰动量的迭代求解公式。

9.根据权利要求1所述的基于精确Zoeppritz方程的多波AVO储层弹性参数反演方法，其特征在于，所述对反演目标函数的求解方法包括：高斯-牛顿法。

10.根据权利要求1所述的基于精确Zoeppritz方程的多波AVO储层弹性参数反演方法，其特征在于，所述步骤f中计算模型参数的扰动量的方法包括：迭代加权最小二乘法。