随着油气勘探开发难度加大,如何从地面地震数据中获取反映储层流体特征参数成为储层地球物理的研究热点.叠前地震反演作为获取储层流体敏感参数的重要途径促进了储层流体识别技术的发展[1-3].Smith和Gidlow 率先提出用叠前地震数据加权叠加得到流体因子和伪泊松比剖面来预测岩性和流体[4].Goodway 等提出了用于流体识别的LMR(Lambda-Mu-Rho,拉梅参数反演)技术[5].但考虑到LMR 技术在实际流体饱和多孔介质中的局限性,Russell等人在前人研究[6-10]的基础上,基于Biot-Gassmann方程得到流体饱和条件下的流体因子,该流体因子可通过叠前地震反演或AVO(AmplitudeVariationwithOffset,振幅随偏移距变化)分析获取纵横波阻抗、密度信息后再间接计算获得[11].然而,实际叠前地震资料特别是深部地震资料受偏移距限制,无法获得对密度信息起决定作用[12]的远道集信息,也就无法通过叠前地震反演或AVO 分析获得合理的密度信息.
叠前地震反演主要致力于反演地下介质弹性参数来指导含油气储层预测,许多学者在利用叠前地震资料反演弹性参数方面做了大量的工作.Fatti[13]用加权叠加的方法得到了纵波阻抗的相对变化量和横波阻抗的相对变化量,进而得到纵横波阻抗体.Simmons和Backus[14]应用线性化的反演方法来求取纵波阻抗、横波阻抗和密度反射系数,反演过程中用到了Gardner经验方程.Buland和Omre[15]研究了关于贝叶斯线性AVO 反演的方法,直接反演纵波速度、横波速度和密度参数,不足之处是该方法未考虑参数之间的关系.Connolly[16]提出弹性阻抗的概念,并有学者认为弹性阻抗反演的抗噪能力优于AVO 反演[17].印兴耀等[18]提出利用Connolly弹性阻抗方程先反演纵、横波速度和密度,进而可以得到纵、横波阻抗,拉梅常数,泊松比等多种弹性参数.王保丽等推导了基于Gray 的弹性阻抗方程,实现了拉梅参数直接提取[19].弹性阻抗反演作为一种比较实用的叠前反演技术,使AVO 信息通过可以直观理解的方式显示出来,在地震解释和岩石属性分析中发挥了重要作用.本文首先基于Biot-Gassmann方程对纵横波速度方程进行了重新推导,基于流体饱和多孔介质模型建立了流体因子与纵横波模量之间的直接关系.新的关系式表明,流体因子仅与纵横波模量信息直接相关,避免了常规求取需要的密度信息无法准确求取的问题.然后由Aki-Richard近似出发,推导得到新的基于纵横波模量的Zeoppritz近似公式,并基于此公式探讨了基于叠前角度叠加地震道集的纵横波模量弹性阻抗反演方法.该方法减少了常规方法先反演纵横波速度和密度再计算纵横波模量带来的累积误差,最后模型和实际应用表明,该方法能够获得合理的流体因子.
2 流体因子与纵横波模量之间的关系均匀各向同性完全弹性无孔介质中纵横波速度方程为
(1) |
(2) |
其中,Vp 为纵波速度,Vs 为横波速度,M为纵波模量,μ 为横波模量,ρ 为介质密度.
而实际地下储层是饱含流体的多孔介质,方程(1)和(2)并不能反映地下介质的真实情况,因此Biot[8]和Gassmann[9]设立了理想的流体饱和多孔介质模型,建立了纵横波速度与拉梅参数的关系,Russell[11]在此基础上建立流体因子与纵横波阻抗、密度的关系如(3)式所示:
(3) |
其中,f为Russell流体因子,Zp、Zs 分别为纵横波阻抗,ρ 为密度,c为参数因子.
笔者基于Biot和Gassmann设立的流体饱和多孔介质模型建立Russell流体因子与纵横波模量之间的关系.纵横波速度关系式可以表示为
(4) |
(5) |
(6) |
式中,Mdry 为干燥岩石的纵波模量,β 为Biot系数,Km 为岩石骨架的体积模量,Km 为孔隙流体的体积模量,Φ 为岩石孔隙度,ρsat 为饱和流体岩石密度.
对比公式(1)和(4),将公式(4)中的β2N记作f,联立公式(4)和(5)可以得到
(7) |
其中,M为流体饱和岩石纵波模量,μ 为流体饱和岩石横波模量,c= (Vp/Vs)dry2 = Kdry/μ +4/3,Murphy[20]测量了大量不同孔隙度下纯石英砂岩的Kdry/μ 值,其测量结果表明Kdry/μ 的平均值为0.9,相应的c值为2.233.如果Kdry/μ 的取值在1.0 附近,对应c值为2.333.实际上,应根据实际工区情况确定最优参数c.
3 基于弹性阻抗的纵横波模量直接反演 3.1 基于纵横波模量的Zeoppritz近似公式1976年,Richards和Frasier研究了性质相近的反射场半空间之间的反射和透射问题,给出了以速度和密度相对变化表示的反射系数近似公式.1980年,Aki和Richards在《定量地震学》经典专著中对Richards 和Frasier 等近似进行了综合整理[21],给出了如下近似公式:
(8) |
其中,α、β 为纵横波速度,ρ 为密度,θ 为入射角,Δα/α、Δβ/β、Δρ/ρ 分别是纵横波速度反射系数和密度反射系数.纵横波模量反射系数与纵横波速度、密度反射系数的关系如下所示:
(9) |
(10) |
将式(9)、(10)代入式(8)可得基于纵横波模量的Zeoppritz近似方程为
(11) |
式中,ΔM/M=2(M1-M2)/(M1+M2),为纵波模量反射系数,Δμ/μ=2(μ1-μ2)/(μ1+μ2),为横波模量反射系数,Δρ/ρ=2(ρ1-ρ2)/(ρ1+ρ2),为密度反射系数,其中,M1、M2、μ1、μ2、ρ1、ρ2 分别是反射界面两侧介质的纵横波模量和密度.
3.2 基于纵横波模量的弹性阻抗方程推导用波阻抗的对数值表示反射系数:
(12) |
式中:EI是弹性阻抗.
由公式(11)得:
(13) |
用K表示β2/α2,并用Δln(x)来替换Δx/x重新整理可得:
(14) |
最后取积分并指数化,可得:
(15) |
(15)式可写成下面的形式:
(16) |
其中:
上面推导出的弹性阻抗(EI(θ))随角度在量纲尺度上有很大变化,不利于进行不同角度的EI(θ)值之间及其与波阻抗(AI)值的对比.为了消除入射角变化对量纲尺度的影响,类比Whitcombe[22]标准化方法,引入3个参考常数M0、μ0 和ρ0,这些常数被定为M、μ 和ρ曲线的平均值,同时用因子(M0ρ0)1/2进一步标定函数,使EI(θ)的量纲尺度变得与AI一样,可以得到基于纵横波模量近似的弹性阻抗公式标准化形式为
(17) |
显然,由(17)式可以得到当θ =0°时,弹性阻抗为常数αρ,即声阻抗值.因此,θ变化时,标准化后的EI(θ)的维数保持为常数,EI函数就不会随着θ 剧变,从而实现不同的EI值间的直观对比,且对于所有的角度θ,EI值变回到常规的AI范围内,克服了不同角度弹性阻抗量纲不统一的不足.
3.4 反射系数对比分析利用Ostrander[23]所构造的含气砂岩模型和Goodway[5]根据实测资料给出的含气砂岩模型,我们对基于纵横波模量的近似方程和弹性阻抗方程进行定量分析,模型参数分别见表 1和表 2.
入射介质波阻抗小于透射介质波阻抗的界面称为正波阻抗界面,反之称为负波阻抗界面.上述两种模型中,上覆页岩、下伏含油气砂岩的反射界面就是负波阻抗界面,而上覆砂岩、下伏页岩的反射界面为正波阻抗界面.以上述模型为基础,分别用精确的Zoeppritz方程、Aki-Richards近似方程、基于纵横波模量的近似方程和基于纵横波模量的弹性阻抗方程计算了不同界面处的反射系数(图 1),由图可知,无论Ostrander所构造的含气砂岩模型还是Goodway根据实测资料给出的含气砂岩模型,基于纵横波模量的近似方程和基于纵横波模量的弹性阻抗方程在入射角为50°左右都与精确Zoeppritz方程有较好的近似,能够适用于大角度入射情况下反射系数的准确求解.
纵横波模量弹性阻抗反演主要包括地震资料处理、测井资料处理、角度子波的提取、弹性阻抗反演和纵横波模量直接提取等,具体反演流程如图 2所示.与常规弹性阻抗反演类似,地震资料处理主要是将地震数据的偏移距数据体转化为部分角度叠加数据体,部分角度叠加数据体考虑振幅随角度变化AVA(Amplitude Variation with Angle)特征,比传统的叠加剖面提供了更多的振幅信息,有利于烃类信息的直接检测;测井曲线分析与处理主要是利用弹性阻抗方程和已知井曲线信息计算井旁道伪弹性阻抗曲线,该曲线除用来作为约束外,还可以弥补地震波传播过程中损失的频率成分.利用提取的各角度子波,通过稀疏脉冲反演可以得到各角度弹性阻抗体.
从EI(θ)中直接提取纵横波模量和密度等弹性参数是叠前反演中重要的一环.我们采用直接提取的方法以减小间接计算带来的累计误差,将方程(10)两边取对数,令A0 = (M0ρ0)$\frac{1}{2}$ 并考虑在角度相同的情况下,同一岩石物性参数
(18) |
通过井旁道弹性阻抗值和已知井纵横波模量和密度计算得到9 个常系数a(θ1)、b(θ1)、c(θ1);a(θ2)、b(θ2)、c(θ2);a(θ3)、b(θ3)、c(θ3).将它们联合反演所得的各角度弹性阻抗体代入方程组(18),从而可获得各道任意一个采样点处的M、μ、ρ.
3.6 模型试算笔者通过模型(图 3)验证基于纵横波模量弹性阻抗方程的反演方法的可行性及其直接提取在抗噪方面的优势.图 3a、图 3b和图 3c分别为纵波模量、横波模量和密度参数的二维地层模型剖面,在1.3s左右有一气层;在2.1s左右之间有一油层.首先由纵横波模量和密度模型剖面正演并抽取获得10°、20°和30°部分角度叠加剖面,随后采用上述基于纵横波模量的弹性阻抗反演方法对3个部分角度叠加剖面进行反演,得到各角度弹性阻抗剖面(如图 4),其中,图 4a、图 4b和图 4c分别为10°、20°和30°弹性阻抗剖面.从图中可以看出,该方法可靠性较强,能够得到合理的不同角度弹性阻抗.
从部分角度叠加剖面和弹性阻抗反演结果可以看出气层有较好的亮点特征,而气层下方,同相轴出现明显下拉现象,且信噪比较低,影响了对油层的识别.我们分别采用常规先提取纵横波速度再间接计算和基于纵横波模量弹性阻抗方程直接提取纵横波模量、密度参数的方法,提取得到的纵横波模量参数剖面如图 5—6所示.其中,图 5a和5b分别是间接计算和直接提取获得的纵波模量参数剖面,图 6a和图 6b分别是间接计算和直接提取获得的横波模量参数剖面,从图中黑框区域可以看出,直接提取获得的纵横波模量剖面的信噪比明显高于间接计算结果.为更精确地分析直接提取方法的优势,我们抽取图 5和图 6中的过气层的单道进行对比分析,如图 7—8所示.图 7a曲线是模型数据与间接计算得到的纵波模量曲线对比,图 7b曲线是模型数据与直接提取得到的纵波模量曲线对比,从该图可以看出,相比间接计算结果,直接提取结果与模型数据吻合程度更高,特别是在气层下方油层附近具有明显的抗噪能力.图 8a曲线是模型数据与间接计算得到的横波模量曲线对比,图 8b曲线是模型数据与直接提取得到的横波模量曲线对比,同样可以得到与图 7 类似的认识.
实际资料来自中国东部某勘探工区,进行叠前反演之前,需要对地震数据进行保幅处理,包括精细的波前扩散补偿、震源组合与检波器组合效应的校正、反Q滤波、地表一致性处理、叠前去噪处理、去除多次波等,并假设处理后的层间多次波、各向异性的影响可以忽略不计.
笔者利用不同角度部分叠加地震剖面实现了基于弹性阻抗的纵横波模量直接反演,图 9a、9b和9c分别是反演得到的过井A 的纵波模量、横波模量和密度剖面.图中的黑色曲线为电阻率测井曲线,井上1.24s附近为含气砂岩层,1.26s和1.3s为含油砂岩层.其中,纵波模量表征介质抵抗正应变的能力,体现介质的抗压缩性和硬度,与骨架和流体性质有关,可以大致区分储层的流体类型,在含油和含气位置处显示低值,其中含气储层位置处更低;相比非储层位置,含油气储层处的孔隙度更大,构成储层的岩石抗剪切能力减弱,因此横波模量和密度含油气储层处也呈低值显示,但三者在1.3s油层处显示较弱甚至无显示.根据本文提出的流体因子计算公式(7),可得流体因子剖面如图 10 所示,其中,根据本区岩石物理和测井数据统计结果,我们取c=2.238.从图 10中可以看出,流体因子对储层流体类型更为敏感,特别在1.3s油层位置处,流体因子剖面有更好的指示.
基于Biot-Gassmann方程建立了流体饱和多孔介质中流体因子与纵横波模量的直接关系,避免流体因子计算需要的密度信息无法准确求取的问题.推导得到新的基于纵横波模量的Zeoppritz近似公式,同时充分利用弹性阻抗反演的优势,基于此公式探讨了基于叠前角度叠加地震道集的纵横波模量弹性阻抗反演方法,进而得到指示储层流体的流体因子,减少了常规纵横波速度和密度弹性阻抗反演间接计算纵横波模量带来的累积误差,模型测试与实际资料表明,该方法能够获得合理有效的流体因子,提供了一种精度更高的基于叠前反演的流体识别方法.
[1] | Gidlow M, Smith G C. The Fluid Factor Angle. EAGE 65th Conference & Exhibition , 2003: 2-5. |
[2] | Smith G C. A comparison of the fluid factor with λ and μ in AVO analysis. Expanded Abstracts of 70th SEG Mtg , 2000: 1940-1945. |
[3] | Quakenbush M, Shang B, Tuttle C. Poisson impedance. The Leading Edge , 2006, 25(2): 128-138. DOI:10.1190/1.2172301 |
[4] | Smith G C, Gidlow P M. Weighted stacking for rock property estimation and detection of gas. Geophys. Prosp. , 1987, 35(9): 993-1014. |
[5] | Goodway B, Chen T W, Jon Downton. Improved AVO fluid detection and lithology discrimination using Lamé petrophysical parameters; "λρ","μρ", &"λ/μ fluid stack", from P and S inversions. Expanded Abstracts of 67th SEG Mtg , 1997: 183-186. |
[6] | Hedlin K. Pore space modulus and extraction using AVO. Expanded Abstracts of 70th SEG Mtg , 2000: 170-173. |
[7] | Hilterman F J. Seismic Amplitude Interpretation. Short Soc. Expl. Geophys, Distinguished Instructor Series 4 , 2001. |
[8] | Biot M A. General theory of three-dimensional consolidation. J. Appl. Physics , 1941, 12(2): 155-164. |
[9] | Gassmann F. Uber die Elastizitat poroser Medien. Vier. Natur Gesellschaft , 1951, 96: 1-23. |
[10] | Garat J, Krief M, Stellingwerff J, et al. A petrophysical interpretation using the velocities of P and S waves. The Log Analyst , 1990, 31(6): 355-369. |
[11] | Russell B H, Hedlin K, Hilterman F J, Lines L R. Fluid-property discrimination with AVO: A Biot-Gassmann perspective. Geophysics , 2003, 68(1): 29-39. DOI:10.1190/1.1543192 |
[12] | Downton J E. Seismic parameter estimation from AVO inversion. Calgary: University of Calgary , 2005. |
[13] | Fatti J L, Smith G C, Vail P J, et al. Detection of gas in sandstone reservoirs using AVO analysis: A 3-D seismic case history using the Geostack technique. Geophysics , 1994, 59(9): 1362-1376. DOI:10.1190/1.1443695 |
[14] | Simmons J L Jr, Backus M M. Waveform-based AVO inversion and AVO prediction-error. Geophysics , 1996, 61(6): 1575-1588. DOI:10.1190/1.1444077 |
[15] | Buland A, Omre H. Bayesian linearized AVO inversion. Geophysics , 2003, 68(1): 185-198. DOI:10.1190/1.1543206 |
[16] | Connolly P. Elastic impedance. The Leading Edge , 1999, 18(4): 438-452. DOI:10.1190/1.1438307 |
[17] | Cambois G. AVO inversion and elastic impedance. Expanded Abstracts of 70th SEG Mtg , 2000: 1953-1956. |
[18] | 印兴耀, 袁世洪, 张繁昌. 从弹性波阻抗中提取岩石物性参数 // 中国石油学会物探专业委员会及美国地球物理学家学会. CPS/SEG 2004国际地球物理会议论文集. 北京: CPS/ SEG 2004国际地球物理会议, 2004: 538-542. Y.in X Y, Yuan S H, Zhang F C. Rock elastic parameters calculated from elastic impedance (in Chinese). Expanded Abstracts CPS/SEG Mtg, 2004: 538-542. |
[19] | Wang B L, Yin X Y, Zhang F C. Lamé parameters inversion based on elastic impedance and its application. Applied Geophysics , 2006, 3(3): 174-178. DOI:10.1007/s11770-006-0026-z |
[20] | Murphy W, Reischer A, Hsu K. Modulus decomposition of compressional and shear velocities in sand bodies. Geophysics , 1993, 58(2): 227-239. DOI:10.1190/1.1443408 |
[21] | Aki K, Richards P G. Quantitative Seismology: Theory and Methods. San Francisco: W H Freeman and Co Cambridge , 1980: 144-154. |
[22] | Whitcombe D N. Elastic impedance normalization. Geophysics , 2002, 67(1): 60-62. DOI:10.1190/1.1451331 |
[23] | Ostrander W J. Plane-wave reflection coefficients for gas sands at nonnormal angles of incidence. Geophysics , 1984, 49(10): 1637-1648. DOI:10.1190/1.1441571 |