作为一种新型天然气资源，页岩气约占非常规天然气资源量的50%，是非常重要的能源.由于页岩地层非均质性强，储层致密，孔隙及流体状况复杂等(董大忠等，2011；张金川等，2008；于炳松，2012；董宁等，2014)，其储层描述和流体识别难度很大.岩石物理模型可将地震数据(纵、横波速度等)与储层参数(孔隙度，流体饱和度，泥质含量等)联系起来，正确适用的岩石物理模型能够为地震正反演提供必要的数据资料，在用地震资料进行反演和解释中有着非常重要的作用(马淑芳等，2010；徐胜峰等，2009；白俊雨等，2012；印兴耀等，2015).对于复杂页岩地层，目前尚无成熟的岩石物理理论.因此，研究适用于页岩的地震岩石物理模型具有重要意义.

国内外一些学者主要针对页岩中的有机质的添加对页岩地层的岩石物理理论进行了研究.目前常用的有机质添加方法有两类，一是将其作为矿物成分，与其他矿物共同混合添加(Bayuk et al.，2008；Wu et al.，2012)，但这种方式不符合电镜照片观察到的有机质在背景介质中填充分布的分布规律；二是将其作为孔隙填充物通过固体替换添加(Zhu et al.，2012；Bandyopadhyay et al.，2012)，但该方式未考虑有机质中存在孔隙及流体的影响.另外，对于孔隙中流体添加，目前一般通过流体替换(Wu et al.，2012；Zhu et al.，2012)进行，未考虑页岩孔隙度较低、存在部分连通性差的微孔隙，流体无法在其中自由流动，不符合常用流体替换模型假设条件的情况.

页岩气藏勘探开发潜力巨大，因其基质具有超低渗透性，裂缝网络的形成成为获得工业性气流的关键，因此开发过程中需要进行大规模水力压裂，而地应力的预测是指导页岩压裂施工的重要基础，其预测的精度直接影响压裂的质量.目前常用的地应力预测方法主要以均匀各向同性介质为基础(Sayers，2010；Goodway et al.，2010；Perez et al.，2011)，而页岩储层具有较强的非均质性和各向异性特征，因此，应力预测存在一定误差.Thiercelin和Plumb(1994)、Suárez-Rivera等(2009)、邓金根等(2013)分别研究了线弹性的横向各向同性地层地应力的计算方法，但其所使用弹性参数较多，无法利用 常规测井等资料直接确定，制约了该方法思路的应用.

本文依据对页岩地层特征的分析，重点针对目前有机质、孔隙流体添加中存在的问题及岩石的各向异性特征，给出了页岩岩石物理等效模型的建立流程，根据此模型，可以实现页岩地层的横波速度的估算及最小水平地应力的求取.在以某页岩气工区井A为例验证了模型构建方法的适用性后，利用该模型求取井A页岩地层最小水平地应力并与已有的各向同性模型估算结果进行比较，结果表明利用本文模型计算的最小水平地应力更加精确，具有较高的实用价值.2 基本原理

本文从页岩岩石物理等效模型的构建出发，提出一种水平地应力有效预测的方法，为页岩地层压裂区域优选提供了一定的理论依据.其基本原理主要分为页岩各向异性岩石物理模型的构建和页岩水平地应力预测两部分.2.1 页岩各向异性岩石物理模型的构建2.1.1 页岩特征分析

有机质是页岩气聚集成藏最重要的控制因素之一，由电镜照片(图 1)观察可知，泥质为页岩中的主要基质矿物，有机质主要在其中填充分布，且其中含有较多微孔隙.依据此分布特点，本文提出将有机质作为泥质背景的包含物，计算等效介质弹性模量，然后向其中添加孔隙、流体的建模思路.

图1

Fig. 1

图 1页岩的扫描电镜照片Fig. 1Scanning electron microphotographs of shale

页岩储层为低孔低渗致密储层，其孔隙度一般小于10%，页岩尤其是其泥质组分中分布有较多不连通的微孔隙，且其中多填充不可动的束缚水和吸附气，而可自由移动的自由水和游离气主要赋存于连通性较好的基质孔隙中(姜文利等，2010；丁文龙等，2012).因此本文中根据页岩储层的特殊性将孔隙分为含有不可动流体的不连通微孔隙和含有可动流体的有一定连通性的微孔隙两类，并假设含有不可动流体的不连通微孔隙主要分布于泥质中.对不连通微孔隙连同其中的不可动流体以包含物的形式向泥质中添加，对可动流体则通过先添加干燥孔隙后进行流体替换的方式添加到基质中.2.1.2 页岩各向异性岩石物理等效模型的构建流程

依据上述页岩特征的分析，可进行页岩各向异性岩石物理等效模型的构建，构建流程主要分为“纯”页岩的获取，基质(VTI介质)等效弹性模量的计算，“干”骨架的弹性模量的求取以及饱和岩石弹性参数的计算四部分.其流程图如图 2所示.

图2

Fig. 2

图 2 页岩岩石物理等效模型构建流程图 Fig. 2 The flow chart of shale rock physics effective model

(1)“纯”页岩的获取：本文中假设含有稀疏分布不连通微孔隙的泥质为“纯”页岩.由于不连通微孔隙微小且多稀疏分布在泥质中，孔隙纵横比较低，其中多充填不可动流体，因此可假设该类孔隙呈针状，连同其中的不可动流体，使用Kuster-Toksöz模型(Kuster and ToksöKuster and Toksöz，1974)加入到背景介质—泥质中，得到“纯”页岩.

应用长波一阶散射理论，Kuster和Toksöz(1974)考虑包含物形状、含量及弹性性质的影响，得 到了应用广泛的Kuster-Toksöz模型，其表达式如下：

其中，分别是饱和岩石、岩石基质和第i种包含物的体积模量和剪切模量，xi为第i种包含物体积分数，P^mi、Q^mi表示在背景基质m中加入包含物材料i后的效果，对于针状包含物，.

(2)基质(VTI介质)等效弹性模量的计算：将 有机质作为椭球形包含物，利用自相容模型(Berryman，1980，1995)加入到“纯”页岩中，得到混合矿物1，计算该等效介质的弹性模量；使用V-R-H平均(Mavko et al.，2009)混合石英和方解石，得到混合矿物2；使用Backus平均(Backus，1962)将两等效介质进行混合，得到基质(VTI介质)的等效弹性模量.

①自相容模型

该模型既考虑了孔隙形状的影响，又考虑了彼此靠近的包含物的相互作用，因此能够适用于包含物含量较高的岩石，可用来添加含量较高的有机质. Berryman(1980，1995)给出了自相容模型的一般形式：

②Voigt-Reuss-Hill平均

Voigt-Reuss-Hill平均可计算各向同性、完全弹性介质的等效弹性模量，由于石英和方解石均为刚性矿物，弹性性质较为接近，可用该模型进行混合，其表达式为

其中，和M_i分别为第i种组分的体积分数和弹性模量，M_VRH为混合矿物的弹性模量.

③Backus平均

层状分布的页岩地层可等效为VTI介质，其基质等效模量可用Backus平均求取，依据Backus(1962)的理论，长波极限下，多层横向各向同性介质材料组成的层状介质等效刚度张量可表示为

其中，

式中C_ij为各组分弹性刚度张量，C_ij^*为求得的VTI介质的弹性刚度张量，运算符号“〈·〉”表示对各组分加权平均.

(3)“干”骨架的弹性模量的求取：利用各向异性SCA-DEM模型(Hornby et al.，1994)向基质中加入有一定连通性的干燥微孔隙，计算“干”骨架的弹性模量，其中假设该孔隙为单一孔隙纵横比的理想椭球体.

由于基质为VTI介质，因此常用添加孔隙的各向同性模型不再适用.Hornby等(1994)将各向同性的SCA模型推广到各向异性介质，得

其中，C ⁿ、v_n分别为第n种组分的弹性刚度张量及体积分数，C ^SCA为利用该模型得到的等效介质弹性刚度张量，K ^{* n}为联系复合材料平均应力与内部各相平均应变的张量，根据本征应变关系可求得K={C[I+Ĝ(C'－C)]}^－1，张量Ĝ是与包含物纵横比有关的参数，其表示形式及计算可参照Mura(1987)的文章.

因此，介质等效弹性刚度张量可表示为

为确保介质在任意孔隙度时的双相连通性，需将SCA模型与DEM模型结合使用(Hornby et al.，1994; Hu et al.，2013).即首先用SCA模型计算孔隙度较高(如50%)的等效介质弹性模量，通过使用DEM模型逐渐加入固体矿物减小孔隙度到指定值，计算新的等效介质的弹性模量，可用差分方程表示为

其中，C ^DEM表示由各向异性DEM模型求得的等效介质刚度张量，v为包含物的体积分数，C ⁱ为包含物的弹性刚度张量.Ĝ为与包含物形状有关的弹性张量，I 是单位张量.该差分方程可通过Runge-Kutta方法求解.

本文中假设岩石孔隙为具有单一纵横比的理想椭球体，将其纵横比看作自由参数，以纵波速度为约束条件进行求取，即将模型计算的纵波速度与测井实测纵波速度相比较，通过最小二乘法来寻找使二者误差足够小的最佳的等效孔隙纵横比.

(4)饱和岩石弹性参数的计算：使用Brown-Korringa方程(Brown and Korringa，1975)将流体加入有一定连通性的干燥微孔隙中，计算饱和岩石弹性参数，进而计算页岩的纵横波速度及其他参数.

Brown和Korringa(1975)研究了各向异性岩石的流体替换，得到饱和流体岩石等效模量：

应力场分布的精确预测对优选有利压裂区有重要的意义.由于页岩通常具有良好的分层结构，可等效为VTI介质，其岩石力学性质也表现出横向各向同性，因此，对页岩地层地应力的预测也可通过横向各向同性模型进行研究(图 3).Suárez-Rivera等(2009)推导了横向各向同性介质地应力计算公式：

图3

Fig. 3

图 3 页岩地层地应力预测模型示意图 Fig. 3 The stress prediction model of shale formation

式中，σ_v为垂直地应力;σ_H和σ_h分别为水平最大和最小地应力;P_p为地层孔隙压力;ε_H为水平方向的构造应变；E_h和E_v分别为水平和垂直方向的杨氏模量；υ_h和υ_v为水平和垂直方向的泊松比.

由于应力应变关系存在两种等效的参数表示形式，即使用工程参数(杨氏模量、泊松比等)表示和使用声学参数(弹性刚度张量)表示，对于VTI介质，二者相互转化关系可表示为

式中，E_h、E_v、υ_h、υ_v分别为水平和垂直方向的杨氏模量和泊松比，C_ij为弹性刚度张量.

由于声学参数更易通过建立的岩石物理模型求 得，因此，依据式(13)—(16)可将式(11)、(12)转化为

国内外学者对于垂直地应力σ_v、地层孔隙压力P_p、水平方向的构造应变ε_H提出了多种计算方式(时贤等，2012；王越等，2014)，如σ_v可通过密度的深度积分进行计算，εH在同一区块内为常数，且不随岩性及位置变化，可通过室内岩石试验获得；P_p可通过Eaton法计算等，但对于各弹性刚度张量(C_ij)难以准确确定.本文利用建立的页岩气岩石物理等效模型可较为准确地计算各弹性刚度张量，进而可实现页岩水平地应力的精确预测.3 实例分析

依据Castagna等(1993)的研究，岩石的纵横波速度存在较为明显的线性或非线性经验关系，我们以某页岩气工区井A为例，分别拟合了其线性与非线性经验公式(图 4)，并将该经验公式应用到井A的横波速度预测中，预测结果如图 5所示.另外，将该井测井解释得到的各矿物含量、含水饱和度及孔隙度信息(图 6)输入本文建立的页岩岩石物理等效模型，计算纵横波速度，与测井结果进行对比(图 7).计算过程中所需矿物及流体参数见表 1.

表1(Table 1)

表 1 所用的矿物组分参数 Table 1 Mineral constituent parameters table 矿物或流体 密度/(g·cm-3) K/GPa μ/GPa 方解石 2.71 76.8 32 石英 2.65 36.6 45 泥质 2.6 21 7 有机质 1.3 2.9 2.7 水 1.02 2.25 0 气 0.00065 0.00013 0

表 1 所用的矿物组分参数 Table 1 Mineral constituent parameters table

图4

Fig. 4

图 4 井A纵横波速度的线性和非线性拟合 (a) 线性拟合； (b) 非线性拟合.Beeler等(1994)实验室摩擦实验结果.实线代表用Dieterich law 对滑动-保持-滑动实验数值模拟的结果. Fig. 4 The linear and nonlinear fitting of well A (a) The linear fitting; (b) The nonlinear fitting.

图5

Fig. 5

图 5 基于线性和非线性拟合的横波速度预测结果 (a)线性拟合；(b)非线性拟合. (a)(b)中自左至右均为纵波速度、横波速度(蓝色曲线为测井值，红色曲线为预测值)、横波速度相对误差. Fig. 5 The prediction of S-wave velocity based on linear and nonlinear fitting (a) The linear fitting; (b) The nonlinear fitting. From left to right of (a) (b) are the P-wave velocity, S-wave velocity (blue curve is logging value and red one is the prediction) and relative error of S-wave velocity.

图6

Fig. 6

图 6 井A矿物体积分数、孔隙度及含水饱和度信息 Fig. 6 The mineral volume fraction, porosity and water saturation of well A

由图 7可知，模型计算得到的纵横波速度与测井结果吻合较好，误差很小，绝大部分小于0.05.与图 5比较发现，利用页岩岩石物理等效模型计算的 横波速度有更高的精度，验证了该模型的可靠适用性.

图7

Fig. 7

图 7 基于页岩岩石物理等效模型的横波速度预测结果 自左至右为：纵波速度(蓝色曲线为测井值，红色曲线为模型计算值)、纵波速度相对误差、横波速度(蓝色曲线为测井值，红色曲线为模型计算值)、横波速度相对误差. Fig. 7 The prediction of S-wave velocity based on shale rock physics equivalent model From left to right: P-wave velocity (the blue curve is logging value and the red one is the prediction), relative error of P-wave velocity, S-wave velocity (the blue curve is logging value and the red one is the prediction) and relative error of S-wave velocity.

由图 8可见，水平方向纵横波速度与垂直方向有一定程度不同，显示地层存在VTI介质特征，这也反映了我们模型的VTI假设具有合理性.

图8

Fig. 8

图 8 不同方向纵横波速度预测结果：两图分别为纵 波速度和横波速度，其中蓝色曲线为水平方向速度，红色曲线为垂直方向速度 Fig. 8 The prediction of P-wave and S-wave velocity in different direction：The two figures are P-wave velocity and S-wave velocity. Among them the blue curves are horizontal velocity and the red curves are vertical velocity

此外，利用建立的岩石物理模型还可获得井A中计算地应力所需的弹性刚度张量C₁₁、C₃₃、C₁₂、C₁₃，如图 9所示.

图9

Fig. 9

图 9 岩石物理等效模型求得的弹性刚度张量(自左至右分别为：C11、C33、C12、C13) Fig. 9 Elastic stiffness tensors calculated by rock physics equivalent model (From left to right: C11、C33、C12、C13)

利用前文提出的方法进行该井垂直地应力σ_v、地层孔隙压力P_p的求取，结果见图 10.另外，实验 室测得的该地区水平方向的构造应变ε_H均值约为 0.3.

图10

Fig. 10

图 10 井A垂直地应力及地层孔隙压力 Fig. 10 Vertical stress and formation pore pressure of well A

至此，我们求得了计算井A最小水平地应力所需的全部参数，依据公式(18)进行最小水平地应力的求取，并将其与各向同性模型求得的最小水平地应力及地层瞬间闭合压力(ISIP)进行对比，结果见图 11.

图11

Fig. 11

图 11 利用各向同性模型和本文模型得到的井A最小水平地应力(标出部分为储层位置) Fig. 11 Minimum horizontal stress of well A calculated by isotropic model and the model in this paper (The marked part is reservoir position)

根据Anderson理论(Anderson，1951)，最小水 平地应力总是小于或等于垂直地应力，同时由于未发生水力压裂，最小水平地应力应大于孔隙压力.由图 11，利用各向同性模型与本文模型计算得到的最小水平地应力数值均小于垂直地应力而大于地层孔隙压力，可见模型构建具有合理性.另外，该井分别测得了2338 m和2373 m处的地层瞬间闭合压力，分别为40.51 MPa和44.3 MPa，将其与模型计算结果比较可知，使用本文中方法求得的最小水平地应力与其吻合较好，且储层与围岩应力差异较大，可较 好指示储层位置(图中标出区域).而各向同性模型计算得到的地应力偏小且储层与围岩应力差异较小，储层位置不明显.因此，基于本文中建立的页岩地层岩石物理模型进行最小水平地应力预测是行之有效的.4 结论和认识

本文通过对研究工区页岩地层的各向异性特征及矿物、孔隙、流体特征的综合分析，将其等效为VTI介质，提出了页岩地层岩石物理等效模型的建立流程，并将其应用于工区横波速度的估算.通过分析其与测井数据的吻合程度以及与经验公式计算结果的对比，验证了此模型在该工区的可靠适用性.同样，对于其他成层性较好、矿物、孔隙、流体等相似的页岩地层，该模型同样是适用的.但因模型中假设矿物种类较少且孔隙为随机分布，该模型对于矿物种类复杂，存在较多垂向裂缝(具有HTI介质特征)的 页岩地层是不适用的，对该类地层的建模方法我们仍需进一步研究.

以建立的岩石物理等效模型为基础，计算了各弹性刚度张量，以VTI介质地应力计算公式为依据，进行了井A最小水平地应力的预测，较各向同 性模型预测的地应力而言，本文中方法预测的最小 水平地应力与地层瞬间闭合压力吻合更好，且储层位置更为明显.最小水平地应力的有效预测可为页岩地层压裂提供一定的依据，但不能为压裂的进行提供完全指导，在实际压裂区域优选过程中，应全面考虑地应力、岩石脆性、含油气性等多方面的因素，提高压裂的有效性，达到页岩气增产的目的.