2. 中国石油大学(北京) 海洋石油勘探国家工程实验室, 北京 102249;
3. 中国石油大学(华东) 地球科学与技术学院, 青岛 266555
2. National Engineering Laboratory for Offshore oil Exploration, China University of Petroleum, Beijing 102249, China;
3. School of Geoscience, China University of Petroleum, Qingdao 266555, China
地下裂缝系统是碳酸盐岩储层和非常规储层(页岩气、致密气和致密油储层)研究的重要部分.裂缝系统不仅可以连通孤立的孔隙,增加储层的有效孔隙度,而且可以为油气运移提供通道,提高储层的渗透率.一般认为页岩、周期性薄互层和定向排列垂直裂隙的地层是横向各向同性(TI)介质.当地层由于构造运动发生褶皱和上冲作用时,TI介质的对称轴将不再是垂直或水平的,而是与观测坐标系存在夹角,这就形成了倾斜横向各向同性(TTI)介质.TTI介质能更真实地反映地层的地质特征,是众多地球物理学家研究的重点内容.
裂缝岩石物理研究是联系裂缝特征参数与地震响应的桥梁.目前常用的两个模型为薄硬币形状裂隙模型(Hudson,1980)和线性滑移模型(Schoenberg,1983).两者的不同在于前者面向研究微裂隙对岩石的影响,后者主要用于描述裂缝对岩石的影响.关于两者的等价关系,Schoenberg和Douma(1988)指出当缝隙面比较平缓、且缝隙包含物模量较小时,分别利用薄硬币形状裂隙模型和线性滑移模型计算的介质弹性模量相等.由于薄硬币形状裂隙模型考虑了缝隙形状、缝隙密度及缝隙充填物的影响,因此鉴于两个模型的等价关系,可以研究线性滑移模型中裂缝岩石物理参数与裂缝形状、裂缝密度及裂缝充填物的关系.Mavko等(2009)分析了薄硬币形状模型的假设条件应用限制,指出在地震频带范围内应先计算含裂缝干岩石骨架的弹性模量,再利用流体替换方法计算饱和岩石的弹性模量.研究表明,定向排列的裂缝可以引起地震波传播的各向异性特征,因此在裂缝型干岩石的流体替换过程中不能忽略岩石各向异性的特征.同时,Mavko等(2009)在岩石物理手册中展示了Gassmann(1951)各向异性岩石的流体替换方程,该方程可以用来实现裂缝型饱和岩石弹性模量、裂缝岩石物理参数的求解.本文从岩石物理等效理论出发,假设各向异性岩石的流体为气体,选取线性滑移模型实现裂缝型岩石物理等效模型的构建.
AVAZ反演方法是预测裂缝的一种有效手段.利用地震反演得到与裂缝特征相关的岩石物理参数,可以指导地下裂缝预测,因此探索裂缝岩石物理参数与地震反射之间的关系式,进而开展裂缝岩石物理参数叠前地震反演研究显得尤为重要.Rüger(1996, 1997)推导了HTI介质的线性近似Zoeppritz方程.基于该方程,Mallick等(1998)利用方位地震道集估计了弹性参数和各向异性参数.Chen等(2014b)在假设小角度入射的情况下进一步简化了该方程,推导了含裂缝柔度的新方程,并利用最小二乘算法直接反演裂缝柔度和相对反射系数.Ivanov和Stovas (2016)推导了弱各向异性近似下两层TTI介质间的二维P波反射系数方程,该方程包含VTI介质的Thomsen三参数.Wang和Peng(2016)推导了两层粘弹TTI介质间的精确六阶Zoeppritz方程.Wang等(2016)推导了两层弱各向异性TTI介质间的三维P波反射系数方程,该方程可以看成shuey三项式和扰动项的叠加.研究表明,裂缝型储层各向异性弹性阻抗表现出随入射角和方位角同时变化的特征,也将其称为方位各向异性弹性阻抗.国内外大量学者研究了弹性阻抗反演以及基于弹性阻抗的参数估测方法(王保丽等,2005;印兴耀等,2013;陈怀震等, 2013, 2014a).
AVO反演通常是病态的,需要对反问题进行正则化约束才能获得稳定的解,贝叶斯理论是建立AVO反演方程的理想框架,它通过联合观测数据与先验信息进行约束,有效提高反演过程的稳定性和反演结果的可靠性.Buland和Omre(2003)基于贝叶斯理论框架,提出了贝叶斯线性AVO反演方法.同时在贝叶斯公式中引入三变量的高斯分布作为先验分布,有效地提高了反演的稳定性.Buland等(2003)基于贝叶斯线性反演框架,将空间相关函数引入到先验分布函数中,用于描述弹性参数的空间相关性,进而提高了多道反演的稳定性和精度.陈建江和印兴耀(2007)基于贝叶斯理论,采用了一种以先验地质信息为约束条件的非线性AVO反演方法,提高了参数估算的分辨率和抗噪声能力, 并且可以获得符合地质规律的反演结果.杨培杰和印兴耀(2008)提出了一种基于非线性二次规划的叠前三参数反演方法,反演方法运算速度快捷,抗噪能力强.王守东和王波(2012)针对时移地震差异数据,提出了一种基于贝叶斯理论的AVO波形反演方法,同时反演出纵波阻抗、横波阻抗和密度的变化.Liu等(2016)基于贝叶斯理论,提出了一种反射率AVO反演方法,相比Zoeppritz方程及其近似式,该方法更能适应透射损失和层间多次波.
本文首先从TTI介质的岩石物理等效模型构建出发,利用测井数据估测可靠的纵横波速度和裂缝岩石物理参数,分析裂缝岩石物理参数的变化特征,直接从测井曲线上识别出裂缝储层位置,指导后续裂缝储层的地震反演预测.其次,通过推导含裂缝岩石物理参数的线性TTI反射系数方程,探索利用不同方位的角度叠加地震道集直接反演裂缝岩石物理参数的方法.最后,本文以包含裂缝柔度的线性TTI反射系数方程为正演模型,基于贝叶斯理论建立反演框架,采用高斯分布约束反演过程,形成AVAZ四个参数的反演算法,实现裂缝岩石物理参数的可靠估测.该流程不仅依靠岩石物理模型从测井上直接估测岩石物理参数,识别裂缝储层位置,而且充分利用不同方位的部分角度叠加地震道集反演裂缝柔度,进而提取整个工区的裂缝岩石物理参数,为储层裂缝预测和描述提供可靠的支撑.
1 方法原理 1.1 裂缝型岩石物理等效模型研究表明,地下裂缝的定向排列将引起地震响应的各向异性特征.当裂缝尺度小于地震波长,等效介质理论可以用来构建裂缝型岩石物理模型.其中,选用线性滑移模型来研究裂缝对各向同性背景岩石的影响,并直接估测裂缝岩石物理参数(ΔN和ΔT),避免了由各向异性参数向裂缝岩石物理参数转换过程中引入的误差.
已知线性滑移模型(Schoenberg和Sayers, 1995)描述TI介质的裂缝型岩石刚度矩阵为
(1) |
其中,M=λ+2μ.λ和μ为是各向同性背景介质的拉梅参数,ΔN为法向柔度,表示垂直裂缝面上裂缝对地震波的影响,ΔT为切向柔度,表示平行于裂缝面上裂缝对地震波的影响.
Bakulin等(2000)根据薄硬币裂隙模型和线性滑移模型的关系,对裂缝岩石刚度矩阵中的ΔN和ΔT进行求解.
(2) |
其中,g=μ/(λ+2μ).e为裂缝密度, α为裂缝高宽比. k'和μ'分别为缝隙充填物的体积模量和剪切模量.研究表明,在地震频带范围内构建裂缝型岩石物理等效模型时,应先利用线性滑移模型添加干裂缝(k'=0,μ'=0)形成裂缝型干岩石骨架,再进行流体替换形成饱和裂缝型岩石.
联合线性滑动模型和Thomsen各向异性参数表达式(Thomsen,1986),TI介质的各向异性参数与裂缝柔度之间的关系可表示为(Bakulin等, 2000)
(3) |
本文假设k'=0和μ'=0,利用线性滑动模型构建裂缝型含气页岩,因此流体替换问题不在本文研究范围内.
1.2 裂缝柔度敏感性分析研究裂缝柔度和裂缝密度、流体类型的关系可以揭示裂缝柔度估计的重要性,为后续的地震反演奠定基础.按照公式(2),本文分析了裂缝柔度随裂缝密度、流体类型的变化特征,如图 1所示.试验中裂缝高宽比为0.0005.图 1为当裂缝饱和流体分别为水、油和气时,裂缝柔度与纵横波速度比(β/α)之间的关系.图 1表明ΔN和ΔT随裂缝密度的增加而增加.当流体类型为气体时,ΔN变化最为明显,其次为油,最后为水.然而对气、油和水填充的不同裂缝,ΔT变化很小.因此,我们可以利用ΔT直接预测岩石的裂缝密度;结合预测的裂缝密度,利用ΔN可估计裂缝填充物的类型.
Wang等(2016)推导了两层弱各向异性TTI介质间的三维P波反射系数方程
(4) |
其中,θ为入射角,φ为方位角.RPPiso(θ)的表达式为
(5) |
Aiso, Biso and Ciso分别为与拉梅参数有关的截距项,梯度项和曲率项
(6) |
其中,Vp和Vs分别为沿对称轴方向的纵横波速度,ρ为密度,Z=ρVp,Δ表示界面两侧参数的变化,横线表明界面两侧参数的均值(在扰动方法中可以作为背景介质参数值).RPPani(θ,φ)同样可以划分为截距项,梯度项,曲率项和扰动项,其表达式为
(7) |
其中,
(8) |
(9) |
(10) |
(11) |
其中,θ01为上层介质的对称轴倾角,(ε1,δ1,γ1)为上层介质的各向异性参数.θ02和(ε2,δ2,γ2)为下层介质的相应参数.
从Wang的TTI介质反射系数公式出发,舍掉密度反射系数项,用裂缝柔度参数替换Thomsen各向异性参数,假设上层为各向同性半空间,下层为TTI半空间,就可得到含裂缝柔度参数的发育定向排列的倾斜裂缝岩石的反射系数近似公式:
(12) |
其中
公式(12)表明利用地震反演方法可以估测裂缝柔度.如果偏移距采样数为m,方位角采样数为n,则模拟的AVAZ角道集可表示为如下矩阵相乘的形式:
(13) |
其中,
(14) |
它们只与背景速度场及角度有关.综上所述,基于线性近似公式的角道集模拟过程可以简写为以下矩阵形式:
(15) |
其中,模型参数为
(16) |
子波矩阵为
(17) |
L=WF为正演算子,n为观测地震数据的噪声.
按照叠前地震反演的贝叶斯理论,已知叠前地震数据,模型参数的反演问题可以通过以下公式求解:
(18) |
其中,P(m|d)是后验概率分布;P(d|m)是观测数据的似然函数;P(m)为先验分布;P(d)是一个常数,称为边缘概率密度,它的定义规则是使后验概率分布的总和等于1.
假设噪声服从高斯分布(Todoeschuck et al., 1990),则观测数据的似然函数
(19) |
其中,Cd表示数据的协方差矩阵.假设地震数据中的噪声相互独立,即Cd=σn2I,其中σn2是噪声协方差,I是Nd×Nd的单位矩阵,Nd是数据的大小.
在贝叶斯反演内,选取不同的先验分布将得到不同性质的解,例如选取高斯分布将得到稳定光滑的解,而选取柯西分布会得到稀疏的解.假设先验分布具有以下的统一形式:
(20) |
其中,Cm为模型参数的协方差矩阵.
从模型参数的后验概率分布获得模型参数预测值的一种方法是直接求最大后验概率解,即求解以下目标函数的极小值,
(21) |
其中,
(22) |
求解目标函数的极小值等价于令其导数等于零,整理得到
(23) |
其中,
为了验证本文含裂缝岩石物理参数的TTI介质近似PP反射系数方程的正确性和有效性,设计一个两层TTI介质模型.对于模型1和2,上层介质为相同的TTI介质,其参数为VP01 = 2.8 km·s-1,VS01= 1.4 km·s-1,ρ=1.8 g·cm-3,δ= 0.05,ε=0.02,γ= 0.01;下层介质分别为强、弱阻抗差和各向异性条件下的TTI介质,参数如表 1所示.分别采用方程(A8)、Wang等(2016)推导的反射系数方程和精确方程(Rokhlin等, 1985)对不同模型分别计算近似和精确PP波反射系数,探讨本文推导的近似方程的准确度.注意,利用方程(A8)计算AVO曲线之前,需要利用公式(3)将Thomsen各向异性参数转化为裂缝柔度ΔN和ΔT.其中,上下两层介质的对称轴倾角都为45°,入射角为0°~50°,方位角为0°, 45°和90°.
图 2为本文近似曲线、Wang等(2016)近似曲线和精确曲线的对比.图中,实线是精确值,点实线是Wang等(2016)近似值,虚线是本文近似值;蓝色、绿色和红色分别表示方位角0°,45°和90°的反射系数曲线.图 2表明:本文近似值的趋势与Wang等(2016)近似值、精确值相同,但存在一定差异.随着入射角度的增大,这种差异逐渐增大,这是由于本文近似方程推导中使用了小角度入射假设条件.注意,在小角度入射情况下,本文近似曲线相比Wang等(2016)近似曲线更加接近精确曲线.此外,图 2a和2b分别显示在强、弱阻抗差和各向异性条件下,三种曲线的变化趋势相同,然而反射系数振幅大小差异较大.最后,图 2显示当方位角为90°时,相比其他方位角,新近似方程的反射系数与精确曲线的差异更大.因此在勘探中,方位效应不能忽略.总体来说,近似解与精确解吻合得比较好,这说明本文的近似公式具有较高的精度,有助于理解各向异性对AVO响应的影响,并且可以为后续方位各向异性的裂缝岩石物理参数反演提供理论基础.
利用模拟测井解释数据合成理论方位各向异性地震数据,为后续的AVAZ裂缝参数反演提供输入. 图 3为模型测井解释数据,分别是纵波速度,横波速度,密度,裂缝密度,裂缝高宽比.图 4为从模型数据上利用岩石物理等效模型估测的纵横波相对反射系数和裂缝岩石物理参数.
对比图 3和图 4可以看出,裂缝岩石物理参数与测井解释信息的变化趋势具有一定的相关性,而且通过分析裂缝岩石物理参数的数值变化可以从测井数据上识别裂缝发育位置.在反演试算前,通过倾角测井方法或者速度分析等方法可以计算对称轴倾角为45°.
为了检验研究反演算法的可行性和抗噪性,利用图 4所示的模型参数进行反演算法试算.利用公式(13)分别计算出不同时间采样点、不同入射角度(以5°为间隔,0°到50°之间的10个角度)、不同方位角度(以30°为间隔,30°到150°之间的5个角度)的纵波反射系数,并与30 Hz的雷克子波进行褶积运算合成方位叠前角道集,如图 5a所示.为了检验反演算法的抗噪性,在合成叠前角道集上加入信噪比为4和1的随机噪声,如图 5b和c所示.利用本文算法进行反演,检验方法的抗噪性.从图 6中的反演结果可以看出,随机噪声对纵横波相对反射系数的估计影响较小,而对裂缝柔度参数的反演结果影响较大,尤其是在模型参数极大(或极小)或者反射较弱的地方,影响尤为严重,但整体还是吻合很好.总体上,可以看出该反演算法具有较高的抗噪能力.
本文主要研究了基于贝叶斯线性AVAZ的TTI介质裂缝柔度参数反演方法,通过构建裂缝岩石物理模型,利用模拟测井解释数据实现裂缝储层的纵横波相对反射系数和裂缝柔度参数估测;通过推导含裂缝柔度参数的TTI介质近似公式,构建了贝叶斯框架下的TTI介质裂缝柔度参数AVAZ反演方法;利用模型测井数据合成的不同方位叠加地震道集反演裂缝柔度参数实现了对裂缝柔度参数的直接反演,减少了间接求解带来的误差,探讨了反演方法的抗噪性.从基于两层TTI介质模型的Wang等(2016)近似方程、精确方程(Rokhlin等, 1985)和新方程的反射系数曲线对比结果可以看出,在小角度入射的情况下,推导的含裂缝岩石物理参数的TTI介质方位各向异性公式与精确反射系数曲线吻合较好,且具有明确的物理意义;从基于模型数据的反演结果中可以看出,纵横波相对反射系数受随机噪声影响较小,而反演的裂缝岩石物理参数受噪声影响较大,但与真实数据吻合较好,验证了反演方法可靠性和抗噪能力.后续研究工作如下:(1)由于本文方法的反演参数不包括TTI介质对称轴倾角,后续拟进行TTI介质对称轴倾角的反演方法研究.(2)将本文方法应用于实际数据测试,进一步明确方法的反演效果.
附录 含裂缝岩石物理参数的TTI介质反射系数方程推导式(4)中各向同性部分可变为
(A1) |
如果
(A2) |
小角度入射情况下,舍掉密度反射系数项,则式(4)的各向同性部分可变为
(A3) |
对式(4)中各向异性部分的下层TTI介质反射系数进行化简,则有
(A4) |
将弱各向异性Thomsen参数与裂缝岩石物理参数的关系式(式3)带入上式可得:
(A5) |
则将含裂缝岩石物理参数的项集中一起,通过进一步化简最终可得
(A6) |
同理,式(4)中各向异性部分的上层TTI介质反射系数可推导为
(A7) |
那么两层弱各向异性TTI介质间的P波反射系数方程可表示为
(A8) |
假设上层为各向同性半空间,下层为TTI半空间,上式即变为公式(12).
Bakulin A, Grechka V, Tsvankin I. 2000. Estimation of fracture parameters from reflection seismic data-part Ⅰ:HTI model due to a single fracture set. Geophysics, 65(6): 1788-1802. DOI:10.1190/1.1444863 |
Buland A, Kolbjórnsen O, Omre H. 2003. Rapid spatially coupled AVO inversion in the Fourier domain. Geophysics, 68(3): 824-836. DOI:10.1190/1.1581035 |
Buland A, Omre H. 2003. Bayesian linearized AVO inversion. Geophysics, 68(1): 185-198. DOI:10.1190/1.1543206 |
Chen H Z, Yin X Y, Gao C G, et al. 2014a. AVAZ inversion for fluid factor based on fracture anisotropic rock physics theory. Chinese Journal of Geophysics (in Chinese), 57(3): 968-978. DOI:10.6038/cjg20140326 |
Chen H Z, Yin X Y, Qu S L, et al. 2014b. AVAZ inversion for fracture weakness parameters based on the rock physics model. Journal of Geophysics and Engineering, 11(6): 065007. DOI:10.1088/1742-2132/11/6/065007 |
Chen J J, Yin X Y. 2007. Three-parameter AVO waveform inversion based on Bayesian theorem. Chinese Journal of Geophysics (in Chinese), 50(4): 1251-1260. |
Gassmann F. 1951. Uber die elastizitat poroser medien. Veirteljahrsschrift der Naturforschenden Gesellschaft in Zzirich, 96: 1-23. |
Hudson J A. 1980. Overall properties of a cracked solid. Mathematical Proceedings of the Cambridge Philosophical Society, 88(2): 371-384. DOI:10.1017/S0305004100057674 |
Ivanov Y, Stovas A. 2016. Weak-anisotropy approximation for P-wave reflection coefficient at the boundary between two tilted transversely isotropic media. Geophysical Prospecting, 65(2): 485-502. |
Liu H X, Li J Y, Chen X H, et al. 2016. Amplitude variation with offset inversion using the reflectivity method. Geophysics, 81(4): R185-R195. DOI:10.1190/geo2015-0332.1 |
Mallick S, Craft K L, Meister L J, et al. 1998. Determination of the principal directions of azimuthal anisotropy from P-wave seismic data. Geophysics, 63(2): 692-706. DOI:10.1190/1.1444369 |
Mavko G, Mukerji T, Dvorkin J. 2009. The Rock Physics Handbook:Tools for Seismic Analysis of Porous Media. 2nd ed. Cambridge: Cambridge University Press.
|
Rüger A. 1996. Reflection coefficients and azimuthal AVO analysis in anisotropic media. Golden: Colorado School of Mines.
|
Rüger A. 1997. P-wave reflection coefficients for transversely isotropic models with vertical and horizontal axis of symmetry. Geophysics, 62(3): 713-722. DOI:10.1190/1.1444181 |
Rokhlin S I, Bolland T K, Adler L. 1985. Reflection and refraction of elastic waves on a plane interface between two generally anisotropic media. The Journal of the Acoustical Society of America, 79(4): 906-918. |
Schoenberg M. 1983. Reflection of elastic waves from periodically stratified media with interfacial slip. Geophysical Prospecting, 31(2): 265-292. DOI:10.1111/gpr.1983.31.issue-2 |
Schoenberg M, Douma J. 1988. Elastic wave propagation in media with parallel fractures and aligned cracks. Geophysical Prospecting, 36(6): 571-590. DOI:10.1111/gpr.1988.36.issue-6 |
Schoenberg M, Sayers C M. 1995. Seismic anisotropy of fractured rock. Geophysics, 60(1): 204-211. DOI:10.1190/1.1443748 |
Thomsen L. 1986. Weak elastic anisotropy. Geophysics, 51(10): 1954-1966. DOI:10.1190/1.1442051 |
Todoeschuck J P, Jensen O G, Labonte S. 1990. Gaussian scaling noise model of seismic reflection sequences:Evidence from well logs. Geophysics, 55(4): 480-484. DOI:10.1190/1.1442857 |
Wang B L, Yin X Y, Zhang F C. 2005. Elastic impedance inversion and its application. Progress in Geophysics (in Chinese), 20(1): 89-92. |
Wang H W, Peng S P. 2016. Reflection coefficient of qP, qS and SH at a plane boundary between viscoelastic TTI media. Geophysical Journal International, 204(1): 555-568. DOI:10.1093/gji/ggv475 |
Wang H W, Peng S P, Du W F. 2016. Azimuthal AVO of P-wave at the boundary between two TTI media. //78th EAGE Conference and Exhibition 2016. Vienna: EAGE.
|
Wang S D, Wang B. 2012. Time-lapse Bayesian AVO waveform inversion. Chinese Journal of Geophysics (in Chinese), 55(7): 2422-2431. DOI:10.6038/j.issn.0001-5733.2012.07.026 |
Yang P J, Yin X Y. 2008. Non-linear quadratic programming Bayesian prestack inversion. Chinese Journal of Geophysics (in Chinese), 51(6): 1876-1882. |
Yin X Y, Zhang S X, Zhang F. 2013. Two-term elastic impedance inversion and Russell fluid factor direct estimation method for deep reservoir fluid identification. Chinese Journal of Geophysics (in Chinese), 56(7): 2378-2390. DOI:10.6038/cjg20130724 |
陈怀震, 印兴耀, 杜炳毅, 等. 2013. 裂缝型碳酸盐岩储层方位各向异性弹性阻抗反演. 地球物理学进展, 28(6): 3073-3079. DOI:10.6038/pg20130631 |
陈怀震, 印兴耀, 高成国, 等. 2014. 基于各向异性岩石物理的缝隙流体因子AVAZ反演. 地球物理学报, 57(3): 968-978. DOI:10.6038/cjg20140326 |
陈建江, 印兴耀. 2007. 基于贝叶斯理论的AVO三参数波形反演. 地球物理学报, 50(4): 1251-1260. |
王保丽, 印兴耀, 张繁昌. 2005. 弹性阻抗反演及应用研究. 地球物理学进展, 20(1): 89-92. |
王守东, 王波. 2012. 时移地震资料贝叶斯AVO波形反演. 地球物理学报, 55(7): 2422-2431. DOI:10.6038/j.issn.0001-5733.2012.07.026 |
杨培杰, 印兴耀. 2008. 非线性二次规划贝叶斯叠前反演. 地球物理学报, 51(6): 1876-1882. |
印兴耀, 张世鑫, 张峰. 2013. 针对深层流体识别的两项弹性阻抗反演与Russell流体因子直接估算方法研究. 地球物理学报, 56(7): 2378-2390. DOI:10.6038/cjg20130724 |