2. 中国石油化工股份有限公司华北分公司, 郑州 450006
2. China Petroleum & Chemical Corporation, Huabei Branch, Zhengzhou 450006, China
AVO理论在过去的30多年里得到长足发展与应用.自从1979年Taner等人(1979)发现与油气有关的低频阴影以来,学者进行大量研究来解释该现象,时频分析与谱分解方法得以深入研究和广泛应用.Castagna等(2003)使用谱分解的方法,将“低频阴影”作为烃类指示因子用于烃类检测.Ebrom(2004)分析了出现“低频阴影”现象的原因,包括地层衰减因素.Chapman等(2006)利用喷射流衰减机制分析了高低频极限与频散时的地震记录,并正演模拟了低频阴影.Odebeatu等(2006)用频谱分解技术对地震剖面进行分频,发现了低频衰减现象可以指示油气,并指出Gassmann理论和衰减理论在不同条件下都会有低频衰减现象.陈学华等(2009)使用黏弹性双相介质的假设,基于波动方程对低频阴影进行了数值模拟,利用广义S变换谱分解检测到了实际资料油气储层的低频阴影.Chen等(2012)通过对物理模型(Wang et al., 2010)的实验模拟所得到的地震数据进行时频分析,验证了第一类AVO的高频亮点现象(Chapman et al., 2005),说明了高频亮点可以用来识别砂泥岩薄互层的气藏.
针对速度频散与地震反射特征之间的关系,一些学者在衰减介质反射系数(Cooper and Henry, 1964;Krebes,1984;Nechtschein and Hron, 1997;Ursin and Stovas, 2002)、利用频变反射特征识别储层(Xu et al., 2011)与裂缝(Maultzsch et al., 2003;Goloshubin and Silin, 2005;Wang and Li, 2006; Wang et al., 2007;Chapman,2009;Al-Harrasi et al., 2011)等方面进行了研究.Chapman等(2006)指出频散介质与弹性介质界面处反射与频率有关.Ren等(2009)讨论了垂直入射时,非频散介质与频散介质界面处的反射系数振幅与相位随频率的变化关系,并总结了三类AVF(Amplitude Variation with Frequency)响应.基于该研究,Liu等(2011)基于Biot孔隙弹性理论模拟了三类AVF响应,结果与理论预测相符.Innanen(2011)推导了吸收反射系 数公式,进而研究了提取品质因子的AVF/AVA方法.
还有学者利用速度频散来定量刻画储层与流体,即频变AVO(Frequency-dependent AVO)反演.室内实验与野外观测已证实了速度频散与地层物性密切相关(Sun,2009),而且频散对地层物性的变化比速度更敏感(Johnston,2005).与流体有关的频散通常产生频变AVO(FAVO)效应.叠后反演与叠前反演得到的频变地震属性如纵波频散属性,已经证实为有效的碳氢指示因子.Wilson等(2010)首先提出了基于Smith-Gidlow近似的频散属性计算方法,建立了基于岩石物理理论、谱分解与AVO分析的FAVO反演框架.吴小羊等(2010)利用Wigner-Ville分布将FAVO反演应用于实际数据.Zhang等(2011)推导了基于纵波频散程度与梯度的频变反射系数近似.程冰洁等(2012)构造了一系列FAVO属性.
然而,以上的频变AVO属性均基于纵波速度频散提出,反映的是介质的整体弹性信息.而基于双相介质理论得到的Gassmann流体项f能够更好地表征孔隙流体的特征(印兴耀等,2013a),将其作为流体因子已经取得了较好的储层流体识别效果(印兴耀等,2010;Zong et al., 2013;印兴耀等,2013b).为此本文基于f-μ-ρ近似(Russell et al., 2011)构建了频变流体项属性.实际上,岩石物理分析已经表明,相比纵波速度,Gassmann流体项对频率更敏感.本文介绍了频变流体项属性的反演流程,并通过模型与实际数据验证其可行性. 2 反演方法 2.1 方程构建
Russell等(2011)提出了AVO反射系数近似
其中f、μ和ρ分别表示界面两侧介质的Gassmann流体项,剪切模量和密度的平均值;Δf、Δμ和Δρ则分别表示界面两侧的Gassmann流体项,剪切模量和密度的差异;γsat和γdry是界面两侧,饱和岩石与干岩石的纵横波速度比的平均值;θ是入射角.忽略密度随频率的变化,可得由于存在速度频散,VP和VS与频率有关,可知ΔVP/VP和ΔVS/VS均与频率有关,进而Δf/f和Δμ/μ也与频率有关.因此我们考虑(1)式的近似频变形式,即关于频率ω进行一阶泰勒展开
其中ω0代表地震记录的主频,ΔR(θ,ω)=RPP(θ,ω)- RPP(θ,ω0).Df=(Δf/f)/ω和Dμ=(Δμ/μ)/ω 分别表示Gassmann流体项(即孔隙流体参数)的频散程度与剪切模量的频散程度,即待反演的频散属性参数.对于两个入射角(θ1和θ2),两个采样点(以上标1和2表示)以及两个频率(ω1和ω2)情形,反演方程为 其中,和表示入射角为θj的第i个采样点的正演算子.此方程可以推广到m个入射角,n个采样点及l个频率的情况.将矩阵进行块化表示,得到 其中,Δ R iii= [ΔR1(θii,ωii)ΔR2(θii,ωii)…ΔRn(θii,ωii)]T(i=1,2,…,l; ii=1,2,…,m)表示的是 第ii个偏移距,频率为ωi所对应的共n个采样点的数据列向量; A iii和 B iii(i=1,2,…,l;ii=1,2,…,m)是第ii个偏移距,频率为ωi所对应斜对角系数矩阵,即为了消除子波的带通滤波作用,从地震振幅信息得到属性的真实值,需要两边引入不同频率对应的子波矩阵,
向量Diii表示入射角为θiii时,主频ωiii剖面振幅与主频ω0剖面振幅的差异. Wi是用ωi主频子波构建的子波矩阵.考虑到子波矩阵的稀疏性,本文采用共轭梯度法求解,它适用于稀疏矩阵线性方程组,克服了最速下降法收敛慢的缺点,又避免了牛顿法需要存 储和计算Hessian矩阵并求逆的缺点.程序框图如图 1.考虑到剪切模量一般与孔隙流体无关,实际当中采用Df来进行流体识别.上述反演方程需要不同频带的地震信息,这需要对叠前数据进行谱分解.本文采用连续小波变换(CWT),因为它具有灵活的时频分辨率调节能力,同时容易实现.具体分频方法如下:
CWT定义式为
给定,在假设地震信号为实信号前提下,(10)式可写作褶积形式这样通过将地震数据与特定时间尺度算子褶积,就可以得到对应频带剖面.实际计算机实现过程中需要将尺度因子a离散化,即表示为a=a0j,从而
其中a0=(ωmin/ωmax)1/Nscale,ωmin、ωmax和Nscale分别 为所需不同频带剖面的最小频率、最大频率和频带个数.由尺度因子含义易知,ω=ω0/a=ω0/(ω0/ωmax)j/Nscale. 因此,为了得到某频带f剖面,我们需将每道地震数据同ωmin、ωmax和Nscale所决定的Morlet小波矩阵的第j列相褶积. 2.3 反演优化 2.3.1 参数去相关与协方差矩阵的建立由于纵横波速度和密度之间的统计相关性,导致流体项频散和剪切模量频散亦存在相关性,因此需要应用两者之间的协方差矩阵对参数进行去相关处理(杨培杰,2008).协方差矩阵的建立一般有三种方法:基于井估计、基于岩石物理关系计算和人为给定.这里用最简单的方法:从附近井资料进行统计估计:
其中:X =[Df Dμ ],N为样本个数.建立协方差矩阵之前首先要由井资料计算流体项和剪切模量,再对井作分频处理,即(i)根据地震数据时间深度范围,对时深转换后包含时间采样点、纵波速度、横波速度和密度的井数据进行时间重采样.并计算流体项f和剪切模量μ.
(ii)根据叠后地震剖面的振幅谱,确定主频ω0及附近两侧两个频率ω-= ω0-Δω与ω+= ω0+Δω.Δω的确定根据频带宽度选择5~10 Hz.利用2.2节的分频方法对井分频,得到主频分别为ω-与ω+的流体项f和剪切模量μ,即f-、μ-、f+ 和μ+.
(iii)最终由井资料估计的Df=(f+- f-)/2Δω,Dμ=(μ+-μ-)/2Δω.
已知Df和Dμ,那么待反演的参数之间的协方差矩阵 Cr 表示为:
将协方差矩阵延伸为N个时间采样点的情况,产生2N×2N的稀疏协方差矩阵 Cx,特征向量分析也可扩展为2N×2N的情况.令
则 经过如下变换 得到 经过变换后的参数的协方差矩阵 C ′ r 等于(20)式非对角线上的元素都为零,即经过变换后的参数之间相互独立. 2.3.2 稀疏先验约束
通常反演是病态的,这是因为子波矩阵是雅克比矩阵,且本文频变AVO反演利用多频带地震资料,必须对大型正演矩阵求逆.因此需要引入先验约束以使求解稳定,反演优化流程图如图 2.首先由已知的VP、VS和ρ井数据,通过上述的井分频处理得到Df和Dμ的井约束信息.由于实际地层中VP、VS和ρ存在统计相关性,因而需要对Df和Dμ进行去相关处理.下一步根据贝叶斯理论将先验信息与似然函数相结合,得到正则化的目标函数.最后通过求解新的目标函数以及反去相关处理,得到最终的频散属性.
假设反演前地震数据已经过去多次波、球面扩散补偿、衰减补偿等处理.为了验证方法的可行性与可靠性,在此设计了三层模型,参数选择参照Chapman等(2006)的模型,如图 3.顶层和底层代表泥 岩,纵波速度、横波速度和密度分别为2743 m·s-1、 1394 m·s-1、2060 kg·m-3.中间层左侧代表水饱和砂岩,中间层右侧代表气饱和砂岩,密度分别为2080 kg·m-3、2060 kg·m-3,孔隙度均为30%,裂隙纵横比为10-3,裂缝密度为0.1,时间尺度参数为2×10-3.10 Hz时的水饱和砂岩的纵波速度、横波速度分别为2790 m·s-1、1463 m·s-1.
根据体积模量与剪切模量表达式(Chapman et al., 2006),分别计算出水饱和与气饱和时,每个频率所对应的纵横波速度与纵横波品质因子,如图 4.可见在1~1000 Hz频率范围内,纵横波速度均随频率增加而增加,同时含水砂岩纵波速度明显高于含气砂岩,而横波速度两者基本相等.当log(2πf)大 约等于2.5时,即频率为50 Hz时,频散与衰减最大.说明该模型参数 能够体现地震频带的衰减与频散.
为了测试频变反射系数近似方程(4)的精度,利 用上述模型计算不同频率(10 Hz、20 Hz、30 Hz、40 Hz、 50 Hz、60 Hz、70 Hz和80 Hz)的含气砂岩与含水砂岩的纵横波速度、密度与Gassmann流体项f,如表 1.可见,砂岩(尤其是含气砂岩)的Gassmann流体项f在低频处较低,与中高频带差异较大.考虑到本文频变AVO反演利用频带差异,为避免低频值对反演的过多影响,实际反演中尽量选取靠近主频两侧的频率进行谱分解,同时可以将谱分解结果作谱均衡处理.
根据表 1中的参数值,分别利用Zoeppritz方程与频变近似方程(4)计算不同入射角(0°、5°、10°、15°、20°、25°、30°、35°和40°)下,含气砂岩与泥岩界面反射系数以及含水砂岩与泥岩界面反射系数,如图 5.由图 5a,5b可知,随着入射角增加,精确反射系数先变小后变大,同时频率越高,反射系数越小.比较图 5c,5d和图 5a,5b可以发现,基于频变AVO近似得到的反射系数变化趋势与前者相同,且在20°入射角范围内,反射系数值比较接近真实值,而高频反射系数在0°至40°范围内均与真实值接近.
利用加拿大卡尔加里大学的CREWES Matlab地震工具箱中的反射率法,分别正演得到砂岩含气与砂岩含水时,三层模型的叠前地震记录以及动校正结果,如图 6和图 7.
可见,随着偏移距的增大,上界面的反射振幅很快变小,且出现极性反转.相比砂岩含水情况,含气砂岩由于衰减较强,底界面的反射能量比较微弱.反演的Df如图 8.可见含气砂岩上界面与下界面的频散程度明显大于含水砂岩.
将本文方法应用于某二维实际叠前数据,叠后剖面如图 9a,其中含气层顶底如圆圈所示.蓝线代表纵波速度测井曲线,可见气层处纵波速度明显降 低.过气层的井旁地震道的时频分析结果如图 9b,可以发现气层主频位于40 Hz左右.
由于叠后资料具有较高的信噪比,故频带选择参考叠后资料的时频分析结果.图 10是第75道的振幅谱.主频约为40 Hz,频带范围为5~60 Hz.因 此,在带宽范围内等间距地选取20 Hz,30 Hz,40 Hz,50 Hz,60 Hz 作为谱分解频率较为合理.
同时在谱分解过程中,鉴于地震信号为实信号,故需要考虑选取小波变换结果的实部还是模值作为反演输入.虽然前人的FAVO反演方法(Wilson,2010;吴小羊,2010;Zhang et al., 2011)采用了模值,但本研究发现相比利用模值的反演结果(图 11),选取小波变换结果的实部可以使分频剖面具有较高分辨率,反演结果也更有效.
图 12a和图 13a分别是20 Hz分频剖面和30 Hz 分频剖面.相应的第75道的时频分析(图 12b和图 13b)结果说明所得分频剖面的主频符合预期的频 带.同时可以看到30 Hz剖面可以更清晰地刻画气层.
图 14a为Df的反演结果,图 14b为Zhang等(2011)提出的DVP频散属性反演结果.可见,Df对气层的定位较准确.
选取另一个实际工区的地震数据做测试,150 m 偏移距剖面如图 15a.反演的流体因子如图 15b.柱状图为测井解释结果,红色表示气层,白色表示泥岩.反演结果上的红色异常与解释结果基本对应.
本文基于f-μ-ρ AVO反射系数方程,推导了频变反射系数近似公式,从而可以利用谱分解所得到的不同频带地震资料,估计新的孔隙流体参数属性Df进行流体预测.反演采用贝叶斯先验约束保证结果合理性与稳定性.基于实际岩石物理参数的模型测试表明,该频变AVO近似在入射角不大的情况下具有较高的精度,单道反演结果进一步验证了该属性可以将含气砂岩与含水砂岩区分开来.实际资料试处理说明了基于新近似方程的频变AVO反演对油气识别的可行性.
[1] | Al-Harrasi O H, Kendall J M, Chapman M. 2011. Fracture characterization using frequency-dependent shear wave anisotropy analysis of microseismic data. Geophysical Journal International, 185(2): 1059-1070, doi: 10.1111/j.1365-246X.2011.04997.x. |
[2] | Castagna J, Sun S, Siegfried R. 2003. Instantaneous spectral analysis: Detection of low-frequency shadows associated with hydrocarbons.The Leading Edge, 22(2): 120-127, doi: 10.1190/1.1559038. |
[3] | Chapman M. 2003. Frequency-dependent anisotropy due to meso-scale fractures in the presence of equant porosity. Geophysical Prospecting, 51(5): 369-379, doi: 10.1046/j.1365-2478.2003.00384.x. |
[4] | Chapman M, Liu E R, Li X Y. 2005. The influence of abnormally high reservoir attenuation on the AVO signature. The Leading Edge, 24(11): 1120-1125, doi: 10.1190/1.2135103. |
[5] | Chapman M, Liu E R, Li X Y. 2006. The influence of fluid sensitive dispersion and attenuation on AVO analysis. Geophysical Journal International, 167(1): 89-105, doi: 10.1111/j.1365-246X.2006.02919.x. |
[6] | Chapman M. 2009. Modeling the effect of multiple sets of mesoscale fractures in porous rock on frequency-dependent anisotropy. Geophysics, 74(6): D97-D103, doi: 10.1190/1.3204779. |
[7] | Chen X H, He Z H, Wen X T, et al. 2009. Numeric simulation and detection of low frequency shadow. Oil Geophysical Prospecting (in Chinese), 44(3): 298-303. |
[8] | Chen S Q, Li X Y, Wang S X. 2012. The analysis of frequency-dependent characteristics for fluid detection: a physical model experiment. Applied Geophysics, 9(2): 195-206, doi: 10.1007/s11770-012-0330-8. |
[9] | Cheng B J, Xu T J, Li S G. 2012. Research and application of frequency dependent AVO analysis for gas recognition. Chinese Journal of Geophysics (in Chinese), 55(2): 608-613. |
[10] | Henry F, Cooper Jr. 1967. Reflection and Transmission of Oblique Plane Waves at a Plane Interface between Viscoelastic Media. The Journal of the Acoustical Society of America, 42(5): 1064-1069, 10. 1121/1.1910691. |
[11] | Ebrom D. 2004. The low-frequency gas shadow on seismic sections. The Leading Edge, 23(8): 772, doi: 10.1190/1.1786898. |
[12] | Goloshubin G, Silin D. 2005. Using frequency-dependent seismic attributes in imaging of a fractured reservoir zone. SEG Technical Program Expanded Abstracts, 1417-1420, doi: 10.1190/1.2147954. |
[13] | Innanen K. 2011. Inversion of the seismic AVF/AVA signatures of highly attenuative targets. Geophysics, 76(1): R1-R14, doi: 10.1190/1.3518816. |
[14] | Johnston D H. 2005. The attenuation of seismic waves in dry and saturated rocks . Cambridge: Massachusetts Institute of Technology. |
[15] | Krebes E S. 1984. On the reflection and transmission of viscoelastic waves-Some numerical results. Geophysics, 49(8): 1374-1380, doi: 10.1190/1.1441765. |
[16] | Liu L F, Cao S Y, Wang L. 2011. Poroelastic analysis of frequency-dependent amplitude-versus-offset variations. Geophysics, 76(3): C31-C40, doi: 10.1190/1.3552702. |
[17] | Maultzsch S, Chapman M, Liu E R, et al. 2003. Modelling frequency-dependent seismic anisotropy in fluid-saturated rock with aligned fractures: Implication of fracture size estimation from anisotropic measurements. Geophysical Prospecting, 51(5): 381-392, doi: 10.1046/j.1365-2478.2003.00386.x. |
[18] | Nechtschein S, Hron F. 1997. Effects of anelasticity on reflection and transmission coefficients. Geophysical Prospecting, 45(5): 775-793, doi: 10.1046/j.1365-2478.1997.590288.x. |
[19] | Odebeatu E, Zhang J, Chapman M. 2006. Application of spectral decomposition to detection of dispersion anomalies associated with gas saturation. The Leading Edge, 25(2): 206-210, doi: 10.1190/1.2172314. |
[20] | Ren H T, Goloshubin G, Hilterman F J. 2009. Poroelastic analysis of amplitude-versus-frequency variations. Geophysics, 74(6):N41-N48, doi: 10.1190/1.3207863. |
[21] | Russell B H, Gray D, Hampson D P. 2011. Linearized AVO and poroelasticity. Geophysics, 76(3): C19-C29, doi: 10.1190/1.3555082. |
[22] | Sun L. 2009. Attenuation and velocity dispersion in the exploration seismic frequency band . Toronto: University of Toronto. |
[23] | Taner M T, Koehler F, Sheriff R. 1979. Complex seismic trace analysis. Geophysics, 44(6): 1041-1063, doi: 10.1190/1.1440994. |
[24] | Ursin B, Stovas A. 2002. Reflection and transmission responses of a layered isotropic viscoelastic medium. Geophysics, 67(1): 307-323, doi: 10.1190/1.1451803. |
[25] | Wang S X, Li X Y. 2006. Layer stripping of azimuthal anisotropy from P-wave reflection moveout in orthogonal survey lines. Journal of Geophysics and Engineering, 3(1): 1, doi: 10.1088/1742-2132/3/1/001. |
[26] | Wang S X, Li X Y, Qian Z P, et al. 2007. Physical modelling studies of 3-D P-wave seismic for fracture detection. Geophysical Journal International, 168(2): 745-756, doi: 10.1111/j.1365-246X.2006. 03215.x. |
[27] | Wang S X, Li X Y, Di B. 2010. Reservoir fluid substitution effects on seismic profile interpretation: A physical modeling experiment. Geophysical Research Letters, 37(10): L10306, doi: 10. 1029/2010GL043090. |
[28] | Wilson A. 2010. Theory and methods of frequency-dependent AVO inversion . Edinburgh: The University of Edinburgh. |
[29] | Wu X Y. 2010. Frequency dependent AVO inversion using spectral decomposition techniques (in Chinese). Wuhan: China University of Geosciences. |
[30] | Wu X Y, Chapman M, Wilson A, et al. 2010. Estimating seismic dispersion from pre-stack data using frequency-dependent AVO inversion. SEG Technical Program Expanded Abstracts, 84(29): 425-429, doi: 10.1190/1.3513759. |
[31] | Xu D, Wang Y H, Gan Q G, et al. 2011. Frequency-dependent seismic reflection coefficient for discriminating gas reservoirs. Journal of Geophysics and Engineering, 8(4): 508, doi: 10.1088/1742-2132/8/4/003. |
[32] | Yang P J. 2008. Seismic wavelet blind extraction and non-linear inversion (in Chinese). Qingdao: China University of Petroleum (Huadong). |
[33] | Yin X Y, Zhang S X, Zhang F C, et al. 2010. Utilizing Russell Approximation-based elastic wave impedance inversion to conduct reservoir description and fluid identification. Oil Geophysical Prospecting (in Chinese), 45(3): 373-380. |
[34] | Yin X Y, Li C, Zhang S X. 2013a. Seismic fluid discrimination based on two-phase media theory. Journal of China University of Petroleum (in Chinese), 37(5): 38-43, doi: 10.3969/j.issn.1673-5005.2013. 05.006. |
[35] | Yin X Y, Zong Z Y, Wu G C. 2013b. Seismic wave scattering inversion for fluid factor of heterogeneous media. Science China: Earth Sciences (in Chinese), 43(12): 1934-1942, doi: 10. 1007/s11430-013-4783-2. |
[36] | Zhang S X, Yin X Y, Zhang G Z. 2011. Dispersion-dependent attribute and application in hydrocarbon detection. Journal of Geophysics and Engineering, 8(4): 498-507, doi: 1742-2140/8/4/002. |
[37] | Zong Z Y, Yin X Y, Wu G C. 2013. Direct inversion for a fluid factor and its application in heterogeneous reservoirs. Geophysical Prospecting, 61(5): 998-1005, doi: 10.1111/1365-2478.12038. |
[38] | 陈学华, 贺振华, 文晓涛等. 2009. 低频阴影的数值模拟与检测. 石油地球物理勘探, 44(3): 298-303. |
[39] | 程冰洁, 徐天吉, 李曙光. 2012. 频变AVO含气性识别技术研究与应用. 地球物理学报, 2012, 55(2): 608-613. |
[40] | 吴小羊. 2010. 基于频谱分析技术的频散AVO反演研究[博士论文]. 武汉: 中国地质大学. |
[41] | 杨培杰. 2008. 地震子波盲提取与非线性反演[博士论文]. 青岛: 中国石油大学(华东). |
[42] | 印兴耀, 张世鑫, 张繁昌等. 2010. 利用基于Russell近似的弹性波阻抗反演进行储层描述和流体识别. 石油地球物理勘探, 45(3): 373-380. |
[43] | 印兴耀, 李超, 张世鑫. 2013a. 基于双相介质的地震流体识别. 中国石油大学学报(自然科学版), 37(5): 38-43, doi: 10.3969/j.issn.1673-5005.2013.05.006. |
[44] | 印兴耀, 宗兆云, 吴国忱. 2013b. 非均匀介质孔隙流体参数地震散射波反演. 中国科学——地球科学, 43(12): 1934-1942, doi: 10.1007/s11430-013-4783-2. |