CN111460602B - 基于岩石物理建模的横观各向同性地层地应力预测方法 - Google Patents

Abstract

本发明公开了一种基于岩石物理建模的横观各向同性地层地应力预测方法，获取页岩物性参数，建立横观各向同性岩石物理模型，计算刚度系数、动静态岩石弹性参数、Biot系数，基于Eaton法利用纵波时差来预测地层孔隙压力，最后计算横观各向同性页岩地层最大、最小水平主应力。本发明构建了基于岩石物理建模的地应力评价体系，提高了静态岩石弹性参数的转换精度和地应力的预测精度，解决了横观各向同性页岩地层的地应力评价难题。利用该发明计算的横观各向同性地层地应力可以为后续的横观各向同性页岩地层井壁稳定性分析、压裂起裂分析等提供基础数据。

Description

基于岩石物理建模的横观各向同性地层地应力预测方法

技术领域

本发明涉及基于岩石物理建模的横观各向同性地层地应力预测方法，属于岩石物理技术领域。

背景技术

地应力是指存在于地壳中的应力，其来源主要包括上覆岩层重力、构造应力、水压附加应力、地温附加应力、化学附加应力等，通常采用垂向地应力、最大水平地应力、最小水平地应力三个参数来表征。对于钻井工程而言，地应力是井壁稳定模型求解的边界条件，其大小直接决定了稳定井壁所需的临界井筒压力或钻井液当量密度，而地应力方向直接影响水平井钻井方位；对于水力压裂而言，地应力直接决定了地层中裂缝起裂、扩展和延伸的力学行为，其大小直接决定了水力压力起裂压力、扩展压力、水力压裂设备选型，其而地应力方向直接影响裂缝起裂与扩展方向，最终直接影响水力压裂效果。

目前，石油工程领域获取地应力的方法较多，可以分为室内测试、矿场测试、测井解释以及其它方法等四个大类：①室内测试主要包括差应变测试、Kaiser声发射测试、波速各向异性测试、古地磁定向测试等；②矿场测试主要包括地破实验、水压致裂法、应力解除法、应力恢复法等；③测井解释主要包括声波测井解释、井壁崩落法、井壁诱导缝解释等；④而其它方法包括地质资料分析、地震资料预测、地应力数值模拟等。其中，矿场测试和室内测试的精度最高、测井解释次之、其它方法最低。但是，矿场测试和室内测试成本高、耗时长，而且只能获得取芯深度点的地应力大小和方向，所测得的数据十分有限；而测井解释具有纵向分辨率高、测井资料相对连续的特点，可以解释得到沿井深相对连续分布的地应力剖面，再结合多种方法和数据进行测井解释剖面的刻度，在一定程度上提升了地应力计算精度，这使得该方法被现场广泛采用。

通常，垂向地应力通过中子密度测井资料进行积分求取，因此，通过测井资料解释求取地应力的难点在于求取最大和最小水平地应力。目前，最大和最小水平地应力的常规测井解释模型主要包括单轴应变模型、黄氏模型、组合弹簧模型和倾斜地层模型，但是，这些常规模型都是以各向同性介质为基础，但实际页岩地层是各向异性介质，因此，常规测井解释无法准确求取地层地应力。对于各向异性地层中地应力的测井解释，Thiercelin和Plumb在各向同性地应力模型的基础上，区分了垂直和水平方向上力学性质的差别，建立了横观各向同性地层地应力模型；Higgins等采用横向各向同性地层地应力模型，对Baxter页岩进行地应力评价，进一步证实了各向异性的显著影响；Sayers、Gholami及邓金根等相继采用横向各向同性地层地应力模型开展了页岩地应力评价。但是，横向各向同性地层地应力测井计算的关键在于计算地层岩石的5个独立刚度系数，若能准确求出刚度系数C₁₁、C₃₃、C₄₄、C₆₆和C₁₃，便可以通过横向各向同性地层地应力模型计算出最小和最大水平地应力。其中，C₃₃和C₄₄可以从阵列声波资料提供的纵波、垂向横波和中子密度测井曲线计算，C₆₆可以通过斯通利波反演水平横波的方法计算，而C₁₁和C₁₃则只能通过室内测试或经验关系计算。然而，对于页岩地层，斯通利波对水平横波的灵敏度降低甚至消失，导致C₆₆的计算精度大幅降低，甚至完全错误，而且，C₁₁和C₁₃的精度也非常低。此外，由于阵列声波测井工具尺寸大，受页岩地层水平井易垮塌、水平段测井工具下入困难等因素影响，通常不能获得阵列声波测井资料，这使得地层岩石5个独立刚度系数的计算受限，这也严重制约了横向各向同性地层地应力测井计算的推广应用。为此，Gray和Goodway等将该模型与地震反演算法相结合，通过地震资料提取的纵横波数据反演计算这5个独立刚度系数，实现了页岩地层区域地应力的评价。但是，地震资料的尺度大、精度低，这种方法通常用于研究区域地应力整体分布情况，通过地震资料提取的单井地应力剖面数据精度不高，难以直接用于页岩气水平井井壁稳定和水力压力设计与施工。因此，发展一种以常规测井资料为基础，准确计算横观各向同性地层地应力的方法非常迫切。

发明内容

针对上述问题，本发明主要是克服现有技术中的不足之处，提出基于岩石物理建模的横观各向同性地层地应力预测方法。

本发明解决上述技术问题所提供的技术方案是：基于岩石物理建模的横观各向同性地层地应力预测方法，包括以下步骤：

步骤S1、通过测井、录井以及室内测试得到页岩物性参数，所述页岩物性参数包括矿物组成及含量、地层流体组成及含量、孔隙度、饱和度；

步骤S2、建立横观各向同性岩石物理模型；

步骤S3、根据横观各向同性岩石物理模型计算横观各向同性的刚度系数；

步骤S4、根据刚度系数与岩石弹性参数的关系式，计算动态岩石弹性参数；

式中：E₁为平行于层理面杨氏模量，GPa；E₃为垂直于层理面的杨氏模量，GPa；ν₁₂为平行于层理面，沿最小水平主应力方向的泊松比，无因次；ν₁₃为平行于层理面，沿垂向方向的泊松比，无因次；ν₃₁为垂直于层理面方向的泊松比，无因次；c₁₁、c₁₂、c₁₃、c₃₃为刚度矩阵中的刚度系数，GPa；

步骤S5、利用动态岩石弹性参数，并结合岩石弹性参数的动静态转换关系来计算静态岩石弹性参数；

步骤S6、根据下式计算横观各向同性岩石物理模型中各向异性的Biot系数；

式中：

为干页岩的刚度系数，GPa；

为背景介质的刚度系数，GPa；α₁₁、α₃₃为平行于层理面和垂直于层理面方向的Biot系数，无因次；

步骤S7、基于Eaton法利用纵波时差来预测地层孔隙压力；

P_p＝σ_v-(σ_v-p_w)(AC_n/AC)^c

式中：σ_v为垂向地应力，MPa；p_w为地层水静液柱压力，MPa；AC_n为该深度点正常压实趋势线纵波时差，μs/ft；AC为实际的纵波时差，μs/ft；P_p为地层孔隙压力，MPa；

步骤S8、利用以下关系式，计算横观各向同性页岩地层最大、最小水平主应力；

式中：σ_H为最大水平主应力，MPa；σ_h为最小水平主应力，MPa；P_p为地层孔隙压力，MPa；E₁为平行于层理面杨氏模量，GPa；E₃为垂直于层理面的杨氏模量，GPa；ν₁₂为平行于层理面，沿最小水平主应力方向的泊松比，无因次；ν₁₃为平行于层理面，沿垂向方向的泊松比，无因次；σ_v为垂向地应力，MPa；ε_H为最大水平主应力方向所对应的最大水平主应变，无因次；ε_h为最小水平主应力方向所对应的最小水平主应变，无因次。

进一步的技术方案是，所述横观各向同性岩石物理模型的具体建立过程为：

步骤S21、利用Hashin-Shtrikman界限计算基质矿物混合物的等效弹性模量；

步骤S22、利用各向同性SCA模型和各向同性DEM模型，将无机质孔隙加入到基质矿物混合物中得到干基质，建立干基质页岩岩石物理模型并计算干基质的等效弹性模量；

步骤S23、利用各向异性SCA模型和各向异性DEM模型计算有机质矿物混合物的等效弹性刚度张量；

步骤S24、利用CL系数描述有机质矿物的成层分布，结合Bond变换得到取向分布的有机质矿物混合物等效弹性刚度张量；

步骤S25、利用各向异性DEM模型，将有机质孔隙加入到取向分布的有机质矿物混合物中得到干有机质，建立干有机质页岩岩石物理模型并计算干有机质的等效弹性刚度张量；

步骤S26、利用各向异性DEM模型，将干基质添加到干有机质中，得到包含空孔隙的干页岩，建立干页岩岩石物理模型并计算干页岩的刚度张量；

步骤S27、根据含气饱和度、含水饱和度，利用Wood公式计算有机质孔隙内气-水混合物的体积模量，并转换成刚度张量；并结合干页岩呈VTI性质的情况，利用Brown-Korringa模型，将混合流体添加到干页岩中得到饱和流体页岩，建立饱和流体页岩岩石物理模型并计算饱和流体页岩的刚度张量。

进一步的技术方案是，所述步骤S21的具体过程是：利用Hashin-Shtrikman界限计算基质矿物混合物体积模量和剪切模量的上下限，并取上下限平均作为基质矿物混合物的体积模量和剪切模量。

进一步的技术方案是，所述Hashin-Shtrikman界限计算基质矿物混合物的公式如下：

式中：K^HS+为基质矿物混合物的等效体积模量上限，GPa；K^HS-为基质矿物混合物的等效体积模量下限，GPa；μ^HS+为基质矿物混合物的等效剪切积模量上限，GPa；μ^HS-为基质矿物混合物的等效剪切模量下限，GPa。

进一步的技术方案是，所述步骤S22中各向同性SCA模型的计算公式如下：

式中：v_i为第i种材料的体积分数，无因次；P^*i为第i种材料的几何因子1，无因次；Q^*i为第i种材料的几何因子2；

为等效的体积模量，GPa；

为等效的剪切模量，GPa；

各向同性DEM模型的计算公式如下：

式中：K₁、μ₁分别为背景介质的体积模量和剪切模量，GPa；K₂、μ₂分别为包裹体的体积模量和剪切模量，GPa；v为包裹体的体积分数，无因次。

进一步的技术方案是，所述步骤S23中各向异性SCA模型的计算公式如下：

式中：

为SCA模型等效的刚度张量，GPa；

为第n相物质的Eshelby刚度张量，无因次；

为四阶单位刚度张量，无因；

为第n相物质的刚度张量，GPa；v_n为第n相物质的体积分数，无因次；

各向异性DEM模型的计算公式如下：

式中：

为背景介质刚度张量，GPa；

为包裹体刚度张量，GPa；

为包裹体的Eshelby刚度张量，无因次；

为四阶单位刚度张量，无因次；v被添加项的体积，小数。

进一步的技术方案是，所述步骤S27中各向异性Schoenberg线性滑动模型计算公式如下：

式中：C为添加裂缝后的干ORT页岩刚度矩阵，GPa；c_11b、c_12b、c_13b、c_33b、c_44b、c_66b为不含垂向裂缝的干VTI页岩刚度系数，GPa；Δ_N、Δ_V、Δ_H描述裂缝特征的弱度，无因次。

进一步的技术方案是，所述步骤S27中Wood公式如下所示：

式中：K_g为气体的体积模量，GPa；K_w为水的体积模量，GPa；K_f为混合流体的等效体积模量，GPa；s_g为含气饱和度，无因次；μ_f为混合流体的等效体积模量，GPa；

Brown-Korringa模型的计算公式如下：

式中：

为饱和岩石的柔度张量，GPa^-1；

为干岩石骨架的柔度张量，GPa^-1；

为岩石基质的柔度张量，GPa^-1；β_fl、β_gr分别为流体和岩石基质的可压缩系数，GPa^-1；φ为孔隙度，无因次；

进一步的技术方案是，所述步骤S5中岩石弹性参数的动静态转换关系通过测量的波速计算岩心的动态岩石弹性参数并和三轴压缩实验得到的静态岩石弹性参数建立得到。

本发明具有以下有益效果：本发明构建了基于岩石物理建模的地应力评价体系，提高了静态岩石弹性参数的转换精度和地应力的预测精度，解决了横观各向同性页岩地层的地应力评价难题；利用该发明计算的横观各向同性地层地应力可以为后续的横观各向同性页岩地层井壁稳定性分析、压裂起裂分析等提供基础数据。

附图说明

图1是本发明的流程框图；

图2是本发明的建模流程图；

图3是x1井矿物及流体的测井解释图；

图4是岩石弹性参数的动、静态转换图；

图5是x1井岩石物理建模预测结果以及同测井测量结果的对比图；

图6是岩石静态弹性参数和各向异性Biot系数图；

图7是孔隙压力和各类模型计算的地应力图。

具体实施方式

下面结合实施例和附图对本发明做更进一步的说明。

如图1所示，本发明的基于岩石物理建模的横观各向同性地层地应力预测方法，具体实施过程如下：

步骤S1、结合测井、录井以及室内测试，得到建模所需的页岩物性参数，其x1井矿物及流体的测井解释图如图3所示，所述页岩物性参数包括矿物组成及含量、地层流体组成及含量、孔隙度、饱和度；

利用测井、录井以及室内测试等手段，获取页岩的矿物组分及相对含量，矿物组分包括基质矿物和有机质矿物两大类，基质矿物主要包括石英、长石、方解石、白云石、黄铁矿等，有机质矿物主要黏土和干酪根；

再利用室内测试或测井，获取页岩孔隙度，以及孔隙中所含的页岩气及水对应的包含度；页岩孔隙中的气体分为吸附气和游离气，基于Xu-White模型的原理，将孔隙分为无机质孔隙和有机质孔隙两类，其中无机质孔隙包含了粒内孔、粒间孔以及微裂缝；利用游离气和吸附气相对含量计算无机质孔隙度和有机质孔隙度；

吸附气的含量可由兰格缪尔等温吸附方程进行计算，具体公式如下：

式中，g_a为吸附气含量，m³/t；V_lc为在油藏温度下经过TOC校正后的兰格缪尔体积，m³/t；P_lt为油藏温度下的兰格缪尔压力，Mpa；p为储层压力，Mpa。

游离气的含量可由以下方程进行计算：

式中，g_f为游离气含量，m³/t；ψ为常数，无因次；ρ为体积密度，g/cm³；B_g为气体体积压缩系数，无因次；φ_e为有效孔隙度，无因次；S_w为有效孔隙度，无因次。

有机质孔隙度和无机质孔隙度的计算方式如下：

式中，φ_o为有机质孔隙度，无因次；φ_g为无机质孔隙度，无因次；φ为总孔隙度，无因次；g_a为吸附气含量，m³/t；g_f为游离气含量，m³/t；

步骤S2、建立横观各向同性岩石物理模型，其具体过程如图2所示；

步骤S21、利用Hashin-Shtrikman界限计算基质矿物混合物的等效弹性模量；

利用Hashin-Shtrikman界限计算基质矿物混合物体积模量和剪切模量的上下限，并取上下限平均作为基质矿物混合物的体积模量和剪切模量；基质矿物是指除了黏土和干酪根以外的矿物，主要包括：石英、长石、方解石、白云石、黄铁矿等；

所述Hashin-Shtrikman界限计算基质矿物混合物的公式如下：

式中：K^HS+为基质矿物混合物的等效体积模量上限，GPa；K^HS-为基质矿物混合物的等效体积模量下限，GPa；μ^HS+为基质矿物混合物的等效剪切积模量上限，GPa；μ^HS-为基质矿物混合物的等效剪切模量下限，GPa；

本实施例中页岩岩石矿物及流体弹性模量如表1所示：

表1

步骤S22、利用各向同性SCA模型和各向同性DEM模型，将无机质孔隙加入到基质矿物混合物中得到干基质，建立干基质页岩岩石物理模型并计算干基质的等效弹性模量；

具体是给定页岩基质岩石的临界孔隙度，利用各向同性SCA模型计算干岩石在临界孔隙度时候的体积模量和剪切模量，然后再利用各向同性DEM模型将孔隙度逐步调整至真实无机质孔隙度φ_g，以计算真实无机质孔隙度时的干页岩基质的体积模量和剪切模量；

各向同性SCA模型的计算公式如下：

式中：v_i为第i种材料的体积分数，无因次；P^*i为第i种材料的几何因子1，无因次；Q^*i为第i种材料的几何因子2，无因次；

为等效的体积模量，GPa；

为等效的剪切模量，GPa；

各向同性DEM模型的计算公式如下：

式中：K₁、μ₁分别为背景介质的体积模量和剪切模量，GPa；K₂、μ₂分布为包裹体的体积模量和剪切模量，GPa；v为包裹体的含量，无因次；初始条件满足K^*(0)＝K₁，μ^*(0)＝μ₁；

步骤S23、利用各向异性SCA模型和各向异性DEM模型计算有机质矿物混合物的等效弹性刚度张量；

具体的是首先利用各向异性SCA模型计算干酪根和黏土含量各占50％情况下的弹性特征，然后用各向异性DEM模型将两者的体积含量调整至对应的百分比，这既保证了干酪根和黏土的相互连通，又避免了由于干酪根和黏土添加顺序的不同而导致混合物弹性模量的不对称；

各向异性SCA模型的计算公式如下：

式中：

为SCA模型等效的刚度张量，GPa；

为第n相物质的Eshelby刚度张量，无因次；

为四阶单位刚度张量，无因；

为第n相物质的刚度张量，GPa；v_n为第n相物质的体积分数，无因次；

各向异性DEM模型的计算公式如下：

式中：

为背景介质刚度张量，GPa；

为包裹体刚度张量，GPa；

为包裹体的Eshelby刚度张量，无因次；

为四阶单位刚度张量，无因次；v被添加项的体积，小数；

步骤S24、利用CL系数描述有机质矿物的成层分布，结合Bond变换正态分布和CL系数形成得到取向分布的有机质矿物混合物等效弹性刚度张量；

步骤S241、利用扫描电镜资料，确定有机质颗粒随对称轴的偏转角度，给出表征有机质颗粒偏转角度的排列分布系数CL；

步骤S242、对于偏转到单个角度下的有机质颗粒，利用Bond变换将完全定向排列的有机质颗粒旋转到对应的偏转角度，得到该角度下对应的有机质颗粒的刚度张量；

步骤S243、将分布函数上所有排列角度的有机质颗粒的刚度矩阵利用V-R-H模型进行叠加，得到不同角度叠加的最终的有机质等效弹性刚度张量；

步骤S25、利用各向异性DEM模型，将有机质孔隙加入到取向分布的有机质矿物混合物中得到干有机质，建立干有机质页岩岩石物理模型并计算干有机质的等效弹性刚度张量；

步骤S26、利用各向异性DEM模型，将干基质添加到干有机质中，得到包含空孔隙的干页岩，建立干VTI页岩石物理模型并计算干页岩的刚度张量；

步骤S27、根据含气饱和度、含水饱和度，利用Wood公式计算孔隙内气-水混合物的体积模量，并转换成刚度张量；并结合干页岩呈VTI性质的情况，利用Brown-Korringa模型，将混合流体添加到干页岩中得到饱和流体页岩，建立饱和流体页岩的岩石物理模型；

Wood公式如下所示：

式中：K_g为气体的体积模量，GPa；K_w为水的体积模量，GPa；K_f为混合流体的等效体积模量，GPa；s_g为含气饱和度，无因次；μ_f为混合流体的等效体积模量，GPa；

Brown-Korringa模型的计算公式如下：

式中：

为饱和岩石的柔度张量，GPa^-1；

为干岩石骨架的柔度张量，GPa^-1；

为岩石基质的柔度张量，GPa^-1；β_fl、β_gr分布为流体和岩石基质的可压缩系数，GPa^-1；φ为孔隙度，无因次；

步骤S3、根据横观各向同性岩石物理模型计算横观各向同性的刚度系数，其结果如图5所示，刚度系数包括：c₁₁、c₁₂、c₁₃、c₃₃、c₄₄、c₆₆；

步骤S4、根据刚度系数与岩石弹性参数的关系式，计算动态岩石弹性参数，其弹性参数包括2个杨氏模量、3个泊松比；

式中：E₁为平行于层理面杨氏模量，GPa；E₃为垂直于层理面的杨氏模量，GPa；ν₁₂为平行于层理面，沿最小水平主应力方向的泊松比，无因次；ν₁₃为平行于层理面，沿垂向方向的泊松比，无因次；ν₃₁为垂直于层理面方向的泊松比，无因次；c₁₁、c₁₂、c₁₃、c₃₃为刚度矩阵中的刚度系数，GPa；

步骤S5、利用动态岩石弹性参数，并结合岩石弹性参数的动静态转换关系来计算静态岩石弹性参数；

具体的是在利用波速测量实验测量页岩岩心波速的同时，对岩心开展了三轴岩石力学压缩实验，其结果如图4所示，用测量的波速计算岩心的动态岩石弹性参数并和三轴压缩实验得到的静态岩石弹性参数建立相应的动静态转换关系；

步骤S6、根据下式计算横观各向同性岩石物理模型中各向异性的Biot系数；

式中：

为干页岩的刚度系数，GPa；

为背景介质的刚度系数，GPa；α₁₁、α₃₃为平行于层理面和垂直于层理面方向Biot系数，无因次；

其岩石静态弹性参数和各向异性Biot系数如图6所示；

步骤S7、基于Eaton法利用纵波时差来预测地层孔隙压力；

P_p＝σ_v-(σ_v-p_w)(AC_n/AC)^c

式中：σ_v为垂向地应力，MPa；p_w为地层水静液柱压力，MPa；AC_n为该深度点正常压实趋势线纵波时差，μs/ft；AC为实际的纵波时差，μs/ft；P_p为地层孔隙压力，MPa；

步骤S8、利用以下关系式，计算横观各向同性页岩地层最大、最小水平主应力，其结果如图7所示；

式中：σ_H为最大水平主应力，MPa；σ_h为最小水平主应力，MPa；ε_H为最大水平主应力方向所对应的最大水平主应变，无因次；ε_h为最小水平主应力方向所对应的最小水平主应变，无因次。

本实施例得到的孔隙压力和各类模型计算的地应力如图7所示；

图7(b)为不同地应力模型计算的最大水平主应力，图中A点为通过压裂施工曲线反算的最大水平主应力，其值为88.26MPa。B点为3161.8米处井下岩心声发射实验(基于Kaiser效应)获取的最大水平主应力，其值为95.21MPa。在B点对应的深度点，VTI地应力模型的预测结果为100.15MPa，ISO地应力模型的预测结果为84.79MPa，两种模型预测的地应力与实测结果的误差分别为5.19％、10.94％。

图7(c)为不同地应力模型计算的最小水平主应力，图中A点为通过压裂施工曲线反算的最小水平主应力，其值为64.70MPa。B点为3161.8米处井下岩心声发射实验(基于Kaiser效应)获取的最小水平主应力，其值为74.13MPa。在B点对应的深度点，VTI地应力模型的预测结果为78.40MPa，ISO地应力模型的预测结果为68.09MPa，两种模型预测的地应力与实测结果的误差分别为：5.76％、8.15％。

以上所述，并非对本发明作任何形式上的限制，虽然本发明已通过上述实施例揭示，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容做出些变动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。

Claims (7)

1.一种基于岩石物理建模的横观各向同性地层地应力预测方法，其特征在于，包括以下步骤：

步骤S1、通过测井、录井以及室内测试得到页岩物性参数，所述页岩物性参数包括矿物组成及含量、地层流体组成及含量、孔隙度和饱和度；

步骤S2、建立横观各向同性岩石物理模型；

步骤S3、根据横观各向同性岩石物理模型计算横观各向同性的刚度系数；

步骤S4、根据刚度系数与岩石弹性参数的关系式，计算动态岩石弹性参数；

式中：E₁为平行于层理面杨氏模量，GPa；E₃为垂直于层理面的杨氏模量，GPa；ν₁₂为平行于层理面，沿最小水平主应力方向的泊松比，无因次；ν₁₃为平行于层理面，沿垂向方向的泊松比，无因次；ν₃₁为垂直于层理面方向的泊松比，无因次；c₁₁、c₁₂、c₁₃、c₃₃为刚度矩阵中的刚度系数，GPa；

步骤S5、利用动态岩石弹性参数，并结合岩石弹性参数的动静态转换关系来计算静态岩石弹性参数；

步骤S6、根据下式计算横观各向同性岩石物理模型中各向异性的Biot系数；

式中：

为干页岩的刚度系数，GPa；

为背景介质的刚度系数，GPa；α₁₁、α₃₃为平行于层理面和垂直于层理面方向的Biot系数，无因次；

步骤S7、基于Eaton法利用纵波时差来预测地层孔隙压力；

P_p＝σ_v-(σ_v-p_w)(AC_n/AC)^c

式中：σ_v为垂向地应力，MPa；p_w为地层水静液柱压力，MPa；AC_n为正常压实趋势线纵波时差，μs/ft；AC为实际的纵波时差，μs/ft；P_p为地层孔隙压力，MPa；

步骤S8、利用以下关系式，计算横观各向同性页岩地层最大、最小水平主应力；

式中：σ_H为最大水平主应力，MPa；σ_h为最小水平主应力，MPa；P_p为地层孔隙压力，MPa；E₁为平行于层理面杨氏模量，GPa；E₃为垂直于层理面的杨氏模量，GPa；ν₁₂为平行于层理面，沿最小水平主应力方向的泊松比，无因次；ν₁₃为平行于层理面，沿垂向方向的泊松比，无因次；σ_v为垂向地应力，MPa；ε_H为最大水平主应力方向所对应的最大水平主应变，无因次；ε_h为最小水平主应力方向所对应的最小水平主应变，无因次；

所述横观各向同性岩石物理模型的具体建立过程为：

步骤S21、利用Hashin-Shtrikman界限计算基质矿物混合物的等效弹性模量；

步骤S22、利用各向同性SCA模型和各向同性DEM模型，将无机质孔隙加入到基质矿物混合物中得到干基质，建立干基质页岩岩石物理模型并计算干基质的等效弹性模量；

步骤S23、利用各向异性SCA模型和各向异性DEM模型计算有机质矿物混合物的等效弹性刚度张量；

步骤S24、利用CL系数描述有机质矿物的成层分布，结合Bond变换得到取向分布的有机质矿物混合物等效弹性刚度张量；

步骤S25、利用各向异性DEM模型，将有机质孔隙加入到取向分布的有机质矿物混合物中得到干有机质，建立干有机质页岩岩石物理模型并计算干有机质的等效弹性刚度张量；

步骤S26、利用各向异性DEM模型，将干基质添加到干有机质中，得到包含空孔隙的干页岩，建立干页岩岩石物理模型并计算干页岩的刚度张量；

步骤S27、根据含气饱和度、含水饱和度，利用Wood公式计算有机质孔隙内气-水混合物的体积模量，并转换成刚度张量；并结合干页岩呈VTI性质的情况，利用Brown-Korringa模型，将混合流体添加到干页岩中得到饱和流体页岩，建立饱和流体页岩岩石物理模型并计算饱和流体页岩的刚度张量。

2.根据权利要求1所述的一种基于岩石物理建模的横观各向同性地层地应力预测方法，其特征在于，所述步骤S21的具体过程是：利用Hashin-Shtrikman界限计算基质矿物混合物体积模量和剪切模量的上下限，并取上下限平均作为基质矿物混合物的体积模量和剪切模量。

3.根据权利要求1所述的一种基于岩石物理建模的横观各向同性地层地应力预测方法，其特征在于，所述Hashin-Shtrikman界限计算基质矿物混合物的公式如下：

式中：K^HS+为基质矿物混合物的等效体积模量上限，GPa；K^HS-为基质矿物混合物的等效体积模量下限，GPa；μ^HS+为基质矿物混合物的等效剪切积模量上限，GPa；μ^HS-为基质矿物混合物的等效剪切模量下限，GPa。

4.根据权利要求1所述的一种基于岩石物理建模的横观各向同性地层地应力预测方法，其特征在于，所述步骤S22中各向同性SCA模型的计算公式如下：

式中：v_i为第i种材料的体积分数，无因次；P^*i为第i种材料的几何因子1，无因次；Q^*i为第i种材料的几何因子2；

为等效的体积模量，GPa；

为等效的剪切模量，GPa；

各向同性DEM模型的计算公式如下：

式中：K₁、μ₁分别为背景介质的体积模量和剪切模量，GPa；K₂、μ₂分别为包裹体的体积模量和剪切模量，GPa；v为包裹体的体积分数，无因次。

5.根据权利要求1所述的一种基于岩石物理建模的横观各向同性地层地应力预测方法，其特征在于，所述步骤S23中各向异性SCA模型的计算公式如下：

式中：

为SCA模型等效的刚度张量，GPa；

为第n相物质的Eshelby刚度张量，无因次；

为四阶单位刚度张量，无因次；

为第n相物质的刚度张量，GPa；v_n为第n相物质的体积分数，无因次；

各向异性DEM模型的计算公式如下：

式中：

为背景介质刚度张量，GPa；

为包裹体刚度张量，GPa；

为包裹体的Eshelby刚度张量，无因次；

为四阶单位刚度张量，无因次；v被添加项的体积，小数。

6.根据权利要求1所述的一种基于岩石物理建模的横观各向同性地层地应力预测方法，其特征在于，所述步骤S27中Wood公式如下所示：

式中：K_g为气体的体积模量，GPa；K_w为水的体积模量，GPa；K_f为混合流体的等效体积模量，GPa；s_g为含气饱和度，无因次；μ_f为混合流体的等效体积模量，GPa；

Brown-Korringa模型的计算公式如下：

式中：

为饱和岩石的柔度张量，GPa^-1；

为干岩石骨架的柔度张量，GPa^-1；

为岩石基质的柔度张量，GPa^-1；β_fl、β_gr分别为流体和岩石基质的可压缩系数，GPa^-1；φ为孔隙度，无因次。

7.根据权利要求1所述的一种基于岩石物理建模的横观各向同性地层地应力预测方法，其特征在于，所述步骤S5中岩石弹性参数的动静态转换关系通过测量的波速计算岩心的动态岩石弹性参数并和三轴压缩实验得到的静态岩石弹性参数建立得到。

Images