CN103792581B - 基于岩石物理模型的页岩气储层地应力测井预测方法 - Google Patents

Abstract

本发明涉及基于岩石物理模型的页岩气储层地应力测井预测方法，建立考虑干酪根颗粒的页岩气储层岩石物理模型，预测测井纵、横波速度；在此基础上计算储层的最大、最小水平主应力，及破裂压力，以在没有横波测井的情况下也能够对页岩气储层进行精确的应力评价；本发明的有益效果是：通过建立符合页岩气储层特征的岩石物理模型，提高了预测速度的精度；在岩石物理模型基础上求取最大、最小水平主应力、破裂压力，使得在没有实测横波测井数据的情况下，依据常规测井曲线就能够对地下应力进行预测，且预测结果精度较高。

Description

基于岩石物理模型的页岩气储层地应力测井预测方法

技术领域

本发明属于页岩气储层的地应力测井预测领域，尤其涉及基于岩石物理模型的页岩气储层地应力测井预测方法。

背景技术

对于页岩气储层，由于压裂改造直接影响着气产量，而对地下应力的预测决定着压裂改造的好坏。储层地应力的测井评价模型有多种，比较成熟的有Mohr Coulomb模型，黄荣樽模型（参见参考文献1），单轴应变模型等。通过这些模型求取最大、最小水平主应力，进而求取破裂压力，应力梯度等，可以对储层进行评价；但是这些模型的应用必须以较为精确的杨氏模量、泊松比等岩石力学参数作为输入参数，而精确的岩石力学参数，必须要有精确的纵、横波速度才能计算得到，通常，纵、横波速度可以由偶极子测井获得。但是，大部分已钻井中没有横波测井资料，这给地应力的评价造成了难度。Mullen（参见参考文献9）等使用一个综合的经验公式求算横波，并使用单轴应力公式计算最小闭合应力以对地下的应力场进行评价。经验公式求算横波受地区具体地质情况影响，很难求取准确。通过建立符合具体岩石物理特征的岩石物理模型，可以计算得到比较准确的横波。页岩气储层岩性复杂，微观孔隙类型多样，富含干酪根，针对这种非常规储层，常规的岩石物理模型不能够求算较为准确的速度，经典的Xu-White模型（参见参考文献2），将砂岩和粘土混合，使用二维的Kuster-模型（参见参考文献6）求取带孔隙的干岩石的弹性模 量，使用Gassmann（参见参考文献10）方程求取饱和流体的岩石的弹性模量，但是这个模型只考虑了粘土和砂岩，岩性不复杂，且孔隙类型也较为简单；在此基础上发展的Xu-Payne岩石物理模型（参见参考文献7）及DEM-Gassmann岩石物理模型（参见参考文献3）都使用三维孔隙谱计算含不同孔隙类型岩石的体积模量和剪切模量，得到较好的效果，但是这两个模型都只适合于比较简单的岩性，且都没有考虑干酪根，不适合于页岩气储层。Bandyopadhyay等（参见参考文献8）使用Brown-Korriga（参见参考文献4）模型，采用固体替换的方式，计算含干酪根岩石的弹性模量，但是根据岩心观察，干酪根是以一定颗粒形态分散存在于泥页岩中的，这种方式模拟的不是有机质在岩石中的真实状态。

发明内容

本发明的目的是提出基于岩石物理模型的页岩气储层地应力测井预测方法的技术方案，建立考虑干酪根颗粒的页岩气储层岩石物理模型，真实的反映在岩石中分散存在的干酪根的的弹性特征；预测测井纵、横波速度；在此基础上计算储层的最大、最小水平主应力，及破裂压力，以在没有横波测井的情况下也能够对页岩气储层进行精确的应力评价；

为了实现上述目的，本发明的技术方案是：基于岩石物理模型的页岩气储层地应力测井预测方法，建立包括干酪根颗粒的页岩气储层岩石物理模型，在常规测井曲线的基础上预测测井纵、横波速度；在此基础上计算储层的最大、最小水平主应力，及破裂压力；所建立的岩石物理模型，包括矿物，干酪根，孔隙及其中所含的流体，所述孔隙中的流体包括油、气、水，其中将干酪根等效为硬币状裂缝形态的 颗粒，使用颗粒纵横比表示其大小和形状，以微分方式计算含干酪根的岩石的弹性参数，得到在岩石中分散存在的干酪根的弹性特征；

所述方法包括以下步骤：

步骤1：通过测井复杂岩性分析程序进行测井解释，获得矿物体积、岩石孔隙度、流体饱和度、干酪根体积，进而求得矿物在岩石中的体积分数，所述矿物包括粘土、石英、方解石；使用交会分析获得每个矿物的弹性模量和密度，使用已知文献中的经验值作为孔隙中流体的弹性模量和密度，所述弹性模量指体积模量、剪切模量；

步骤2：根据步骤1中所得的矿物的体积模量和剪切模量，矿物体积分数，使用等应力和等应变平均的算术平均估算岩石基质的体积模量和剪切模量，所述岩石基质指不含孔隙和流体的固体岩石部分；在此基础上，将岩石基质与干酪根颗粒混合，即计算含干酪根颗粒的岩石的体积模量和剪切模量，将干酪根等效为硬币状裂缝形态的颗粒，使用颗粒纵横比表示干酪根颗粒的形状，使用微分等效介质方法计算含所述干酪根颗粒的岩石的体积模量和剪切模量；在所述的微分等效介质方法中，形状因子使用Berryman的三维包含物理论中的硬币状裂缝公式进行计算，所述颗粒纵横比指硬币状裂缝形状颗粒的短轴和长轴之比，在岩石物理模型中颗粒纵横比的取值范围为0～1，作为迭代计算的变量参数；

步骤3：将步骤1所得的孔隙度划分为砂岩矿物孔隙度、灰岩矿物孔隙度、粘土矿物孔隙度三种；将灰岩矿物孔隙度、砂岩矿物孔隙度等效为刚性孔隙占多数，柔性孔隙占少数的孔隙谱；将粘土矿物孔隙度等效为柔性孔隙占多数，刚性孔隙占少数的孔隙谱；将不同的孔隙度与步骤2所得到的含干酪根的岩石基质相混合，这个过程使用微分 等效介质方法计算，得到的结果为包括孔隙、干酪根、岩石基质的干岩石的体积模量、剪切模量；所述干岩石为不含流体的岩石固体部分。

步骤4：将流体与步骤3所得到的干岩石相混合，该过程使用流体替换方程计算，所得到的结果为饱和流体岩石的体积模量和剪切模量；所述的流体指混合流体，是在步骤1得到的饱和度及流体的弹性模量的基础上，使用等应力平均公式计算不同饱和度的油、气、水混合时的体积模量，混合流体的剪切模量为0；根据步骤1所得的矿物体积、孔隙度、饱和度、干酪根体积，使用算术平均计算饱含流体岩石的密度；所述流体替换方程为Gassmann方程，所述算术平均计算饱含流体岩石的密度的方法采用Wood公式；

步骤5：由步骤4所得的饱和流体岩石的体积模量、剪切模量、密度计算饱和流体岩石的纵、横波速度，即预测纵波速度和预测横波速度，对比实测纵波速度与预测纵波速度，相减求误差，当误差大于限定值时，返回步骤2，调整所述颗粒纵横比的取值范围，直到误差小于限定值或误差达到极小值，当误差小于限定值或达到极小值时，输出预测的纵、横波速度，这是一个迭代的过程；当测井有实测横波速度时，也需要对预测横波速度与实测横波速度进行对比求误差并用于判别迭代是否结束。

步骤6：使用步骤5所得到的预测纵、横波速度计算动态杨氏模量和泊松比，使用地区经验公式将动态杨氏模量和泊松比转换为静态杨氏模量和泊松比；所述的地区经验公式是张晋言和孙建孟关系式；使用测井数据求取上覆地层压力和孔隙压力，所述测井数据特指密度。

步骤7：将步骤6所得的静态杨氏模量和泊松比，上覆地层压力、 孔隙压力代入黄荣樽应力模型计算最大水平主应力、最小水平主应力，在这个过程中，需要实际工区岩心的应力测试数据对黄荣樽应力模型中的参数进行标定；根据应力参数之间的关系，计算破裂压力；综合最大水平主应力、最小水平主应力、破裂压力，进行地层应力评价。

更进一步，步骤1中所述复杂岩性分析程序，是Atlas公司的复杂岩性程序。

更进一步，步骤1中，所述的交会分析是指，在测量有纵波、横波、密度测井数据的井中，根据纵横波、密度计算体积模量、剪切模量，根据测井、录井、钻井及地质信息综合判断出井筒中不同矿物较纯的岩性段，分别将这些层段的体积模量、剪切模量、密度数据使用交会图的方法进行显示，选择数据点较为集中的位置作为纯的矿物的体积模量、剪切模量和密度数据。

本发明的有益效果是：所建立岩石物理模型可更真实的反映页岩气储层的复杂矿物以及复杂矿物造成的多孔隙类型，即以一定形态分散存在的干酪根的岩石物理特征，从而在使用该岩石物理模型进行储层速度预测的过程中提高预测速度的精度；在此基础上的应力预测使得在没有测井横波速度的情况下，根据常规测井曲线可以计算得到与测量值吻合的地下的最大、最小水平主应力，以及破裂压力，是一种可靠的页岩气储层地应力测井预测方法。

下面结合附图和实施例对本发明作一详细描述。

附图说明

图1是基于岩石物理模型的页岩气储层地应力测井预测技术流程图；

图2是所建立包含干酪根的岩石物理模型的流程示意图；

图3是建立的新岩石物理模型与经典的Xu-White模型速度预测结果对比图，从左往右，图道依次为地层，深度，录井岩性，气测甲烷与烃类，新岩石物理模型预测纵波时差与实测纵波时差，Xu-White岩石物理模型预测纵波时差与实测纵波时差，新岩石物理模型预测横波时差与实测横波时差，Xu-White岩石物理模型预测横波时差与实测横波时差；岩心图片为箭头所指位置的储层特征；图道中，CH4指甲烷，Hydrocarbon指烃类，Acpre_new指本发明预测的纵波速度，Dspre_new指本发明预测的横波速度，Acpre_XW指Xu-White模型预测的纵波速度，Dspre_XW指Xu-White模型预测的横波速度，AC指测量纵波时差，DTS指测量横波时差；所述时差为速度的导数；

图4为建立的新岩石物理模型与经典的Xu-White岩石物理模型预测速度与实测速度之间的相对误差的对比图，a),b)图为新模型预测速度的纵、横波误差；c),d)为Xu-White模型预测速度的纵、横波误差；

图5是页岩气储层应力计算结果，从左往右，图道依次为地层，伽马与井径曲线，深度，录井岩性，气测甲烷与烃类,静态泊松比，静态杨氏模量，实测与预测上覆地层压力、孔隙压力，实测与预测最小水平主应力，实测与预测最大水平主应力，破裂压力；其中，实测的最大、最小水平主应力，上覆地层压力为杆状离散点，预测最大、最小水平主应力，上覆地层压力为曲线。

具体实施方式

基于岩石物理模型的页岩气储层地应力测井预测方法，建立考虑干酪根颗粒的页岩气储层岩石物理模型，在常规测井曲线的基础上预 测测井纵、横波速度；在此基础上计算储层的最大、最小水平主应力，及破裂压力，以在没有横波测井的情况下根据常规测井曲线对页岩气储层进行精确的应力评价；所建立的岩石物理模型，包括矿物，干酪根，干酪根，孔隙及其中所含的流体，所述孔隙中的流体包括油、气、水，其中将干酪根等效为硬币状裂缝形态的颗粒，使用颗粒纵横比表示其大小和形状，以微分方式计算含干酪根的岩石的弹性参数，得到在岩石中分散存在的干酪根的弹性特征。

以下实施例结合中国南方某页岩气盆地侏罗系下统页岩气储层为例来阐述本发明的具体实施方式：

如图1所示，为基于岩石物理模型的页岩气储层地应力测井预测方法的流程图，首先通过测井解释得到地层中的矿物体积、孔隙度、饱和度、干酪根体积；通过交会分析得到主要矿物的弹性模量和密度，通过文献调研得到流体的弹性模量和密度，此处弹性参数特指体积模量、剪切模量；将测井解释结果引入新建立的页岩气储层岩石物理模型，通过调整岩石物理模型中的孔隙类型分布谱参数，迭代计算纵波速度与横波速度。使用预测纵、横波速度计算动态杨氏模量和泊松比，并根据地区实测数据总结的经验关系式将动态的杨氏模量和泊松比转换为静态杨氏模量和泊松比；由测井数据计算上覆地层压力和孔隙压力，此处测井数据特指密度和声波时差；将静态的杨氏模量、泊松比，上覆地层压力、孔隙压力代入黄荣樽应力模型中计算最大水平主应力和最小水平主应力，在这个过程中需要根据实际工区岩心应力测试数据标定黄荣樽应力模型中的参数；在此基础上计算破裂压力；在最大水平主应力，最小水平主应力、破裂压力的基础上就可以对地层应力进行评价。

对图1所述流程共分7个步骤进行详细说明，包括两部分主要内容，第一部分为岩石物理模型的构建和基于岩石物理模型的速度预测，包括下述步骤1、2、3、4、5，图2为岩石物理模型的流程图；第二部分为应力的预测，包括下述步骤6、7。以下结合公式即实例进行详细说明。

步骤1：通过测井复杂岩性分析程序进行测井解释，获得矿物体积、岩石孔隙度、流体饱和度、干酪根体积，进而求得矿物在岩石中的体积分数，所述矿物包括粘土、石英、方解石；使用交会分析获得每个矿物的弹性模量和密度，使用已知文献中的经验值作为孔隙中流体的弹性模量和密度，所述弹性模量指体积模量、剪切模量。

所述复杂岩性分析程序，是Atlas公司的复杂岩性程序；通过该程序，可以根据常规测井曲线计算较为精细的泥质含量，矿物体积分数，孔隙度，流体饱和度等参数。其具体公式记载在雍世和、张超谟主编的于2002年由中国石油大学出版社出版的《测井数据处理与综合解释》书的313-319页，参见参考文献12。

所述的交会分析是指，在测量有纵波、横波、密度测井数据的井中，根据纵横波、密度计算体积模量、剪切模量，根据测井、录井、钻井及地质信息综合判断出井筒中不同矿物较纯的岩性段，分别将这些层段的体积模量、剪切模量、密度数据使用交会图的方法进行显示，选择数据点较为集中的位置作为纯的矿物的体积模量、剪切模量和密度数据。

本实施例中，得到的矿物、流体的弹性模量和密度如表1所示。

表1本发明中矿物及流体的性质参数

表1中部分数据参考Mavko等人于1999年由Cambridge University出版的The rock physics handbook:tools for seismic analysis in porous media书中的306-309页。

步骤2：根据步骤1中所得的矿物的体积模量和剪切模量，矿物体积分数，使用等应力和等应变平均的算术平均估算岩石基质的体积模量和剪切模量，所述岩石基质指不含孔隙和流体的固体岩石部分；在此基础上，将岩石基质与干酪根颗粒混合，即计算含干酪根颗粒的岩石的体积模量和剪切模量，将干酪根等效为硬币状裂缝形态的颗粒，使用颗粒纵横比表示干酪根颗粒的形状，使用微分等效介质方法计算含这种干酪根颗粒的岩石的体积模量和剪切模量；在所述的微分等效介质方法中，形状因子使用Berryman的三维包含物理论中的硬币状裂缝公式进行计算，而颗粒纵横比指硬币状裂缝的短轴和长轴之比，在岩石物理模型中颗粒纵横比的取值范围为0～1，作为迭代计算的变量参数之一；

所述等应变平均指Voigt平均，由Voigt在1928年提出，使用公式（1）表示；等应力平均指Reuss平均，由Reuss在1929年提出，使用公式（3）表示；两者的算术平均指Hill平均，由Hill在1952年提出，使用公式（1）表示，公式（1），（2），（3）分别在Mavko等于1999年由Cambridge University出版的The rock physics handbook:tools for seismic analysis in porous media书中的110,110,114页。

$M_{\mathrm{VRH}}=\frac{M_V+M_R}{2}---\left(1\right)$

$M_V=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{V}_i{M}_i---\left(2\right)$

$\frac{1}{M_R}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{V_i}{M_i}---\left(3\right)$

式中，M_VRH为要求的弹性模量，MV为Voigt模型所求的弹性模量,MR为Reuss模型所求的弹性模量,Vi和Mi为第i个组成成分的体积分数和弹性模量。此时，这里使用的粘土、石英、方解石的体积含量分别为：

$V_{\mathrm{shale}}^{\′}=\frac{V_{\mathrm{shale}}}{1-{\mathrm{\φ}}_{\mathrm{total}}-V_{\mathrm{kerogen}}},V_{\mathrm{quartz}}^{\′}=\frac{V_{\mathrm{quartz}}}{1-{\mathrm{\φ}}_{\mathrm{total}}-V_{\mathrm{kerogen}}},V_{\mathrm{calcite}}^{\′}=\frac{V_{\mathrm{calcite}}}{1-{\mathrm{\φ}}_{\mathrm{total}}-V_{\mathrm{kerogen}}}.$

其中φtotal为总孔隙度，Vshale,Vquartz,Vcalcite,Vkerogen分别为粘土，石英，方解石，干酪根占岩石总体积的百分比，V’shale,V’quartz,V’calcite分别为粘土，石英，方解石占岩石基质的百分比。

所述微分等效介质方法具体公式为：

$(1-y)\frac{d}{\mathrm{dy}}\left[K^{*}\left(y\right)\right]=(K_{\mathrm{kerogen}}-K^{*})P^{(*2)}\left(y\right)---\left(4\right)$

$(1-y)\frac{d}{\mathrm{dy}}\left[{\mathrm{\μ}}^{*}\left(y\right)\right]=({\mathrm{\μ}}_{\mathrm{kerogen}}-{\mathrm{\μ}}^{*})Q^{(*2)}\left(y\right)---\left(5\right)$

式中，K*、μ*分别代表所求的含干酪根颗粒的岩石的体积模量和剪切模量；Kkerogen、μkerogen分别代表干酪根的体积模量和剪切模量；y代表干酪根所占的体积分数，P、Q为Berryman三维等效包含物理论中的孔隙类型形状因子，具体公式见Berryman于1995年发表在America Geophysical Union上205-228的《Mixture theories for rock properties,in A handbook of Physical Constants》一文中,见参考文献（5），此处P、Q代表干酪根颗粒的形状因子，代表硬币状裂缝形态，是关于颗粒纵横比的函数。

步骤3：将步骤1所得的孔隙度划分为砂岩矿物孔隙度、灰岩矿物孔隙度、粘土矿物孔隙度三种；将灰岩矿物孔隙度、砂岩矿物孔隙度等效为刚性孔隙占多数，柔性孔隙占少数的孔隙谱；将粘土矿物孔隙度等效为柔性孔隙占多数，刚性孔隙占少数的孔隙谱；将不同的孔隙度与步骤2所得到的含干酪根的岩石基质相混合，这个过程使用微分等效介质方法计算，得到的结果为包括孔隙、干酪根、岩石基质的干岩石的体积模量、剪切模量；所述干岩石为不含流体的岩石固体部分。

本步骤中，岩石中的孔隙度划分为砂岩矿物孔隙度、灰岩矿物孔隙度、粘土矿物孔隙度，其划分方法用公式（6）表示

式中，Vshale,Vquartz,Vcalcite,Vkerogen分别为粘土，石英，方解石，干酪根的体积分数，单位为%；φtotal为孔隙度；φquartz，φcalcite，φclay分别为石英，方解石，粘土矿物孔隙度。

将不同的孔隙度与含干酪根的岩石基质相混合，过程使用微分等效介质方法计算。在计算的过程中，将粘土孔隙度、砂岩孔隙度、灰岩孔隙度依次进行计算，计算通用的公式为：

$(1-y)\frac{d}{\mathrm{dy}}\left[{K_{\mathrm{dry}}}^{*}\left(y\right)\right]=(K_i-{K_{\mathrm{dry}}}^{*})P^{(*2)}\left(y\right)---\left(7\right)$

$(1-y)\frac{d}{\mathrm{dy}}\left[{{\mathrm{\μ}}_{\mathrm{dry}}}^{*}\left(y\right)\right]=({\mathrm{\μ}}_i-{{\mathrm{\μ}}_{\mathrm{dry}}}^{*})Q^{(*2)}\left(y\right)---\left(7\right)$

式中，Kdry*、μdry*分别代表所求的干岩石的体积模量和剪切 模量；Ki、μi分别代表孔隙的体积模量和剪切模量；y代表包含物所占的百分比，P、Q为每种孔隙的形状因子。

步骤4：将流体与步骤3所得到的干岩石相混合，该过程使用流体替换方程计算，所得到的结果为饱和流体岩石的体积模量和剪切模量；所述的流体指混合流体，是在步骤1得到的饱和度及流体的弹性模量的基础上，使用等应力平均公式计算不同饱和度的油、气、水混合时的体积模量，混合流体的剪切模量为0；根据步骤1所得的矿物体积、孔隙度、饱和度、干酪根体积，使用算术平均计算饱含流体岩石的密度；所述算术平均计算饱含流体岩石的密度的方法即Wood公式；

所述等应力平均公式见公式（3）；

所述流体替换方程为Gassmann方程，由Gassmann在1951年提出，具体公式在Mavko等人于1999年由Cambridge University出版的The rock physics handbook:tools for seismic analysis in porous media书中的168页，公式为

$\frac{K_{\mathrm{sat}}}{K_0-K_{\mathrm{sat}}}=\frac{K_{\mathrm{dry}}}{K_0-K_{\mathrm{dry}}}+\frac{K_{\mathrm{fl}}}{{\mathrm{\φ}}_{\mathrm{total}}(K_0-K_{\mathrm{fl}})},{\mathrm{\μ}}_{\mathrm{sat}}={\mathrm{\μ}}_{\mathrm{dry}}---\left(9\right)$

其中，Kdry,Ksat,K0,Kfl分别是干岩石，饱和流体岩石，基质岩石，流体的体积模量,φtotal为总孔隙度,μsat,μdry为所求的饱和流体岩石，干岩石的剪切模量。

所述计算饱和流体岩石密度的Wood公式为

$\mathrm{\ρ}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{f}_i{\mathrm{\ρ}}_i---\left(10\right)$

式中，fi指油、气、水的饱和度、基质矿物的体积比，ρi为油、气、水等流体、基质矿物的密度，ρ为所求饱和流体岩石的密度。

步骤5：由步骤4所得的饱和流体岩石的体积模量、剪切模量、密 度计算饱和流体岩石的纵、横波速度，即预测纵波速度和预测横波速度，对比实测纵波速度与预测纵波速度，相减求误差，当误差大于限定值时，返回步骤2，调整所述颗粒纵横比的取值范围，其调整范围为步骤2中所述的0～1，直到误差小于限定值或误差达到极小值，当误差小于限定值或达到极小值时，输出预测的纵、横波速度，这是一个迭代的过程；当测井有实测横波速度时，也需要对预测横波速度与实测横波速度进行对比求误差并用于判别迭代是否结束。

本实施例中，所述误差的限定值为0.0001。

由体积模量、剪切模量、密度计算饱和流体岩石的纵、横波速度，公式为

$V_p=\sqrt{\frac{K_{\mathrm{sat}}+4/3{\mathrm{\μ}}_{\mathrm{sat}}}{\mathrm{\ρ}}},V_s=\sqrt{\frac{{\mathrm{\μ}}_{\mathrm{sat}}}{\mathrm{\ρ}}}---\left(11\right)$

其中，Vp，Vs为计算的纵波速度和横波速度，Ksat与μsat为饱和流体岩石的体积模量、剪切模量，ρ为饱和流体岩石的密度。

所述误差为相对误差，使用Error=(Vp实测-Vp预测)/Vp实测计算。 

图3为利用上述5个步骤得到的速度预测结果图，图中对比了工业界常用的Xu-White岩石物理模型速度预测结果和本发明的岩石物理模型所得速度预测结果，图中从左往右的图道依次为地层，深度，录井岩性，气测甲烷与烃类，新岩石物理模型预测纵波时差与实测纵波时差，Xu-White岩石物理模型预测纵波时差与实测纵波时差，新岩石物理模型预测横波时差与实测横波时差，Xu-White岩石物理模型预测横波时差与实测横波时差；岩心图片为箭头所指位置的储层特征；图道中，Formation指地层，Depth指深度，Lithology指岩性，CH4指甲烷，Hydrocarbon指烃类，Acpre_new指本发明预测的纵波速度， Dspre_new指本发明预测的横波速度，Acpre_XW指Xu-White模型预测的纵波速度，Dspre_XW指Xu-White模型预测的横波速度，AC指测量纵波时差，DTS指测量横波时差。从图中可以看出，新的岩石物理模型在页岩气储层段的预测纵、横波速度与实测纵、横波速度之间的距离相比较于Xu-White模型预测的纵、横波速度和实测纵、横波速度之间的距离更加靠近，表明新的岩石物理模型预测速度更加精确。图4为新岩石物理模型的预测速度的相对误差与经典的Xu-White岩石物理模型预测速度的相对误差的直方图对比，a,b为新模型的预测纵、横波相对误差直方图，c,d为Xu-White模型的预测纵、横波相对误差直方图，新模型的纵、横波相对误差明显小于Xu-White模型的纵、横波误差，证明了这种岩石物理模型求取速度的有效性。

步骤6：使用步骤5所得到的预测纵、横波速度计算动态杨氏模量和泊松比，使用地区经验公式将动态杨氏模量和泊松比转换为静态杨氏模量和泊松比；

所述动态杨氏模量和泊松比的计算公式为：

$E_d=\frac{\mathrm{\ρ}{V_s}^2({3V}_p^2-4{V_s}^2)}{(V_p^2-{V_s}^2)}---\left(12\right)$

$v_d=\frac{V_p^2-{{2V}_s}^2}{2(V_p^2-{V_s}^2)}---\left(13\right)$

式中，Ed为动态杨氏模量，νd为动态泊松比，Vp为纵波速度，Vs为横波速度，ρ为饱和流体岩石的密度。

所述的地区经验公式是张晋言和孙建孟关系式,记载在张晋言、孙建孟发表在《测井技术》2012,36(2):146-153页的《利用测井资料评价泥页岩油气“五性”指标》文章中，见参考文献（1）；

E_s＝0.37E_d+0.655    (14)

ν_s＝0.44ν_d+0.446    (15) 

式中Es为静态弹性模量;νs为静态泊松比。

使用测井数据求取上覆地层压力和孔隙压力，公式为（16）（17），所述测井数据特指声波时差和密度。

${\mathrm{\σ}}_v={\∫}_0^z{\mathrm{\ρ}}_z\mathrm{gdz}---\left(16\right)$

$p_p={\∫}_0^z{\mathrm{\ρ}}_f\mathrm{gdz}---\left(17\right)$

σv为上覆地层压力,ρz，ρf为密度、流体密度曲线，g为重力加速度。

步骤7：将步骤6所得的静态杨氏模量和泊松比，上覆地层压力、孔隙压力代入黄荣樽应力模型计算最大水平主应力、最小水平主应力，在这个过程中，需要实际工区岩心的应力测试数据对黄荣樽应力模型中的参数进行标定；根据应力参数之间的关系，计算破裂压力；综合最大水平主应力、最小水平主应力、破裂压力，进行地层应力评价。

所述黄荣樽应力模型公式为（18）（19），指黄荣樽1984年提出的计算应力模型，本处所使用的公式为张晋言和孙建孟在2012年发表于测井技术上的文章《利用测井资料评价泥页岩油气“五性”指标》，见参考文献（1）：

${\mathrm{\σ}}_h=\frac{v_s}{1.0-v_s}({\mathrm{\σ}}_v-{\mathrm{\αp}}_p)+B_1({\mathrm{\σ}}_v-{\mathrm{\αp}}_p)+{\mathrm{\αp}}_p---\left(18\right)$

${\mathrm{\σ}}_H=\frac{v_s}{1.0-v_s}({\mathrm{\σ}}_v-{\mathrm{\αp}}_p)+B_2({\mathrm{\σ}}_v-{\mathrm{\αp}}_p)+{\mathrm{\αp}}_p---\left(19\right)$

式中，B1,B2为地质构造应力系数，由工区目标实际岩心测量的应力数据根据公式反算标定得到；α为Biot系数；Pp为孔隙压力，σH为最大水平主应力；σh为最小水平主应力；vs为静态泊松比。

根据最大、最小水平主应力计算破裂压力公式为（20），参考于朱玉林于2007年的在中国石油大学（华东）的学位论文《测井资料在地 应力研究中的应用》一文中，见参考文献（11），所述破裂压力指井内钻井液柱压力过大时，井壁会出现张性应力，导致井壁岩石发生拉张破裂，出现井漏事故。这种拉张破裂时的井内液柱压力称为地层破裂压力。

p_F＝3σ_h-D_H-αp_p+F_co/4.0    (20) 

式中，PF为破裂压力，σh为最小水平主应力，σH为最大水平主应力，Pp为孔隙压力，α为Biot系数，Fco为内聚力。 

图5为页岩气储层应力计算结果，图中给出了步骤6所得静态杨氏模量和泊松比，以及步骤7中所得的最大、最小水平主应力、破裂压力。图中，从左往右，图道依次为地层，伽马与井径曲线，深度，录井岩性，气测甲烷与烃类,静态泊松比，静态杨氏模量，实测与预测上覆地层压力、孔隙压力，实测与预测最小水平主应力，实测与预测最大水平主应力，破裂压力；图道中，Formation指地层，GR指伽马测井曲线，CAL指井径曲线，Depth指深度，Lithology指岩性，CH4指甲烷，Hydrocarbon指烃类，图中，上覆地层压力，最大与最小水平主应力道中有一个实际测量值，这个测量的数据点，证明了预测结果与实测结果是比较吻合的，说明这种基于岩石物理模型的页岩气储层地应力测井预测技术是可靠的。

步骤1中的矿物类型包括粘土、方解石、石英，考虑了页岩气储层中最主要的三种岩性，也是根据测井数据可以比较准确求得的。

步骤2中的干酪根颗粒等效方式，使用Berryman的三维包含物形态中的裂缝形态来描述与计算，因为干酪根为柔性物质，不同于固体矿物，又存在剪切模量，不同于流体，使用这种表达方式，与镜下观察的有机质存在方式更为符合。

步骤3中的孔隙度表述方式，使用Berryman的三维包含物形态中 球形孔隙，针形孔隙与裂缝孔隙来表述与计算，其中方解石、石英、白云石孔隙为以球形孔隙和针形孔隙等刚性孔隙占主要部分，裂缝型孔隙等柔性孔隙占少数部分的孔隙谱表述，粘土孔隙为以裂缝型孔隙占多数部分，球形和针形孔隙占少数部分的孔隙谱表述，具体表述公式见Berryman于1995年发表在America Geophysical Union上205-228的《Mixture theories for rock properties,in A handbook of Physical Constants》一文中。

以上三点使得新建的岩石物理模型优于经典的Xu-White模型，即包括的矿物类型较多；较为合理的表述并计算了岩石中有机质对弹性性质的影响；分开使用较为复杂的孔隙类型表述了岩石中石英，方解石、粘土成分中的孔隙类型，与岩石中的实际情况更加符合，使得速度预测的精度相比较于Xu-White模型有较大的改善。

步骤6、7中，将符合页岩气储层岩石物理特征的岩石物理模型预测的弹性参数引入应力计算模型，得到较为精确的评价页岩储层地下应力的数据，即将岩石物理模型与应力计算模型相结合，在没有测量横波的常规测井曲线的基础上，根据常规测井曲线也可以进行较为精确的地下应力场评价。

参考文献： 

1.张晋言,孙建孟.利用测井资料评价泥页岩油气“五性”指标.测井技术,2012,36(2):146-153；

2.Xu,S.,R.E.White.A new velocity model for clay-sand mixtures:Geophysical Prospecting,1995,43:91-118；

3.Sun,S.Z,H.Wang,Z.Liu,et al.The theory and  application of DEM-Gassmann rock physics model for complex carbonate reservoirs.The leading edge,2012,31(2):152-158. 

4.Brown,R.,and J.Korringa.On the dependence of the elastic properties of a porous rock on the compressibility of the pore fluid.Geophysics,1975,40,608-616；

5.Berryman,J.G..Mixture theories for rock properties,in A handbook of Physical Constants.America Geophysical Union,1995,205-228；

6.Kuster,G.T.,M.N.Velocity and attenuation of seismic waves in two-phase media.Geophysics,1974,39:587-618；

7.Xu,S.,M.A.Payne.Modeling elastic properties in carbonate rocks:The leading edge,2009,28(1):66-74；

8.Bandyopadhyay,K.,R.Sain,E.Liu,and et al.,Rock property inversion in organic-rich shale:uncerties,ambiguities,and pitfalls:SEG annual meeting,2012；

9.Mullen,M.,R.Roundtree.,and B.Barree.A composite determination of mechanical rock properties for stimulation design(what to do when you don’t have a sonic log).Society of Petroleum Engineers.2007:SPE108139；

10.Mavko,G.,T.Mukerji.,and J.Dvorkin.The rock physics handbook:tools for seismic analysis in porous media.Cambridge University press.1999；

11.朱玉林.测井资料在地应力研究中的应用.东营:中国石油大 学,2007；

12.雍世和，张超谟.测井数据处理与综合解释.东营.中国石油大学出版社，2002：313-319。

Claims (3)

1.基于岩石物理模型的页岩气储层地应力测井预测方法，其特征在于，建立包括干酪根颗粒的页岩气储层岩石物理模型，在常规测井曲线的基础上预测测井纵、横波速度；在此基础上计算储层的最大、最小水平主应力，及破裂压力；所建立的岩石物理模型，包括矿物，干酪根，孔隙及其中所含的流体，所述孔隙中的流体包括油、气、水，其中将干酪根等效为硬币状裂缝形态的颗粒，使用颗粒纵横比表示其大小和形状，以微分方式计算含干酪根的岩石的弹性参数，得到在岩石中分散存在的干酪根的弹性特征；

所述方法包括以下步骤：

步骤1：通过测井复杂岩性分析程序进行测井解释，获得矿物体积、岩石孔隙度、流体饱和度、干酪根体积，进而求得矿物在岩石中的体积分数，所述矿物包括粘土、石英、方解石；使用交会分析获得每个矿物的弹性模量和密度，使用已知文献中的经验值作为孔隙中流体的弹性模量和密度，所述弹性模量指体积模量、剪切模量；

步骤2：根据步骤1中所得的矿物的体积模量和剪切模量，矿物体积分数，使用等应力和等应变平均的算术平均估算岩石基质的体积模量和剪切模量，所述岩石基质指不含孔隙和流体的固体岩石部分；在此基础上，将岩石基质与干酪根颗粒混合，即计算含干酪根颗粒的岩石的体积模量和剪切模量，将干酪根等效为硬币状裂缝形态的颗粒，使用颗粒纵横比表示干酪根颗粒的形状，使用微分等效介质方法计算含所述干酪根颗粒的岩石的体积模量和剪切模量；在所述的微分等效介质方法中，形状因子使用Berryman的三维包含物理论中的硬币状裂缝公式进行计算，所述颗粒纵横比指硬币状裂缝形状颗粒的短轴和长轴之比，在岩石物理模型中颗粒纵横比的取值范围为0～1，作为迭代计算的变量参数；

步骤3：将步骤1所得的孔隙度划分为砂岩矿物孔隙度、灰岩矿物孔隙度、粘土矿物孔隙度三种；将灰岩矿物孔隙度、砂岩矿物孔隙度等效为刚性孔隙占多数，柔性孔隙占少数的孔隙谱；将粘土矿物孔隙度等效为柔性孔隙占多数，刚性孔隙占少数的孔隙谱；将不同的孔隙度与步骤2所得到的含干酪根的岩石基质相混合，这个过程使用微分等效介质方法计算，得到的结果为包括孔隙、干酪根、岩石基质的干岩石的体积模量、剪切模量；所述干岩石为不含流体的岩石固体部分；

步骤4：将流体与步骤3所得到的干岩石相混合，该过程使用流体替换方程计算，所得到的结果为饱和流体岩石的体积模量和剪切模量；所述的流体指混合流体，是在步骤1得到的饱和度及流体的弹性模量的基础上，使用等应力平均公式计算不同饱和度的油、气、水混合时的体积模量，混合流体的剪切模量为0；根据步骤1所得的矿物体积、孔隙度、饱和度、干酪根体积，使用算术平均计算饱含流体岩石的密度；所述流体替换方程为Gassmann方程，所述算术平均计算饱含流体岩石的密度的方法采用Wood公式；

步骤5：由步骤4所得的饱和流体岩石的体积模量、剪切模量、密度计算饱和流体岩石的纵、横波速度，即预测纵波速度和预测横波速度，对比实测纵波速度与预测纵波速度，相减求误差，当误差大于限定值时，返回步骤2，调整所述颗粒纵横比的取值范围，直到误差小于限定值或误差达到极小值，当误差小于限定值或达到极小值时，输出预测的纵、横波速度，这是一个迭代的过程；当测井有实测横波速度时，也需要对预测横波速度与实测横波速度进行对比求误差并用于判别迭代是否结束；

步骤6：使用步骤5所得到的预测纵、横波速度计算动态杨氏模量和泊松比，使用地区经验公式将动态杨氏模量和泊松比转换为静态杨氏模量和泊松比，所述动态杨氏模量和泊松比的计算公式为：

$E_d=\frac{\mathrm{\ρ}{V}_s^2(3V_p^2-4V_s^2)}{(V_p^2-V_s^2)},$

$v_d=\frac{V_p^2-2V_s^2}{2(V_p^2-V_s^2)},$

式中，Ed为动态杨氏模量，νd为动态泊松比，Vp为纵波速度，Vs为横波速度，ρ为饱和流体岩石的密度；所述的地区经验公式是张晋言和孙建孟关系式：

E_s＝0.37E_d+0.655，

ν_s＝0.44ν_d+0.446，

式中Es为静态杨氏模量；νs为静态泊松比；

使用测井数据求取上覆地层压力和孔隙压力，所述测井数据特指密度；

步骤7：将步骤6所得的静态杨氏模量和泊松比，上覆地层压力、孔隙压力代入黄荣樽应力模型计算最大水平主应力、最小水平主应力，在这个过程中，需要实际工区岩心的应力测试数据对黄荣樽应力模型中的参数进行标定；根据应力参数之间的关系，计算破裂压力；综合最大水平主应力、最小水平主应力、破裂压力，进行地层应力评价；所述黄荣樽应力模型公式为：

${\mathrm{\σ}}_h=\frac{v_s}{1.0-v_s}({\mathrm{\σ}}_v-{\mathrm{\αp}}_p)+B_1({\mathrm{\σ}}_v-{\mathrm{\αp}}_p)+{\mathrm{\αp}}_p,$

${\mathrm{\σ}}_H=\frac{v_s}{1.0-v_s}({\mathrm{\σ}}_v-{\mathrm{\αp}}_p)+B_2({\mathrm{\σ}}_v-{\mathrm{\αp}}_p)+{\mathrm{\αp}}_p,$

式中，B1,B2为地质构造应力系数，由工区目标实际岩心测量的应力数据根据公式反算标定得到；α为Biot系数；Pp为孔隙压力，σH为最大水平主应力；σh为最小水平主应力；vs为静态泊松比。

2.根据权利要求1所述的基于岩石物理模型的页岩气储层地应力测井预测方法，其特征在于，步骤1中所述复杂岩性分析程序，是Atlas公司的复杂岩性程序。

3.根据权利要求1所述的基于岩石物理模型的页岩气储层地应力测井预测方法，其特征在于，步骤1中，所述的交会分析是指，在测量有纵波、横波、密度测井数据的井中，根据纵横波、密度计算体积模量、剪切模量，根据测井、录井、钻井及地质信息综合判断出井筒中不同矿物较纯的岩性段，分别将这些层段的体积模量、剪切模量、密度数据使用交会图的方法进行显示，选择数据点较为集中的位置作为纯的矿物的体积模量、剪切模量和密度数据。