基于黏弹性Chapman-Kelvin模型的裂缝储层频变AVAZ反演研究

2024-02-04张青廖建平刘和秀周林张学娟

地球物理学报 2024年2期

张青,廖建平,*,刘和秀,周林,张学娟

1 湖南科技大学地球科学与空间信息工程学院,湘潭 411201 2 中国石化地球物理重点实验室,中石化石油物探技术研究院有限公司,南京 211103 3 重庆科技学院非常规油气开发研究院,重庆 401331

0 引言

随着油气勘探逐渐从构造油气藏转向非常规油气藏,地下裂缝在油气生产中发挥着重要作用(Zeng et al.,2013;Ghasemi et al.,2020;Lei et al.,2022).Sullivan等(2006)表明全球三分之一以上的油气储存在碳酸盐岩储层中,并且大部分为裂缝储层.对于致密储层,裂缝不仅可以连通孤立孔隙,增加储层有效孔隙度,提高渗透率,更是储层中主要的渗流通道(Lei and Wang,2016;Shi et al.,2018;Hou et al.,2021).因此,裂缝描述是包括碳氢化合物和地下水的生产和监测二氧化碳的地质储存等各种应用的一个关键问题(Iding and Ringrose,2010;刘敬寿等,2019;Sheng et al.,2021;Yasin et al.,2022).Liu和Martinez(2012)表明,可以通过露头、岩心或钻孔测井观测和地震资料解释等方法来描述裂缝的特征.通常需要将所有这些方法与不同尺度的测量数据结合起来,对裂缝储层进行准确的描述.

裂缝储层最基本的特征是各向异性(Cao et al.,2018;Song et al.,2019,2020;Jiang et al.,2021).地震各向异性是描述裂缝特征最常用的方法,并且为了建立各向异性参数与裂缝系统特征之间的联系,使用裂缝储层的等效介质理论来描述裂缝岩石(Pan et al.,2018;Ding et al.,2020;Wang et al.,2022).对裂缝的假设不同,目前常用的等效介质理论方法有两种.其一是基于夹杂物的Hudson模型(Hudson,1981),将裂缝视为嵌入各向同性固体中的椭球体孔隙;另一种是线性滑动模型(Schoenberg,1980),认为裂缝具有线性滑动边界条件的薄弱面.Hudson模型和Schoenberg(1980)提出的线性滑动模型是常用的静态等效介质模型.当储层中存在一组垂直排列的裂缝时,该储层可以看作具有水平对称轴的横向各向同性介质,基于这种静态等效介质模型的裂缝预测方法无法获取裂缝长度、裂缝密度等裂缝参数信息,也不能预测出储层的孔隙度、含水饱和度等物性信息.同样,这类模型也不能解释在实际裂缝储层岩石中地震波的速度频散和衰减现象,这将与真实裂缝储层岩石中的频变各向异性理论相违背(李博南等,2017;Kumar et al.,2018;潘新朋,2019).流体在裂缝与孔隙之间流动的动态等效介质理论模型,更加全面地描述出实际裂缝储层中的频变各向异性特征(Kobayashi and Mavko,2016;Gao and Deng,2022).这类模型通过地震波致流作用,可以反映出真实裂缝储层岩石中的地震波频散和衰减现象.Thomsen(1995)通过考虑裂缝和孔隙之间的流体传递,认识到频散对解释地震各向异性的重要性;Hudson等(1996)建立了一个动态等效介质模型,该模型考虑了整个频率范围内地震波传播特性,但在低频极限时与Thomsen(1995)等径孔隙模型的相应结论不一致,所以不能在全频段内有效;Chapman(2003)基于Chapman等(2002)的颗粒喷射流的孔弹性理论,结合一组定向排列的介观尺度硬币型裂缝,考虑二者之间的相互作用,提出了各向异性模型,通过对比各向异性Gassmann方程和Hudson模型在低频和高频极限下的结论,发现与其一致,说明Chapman(2003)动态等效介质模型可以反映出实际裂缝储层岩石中的地震波的频散和衰减现象.

近年来,一些学者开展了基于动态等效介质模型的裂缝储层参数反演的研究(刘倩等,2016;潘新朋等,2018;Li et al.,2017;Zhang et al.,2022).对于裂缝储层的反演,主要有定性反演和定量反演两种类型.其中,定量反演是指基于地震各向异性来反演裂缝密度和裂缝走向等裂缝信息,包括叠前AVAZ反演等方法(陈怀震等,2014;刘宇巍等,2016;葛子建等,2018;陈珂磷等,2022).Maultzsch等(2003)根据Chapman(2003)模型,反演了裂缝密度和裂缝半径,通过反演得到的裂缝参数信息与裂缝钻孔成像的观测结果一致,说明Chapman(2003)模型的有效性;Ali和Jakobsen(2011)通过T矩阵方法求取了裂缝孔隙介质的等效刚度矩阵系数,并且通过求解Christoffel方程,建立裂缝密度和裂缝走向与地震波速度频散和衰减之间的关系,同时反演了裂缝密度和裂缝走向信息;基于Chapman(2003)模型,Al-Harrasi等(2011)使用网格搜索法对裂缝走向、裂缝密度和裂缝长度进行反演,反演得到的裂缝尺寸与地层露头观测结果一致;Ali和Jakobsen(2014)结合T矩阵法和HTI介质PP波反射系数的Rüger近似方程(Rüger,1998),计算了含有一组垂直排列的介观尺度裂缝的裂缝孔隙介质反射PP波频变AVAZ数据,通过T矩阵法所计算的等效刚度系数建立了裂缝参数与频变AVAZ数据的关系,结合贝叶斯方法和蒙特卡罗随机采样方法对裂缝密度、裂缝走向、裂缝开度和裂缝半径进行了同时反演,并表明只要在反演过程中采用动态等效介质模型,并提供有效的裂缝长度先验信息,可以改善裂缝储层的特征;Liu等(2021)将频变各向异性用于裂缝储层的检测和评价,并研究了裂缝储层随方位变化的规律.由于可以通过分析纵波反射率随方位角和偏移距的变化来确定裂缝方向和裂缝密度,引起了人们的广泛关注(Rüger,1998;Hall and Kendall,2003;刘喜武等,2015;Xie et al.,2020;Ma et al.,2022).然而,由于传统的地震各向异性对裂缝长度不敏感(Hudson,1981;Liu等,2000),因此,从纵波反射中估计裂缝大小仍然很具有挑战性.

本文研究目的是基于已提出的黏弹性Chapman-Kelvin的动态等效介质模型,深入研究在黏弹性裂缝储层中地震波的传播特性,分析裂缝储层参数(主要为裂缝密度、裂缝长度、孔隙度和含水饱和度)对地震波传播特性的影响,以便为复杂裂缝的储层参数反演问题研究提供有力的理论依据.我们考虑到传统的等效介质模型未充分考虑非完全弹性介质理论和频变各向异性理论的双相或多相流体假设,也无法对地震波在实际裂缝储层中的频散和衰减给出准确合理的解释.基于黏弹性Chapman-Kelvin的动态等效介质模型(廖建平等,2023a;Liao et al.,2023a,b),我们将关于黏弹性、喷射流以及斑块效应耦合的频变弹性系数结合到Schoenberg和Protázio(1992)概括的Zoeppritz方程,计算出反射系数,同时考虑到发生地震频散时,反射系数和频率产生关系,构建了在角度、方位和时间域内的新型正演方程.最后,基于PP波频变反射系数对裂缝储层参数变化有较好的敏感性,开展了两种反演方法研究:其一是基于贝叶斯理论直接反演方法;其二是基于MCMC的随机反演方法.通过频变AVAZ进行储层参数的反演研究,结果表明,MCMC随机反演方法在缺失先验的储层参数信息时,反演结果的不确定性较强.当存在有效且足够的先验信息时,反演结果的可靠性进一步提升.而基于贝叶斯理论直接反演方法,由于在方位角AVO分析中对地震频散的考虑,则显示出区别大尺度裂缝和微尺度裂隙的潜力,进而说明基于频变反射系数进行储层参数反演的可行性.

1 裂缝储层的等效介质建模

1.1 基于黏弹性Chapman-Kelvin动态等效介质模型

为了描述实际的裂缝储层介质在全频带范围内的地震响应规律,Chapman等(2002)提出一种微观结构的喷射理论模型,该理论模型在低频时符合Gassmann方程,在高频时符合与喷射相关的频散,并且描述了在含有微裂隙和孔隙的岩石中,流体在微裂隙和孔隙之间,以及微裂隙之间的流动.Chapman(2003)认为,频变各向异性理论的有效性不仅取决于裂缝尺度,还取决于颗粒尺度的流体效应.所以,在基于Chapman等(2002)的喷射理论基础上,描述了流体流动引起的喷射流产生的各向异性频散效应,引入了一套垂直排列的介观尺度裂缝.该模型考虑了介观和微观两种尺度的流体流动,能够对地震频段地震波的频散和衰减给出合理的解释.为了准确描述模型,需要考虑一系列参数,包括孔隙度、裂缝密度、流体参数、裂缝参数以及时间尺度等.Chapman(2003)模型中存在椭球状微裂隙、球形孔隙和定向排列的椭球状介观裂缝.图1为Chapman(2003)模型的示意图,其中微裂隙、孔隙与岩石颗粒均为微观尺度,而裂缝为介观尺度.

图1 Chapman裂缝孔隙介质示意

考虑流体在整个孔隙空间内交换,Chapman(2003)模型的频变刚度矩阵形式为:

(1)

(2)

式(2)的矩阵形式,如式(3)所示:

(3)

其中,C′11、C′33、C′44、C′23、C′13、C′66为两种非混相流体饱和裂缝岩石的有关喷射流和斑块效应的频变各向异性弹性刚度系数,且黏弹性Chapman-Kelvin模型的等效刚度张量在本文的研究中主要取决于以下参数:频率、裂缝密度、裂缝长度、孔隙度和含水饱和度.

1.2 等效黏弹性刚度张量

五个独立的等效黏弹性刚度系数C11、C13、C33、C44、C66是储层参数的函数,储层参数是由基于黏弹性Chapman-Kelvin动态等效介质模型给出.我们对基于Chapman-Kelvin模型的等效黏弹性刚度参数的响应特征进行分析,这里的刚度参数分别是频率f、裂缝密度d、裂缝长度l、孔隙度por和含水饱和度sw的函数.图2为各刚度参数的实部曲线和虚部曲线,其中设定裂缝密度为0.05和孔隙度为0.1,可以观察到刚度张量C11、C13、C33以及C66的实部分量随频率的增大,先增大,当达到高频104Hz以后,逐渐趋于平缓,达到最大值(图2a、c、e、i),而C44的实部分量随频率逐渐减小,当达到高频104Hz以后,逐渐趋于平缓,达到最小值(图2g);刚度张量C11、C13的虚部分量表现为先随频率的增大而增大,但在100Hz以后出现了显著的下降现象,如图2b、d所示.观察到地震频段为0到100 Hz的频散和衰减现象,其中较为明显的是较大尺度裂缝,即l为1 m和较高含水饱和度,即sw为0.6,也反映出基于黏弹性Chapman-Kelvin模型的频变AVAZ裂缝储层参数反演,与裂缝尺寸和含水饱和度关系不唯一.

图2 刚度参数的实部和虚部随频率的变化

1.3 基于黏弹性Chapman-Kelvin模型的裂缝储层介质地震波频变特征的分析

裂缝储层中地震波的频散、衰减和频变各向异性,以及分析出重要的裂缝储层参数对地震波在裂缝介质中传播特性的影响,将为裂缝储层的储层参数反演提供有力的理论基础.在基于黏弹性Chapman-Kelvin模型中,我们研究了裂缝密度、裂缝长度、孔隙度、含水饱和度对地震波频变特征的影响.试验中考虑了一块被水和CO2气体饱和的岩石,体积模量和剪切模量代表不存在裂缝时的测量值,表1给出了各个参数的测量值.其中,附录A中给出了横波速度、横波分裂、横波衰减和准纵波速度的计算公式.

表1 模型计算中输入的参数值

裂缝密度是裂缝储层中的一个重要参数,它描述了储层中裂缝发育的程度.图3显示了裂缝密度对准纵波速度的影响.首先分别取频率为1 Hz、10 Hz、100 Hz、1000 Hz,从图3a可以看出:准纵波速度随裂缝密度的增大而增大.然后分别取裂缝密度为0.02、0.04、0.06和0.08,研究了在不同的裂缝密度下准纵波速度与频率的关系,从图3b可以看出,准纵波速度随频率的增大,先增大,后缓慢减小;准纵波产生的地震波速度频散的频率范围与裂缝密度的变化无关,而频散的幅度随着裂缝密度的增大而增大.

图3 准纵波速度随裂缝密度(a)和频率(b)的变化

裂缝与孔隙之间流体流动的弛豫时间τf与裂缝长度l成正比关系,因而裂缝长度将会影响频散出现的频段(刘喜武等,2015).图4显示了裂缝长度对横波分裂的影响.首先从图4a可以看出,横波分裂随裂缝长度的增大,而非线性增大,其中在低频即f=1 Hz和10 Hz时,变化的幅度小,但是在高频即f=100 Hz和1000 Hz时,变化的幅度变大.然后分别取裂缝长度为0.15 m、0.25 m、0.35 m和0.45 m,研究了不同的裂缝长度下横波分裂与频率的关系,从图4b可以看出,横波分裂产生的速度频散的频率范围不随裂缝长度变化,而频散的幅度随着裂缝长度的增大而增大.

图4 横波分裂随裂缝长度(a)和频率(b)的变化

图5显示了孔隙度对横波速度的影响.从图5a中,可以看出,横波速度随孔隙度的增大而增大,还可以发现横波速度随频率变化的敏感性较低.分别取孔隙度为0.05、0.1、0.15和0.2,研究了不同的孔隙度下横波速度与频率的关系,从图5b可以看出,横波产生的速度频散的频率范围不随孔隙度变化,而频散的幅度随着孔隙度的增大而增大.

图5 横波速度随孔隙度(a)和频率(b)的变化

图6显示了含水饱和度对横波衰减的影响.从图6a中研究发现,横波衰减在频率为1000 Hz时,随着含水饱和度的增大而减小,在其他频率下(1 Hz、10 Hz、1000 Hz),横波衰减随含水饱和度的增大,先增大后减小.为了研究不同含水饱和度下横波衰减随频率的变化,分别取含水饱和度为0、0.25、0.5、0.75和1,从图6b可以看出,含水饱和度为零时,横波衰减随频率先增大后减小;在含水饱和度不为零时,横波产生衰减的频率范围也与含水饱和度的变化无关,而衰减的幅度随含水饱和度的增大而增大.

2 基于黏弹性Chapman-Kelvin模型的AVAZ正问题与分析

可以利用地震属性检测裂缝导致的岩石弹性性质变化,如振幅随方位角和入射角的变化、地震波传播速度随方位的变化等(Liu et al.,2012;Zhang,2020).概况地说,目前可以分为两类,其一是根据AVO属性随方位角的变化进行椭圆拟合,定性估计裂缝发育方位和强度;另一种则是根据相应的裂缝等效介质理论,将反射的地震波振幅与偏移距和方位角联系起来,从而建立出反射系数与裂缝储层参数的关系,进而反演出裂缝储层参数(刘喜武等,2015).上述研究均是基于静态等效介质理论,其刚度系数矩阵与频率无关.并且在基于静态等效介质的AVO或AVAZ分析中,忽略了频散和衰减的影响.而动态等效介质模型可以描述具有频变特性的裂缝介质,如在1.1小节中的基于黏弹性Chapman-Kelvin动态等效介质模型.速度频散和衰减的出现将使得上覆地层与裂缝储层分界面的反射系数具有频变性,这种频变性与裂缝储层参数密切相关.因此,基于动态等效介质理论计算的频变反射系数与基于静态等效介质理论计算的不具有频变的反射系数相比,具有更为接近实际反射系数的值,即包含了更多裂缝储层信息.

本部分采用黏弹性Chapman-Kelvin动态等效介质模型描述裂缝孔隙介质,结合Schoenberg和Protázio(1992)概括的Zoeppritz方程,给出了频变反射系数的计算方法,并分析了重要的储层参数与频变反射系数的关系.最后,给出了正演方程的建立过程.

2.1 PP波和PS波反射系数的计算与频变AVAZ分析

R=(Y′-1Y+X′-1X)-1(Y′-1Y-X′-1X)

(4)

从而得到PP波和PS波的反射系数分别为R(1,1)和R(2,1),这个反射系数是入射角、方位角、频率以及储层参数矢量m的函数:

m=[l,fd,por,sw],

(5)

其中,l为裂缝长度;fd为裂缝密度;por为孔隙度;sw为含水饱和度.

为了研究裂缝孔隙介质与上覆介质分界面处的PP波和PS波反射系数具有频变性的同时,也随着入射角和方位角的变化,我们设计一个两层介质模型,即上层为页岩,并且假定为各向同性;下层为带有垂直裂缝排列的碳酸盐岩(各向异性,被80%水和20%CO2所饱和,体积模量、剪切模量等其他参数与表1中的参数一致),该模量的具体参数值,如表2所示.

表2 水和气饱和模型的参数

采用表1和表2中的参数,通过式(4)计算PP波和PS波的反射系数,并分析PP波和PS波的频变AVAZ特性.

图7给出了不同入射角情况下,PP波、PS波反射系数随方位角和频率的变化.当频率一定时,从图7a、b可以看出,PP波和PS波的反射系数幅值均随方位角的变化,而且这种变化随着入射角的增大变得更为显著;当固定方位角为50°,从图7c、d我们可以发现PP波反射系数幅值随频率的增大,先减小后增大;且PP波反射系数幅值在地震频段(0～100 Hz)内随频率的变化最为显著,尤其在大角度(40°)入射情况下;PS波反射系数幅值随频率的增大,先增大后减小;PS波反射系数幅值随入射角的增大而增大.

图7 不同入射角情况下PP波和PS波反射系数随频率和方位角的变化

2.2 储层参数变化对PP波频变反射系数的影响

研究裂缝储层参数(裂缝密度fd、裂缝长度l、孔隙度por和含水饱和度sw)的变化对PP波频变反射系数的影响.分析PP波频变AVO特性对这些储层参数的敏感性,以便于我们探究利用频变AVO方法反演这些参数的可行性.其中,选择表2所示的两层模型,并利用式(4)计算反射系数.

图8显示了裂缝密度变化对PP波频变反射系数的影响.分别取裂缝密度为0.02、0.05、0.08和0.1,其他模型参数保持不变.不同裂缝密度情况下,反射系数幅值均随着入射角的增大而增大,在入射角大于20°的范围内,这一现象较为显著;反射系数幅值随着裂缝密度的增大而增大.反映出PP波频变AVO特性对裂缝密度的变化较为敏感.

图8 不同裂缝密度的PP波反射系数幅值随入射角和频率的变化

图9显示了裂缝长度变化对PP波频变反射系数的影响.分别取裂缝长度为0.2 m、0.5 m、1 m和2 m,其他模型参数保持不变.我们可以发现频变反射系数的变化主要在入射角30°到45°之间,随着裂缝长度的增大,反射系数幅值随频率的变化幅度减小,而不同频率的反射系数幅值均有所增大(如图9所示);裂缝长度的增大导致不同入射角的反射系数幅值均有所增大,而反射系数幅值随入射角的变化幅度减小(如图9所示).反映出PP波频变AVO特性对裂缝长度的变化具有一定的敏感性,尤其在大角度即30°到45°、低频即10 Hz到40 Hz情况下,但其敏感性不如裂缝密度.

图9 不同裂缝长度的PP波反射系数幅值随入射角和频率的变化

图10显示了PP波频变反射系数随孔隙度的变化.分别取孔隙度为0.05、0.1、0.15和0.2,其他模型参数保持不变.从图10中可以看出,孔隙度的变化引起PP波频变AVO特性的显著变化,随着孔隙度的增大,不同频率、不同入射角对应的反射系数幅值均明显减小.从而说明,PP波频变AVO特性对孔隙度的变化非常敏感.因而,可以利用频变AVO反演技术较为准确地反演出孔隙度.

图10 不同孔隙度的PP波反射系数幅值随入射角和频率的变化

图11显示了PP波频变反射系数随含水饱和度的变化.分别取含水饱和度为0.2、0.4和0.6,其他模型参数保持不变.从图11可以看出,只有在低频时,反射系数幅值随着入射角的增大而增大,并且随着含水饱和度的增大,反射系数的幅值在不断地变小.所以,PP波频变AVO特性对含水饱和度的变化具有一定的敏感性,但要远远弱于裂缝密度、裂缝长度和孔隙度.

图11 (a) sw=0.2; (b) sw=0.4; (c) sw=0.6; (d) 不同含水饱和度的PP波反射系数幅值随入射角和频率的变化

2.3 正演方程的建立

传统的合成地震记录是通过震源与时间序列的反射率的褶积来生成,但是当地震频散发生时,反射系数与频率产生关系,使得我们在时间域上难以应用褶积.于是,我们考虑频率域中的地震子波(雷克子波)和傅里叶反变换方法来合成叠前地震记录,模型参数m和叠前地震数据的关系为:

d(θ,φ,t)=F-1[W(f)·R(θ,φ,f;m)]+n,

(6)

式(6)可以简化为:

d=F-1G(m)+n,

(7)

其中,dT=(d1,d2,…,dM)T为叠前观测地震数据;mT=(m1,m2,…,ml)T为模型参数;F-1表示傅里叶反变换;G为正演算子;W(f)为频率域中的地震子波;R为反射系数,如式(4)所示;f为频率;θ为入射角;φ为方位角;n为观测地震记录中的噪声.

3 裂缝孔隙介质中储层参数的反演

裂缝储层参数反演是裂缝储层评价的重要内容,对于裂缝油气藏的勘探和开发有着重要的意义.用于描述裂缝储层介质的传统裂缝等效介质理论是基于静态等效介质.由于没有考虑裂缝储层介质弹性响应的频变性,在含有少量较大尺寸裂缝的介质与含有较多小尺寸微裂隙的介质中,它们的弹性常数和各向异性参数可能相等,所以静态等效介质理论无法区别出地震各向异性是由微裂隙引起,还是由较大尺寸的裂缝引起(Maultzsch et al.,2003;Al-Harrasi et al.,2011).与微裂隙相比,一般来说,大尺度裂缝在裂缝储层中主导着流体流动.因此,裂缝尺寸是裂缝储层的一个重要参数,也是储层评价标准的一个重要参考.又根据2.2小节可以得知,动态等效介质模型所描述的地震响应的频变特征对裂缝密度、裂缝长度、孔隙度、含水饱和度等裂缝储层参数具有很好的敏感性.

在已知裂缝走向的条件下,利用贝叶斯理论,通过叠前频变AVAZ反演得到裂缝储层参数,即裂缝密度、裂缝长度、孔隙度和含水饱和度,以及显示出区别大尺度裂缝和微尺度裂隙的潜力.

3.1 贝叶斯理论

式(7)所示的反问题在地球物理领域中总是不适定的非线性问题,这种不适定性是由噪声干扰或不适当的正演算子等因素引起的(黄广谭,2019).因此,求解相应的反问题需要引入尽可能多的先验约束信息来缓解这种不适定性,这通常被称为正则化,本文我们采用贝叶斯理论来减轻这种不适定性问题.

贝叶斯反演基于概率和统计的思想,即利用先验信息和似然函数获得后验信息,从而把求取反演目标函数的极小值问题,转化为获得模型参数的后验概率分布的最大化问题.贝叶斯条件概率的公式为:

(8)

P(m|d)∝L(d|m)P(m).

(9)

3.2 裂缝孔隙介质的储层参数反演

基于统计理论的贝叶斯反演,通过引入模型参数的先验信息作为反演的正则化约束条件,能够有效地降低反演的非唯一性.模型参数的先验分布可以服从不同的分布函数(如均匀分布、高斯分布和柯西分布等),进而可以获得不同的反演结果(黄广谭,2019).本文选用多元高斯分布来构建先验模型,则有:

(10)

其中,Nd为模型参数的个数;m=[l,fd,por,sw]T～P(μm,Cm),l为裂缝长度,fd为裂缝密度,por为孔隙度,μm为模型参数的期望,即μm=E{m}=[μl,μfd,μpor,μsw]T;Cm为模型参数m的协方差矩阵,如式(11)所示:

(11)

其中,Cov表示协方差.

似然函数,即失配函数,是指由模型向数据的映射,描述了观测数据与合成数据之间的误差,表3给出了5种常用的失配函数(Huang et al.,2020).

表3 5种常用的失配函数

L1范数和L2范数度量作为地震反演中最常用的两种方法,常被用作失配函数.但是这两种方法都有其局限性.L1范数方法在零点处存在奇异性问题,L2范数方法因收敛性差而受到限制(黄广谭,2019).而Huber范数、混合范数和对数绝对范数能有效地克服常规范数的缺陷,其中对数绝对范数最为突出,对数绝对范数失配函数既能解决L1范数零点处的奇异性问题,又能尽可能保持与L1范数的收敛速度.

所以,我们在计算观测数据d(叠前地震数据)与正演模拟响应的残差时,引入了对数绝对范数,则残差表示为:

(12)

我们从残差中计算出似然函数,则似然函数表示为:

P(d|m)=aexp(-bΔE),

(13)

其中,a为归一化系数;b为常数.

式(13)是由Ulrych等(2001)分析得出的;Ulrych等(2001)和Mavko和Mukerji(1998)从理论上解释b的确定如何与观测的标准误差有关.然而Kirkpatrick等(1983)解释式(13)作为一个经验关系.

最终的似然函数表示为:

(14)

略去无关的归一化常量,可以得到m的后验分布如下:

P(m|d)∝exp[-b·[‖d-G(m)‖1

(15)

因此,目标函数J(m)可以简化为:

(16)

在得到后验概率分布之后,为了量化与裂缝相关的反演模型参数中的不确定性,可以通过蒙特卡罗随机采样方法得到反演结果.蒙特卡罗随机采样是通过大量重复试验得到某一事件的出现概率.在地震勘探领域中,一般认为地层参数具有马尔科夫性质,即在某一时刻t的模型参数的概率分布只与前一时刻(t-1)的模型参数有关.马尔科夫链是不同状态之间的转移,概率矩阵P决定转移的过程.其中,不同的转移矩阵的构建方法将产生不同的MCMC方法(Ali and Jakobsen,2014).因此,一般将蒙特卡罗随机采样方法和马尔科夫链结合起来进行采样.

为了有效地避免陷入局部极小的情况,在贝叶斯理论下,采用Metropolis准则对后验概率密度分布进行采样.Metropolis准则对任何样本都是以一定的概率接受或拒绝(Rubino et al.,2013).假设当前所在点为i,随机选择的下一个点为j,由i向j传递所遵循的Metropolis采样的接收准则为:

(17)

4 数值试验

4.1 基于贝叶斯理论的反演试验

本文使用表4中的参数进行合成研究,并关注于裂缝长度、裂缝密度、孔隙度和含水饱和度的估计以及在方位角AVO分析中,能够通过对地震频散的考虑,区分出大尺度裂缝与微尺度裂隙.本文的合成地震记录考虑一个具有两个界面的例子,即上层为砂岩,并且假定为各向同性;中间层为带有垂直排列裂缝的碳酸盐岩被60%的水和40%的CO2饱和,体积模量、剪切模量等其他参数与表1中的参数一致;下层为砂岩,并且假定为各向同性介质,具体的参数设定如表4所示.本文考虑了两种裂缝的情况,第一种是微尺度裂隙,即裂隙长度为1 mm,第二种是大尺度裂缝,即裂缝长度为1 m.利用公式(7)合成了角度域和方位角域的两种裂缝情况的地震道集记录,并在道集上加入了信噪比为4的随机噪声,信噪比SNR的定义为(Chen and Fomel,2015):

表4 水和气体饱和模型的参数 (合成示例)

(18)

合成的角度域和方位域的地震道集,如图12所示.其中标注框部分为考虑微尺度裂隙和考虑大尺度裂缝所合成的地震道集中主要的不同部分.

按照3.2小节中的反演步骤,本文首先假设裂缝长度l、裂缝密度fd、孔隙度por和含水饱和度sw满足高斯分布,均值分别取5、0.04、0.12、0.5,考虑到裂缝长度、裂缝密度、孔隙度以及含水饱和度在不同地区随岩性变化较大,因此选取较大的方差,即分别取2.5、0.25、0.4、0.35,并将裂缝长度、裂缝密度、孔隙度和含水饱和度控制在一定范围内,即l∈[0,10],fd∈[0,0.1],por∈[0,0.2],sw∈[0,1],先验概率分布函数如图13所示,并假定所有其他参数都已知.在实际应用中,这些先验信息可以从测井分析或地震解释中得到(Aster et al.,2005).

通过使用式(12)求得观测数据和正演模拟数据的残差,再通过式(14)得到似然函数.然后根据式(15)计算出后验概率分布P(fd,l,por,sw|d).对于给定问题的概率分布P(x,y,z),每个变量的边际概率分布为:

(19)

(20)

(21)

因此,可以从边际概率的角度出发,直接从后验分布P(fd,l,por,sw|d)中,计算出的裂缝密度fd、裂缝长度l、孔隙度por、含水饱和度sw的边际概率分布的表达式,分别为:

(22)

(23)

(24)

(25)

我们分别应用于大尺度裂缝和微尺度裂隙情况,即微尺度裂隙,取裂隙密度为0.05,裂隙长度为1 mm,并且在该种情况下孔隙度取0.1,含水饱和度取0.6;大尺度裂缝,取裂缝密度为0.05,裂缝长度为1 m,并且在该种情况下孔隙度取0.1,含水饱和度取0.6.图14、图15、图16、图17分别为后验概率P(fd,l,por,sw|d)的裂缝密度、裂缝长度、含水饱和度和孔隙度的边际概率分布.从图14的裂缝密度边际概率分布来看,裂缝密度在大尺度裂缝下反演结果的最大概率分布集中在0.05处,与真实的裂缝密度很吻合,而在微尺度裂隙下的反演结果的最大概率分布集中在0.04,很显然与真实的裂缝密度相差很大.从图15的裂缝长度边际概率分布来看,微尺度裂隙下裂缝长度的最大概率主要集中在10-4m和10-2m之间,不能准确地估计出真实微尺度裂隙的长度,而大尺度裂缝的最大概率在100m处,与真实大尺度裂缝的长度相符合.从图16的含水饱和度边际概率分布来看,大尺度裂缝下含水饱和度的最大概率主要集中在0.6处,与真实含水饱和度值相吻合,而微尺度裂隙下含水饱和度的最大概率集中在0.82处,很显然不符合真实值.最后,从图17的孔隙度边际概率分布来看,大尺度裂缝下孔隙度的最大概率分布集中在0.1和0.12区间,而微尺度裂隙下孔隙度的最大概率集中在0.08处,相对大尺度裂缝下的结果来说,精确度不够.

图14 微尺度裂隙(a)、大尺度裂缝(b)下裂缝密度的后验概率分布

图15 微尺度裂隙(a)、大尺度裂缝(b)下裂缝长度的后验概率分布

图16 微尺度裂隙(a)、大尺度裂缝(b)下含水饱和度的后验概率分布

图17 微尺度裂隙(a)、大尺度裂缝(b)下孔隙度的后验概率分布

通过求解裂缝密度、裂缝长度、含水饱和度和孔隙度等储层参数的后验概率分布,反演出各个参数最大概率的分布情况,可以很好地区分开大尺度裂缝和微尺度裂隙两种情况.也可以发现,我们可以精确地反演出大尺度裂缝的长度,而不能够准确地恢复微尺度裂隙的长度,这是因为地震频散受到了裂缝长度的控制,当裂缝长度足够小时,即微尺度裂隙,地震频散可以忽略不计,所以对于地震频散可忽略不计的微尺度裂隙情况,不能准确地反演出裂缝密度、含水饱和度和孔隙度.因此,在方位角AVO分析中对地震频散的考虑,显示出将大尺度裂缝和微尺度裂隙区分开的潜力.

4.2 基于频变反射系数的MCMC反演试验

基于公式(7)的正演方程,若不考虑地震子波的作用,直接利用反射系数,并增加约束,对储层参数进行反演.我们将储层参数集M设置为:

M=[fd,l,por,sw],

(26)

即我们关注的储层参数主要为裂缝密度fd、裂缝长度l、孔隙度por和含水饱和度sw.因此,随频率变化的P波反射系数可以表示成:

R(θ,φ,f;M),

(27)

其中,θ和φ分别为入射角和方位角,f为频率.

对于K个入射角(k=1,2,…,K),Z个方位角(z=1,2,…,Z),F个频率(f=1,2,…,F),反射系数R可以写作为:

Rkzf(M)=R(θk,φz,ff;M),

(28)

将反射系数进一步写成矩阵形式:

R=GM,

(29)

其中,G为基于Chapman-Kelvin模型和Schoenberg和Protázio(1992)概括的Zoeppritz方程所计算的频变反射系数的核函数.

基于表4中的模型参数,且此时考虑大尺度裂缝的情况,即我们假定真实裂缝密度为0.05,真实裂缝长度为1 m,真实孔隙度为0.1,以及真实含水饱和度为0.8.然后,我们尝试使用计算得到的合成地震AVAZ资料来估计真实值,即反射系数的大小作为不同入射角的地震频率(≤100 Hz)和方位角的函数.使用1.1小节中描述的黏弹性Chapman-Kelvin模型和附录B中的阻抗矩阵(Schoenberg and Protázio,1992)来计算反射系数的大小.假设入射角为10°、20°、30°、40°,方位角为20°、40°、60°、80°,频率为10 Hz、20 Hz、30 Hz、40 Hz、50 Hz、60 Hz、70 Hz、80 Hz、90 Hz、100 Hz.当缺少先验信息时,式(16)中的目标函数J(m)变为:

(30)

图18 不同的入射角下添加信噪比为5的随机噪声的PP波反射系数

图19显示了在没有关于任何储层参数的先验信息的情况下,使用频变的地震AVAZ数据,以储层参数的后验概率密度函数形式得到的MCMC采样的反演结果.从图19可以看出反演结果的不确定性很强,无法得到储层参数的准确解.

图19 无先验信息的储层参数的反演结果

我们接下来通过假设一个关于裂缝密度的先验信息来进行反演.本文我们将裂缝密度、裂缝长度、孔隙度和含水饱和度的目标函数J(m),以及关于裂缝密度的先验信息定义为:

(31)

其中,第一个平方项是在4个不同的入射角、4个不同的方位角和10个频率(≤100 Hz)上反射系数的观测值和计算值之间的不匹配;第二个平方项表示裂缝密度的先验信息,fd0为裂缝密度先验信息的均值;Δfd为裂缝密度先验分布的标准差.

图20显示了利用频变的地震AVAZ数据和裂缝密度的先验信息对储层参数的后验概率密度函数形式的MCMC反演结果.我们可以清楚地看到,包含裂缝密度的先验信息可以在一定程度上提高裂缝长度、孔隙度和含水饱和度的确定性,同时裂缝密度的反演结果的不确定性也大大减小,但含水饱和度的确定性仍然较差.

同样,假定已知裂缝长度和孔隙度的先验信息.可以通过观测和钻孔数据的插值得到裂缝长度的先验信息.Odling等(1999)讨论了一种从露头工区中求取裂缝长度的方法,在这种方法中,裂缝长度总体绘制为归一化的累积频率分布,即T(l),其中l为裂缝长度.当使用对数轴时,指数为c,裂缝长度的分布如下:

T(l)∝l-c.

(32)

Maultzsch等(2003)通过使用频率相关的剪切波各向异性和P波衰减的反演也可以获得关于裂缝长度的先验信息.

我们将裂缝密度、裂缝长度、孔隙度和含水饱和度的目标函数J(m),以及关于裂缝长度和孔隙度的先验信息定义为:

(33)

其中,第一个平方项与式(31)中的定义一致;第二个平方项分别表示裂缝长度和孔隙度的先验信息,l0、por0分别为裂缝长度和孔隙度先验信息的均值;Δl、Δpor分别为裂缝长度和孔隙度先验分布的标准差.

图21显示了利用频变地震AVAZ数据和裂缝长度和孔隙度的先验信息,对储层参数的后验概率密度函数形式的MCMC反演结果.从图21a的裂缝密度反演结果来看,相对于没有考虑先验信息的图19a的裂缝密度反演的准确度有了进一步的提高.但是我们对比图20d、d,发现关于裂缝长度和孔隙度的先验信息对含水饱和度的估计没有太大的改善作用.

图21 具有裂缝长度和孔隙度先验信息的储层参数的反演结果

最后,我们针对具有裂缝长度和孔隙度先验信息的含水饱和度反演结果确定性仍然不是很理想,所以我们考虑了足够多的先验信息,即假定裂缝长度、裂缝密度和孔隙度的先验信息是已知的.本文将裂缝密度、裂缝长度、孔隙度和含水饱和度的目标函数J(m),以及关于裂缝长度、裂缝密度和孔隙度的先验信息定义为:

(34)

其中,第一个平方项与式(31)中的定义一致;第二个平方项分别表示裂缝长度、裂缝密度和孔隙度的先验信息,l0、fd0、por0分别为裂缝长度、裂缝密度和孔隙度先验信息的均值;Δl、Δfd、Δpor分别为裂缝长度、裂缝密度和孔隙度先验分布的标准差.

图22显示了利用频变地震AVAZ数据和裂缝长度、裂缝密度和孔隙度的先验信息对储层参数的后验概率密度函数形式的MCMC反演结果.从图22d含水饱和度的反演结果来看,相比之前,反演结果的确定性有了一定的加强.并且通过对比图21和图22中含水饱和度的反演结果,我们也不难发现关于裂缝密度、裂缝长度、孔隙度三者的先验信息最终对含水饱和度的估计变得更加关键.

图22 具有裂缝长度、裂缝密度和孔隙度先验信息的储层参数的反演结果

5 讨论

裂缝等效介质模型能够更精确地模拟裂缝储层介质的性质,本文根据所提出的基于黏弹性Chapman-Kelvin动态等效介质模型,其能够适应复杂的储层情况,并初步讨论了裂缝储层的储层参数反演问题,但是离实际裂缝储层中工业应用还有很远的距离.综合国内外对裂缝储层相关问题的研究现状,对其进一步的研究应包括以下几个方面:

(1)完善裂缝储层动态等效介质模型理论.在Chapman(2003)模型的基础上,考虑了实际地下介质的黏弹性和地震频散,引入了线性流变Kelvin黏弹性理论模型,构建出黏弹性Chapman-Kelvin动态等效介质模型.但是能否将线性流变理论模型与基于Boltzmann理论的黏弹性模型统一结合到Chapman(2003)模型中,以便适应更加复杂的黏弹性储层情况还有待后续研究.

(2)裂缝储层参数反演由于各向异性、地震资料的处理、参数的选择以及反演策略的原因具有很强的多解性,尤其在考虑地震频散的情况下,表现得更为显著.本文基于贝叶斯理论进行裂缝储层参数的反演,在缺失先验的储层参数信息时,反演结果的不确定性较强.当存在有效且足够的先验信息时,反演结果的可靠性进一步提升.但是,关于储层参数的先验信息很难获取.

(3)本文所提出的频变AVAZ反演技术,为裂缝储层参数的反演和区别大尺度裂缝和微尺度裂隙提供了新的思路.在实际应用中,需要宽方位地震采集数据和高保真去噪技术来获得高质量的叠前数据,采用多尺度地震分频技术,按照本文的反演流程进行大尺度裂缝和微尺度裂隙的区分,以及裂缝储层参数的反演,并结合实际的测井资料,验证所提出的反演方法的实际应用效果.

6 结论

裂缝储层介质的等效介质模型建立了地震数据与储层特征的联系.考虑到黏弹性介质比弹性介质更能代表地球的性质,以及传统的裂缝等效介质模型无法对实际裂缝储层中地震波的频散和衰减给出合理的解释.因此,本文基于已经提出的黏弹性Chapman-Kelvin动态等效介质模型,并在此基础上就裂缝密度、裂缝长度、孔隙度以及含水饱和度对地震波频变特征的影响做了分析.通过理论分析和数值试验,我们可得到以下结论和认识:

(1)利用黏弹性Chapman-Kelvin动态等效介质模型与Schoenberg和Protazio概括的Zoeppritz方程计算了反射系数,分析了反射PP波和PS波的频变AVAZ特性,即当频率一定时,PP波和PS波的反射系数幅值均随方位角的变化,而且这种变化随着入射角的增大变得更为显著;当固定方位角为50°,我们可以发现PP波反射系数幅值随频率的增大,先减小后增大;且PP波反射系数幅值在地震频段(0～100 Hz)内随频率的变化最为显著,尤其在大角度入射情况下;PS波反射系数幅值随频率的增大,先增大后减小,而随入射角的增大而增大.并就重要的储层参数对PP波频变反射系数的影响进行了分析,发现孔隙度的敏感性最强,其他依次为裂缝密度、裂缝长度和含水饱和度.

(2)基于PP波频变反射系数对裂缝储层参数变化有较好的敏感性,开展了两种反演方法研究,其一是基于贝叶斯理论直接反演方法,其中以对数绝对范数作为似然函数和高斯分布即L2范数度量作为先验约束;其二是基于MCMC的随机反演方法.通过频变AVAZ进行储层参数的反演研究,结果表明,MCMC随机反演方法在缺失先验信息时,反演解的不确定性较强;包含裂缝密度的先验信息在一定程度上提高裂缝长度、孔隙度和含水饱和度的确定性,但含水饱和度的确定性仍然很差;假定裂缝长度、裂缝密度和孔隙度三者先验信息是已知时,含水饱和度反演的确定性有了进一步的提高.而基于贝叶斯理论直接反演方法,由于考虑到地震频散的影响,通过后验概率的边际概率分布情况,区别出了大尺度裂缝和微尺度裂隙,使得利用频变方位AVO来区分微尺度裂隙和大尺度裂缝成为可能.

致谢感谢审稿专家提出的宝贵意见.

附录A 基于黏弹性Chapman-Kelvin模型的横波速度、横波分裂、横波衰减和准纵波速度的计算

基于黏弹性Chapman-Kelvin地震波传播的模型,如式(A1)所示:

(A1)

式(A1)写成矩阵形式为:

(A2)

其中,子块矩阵:

(A3)

mfw=C44+Cmu,

(A4)

其中:

Cmu=2Udry-real(C44),

(A5)

其中,Udry为干燥岩石的剪切模量;real指复数的实部部分.

则横波速度vs的表达式为:

(A6)

其中,r为岩石密度,横波衰减qis的表达式为:

(A7)

其中,imag指复数的虚部部分,横波分裂sws的表达式为:

(A8)

其中,vps为纯横波的速度,表达式为:

(A9)

其中,α为地震波传播和对称轴之间的角度(0°对应于沿裂缝的法线传播);且vqs为准横波速度,其表达式为:

(A10)

其中,H为常数,表达式为:

H=[(C11-C44)sin2α-(C33-C44)cos2α]2

+(C13+C44)2sin2(2α).

(A11)

最后,准纵波速度vqp的计算表达式为:

(A12)

附录B 基于Schoenberg和Protázio(1992)阻抗矩阵的计算

Schoenberg和Protázio(1992)利用Zoeppritz方程系数矩阵的子矩阵给出了平面波反射和透射问题的显示解;Chapman和Liu(2003)在Schoenberg和Protázio(1992)的研究基础上,推导了与入射角和方位角相关的反射率;Jin等(2017)假设各向异性介质有一个平行于X1-X2平面的水平对称平面,设X3=0是两个半空间之间的水平界面.设定上部弹性介质中入射波的入射角为θ,方位角为φ,传播方向可以表示为n=(sinθcosφ,sinθsinφ,cosθ).当地震波在垂直X1-X3平面内传播时,位移可以表示成:

(B1)

其中,e1、e2和e3表示极化,s1为X1方向上的水平慢度,s2为X2方向上的水平慢度,s3为垂直慢度;水平慢度s1和s2分别为式(B2)、(B3)所示:

s1=ξsinθcosφ,

(B2)

s2=ξsinθsinφ,

(B3)

Tsvankin(2012)给出了平面波的解析表达:

(B4)

其中,μ为位移;Cijkl为弹性刚度张量.将式(B1)代入式(B4),可以得出:

ρμi=Cijklslsjμk,

(B5)

其中,Cijklslsj定义为Christoffel矩阵(Mavko et al.,2020),该矩阵建立了弹性模量与水平慢度和垂直慢度之间的联系.

将式(B5)进行等价改写成:

(ρδik-Cijklslsj)μk=0,

(B6)

其中,若i等于k,则δik为1;若i不等于k,则δik为0.

通过求解式(B7)中的方程,可以得到垂直慢度s3与弹性模量Cijkl之间的关系:

|ρδik-Cijklslsj|=0,

(B7)

其中,特征值分别是准P波、第一个到达的准S波和最后一个到达的准S波的垂直慢度s3p、s3s和s3t的平方解,相关的特征向量为相应的极化ep、es和et.

上部介质中的地震波场包括入射和反射的P波和S波,而下部介质只有透射波(Jin et al.,2018).上部介质的位移场可以写成:

(B8)

其中,ip,rp,is,rs,it和rt分别是入射准纵波、反射准纵波、入射快准横波、反射快准横波、入射慢准横波和反射慢横波的幅值;ep1,ep2,ep3;es1,es2,es3;et1,et2,et3分别为准P波、第一个到达准S波和最后一个到达准S波的极化;s3p,s3s,s3t分别为准P波、第一个到达准S波和最后一个到达准S波的垂直慢度;ω为角频率.

而在下层介质中只有透射的P波和S波,其位移场表示为:

(B9)

其中,tp,ts和tt分别为透射的准纵波、透射的快准横波和透射的慢准横波的幅值;所有其他带撇号的参数具有与式(B8)中不带撇号的参数相同的含义.

可以看出反射系数计算的关键是求出垂直慢度相应的极化值.由Christoffel方程可以从给定的水平慢度条件下,计算出垂直慢度和相应的极化.Zoeppritz方程给出了在两个弹性介质之间的界面上,入射的纵波和横波的反射和透射平面波振幅的精确计算.Schoenberg和Protázio(1992)通过引入两个阻抗矩阵来求解Zoeppritz方程,得到了水平和垂直慢度值,以及相应的极化值,两个阻抗矩阵的表达式为:

(B10)

(B11)

其中,Cij表示入射介质实数域的刚度矩阵元素.