一种页岩气地层地应力计算方法与流程

本发明涉及非常规油气地震勘探领域，特别涉及一种页岩气地层地应力计算方法，具体公开了一个评价非常规页岩气地层是否可压裂成网的重要指示因子，即正交各向异性介质水平应力差异比(ODHSR，OrthorhombicDifferentialHorizontalStressRatio)。
背景技术：
：地应力是地球固体介质受重力、多种地球构造动力和天文动力以及探掘工程附加动力的作用,在介质内部单元引起响应变形的力学参数，对非常规页岩油气勘探开发具有重要的影响和意义。起初，由于地应力成因复杂且状态多变，通过数学计算和理论推导计算地应力的大小和方向较为困难，开始采用测量的方法得到地应力。但地应力测量只能得到岩体某点的应力大小和方向，所测得数据是有限的，不能得到连续的地应力剖面。因此，另一种获取地应力数据的手段应运而生，即地应力计算及估算方法。如利用测井资料计算地应力、数值模拟地应力以及刚刚兴起的基于地震资料估算地应力等方法，可以得到沿井深连续分布的分层地应力剖面。为了使计算结果准确，可以用地应力实测的数据进行检验和标定地应力计算结果。现有的计算水平应力差异比主要基于具有对称轴的横向各向异性介质理论，即HTI介质假设。首先，确定地层地震数据中的岩石力学参数，利用方位地震数据反演得到岩石力学参数，然后用反演得到的岩石力学参数计算地层的各向异性弹性参数，最后利用计算得到的地层各向异性参数和岩石力学参数估算最大水平应力、最小水平应力以及水平应力差异比DHSR(即，DHSR，DifferentialHorizontalStressRatio)。该方法没有考虑水平层理对地应力的影响，而实际页岩气地层具有很强的垂向各向异性特征。技术实现要素：本发明综合考虑了VTI介质的水平层理特征和HTI介质的裂缝扰动性质，提出了基于正交各向异性介质的水平应力差异比(即，ODHSR，OrthorhombicDifferentialHorizontalStressRatio)的计算公式，实现了ODHSR的求解，这是一种新的页岩气地层地应力评价方法，可有效评价页岩油气储层是否易于压裂成网。为了解决上述技术问题，本发明采用的技术方案是：一种页岩气地层地应力计算方法，包括如下步骤：第一步，基于岩石物理理论的正交各向异性介质地应力方程推导：设正交各向异性介质的刚度矩阵为：c11c12c13000c12c22c23000c13c23c33000000c44000000c55000000c66---(1),]]>式(1)中cij为正交各向异性介质的弹性刚度张量。在弹性形变范围内，利用广义虎克定律描述弹性介质应力和应变的关系，将广义虎克定律的一般形式进行变换，得到应变是应力的函数形式，写成矩阵形式为：ϵ1ϵ2ϵ3ϵ4ϵ5ϵ6=ABC000BDE000CEF000000H000000M000000Nσ1σ2σ3σ4σ5σ6---(12),]]>式(12)中，σi是应力张量，εi是应变张量，A=c22c33-c232c11c22c33+2c12c23c13-c11c232-c33c122-c22c132---(3),]]>B=c23c13-c12c33c11c22c33+2c12c23c13-c11c232-c33c122-c22c132---(4),]]>C=c12c23-c22c13c11c22c33+2c12c23c13-c11c232-c33c122-c22c132---(5),]]>D=c11c33-c132c11c22c33+2c12c23c13-c11c232-c33c122-c22c132---(6),]]>E=c12c13-c11c23c11c22c33+2c12c23c13-c11c232-c33c122-c222c13---(7),]]>F=c11c22-c122c11c22c33+2c12c23c13-c11c232-c33c122-c22c132---(8),]]>H=1c44---(9),]]>M=1c55---(10),]]>N=1c66---(11),]]>假设存在垂直方向的主应力和两个水平方向的应力，且假设地下岩石是有界的，不能移动，则水平方向上的应变为零，εx＝ε1＝Aσx+Bσy+Cσz＝0(13)，εy＝ε2＝Bσx+Dσy+Eσz＝0(14)，由方程(13)和(14)得到最小水平应力与最大水平应力分别为：σx=σzBE-CDAD-B2---(15),]]>σy=σzBC-AEAD-B2---(16),]]>垂直地应力表示为：σz=∫0Hgρ(h)d(h)---(17),]]>方程(17)中，h为地层厚度，H为深度，g为重力加速度，ρ为密度；利用方程(15)和(16)从地震数据直接计算正交各向异性介质的水平应力差异比ODHSR，ODHSR的计算公式如下所示：ODHSR=σy-σxσy=C(B+D)-E(A+B)BC-AE---(18),]]>正交各向异性弹性介质的矩阵有九个相对独立的弹性参数来确定，设用两个垂向速度和七个无量纲的参数来表征正交各向异性介质弹性性质，具体定义为：VP0=c33ρ---(19),]]>VS0=c55ρ---(20),]]>ϵ(1)=c22-c332c33---(21),]]>ϵ(2)=c11-c332c33---(22),]]>γ(1)=c66-c552c55---(23),]]>γ(2)=c66-c442c44---(24),]]>δ(1)=(c23+c44)2-(c33-c44)22c33(c33-c44)---(25),]]>δ(2)=(c13+c55)2-(c33-c55)22c33(c33-c55)---(26),]]>δ(3)=(c12+c66)2-(c11-c66)22c11(c11-c66)---(27),]]>其中，VP0表示为准纵波的垂直速度；VS0表示为极化方向为x方向的准横波垂直速度；ε、γ和δ表示介质各向异性强度的三个无量纲因子，它们与TI介质Thomsen参数意义相同，上标(1)、(2)、(3)分别表示对应面的法线方向为x、y、z方向，将正交各向异性介质的弹性刚度张量表示为：c11=VP02ρ(1+2ϵ(2))---(28),]]>c22=VP02ρ(1+2ϵ(1))---(29),]]>c33=VP02ρ---(30),]]>c44=VS02ρ(1+2γ(1))1+2γ(2)---(31),]]>c55=VS02ρ---(32),]]>c66=VS02ρ(1+2γ(1))---(33),]]>f=1-VS02VP02=c33-c55c33---(34),]]>c12=VP02ρ2(1+2ϵ(2))[(1+2ϵ(2))-(1-f)(1+2γ(1))]δ(3)+[(1+2ϵ(2))-(1-f)(1+2γ(1))]2-VS02ρ(1+2γ(1))---(35),]]>c13=VP02ρ2fδ(2)+f2-VS02ρ---(36),]]>c23=VP02ρ2[1-(1-f)(1+2γ(1))1+2γ(2)]δ(1)+[1-(1-f)(1+2γ(1))1+2γ(2)]2-VS02ρ(1+2γ(1))1+2γ(2)---(37),]]>将用各向异性参数表示的刚度张量cij代入方程(3)到(11)，在将用各向异性参数表示的方程(3)到(11)代入到方程(18)中，就得到用各向异性参数表示的正交各向异性介质的水平应力差异比ODHSR：ODHSR=C(B+D)-E(A+B)BC-AE=(c12c23-c22c13)[c13(c23-c13)+c33(c11-c12)]-(c12c13-c11c23)[c33(c22-c12)+c23(c13-c23)(c23c13-c12c33)(c12c23-c22c13)-(c22c33-c232)(c12c13-c11c23)---(38);]]>第二步，正交各向异性介质地应力近似表征：正交各向异性介质主应力的近似表达式：σx=σzVP022fδ+f2-VS02ZN[VP04ρ(1+2ϵ)-ρ(VP022fδ+f2-VS02)2]+VP02---(42),]]>σy=σz(VP022fδ+f2-VS02)[1+2VS02ρ(2γ+1)ZN]ZN[VP04ρ(1+2ϵ)-ρ(VP022fδ+f2-VS02)2]+VP02---(43),]]>利用最大水平地应力σy和最小水平地应力σx得到正交各向异性介质的水平应力差异比ODHSR近似式为：ODHSR=2VS02ρ(2γ+1)ZN1+2VS02ρ(2γ+1)ZN---(44),]]>其中，M＝λ+2μ(45)，g=VS02VP02---(46),]]>ΔN=-ϵ(V)2g(1-g)---(47),]]>ZN=ΔNM(1-ΔN)---(48);]]>第三步，利用方位叠前3D地震数据反演得到的岩石力学参数和各向异性参数，然后根据公式(44)，估算正交各向异性介质的水平应力差异比ODHSR与HTI假设条件下计算的DHSR相比，本发明方法同时考虑了各向异性垂向与水平对称轴的作用，基于正交各向异性介质理论，推导出了水平应力差异比岩石物理参数与弹性参数及各向异性参数之间的良好关系，并提出了正交各向异性介质的水平应力差异比这个概念(即ODHSR，OrthorhombicDifferentialHorizontalStressRatio)，使页岩的地应力评价具有符合实际地质意义，对页岩油气储层的勘探开发具有重要的意义附图说明图1本发明一种页岩气地层地应力计算方法的示意图。图2A井的弹性参数和各向异性参数。图3a未消去VTI特征影响的DHSR和ODHSR对比。图3b消去VTI特征影响的DHSR和ODHSR对比。图4a计算HTI介质的DHSR。图4b计算OA介质的ODHSR。具体实施方式本发明提出的一种页岩气地层地应力计算方法的如图1所示：第一步：基于岩石物理理论的正交各向异性介质地应力方程推导正交各向异性介质是由水平层状互层矿物组成的薄层和沿水平对称轴方向排列的裂缝共同构成的各向异性介质，它由三个相互垂直的对称面组成，其刚度矩阵为：c11c12c13000c12c22c23000c13c23c33000000c44000000c55000000c66---(1)]]>(1)式中cij为正交各向异性介质(简称OA介质)的弹性刚度张量。在弹性形变范围内，可以利用广义虎克定律描述弹性介质应力和应变的关系，将广义虎克定律的一般形式进行变换，得到应变是应力的函数形式，写成矩阵形式为：ϵ1ϵ2ϵ3ϵ4ϵ5ϵ6=c22c33-c232c11c22c33+2c12c23c13-c11c232-c33c122-c22c132c23c13-c12c33c11c22c33+2c12c23c13-c11c232-c33c122-c22c132c12c23-c22c13c11c22c33+2c12c23c13-c11c232-c33c122-c22c132000c23c13-c12c33c11c22c33+2c12c23c13-c11c232-c33c122-c22c132c11c33-c132c11c22c33+2c12c23c13-c11c232-c33c122-c22c132c12c13-c11c23c11c22c33+2c12c23c13-c11c232-c33c122-c22c132000c12c23-c22c13c11c22c33+2c12c23c13-c11c232-c33c122-c22c132c12c13-c11c23c11c22c33+2c12c23c13-c11c232-c33c122-c22c132c11c22-c122c11c22c33+2c12c23c13-c11c232-c33c122-c22c1320000001c440000001c550000001c66σ1σ2σ3σ4σ5σ6---(2)]]>其中，σi是应力张量，εi是应变张量。由于该方程的矩阵中元素过于繁冗，为了克服数学上的复杂性和物理上的非直观性，因此设：A=c22c33-c232c11c22c33+2c12c23c13-c11c232-c33c122-c22c132---(3)]]>B=c23c13-c12c33c11c22c33+2c12c23c13-c11c232-c33c122-c22c132---(4)]]>C=c12c23-c22c13c11c22c33+2c12c23c13-c11c232-c33c122-c22c132---(5)]]>D=c11c33-c132c11c22c33+2c12c23c13-c11c232-c33c122-c22c132---(6)]]>E=c12c13-c11c23c11c22c33+2c12c23c13-c11c232-c33c122-c22c132---(7)]]>F=c11c22-c122c11c22c33+2c12c23c13-c11c232-c33c122-c22c132---(8)]]>H=1c44---(9)]]>M=1c55---(10)]]>N=1c66---(11)]]>于是，得到简化后的虎克定律的矩阵形式为：ϵ1ϵ2ϵ3ϵ4ϵ5ϵ6=ABC000BDE000CEF000000H000000M000000Nσ1σ2σ3σ4σ5σ6---(12)]]>假设存在垂直方向的主应力和两个水平方向的应力，且假设地下岩石是有界的，不能移动，则水平方向上的应变为零。εx＝ε1＝Aσx+Bσy+Cσz＝0(13)εy＝ε2＝Bσx+Dσy+Eσz＝0(14)由方程(13)和(14)可以得到最小水平应力σx与最大水平应力σy分别为：σx=σzBE-CDAD-B2---(15)]]>σy=σzBC-AEAD-B2---(16)]]>垂直地应力(σz或σv)表示为：σz=∫0Hgρ(h)d(h)---(17)]]>其中，h为地层厚度，H为深度，g为重力加速度，ρ为密度。从公式可以看出垂直地应力是利用密度测井曲线进行积分估算得到的。此外，利用方程(15)和(16)可以从地震数据直接计算正交各向异性介质的水平应力差异比(ODHSR，OrthorhombicDifferentialHorizontalStressRatio)，它是评价储层是否易于进行水力压裂的重要因子，而不需要知道垂直地应力，低ODHSR值表明其所在区域易于压裂成网，ODHSR的计算公式如下所示：ODHSR=σy-σxσy=C(B+D)-E(A+B)BC-AE---(18)]]>正交各向异性弹性介质的矩阵有九个相对独立的弹性参数来确定，Tsvankin提出用两个垂向速度和七个无量纲的参数来表征OA介质弹性性质。具体定义为：VP0=c33ρ---(19)]]>VS0=c55ρ---(20)]]>ϵ(1)=c22-c332c33---(21)]]>ϵ(2)=c11-c332c33---(22)]]>γ(1)=c66-c552c55---(23)]]>γ(2)=c66-c442c44---(24)]]>δ(1)=(c23+c44)2-(c33-c44)22c33(c33-c44)---(25)]]>δ(2)=(c13+c55)2-(c33-c55)22c33(c33-c55)---(26)]]>δ(3)=(c12+c66)2-(c11-c66)22c11(c11-c66)---(27)]]>其中，VP0表示为准纵波的垂直速度；VS0表示为极化方向为x方向的准横波垂直速度；ε、γ和δ表示介质各向异性强度的三个无量纲因子，它们与TI介质Thomsen参数意义相同，上标(1)、(2)、(3)分别表示对应面的法线方向为x、y、z方向。根据OA介质中的刚度张量与Tsvankin速度各向异性参数的关系，可将OA介质的弹性刚度张量表示为：c11=VP02ρ(1+2ϵ(2))---(28)]]>c22=VP02ρ(1+2ϵ(1))---(29)]]>c33=VP02ρ---(30)]]>c44=VS02ρ(1+2γ(1))1+2γ(2)---(31)]]>c55=VS02ρ---(32)]]>c66=VS02ρ(1+2γ(1))---(33)]]>f=1-VS02VP02=c33-c55c33---(34)]]>c12=VP02ρ2(1+2ϵ(2))[(1+2ϵ(2))-(1-f)(1+2γ(1))]δ(3)+[(1+2ϵ(2))-(1-f)(1+2γ(1))]2-VS02ρ(1+2γ(1))---(35)]]>c13=VP02ρ2fδ(2)+f2-VS02ρ---(36)]]>c23=VP02ρ2[1-(1-f)(1+2γ(1))1+2γ(2)]δ(1)+[1-(1-f)(1+2γ(1))1+2γ(2)]2-VS02ρ(1+2γ(1))1+2γ(2)---(37)]]>最后，将用各向异性参数表示的刚度张量cij代入方程(3)到(11)，在将用各向异性参数表示的方程(3)到(11)代入到方程(18)中，就得到用各向异性参数表示的OA介质的水平应力差异比(ODHSR)：ODHSR=C(B+D)-E(A+B)BC-AE=(c12c23-c22c13)[c13(c23-c13)+c33(c11-c12)]-(c12c13-c11c23)[c33(c22-c12)+c23(c13-c23)](c23c13-c12c33)(c12c23-c22c13)-(c22c33-c232)(c12c13-c11c23)=]]>第二步，正交各向异性介质地应力近似表征由于上面推导的OA介质的主应力和ODHSR公式复杂，不便于实际应用，所以将OA介质的一般公式进行近似简化。由于正交各向异性介质是由水平层状互层矿物组成的薄层和沿水平对称轴方向排列的裂缝共同构成的各向异性介质，所以根据Schoenberg和Sayers的线性滑动理论，可以将OA介质的柔度矩阵近似等价于具有垂直对称轴的横向各向同性背景介质的柔度矩阵与扰动裂缝的柔度矩阵之和。根据Schoenberg和Sayers，扰动裂缝的柔度张量Sf可以写成如下形式：Sf=ZN000000000000000000000000000ZT000000ZT---(39)]]>其中，ZN为法向柔度，ZT为切向柔度。用刚度张量表示的VTI背景介质的刚度矩阵为：Cv=c11bc11b-2c66bc13b000c11b-2c66bc11bc13b000c13bc13bc33b000000c44b000000c44b000000c66b---(40)]]>其中，cijb为VTI介质的刚度张量。根据Schoenberg和Sayers提出的线性滑动理论，可以得到OA介质柔度矩阵的近似式：S=Sv+Sf=c11bc33b-c13b24c66b(c11bc33b-c33bc66b-c13b2)+ZN2c66bc33b+c13b2-c11bc33b4c66b(c11bc33b-c33bc66b-c13b2)-c13bc66b2c66b(c11bc33b-c33bc66b-c13b2)0002c66bc33b+c13b2-c11bc33b4c66b(c11bc33b-c33bc66b-c13b2)c11bc33b-c13b24c66b(c11bc33b-c33bc66b-c13b2)-c13bc66b2c66b(c11bc33b-c33bc66b-c13b2)000-c13bc66b2c66b(c11bc33b-c33bc66b-c13b2)-c13bc66b2c66b(c11bc33b-c33bc66b-c13b2)c11bc66b-c66b2c66b(c11bc33b-c33bc66b-c13b2)0000001c44b0000001c44b+ZT0000001c66bZT---(41)]]>则根据虎克定律，经过一系列推导，可以得到OA介质主应力的近似表达式：σx=σzVP022fδ+f2-VS02ZN[VP04ρ(1+2ϵ)-ρ(VP022fδ+f2-VS02)2]+VP02---(42)]]>σy=σz(VP022fδ+f2-VS02)[1+2VS02ρ(2γ+1)ZN]ZN[VP04ρ(1+2ϵ)-ρ(VP022fδ+f2-VS02)2]+VP02---(43)]]>其中，VP0表示为准纵波的垂直速度；VS0表示为极化方向为x方向的准横波垂直速度；ε、γ和δ表示VTI介质各向异性强度的三个无量纲因子。利用最大水平地应力σy和最小水平地应力σx可以得到OA介质的水平应力差异比(ODHSR)近似式为：ODHSR=2VS02ρ(2γ+1)ZN1+2VS02ρ(2γ+1)ZN---(44)]]>其中，M＝λ+2μ(45)g=VS02VP02---(46)]]>ΔN=-ϵ(V)2g(1-g)---(47)]]>ZN=ΔNM(1-ΔN)---(48)]]>最后利用方位叠前3D地震数据反演得到的岩石力学参数和各向异性参数估算水平应力差异比(ODHSR)。正交各向异性介质求得的ODHSR既有HTI介质的各向异性参数也有VTI介质的各向异性参数，为了验证近似方程的合理性，进一步将ODHSR方程中的VTI介质的各向异性参数消去，我们可以得到HTI情况下DHSR表示形式：ODHSR=2VS02ρ(2γ+1)ZN1+2VS02ρ(2γ+1)ZN---(49)]]>消去VTI介质的各向异性参数，即令VTI介质的各向异性参数γ等于零，可以得到：DHSR=2VS02ρZN1+2VS02ρZN---(50)]]>由于且将方程(50)进行变换，得到：DHSR=EZN1+EZN+v---(51)]]>其中，μ为剪切模量，E为杨氏模量，v为泊松比。方程(51)与Gray(2011)提出的HTI介质计算的水平应力相对变换(DHSR)一致，因此可以证明ODHSR公式的合理性。选取中国某页岩储层工区A井，分别对HTI介质的水平应力差异比(DHSR)和OA介质的水平应力差异比进行计算。A井的纵波速度VP0、横波速度VS0、密度ρ以及各向异性参数ε(V)、δ(V)、γ(V)、γ如图2所示，其中ε(V)、δ(V)、γ(V)是HTI介质的各向异性参数，γ是VTI介质的各向异性参数，它可以通过HTI介质与VTI介质的各向异性参数的关系，由HTI的各向异性参数γ(V)经过转化得到。先计算HTI介质的水平应力差异比(DHSR)，利用A井的弹性参数计算得到岩石力学参数，然后利用计算得到的岩石力学参数和各向异性参数得到地层的DHSR。然后再对OA介质的水平应力差异比(ODHSR)进行计算，利用A井的弹性参数和各向异性参数即可计算地层的ODHSR，结果如图3a、图3b所示，不同介质估算的DHSR的对比，其中虚线为DHSR，实线为ODHSR。从图3aHTI介质DHSR与OA介质ODHSR的对比可以看出，两者的趋势基本一致，ODHSR低值部分与DHSR低值部分相同，ODHSR的高值部分则略高于DHSR高值部分，这是由于VTI介质的各向异性参数γ的影响导致。消去VTI介质的各向异性参数γ的影响，得到的ODHSR与DHSR是完全一致的，如图3b所示，这说明了新建立的正交各向异性介质水平应力差异比的合理性。进一步对实际地震资料进行计算，通过方位叠前地震资料反演可得到HTI介质的岩石力学参数和各向异性参数，再对不同入射角的3D地震数据的地震振幅进行反演得到VTI介质的弹性参数。根据HTI介质与VTI介质的各向异性参数的精确关系，将HTI介质的各向异性参数γ(V)转化为VTI介质的各向异性参数γ，分别计算HTI介质的DHSR和OA介质的ODHSR，通过对比可以看出(如图4a、图4b所示)，本发明估算的ODHSR和DHSR的趋势一致，但因考虑VTI的影响，局部细节有些差异。因此，本发明同时考虑各向异性垂向与水平对称轴的作用，基于正交各向异性介质理论，推导出了水平应力差异比岩石物理参数与弹性参数及各向异性参数之间的关系，实现了对正交各向异性介质水平应力差异比的求解，这是一种新的页岩气地层地应力评价方法，并提出了一个正交各向异性介质水平应力差异比(即，ODHSR，OrthorhombicDifferentialHorizontalStressRatio)的概念。通过将正交各向异性介质估算的ODHSR与具有水平对称轴的横向各向异性介质的DHSR进行对比，证明了本发明的正交各向异性介质计算地应力的可行性。当前第1页1 2 3