CN105572727A - 基于孔隙流体参数频变反演的储层流体识别方法 - Google Patents

Abstract

基于孔隙流体参数频变反演的储层流体识别方法，包括对叠前地震数据进行保幅、去噪及偏移处理；进行分角度叠加得到至少三个角度的分角度叠加地震资料；建立孔隙流体参数以及剪切模量模型，提取角度子波，对不同角度的叠加地震资料进行小波多尺度分频处理；进行孔隙流体参数的频变属性反演；根据孔隙流体参数频变属性演变值高低判断储层流体发育情况。本方法更多地利用接收的地震信息，能够较好地解决常规流体识别中的假象问题，孔隙流体敏感参数的频变属性对储层流体的敏感性更高，对流体识别的精度更高，并且利用贝叶斯反演理论进行反演，稳定性强，可信度高。

Description

基于孔隙流体参数频变反演的储层流体识别方法

技术领域

本发明用于地球物理勘探资料处理方法领域，特别是一种利用孔隙流体参数频散属性进行砂岩、泥页岩、碳酸盐岩和火成岩等各类岩石储层流体识别的方法。

背景技术

获取高精度的储层流体识别结果是地球物理勘探的最终目标。随着勘探程度的加深，我们所面临的地下条件也变得日益复杂，以地震资料为主体进行储层流体识别可以在一定程度上增加勘探成功率，提高生产效率以及减少勘探开发成本，是现阶段油气勘探领域的研究热点之一。

地球物理勘探中所接收到的反射地震波是地下地层的综合响应，反映了地下储层的物性、流体等发育情况。常规的流体识别方法是通过地下地层的介质弹性参数(如纵波阻抗、拉梅参数、泊松比等)进行储层预测以及流体识别。基于地震资料的储层流体识别是一项在岩石物理理论指导下将与孔隙流体有关的异常特性表征为流体因子，利用流体因子实现介质孔隙流体类型判别的技术。20世纪七八十年代，利用“亮点”现象进行储层预测的技术得到广泛应用。1987年，Smith与Gidlow率先提出了流体因子的概念，即利用纵横波速度相对变化量的加权叠加判别储层含烃异常。Goodway等(1997)提出了LMR技术，使用拉梅常数和密度进行流体指示的方法。George(2003)提出了流体因子角和交会图角两种流体因子概念。Gassmann(1951)在一定的假设前提下，推导了表征双相介质弹性模量参数的Gassmann方程。Russell等人(2003)在多孔弹性介质岩石物理理论的指导下，基于Biot-Gassmann方程推导出了可以反映孔隙流体类型的ρf参数，并于2006年进一步研究了孔隙流体参数(Gassmann流体项f)，指出可以将f直接作为一项流体因子参与流体检测。印兴耀等(2010)提出了包含Gassmann流体项的弹性阻抗方程，可以直接反演高精度的流体项。这些方法仅仅只是利用了地震波的振幅信息，而接收的地震波中不仅包含了振幅信息，还包含了地震波的频率以及相位信息。实际观测数据与岩石物理实验均证实，当地震波穿过含烃储层时，在地震频带内会发生较为明显振幅频散现象，主要是因为介质弹性模量在不同的频率响应下具有不同数值，从而振幅在不同的频率相应下也具有不同的数值。因此，利用振幅随频率变化的信息有助于我们进一步明确储层流体的发育状况，对勘探开发油气藏方面具有理论意义和实用价值，是储层流体识别技术发展的必然趋势。考虑到地震波经过油气储集层后表现出吸收衰减、速度频散等异常现象，可以利用地震资料中所蕴含的频率信息进行储层描述与流体识别。利用低频阴影识别以及优势频带识别储层流体是目前利用与频率有关的属性进行储层描述的主要方法。Taner(1979)在通过复地震道分析提取瞬时频率和瞬时振幅等属性，并在研究中发现了与储层相关的“低频阴影”现象。谱分解方法是低频阴影检测的重要手段。Castagna(2003)通过瞬时谱分析，进而检测与储层相关的低频阴影。陈学华等(2009，2011)利用黏滞型波动方程对饱含流体地质模型进行低频阴影的数值模拟，发现它能刻画由于流体弥散性、黏滞性以及速度频散引起的反射同相轴时移、相位畸变、主频降低和振幅衰减及低频阴影现象。张固澜等(2010)基于改进的广义S变换检测与油气相关的低频阴影，提高了储层纵横向位置的准确度，直接指示油气的存在。张波等(2011)利用时频连续小波变换，结合低频阴影对储层进行了含气性检测。但造成“低频阴影”原因有很多，Ebrom(2004)总结了叠前和叠后地震资料中的可能造成“低频阴影”的10种机理，可见利用低频阴影检测油气不具有普适性。Chen等(2008)利用谱分解技术对大角度角道集的深部储层流体的频率响应的研究，发现在低频剖面中能明显区分油水界面。Goloshubin(2000，2002，2006)讨论了储层流体饱和的低频效应，并展示了利用低频信息进行流体识别的实例。陈学华等(2012)认为在不同角度叠加道集中，含流体地层地震反射中不同频率成分存在差异，通过广义S变换进行谱分解，再利用低频分量构建瞬时谱差异信息提取公式进而进行流体识别，与实际钻井结果相符。上述利用频率信息进行储层流体识别的方法主要是定性的，要实现利用频率属性进行定量储层流体识别，最为关键的问题是频散属性与地震反射特征之间的数学关系表示。Cooper(1967)推导了线性粘弹性介质界面处的反射与透射系数；Krebes(1984)和Nechtschein等人(1997)则推导了粘弹性介质的反射与透射系数；Ursin(2002)推导了粘弹各向同性薄层介质的反射和透射系数；Sidler和Carcione(2007)分析了VTI粘弹介质中反射和透射特征；而Ren(2009)结合中观White模型，利用数值模型的方法研究了频散界面处的法向入射振幅随频率的变化规律，并将其分为三种类型，这对AVF(AmplitudeVariationwithFrequency)的研究提供了理论指导，以上这些研究主要是以基于各种假设的理论分析为主，围绕实际资料实现频散属性提取的难度较大，直到2009年Wilson根据Chapman岩石物理模型Smith与Gidlow两项AVO反射系数近似公式看作与频率的函数，给出了纵、横波频散属性表达式，并利用了最小二乘方法进行反演，定量地利用频散属性进行流体识别才得以实现。吴小羊(2010)、SunZD(2012)分别基于伪平滑Wigner-Ville分布的频谱分析技术和广义S变换技术实现了频变AVO反演，得到纵横波速度的频散属性，得到较好的储层流体识别效果。ZhangSX(2011)等推导了基于AVO梯度频散属性的频变AVO近似方程，并利用最小二乘方法进行了反演。

虽然以上方法能够在一定程度上识别流体，但仍然存在流体识别的假象。不同的弹性参数的频率响应特征也不相同，而速度频散属性并不是对储层流体敏感性最强的流体因子,并且最小二乘反演方法缺乏稳定性。因此，在常规的AVO近似方程中加入频率信息，可以更好地利用接收得到的资料，结合测井信息，研究储层流体识别方法，对减少流体识别假象具有重要意义。

发明内容

本发明的目的是针对实际生产中现存技术的不足，给出了一种基于孔隙流体参数频变反演的储层流体识别方法，旨在岩石物理分析不同弹性参数的频散属性对储层流体的敏感性，提出将对储层流体更为敏感的孔隙流体参数频散属性作为流体因子，从而推导基于孔隙流体参数频散属性的频变反射系数近似方程，在贝叶斯反演理论指导下加入模型约束，增强反演的稳定性，以达到提高储层流体识别精度、减少流体识别假象的问题。

本发明的总体技术路线是利用地震岩石物理频散特征分析，从所有参数中选择频散最为严重的弹性参数，反演其频散属性，进行相应的储层预测和流体识别。利用平面反射方程，从孔隙流体参数AVO近似方程出发，考虑孔隙流体参数的频变特征，推导包含弹性参数频散程度和频散程度随偏移距变化梯度的两个频散属性的地震波AVO近似频变反射特征公式；利用角度叠加地震数据，分别进行小波多尺度分解，基于模型约束利用贝叶斯反演方法直接反演孔隙流体参数的频散属性，用于储层预测和流体识别。其中：根据地震岩石物理理论分析了在含流体相同情况下岩石不同弹性参数的频变特征，选取了频散最为严重(对流体最为敏感)的弹性参数的频散属性用于储层预测和流体识别。考虑孔隙流体参数的频变特征，从Russell孔隙流体参数近似公式出发，推导了包含弹性参数频散程度和频散程度随偏移距变化梯度的两个频散属性的地震波AVO近似频变反射特征公式。

基于孔隙流体参数频变反演的储层流体识别方法，包括：

(1)对叠前地震数据进行保幅、去噪及偏移处理，即将叠前CMP道集数据转化叠前CRP道集数据；

(2)在对叠前地震数据进行保幅、去噪及偏移处理的基础上，进行分角度叠加得到至少三个角度的分角度叠加地震资料，即根据叠前CRP道集提取至少三个角度叠加道集数据；

(3)利用测井曲线计算孔隙流体参数以及剪切模量的伪井曲线，并建立孔隙流体参数以及剪切模量模型，根据不同角度的叠加地震资料以及测井资料提取角度子波，对不同角度的叠加地震资料进行小波多尺度分频处理，转化为不同主频的地震资料；

(4)结合不同频率不同角度的地震资料以及贝叶斯最大后验概率反演方法，进行孔隙流体参数的频变属性反演；

(5)根据孔隙流体参数频变属性演变值高低判断储层流体发育情况。

上述方案还包括：

所述步骤(4)中，考虑孔隙流体参数的频变特征，从Russell孔隙流体参数近似公式出发，得出包含弹性参数频散程度和频散程度随偏移距变化梯度的两个频散属性的地震波AVO近似频变反射特征公式。

所述近似频变反射特征公式如下：

ΔR(θ,ω)＝(ω-ω₀)a(θ)D_f+(ω-ω₀)b(θ)D_μ(1)

其中，θ为入射角的角度，ω代表频率，ω₀为参考频率； D_f为孔隙流体参数的频散属性，表达式为D_μ为剪切模量的频散属性，表达式为ΔR(θ,ω)为入射角为θ、频率为ω时的反射系数与参考频率ω₀下的反射系数之差，表达式为ΔR(θ,ω)＝R(θ,ω)-R(θ,ω₀)，R(θ,ω)表示入射角为θ、频率为ω时的反射系数。

上述方案进一步包括：

所述步骤(4)中结合不同频率不同角度的地震资料以及贝叶斯最大后验概率反演方法后，反演目标函数加入模型约束进行孔隙流体参数的频变属性反演，最终目标函数如(2)式所示：

$F\left(m\right)={\left|\right|\mathrm{Gm}-d\left|\right|}_2^2+\mathrm{\λQm}+{\mathrm{\α}}_1{\left|\right|m_{D_f}-{\mathrm{\ξ}}_{D_f}\left|\right|}_2^2+{\mathrm{\α}}_2{\left|\right|m_{D_{\mathrm{\μ}}}-{\mathrm{\ξ}}_{D_{\mathrm{\μ}}}\left|\right|}_2^2---\left(2\right)$

通过对上式进行求导取极值可得最终反演方程：

$(G^{T}G+\mathrm{\λQ}+{\mathrm{\α}}_{1}I+{\mathrm{\α}}_{2}I)m=G^{T}d+{\left[{\mathrm{\α}}_1{{\mathrm{\ξ}}_{D_f}}^{\text{T}}{\mathrm{\α}}_2{{\mathrm{\ξ}}_{D_{\mathrm{\μ}}}}^T\right]}^T---\left(3\right)$

(2)与(3)式中，G为正演算子，m为待反演参数，有孔隙流体参数频散属性D_f以及剪切模量频散属性D_μ构成的列向量，d为经过分角度叠加、分频处理后的地震数据组成的列向量，Q为柯西约束矩阵，为待反演的孔隙流体参数频散属性向量，m_Dμ为待反演剪切模量频散属性向量，与分别为孔隙流体参数频散属性与剪切模量频散属性模型向量，G,d，m以及Q的具体表达式如下：

$G=\left(\begin{array}{c}(W_1^1{\mathrm{\ω}}_1-W_1^0{\mathrm{\ω}}_0)\left[\begin{array}{cc}{A}_1& B_1\end{array}\right]\\ (W_2^1{\mathrm{\ω}}_2-W_2^0{\mathrm{\ω}}_0)\left[\begin{array}{cc}{A}_1& B_1\end{array}\right]\\ \·\\ \·\\ \·\\ (W_N^1{\mathrm{\ω}}_N-W_N^0{\mathrm{\ω}}_0)\left[\begin{array}{cc}{A}_1& B_1\end{array}\right]\\ \·\\ \·\\ \·\\ (W_j^p{\mathrm{\ω}}_j-W_j^0{\mathrm{\ω}}_0)\left[\begin{array}{cc}{A}_p& B_p\end{array}\right]\\ \·\\ \·\\ \·\\ (W_1^n{\mathrm{\ω}}_1-W_N^0{\mathrm{\ω}}_0)\left[\begin{array}{cc}{A}_n& B_n\end{array}\right]\\ (W_2^n{\mathrm{\ω}}_2-W_2^0{\mathrm{\ω}}_0)\left[\begin{array}{cc}{A}_n& B_n\end{array}\right]\\ \·\\ \·\\ \·\\ (W_N^n{\mathrm{\ω}}_N-W_N^0{\mathrm{\ω}}_0)\left[\begin{array}{cc}{A}_n& B_n\end{array}\right]\end{array}\right]d=\left[\begin{array}{c}\mathrm{\Δ}{d}_1^1\left({\mathrm{\θ}}_1\right)\\ \mathrm{\Δ}{d}_2^1\left({\mathrm{\θ}}_1\right)\\ \·\\ \·\\ \·\\ \mathrm{\Δ}{d}_i^j\left({\mathrm{\θ}}_p\right)\\ \·\\ \·\\ \·\\ \mathrm{\Δ}{d}_M^N\left({\mathrm{\θ}}_n\right)\end{array}\right]m=\left[\begin{array}{c}{D}_f\\ D_{\mathrm{\μ}}\end{array}\right]---\left(4\right)$

上式中，i表示时间采样点，j表示为频率，p表示角度，N表示小波分频的个数，M表示时间采样点的个数，

则最终反演结果表示为：

$m={(G^{T}G+\mathrm{\λQ}+{\mathrm{\α}}_{1}I+{\mathrm{\α}}_{2}I)}^{-1}(G^{T}d+{\left[{\mathrm{\α}}_1{{\mathrm{\ξ}}_{D_f}}^T{\mathrm{\α}}_2{{\mathrm{\ξ}}_{D_{\mathrm{\μ}}}}^T\right]}^T)---\left(6\right)$

公式(2)中柯西约束项矩阵Q是与待反演参数m的函数，公式(2)通过迭代最小二乘方法进行求解，提高利用(6)式计算反演得到反演结果的精度。

所述迭代最小二乘反演步骤为：

步骤1：给定初始值m₀，迭代次数i＝0，m_i＝m₀；

步骤2：令m＝m_i，

步骤3：利用m计算柯西约束矩阵Q；

步骤4：利用(6)式计算反演得到反演结果m_i

步骤5：计算精度若不满足要求，返回步骤2，若满足精度，输出结果。

本发明首先根据地震岩石物理分析了各弹性参数频散属性对储层流体的敏感性，确立了孔隙流体参数频散属性作为流体因子用于检测储层含油气性的可行性。基于弹性模量对变量变化特征推导了基于频变孔隙流体参数的反射特征近似方程，然后结合角度部分叠加地震数据和测井资料在贝叶斯理论框架下反演孔隙流体参数的频变特征，最终将反演得到的孔隙流体参数频变属性用于储层流体识别。孔隙流体参数频散属性的高值对应于储层流体发育的位置。与常规地震反演方法相比该储层流体识别技术的优势体现在：不仅仅利用了接收到的反射地震振幅随偏移距的变化，还应用了反射振幅随频率的变化信息，孔隙流体敏感参数的频变属性对储层流体的敏感性更高，对流体识别的精度更高，并且利用贝叶斯反演理论进行反演，稳定性强，可信度高。

本发明提出的基于孔隙流体参数频变反演的储层流体识别方法更多地利用接收的地震信息，能够较好地解决常规流体识别中的假象问题。

附图说明

图1基于贝叶斯方法反演孔隙流体参数频散属性流程图

图2不同弹性参数频散属性标准化比较

图3三个角度部分叠加数据剖面：(a)小角度，(b)中角度，(c)大角度

图4短时傅立叶变换与小波变换时频窗口比较

图5小角度地震数据A井井旁道频谱分析

图6小角度地震数据分频(a)10Hz，(b)15Hz，(c)20Hz，(d)25Hz，(e)30Hz

图7孔隙流体参数频散属性预测储层含油气性(高值代表储层油气发育)

图8孔隙流体参数频散属性预测储层含油气平面分布(高值代表储层油气发育良好)

图9本发明一种实施例的识别流程图

具体实施方式

下面结合附图和具体实施例对本发明的技术方案做进一步说明。

如图9所示，基于孔隙流体参数频变反演的储层流体识别方法主要分为以下几个步骤：

(1)将叠前CMP道集数据转化叠前CRP道集数据，再根据叠前CRP道集提取至少3个部分角度叠加数据；

(2)结合角度部分叠加数据和地震数据多尺度分解方法，将各角度的地震数据转化为5个不同主频的地震资料；

(3)结合不同频率不同角度的地震资料以及贝叶斯理论反演方法可以反演孔隙流体参数的频变属性；

(4)孔隙流体参数频变属性的高值代表了储层流体的发育良好，图7为孔隙流体参数频散属性预测储层含油气性，预测气层发育情况与测井解释结果相一致；图8为孔隙流体参数频散属性预测储层含油气平面分布图。

下面通过优化的具体实施例做更进一步说明。

参照图9和图1所示。

步骤1：

根据地震岩石物理分析、计算各介质弹性参数(纵波阻抗I_p、纵波速度V_p、泊松比σ、体积模量K、拉梅参数λ、孔隙流体参数f)的频变属性，为了能够更好地进行对比，对各个频变弹性参数进行标准化处理，将低频段数值校正为1，如图2所示。各参数可近似为线性频散，从图2中可以得到如下关系式：

$\frac{\mathrm{df}}{\mathrm{d\ω}}>\frac{\mathrm{d\λ}}{\mathrm{d\ω}}>\frac{\mathrm{dK}}{\mathrm{d\ω}}>\frac{\mathrm{d\σ}}{\mathrm{d\ω}}>\frac{d{V}_p}{\mathrm{d\ω}}=\frac{d{I}_p}{\mathrm{d\ω}}---\left(1\right)$

从图2中可以看出，各弹性参数都存在频散特征，其中纵波速度频散曲线与纵波阻抗频散曲线重合，从侧面表明了岩石的密度不存在频散现象，在众多的频散曲线中，孔隙流体参数的频散程度最大，对储层流体的敏感性最强，可用于储层流体识别，提高识别精度。

步骤2：

基于孔隙流体参数的Zoeppritz方程近似式如：

$R_{\mathrm{pp}}\left(\mathrm{\θ}\right)\left[(\frac{1}{4}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{4{\mathrm{\γ}}_{\mathrm{sat}}^2}){\sec }^2\mathrm{\θ}\right]\frac{\mathrm{\Δf}}{f}+[\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{4{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}]\frac{\mathrm{\Δ\μ}}{\mathrm{\μ}}+[\frac{1}{2}-\frac{{\sec }^2\mathrm{\θ}}{4}]\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}---\left(2\right)$

上式中，θ为入射角，γ_dry与γ_sat分别为干岩石与饱和岩石的纵横波速比，R_pp(θ)为纵波勘探中入射角为θ时的反射系数。

由于孔隙流体参数f以及剪切模量均是与频率相关的函数，密度项ρ一般假设不发生频散，因而，反射系数R_pp也是与频率有关的函数，因此，我们可以将(2)式写作频率的函数：

$R_{\mathrm{pp}}(\mathrm{\θ},\mathrm{\ω})\≈a\left(\mathrm{\θ}\right)\frac{\mathrm{\Δf}}{f}\left(\mathrm{\ω}\right)+b\left(\mathrm{\θ}\right)\frac{\mathrm{\Δ\μ}}{\mathrm{\μ}}\left(\mathrm{\ω}\right)+c\left(\mathrm{\θ}\right)\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}---\left(3\right)$

上式中，ω表示频率，a(θ)、b(θ)与c(θ)分别为：

$a\left(\mathrm{\θ}\right)=(\frac{1}{4}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{4{\mathrm{\γ}}_{\mathrm{sat}}^2}){\sec }^2\mathrm{\θ}$ $b\left(\mathrm{\θ}\right)=\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{4{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}$ $c\left(\mathrm{\θ}\right)=\frac{1}{2}-\frac{{\sec }^2\mathrm{\θ}}{4}$

从图2中可以看出，频散曲线在频散频带内可近似为线性函数，因而可以对与在频散频带内进行泰勒展开并去一阶近似项，设参考频率为ω₀，

$\frac{\mathrm{\Δf}}{f}\left(\mathrm{\ω}\right)=\frac{\mathrm{\Δf}}{f}\left({\mathrm{\ω}}_0\right)+(\mathrm{\ω}-{\mathrm{\ω}}_0)\frac{\∂}{\∂\mathrm{\ω}}\left[\frac{\mathrm{\Δf}}{f}\left(\mathrm{\ω}\right)\right]---\left(4\right)$

$\frac{\mathrm{\Δ\μ}}{\mathrm{\μ}}\left(\mathrm{\ω}\right)=\frac{\mathrm{\Δ\μ}}{\mathrm{\μ}}\left({\mathrm{\ω}}_0\right)+(\mathrm{\ω}-{\mathrm{\ω}}_0)\frac{\∂}{\∂\mathrm{\ω}}\left[\frac{\mathrm{\Δ\μ}}{\mathrm{\μ}}\left(\mathrm{\ω}\right)\right]---\left(5\right)$

(4)与(5)两式中分别为孔隙流体参数频散属性与剪切模量频散属性，分别即为D_f、D_μ，将(4)、(5)两式代入(3)可以得到

$\begin{array}{c}{R}_{\mathrm{pp}}(\mathrm{\θ},\mathrm{\ω})\≈a\left(\mathrm{\θ}\right)\{\frac{\mathrm{\Δf}}{f}\left({\mathrm{\ω}}_0\right)+(\mathrm{\ω}-{\mathrm{\ω}}_0)D_f\}\\ +b\left(\mathrm{\θ}\right)\{\frac{\mathrm{\Δ\μ}}{\mathrm{\μ}}\left({\mathrm{\ω}}_0\right)+(\mathrm{\ω}-{\mathrm{\ω}}_0)D_{\mathrm{\μ}}\}+c\left(\text{\θ}\right)\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}\end{array}---\left(6\right)$

令频率为ω₀时的反射系数为R_pp0(θ,ω₀)：

$R_{\mathrm{pp}0}(\mathrm{\θ},{\mathrm{\ω}}_0)=a\left(\mathrm{\θ}\right)\frac{\mathrm{\Δf}}{f}\left({\mathrm{\ω}}_0\right)+b\left(\mathrm{\θ}\right)\frac{\mathrm{\Δ\μ}}{\mathrm{\μ}}\left({\mathrm{\ω}}_0\right)+c\left(\mathrm{\θ}\right)\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}---\left(7\right)$

(7)式代入到(6)是可以得到：

$\begin{array}{c}\mathrm{\Δ}{R}_{\mathrm{pp}}(\mathrm{\θ},\mathrm{\ω})=R_{\mathrm{pp}}(\mathrm{\θ},\mathrm{\ω})-R_{\mathrm{pp}}(\mathrm{\θ},{\mathrm{\ω}}_0)\\ =a\left(\mathrm{\θ}\right)(\mathrm{\ω}-{\mathrm{\ω}}_0)D_f+b\left(\mathrm{\θ}\right)(\mathrm{\ω}-{\mathrm{\ω}}_0)D_{\mathrm{\μ}}\end{array}---\left(8\right)$

公式(8)即为本发明推导的包含孔隙流体参数频散属性的频变反射系数方程。

步骤3：

由公式(8)可知，反演孔隙流体参数频散属性D_f之前需要获得不同角度的不同频率的地震信息。为了提高地震数据的信噪比进而增强反演稳定性，首先对叠前道集进行分角度叠加，分角度部分叠加数据如图3所示，图中椭圆中部井位上深色部分为测井解释气层。在分角度叠加地震数据的基础之上进行分频处理。高频率信号的分辨率比低频信号的分辨高，因而频率越高窗口应越小。然而如图4a所示，短时傅立叶变换窗口大小不随频率变化，并不满足这一要求，为了适应窗口大小随频率的变化，选择如图4b所示的连续小波变换。针对目的层的频谱关系选择合适尺度的Ricker小波，通过分频处理将地震记录分解为不同频带的地震信息，与常规的梯形窄带滤波器分频方法相比，基于小波变换的分频方法具有一定的自适应分辨分析特性，可以将地震记录分解为一系列具有质心频率的窄带剖面，实现了不同频带范围的地震振幅信息分离，较好的克服了常规梯形滤波器分频方法对信号的较强正弦调谐影响。对地震数据进行分频处理需要在地震数据的有效频带内进行，并要选择合适的参考频率。因此，首先对图3a小角度地震数据井旁道进行频谱分析，如图5所示，图5a为井旁道地震数据，图5b为利用连续小波变换方法对地震数据的井旁道地震记录时频分析图，图5c为井旁道的频谱，从图中可以得到井旁道的有效频带范围为7Hz-35Hz，根据时频图以及振幅谱选择优势频率25Hz作为参考频率ω₀，进一步在有效频带范围选择10Hz、15Hz、20Hz、25Hz与30Hz对分角度叠加的地震数据进行小波分频处理得到相应频带范围的地震数据，如图6所示。利用同样的方法对中角度和大角度部分叠加地震数据进行分频处理。

步骤4：

进一步考虑基本的褶积模型假设，将不同频带地震记录对应的子波信息(此时的子波信息是原始子波与相应频带Ricker小波的褶积结果)以子波矩阵的形式引入公式(8)，最终得到如下形式的方程

上式中，i表示时间采样点，j表示为频率，p表示角度，N表示小波分频的个数，M表示时间采样点的个数；A_p与B_p分别是由a(θ_p)与b(θ_p)组成的对斜对角矩阵，矩阵大小为M×M。

为了更为稳定的反演得到孔隙流体参数频散属性，本发明采用基于贝叶斯反演框架下的模型约束反演方法。一般情况下，地震记录中的噪音信息服从正态分布，且相互独立，即贝叶斯理论中的似然函数是高斯分布，为了保护地下的弱反射信息，在此假设反映地下的弹性参数频散属性的先验信息服从柯西分布，在此假设待反演参数孔隙流体参数频散属性D_f与剪切模量频散属性D_μ是相互独立的，以单变量柯西分布进行正则先验约束，最终利用贝叶斯公式结合似然函数与先验函数，得到包含反演参数m的后验概率函数：

$P\left(m\right|d)=\∝K_1{K}_2\·\exp \left[\frac{-{(\mathrm{Gm}-d)}^T(\mathrm{Gm}-d)}{2{\mathrm{\σ}}^2}\right]\·\underset{i=1}{\overset{2*M}{\mathrm{\Π}}}\left[\frac{1}{1+m_i^{\′2}/{\mathrm{\σ}}_m^2}\right]---\left(10\right)$

上式中，σ表示地震数据中噪音的方差；σ_m表示代反演参数的方差；通过求取最大后验概率，可以得到目标函数：

$F\left(m\right)={\left|\right|\mathrm{Gm}-d\left|\right|}_2^2+\mathrm{\λQm}---\left(11\right)$

其中，λ为稀疏约束系数，Q为斜对角加权矩阵，定义如下：

进一步考虑反演稳定性问题，加入模型约束，得到最终的目标函数如(12)所示：

$F\left(m\right)={\left|\right|\mathrm{Gm}-d\left|\right|}_2^2+\mathrm{\λQm}+{\mathrm{\α}}_1{\left|\right|m_{D_f}-{\mathrm{\ξ}}_{D_f}\left|\right|}_2^2+{\mathrm{\α}}_2{\left|\right|m_{D_{\mathrm{\μ}}}-{\mathrm{\ξ}}_{D_{\mathrm{\μ}}}\left|\right|}_2^2---\left(12\right)$

上式中，与为孔隙流体参数频散属性与剪切模量的频散属性模型，α₁与α₂为模型约束系数，对(12)式进行求导取零值，可以得到反演方程：

$m={(G^{T}G+\mathrm{\λQ}+{\mathrm{\α}}_{1}I+{\mathrm{\α}}_{2}I)}^{-1}(G^{T}d+{\left[{\mathrm{\α}}_1{{\mathrm{\ξ}}_{D_f}}^T{\mathrm{\α}}_2{{\mathrm{\ξ}}_{D_{\mathrm{\μ}}}}^T\right]}^T)---\left(13\right)$

公式(13)中柯西约束项矩阵Q是与待反演参数m的函数，并且公式(13)具有较弱的非线性，通过迭代最小二乘方法进行求解可较稳定、准确地获取孔隙流体参数频散属性，利用提取的孔隙流体参数频散属性进行储层流体识别。

迭代最小二乘反演步骤为：

步骤1：给定初始值m₀，迭代次数i＝0，m_i＝m₀；

步骤2：令m＝m_i，

步骤3：利用m计算柯西约束矩阵Q；

步骤4：利用(13)式计算反演得到反演结果m_i

步骤5：计算精度若不满足要求，返回步骤2，若满足精度，输出结果。

利用上述方法对图3的地震数据进行反演得到孔隙流体参数频散属性，如图7所示，由岩石物理分析可知，图7中的高值指示了储层流体发育良好。如图8所示，从反演得到的三维孔隙流体参数频散属性体提取沿层切片可以识别和预测储层含油气性的平面展布。

以上便是本发明基于孔隙流体参数频变反演的储层流体识别方法以及应用,对实际地震资料处理的具体流程如图9。

通过上述实施例进一步说明，该方法首先对地震资料进行分角度叠加处理，并利用测井曲线计算孔隙流体参数以及剪切模量的伪井曲线，并建立孔隙流体参数以及剪切模量模型，利用角度部分叠加道集提取角度子波，然后对分角度叠加地震数据进行小波多尺度分频处理，最后在模型约束条件下利用贝叶斯最大后验概率反演方法进行孔隙流体参数频散属性反演。利用小波多尺度分频方法能够自适应窗口随频率增高而缩小的要求。地震波穿过含烃储层会发生较为明显的频散现象，因此，孔隙流体参数频散属性反演结果的高值指示储层流体含油气。

Claims (6)

1.基于孔隙流体参数频变反演的储层流体识别方法，其特征在于包括：

(1)对叠前地震数据进行保幅、去噪及偏移处理；

(2)在对叠前地震数据进行保幅、去噪及偏移处理的基础上，进行分角度叠加得到至少三个角度的分角度叠加地震资料；

(3)利用测井曲线计算孔隙流体参数以及剪切模量的伪井曲线，并建立孔隙流体参数以及剪切模量模型，根据不同角度的叠加地震资料以及测井资料提取角度子波，对不同角度的叠加地震资料进行小波多尺度分频处理，转化为不同主频的地震资料；

(4)结合不同频率不同角度的地震资料以及贝叶斯最大后验概率反演方法，进行孔隙流体参数的频变属性反演；

(5)根据孔隙流体参数频变属性演变值高低判断储层流体发育情况。

2.根据权利要求1所述的基于孔隙流体参数频变反演的储层流体识别方法，其特征在于：所述步骤(4)中，考虑孔隙流体参数的频变特征，从Russell孔隙流体参数近似公式出发，得出包含弹性参数频散程度和频散程度随偏移距变化梯度的两个频散属性的地震波AVO近似频变反射特征公式。

3.根据权利要求2所述的基于孔隙流体参数频变反演的储层流体识别方法，其特征在于：所述近似频变反射特征公式如下：

$\mathrm{\ΔR}(\mathrm{\θ},\mathrm{\ω})=(\mathrm{\ω}-{\mathrm{\ω}}_0)a\left(\mathrm{\θ}\right)D_f+(\mathrm{\ω}-{\mathrm{\ω}}_0)b\left(\mathrm{\θ}\right)D_{\mathrm{\μ}}---\left(1\right)$

其中，θ为入射角的角度，ω代表频率，ω₀为参考频率； $b\left(\mathrm{\θ}\right)=\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{4{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ};$ D_f为孔隙流体参数的频散属性，表达式为 $D_f=\frac{\∂}{\∂\mathrm{\ω}}\left(\frac{\mathrm{\Δf}}{f}\right),$ D_μ为剪切模量的频散属性，表达式为ΔR(θ,ω)为入射角为θ、频率为ω时的反射系数与参考频率ω₀下的反射系数之差，表达式为ΔR(θ,ω)＝R(θ,ω)-R(θ,ω₀)，R(θ,ω)表示入射角为θ、频率为ω时的反射系数。

4.根据权利要求1所述的基于孔隙流体参数频变反演的储层流体识别方法，其特征在于：所述步骤(4)中结合不同频率不同角度的地震资料以及贝叶斯最大后验概率反演方法后，反演目标函数加入模型约束进行孔隙流体参数的频变属性反演，最终目标函数如(2)式所示：

$F\left(m\right)={\left|\right|\mathrm{Gm}-d\left|\right|}_2^2+\mathrm{\λQm}+{\mathrm{\α}}_1{\left|\right|m_{D_f}-{\mathrm{\ξ}}_{D_f}\left|\right|}_2^2+{\mathrm{\α}}_2{\left|\right|m_{D_{\mathrm{\μ}}}-{\mathrm{\ξ}}_{D_{\mathrm{\μ}}}\left|\right|}_2^2---\left(2\right)$

通过对上式进行求导取极值可得最终反演方程：

$(G^{T}G+\mathrm{\λQ}+{\mathrm{\α}}_{1}I+{\mathrm{\α}}_{2}I)m=G^{T}d+{\left[\begin{array}{cc}{\mathrm{\α}}_1{{\mathrm{\ξ}}_{D_f}}^T& {\mathrm{\α}}_2{{\mathrm{\ξ}}_{D_{\mathrm{\μ}}}}^T\end{array}\right]}^T---\left(3\right)$

(2)与(3)式中，G为正演算子，m为待反演参数，有孔隙流体参数频散属性D_f以及剪切模量频散属性D_μ构成的列向量，d为经过分角度叠加、分频处理后的地震数据组成的列向量，Q为柯西约束矩阵，为待反演的孔隙流体参数频散属性向量，m_Dμ为待反演剪切模量频散属性向量，与分别为孔隙流体参数频散属性与剪切模量频散属性模型向量，G,d，m以及Q的具体表达式如下：

$G=\left[\begin{array}{c}({W_1}^1{\mathrm{\ω}}_1-{W_1}^0{\mathrm{\ω}}_0)\left[\begin{array}{cc}{A}_1& B_1\end{array}\right]\\ W_2^1{\mathrm{\ω}}_2-W_2^0{\mathrm{\ω}}_0\left[\begin{array}{cc}{A}_1& B_1\end{array}\right]\\ .\\ .\\ .\\ (W_N^1{\mathrm{\ω}}_N-W_N^0{\mathrm{\ω}}_0)\left[\begin{array}{cc}{A}_1& B_1\end{array}\right]\\ .\\ .\\ .\\ (W_j^p{\mathrm{\ω}}_j-W_j^0{\mathrm{\ω}}_0)\left[\begin{array}{cc}{A}_p& B_p\end{array}\right]\\ .\\ .\\ .\\ (W_1^n{\mathrm{\ω}}_1-W_N^0{\mathrm{\ω}}_0)\left[\begin{array}{cc}{A}_n& B_n\end{array}\right]\\ (W_2^n{\mathrm{\ω}}_2-W_2^0{\mathrm{\ω}}_0)\left[\begin{array}{cc}{A}_n& B_n\end{array}\right]\\ .\\ .\\ .\\ (W_N^n{\mathrm{\ω}}_N-W_N^0{\mathrm{\ω}}_0)\left[\begin{array}{cc}{A}_n& B_n\end{array}\right]\end{array}\right]$ $d=\left[\begin{array}{c}{\mathrm{\Δd}}_1^1\left({\mathrm{\θ}}_1\right)\\ {\mathrm{\Δd}}_2^1\left({\mathrm{\θ}}_1\right)\\ .\\ .\\ .\\ {\mathrm{\Δd}}_i^j\left({\mathrm{\θ}}_p\right)\\ .\\ .\\ .\\ {\mathrm{\Δd}}_M^N\left({\mathrm{\θ}}_n\right)\end{array}\right]$ $m=\left[\begin{array}{c}{D}_f\\ D_{\mathrm{\μ}}\end{array}\right]---\left(4\right)$

上式中，i表示时间采样点，j表示为频率，p表示角度，N表示小波分频的个数，M表示时间采样点的个数，

则最终反演结果表示为：

$m={(G^{T}G+\mathrm{\λQ}+{\mathrm{\α}}_{1}I+{\mathrm{\α}}_{2}I)}^{-1}(G^{T}d+{\left[\begin{array}{cc}{\mathrm{\α}}_1{{\mathrm{\ξ}}_{D_f}}^T& {\mathrm{\α}}_2{{\mathrm{\ξ}}_{\mathrm{\μ}}}^T\end{array}\right]}^T)---\left(6\right)$

5.根据权利要求4所述的基于孔隙流体参数频变反演的储层流体识别方法，其特征在于：公式(2)中柯西约束项矩阵Q是与待反演参数m的函数，公式(2)通过迭代最小二乘方法进行求解，提高利用(6)式计算反演得到反演结果的精度。

6.根据权利要求5所述的基于孔隙流体参数频变反演的储层流体识别方法，其特征在于迭代最小二乘反演步骤为：

步骤1：给定初始值m₀，迭代次数i＝0，m_i＝m₀；

步骤2：令m＝m_i，

步骤3：利用m计算柯西约束矩阵Q；

步骤4：利用(6)式计算反演得到反演结果m_i

步骤5：计算精度若不满足要求，返回步骤2，若满足精度，输出结果。