0 引言

钻孔中的电学或超声波井壁图像测井提供了分辨率高达5 mm的井壁图像，对于复杂的裂缝型地层，研究人员多是通过观察井壁成像测井资料描述地层信息，如岩层面、天然或钻井诱导裂缝及井壁破裂、坍塌等；在此基础上结合岩心观察建立地质特征解释模式^[1-2]。然而，钻孔成像的探测深度很小，电阻率成像探测深度通常为几厘米，超声波成像数据也只能提供井壁的声幅或到时图像^[3]。因此，很难评估井壁张开裂缝在地层中的分布情况。如当最大水平应力方向与天然裂缝走向近一致时(通常称为“Ⅰ型裂缝”)，钻井诱导裂缝和天然裂缝的井壁图像特征非常相似^[4]，单独使用成像测井通常难以鉴别；此外，通常岩层界面或风化节理的井壁图像特征类似于裂缝^[5]。

声波测井的探测深度比井壁超声波或微电阻率扫描成像(FMI)大得多^[6]，因此可以反映相对远场地层信息^[7-8]。声波测井一般应用弹性各向异性表征地应力场和裂缝特性等^[9-11]，如Esmersoy等^[12]最早研究了由方位弹性各向异性引起的井眼挠曲横波分裂现象。通常由于复杂的地质原因，如岩石的固有各向异性、高角度地层层面、裂缝或不平衡应力挤压等都可以导致弹性波表观各向异性，因此往往难以判断各向异性的成因^[13]。目前业界一般在频率—慢度域分析纵波、横波、斯通利波波形频散推断声波各向异性成因^[14]。Sinha^[9]、Winkler等^[11]认为地层受不均衡地应力挤压导致快、慢横波的慢度频散曲线呈“交叉”特征，固有各向异性(如裂缝介质)导致快、慢横波频散曲线呈“平行”特征。上述理论模拟结果都是在假定井轴与各向异性介质对称轴一致的情况下得到的^[15]。斜井或井眼与地层倾斜相交时，波传播的复杂性常常使“交叉”或“平行”特征混杂、模糊。通过四分量正交偶极子波形数据提取快、慢波阵列波形可反演频散曲线，因此频散曲线的可靠性依赖于快、慢波阵列波形的准确性^[16]。

当声波各向异性成因存在多解性时，有必要建立结合成像测井信息与声波弹性属性的数学物理定量研究模型，利用井眼地质信息模拟和解释裂缝、应力等引起的声波各向异性，进而确认井壁图像信息的真实性及其地质意义。

前人针对单一方位二维裂缝模型进行了声波正演模拟，但不适用于井筒含多组任意方位裂缝介质的弹性波模拟^[17-19]。此外，现有应力—速度线性关系模型存在必须依赖实验建模验证模型合理性并确定耦合系数的应用局限性，而且对于孔隙度较小和裂缝较稀疏的岩石，在无损样本实验条件下，应力—速度耦合效应并不明显。由于实验需记录两个正交方向的声波速度，因此只能在另一个空间方向施加单轴应力，从而限制了应力—速度模型的普适性。

为此，本文首先拾取井壁成像测井中导致声波表观各向异性的三类结构特征(天然裂缝、层面缝(或薄层)和钻井应力诱导缝)信息，同时定量化层理结构信息(如倾向、方位等)，并将裂缝分层、分组；然后从声波全波列中提取地层的快、慢横波慢度及各向异性方位；接着在岩石和裂缝信息井壁成像建模基础上，通过裂缝附加柔度模型、层面伪裂缝模型和非线性应力场弹性波动理论正演模拟岩石的声波属性；最后通过对比声波速度及各向异性方位等信息的模型预测结果与实测结果，利用最小误差拟合反演确定井壁图像表征的裂缝类型，并对比、分析了频散曲线特征。

1 成像测井裂缝提取与分组

高分辨率电成像测井或超声波成像测井通常用于钻孔地质解释，图像分辨率通常为5 mm，探测深度通常为几厘米，具体取决于所用工具。对比、分析电阻率(或声学参数)与图像形态，可识别主要的平面或线性构造特征(层理面、褶皱面、裂缝面和节理)，并测量平面结构的倾角和方位角。本文主要关注在图像特征上容易混淆的天然裂缝、层节理缝及反映地应力特征的钻井诱导缝(图 1左)。天然裂缝中的开启高导缝主要由构造作用形成，图像呈黑色条带状，可拟合为正弦曲线，对于储层的形成和改造具有重要作用，并决定油气的储、渗性质。诱导缝为钻井过程中产生的人工缝，是由钻具振动、应力释放和钻井液压裂等因素形成的，对储层原始储、渗空间没有贡献，图像上往往呈羽状排列的不连续短裂缝或者在对称位置出现的两条长直缝，这组裂缝的走向即为现今最大水平主应力方向；层理或风化节理缝在图像上呈顺岩石层面或产状规律分布的似裂缝暗色条带状，且不穿过岩层，在砂泥岩纹层或花岗岩风化节理地层中较常见。

需要注意的是成像测井和声波全波列测井的纵向分辨率差异较大，前者为毫米级，后者通常为米级。井筒裂缝地层的声波模拟首先在声波仪器长度范围内抽化成像测井地质信息(图 1左)，即先分层段再按方位分组，将某一层段的图像地层特征按方位相近的原则抽化为几组裂缝或地层特征，形成正演模型(图 1中)。既体现了裂缝的地质成因，又可以按裂缝方位分组建立裂缝的附加柔度模型。文中结合网状裂缝说明裂缝分组的方法及效果(图 2)。对处理窗口深度段内的成像测井信息按照相对方位分为1~3组(1组即不分组)。首先根据井壁电阻率图像拾取裂缝方位角(图 2a)，得到网状缝分布(图 2b)；计算全部裂缝方位向量间的点积，搜索点积最大的两个方位作为计算窗口内具有极端方向的主裂缝方位，其他裂缝围绕这两个主裂缝方位分组。通过计算剩余裂缝方位和主裂缝方位之间的点积，将点积较小者归入相应主裂缝方位，分为两组网状缝(图 2c、图 2d)。对三分组而言，采用海伦三角剖分公式计算给定裂缝方位三元组之间的“面积”，搜索其中最大面积的三联线组为三个主裂缝方位。一旦确定主裂缝方位，计算每个非主裂缝方位和主裂缝方位之间的点积，将点积较小者分别归入相应三个主裂缝方位。

2 裂缝地层横波分裂和方位指示

偶极声波测井的主要目的是测量地层的横波各向异性。由构造应力或其他地质因素形成的裂缝型地层的横波速度通常呈方位各向异性，如沿井轴方向传播的横波，当质点振动方向平行于裂缝走向时速度高，当质点振动方向垂直于裂缝走向时速度低。如果横波的质点振动方向与裂缝走向成一角度，则入射横波分裂为质点振动方向平行和垂直于裂缝走向并以不同速度传播的快横波和慢横波——横波分裂现象(图 1右)。值得指出的是，这种横波分裂现象存在于裂缝型地层及地应力不平衡的非裂缝型地层中，即存在于强各向异性地层中；当地层倾角较大且层理发育时会表现出较强的各向异性，可结合频散图及本文研究成果进一步分析。

在各向异性介质中存在三种体波(一种纵波和两种横波)，其对应的偏振方向相互正交，并且沿着Christoffel矩阵的特征向量方向。配备有偶极子源和交叉接收阵列的声波工具可记录平行和垂直于源偶极子轴向的剪切波场。这些剪切波形被输入到Alford旋转算法中，其目标是沿着仪器轴旋转波形，使其旋转至与各向异性介质的剪切波场偏振方向对应。然后分析旋转剪切波形，以确定相应的快、慢横波的慢度及快横波方位(FSA)。为此，本文基于井眼图像获得的裂缝数据及岩石基质的弹性属性建立正演模型(图 1中)及相应的刚度矩阵，然后通过上述方法模拟各模式波的速度及FSA，并与实际声波测井资料处理结果对比，以确定声波弹性各向异性的成因。

3 井壁成像和声波联合模拟 3.1 天然裂缝模拟

采用Hudson裂缝介质等效理论与线性滑动理论模拟裂缝型地层的声弹性。在声波传播过程中，裂缝被看作是一个位移不连续的线性滑动界面。在线性滑动模型中，假定裂缝可以用一个柔性面表示，穿过该界面，由声波引起的位移是不连续的，而应力是连续的。位移矢量的位错与应力矢量之间的线性关系由裂缝切向和法向柔度张量等决定。据Schoenberg等^[20]的裂缝附加柔度理论，将已知背景岩石中包含的裂缝表示为附加应变源，将裂缝产生的柔度张量S_ijkl^f加到背景柔度张量S_ijkl^b上，得到裂隙介质柔度张量

$ \boldsymbol{S}_{i j k l}=\boldsymbol{S}_{i j k l}^{\mathrm{b}}+\boldsymbol{S}_{i j k l}^{\mathrm{f}} $

(1)

$ \boldsymbol{S}_{i j k l}^{\mathrm{f}}=\delta_{i k} \boldsymbol{\alpha}_{j l}+\delta_{i l} \boldsymbol{\alpha}_{j k}+\delta_{j k} \boldsymbol{\alpha}_{i l}+\delta_{j l} \boldsymbol{\alpha}_{i k}+\boldsymbol{\beta}_{i j k l} $

(2)

式中: δ为克罗内克函数；α、β分别为二阶和四阶正交张量，分别表征裂缝的切向和法向柔度，张量下标索引i、j、k、l取值为1~3，即

$ \boldsymbol{\alpha}=\frac{1}{4 V} \sum\limits_r Z_{\mathrm{T}}^{(r)} \boldsymbol{n}_i^{(r)} \boldsymbol{n}_j^{(r)} A^{(r)} $

(3)

$ \boldsymbol{\beta}=\frac{1}{V} \sum\limits_r\left[Z_{\mathrm{N}}^{(r)}-Z_{\mathrm{T}}^{(r)}\right] \boldsymbol{n}_i^{(r)} \boldsymbol{n}_j^{(r)} \boldsymbol{n}_k^{(r)} \boldsymbol{n}_l^{(r)} A^{(r)} $

(4)

式中：n^(r)为第r条裂缝的法向向量；A^(r)为第r条裂缝的面积；V为裂缝所在岩石体积，对于板状平行裂缝该项可省略；Z_N^(r)、Z_T^(r)分别为第r条裂缝的法向、切向柔度，对其调整可使实测和模拟横波慢度之间匹配更好。

对于井壁成像测井获取的任意一组层面倾角为θ、方位角为φ的平行板状裂缝(图 3虚线所示)，在直角坐标系下该组裂缝面的法向向量为

$ \boldsymbol{n}^{(r)}=(\cos \varphi \sin \theta, \sin \varphi \cos \theta, -\cos \theta) $

(5)

对于VTI介质(对称轴沿z方向)弹性刚度矩阵形式最简单，对于具有任意方向对称轴的TTI介质，矩阵形式的复杂性取决于TTI介质的对称轴方向。在实际地层和井眼情况下，可以按前文方法对成像测井拾取的裂缝分组，然后进行裂缝柔度矢量叠加。由于裂缝介质的对称轴和VTI介质刚度矩阵对称轴之间存在一定角度(图 3)，需要分别描述各组裂缝体系、叠加背景并旋转坐标。TTI介质相对于观测系统的方向由θ和φ定义。为了研究xyz系统的波动，必须将矩阵从x'y'z'系统旋转到xyz系统。这里，可以把旋转看作两个单独旋转的乘积。第一个是x'Oy'面旋转角度θ到xOy面，旋转矩阵为

$ \boldsymbol{R}_1=\left(\begin{array}{ccc} \cos \theta & 0 & \sin \theta \\ 0 & 1 & 0 \\ -\sin \theta & 0 & \cos \theta \end{array}\right) $

(6-1)

第一次旋转后，z'轴和z轴重合。第二个是绕z轴以90°－φ的角度顺时针旋转，旋转矩阵为

$ \boldsymbol{R}_2=\left(\begin{array}{ccc} \sin \varphi & \cos \varphi & 0 \\ -\cos \varphi & \sin \varphi & 0 \\ 0 & 0 & 1 \end{array}\right) $

(6-2)

因此，从x'y'z'系统到xyz系统的旋转矩阵为

$ \boldsymbol{R}=\boldsymbol{R}_2 \boldsymbol{R}_1=\left(\begin{array}{ccc} \cos \theta \sin \varphi & \cos \varphi & \sin \theta \sin \varphi \\ -\cos \theta \sin \varphi & \sin \varphi & -\sin \theta \cos \varphi \\ -\sin \theta & 0 & \cos \theta \end{array}\right) $

(6-3)

联合式(1)~式(6)，根据Bond变换^[21]得到观测坐标系的柔度，从而得到一组任意方位裂缝的附加柔度模型，为一个6×6阶的矩阵Z_r=(Z₁ Z₂), 其中

$ Z_1=\left(\begin{array}{ccc} Z_{\mathrm{T}} \sin ^2 \theta_r \sin ^2 \varphi_r & 0 & 0 \\ 0 & Z_{\mathrm{T}} \sin ^2 \theta_r \cos ^2 \varphi_r & 0 \\ 0 & 0 & Z_{\mathrm{N}} \cos ^2 \theta_r \\ 0 & -Z_{\mathrm{T}} \sin \theta_r \cos \theta_r \cos \varphi_r & -Z_{\mathrm{N}} \sin \theta_r \cos \theta_r \cos \varphi_r \\ -Z_{\mathrm{T}} \sin \theta_r \cos \theta_r \sin \varphi_r & 0 & -Z_{\mathrm{N}} \sin \theta_r \cos \theta_r \sin \varphi_r \\ Z_{\mathrm{T}} \sin ^2 \theta_r \sin \varphi_r \cos \varphi_r & Z_{\mathrm{T}} \sin ^2 \theta_r \sin \varphi_r \cos \varphi_r & 0 \end{array}\right) $

(7-1)

$ Z_2=\left(\begin{array}{ccc} 0 & -Z_{\mathrm{T}} \sin \theta_r \cos \theta_r \sin \varphi_r & Z_{\mathrm{T}} \sin ^2 \theta_r \sin \varphi_r \cos \varphi_r \\ -Z_{\mathrm{T}} \sin \theta_r \cos \theta_r \cos \varphi_r &0 &Z_{\mathrm{T}} \sin ^2 \theta_r \sin \varphi_r \cos \varphi_r \\ -Z_{\mathrm{N}} \sin \theta_r \cos \theta_r \cos \varphi_r & -Z_{\mathrm{N}} \sin \theta_r \cos \theta_r \sin \varphi_r & 0 \\ Z_{\mathrm{N}} \sin ^2 \theta_r \cos ^2 \varphi_r+Z_{\mathrm{T}} \cos ^2 \theta_r & Z_{\mathrm{N}} \sin ^2 \theta_r \sin \varphi_r \cos \varphi_r & -Z_{\mathrm{T}} \sin \theta_r \cos \theta_r \sin \varphi_r \\ Z_{\mathrm{N}} \sin ^2 \theta_r \sin \varphi_r \cos \varphi_r & Z_{\mathrm{N}} \sin ^2 \theta_r \cos ^2 \varphi_r+Z_{\mathrm{T}} \cos ^2 \theta_r & -Z_{\mathrm{T}} \sin \theta_r \cos \theta_r \cos \varphi_r \\ -Z_{\mathrm{T}} \sin \theta_r \cos \theta_r \sin \varphi_r & -Z_{\mathrm{T}} \sin \theta_r \cos \theta_r \cos \varphi_r & Z_{\mathrm{T}} \sin ^2 \theta_r \end{array}\right) $

(7-2)

式中：Z_N、Z_T分别为裂缝的法向和切向柔度；下标r代表第r条裂缝。根据Hudson裂缝介质等效理论与线性滑动理论，可以计算Z_N、Z_T，即

$ Z_{\mathrm{N}}=\frac{4 e}{3 g(1-g)\left[1+\frac{\left(k_{\mathrm{f}}+4 \mu_{\mathrm{f}} / 3\right) a}{\pi(1-g) \mu c}\right]} $

(8)

$ Z_{\mathrm{T}}=\frac{16 e}{3(3-2 g)\left[1+\frac{4 \mu_{\mathrm{f}} a}{\pi(3-2 g) \mu c}\right]} $

(9)

式中：e为裂缝密度；k_f、μ_f分别为裂缝充填物的体积模量、切变模量；μ为背景介质的体积模量；a为裂缝长度；c为裂缝宽度；g= v_S²/v_P²为横、纵波速度比的平方。这些信息可由成像和声波测井获取。

含多组任意方位裂缝的三维岩石柔度模型可由式(1)和式(7)得到。刚度矩阵C为柔度矩阵S的逆，即C＝S^－1，它确定了应力与应变之间的关系，下文讨论求解声波属性的方法。

通过将具有不同倾角和方位角裂缝的地层添加到TI背景中，推导TTI介质的弹性矩阵。通过坐标旋转使裂缝角度交替变化，坐标旋转法只能在裂缝方位一致时工作。考虑到TTI介质也是由多期次应力作用挤压产生的，因此本文预先对裂缝按方位相近的原则分组。

由于声波工具的发射器和接收器的探测距离一般为9 m，而由成像测井和斯通利波观察到的裂缝宽度通常为10 μm~2 mm。假设声波频率为1~5 kHz，地层慢度为100~800 μs/ft，则波长约为0.25 ~10 ft(0.08~3.05 m)。因此，低频声波波长总是大于裂缝宽度而小于测井仪器探测范围，现代声波测井技术满足长波长覆盖范围条件。式(1)~式(9)适用于测井声波裂缝模拟。利用时域有限差分方法^[22-23]模拟井眼周围含裂缝地层的偶极声波传播。首先由岩石基质速度计算地层的拉梅常数λ= ρ(v_P²－2v_S²)和μ=ρv_S²，联合刚度矩阵C计算弹性模量，然后进行声波测井模拟计算。图 4为含裂缝地层中偶极子测井模拟的波形数据。由图可以清晰地看到横波速度各向异性，沿裂缝占优方向偏振的横波(红线)速度大于次级裂缝方向偏振的横波(黑线)，二者相差约为8%，这是由模型中的弹性模量及裂缝柔度参数决定的。

通过解释测井图像了解了裂缝性质，正演模拟前必须确定基岩和裂缝的柔度，基岩柔度与介质的弹性性质有关，即无裂缝岩石的弹性参数(可直接由声波测井获取)。实际模拟时可以将基岩视为各向同性，利用高频单极子(不反映各向异性)的纵波时差、低频偶极子的慢横波时差(偏振方向垂直于裂缝面)及测井密度计算初始弹性参数。

将含多组任意方位裂缝的三维岩石柔度模型应用于实际花岗岩测井资料处理，背景介质为各向同性的均质岩石，由FMI观察到的每条裂缝都视为具有相同柔度的裂缝集的一部分。首先，根据FMI地质解释，确定每条裂缝的类型。在模拟迭代过程中，在某一特定深度，通过在深度域开窗(窗长为声波工具长度)，根据裂缝的方位及类型对裂缝分组，将不同组的所有裂缝及基岩的柔度添加到复合柔度张量中。对得到的柔度张量求逆，得到等效刚度张量，由刚度张量求解的声波速度及快横波方位与实际测井资料匹配良好(图 5中第5、第6列)。图 5表明模拟过程中纵波和快、慢横波速度拟合误差一般较小且易于收敛，容易和实测资料匹配，而FSA对裂缝较敏感，主导了模拟匹配优化过程。

3.2 层面裂缝模拟

在电阻率图像上易将层面裂缝或页岩层理面误判为裂缝，其声波测井数据也呈各向异性。一般认为层面各向异性介质是由颗粒定向排列而形成的横观各向同性介质，由电阻率成像测井可识别TI介质的对称轴。因此可以用VTI介质正演模型模拟层状地层和页岩的声波响应。在每个声波计算窗口内对由井壁图像获取的地层倾角和方位角进行平均，平均层面法向可以作为TI介质对称轴。由于横波在各向同性面发生极化，不存在方位各向异性，通常观察到的是仪器的旋转方位，而不是各向异性方位，因此FSA没有意义；斜井的介质或者倾斜层状介质为TTI介质，FSA与层状介质走向平行或垂直。因此声波模拟重点为波速，而不必考虑FSA。本文针对层面缝岩石应力—应变模型(图 6)，利用弹性参数的基本定义推导含层面缝或页岩层理面岩石的基本模量公式。

假设模型由q层地层组成，每层厚度为t_w(w为层号)，杨氏模量为E_w，地层中裂缝介质杨氏模量为E_fw(f代表裂缝，下同)，裂缝剪切模量为G_f，裂缝密度为λ_w。与地层骨架的总厚度$ L=\sum\limits_w t_w$相比，每条裂缝的宽度可以忽略, 每层地层骨架的应变$\varepsilon_w$ $=\frac{\sigma}{E_w}$, 位移$\delta_w=\frac{t_w \sigma}{E_w}, \sigma$为应力。对任一裂缝, 法向位移$\delta_{\mathrm{d}}=\frac{\sigma}{E_{\mathrm{f} w}}$, 因此对于第$w$层的$m$条裂缝, 位移为$\delta_{\mathrm{f}w}=\frac{\lambda_w t_w \sigma}{E_{\mathrm{f}w}}$, 总法向位移$\delta_{\mathrm{N}}=\sum\limits_w \delta_w+\sum\limits_w \delta_{\mathrm{f}w}$是岩层骨架和裂缝引起的位移之和, 由此可得层面缝岩层总应变为

$ \varepsilon=\frac{\delta_{\mathrm{N}}}{L}=\frac{\sigma \sum\limits_w t_w\left(\frac{1}{E_w}+\frac{\lambda_w}{E_{\mathrm{f}w}}\right)}{\sum\limits_w t_w} $

(10-1)

层面缝岩层杨氏模量为

$ E_m=\frac{\sigma}{\varepsilon}=\frac{\sum\limits_w t_w}{\sum\limits_w t_w\left(\frac{1}{E_w}+\frac{\lambda_w}{E_{\mathrm{f}w}}\right)} $

(10-2)

同理，可得层面缝岩层剪切模量为

$ G_m=\frac{\sum\limits_w t_w}{\sum\limits_w t_w\left(\frac{1}{G_w}+\frac{\lambda_w}{G_{\mathrm{f}w}}\right)} $

(10-3)

上述公式中若假设裂缝密度λ_w为零，则变为描述简单层状模型的公式。

3.3 应力诱导裂缝模拟

利用井眼声波测井研究地应力的依据是应力造成的横波各向异性效应。利用基于非线性弹性理论的岩石物理模型，构建了各向异性地层三维应力场函数的有效刚度张量，以预测任意三维应力状态下任意方向的纵、横波运动。因此，该模型克服了现有应力—速度线性关系模型^[24]必须依赖实验建模验证模型正确性和确定耦合系数的应用局限性。

由非线性弹性波动理论模拟应力状态下弹性波的传播^[25]，固体介质的弹性性质依赖于施加的应力，包括外加应力(或称静态应力)和波动造成的扰动压力(或称动态应力)，这种方法采用非线性应力—应变关系。非线性弹性理论提供了有效弹性刚度张量C_ijkl、主应力σ_ik及应变ε_mn之间的关系

$ \left\{\begin{array}{l} \boldsymbol{C}_{i j k l}=\boldsymbol{A}_{i j k l}+\boldsymbol{A}_{i j k l m n} \varepsilon_{m n} \\ \boldsymbol{\varepsilon}_{i j}=\boldsymbol{A}_{i j k l}^{-1} \boldsymbol{\sigma}_{k l} \end{array}\right. $

(11)

式中：A_ijkl为无应力状态的四阶刚度张量，相当于Voigt标记法二阶弹性常数C_ij⁰符号；A_ijklmn为六阶张量，相当于Voigt标记法三阶弹性常数C_ijk符号；ε_ij为剪切应变；重复下标(如m, n)表示爱因斯坦求和约定。

对于固结和中等应力范围岩石，上述非线性弹性模型的刚度和应力之间符合准线性关系。无应力状态的各向同性岩石，施加应力后其弹性特征具正交对称性，可用两个二阶弹性常数(C₃₃和C₅₅)以及三个三阶弹性常数(C₁₁₁、C₁₁₂和C₁₂₃)描述^[26]

$ \left\{\begin{array}{l} C_{11}=C_{33}^0+C_{111} \varepsilon_{11}+C_{112} \varepsilon_{22}+C_{112} \varepsilon_{33} \\ C_{22}=C_{33}^0+C_{112} \varepsilon_{11}+C_{111} \varepsilon_{22}+C_{112} \varepsilon_{33} \\ C_{33}=C_{33}^0+C_{112} \varepsilon_{11}+C_{112} \varepsilon_{22}+C_{111} \varepsilon_{33} \\ C_{23}=C_{12}^0+C_{123} \varepsilon_{11}+C_{112} \varepsilon_{22}+C_{112} \varepsilon_{33} \\ C_{13}=C_{12}^0+C_{112} \varepsilon_{11}+C_{123} \varepsilon_{22}+C_{122} \varepsilon_{33} \\ C_{12}=C_{12}^0+C_{112} \varepsilon_{11}+C_{112} \varepsilon_{22}+C_{123} \varepsilon_{33} \\ C_{44}=C_{55}^0+C_{144} \varepsilon_{11}+C_{155} \varepsilon_{22}+C_{155} \varepsilon_{33} \\ C_{55}=C_{55}^0+C_{155} \varepsilon_{11}+C_{144} \varepsilon_{22}+C_{155} \varepsilon_{33} \\ C_{66}=C_{55}^0+C_{155} \varepsilon_{11}+C_{155} \varepsilon_{22}+C_{144} \varepsilon_{33} \end{array}\right. $

(12)

式中：C₁₁₁、C₁₁₂和C₁₂₃为独立三阶弹性参数，C₁₄₄=(C₁₁₂－C₁₂₃)/2，C₁₅₅=(C₁₁₁－C₁₁₂)/4；C₃₃⁰和C₅₅⁰为独立二阶弹性参数(无应力状态)，因为C₁₂⁰=C₃₃⁰－2C₅₅⁰，又由于C₃₃⁰=λ+2μ，C₅₅⁰=μ，所以C₁₂⁰=λ。

纵波速度(与C₁₁、C₂₂和C₃₃相关)的应力敏感性由非线性系数C₁₁₁、C₁₁₂控制，而横波速度(与C₄₄、C₅₅和C₆₆相关)的应力敏感性由非线性系数C₁₄₄、C₁₅₅控制。

对于无应力状态的各向同性岩石，根据虎克定律，有

$ \boldsymbol{\sigma}_{i j}=\lambda \boldsymbol{\delta}_{i j} \boldsymbol{\varepsilon}_{a a}+2 \mu \boldsymbol{\varepsilon}_{i j} $

(13)

因此

$ \boldsymbol{\varepsilon}_{i j}=\frac{1}{E}\left[(1+\nu) \boldsymbol{\sigma}_{i j}-\nu \boldsymbol{\delta}_{i j} \boldsymbol{\sigma}_{a a}\right] $

(14)

式中：σ_aa与ε_aa分别为体应力与应变；δ_ij为克罗内克函数；$ E=2 \mu(1+\nu)=C_{55}^0\left(3 C_{33}^0-4 C_{55}^0\right) /\left(C_{33}^0-C_{55}^0\right)$为杨氏模量，$\nu=\frac{\lambda}{2(\lambda+\mu)}=\left(C_{33}^0-2 C_{55}^0\right) /\left[2 \left(C_{33}^0-\right.\right. $$ \left.\left(C_{55}^0\right)\right]$为泊松比。

将式(14)代入式(12)，得

$ \left\{\begin{array}{l} C_{11}-C_{33}=\left(C_{111}-C_{112}\right) \frac{1+\nu}{E}\left(\sigma_{11}-\sigma_{33}\right) \\ C_{22}-C_{33}=\left(C_{111}-C_{112}\right) \frac{1+\nu}{E}\left(\sigma_{22}-\sigma_{33}\right) \\ C_{44}-C_{66}=\left(C_{144}-C_{155}\right) \frac{1+\nu}{E}\left(\sigma_{11}-\sigma_{33}\right) \\ C_{55}-C_{66}=\left(C_{144}-C_{155}\right) \frac{1+\nu}{E}\left(\sigma_{22}-\sigma_{33}\right) \end{array}\right. $

(15)

式中σ₁₁、σ₂₂、σ₃₃分别为x、y、z方向的正应力。

弹性张量矩阵中的其他元素(如C₁₁、C₁₂等)也受裂隙影响，但在偶极子声波测井中对横波的影响不大。与声波测井横波各向异性研究密切相关的主要是C₄₄、C₅₅、C₆₆。式(15)中与横波相关的未知数为C₁₄₄、C₁₅₅、C₄₄、C₆₆。为此，还需要通过基尔施方程描述钻孔附近地层应力分布及弹性模量分布，再利用扰动理论建立三个剪切模量的差分方程，即

$ \begin{aligned} \Delta C_{55}(d, \theta)= & {\left[C_{55}^0-\nu C_{144}+(1-\nu) C_{155}\right] \frac{\Delta \sigma_{11}}{E}+} \\ & {\left[C_{144}-(1+2 \nu) C_{55}^0-2 \nu C_{155}\right] \frac{\Delta \sigma_{22}}{E}+} \\ & {\left[2 C_{55}^0(1+\nu)+C_{55}^0-\nu C_{144}+\right.} \\ & \left.(1-\nu) C_{155}\right] \frac{\Delta \sigma_{33}}{E} \end{aligned} $

(16-1)

$ \begin{aligned} \Delta C_{44}(d, \theta)= & {\left[-(1+2 \nu) C_{44}^0+C_{144}-2 \nu C_{155}\right] \frac{\Delta \sigma_{11}}{E}+} \\ & {\left[-\nu C_{144}+C_{44}^0+(1-\nu) C_{155}\right] \frac{\Delta \sigma_{22}}{E}+} \\ & {\left[2 C_{44}(1+\nu)+C_{44}^0-\nu C_{144}+\right.} \\ & \left.(1-\nu) C_{155}\right] \frac{\Delta \sigma_{33}}{E} \end{aligned} $

(16-2)

$ \begin{aligned} \Delta C_{66}(d, \theta)= & {\left[(2+\nu) C_{66}^0-\nu C_{144}+(1-\nu) C_{155}\right] \times } \\ & \frac{\Delta \sigma_{11}+\Delta \sigma_{22}}{E}+\left[C_{144}-(1+2 \nu) C_{66}^0-\right. \\ & \left.2 \nu C_{155}\right] \frac{\Delta \sigma_{33}}{E} \end{aligned} $

(16-3)

式中：d为径向距离；Δσ₁₁、Δσ₂₂、Δσ₃₃为正应力的扰动量；ΔC₅₅、ΔC₄₄由偶极子测井快、慢横波径向速度反演剖面求取^[27-29]；ΔC₆₆由斯通利波频散反演横波获取(模拟偶极子快、慢横波不需要C₆₆)；C₆₆¹⁰为无外加应力岩石切变模量。至此，方程有完备解析解，可计算井眼周围诱导缝地层弹性参数。

图 7为根据上述理论模拟的井壁附近径向和环向应力诱发的声波速度分布。由图可见：由于井壁附近应力集中效应及应力对波速的影响，使井壁附近的波速变化很大，沿最大主应力方向横波速度较低，而在距井壁2~3倍半径处速度明显升高；沿最小主应力方向井壁附近横波速度较高，随着距井壁的距离增加速度降低，受不均衡地应力影响，在距井壁较远处沿最大主应力方向的横波速度大于沿速度最小主应力方向。因此，应力诱发的各向异性速度在径向和环向均发生变化，而延伸较远的天然裂缝以及薄互层等导致的层面各向异性只和方位有关。利用上述特征可更好地指示各向异性成因及井壁裂缝类型。

4 声波属性计算与裂缝类型识别

在上述三类裂缝性地层(天然裂缝、层面伪裂缝和井壁诱导缝)声波弹性模型基础上，计算声波属性与反演裂缝类型修正成像测井人工观察裂缝类型的初始认识。

各向同性介质的应力—应变原理、运动微分方程均适用于各向异性介质，不同的是应力—应变关系。一般情况下由广义虎克定律表达应力—应变关系，其中刚性系数因介质性质不同而变化。上述三个模型均产生各向异性刚度张量，将广义虎克定律代入弹性介质运动微分方程可以得到以位移表示的运动微分方程

$ \rho \ddot{\boldsymbol{u}}_j=\boldsymbol{C}_{i j k l} \boldsymbol{u}_{k, j l} $

(17)

式中：ρ为介质密度；ü为对位移u求偏微分；C_ijkl为刚度张量，下标i、j取值为1~3，重复下标k、l表示求和。

在计算窗口内用前述刚度模型构造Christoffel张量，给定沿钻孔轴线的传播方向，求解Christoffel张量的特征值和特征向量，得到三种弹性波传播模式(准纵波qP和准剪切波qSH、qSV^[30])的极化矢量。在与钻孔正交的平面内，分析剪切波速的方位变化得到反映声波各向异性的FSA以及快速、慢速剪切波慢度。

文中对比了测量和模拟的FSA、声波慢度计算纵波时差(DTCO)、快横波时差(DTSM_Fast)、快慢横波时差的差值(DTSM)及FSA测量值方差的倒数。方差值的设置考虑了各物理量的取值范围及误差权重协调性，通过在模拟过程中设置DTCO和DTSM_Fast的方差，使其在优化过程中的影响最小化(即这两个值很大)，因此DTSM及FSA误差主导了优化过程。在每一个裂缝发育深度段，通过最小二乘误差确定该层段的裂缝类型，进而判断导致声波各向异性的成因机理。

5 应用实例

将本文方法应用于花岗岩潜山裂缝性气田生产测试井以验证方法合理性，岩心岩性为层状花岗岩、裂缝型花岗岩。层状花岗岩在FMI图像上呈明显的成层性，竖直诱导缝及羽状诱导缝发育，特别是在FMI图像上常见暗色条带，极易和天然裂缝混淆；裂缝型花岗岩在FMI图像上可见多种裂缝特征(图 8)。正交偶极子声波测井数据的横波各向异性表明，局部井段各向异性较强，各向异性受层界面、应力及裂缝的共同影响。利用声波测井可分析快、慢横波的各向异性(图 8中第5、第9、第11列)，并结合FMI模拟井壁地质现象，综合分析成像地质特征和声波各向异性的成因，进而评价裂缝有效性。可见：①在2960~2992 m井段，FMI图像疑似反映了层面裂缝和钻井诱导缝，声波测井数据呈中等各向异性，难以判断裂缝的有效性；综合模拟表明，声波FSA与层界面走向总体接近，与基于层面模型和应力模型正演模拟的FSA结果较接近，因此判定此段的顺层面发育裂缝为无效层面缝，声波各向异性主要由地应力引起，且快、慢横波频散曲线呈“交叉”特征。②在2940~2950 m井段，FMI图像存在疑似裂缝特征，FSA总体与裂缝模型正演结果较一致；频散特征表明，低频段快、慢横波分离且呈平行状，因此认为本段成像裂缝为开启有效缝。③在2905~2940 m井段，FMI图像反映了网状裂缝及层状花岗岩特征，FSA介于层面模型和裂缝模型正演结果之间，因此层界面和天然裂缝为此段各向异性成因。④根据不同井段的测井数据分析，可以自动划分近井裂缝各向异性成因分类标识(图 8中第12列)。B井地层分上、下两层，下层(2950~3000 m)发育层面伪裂缝、钻井诱导缝及少量有效缝，测试产能较低，储层无阻流量为3620 m³/d；上层(2900~2940 m)有效缝发育，测试产能较高，储层无阻流量为23950 m³/d，验证了研究结果的合理性。

6 结论

本文通过将井眼图像获得的裂缝等信息与声波各向异性有效整合，创建并分析了不同弹性各向异性模型(天然裂缝(裂缝各向异性)、诱导缝(应力各向异性)和层理缝(层面各向异性))及其特征，利用裂缝附加柔度模型推导了含多组不同方位裂缝地层的刚度矩阵，并结合应力声场模型联合正、反演技术评估裂缝类型及成因。实际模拟结果表明，波速对各向异性成因的敏感度不如快横波方位，波速拟合误差一般较小且易于收敛，而快横波方位更能区分和判断地层的各向异性类型，从而有效、可靠地指示裂缝类型。