碳酸盐岩缝洞型储层波场数值模拟

刘炯魏修成陈天胜（中国石化石油勘探开发研究院，北京100083）

摘要缝洞型碳酸盐岩储层作为当今最重要的储层类型之一，一直是油气勘探开发领域的研究重点。但是这类储层往往结构比较复杂，因而地震预测比较困难。本文利用裂缝等效介质理论和随机介质理论建立了在定向裂缝介质中含有随机分布孔洞的各向异性随机模型来描述碳酸盐岩缝洞型储层，并采用伪谱法模拟地震波在模型中的传播。模拟结果显示：由于定向裂缝的作用地震波在缝洞型储层中传播时，地震波会出现横波分裂的现象；另外地震波会因储层中的孔洞发生散射，使得空间波场变得复杂。

关键词裂缝孔洞各向异性随机介质

WavefieldSimulationinCarbonateKarstReservoirs

LIUJiong，WEIXiucheng，CHENTiansheng（SINOPECExploration＆ProductionResearchInstitute，Beijing100083，China）

AbstractCarbonatekarstreservoirisoneofthemostimportantreservoirtypesintheworld，andithasbeenaresearchcenterinoilandgasexplorationanddevelopment.Becauseofthecomplexstructure，karstreservoirisdifficulttobepredictedbytraditionalseismicexplorationtechnology.Basedonfractureequivalentmediatheoryandrandommediatheory，theanisotropicrandommodelissetuptodepictcarbonatekarstreservoir，inwhichcavesandalignedfracturesarecommon.Thenpseudospectralmethodisusedtosimulateseismicwavepropagationinthistypeofreservoir.Resultsshowthatwhenseismicwavepropagateinkarstreservoirshearwaveswillsplitbecauseofalignedfractures.Andonthesurfacesofcaves，seismicwavewillscatter，whichmakesthewavefieldcomplicatedinspace.

Keywordsfracture；cave；anisotropy；randommedia

随着勘探程度的提高，大型构造型油气藏越来越少，勘探目标开始转向复杂储层，如向含裂缝型储层以及裂缝、孔洞结构并存的碳酸盐岩介质转移。

裂缝作为一种复杂的空间结构，大量存在于岩石、地层中。大量的油气勘探实践表明，在储存空间中的裂缝是流体运移的通道，直接关系到油气的产量，同时裂缝在许多储层中也是油气储层的空间，影响储层的油气含量。许多学者对裂缝进行了大量的研究。20世纪80年代，Crampint［1］通过研究发现，地震波在定向裂隙介质中传播时和波在各向异性介质中的传播等效，都会出现快横波和慢横波分裂的现象，并将含定向裂隙的介质称为广泛扩容性各向异性EDA（extensivedilatancyanisotropy）介质。对于一般岩石EDA介质中的众多小裂缝，Hudson［2，3］将它们看成是一个个非常扁的椭球体，并用弹性扰动理论推导出裂缝等效各向异性介质的弹性系数与各向同性背景介质的弹性系数、裂缝参数之间的关系，还给出了裂缝中不同充填物对弹性常数的影响。Schoenberg和Sayers［4］将裂缝看成是具有线性滑动边界条件的柔性边界，推导出了裂缝等效各向异性介质的柔性矩阵。随后很多学者运用这两种等效介质理论研究了地震波在裂缝介质中的传播特点。在这些研究中，人们主要从裂缝角度来考虑对地震波的影响，但是实际地层不仅包含裂缝，还可能含有大量尺度不等的其他结构，如孔、洞等。

地层中的孔洞作为区别于背景介质的不均匀结构，往往使得地震波在界面上发生散射、衰减，从而使得波场变得复杂。许多学者对地层中不规则孔洞做了许多工作，研究了它们对地震波传播的影响［5，6］。在以往的孔洞研究中，学者往往用确定性的方法来描述孔洞在空间中的位置，但是对于实际中大量存在的尺度比波长小很多的孔洞，用空间随机分布的方法去表述更为合理。随机模型起源于20世纪60年代。Aki［7］提出了地下介质的随机不均匀性引起的尾波是导致地下振动长久持续的主要原因。Berteusen［8］将随机介质中标量波散射理论应用于大孔径台阵远震P波记录的相位、振幅涨落问题研究。Frankel和Clayton［9］的方法是用数值研究随机弹性介质中纵波的散射衰减。Liu［10］等运用孔隙介质理论和随机理论建立了随机孔隙介质模型，并研究了地震波在其间传播的能量衰减。

前人对单独的裂缝结构和孔洞结构储层已经做了许多研究。但是在实际地层中，地质结构往往不是单一的，如海相碳酸盐岩储层中，由于地质作用，裂缝、孔洞同时存在。为此本文运用Hudson裂缝等效介质理论和随机介质理论建立裂缝介质中含有空间分布小孔洞的各向异性随机介质模型，并采用伪谱法来模拟地震波在该模型中的传播，以此来认识地震波在含微小裂缝和宏观小孔洞的碳酸盐岩储层中的传播规律。

1各向异性随机介质模型的原理和方法

本部分首先介绍裂缝等效各向异性介质的Hudson理论，然后对弹性波在二维方位各向异性介质中的传播方程做了推导，最后阐述了空间随机介质的建模方法。

1.1Hudson裂缝等效介质理论

20世纪80年代，Hudson在长波近似、地震波场范围内裂纹位置分布均匀、裂纹在岩石空间中稀疏且彼此不连通的假设前提下，得到了小纵横比扁球体裂缝性质同岩石整体宏观性质之间的关系。在Hudson理论中，含小裂缝的岩石等效弹性常数可以表示成如下形式：

油气成藏理论与勘探开发技术：中国石化石油勘探开发研究院2011年博士后学术论坛文集.4

式中：是各向同性背景介质的弹性常数；是由于裂缝存在而产生的一阶、二阶修正。

对于垂直裂缝组，裂缝介质显示出横向各向同性的对称性，其总体弹性参数矩阵可以表示为

油气成藏理论与勘探开发技术：中国石化石油勘探开发研究院2011年博士后学术论坛文集.4

在Hudson理论中，式（2）的各弹性常数的表达式如下：

油气成藏理论与勘探开发技术：中国石化石油勘探开发研究院2011年博士后学术论坛文集.4

式中：λ、μ是各向同性背景介质的拉梅常数；ε表示裂缝密度；参数可以表示成

油气成藏理论与勘探开发技术：中国石化石油勘探开发研究院2011年博士后学术论坛文集.4

垂直裂缝弹性常数表达式中的U1、U3依赖于裂缝内的充填，本文考虑裂缝干燥含气时的状态，表达式如下：

油气成藏理论与勘探开发技术：中国石化石油勘探开发研究院2011年博士后学术论坛文集.4

1.2二维方位各向异性介质中地震波的传播

裂缝介质总体弹性系数矩阵（2）是在自身本构坐标系下的表达形式。当裂缝在地质构造作用下发生变化，如层位发生倾斜，即裂缝所在的本构坐标系和观测坐标系存在一定的倾角时，观测坐标系下的弹性矩阵形式会有变化，可以由本构坐标系下的弹性矩阵通过Bond变换获得，其关系如下：

油气成藏理论与勘探开发技术：中国石化石油勘探开发研究院2011年博士后学术论坛文集.4

式中：C、C′分别是裂缝介质在本构坐标系和观测坐标系下的弹性系数矩阵；M是坐标转换的Bond矩阵，M′是M的转置。

由式（12）得到裂缝方位各向异性介质的弹性矩阵，其元素一般都不为0。在此基础上运用弹性介质理论，将应变-位移关系代入关于C′的本构方程中，并将е/еy取为0，可以得到x-z平面内的应力-位移关系如下：

油气成藏理论与勘探开发技术：中国石化石油勘探开发研究院2011年博士后学术论坛文集.4

油气成藏理论与勘探开发技术：中国石化石油勘探开发研究院2011年博士后学术论坛文集.4

式中：σx、σz分别是介质沿x、z方向的正应力；τyz、τxz、τxy是介质的剪应力；ux、uz是空间质点沿x、z轴的位移分量；是一般方位各向同性介质的弹性常数。

一般方位各向同性介质中，x-z面内质点的运动方程如下：

油气成藏理论与勘探开发技术：中国石化石油勘探开发研究院2011年博士后学术论坛文集.4

式中：ρ表示介质的密度。

方程（13）～（20）构成了方位各向同性介质x-z平面内弹性波传播的控制方程。

1.3随机介质理论

以往波传播问题的研究多考虑均匀或分层均匀的情况，而实际介质往往是非均匀的。对于实际地下大量存在而且分布不规则的异常介质往往用随机介质模型来描述更接近真实情况。

根据随机过程理论，任意二维空间随机分布量f可以表示成如下平均值和扰动量之和的形式［11］：

油气成藏理论与勘探开发技术：中国石化石油勘探开发研究院2011年博士后学术论坛文集.4

式中：f0表示f的空间平均值，它是常数；γ（x，z）是在点（x，z）处f相对于平均值的扰动。为了数学上的处理方便，假设空间随机扰动γ（x，z）是均值为0的空间平稳随机过程，即：

油气成藏理论与勘探开发技术：中国石化石油勘探开发研究院2011年博士后学术论坛文集.4

除了均值人们还往往用方差σ2和自相关函数φ来描述平稳随机过程，它们的表达式如下：

油气成藏理论与勘探开发技术：中国石化石油勘探开发研究院2011年博士后学术论坛文集.4

上述方程中〈·〉表示空间平均算子。

根据随机过程理论，γ（x，z）的功率谱就是其自相关函数φ（x，z）的傅立叶变换，所以可以用随机过程的自相关函数用谱展开的方法来构建γ（x，z）空间随机分布。在构建随机介质的过程中自相关函数φ（x，z）的选择有多种，如高斯型、指数型、VonKarman型，本文采用指数型自相关函数，其形式如下：

油气成藏理论与勘探开发技术：中国石化石油勘探开发研究院2011年博士后学术论坛文集.4

其中a和b分别是随机介质在x和z方向上的自相关长度。

对于均匀背景介质中含孔洞随机分布的模型，可以按以下的方法建立：首先选取指数型函数（25）作为某空间分布量的自相关函数，并选取自相关函数中相关长度，然后用谱展开方法得到函数值的空间随机分布。选取某值作为阈值，当空间某点的分布量大于阈值时认为该点为背景介质区域，否则就是孔洞区域。

2平面各向异性随机模型的建立和地震波的数值模拟

本节首先建立同时反映裂缝、孔洞性质的各向异性随机模型，并通过地震波波动方程数值模拟的方法，来直观认识地震波在含裂缝、孔洞复杂碳酸盐岩储层中的传播规律。然后建立含碳酸盐岩储层随机各向异性介质的两层地质模型，通过记录地表地震的方法来认识随机各向异性介质模型在地震反射记录上的特点。为了方便认识裂缝、孔洞对地震波的作用，我们对两种模型对应的只含裂缝的情况也做了相应的数值模拟。

2.1二维各向异性随机模型及其波场模拟

根据前面介绍的Hudson理论建立含微小裂缝的等效各向异性介质，然后以等效的各向异性介质为背景，用谱展开理论建立了含孔洞的各向异性随机介质模型。

为了描述孔洞空间的大小，本文借用孔隙度的概念定义孔洞孔隙度φ如下：

油气成藏理论与勘探开发技术：中国石化石油勘探开发研究院2011年博士后学术论坛文集.4

式中：Vc表示孔洞的体积；Vt是整个模型空间的体积。

2.1.1二维各向异性随机模型

取自相关函数（25）中的相关尺寸a＝b＝15m，建立二维各向异性随机模型。如图1所示，背景是含裂缝的等效各向异性介质，其中裂缝主轴水平，偏离观测坐标系x轴-45°。在含裂缝的介质中，背景各向同性介质的密度ρ＝2510kg/m3，拉梅常数λ＝7.7351×109Pa，μ＝23.044×109Pa，裂缝体积密度ε＝0.05，其内部假设为干的状态。图1中散布的点代表小尺度孔洞，孔洞中包含物设定为各向同性的泥浆，密度ρc＝1580kg/m3，其弹性参数为λc＝0.5432×109Pa，μc＝1.3665×109Pa，孔洞所占总空间的孔隙度为0.1％。

图1各向异性随机模型

2.1.2各向异性随机模型的波场模拟

伪谱法作为波动方程数值模拟的方法之一，主要通过傅立叶变换的方法把物理变量对空间的微分转化为空间频率域中的代数运算，然后再把结果通过傅立叶逆变换转换到物理空间，从而求得对应量的空间微分值［12，13］。理论上其精度可以和有限差分、空间差分精度达到无穷时的情况相当［14］。在本文研究中我们采用伪谱法来模拟地震波在各向异性随机模型中的传播。

首先对图1所示的二维各向异性随机介质空间采用dx＝dz＝15m的正方形网格进行离散，然后采用伪谱法数值求解。模拟中时间步长取为dt＝0.5ms，震源采用心频率取为50Hz的雷克子波，以x方向集中力源的形式安置在模型的中心点。图2是地震波在各向异性随机模型中传播0.325s时x、y、z3个方向位移分量的波场快照。为了更好地显示孔洞结构对整体波场的影响，我们对只含方位水平裂缝介质的情况也进行了模拟，对应时刻的三分量位移快照如图3所示。

图2各向异性随机介质中的位移波场快照

从图2、图3的波场快照可以看出准纵波qP、快横波qS1、慢横波qS2的传播，这表明伪谱法可以很好地模拟地震波在各向异性随机介质和裂缝介质中的传播。模拟的结果也表明，由于不同偏振方向的横波传播速度不同，像在裂缝介质中一样，横波在各向异性随机介质中也会产生分裂的现象。与在裂缝介质中不同的是，地震波在各向异性随机模型传播过程中还会出现很多杂乱的散射波。这是因为波遇到孔洞会发生散射，此时每一个小孔洞等效于次生震源，并以此为中心产生与均匀各向异性介质中相似的次生波场。由于模型中含有多个孔洞，多个孔洞产生的次生波场与原波场相互叠加，导致空间总波场变得复杂。

图3裂缝介质的位移波场快照

2.2碳酸盐储层各向异性随机模型的地震记录模拟

为了观测缝洞型碳酸盐岩储层在地震记录上的特点，我们设计了含碳酸盐岩储层的两层地质模型，通过伪谱法来模拟地震波在其间的传播，并在地表布置检波器来接收三分量地震位移记录。

2.2.1含碳酸盐储层的各向异性随机模型

设计含孔洞型碳酸盐岩储层的两层介质模型如图4所示。其中第一层从0至960m的深度，为碳酸盐岩的各向异性随机介质，其裂缝的对称轴水平，偏离观测坐标系-45°夹角。背景介质的密度ρ1＝2400kg/m3，λ1＝3.287×109Pa，μ1＝17.496×109Pa，裂缝体积密度ε＝0.05，孔洞内物质参数同2.2中的各向异性模型，但随机孔洞的孔隙度取为0.05％。第二层从690m至1755m深度，为均匀各向同性介质，其密度ρ2＝2650kg/m3，λ2＝6.731×109Pa，μ2＝32.463×109Pa。

图4含碳酸盐岩储层的两层地质模型

图5各向异性随机两层模型的三分量地震记录

2.2.2碳酸盐储层各向异性随机模型的地震记录模拟

用伪谱法对地震波的传播进行数值模拟。模拟中用dx＝dz＝15m的正方形网格进行空间离散，时间步长dt＝0.5ms；纵波震源采用50Hz的雷克子波，其中心位于x＝870m，z＝150m的空间点上。在地表每隔15m布置一个检波器接收位移三分量地震记录。整个模拟时间取为0.7s。为了消除人工边界的影响，在模型四周加了完全匹配层吸收边界条件。得到的位移地震记录如图5所示。为了更好地显示各向异性随机模型的反射特征，我们对同样两层模型但第一层是不含孔洞的裂缝各向异性介质的情况也做了模拟，其地震记录如图6所示。

图6各向异性两层模型的三分量地震记录

从图5、图6的直达准纵波记录qP和直达准横波qS可以看出，虽然在模拟中采用纵波震源激发，但激发的纵波也会产生横波，这是由于在方位裂缝各向异性介质中不同方向的运动是相互耦合的。从模拟结果还可以看出各向异性模型记录中出现了反射准纵波qPqP、反射快横波qPqS1和反射慢横波qPqS2，而碳酸盐岩各向异性随机模型的地震记录上也出现了同样的现象。这是因为准纵波在向下传播时遇到不同介质的分界面，一部分准纵波发生反射并向上传播至地表形成反射准纵波qPqP，一部分准纵波在下行传播过程中遭遇界面时会发生波型转换产生反射横波。由于裂缝结构的作用，反射的横波会分裂成快横波qPqS1和慢横波qPqS2，并先后被地表检波器所接收。从模拟结果还可以看出，在各向异性随机介质的地震记录中出现明显散射现象，而裂缝各向异性介质的地震波记录却没有出现这种现象。这是因为准纵波在传播遭遇到孔洞时，会发生反射和波型转换，产生反射准纵波和反射准横波，这些叠加在裂缝介质的地震记录上，形成各向异性随机介质的复杂地震记录。当孔洞的反射很强时，孔洞的反射波可能会掩盖其他的反射波，如在图5的z方向位移记录上，由于孔洞的影响，反射慢横波变得不明显。

3结论

本文运用Hudson裂缝等效介质理论和随机介质理论建立裂缝介质中含有空间分布小孔洞的各向异性随机介质模型来描述海相碳酸盐岩储层中尺度比波长小很多的裂缝和宏观孔洞，采用伪谱法来模拟地震波在该模型中的传播波场；并建立了含碳酸盐岩储层各向异性随机介质的两层模型，模拟了地震波在该模型中的地震反射记录。为了方便认识，同时对不含孔洞的裂缝各向异性模型也做了相应模拟。

两种数值结果都表明，由于介质中定向裂缝的存在，地震波在各向异性随机模型中传播时，会发生横波分裂的现象；地震波会在空间随机分布的孔洞上发生散射，使得整个空间的波场变得复杂，当孔洞的反射很强时，这种散射波可能会在地表记录干扰其他的反射波型，实际解释时要特别加以注意。

参考文献

［1］CrapmpinS.Effectiveanisotropicelasticconstantsforwavepropagationthroughcrackedsolids.GeophysicalJournaloftheRoyalAstronomical，1984，76（1）：135～145.

［2］HudsonJA.Overallpropertiesofacrackedsolid.MathematicalProceedingoftheCambridge，1980，88：371～384.

［3］HudsonJA.Wavespeedsandattenuationofelasticwavesinmaterialcontainingcracks.GeophysicalJournaloftheRoyalAstronomical，64（1）：133～150.

［4］SchoenbergM，SayersCM.Seismicanisotropyanisotropyoffrctredrock.Geophysics，1995，60：204～211.

［5］HuangBS.Aprogramfortwo-dimensionalseismicwavepropagationbythepseudospectrummethod.Computers＆Geosciences，1992，32：289～307.

［6］董良国，黄超，刘玉柱，等.溶洞地震反射波特征数值模拟研究.石油物探，2010，49（2）：121～124.

［7］AkiK.Analysisoftheseismiccodaoflocalearthquakesasscatteredwaves.JournalofGeophysicalResearch，1969，74：615～631.

［8］BerteussenKA，ChristofferssonAES，HusebyeES，andDahleA.WavescatteringtheoryinanalysisofPwaveanomaliesatNORSARandLASA：GeophysicalJournal.RoyalAstronomicalSociety，1975，42：403～417.

［9］FrankelA，ClaytonRW.Finite-differencesimulationsofseismicscatteringimplicationsforthepropagationofshort-periodseismicwavesinthecrustandmodelsofcrustalheterogeneity.JournalofGeophysicsResearch，1986，91，6465～6489.

［10］LiuJ，BaJ，MaJW，YangHZ.Ananalysisofseismicattenuationinrandomporousmedia：ScienceChina-Physics，Mechanics＆Astronomy，2010，53：628～637.

［11］姚姚，奚先.随机介质模型正演模拟及其地震波场分析.石油物探，2002，41（1）：31～36.

［12］GazdagJ.Modelingoftheacousticwaveequationwithtransformmethods.Geophysics，1981，46（6）：854～859.

［13］KosloffDD，BaysalE.ForwardmodelingbyaFouriermethod.Geophysics，1982，47（10）：1402～1412.

［14］ForbergB.Thepseudospectralmethod：Comparisonswithfinitedifferencesfortheelasticwaveequation.Geophysics，1987，52（4）：483～501.