CN104007467A - 一种基于混合范数正则化的叠前三参数反演实现的储层与流体预测方法 - Google Patents

Abstract

本发明涉及一种基于混合范数正则化的叠前三参数反演实现的储层与流体预测方法，基于贝叶斯参数估计理论和AVO基本原理，综合稀疏理论、剩余残差分布、地质意义等因素建立了混合罚函数方程，对叠前三参数反演过程进行正则化约束。模型与实际资料应用证实：混合范数正则化可进增加结果稳定性与地质意义，特别是使得相同偏移距和信噪比条件下的密度反演结果更加准确。本发明优势在于可直接应用于实际资料，获得合理的反演结果，提高了储层与流体预测精度。

Description

一种基于混合范数正则化的叠前三参数反演实现的储层与流体预测方法

技术领域

本发明属于石油勘探领域，涉及一种基于混合范数正则化的叠前三参数反演实现的储层与流体预测方法，适应于实际不适定AVO（Amplitude versusoffset）问题的弹性参数反演方法。

背景技术

叠前反演可以同时获得纵波、横波、纵横波速度比、泊松比和拉梅常数等多种弹性参数，从而从弹性域对储层与流体进行更为精确的综合预测（Zhang etal.,2013）。实际连续多界面的纵波AVO反射系数方程近似表示为（Buland和More，2003）：

$R_{\mathrm{pp}}(t,\mathrm{\θ})\≈\frac{{\sec }^2\mathrm{\θ}}{2}\frac{\∂}{\∂t}\ln I_p\left(t\right)-4\frac{v_s^2}{v_p^2}{\sin }^2\mathrm{\θ}\frac{\∂}{\∂t}\ln I_s\left(t\right)+\frac{1}{2}(4\frac{v_s^2}{v_p^2}{\sin }^2\mathrm{\θ}-{\tan }^2\mathrm{\θ})\frac{\∂}{\∂t}\ln \mathrm{\ρ}\left(t\right)---\left(1\right)$

其中，θ=(θ₁+θ₂)/2，θ₂=arcsin(sinθ₁/v_p1 ^*v_p2)，θ₁和θ₂分别为纵波入射与透射角度，I_p与I_s分别为上下层介质纵波与横波阻抗，v_p1，v_p2，v_p；v_s1，v_s2，v_s；ρ₁，ρ₂，ρ为别为上下层纵波、横波速度与密度及其均值。其加入子波矩阵以考虑地震资料的带限特征及调谐效应，综合噪音因素后，可得到如下K个界面N个入射角度的单道反演公式的离散形式为：

$d(t,\mathrm{\θ})=\frac{1}{2}A\left(\mathrm{\θ}\right)W_{\mathrm{\θ}}\mathrm{gD}{L}_P+\frac{1}{2}B\left(\mathrm{\θ}\right)W_{\mathrm{\θ}}g{\mathrm{DL}}_S+C\left(\mathrm{\θ}\right)W_{\mathrm{\θ}}g{\mathrm{DL}}_D+n$

即， $d_{\mathrm{KN}*1}=G_{\mathrm{KN}*3K}{m}_{3K*1}+n---\left(2\right)$

其中，d为KN行实际地震观测数据；m为3K行待求反演的参数（纵波阻抗，横波阻抗和密度）；G为KN行3K列的正演算子，由子波矩阵W，一阶差分矩阵D和系数矩阵组成；n为KN行地震数据包含的噪音；L_p=Ln(I_p)，L_s=Ln(I_s)，L_D=Ln(I_D)，

$d=\left[\begin{array}{c}{d}_{11}\\ M\\ d_{1N}\\ M\\ d_{K1}\\ M\\ d_N\end{array}\right],G=\left[\begin{array}{ccc}{\mathrm{WDA}}_1& {\mathrm{WDB}}_1& {\mathrm{WDC}}_1\\ M& M& M\\ W{\mathrm{DA}}_N& {\mathrm{WDB}}_N& {\mathrm{WDC}}_N\end{array}\right],m=\left[\begin{array}{c}\ln I_{p1}\\ M\\ \ln I_{\mathrm{pK}}\\ \ln I_{s1}\\ M\\ \ln I_{\mathrm{sK}}\\ \ln {\mathrm{\ρ}}_1\\ M\\ \ln {\mathrm{\ρ}}_K\end{array}\right]$

适定反问题的实质是要求其正演算子的逆存在且连续（王彦飞，2007），而实际情况下多数反问题都不能满足此条件，具有明显的不适定性，即反问题的解不连续依赖于原始数据（肖庭延等，2003）。从AVO的基本理论可以看出：矩问题（2）本身是对连续正问题（1）在一定精度下的离散表达，而平面波近似理论下对Zoeppritz方程的近似的正演方程本身并非是地震波的精确表达，加上实际情况下数据有效带宽的固有特征和各类噪音的影响使得地球物理AVO反演问题严重病态，存在严重的不适定性，其解存在非唯一和不稳定性（陈建江等，2007；杨培杰，2008）。

对于不适定反演问题的求解，通常都是求取最小平方误差意义下的解，即求m^*使得min||Gm^*-d||，但是求得的结果非唯一并且不稳定。为此必须在求解目标函数的过程中加入先验信息来约束反演过程，以获得唯一、稳定的解，这种方法称为正则化方法。正则化方法主要用于解决第一类不适定性积分方程，最早由苏联院士Tikhonov A N在1963年提出，将不适定问题的求解转换为对变分问题min||Gm^*-d||+α·Ω[m^*]的求解，即

(G^TG+αM)m=G^Td   （3）

其中α为正则化参数，Ω[]为罚函数。M为一相应于Ω[]的（半）正定可微因子。对罚函数的类型和正则化参数的不同选择就发展了不同的正则化方法，比如针对超定问题的最小二乘解，针对欠定问题的最小长度解，针对实际离散欠定问题考虑“平缓度”或“粗糙度”的加权最小长度解，以及针对带噪问题的阻尼最小二乘和加权阻尼最小二乘解等。近年来广泛地应用于地球物理领域的贝叶斯参数估计理论能够有效地利用各方面资料作为先验信息对反演过程进行约束,提高反演结果精度，增加实际地质意义，其本质也是正则化方法的一种。据其构建反演的方程为：

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

其中，P(d|m)为观测数据分布的条件概率，P(m)为待求参数的先验分布函数，P(m|d)为待求参数的后验概率分布函数，如果仅关心P(m|d)的形状，则对于给定数据d，分母P(d)可以被忽略，即待求参数最大后验估计等价于观测数据的最大似然估计与待求参数的先验分布的乘积。典型的基于Cauchy分布假设的贝叶斯框架下的AVO反演正则化方程为：

$(G^{T}G+2\frac{{\mathrm{\σ}}_n^2}{{\mathrm{\σ}}_m^2}Q)m=G^{T}d---\left(5\right)$

其中，μ=2σ_n ²/σ_m ²为Cauchy分布因子，Q为对角矩阵，即 $Q=\mathrm{diag}\left(\frac{1}{{(1+m_i^2/{\mathrm{\σ}}_m^2)}^2}\right).$

先验分布P(m)及正则化约束条件的引入可以提高反演的稳定性，对先验信息不同的假设条件会影响最大后验解的求取结果。对于实际AVO反演问题，Buland和More(2003)提出贝叶斯框架下的叠前反演正则化方法，通过对模型参数的先验分布进行正态分布假设，并借此对反演过程引入罚函数约束，以减小解的多解性。在此基础上许多学者进行了不同的研究，如修改反演的目标参数以提高稳定性（陈建江等，2007）；考虑待求参数之间相关性，将模型参数由单变量概率分布改良为多变量联合分布（Alemie，2010）；假设噪音为非高斯分布以及对比模型在不同先验分布的假设下对结果的影响（Theune，2010）；基于Cauchy分布消除高斯分布平滑作用（Theune，2010）的自适应边缘保持叠前反演方法（张赛民，2011；王康宁等，2013）；以及获取基于贝叶斯框架下不同正演参数的反演方法（宗兆云等，2011）等。

虽然通过先验分布P(m)的引入可以提高反演的稳定性，但通常情况下对待求参数纵波阻抗、横波阻抗和密度的扰动变化量的特征呈正态分布的假设忽略了三者之间的关系。而AVO三参数反演问题本身受高度病态性的影响，其反演结果（尤其是密度项）往往极不稳定。前人曾尝试对三参数方程进行转换，去除方程中的密度项以稳定反演，但不可避免地要要权衡提高反问题稳定性和损失信息程度之间的矛盾。此外，Theune（2010）在对比不同先验分布的假设下反演结果时，发现Cauchy约束反演时有时会出现人为造成的横轴假象。实际上，现有实际AVO反演问题的正则化约束条件普遍只关注初始模型、与地震数据的吻合程度等，对实际弹性参数的地质意义考虑不足，因而最终反演的结果常受控于初始不适定问题的固有限制，表现为反演的纵波、横波与密度等弹性参数质量不一，相关性不足，大大制约了叠前反演在储层与流体预测中的应用。

发明内容

本发明提出一种基于混合范数正则化的叠前三参数反演实现的储层与流体预测方法，该方法能够用于实际不适定AVO反演问题的弹性参数求解。通过本发明提供的方法，能够在实际情况下获得更加准确的弹性参数，特别是使得相同偏移距和信噪比条件下的密度反演结果更加准确，提高储层与流体预测精度。

本发明实现上述目的的具体实施方案如下：一种提高储层与流体预测精度方法，包括以下的步骤：

步骤1：对地震资料进行保幅处理，抽取高分辨率、高信噪比，高保真的共反射点（CRP，Common reflection point）道集；

步骤2：基于高分辨率、高信噪比，高保真的CRP道集获得分角度叠加数据体，结合测井数据估算各数据体地震子波；

步骤3：通过测井数据与地质层位建立纵波、横波与密度的初始模型，统计纵波与横波速度、纵波速度与密度之间的关系；

步骤4：基于估算的子波、建立的初始模型与弹性参数关系，采用本发明的混合范数正则化方法，建立叠前AVO反演的混合范数正则化反演方程，

(G'^TG'+μQ'+βI)m'=βm₀'+G'^Td

其中，μ为Cauchy分布约束因子，β为初始模型约束因子，Ⅰ为单位矩阵，Q为对角矩阵：即 $Q^{\′}=\mathrm{diag}\left(\frac{1}{{(1+m_i^{\′2}/{\mathrm{\σ}}_{m^{\′}}^2)}^2}\right),$ $G^{\′}=\left[\begin{array}{ccc}\mathrm{WD}{A_1}^{\′}& \frac{1}{2}\mathrm{WD}{B}_1& \mathrm{WD}{C}_1\\ M& M& M\\ \mathrm{WD}{A_N}^{\′}& \frac{1}{2}\mathrm{WD}{B}_N& \mathrm{WD}{C}_N\end{array}\right];$ $m^{\′}=\left[\begin{array}{c}\ln I_{p1}\\ M\\ \ln I_{\mathrm{pK}}\\ \mathrm{\Δ}\ln I_{s1}\\ M\\ \mathrm{\Δ}\ln I_{\mathrm{sK}}\\ \mathrm{\Δ}\ln {\mathrm{\ρ}}_1\\ M\\ \mathrm{\Δ}\ln {\mathrm{\ρ}}_K\end{array}\right];{m_0}^{\′}=\left[\begin{array}{c}\ln I_{p01}\\ M\\ \ln I_{p0K}\\ 0\\ M\\ 0\\ 0\\ M\\ 0\end{array}\right];$

步骤5：基于本发明建立的反演方程进行叠前反演，算出每个CRP点上每个时间采样点上的纵波阻抗、横波阻抗、密度等弹性参数；

步骤6：输出纵波阻抗、横波阻抗、纵横波速度比和密度等数据体；

步骤7：基于反演获得的弹性参数进行储层或流体的综合预测，并进行地质意义分析。

本发明基于贝叶斯参数估计理论和AVO基本原理，综合稀疏理论、剩余残差分布、地质意义等因素建立了混合罚函数方程，对叠前三参数反演过程进行正则化约束。模型与实际资料应用证实：混合范数正则化可进增加结果稳定性与地质意义，特别是使得相同偏移距和信噪比条件下的密度反演结果更加准确。本发明优势在于可直接应用于实际资料，获得合理的反演结果，提高了储层与流体预测精度。

附图说明

图1是二维地质模型(a,e,i)、反演初始值(b,f,j)、正演子波(c,d)及合成地震记录(g,h,k,l)示意图。

图2是无正则化约束条件下不同信噪比资料反演结果对比(模型值：a,e,i；无噪音下反演结果：b,f,j；S/N=10时反演结果：c,g,k；S/N=1时反演结果：d,h,l)示意图。

图3是基于Cauchy分布贝叶斯框架下不同信噪比资料反演结果对比(模型值：a,e,i；无噪音下反演结果：b,f,j；S/N=10时反演结果：c,g,k；S/N=1时反演结果：d,h,l)示意图。

图4是基于混合范数正则化约束条件下不同信噪比资料反演结果对比(模型值：a,e,i；无噪音下反演结果：b,f,j；S/N=10时反演结果：c,g,k；S/N=1时反演结果：d,h,l)示意图。

图5是CDP141位置处无正则化约束条件下不同信噪比资料反演结果对比(无噪音下反演结果：a,d,g；S/N=10时反演结果：b,e,h；S/N=1时反演结果：c,f,i)示意图。

图6是CDP141位置处基于Cauchy分布贝叶斯框架下不同信噪比资料反演结果对比(无噪音下反演结果：a,d,g；S/N=10时反演结果：b,e,h；S/N=1时反演结果：c,f,i)示意图。

图7是CDP141位置处基于混合范数正则化约束条件下不同信噪比资料反演结果对比(无噪音下反演结果：a,d,g；S/N=10时反演结果：b,e,h；S/N=1时反演结果：c,f,i)示意图。

图8是基于混合范数正则化的叠前三参数反演方法流程图。

图9是中古8井区某过井分角度叠加剖面（a-c）和密度初始值（d）示意图。

图10是中古8井区无正则化约束条件下的反演结果（a：纵波阻抗；b：纵横波速度比；c：密度）示意图。

图11是中古8井区基于Cauchy分布贝叶斯框架下的反演结果（a：纵波阻抗；b：纵横波速度比；c：密度）示意图。

图12是中古8井区基于混合范数正则化约束条件下的反演结果（a：纵波阻抗；b：纵横波速度比；c：密度）示意图。

具体实施方式

以下结合实例与附图说明本发明具体实施方式。

一种基于混合范数正则化的叠前三参数反演实现的储层与流体预测方法，包括如下步骤：

步骤1：对地震资料进行保幅处理，抽取高分辨率、高信噪比，高保真共反射点道集；

步骤2：基于共反射点道集获得分角度叠加数据体，结合测井数据估算各数据体地震子波；

步骤3：通过测井数据与地质层位建立纵波、横波与密度的初始模型，统计纵波与横波速度、纵波速度与密度之间的关系；

步骤4：基于估算的地震子波、建立的初始模型与弹性参数关系，采用混

合范数正则化方法，建立叠前AVO反演的混合范数正则化反演方程，

(G'^TG'+μQ'+βI)m'=βm₀'+G'^Td

其中，μ为Cauchy分布约束因子，β为初始模型约束因子，Ⅰ为单位矩阵，Q为对角矩阵：即 $Q^{\′}=\mathrm{diag}\left(\frac{1}{{(1+m_i^{\′2}/{\mathrm{\σ}}_{m^{\′}}^2)}^2}\right),$ $G^{\′}=\left[\begin{array}{ccc}\mathrm{WD}{A_1}^{\′}& \frac{1}{2}\mathrm{WD}{B}_1& \mathrm{WD}{C}_1\\ M& M& M\\ \mathrm{WD}{A_N}^{\′}& \frac{1}{2}\mathrm{WD}{B}_N& \mathrm{WD}{C}_N\end{array}\right];$ $m^{\′}=\left[\begin{array}{c}\ln I_{p1}\\ M\\ \ln I_{\mathrm{pK}}\\ \mathrm{\Δ}\ln I_{s1}\\ M\\ \mathrm{\Δ}\ln I_{\mathrm{sK}}\\ \mathrm{\Δ}\ln {\mathrm{\ρ}}_1\\ M\\ \mathrm{\Δ}\ln {\mathrm{\ρ}}_K\end{array}\right];{m_0}^{\′}=\left[\begin{array}{c}\ln I_{p01}\\ M\\ \ln I_{p0K}\\ 0\\ M\\ 0\\ 0\\ M\\ 0\end{array}\right];$

步骤5：基于所述反演方程进行叠前反演，算出每个共反射点上每个时间采样点上的纵波阻抗、横波阻抗、密度等弹性参数；

步骤6：输出纵波阻抗、横波阻抗、纵横波速度比和密度等数据体；

步骤7：基于反演获得的弹性参数进行储层或流体的综合预测，并进行地质意义分析。

其中的高分辨率、高信噪比，高保真是针对地震资料进行保幅处理后数据中相对的高分辨率、高信噪比，高保真。

实施例中，所述混合范数正则化方法，建立叠前AVO反演的混合范数正则化反演方程是：

步骤4中的正则化方程，是本发明综合稀疏理论、剩余残差分布、地质意义等因素建立了混合罚函数方程。包括稀疏性约束条件、反演结果与实际地震数据的匹配程度、空间连续性、先验弹性参数模型、不同弹性参数之间的内在联系等具体条件。通过步骤5中混合范数正则化约束条件对叠前反演的约束，可以获得更加符合地质意义的三参数反演结果，特别是使得相同偏移距和信噪比条件下的密度反演结果更加准确。

本发明所提供的基于混合范数正则化的叠前三参数反演方法求解不适定AVO反演问题的基本原理如下：

该方法综合考虑各种先验信息和地质条件以建立合理的反演方程，将实际AVO反演问题的目标函数表征为稀疏性、剩余残差和地质相关性三大类，建立综合目标函数为：

F(m)=αF_Contrast+λF_Seismic+βF_Geology   （6）

其中，α、λ和β为各项权重；稀疏性约束条件（F_Contrast）主要基于稀疏理论，即反射系数是叠加于高斯背景之上，因而过度稀疏的结果实际是噪音；剩余残差（F_Seismic）则描述了结果与输入地震数据之间的匹配程度；而地质条件的约束（F_Geology）则更多考虑和结果的实际意义，包括空间连续性、先验弹性参数模型、不同弹性参数之间的内在联系等。则混合范数的目标函数为：

$\begin{array}{c}F\left(m\right)=\mathrm{\α}{F}_{\mathrm{Contrast}}+\mathrm{\λ}{F}_{\mathrm{Seismic}}+{\mathrm{\β}}_1{F}_{\mathrm{trend}}+{\mathrm{\β}}_2{F}_{\mathrm{spatial}}+{\mathrm{\β}}_3{F}_{\mathrm{Gardner}}+{\mathrm{\β}}_4{F}_{\mathrm{Mudrock}}\\ =\mathrm{\α\Σ}\left|\right|\mathrm{\Δ}{m}_{\mathrm{Li}}\left|\right|+\mathrm{\λ\Σ}{\left|\right|{\mathrm{Gm}}_i-d_i\left|\right|}^2+{\mathrm{\β}}_1\mathrm{\Σ}\left|\right|m_i-m_{0i}\left|\right|+{\mathrm{\β}}_2\underset{i=x,y,z}{\mathrm{\Σ}}\left|\right|D_{\mathrm{im}}\left|\right|+{\mathrm{\β}}_3\mathrm{\Σ}\left|\right|\frac{\ln {\mathrm{\ρ}}_i-\ln {\mathrm{Ip}}_i}{1+\mathrm{Gard}}\left|\right|+{\mathrm{\β}}_4\mathrm{\Σ}\left|\right|\frac{\ln v_{\mathrm{si}}-\ln v_{\mathrm{pi}}}{\mathrm{Mudrock}}\left|\right|\end{array}---\left(7\right)$

其中，Δm_L表示弹性参数低通滤波后在纵向上的变化；m₀为初始值；D_im代表弹性参数在x,y,z三维空间三个方向上的变化量；Gard和Mudrock分别代表纵波速度和密度，纵波速度和横波速度的斜率（即ρ=a₁v_p ^Gard;Lv_s=Mudrock^*Lv_p+b₂）。

上式虽然涵盖了地震反演问题的所有必要约束条件，但无论对于从自变量出发还是求目标函数导数的何种算法，其计算都过于复杂。为此将混合范数目标函数项进行分类整理。首先，借助Hampson及Russell(2005)提出的矩阵转换思想，可以将纵波与密度关系和纵横波速度关系约束条件融入正演算子G中，借助纵波的高可解性能提高横波和密度反演精度。令纵波速度和密度，纵波速度和横波速度之间的关系为：

L_D=a₁L_p+a_d+△L_D

L_S=a₂L_P+a_s+△L_S   （8）

其中，式中的a₁，a_d，a₂和a_s可以通过分析先验测井信息获得，ΔL_D和ΔL_S分别表示岩石内流体的影响因子，也是反演需要求取的项，对围岩等于0。将一阶差分因子D引入上式，则有：

DL_D=a₁DL_p+D△L_D

DL_S=a₂DL_P+D△L_S   （9）

则公式（2）可以改写为：

$\begin{array}{c}d(t,\mathrm{\θ})=[\frac{A\left(\mathrm{\θ}\right)}{2}+a_2\frac{B\left(\mathrm{\θ}\right)}{2}+a_{1}C\left(\mathrm{\θ}\right)]W_{\mathrm{\θ}}g{\mathrm{DL}}_P+\frac{1}{2}B\left(\mathrm{\θ}\right)W_{\mathrm{\θ}}g{\mathrm{D\ΔL}}_S+C\left(\mathrm{\θ}\right)W_{\mathrm{\θ}}g{\mathrm{D\ΔL}}_D+n\\ =A^{\′}\left(\mathrm{\θ}\right)W_{\mathrm{\θ}}g{\mathrm{DL}}_P+\frac{1}{2}B\left(\mathrm{\θ}\right)W_{\mathrm{\θ}}g{\mathrm{D\ΔL}}_S+C\left(\mathrm{\θ}\right)W_{\mathrm{\θ}}g{\mathrm{D\ΔL}}_D+n\end{array},---\left(10\right)$

其中，A(θ)=A(θ)/2+α₁ ^*B(θ)/2+α₂ ^*C(θ)。据上式建立的反演方程已经采用了待求参数之间的关系进行约束，可以求得更加稳定的结果，但其各个参数扰动量仍存在正态分布假设，相当于对反演过程施加了一个均匀的正则化约束条件，难于在参数变化较大的边界处难以取得准确的反演结果。因而将公式（7）中的第一和第四项融合，参考边界保持约束反演方法在纵向和横向上加入扰动量符合柯西先验分布的假设（张赛民，2011；王康宁，2013等）以获得具有清楚边界的呈块状的解，并将基于测井数据的初始模型作为约束项加入目标函数以获得丰富的频率信息，则混合范数正则化下的叠前三参数反演的相应正则化方程为：

(G'^TG'+μQ'+βI)m'=βm₀'+G'Td   （11）

其中，μ为Cauchy分布约束因子，β为初始模型约束因子，Ⅰ为单位矩阵，Q为对角矩阵：即 $Q^{\′}=\mathrm{diag}\left(\frac{1}{{(1+m_i^{\′2}/{\mathrm{\σ}}_{m^{\′}}^2)}^2}\right),$ $G^{\′}=\left[\begin{array}{ccc}\mathrm{WD}{A_1}^{\′}& \frac{1}{2}\mathrm{WD}{B}_1& \mathrm{WD}{C}_1\\ M& M& M\\ \mathrm{WD}{A_N}^{\′}& \frac{1}{2}\mathrm{WD}{B}_N& \mathrm{WD}{C}_N\end{array}\right];$ $m^{\′}=\left[\begin{array}{c}\ln I_{p1}\\ M\\ \ln I_{\mathrm{pK}}\\ \mathrm{\Δ}\ln I_{s1}\\ M\\ \mathrm{\Δ}\ln I_{\mathrm{sK}}\\ \mathrm{\Δ}\ln {\mathrm{\ρ}}_1\\ M\\ \mathrm{\Δ}\ln {\mathrm{\ρ}}_K\end{array}\right];{m_0}^{\′}=\left[\begin{array}{c}\ln I_{p01}\\ M\\ \ln I_{p0K}\\ 0\\ M\\ 0\\ 0\\ M\\ 0\end{array}\right].$

通过对式（11）的求解，即可求得叠前弹性参数。加入约束信息的实际反演问题的非线性程度较弱，使用梯度类（如共轭梯度法）算法一般迭代三次即可求得最终解。

图1～图7用数值模型正演模拟，并据正演的资料分别对无正则化、不考虑参数相关性的正则化和提出的基于混合范数正则化叠前三参数反演方法进行了对比研究说明后者的有效性。仿照塔里木盆地次生碳酸盐岩储层特征设计了一个二维模型，如图1a，e，i所示；而后通过Zoeppritz方程正演计算该模型在入射角分别为5°、15°、25°、35°时的反射系数序列，选择主频为40Hz的雷克子波（图1c和d）分别正演形成四个角度的地震记录；最后基于合成的地震记录加入随机噪音分别形成信噪比为10和2的资料（图1g，k和h，i等）用于叠前反演对比分析。

分别采用无正则化约束条件（公式5去掉左端Cauchy因子），基于Cauchy分布正则化约束条件（公式5）、及发明的基于混合范数正则化（公式11）叠前三参数反演方法进行合成的不同信噪比模型资料的叠前反演，所得结果显示分别如图2-4所示。在实际地质条件下，纵波反射系数R_p，横波反射系数R_s，和密度梯度R_D三个未知量前面的系数数值常常呈现逐渐减小的趋势（特殊情况下，R_s的系数可能大于R_p），R_p的系数数值通常数倍于R_s，百倍于R_D；R_D的系数在实际地质情况下，入射角不大时常常近似为零。所以，对于实际含有噪音的资料而言，纵波信息最容易准确求解，横波其次；R_D及密度信息相对纵波和横波抗噪性最差，需要大偏移距资料求解。可以看出：无正则化约束条件下的反演结果精度最低，抗噪性和稳定性最差，各项反演结果随着资料噪音程度的增加而急剧变差；而加入Cauchy分布正则化约束条件后的三项反演结果抗噪性则大大增加，在不同信噪比条件下均能获得较为稳定的结果。然而，这两种方法获得的密度结果精度均不高（图2j-l和图3j-l），为初始模型（图1j）附近的少量扰动。

对比图2-4可以发现：基于混合范数正则化约束下的反演结果具备较好的抗噪性，特别是利用了纵波与横波，纵波与密度之间的关系对反演结果进行约束，获得了更为稳定的横波和密度反演结果。

图5-7为图2-4中从左至右第141个点处的各项反演结果对比。可以看出无正则化约束条件下的反演结果精度最低，抗噪性和稳定性最差，而加入Cauchy分布正则化约束条件后可以获得稳定的反演结果，但二者反演的密度结果精度均不高（图5和图6c，f，i）。基于混合范数正则化约束下的反演结果具备较好的抗噪性，并且反演的横波（图7b，g，h）和尤其是密度（图7c，e，h）结果更加准确。

图8是一种基于混合范数正则化的叠前三参数反演实现的储层与流体预测方法流程图：

步骤1：对地震资料进行保幅处理，抽取高分辨率、高信噪比，高保真的共反射点（CRP，Common reflection point）道集；

步骤2：基于高分辨率、高信噪比，高保真的CRP道集获得分角度叠加数据体，结合测井数据估算各数据体地震子波；

步骤3：通过测井数据与地质层位建立纵波、横波与密度的初始模型，统计纵波与横波速度、纵波速度与密度之间的关系；

步骤4：基于估算的子波、建立的初始模型与弹性参数关系，采用本发明的混合范数正则化方法，建立叠前AVO反演的混合范数正则化反演方程，

(G'^TG'+μQ'+βI)m'=βm₀'+G'^Td

其中，μ为Cauchy分布约束因子，β为初始模型约束因子，Ⅰ为单位矩阵，Q为对角矩阵：即 $Q^{\′}=\mathrm{diag}\left(\frac{1}{{(1+m_i^{\′2}/{\mathrm{\σ}}_{m^{\′}}^2)}^2}\right),$ $G^{\′}=\left[\begin{array}{ccc}\mathrm{WD}{A_1}^{\′}& \frac{1}{2}\mathrm{WD}{B}_1& \mathrm{WD}{C}_1\\ M& M& M\\ \mathrm{WD}{A_N}^{\′}& \frac{1}{2}\mathrm{WD}{B}_N& \mathrm{WD}{C}_N\end{array}\right];$ $m^{\′}=\left[\begin{array}{c}\ln I_{p1}\\ M\\ \ln I_{\mathrm{pK}}\\ \mathrm{\Δ}\ln I_{s1}\\ M\\ \mathrm{\Δ}\ln I_{\mathrm{sK}}\\ \mathrm{\Δ}\ln {\mathrm{\ρ}}_1\\ M\\ \mathrm{\Δ}\ln {\mathrm{\ρ}}_K\end{array}\right];{m_0}^{\′}=\left[\begin{array}{c}\ln I_{p01}\\ M\\ \ln I_{p0K}\\ 0\\ M\\ 0\\ 0\\ M\\ 0\end{array}\right];$

步骤5：基于本发明建立的反演方程进行叠前反演，算出每个CRP点上每个时间采样点上的纵波阻抗、横波阻抗、密度等弹性参数；

步骤6：输出纵波阻抗、横波阻抗、纵横波速度比和密度等数据体；

步骤7：基于反演获得的弹性参数进行储层或流体的综合预测，并进行地质意义分析。

图9a-c为过中古8井区典型井剖面的三个分角度叠加数据，可以看出在鹰山组（图上天蓝色曲线）存在非均质性分布极强的强反射，表现为由近致远振幅逐渐增加的第三类AVO特征，但是随着偏移距的增大资料品质有所降低。ZG111和ZG11井均在鹰山组钻遇了优质储层，如图上孔隙度曲线所示。特别地，为突出不同方法的差别，密度初始值选作常数，如图8d所示。

图10-12为分别采用无正则化约束条件，基于Cauchy分布正则化约束条件和发明的基于混合范数正则化叠前三参数反演方法反演的中古8井区实际结果对比。可以看出：三种反演方法获得的纵波结果相比其他横波，密度等参数更加稳定，差别最小，但是无正则化约束条件下的反演结果信噪比明显较低，特别是密度反演结果非常凌乱，与实际井区线完全不吻合；而基于Cauchy分布正则化约束和基于混合范数正则化约束下的反演结果抗噪性较好，纵波阻抗和纵横波速度比反演结果相差不大，但由于忽略了参数相关性，加上实际资料角度范围较小，前者反演的密度结果只是在初始模型附近的扰动，与实际井信息相差较大。显然，基于混合范数正则化约束下的反演结果具备较好的抗噪性，特别是利用了纵波与横波，纵波与密度之间的关系对反演结果进行约束，获得了更为稳定的横波和密度反演结果。

Claims (4)

1.一种基于各向异性有效场的泥页岩岩石物理模型方法，其特征在于，包括以下步骤：

步骤1，通过测井分析，得到弹性张量计算所需的岩性和基本物性参数，包括岩石基质、泥质含量、孔隙度、有机质含量；

步骤2，通过Voigt-Reuss-Hill平均公式（1952）得到岩石骨架的体积模量和剪切模量，通过Wood方程（1955）得到混合流体的弹性参数，并根据各组分的密度和体积分数求算饱和流体介质的总密度；

步骤3，利用Kuster-Toksoz模型，给定干酪根成熟度初值或者依照化学分析的数据，将油气与干酪根混合，得到干酪根与油气混合物的弹性张量；

步骤4，将干酪根油气混合物作为背景介质，给定岩石基质的纵横比（范围0～1），利用各向异性有效场理论计算基质与干酪根混合物的复合岩石的弹性张量；

步骤5，利用各向异性有效场理论将粒内孔隙、粒间孔隙、裂缝加入到混合岩石中形成干岩石，计算其弹性张量参数；

步骤6，通过Brown-Korringa（1975）公式进行流体替代，计算饱和流体介质的弹性张量参数。

2.根据权利要求1所述的方法，其特征在于，步骤3所述的干酪根与油气混合物的弹性张量通过Kuster-Toksoz模型进行计算，设S为流体的饱和度，也即为干酪根的成熟度，且S=φ_f/(φ_f+φ_kerogen)，φ_f和φ_kerogen分别表示流体与干酪根的体积分数，且混合物的弹性张量为：

$\frac{c_{13}^{\mathrm{if}}+\frac{2}{3}{c}_{55}^{\mathrm{if}}}{K_k}=\frac{1+[4{\mathrm{\μ}}_k(K_f-K_k)/(3K_f+4{\mathrm{\μ}}_k)/K_k]S}{1-[3(K_f-K_k)/(3K_f+4{\mathrm{\μ}}_k)]S}---\left(1\right)$

$\frac{c_{55}^{\mathrm{if}}}{{\mathrm{\μ}}_k}=\frac{(1-S)(9K_f+8{\mathrm{\μ}}_k)}{9K_f+8{\mathrm{\μ}}_k+S(6K_k+12{\mathrm{\μ}}_k)}---\left(2\right)$ 其中符号k，f，if分别表示干酪根，流体，干酪根与油气混合物；K，μ，S，c分别表示体积模量，剪切模量，流体饱和度即干酪根成熟度，复合岩石的弹性张量。

3.根据权利要求1所述的方法，其特征在于，步骤4所述各向异性有效场理论通过以下方式计算复合介质的有效弹性张量：

$C_{\mathrm{ijkl}}^{\mathrm{eff}}=C_{\mathrm{ijkl}}^{\left(0\right)}+v[{\left(C_{\mathrm{ijkl}}^{\left(1\right)}\right)}^{-1}+(1-v)P_{\mathrm{ijkl}}{]}^{-1}---\left(3\right)$ 其中是介质的有效弹性张量，是基质的弹性张量，是包含物与基质弹性张量差异，v是包含物的体积分数，且这里的G(x)是各项异性无限大介质的格林函数，且符号(ij)代表包含的符号ij的对称。

4.根据权利要求1所述的方法，其特征在于，步骤5所述的粒间孔隙、粒内孔隙、裂缝，利用不同的孔隙纵横比0.8、0.15、0.01分别表征。