引言

页岩气是一种重要的非常规能源，受北美页岩气革命的影响，中国也加入到了页岩气的勘探与开发热潮中^[1-2]。页岩气具有低孔低渗的特点，水平井钻井技术和体积压裂改造是页岩气成功开发的两大关键技术^[3]岩石的力学特性是影响井壁稳定，造成井眼垮塌、漏失的重要因素^[4]在水力压裂过程中，不同空间位置的岩石力学特性展布会影响裂缝的形成与扩展，例如，在胶结程度较弱的层理面往往先于页岩本体开裂，使得水力裂缝在延伸过程中沿层理面优先扩展，从而抑制缝网的形成，降低水力压裂效率^[5]因此，研究岩石的力学特性对油气勘探开发有着重要的意义。

相比横向各向同性(VTI)介质，在地下岩层中正交各向异性(ORT)介质是更为广泛的存在，尤其在沉积盆地中是最普通的情况之一。目前与沉积岩，如页岩相关的岩石力学研究中，主流观点是将地层岩石视为宏观同性介质^[6-8]。海、河水的方向性流动使得岩石颗粒的定向排列，水平最大最小主应力的差异造成两个方向上的裂缝、孔隙具有不同程度的挤压或新裂缝的产生，以及钻井过程中钻遇背(向)斜地层等是造成地层表现出正交各向异性的主要因素^[9-10]。忽视地层的正交各向异性，利用VTI模型获取的泊松比、杨氏模量等岩石力学参数计算地应力，并用于进一步的井壁稳定分析和压裂施工指导等，有时会导致错误的结果。

与VTI介质不同，ORT介质研究相对较少。Hudson给出裂隙介质的本构关系^[11]，希望通过研究地震波在裂缝介质中的响应特征预测储层目标区域的裂缝方位、填充物等微观裂缝参数。Crampin通过将定向排列的充满液体的垂直裂隙加入到周期性薄层的横向各向同性介质中，研究了正交各向异性的横波偏振特征^[12]Bush在巴黎盆地获取的资料中观察到了横波分裂现象，使人们认识到沉积盆地中广泛存在着正交各向异性^[13]Tsvankin将VTI介质的Thomsen参数推广到了ORT介质，便于研究裂缝引起的各向异性^[14]Sehoenberg和Helbig从垂直裂缝介质的弹性方程出发，模拟了ORT介质中3个正交对称面上纵横波的波速面^[15]Bakulin等对裂缝诱导的正交各向异性进行了研究，利用地震波反射数据估计裂缝参数^[16]Cheng改进了Hudson的一阶扰动模型，在Eshelby的基础上建立了Eshelby-Cheng模型^[17]刘恩儒和曾新吾将多种裂缝模型总结为3类，推导了3种模型对应的柔度解析表达式，给出统一的解析形式^[18]Franquet和Rodriguez利用测井资料，基于刚度矩阵参数，建立了正交各向异性地层的地应力计算方法^[19]Alkhalifah利用在声学介质假设下导出的色散关系，获得了正交介质的声波方程^[20]马妮等利用地震资料根据水平应力差异比计算了正交各向异性地层的地应力^[21]物理模型实验方面，Chaedle等最先对正交各向异性样品进行了速度测试，测得了3个主轴和6个45°对角方向上的纵横波速度^[22]李跃等用铝板和有机玻璃板构成一个薄互层等效各向异性模型，讨论了薄层厚度对各向异性的影响^[23]魏建新通过实验观测了正交各向异性介质的声波特征^[10-24]。李军等利用现场实钻岩芯进行了岩石力学性质正交各向异性实验，并建立了正交各向异性地层井壁围岩应力模型^[25-26]。水平井钻井和体积压裂是目前页岩气勘探开发的主要方式，井筒与地层的夹角(井斜角和方位角)随井眼轨迹的变化而变化。而以上研究成果都是基于岩石本构坐标系而得到的，没有考虑岩石力学特性随观测坐标系的变化，得到的结果不便于水平井钻井和压裂增产改造的直接应用。对于VTI介质，由于观测坐标系的不同，造成的岩石力学特性差异已有一些相关研究^[27-28]，而ORT介质则无相应报道。

针对上述问题，考虑实际页岩地层中的各向异性特点，从岩石本构方程出发，推导出ORT介质岩石力学参数的弹性表达式。然后考虑井眼观测坐标系与岩石本构坐标系的不同交角，根据岩石的力学参数定义，建立了观测坐标系下ORT介质的岩石力学参数计算方法，并进行了分析。最后，考虑页岩与砂岩的互相沉积成岩，研究了不同砂岩含量下页岩-砂岩互层在观测坐标系下的岩石力学特性。

1 岩石本构方程的各向异性表征

对于一个一般各向异性的线性弹性固体，应力σ与应变ε之间存在如下线性关系^[29]

$ \mathit{\boldsymbol{\sigma }} = \mathit{\boldsymbol{C\varepsilon }} $

(1)

式中：

σ-应力张量，GPa；

ε-应变张量，无因次；

C-岩石刚度矩阵，GPa。

由于对称性，C最多具有21个独立常数。反过来，应变也可以表示为应力的线性组合

$ \mathit{\boldsymbol{\varepsilon }} = \mathit{\boldsymbol{S\sigma }} $

(2)

式中：S-岩石的柔度矩阵，GPa^-1。

柔度矩阵S和刚度矩阵C互为逆矩阵。对于如图 1所示的正交各向异性岩石，其具有$C_{11}$、$C_{12}$、$C_{13}$、$C_{22}$、$C_{23}$、$C_{33}$、$C_{44}$、$C_{55}$、$C_{66}$共9个独立的弹性刚度参数。

在正交各向异性条件下，刚度矩阵C₀和柔度矩阵S₀可以分别表示如下

$ {{\boldsymbol{C}}}_0 = \left[ {\begin{array}{*{20}{c}} {{C_{11}}}&{{C_{12}}}&{{C_{13}}}&0&0&0\\ {{C_{12}}}&{{C_{22}}}&{{C_{23}}}&0&0&0\\ {{C_{13}}}&{{C_{23}}}&{{C_{33}}}&0&0&0\\ 0&0&0&{{C_{44}}}&0&0\\ 0&0&0&0&{{C_{55}}}&0\\ 0&0&0&0&0&{{C_{66}}} \end{array}} \right] $

(3)

式中：

C₀-正交各向异性条件下岩石刚度矩阵，GPa

$C_{ij}$-岩石刚度系数($i$=1, 2, 3, 4, 5, 6，$j$=1, 2, 3, 4, 5, 6)，GPa。

$ \boldsymbol{S}_{0} = \left[ \begin{array}{l} \frac{1}{{{E_1}}}\quad \frac{{ - {\mu _{21}}}}{{{E_2}}}\frac{{ - {\mu _{31}}}}{{{E_3}}}\quad 0\quad \; 0\quad \; 0\\ \frac{{ - {\mu _{12}}}}{{{E_1}}}\;\;\frac{1}{{{E_2}}}\;\;\frac{{ - {\mu _{32}}}}{{{E_3}}}\quad 0\quad 0\quad \;0\\ \frac{{ - {\mu _{13}}}}{{{E_1}}}\;\frac{{ - {\mu _{23}}}}{{{E_2}}}\;\;\;\frac{1}{{{E_3}}}\quad 0\quad 0\quad \;0\\ 0\qquad 0\qquad 0\quad \frac{1}{{{G_{23}}}}\quad 0\quad 0\\ 0\qquad 0\;\;\;\;\;\;\;0\;\;\;\;\;0\;\;\;\frac{1}{{{G_{13}}}}\;\;0\\ 0\qquad 0\;\;\;\;\;\;\;0\;\;\;\;\;0\;\;\;\;0\;\;\frac{1}{{{G_{12}}}} \end{array} \right] $

(4)

式中：

S₀-正交各向异性条件下岩石柔度矩阵，GPa^-1；

$E_1$-沿$x$轴方向的杨氏模量，GPa；

$E_2$-沿$y$轴方向的杨氏模量，GPa；

$E_3$-沿$z$轴方向的杨氏模量，GPa；

$\mu_{12}$-垂直于$x$轴，沿$y$轴的泊松比，无因次；

$\mu_{13}$-垂直于$x$轴，沿$z$轴的泊松比，无因次；

$\mu_{21}$-垂直于$y$轴，沿$x$轴的泊松比，无因次；

$\mu_{23}$-垂直于$y$轴，沿$z$轴的泊松比，无因次；

$\mu_{31}$-垂直于$z$轴，沿$x$轴的泊松比，无因次；

$\mu_{32}$-垂直于$z$轴，沿$y$轴的泊松比，无因次；

$G_{12}$-$xy$平面的剪切模量，GPa；

$G_{13}$-$xz$平面的剪切模量，GPa；

$G_{23}$-$yz$平面的剪切模量，GPa。

利用Sonic Scanner^[30]测井平台能够获取包括：$C_{33}$、$C_{44}$、$C_{55}$、$C_{66}$在内的4个参数，其余的参数需要通过实验获取。目前，有两类方法可以获取这些参数，岩芯的超声波实验测试和室内三轴岩石力学压缩实验。

用超声波波速测试确定这9个参数需要测量6个方向的纵波波速和3个方向的横波波速^{[22, 31]}，分别包括：(1)沿$x$、$y$、$z$轴方向传播纵波波速$v_{{\rm p}-x}$、$v_{{\rm p}-y}$、$v_{{\rm p}-z}$，如图 2a；(2)沿$z$轴传播，偏振方向为$x$($y$)轴的横波波速$v_{{\rm{sv}}-x}$、$v_{{\rm{sv}}-y}$和沿$y$轴方向传播、偏振方向为$x$轴的横波波速$v_{\rm{sh}}$，如图 2b；(3)在$xy$($zx$、$yz$)平面与$x$($z$、$y$)轴成45°角方向传播的纵波波速$v_{{\rm p}-xy45{°}}$、$v_{{\rm p}-zx45{°}}$、$v_{{\rm p}-yz45{°}}$，如图 2c。

图 2为纵横波波速测量示意图，图中长虚线单向箭头表示波的传播方向，短虚线双向箭头表示波的振动方向。结合岩石的密度$\rho$和测量的波速，9个独立的弹性刚度参数可表示如式(5)、式(6)、式(7)所示

$ \left\{ \begin{array}{l} {C_{11}} = \rho v_{{\rm{p}} - x}^{\rm{2}}\\ {C_{22}} = \rho v_{{\rm{p}} - y}^{\rm{2}}\\ {C_{33}} = \rho v_{{\rm{p}} - z}^{\rm{2}} \end{array} \right. $

(5)

式中：

$C_{11}$，$C_{22}$，$C_{33}$-刚度系数，GPa；

$\rho$-岩石密度，g/cm$^3$；

$v_{{\rm p}-x}$-沿$x$轴传播纵波波速，km/s；

$v_{{\rm p}-y}$-沿$y$轴传播纵波波速，km/s；

$v_{{\rm p}-z}$-沿$z$轴传播纵波波速，km/s。

$ \left\{ \begin{array}{l} {C_{44}} = \rho v_{{\rm{sv}} - x}^{\rm{2}}\\ {C_{55}} = \rho v_{{\rm{sv}} - y}^{\rm{2}}\\ {C_{66}} = \rho v_{{\rm{sh}}}^{\rm{2}} \end{array} \right. $

(6)

式中：

$C_{44}$，$C_{55}$，$C_{66}$-刚度系数，GPa；

$v_{{\rm{sv}}-x}$-沿$z$轴传播，偏振方向为$x$轴的横波波速，km/s；

$v_{{\rm{sv}}-y}$-沿$z$轴方向传播，偏振方向为$y$轴的横波波速，km/s；

$v_{\rm{sh}}$-沿$y$轴方向传播，偏振方向为$x$轴的横波波速，km/s。

$ \left\{ \begin{array}{l} {C_{12}} = \sqrt {{{\left( {2\rho v_{{{\rm{p}} - xy}{{45°}}}^2 - {C_{66}} - \dfrac{{{C_{22}}}}{2} - \dfrac{{{C_{11}}}}{2}} \right)}^2} - \dfrac{{{{\left( {{C_{22}} - {C_{11}}} \right)}^2}}}{4}} - {C_{66}}\\ {C_{13}} = \sqrt {{{\left( {2\rho v_{{{\rm{p}} - zx}{{45°}}}^2 - {C_{55}} - \dfrac{{{C_{11}}}}{2} - \dfrac{{{C_{33}}}}{2}} \right)}^2} - \dfrac{{{{\left( {{C_{11}} - {C_{33}}} \right)}^2}}}{4}} - {C_{55}}\\ {C_{23}} = \sqrt {{{\left( {2\rho v_{{{\rm{p}} - yz}{{45°}}}^2 - {C_{44}} - \dfrac{{{C_{22}}}}{2} - \dfrac{{{C_{33}}}}{2}} \right)}^2} - \dfrac{{{{\left( {{C_{22}} - {C_{33}}} \right)}^2}}}{4}} - {C_{44}} \end{array} \right. $

(7)

式中：

$C_{12}$，$C_{13}$，$C_{23}$-刚度系数，GPa；

$v_{{\rm p}-xy45{°}}$-在$xy$平面，与$x$轴夹角成45°的纵波波速，km/s；

$v_{{\rm p}-zx45{°}}$-在$zx$平面，与$z$轴夹角成45°的纵波波速，km/s；

$v_{{\rm p}-yz45{°}}$-在$yz$平面，与$y$轴夹角成45°的纵波波速，km/s。

如果采用三轴室内岩石力学实验确定这些参数，则需要测量沿上述6个方向取芯岩样的应力应变曲线就能完全确定岩石中的各向异性^[32]

2 岩石力学参数的弹性表征

由上述可知，由室内实验或现场数据可以得到相应的弹性刚度系数，而要获取对应的岩石力学参数(泊松比、杨氏模量)，则需要建立弹性刚度参数与岩石力学参数的关系。将式(3)取逆，并结合式(4)，可以推出岩石力学参数关于岩石刚度系数的表达式，包括3个杨氏模量和6个泊松比。3个杨氏模量的表达式为

$ \left\{ \begin{array}{l} {E_1} = \Delta /(C_{23}^2 - {C_{22}}{C_{33}})\\ {E_2} = \Delta /(C_{13}^2 - {C_{11}}{C_{33}})\\ {E_3} = \Delta /(C_{12}^2 - {C_{11}}{C_{22}}) \end{array} \right. $

(8)

其中：

$ \Delta = {C_{33}}C_{12}^2 + {C_{22}}C_{13}^2 + {C_{11}}C_{23}^2 - \\ \;\;\;\;\;\;\;\;{C_{11}}{C_{22}}{C_{33}} - 2{C_{12}}{C_{13}}{C_{23}} $

3个方向6个泊松比为

$ \left\{ \begin{array}{l} {\mu _{12}} = \left( {{C_{13}}{C_{23}}{\rm{ }} - {\rm{ }}{C_{12}}{C_{33}}} \right)/(C_{23}^2 - {C_{22}}{C_{33}})\\ {\mu _{13}} = \left( {{C_{12}}{C_{23}}{\rm{ }} - {\rm{ }}{C_{13}}{C_{22}}} \right)/(C_{23}^2 - {C_{22}}{C_{33}})\\ {\mu _{21}} = \left( {{C_{13}}{C_{23}}{\rm{ }} - {\rm{ }}{C_{12}}{C_{33}}} \right)/(C_{13}^2 - {C_{11}}{C_{33}})\\ {\mu _{23}} = \left( {{C_{12}}{C_{13}}{\rm{ }} - {\rm{ }}{C_{11}}{C_{23}}} \right)/(C_{13}^2 - {C_{11}}{C_{33}})\\ {\mu _{31}} = \left( {{C_{12}}{C_{23}}{\rm{ }} - {\rm{ }}{C_{13}}{C_{22}}} \right)/(C_{12}^2 - {C_{11}}{C_{22}})\\ {\mu _{32}} = \left( {{C_{12}}{C_{13}}{\rm{ }} - {\rm{ }}{C_{11}}{C_{23}}} \right)/(C_{12}^2 - {C_{11}}{C_{22}}) \end{array} \right. $

(9)

杨氏模量和泊松比的关系满足

$E_2$$\mu_{12}$=$E_1$$\mu_{21}$、$E_3$$\mu_{13}$=$E_1$$\mu_{31}$、$E_2$$\mu_{32}$=$E_3$$\mu_{23}$

由式(9)可知

$\mu_{12}$≠$\mu_{21}$、$\mu_{13}$≠$\mu_{31}$、$\mu_{23}$≠$\mu_{32}$。

当获得了9个独立的弹性刚度参数$C_{11}$、$C_{12}$、$C_{13}$、$C_{22}$、$C_{23}$、$C_{33}$、$C_{44}$、$C_{55}$、$C_{66}$之后，将其代入式(8)、式(9)，便可得到对应的岩石力学参数，这样便建立起了岩石力学参数的弹性刚度计算方法。

对于VTI介质，弹性刚度参数满足

$C_{11}$=$C_{22}$、$C_{13}$=$C_{23}$、$C_{44}$=$C_{55}$

则式(8)、式(9)退化为

$ \left\{ \begin{array}{l} {E_{\rm{v}}}{\rm{ = }}{C_{33}} - 2\dfrac{{{C_{13}}^2}}{{{C_{11}} + {C_{12}}}}\\ {E_{\rm{h}}}{\rm{ = }}\dfrac{{\left( {{C_{11}}{\rm{ - }}{C_{12}}} \right)\left( {{C_{11}}{C_{33}} - 2{C_{13}}^2 + {C_{12}}{C_{33}}} \right)}}{{{C_{11}}{C_{33}} - {C_{13}}^2}}\\ {\mu _{\rm{v}}}{\rm{ = }}\dfrac{{{C_{13}}}}{{{C_{11}} + {C_{12}}}}\\ {\mu _{\rm{h}}}{\rm{ = }}\dfrac{{{C_{12}}{C_{33}} - C_{13}^2}}{{{C_{11}}{C_{33}} - {C_{13}}^2}} \end{array} \right. $

(10)

式中：$E_{\rm{v}}$-垂直于层理面的杨氏模量，GPa；

$E_{\rm{h}}$-平行于层理面的杨氏模量，GPa；

$\mu_{\rm{v}}$-垂直于层理面的泊松比，无因次；

$\mu_{\rm{h}}$-平行于层理面的泊松比，无因次。

3 岩石力学参数的空间变换

在钻井过程中，井眼的方位角和井斜角随时都在变化，这导致了观测坐标系和岩石本构坐标系不一致。尽管室内实验可以获取任意对应方位角和井斜角条件下的岩石力学参数，但岩芯的获取较为困难，测试费用高，并且钻取的角度有限。因此，如何利用测量到的岩石本构坐标系下弹性刚度矩阵，去计算任意观测坐标系下的岩石力学参数，具有很强现实意义。本文利用Bond变换，首先，将岩石本构坐标系下的弹性刚度矩阵转换到观测坐标系，然后，根据岩石力学参数的定义，计算任意方位井斜角条件下的岩石力学参数，具体流程如下。

假定地层岩石的本构坐标系为$O-x_0y_0z_0$，且$Ox_0$轴，$Oy_0$轴和$Oz_0$轴分别与主应力$\sigma_{\rm H}$、$\sigma_{\rm h}$和$\sigma_{\rm v}$方向一致。本构坐标系$O-x_0y_0z_0$先以$Oz_0$为轴旋转方位角$\alpha$，变为$O-x_1y_1z_1$坐标系，再以$Oy_1$为轴旋转井斜角$\beta$，则得到观测坐标系(井筒坐标系)$O-xyz$，如图 3所示。

岩石本构坐标系($O-x_0y_0z_0$)和观测坐标系($O-xyz$)之间的坐标转换关系可由坐标之间的方向余弦$a_{ij}$表示

$ {a_{ij}} = \left[ {\begin{array}{*{20}{c}} {\cos \beta \cos \alpha }&{\cos \beta \sin \alpha }&{ - \sin\beta }\\ { - \sin\alpha }&{\cos \alpha }&0\\ {\sin \beta \cos \alpha }&{\sin \beta \sin \alpha }&{\cos \beta } \end{array}} \right] $

(11)

式中：

$a_{ij}$-方向余弦($i, j$=1, 2, 3)，无因次；

$\alpha$-方位角，弧度；

$\beta$-井斜角，弧度。

假设$\sigma$、$\varepsilon$、C、S和$\sigma_0$、$\varepsilon_0$、C₀、S₀分别为观测坐标系和本构坐标系下的应力、应变、刚度矩阵和柔度矩阵，则观测坐标系和本构坐标系下应力、应变满足[33, 34}

$ \left\{ {\begin{array}{*{20}{l}} {\sigma = \mathit{\boldsymbol{M}}{\sigma _0}}\\ {\varepsilon = \mathit{\boldsymbol{N}}{\varepsilon _0}} \end{array}} \right. $

(12)

式中：M，N-过渡矩阵。

过渡矩阵M、N满足N^-1=M，M矩阵可由$a_{ij}$的表达式计算。

$ {\boldsymbol{M}} = \\ \left[ {\begin{array}{*{20}{c}} {a_{11}^2}&{a_{12}^2}&{a_{13}^2}&{2{a_{12}}{a_{13}}}&{2{a_{11}}{a_{13}}}&{2{a_{11}}{a_{12}}}\\ {a_{21}^2}&{a_{22}^2}&{a_{23}^2}&{2{a_{22}}{a_{23}}}&{2{a_{21}}{a_{13}}}&{2{a_{21}}{a_{22}}}\\ {a_{31}^2}&{a_{32}^2}&{a_{33}^2}&{2{a_{32}}{a_{33}}}&{2{a_{31}}{a_{33}}}&{2{a_{31}}{a_{32}}}\\ {{a_{21}}{a_{31}}}&{{a_{22}}{a_{32}}}&{{a_{23}}{a_{33}}}&{{a_{22}}{a_{33}} + {a_{32}}{a_{23}}}&{{a_{21}}{a_{33}} + {a_{31}}{a_{23}}}&{{a_{21}}{a_{32}} + {a_{31}}{a_{22}}}\\ {{a_{11}}{a_{31}}}&{{a_{12}}{a_{32}}}&{{a_{13}}{a_{33}}}&{{a_{12}}{a_{33}} + {a_{32}}{a_{13}}}&{{a_{11}}{a_{33}} + {a_{31}}{a_{13}}}&{{a_{11}}{a_{32}} + {a_{31}}{a_{12}}}\\ {{a_{11}}{a_{21}}}&{{a_{12}}{a_{22}}}&{{a_{13}}{a_{23}}}&{{a_{12}}{a_{23}} + {a_{22}}{a_{13}}}&{{a_{11}}{a_{23}} + {a_{21}}{a_{13}}}&{{a_{11}}{a_{22}} + {a_{21}}{a_{12}}} \end{array}} \right] $

(13)

将式(11)代入到上述表达式，可得

$ {\boldsymbol{M}} = \\ \left[ {\begin{array}{*{20}{c}} {{{\left( {\cos \alpha \cos \beta } \right)}^2}}&{{{\left( {\sin \alpha \cos \beta } \right)}^2}}&{{{\left( {\sin \beta } \right)}^2}}&{ - \sin \alpha \sin 2\beta }&{ - \cos \alpha \sin 2\beta }&{\sin 2\alpha {{\left( {\cos \beta } \right)}^2}}\\ {{{\left( {\sin \alpha } \right)}^2}}&{{{\left( {\cos \alpha } \right)}^2}}&0&0&0&{ - \sin 2\alpha }\\ {{{\left( {\cos \alpha \sin \beta } \right)}^2}}&{{{\left( {\sin \alpha \sin \beta } \right)}^2}}&{{{\left( {\cos \beta } \right)}^2}}&{\sin \alpha \sin 2\beta }&{\cos \alpha \sin 2\beta }&{\sin 2\alpha {{\left( {\sin \beta } \right)}^2}}\\ { - 0.5\sin 2\alpha \sin \beta }&{0.5\sin 2\alpha \sin \beta }&0&{\cos \alpha \cos \beta }&{ - \sin \alpha \cos \beta }&{\cos 2\alpha \sin \beta }\\ {0.5{{\left( {\cos \alpha } \right)}^2}\sin 2\beta }&{0.5{{\left( {\sin \alpha } \right)}^2}\sin 2\beta }&{ - 0.5\sin 2\beta }&{\sin \alpha \cos 2\beta }&{\cos \alpha \cos 2\beta }&{0.5\sin 2\alpha \sin 2\beta }\\ { - 0.5\sin 2\alpha \cos \beta }&{0.5\sin 2\alpha \cos \beta }&0&{ - \cos\alpha \sin \beta }&{\sin \alpha \sin \beta }&{\cos 2\alpha \cos \beta } \end{array}} \right] $

(14)

观测坐标系及本构坐标系下的刚度矩阵C和C₀，柔度矩阵S和S₀分别满足以下关系

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{C}} = \mathit{\boldsymbol{M}} \cdot {\mathit{\boldsymbol{C}}_0} \cdot {\mathit{\boldsymbol{M}}^{\rm{T}}}\\ \mathit{\boldsymbol{S}} = \mathit{\boldsymbol{N}} \cdot {\mathit{\boldsymbol{S}}_0} \cdot {\mathit{\boldsymbol{N}}^{\rm{T}}} \end{array} \right. $

(15)

式中：

M^T-矩阵M的转置矩阵，无因次；

N^T-矩阵N的转置矩阵，无因次。

根据式(1)、式(2)，可以得到观测坐标系下应力和应变的关系式

$ \sigma = {\boldsymbol{M}} \cdot {{\boldsymbol{C}}_0} \cdot {{\boldsymbol{M}}^{\rm T}} \cdot \varepsilon $

(16)

$ \varepsilon = {\boldsymbol{N}} \cdot {{\boldsymbol{S}}_0} \cdot {{\boldsymbol{N}}^{\rm T}} \cdot \sigma $

(17)

通过变换，能够求得任意方位角和井斜角条件下的柔度矩阵S。利用式(17)，根据杨氏模量和泊松比的定义，即可计算任意方位角度下的杨氏模量和泊松比。

以观测坐标系下$z$轴方向的岩石力学参数为例。根据定义^[35]，杨氏模量为岩石受到单轴应力时，应力相对于应变的变化率。而泊松比的定义则是岩石在单轴受压条件下横向应变与纵向应变之比。在观测坐标系下，假定岩石只受到单轴应力$\sigma_z$的作用，即

$ \sigma = {\left[ {0, 0, {\rm{ }}{\sigma _z}, 0, 0, 0} \right]^{\rm T}} $

(18)

式中：$\sigma_z$-沿$z$轴的单轴应力，GPa。

将式(18)代入式(17)，可以得到观测坐标系下与应力$\sigma$相关的应变$\varepsilon$表达式，假设为

$ \varepsilon = {\left[ {{\varepsilon _{xx}}, {\rm{ }}{\varepsilon _{{yy}}}, {\rm{ }}{\varepsilon _{{zz}}}, 0, 0, 0} \right]^{\rm T}} $

(19)

式中：$\varepsilon_{xx}$-沿$x$轴的应变，无因次；

$\varepsilon_{yy}$-沿$y$轴的应变，无因次；

$\varepsilon_{zz}$-沿$z$轴的应变，无因次。

则根据定义，可以求取该方向的杨氏模量和泊松比

$ \left\{ \begin{array}{l} {E_3} = \frac{{{\sigma _z}}}{{{\varepsilon _{zz}}}}\\ {\mu _{13}} = - \frac{{{\varepsilon _{xx}}}}{{{\varepsilon _{zz}}}}\\ {\mu _{23}} = - \frac{{{\varepsilon _{yy}}}}{{{\varepsilon _{zz}}}} \end{array} \right. $

(20)

若式(17)中的S₀用C₀表示，则式(20)中解出的杨氏模量和泊松比只与C₀、方位角$\alpha$和井斜角$\beta$有关。类似地，可以求取观测坐标系下$x$轴和$y$轴方向的杨氏模量和泊松比。在正交各向异性条件下，观测坐标系中关于C₀、$\alpha$和$\beta$的杨氏模量和泊松比的解析表达式非常复杂，这里不给出解析表达式，下面进行数值分析。

在本构坐标系下，本文测量了页岩的刚度系数：$C_{11}$=31.89GPa、$C_{12}$=10.29GPa、$C_{13}$=9.94GPa、$C_{22}$=28.18GPa、$C_{23}$=10.11GPa、$C_{33}$=24.48GPa、$C_{44}$=10.64GPa、$C_{55}$=10.72GPa、$C_{66}$=10.80GPa。由于测量过程中岩芯露头有限，没有测量到$C_{22}$和$C_{23}$的数据，这里$C_{22}$取$C_{11}$和$C_{33}$的平均值，$C_{23}$取$C_{12}$和$C_{13}$的平均值。利用式(8)、式(9)计算出本构坐标系下页岩的岩石力学参数：$E_1$=26.26GPa，$E_2$=22.63GPa，$E_3$=19.47GPa，$\mu_{12}$=0.258，$\mu_{13}$=0.300，$\mu_{21}$=0.222，$\mu_{23}$=0.323，$\mu_{31}$=0.222，$\mu_{32}$=0.278。为方便表述，取$E_{\rm{max}}$=26.26GPa，$E_{\rm{min}}$=19.47GPa，$\mu_{\rm{max}}$= 0.323，$\mu_{\rm{min}}$=0.222。

方位角和井斜角对岩石力学参数的影响均是关于$\Pi$的周期函数。因此，取方位角$\alpha$和井斜角$\beta$的变化范围均为180°，间隔为15°。这里解出了一个周期内所有方位角、井斜角所对应的岩石力学参数，为节省篇幅，分析了由VTI模型和ORT模型计算的杨氏模量所产生的相对误差。

图 4a为不同观测坐标系下沿$x$轴方向的杨氏模量$E_1$，图 4b为采用ORT模型和VTI模型计算的$E_1$所产生的相对误差。$E_1$在井斜角30°(150°)和方位角30°(150°)时取最大值26.67 GPa，超过$E_{\rm{max}}$，井斜角为90°时取最小值19.47GPa。用VTI模型计算的杨氏模量在方位角30°$\sim$150°和井斜角0°$\sim$60°、120°$\sim$180°所产生的误差较大，最大值超过了15%。

图 5a为不同观测坐标系下沿$y$轴方向的杨氏模量$E_2$，图 5b为采用ORT模型和VTI模型计算的$E_2$所产生的相对误差。$E_2$随着方位角角度的增加呈现增、降、增、降的趋势，并在方位角为60°(120°)时取得最大值26.41GPa，超过$E_{\rm{max}}$，井斜角的变化对$E_2$没有影响。用VTI模型计算的杨氏模量在方位角0°$\sim$60°和120°$\sim$180°所产生的误差较大，最大值超过了15%。

图 6a为不同观测坐标系下沿$z$轴方向的杨氏模量$E_3$，图 6b为采用ORT模型和VTI模型计算$E_3$所产生的相对误差。$E_3$随着方位角和井斜角角度的变化趋势复杂，在井斜角为0°(180°)时取最小值19.47GPa，井斜角为60°(120°)和方位角为30°(150°)时取最大值26.67GPa，超过$E_{\rm{max}}$。VTI模型计算的杨氏模量在方位角45°$\sim$135°和井斜角45°$\sim$135°所产生的误差较大，最大值超过了15%。

图 7为不同坐标系下沿$x$轴方向的两个泊松比，可以看出同一方向的两个泊松比具有大致类似的变化趋势。$\mu_{21}$随井斜角角度的增加呈现增、降的趋势，在井斜角0°(180°)和方位角45°(135°)时取最小值0.193，小于$\mu_{\rm{min}}$，在井斜角30°(150°)和方位角90°时取最大值0.331，大于$\mu_{\rm{max}}$。$\mu_{31}$随井斜角角度的增加呈现降、增、降、增的趋势，在井斜角45°(135°)和方位角0°(180°)时取最小值0.167，小于$\mu_{\rm{min}}$，在井斜角90°和方位角60°(120°)时取最大值0.331，大于$\mu_{\rm{max}}$。

图 8为不同观测坐标系下沿$y$轴方向的两个泊松比，这两个泊松比随方位角和井斜角的变化并无明显的对应关系，但具有类似的变化趋势。$\mu_{12}$随井斜角和方位角角度的变化关系复杂，在井斜角0°(180°)和方位角45°(135°)时取最小值0.193，小于$\mu_{\rm{min}}$，在井斜角0°(180°)和方位角45°(135°)时取最大值0.301。$\mu_{32}$随井斜角和方位角角度的变化复杂，在井斜角90°和方位角45°(135°)时取最小值0.193，在井斜角45°(135°)和方位角0°(180°)时取最大值0.301。

不同观测坐标系下沿$z$轴方向的两个泊松比随井斜角和方位角的变化具有一定的相似趋势，但又有所不同(图 9)。$\mu_{13}$随井斜角的增加呈现降、增、降、增的趋势，在井斜角45°(135°)和方位角0°(180°)时取最小值0.167，小于$\mu_{\rm{min}}$，在井斜角0°(180°)和方位角60°(120°)时取最大值0.331，大于$\mu_{\rm{max}}$。$\mu_{23}$随井斜角的增加呈现降、增的趋势，在井斜角90°和方位角45°(135°)时取最小值0.193，小于$\mu_{\rm{min}}$，在井斜角0°(180°)和方位角30°(150°)时取最大值0.331，大于$\mu_{\rm{max}}$。

4 层状岩石的等效弹性性质

沉积岩在成岩过程中，通常都是两种或两种以上的岩石互相沉积而成。以页岩为例，浅色的粉砂岩和暗色的泥页岩互层(如图 1所示)，韵律交替发育是造成页岩成层分布的主要原因^{[36, 37]}。不同的岩性具有不同的弹性力学性质，研究砂-页岩互层组成的岩石的弹性性质，及其对应的岩石力学特性，为钻井和压裂提供基础数据和指导依据，具有重要的现实意义和工程意义。VTI介质的层状岩石的等效弹性性质可以由Backus理论进行模拟^[38]，对于ORT介质，则可以用Gerrard提出的模型来分析多层正交各向异性岩石的弹性性质^[39]假定岩石是由页岩和砂岩两种岩体互相沉积而成，两种岩体的体积分别为$V_1$和$V_2(V_1+V_2=1)$，并且两种岩石本构坐标系的对称轴互相重合。

该模型需要满足两点假设：(1)与任何单个层的厚度相比，层系统的规模较大；(2)层的厚度和材料特性相对于它们在系统内的各自位置而随机变化。模型的具体表达式如下^{[39, 40]}

$ \left\{ \begin{array}{l} {\mu _{12}} = \dfrac{\zeta }{a};\begin{array}{*{20}{c}} {} \end{array}{\mu _{13}} = \chi - \dfrac{{\lambda \zeta }}{a};\begin{array}{*{20}{c}} {} \end{array}{\mu _{23}} = \lambda - \dfrac{{\chi \zeta }}{b}\\ {E_1} = \dfrac{{ab - {\zeta ^2}}}{a}\begin{array}{*{20}{c}} ; \end{array}{E_2} = \dfrac{{ab - {\zeta ^2}}}{b}\\ \dfrac{1}{{{E_3}}} = \sum\limits_i {\dfrac{{{v_i}}}{{{E_{3i}}}} + \sum\limits_i {\left( {\dfrac{{{\mu _{13}}}}{{{E_1}}} - \dfrac{{{\mu _{13i}}}}{{{E_{1i}}}}} \right){\chi _i} + \sum\limits_i {\left( {\dfrac{{{\mu _{23}}}}{{{E_2}}} - \dfrac{{{\mu _{23i}}}}{{{E_{2i}}}}} \right){\lambda _i}} } } \\ {G_{12}} = \sum\limits_i {{v_i}{G_{12i}}} ;\begin{array}{*{20}{c}} {} \end{array}\dfrac{1}{{{G_{13}}}} = \sum\limits_i {\dfrac{{{v_i}}}{{{G_{13i}}}}} ;\begin{array}{*{20}{c}} {} \end{array}\dfrac{1}{{{G_{23}}}} = \sum\limits_i {\dfrac{{{v_i}}}{{{G_{23i}}}}} \end{array} \right. $

(21)

其中：

$ \begin{array}{l} a = \sum\limits_i {\frac{{{v_i}{E_{2i}}}}{{1 - {\mu _{12i}}{\mu _{21i}}}}} ;\;\;\;\;\;\;\;\;\;b = \sum\limits_i {\frac{{{v_i}{E_{1i}}}}{{1 - {\mu _{12i}}{\mu _{21i}}}}} \\ \zeta = \sum\limits_i {\frac{{{v_i}{E_{2i}}{\mu _{12i}}}}{{1 - {\mu _{12i}}{\mu _{21i}}}}} = \sum\limits_i {\frac{{{v_i}{E_{1i}}{\mu _{21i}}}}{{1 - {\mu _{12i}}{\mu _{21i}}}}} \\ {\chi _i} = \frac{{{v_i}\left( {{\mu _{13i}} + {\mu _{12i}}{\mu _{23i}}} \right)}}{{1 - {\mu _{12i}}{\mu _{21i}}}}\;;\;\;\;\;\;\;\;{\lambda _i} = \frac{{{v_i}\left( {{\mu _{23i}} + {\mu _{13i}}{\mu _{21i}}} \right)}}{{1 - {\mu _{12i}}{\mu _{21i}}}}\\ \chi = \sum\limits_i {{\chi _i}} {\rm{;}}\;\;\;\;\;\;\;\lambda = \sum\limits_i {{\lambda _i}} \end{array} $

5 算例分析

假定砂岩的刚度系数为$C_{11}$=49.27 GPa、$C_{12}$=18.00 GPa、$C_{13}$=16.24 GPa、$C_{22}$=45.38 GPa、$C_{23}$=17.12 GPa、$C_{33}$=41.49 GPa、$C_{44}$=14.94 GPa、$C_{55}$=15.29 GPa、$C_{66}$=15.63 GPa。由式(8)、式(9)计算出本构坐标系下砂岩的岩石弹性参数：$E_1$=39.58 GPa，$E_2$=35.34 GPa，$E_3$=32.91 GPa，$\mu_{12}$= 0.295，$\mu_{13}$=0.270，$\mu_{21}$=0.263，$\mu_{23}$=0.310，$\mu_{31}$= 0.224，$\mu_{32}$=0.288。取砂岩体积$V_2$分别为0，0.1，…，1.0，计算出混合岩石的等效岩石力学参数。

图 10为杨氏模量随砂岩体积含量的变化关系，随着砂岩体积含量的增加，杨氏模量呈现增加的趋势，其中，$E_1$和$E_2$呈线性增加，而$E_3$以一种近线性的方式增加。

图 11为泊松比随砂岩体积含量的变化关系，与杨氏模量不同，泊松比随砂岩体积含量的变化关系更为复杂。$\mu_{12}$和$\mu_{21}$随砂岩体积含量的增加以近似线性的方式增加，$\mu_{13}$和$\mu_{23}$随砂岩体积含量的增加以近似线性的方式降低，$\mu_{31}$和$\mu_{32}$随砂岩体积含量的增加先降低后增加。

为研究不同砂岩含量条件下页岩和砂岩互层的岩石力学参数在观测坐标系下的变化特征，取砂岩体积$V_2$为0.2、0.4、0.6、0.8，利用上述流程，模拟了不同方位角和井斜角条件下的岩石力学参数的分布特征。考虑到$E_2$$\mu_{12}$= $E_1$$\mu_{21}$、$E_3$$\mu_{13}$= $E_1$$\mu_{31}$、$E_2$$\mu_{32}$= $E_3$$\mu_{23}$。本文只给出了3个杨氏模量$E_1$、$E_2$、$E_3$和3个泊松比$\mu_{12}$、$\mu_{23}$、$\mu_{31}$在不同砂岩体积含量下随方位角和井斜角的变化关系。

图 12为不同观测坐标系下杨氏模量$E_1$在不同砂岩体积含量条件下的变化图。可以看出，随着砂岩体积含量的增加，杨氏模量随方位角和井斜角变化的曲面形态基本不变，数值上逐渐增加。图 12a中红色圆圈部分相比于图 12d中的对应部分的杨氏模量更为平缓。

图 13为不同观测坐标系下杨氏模量$E_2$在不同砂岩体积含量条件下的变化图。随着砂岩体积含量的增加，$E_2$随方位角和井斜角变化的曲面形态完全一致，砂岩体积含量越高，$E_2$越大。

图 14为不同观测坐标系下杨氏模量$E_3$在不同砂岩体积含量下的变化图。杨氏模量随方位角和井斜角变化的曲面形态随着砂岩体积含量的增加基本保持不变，数值上逐渐增加。相比于图 14a中红色圆圈部分，图 14d中的对应部分的杨氏模量更为陡峭。

图 15为在不同砂岩体积含量条件下泊松比$\mu_{12}$随方位角和井斜角的变化，不同砂岩体积含量下$\mu_{12}$曲面形态基本一致，数值上略有增加。砂岩体积含量越高，$\mu_{12}$曲面的变化越趋于平缓。图 14a中红色圆圈部分，4个曲面角向下折，而对应的图 15d中4个曲面角向上翘。

图 16为在不同砂岩体积含量条件下泊松比$\mu_{23}$随方位角和井斜角的变化，不同砂岩体积含量下$\mu_{23}$曲面形态基本一致。砂岩体积含量越高，$\mu_{23}$曲面的变化更趋于平缓，并且在方位角90°(图 16a中红色椭圆部分)附近出现较为明显的增加，其余部分则变化不大。

图 17为泊松比$\mu_{31}$在不同砂岩体积含量条件下随方位角和井斜角的变化，不同砂岩体积含量下$\mu_{31}$曲面形态基本一致。砂岩体积含量的增加使得$\mu_{31}$曲面的变化更趋于平缓，这种变化主要体现在相同方位角不同井斜角的剖面上。而相同井斜角下，不同方位角对曲面形态的影响较小。

6 结论

(1) 页岩的杨氏模量随方位角和井斜角的变化趋势复杂，3个杨氏模量$E_1$、$E_2$、$E_3$随方位角和井斜角变化的最大值均超过了本构坐标系下杨氏模量的最大值26.26GPa，通过比较由VTI模型和ORT模型计算的杨氏模量所产生的相对误差，发现误差的最大值均超过了15%，正交各向异性条件下的岩石力学参数不能用VTI模型近似。

(2) 页岩的泊松比随方位角和井斜角的变化并无明显的对应关系，但同一方向的泊松比随方位角和井斜角的变化具有大致相似的变化趋势。在观测坐标系下，所有泊松比的最小值均小于本构坐标系下的最小值0.222，除了$\mu_{12}$、$\mu_{32}$以外，其余的4个泊松比随方位角和井斜角变化的最大值均大于本构坐标系下的最大值0.323。

(3) 在本构坐标系下，随砂岩和页岩互层砂岩体积含量的增加，杨氏模量增加，泊松比$\mu_{12}$和$\mu_{21}$增加，$\mu_{13}$和$\mu_{23}$降低，$\mu_{31}$和$\mu_{32}$则先降低后增加。

(4) 相同井斜角和方位角条件下，砂岩和页岩互层杨氏模量随砂岩含量的增加而有所增加，而泊松比则变化不一。岩石力学参数的曲面形态除了在某些方位角和井斜角有所变化外，基本保持一致。砂岩体积含量越高，则岩石力学参数曲面变化越趋于平缓。