0 引言

非常规油气，特别是页岩油气的勘探开发在美国获得巨大成功，成为改变化石能源短缺现状的重要后备资源，受到社会广泛重视^[1-2]。页岩是具有大量片状或片状层理分布的沉积岩，岩石物理学研究表明，页岩具有不可忽视的各向异性特征(不同方向的速度差异可达35%)^[3-4]。目前开发的大多数页岩油气储层具有明显的各向异性特征^[5]。

页岩油气的开采需要实施水力压裂改造储层。储层压裂产生微地震，受储层强各向异性的影响，这种被动源的微地震波在页岩中传播时，信号的相应特征(包括传播时间、振幅、相位和频率等)可能会发生改变。波场特征信息的变化，会直接影响微地震监测(定位成像、震源机制反演等)的结果^[6-8]，进而影响对储层的评价。因此，定量评估介质各向异性对被动源微地震波场的影响，对非常规油气的勘探开发有着重要意义。Carcione等^[9]研究了介质各向异性对反射震相(PP和PS)振幅的影响，认为在AVO分析中考虑各向异性是非常必要的。被动源微地震监测的研究表明，如果忽略页岩各向异性，会导致微地震事件的定位出现偏差，因此需对各向异性造成的影响以及各向异性参数的反演开展进一步的研究^[7]。Tsvankin等^[10]总结了各向异性介质中地震波数值模拟、数据处理及反演的最新进展，认为各向异性的影响对微地震波(特别是剪切波)的影响是不可忽视的。通过基于VTI声波方程的微地震波场正演研究，孙伟家等^[11]发现，与各向同性介质相比，页岩的强各向异性会导致地震波的旅行时和振幅出现显著偏差。Li等^[12]研究了页岩各向异性引起的反射地震波振幅及相位偏差，定量分析了不同程度各向异性对反射波波场特征的影响。

上述关于页岩各向异性对地震波影响的研究，更多是利用声波方程，或者是利用弹性波方程研究反射波。基于三维弹性波动方程、对被动源微地震波场的研究相对较少。相比声波，弹性波更接近真实的地震波传播特性。本文基于三维介质弹性波动方程，引入Thomsen参数^[13]表征页岩各向异性，并采用交错网格有限差分算法^[14-19]，正演模拟不同各向异性模型、不同震源机制^[20]的微地震波场记录；通过计算不同波场记录的旅行时偏差和振幅相对偏差，定性分析和定量评估页岩各向异性对被动源微地震波场特征的影响。

1 方法原理 1.1 三维各向异性介质微地震波正演理论

一般认为，弹性波模拟更符合实际地层中地震波的传播。三维介质弹性波波动方程的表达式为

$ \rho \frac{{{\partial ^2}{u_i}}}{{\partial {t^2}}} = \frac{\partial }{{\partial {x_j}}}\left( {{c_{ijkl}}\frac{{\partial {u_k}}}{{\partial {x_l}}}} \right) + \rho {F_i} $

(1)

式中：u代表质点的位移场；ρ和F分别表示介质的密度和所受外力；x、t分别表示空间和时间变量；c_ijkl表示4阶Hooke定律中的刚度矩阵元素，描述了应力与应变的对应关系，下标i, j, k, l=1、2、3，分别表示空间的x、y、z方向。由于弹性系数张量的对称性，刚度矩阵中一共只有21个独立变量。如果是VTI介质，弹性刚度矩阵可简化为

$ \left( {\begin{array}{*{20}{c}} {{c_{11}}}&{{c_{12}}}&{{c_{13}}}&0&0&0\\ {{c_{12}}}&{{c_{11}}}&{{c_{13}}}&0&0&0\\ {{c_{13}}}&{{c_{13}}}&{{c_{33}}}&0&0&0\\ 0&0&0&{{c_{44}}}&0&0\\ 0&0&0&0&{{c_{44}}}&0\\ 0&0&0&0&0&{{c_{66}}} \end{array}} \right) $

(2)

式(2)中的弹性参数无明确的物理意义。因此引入Thomsen参数^[13]

$ \left\{ {\begin{array}{*{20}{l}} {\varepsilon = \frac{{{c_{11}} - {c_{33}}}}{{2{c_{33}}}}}\\ {\gamma = \frac{{{c_{66}} - {c_{44}}}}{{2{c_{44}}}}}\\ {\delta = \frac{{{{({c_{13}} + {c_{44}})}^2} - {{({c_{33}} - {c_{44}})}^2}}}{{2{c_{33}}({c_{33}} - {c_{44}})}}} \end{array}} \right. $

(3a)

$ \left\{ {\begin{array}{*{20}{l}} {{V_{{\rm{Pz}}}} = \sqrt {\frac{{{c_{33}}}}{\rho }} }\\ {{V_{{\rm{Sz}}}} = \sqrt {\frac{{{c_{44}}}}{\rho }} }\\ {{V_{{\rm{Pn}}}} = {V_{{\rm{Pz}}}}\sqrt {1 + 2\delta } } \end{array}} \right. $

(3b)

式中：ε表征纵波各向异性强度；γ表征横波各向异性的强度，同时也表征横波分裂的强度；δ为变异系数，表征地层纵波速度非对称轴方向的变化；V_Pz和V_Sz分别代表纵波与横波的相速度；V_Pn表示地层内纵波的NMO速度。由式(3)可得^[18]

$ \left\{ {\begin{array}{*{20}{l}} {{c_{11}} = \rho (1 + 2\varepsilon )V_{{\rm{Pz}}}^2}\\ {{c_{33}} = \rho V_{{\rm{Pz}}}^2}\\ {{c_{44}} = \rho V_{{\rm{Sz}}}^2}\\ {{c_{12}} = \rho V_{{\rm{Pz}}}^2[1 + 2\varepsilon - 2(1 - {f^\prime })(1 + 2\gamma )]}\\ {{c_{13}} = \rho V_{{\rm{Pz}}}^2\sqrt {{f^\prime }({f^\prime } + 2\delta )} - \rho V_{{\rm{Sz}}}^2}\\ {{c_{66}} = \rho (1 + 2\gamma )V_{{\rm{Sz}}}^2 = \frac{1}{2}({c_{11}} - {c_{12}})} \end{array}} \right. $

(4)

其中f′=$1-\frac{V_{\mathrm{Sz}}^{2}}{V_{\mathrm{Pz}}^{2}}=\frac{c_{33}-c_{44}}{c_{33}}$表示横纵波速比。

由式(4)可看出，VTI介质的弹性参数矩阵一共有6个参数，其中有一个(c₆₆)是非独立的，所以只需要5个Thomsen参数就可以完全描述VTI介质。

1.2 三维微地震波交错网格有限差分算法

由式(1)~式(4)建立了VTI介质的弹性波正演理论。具体模拟时，需将波场分量沿时间和空间方向离散化到相应数值网格上。本文采用经典的交错网格有限差分法^{[5, 21-24]}，该方法具有高效率、低内存需求的特点，已被广泛用于勘探地球物理的地震波场模拟研究中。

交错网格有限差分算法原理为，计算一组空间点的压力及另一组空间点的速度^[21]，用网格点上的差商代替偏微分方程的微商建立差分方程，通过解方程组的形式得到波动方程的数值解。交错网格方法的特点是速度独立于应力而更新，因而数值模拟效率更高效^[22]。

弹性波波动方程(式(1))用速度—应力表示为

$ \rho (\mathit{\boldsymbol{x}})\frac{{\partial {v_i}}}{{\partial t}}(x,t) = \sum\limits_{j = 1}^3 {\frac{{\partial {T_{ij}}}}{{\partial {x_j}}}} (x,t) + \rho (x){F_i}(x,t) $

(5)

式中v_i(x, t)、T_ij(x, t)分别表示质点速度和应力张量。VTI介质中应力—应变关系式为

$ \left[ {\begin{array}{*{20}{c}} {\frac{{\partial {T_{11}}}}{{\partial t}}}\\ {\frac{{\partial {T_{22}}}}{{\partial t}}}\\ {\frac{{\partial {T_{33}}}}{{\partial t}}}\\ {\frac{{\partial {T_{23}}}}{{\partial t}}}\\ {\frac{{\partial {T_{13}}}}{{\partial t}}}\\ {\frac{{\partial {T_{12}}}}{{\partial t}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{c_{11}}}&{{c_{12}}}&{{c_{13}}}&0&0&0\\ {{c_{12}}}&{{c_{11}}}&{{c_{13}}}&0&0&0\\ {{c_{13}}}&{{c_{13}}}&{{c_{33}}}&0&0&0\\ 0&0&0&{{c_{44}}}&0&0\\ 0&0&0&0&{{c_{44}}}&0\\ 0&0&0&0&0&{{c_{66}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\frac{{\partial {v_1}}}{{\partial {x_1}}}}\\ {\frac{{\partial {v_2}}}{{\partial {x_2}}}}\\ {\frac{{\partial {v_3}}}{{\partial {x_3}}}}\\ {\frac{{\partial {v_2}}}{{\partial {x_3}}} + \frac{{\partial {v_3}}}{{\partial {x_2}}}}\\ {\frac{{\partial {v_1}}}{{\partial {x_3}}} + \frac{{\partial {v_3}}}{{\partial {x_1}}}}\\ {\frac{{\partial {v_1}}}{{\partial {x_2}}} + \frac{{\partial {v_2}}}{{\partial {x_1}}}} \end{array}} \right] $

(6)

其中刚度矩阵可由式(4)给出。

交错网格有限差分算法中偏导数可基于以下原理展开

$ \begin{array}{l} \frac{{\partial f}}{{\partial a}}(a) \approx \frac{1}{{\Delta a}}\sum\limits_{m = 1}^N {{k_m}} \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left[ {f\left( {a + \frac{{2m - 1}}{2}\Delta a} \right) - f\left( {a - \frac{{2m - 1}}{2}\Delta a} \right)} \right] \end{array} $

(7)

式中：k_m是有限差分系数(通过求解范德蒙矩阵来计算^[23])；m=1, 2, …, N, 差分阶数为2N(在本文的正演模拟中，时间采用二阶，空间使用十阶)。

根据式(7)，式(5)可以离散化为(这里示例给出v₁)

$ \begin{array}{*{20}{l}} \begin{array}{l} {v_1}\left( {{x_1},{x_2},{x_3},t + \Delta t} \right)\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = {v_1}\left( {{x_1},{x_2},{x_3},t} \right) + {F_1}\Delta t + \frac{{\Delta t}}{{\Delta x}}\frac{1}{\rho }\sum\limits_{m = 1}^N {{k_m}} \times \\ \qquad \begin{array}{*{20}{l}} {\left[ {{T_{11}}\left( {{x_1} + \frac{{2m - 1}}{2}\Delta x,{x_2},{x_3},t + \frac{1}{2}\Delta t} \right) - } \right.}\\ {{T_{11}}\left( {{x_1} - \frac{{2m - 1}}{2}\Delta x,{x_2},{x_3},t + \frac{1}{2}\Delta t} \right) + }\\ {{T_{12}}\left( {{x_1},{x_2} + \frac{{2m - 1}}{2}\Delta x,{x_3},t + \frac{1}{2}\Delta t} \right) - }\\ {{T_{12}}\left( {{x_1},{x_2} - \frac{{2m - 1}}{2}\Delta x,{x_3},t + \frac{1}{2}\Delta t} \right) + }\\ {{T_{13}}\left( {{x_1},{x_2},{x_3} + \frac{{2m - 1}}{2}\Delta x,t + \frac{1}{2}\Delta t} \right) - }\\ {\left. {{T_{13}}\left( {{x_1},{x_2},{x_3} - \frac{{2m - 1}}{2}\Delta x,t + \frac{1}{2}\Delta t} \right)} \right]} \end{array} \end{array} \end{array} $

(8)

根据式(7)，式(6)可以离散化为(这里示例给出T₁₁)

$ \begin{array}{l} \begin{array}{*{20}{l}} {{T_{11}}({x_1},{x_2},{x_3},t + \Delta t)}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = {T_{11}}({x_1},{x_2},{x_3},t) + \frac{{\Delta t}}{{\Delta x}}\sum\limits_{m = 1}^N {{k_m}} \times } \end{array}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{c}} {\left\{ {{c_{11}}\left[ {{v_1}\left( {{x_1} + \frac{{2m - 1}}{2}\Delta x,{x_2},{x_3},t + \frac{1}{2}\Delta t} \right) - } \right.} \right.}\\ {\left. {{v_1}\left( {{x_1} - \frac{{2m - 1}}{2}\Delta x,{x_2},{x_3},t + \frac{1}{2}\Delta t} \right)} \right] + }\\ {{c_{12}}\left[ {{v_2}\left( {{x_1},{x_2} + \frac{{2m - 1}}{2}\Delta x,{x_3},t + \frac{1}{2}\Delta t} \right) - } \right.}\\ {\left. {{v_2}\left( {{x_1},{x_2} - \frac{{2m - 1}}{2}\Delta x,{x_3},t + \frac{1}{2}\Delta t} \right)} \right] + }\\ {{c_{13}}\left[ {{v_3}\left( {{x_1},{x_2},{x_3} + \frac{{2m - 1}}{2}\Delta x,t + \frac{1}{2}\Delta t} \right) - } \right.}\\ {\left. {\left. {{v_3}\left( {{x_1},{x_2},{x_3} - \frac{{2m - 1}}{2}\Delta x,t + \frac{1}{2}\Delta t} \right)} \right]} \right\}} \end{array} \end{array} $

(9)

式中的Δx和Δt分别是空间和时间步长。交错网格有限差分算法中，速度场和应力场在时间和空间上是交错的^{[5, 22, 24]}。

本文正演模拟的稳定性条件为

$ \left\{ {\begin{array}{*{20}{l}} {\Delta t \le \frac{1}{{\sqrt 3 }} \cdot \frac{{\Delta x}}{{{v_{\max }}}}}\\ {\Delta x \le \frac{{{\lambda _{\min }}}}{5}} \end{array}} \right. $

(10)

式中：v_max是介质中的最大速度；λ_min是地震波场的最小波长。此外，在处理边界问题时，使用了吸收边界条件，即对边界区域应用以下衰减函数

$ A(i) = {\rm{exp}}\left[ { - d{{\left( {\frac{i}{M}} \right)}^2}} \right] $

(11)

式中：M和i分别为是衰减区域中网格点总数和序号; d是衰减系数(本文中，M=100，d=0.1)。

1.3 被动源微地震波场受各向异性影响的评估

微地震信号易受噪声干扰，实际一般只能清晰识别出直达波震相。据此，本文主要研究各向异性对直达纵波(P波)、直达横波(S波)的影响。考虑到Thomsen参数中ε、δ是纵波各向异性强度的度量，γ表征横波各向异性强度，本文的研究思路为：首先基于均匀及层状介质模型，建立符合页岩特征的不同程度各向异性介质模型，分别分析ε和δ对直达P波、γ对直达S波旅行时和振幅的影响。

计算各向异性模型与各向同性模型地震记录的互相关函数，表征各向异性对微地震波旅行时的影响，即

$ {r_{i,j}}(k) = \frac{{\sum\limits_{n = 1}^N {{D_i}} (n){D_j}(n + k)}}{{\sqrt {\sum\limits_{n = 1}^N {D_i^2} (n)\sum\limits_{n = 1}^N {D_j^2} (n)} }} $

(12)

式中：D_i表示选定的参考道，即各向同性模型的地震波记录；D_j为加入各向异性后的地震记录。假设互相关函数r_{i, j}(k)的最大值对应的时刻为Δt_{i, j}，则表明当D_j平移Δt_{i, j}时，D_i和D_j波形相似度最大，故Δt_{i, j}即为求得的各向异性模型与各向同性模型地震记录之间的旅行时偏差。

表征各向异性对微地震波振幅的影响时，使用每个检波器地震信号的最大振幅(含正负，符号表示地震波的极性)，以及各向异性与各向同性模型微地震信号的最大振幅偏差百分比(Maximum Amplitude Difference Percentage，MADP)

$ {\rm{MADP}} = \left| {\frac{{{A_{{\rm{ani}}}}(n) - {A_{{\rm{iso}}}}(n)}}{{{A_{{\rm{iso}}}}(n)}}} \right| \times 100\% $

(13)

式中：A_ani(n)为各向异性介质模型第n个检波器中的最大振幅(含正负)；A_iso(n)为各向同性介质模型第n个检波器中的最大振幅。MADP可定量描述各向异性对微地震信号振幅的影响程度。

页岩各向异性对被动源微地震波场的影响评估按如下步骤进行：

(1) 根据页岩储层的地质特征，建立均匀各向异性介质模型M₀、M_E、M_D和M_G(具体参数见表 1)，各向异性参数的选择借鉴了其他页岩各向异性研究成果^[25-26]，基于控制变量法分析不同各向异性参数的影响；

(2) 使用交错网格有限差分正演方法，分别模拟模型M₀、M_E、M_D和M_G的地震波场；

(3) 计算模型M_E、M_D和M_G相对于模型M₀的直达P波、直达S波旅行时偏差和最大振幅偏差。

2 数值模拟分析

分别针对均匀和层状VTI介质，建立不同各向异性参数模型，分析不同各向异性对微地震波场特征的影响。

2.1 均匀介质模型

设计一个网格节点数为201×201×201的均匀模型：网格大小为10m×10m×10m，V_P0为3500m/s，V_S0为2020m/s，密度为2.3g/cm³。震源采用主频为20Hz的雷克子波，置于模型内部(1000m, 1000m, 1250m)。为了模拟不同震源机制的影响，设立一组爆炸源(简写为ISO源)及两组剪切源(分别简写为DC1、DC2震源)，矩张量及机制球表现形式见表 2。模型介质各向异性参数同表 1。在地表平面内以震中为中心，沿x、y方向分别布设201个检波器，检波器间距均为10m，沿x方向检波器序号为1~201，沿y方向检波器序号为202~402(图 1)。时间采样间隔为0.5ms，采样时长为1.25s。

2.1.1 各向异性对旅行时的影响

考虑到地震勘探通常使用单分量(一般为z分量)检波器，本文主要分析z分量微地震波场。图 2为t=0.3s时三种震源机制x-z切面的波场快照；图 3为各向异性模型M_E相对于各向同性模型M₀的z分量初至P波旅行时的偏差图形；图 4为不同各向异性模型地面接收波形图。分析图 2~图 4，可以看出：

(1) 各向异性对微地震波旅行时的影响非常显著。在图 3中，旅行时偏差最大为0.015s(30个采样点)，相当于子波周期(0.05s)的30%。

(2) 各向异性影响直达P波在不同方向上的传播。同一震源机制下，不同各向异性介质模型微地震波z分量初至P波相位不同(图 4箭头处)。此外，M_E模型的波场快照(图 2b)中P波在垂直各向异性对称轴方向(x方向)上传播最快，在平行各向异性对称轴方向(z方向)上传播最慢。Thomsen各向异性假设下，纵、横波速度可表示为

$ \left\{ {\begin{array}{*{20}{l}} {{V_{\rm{P}}}(\theta ) = {V_{{\rm{Pz}}}}(1 + \delta {{\sin }^2}\theta {{\cos }^2}\theta + \varepsilon {{\sin }^4}\theta )}\\ {{V_{{\rm{SV}}}}(\theta ) = {V_{{\rm{Sz}}}}\left[ {1 + \frac{{V_{{\rm{Pz}}}^2}}{{V_{{\rm{Sz}}}^2}}(\varepsilon - \delta ){{\sin }^2}\theta {{\cos }^2}\theta } \right]}\\ {{V_{{\rm{SH}}}}(\theta ) = {V_{{\rm{Sz}}}}(1 + \gamma {{\sin }^2}\theta )} \end{array}} \right. $

(14)

式中θ表示波前法向与垂直方向的夹角。沿z方向传播时，θ=0°, V_P(0°)=V_Pz; 沿x方向传播，θ=90°, V_P(90°)=V_Pz(1+ε)。本文基于页岩储层实际情况，M_E模型ε取值为0.25，所以P波在x方向上传播更快。

(3) 各向异性对直达S波的影响与P波不同。从图 2中可以看到，震源机制相同、不同各向异性介质模型(每一行)得到的z分量直达S波，在平行和垂直各向异性对称轴上旅行时相同。通过分析检波器三分量记录的直达S波属性可知，z分量对应SV震相(表 3所示)。对于M₀、M_E、M_D和M_G四种模型，沿z方向传播时，θ=0°, V_SV(0°)=V_Sz；沿x方向传播时，θ=90°, V_SV(90°)=V_Sz，因而在这两个方向上的SV波传播速率相同。

(4) 各向异性对微地震波旅行时的影响与炮检距有关。从图 3可看出，同一震源机制下，随着炮检距的增加，M_E与M₀的地震波旅行时偏差随之增大。由于炮检距增大，地震波的传播路径变长，各向异性的影响累积增加；其次，随着炮检距增大，检波器接收到地震波场的水平分量能量增强，导致总波场的能量分配从垂直分量占绝对优势变为垂直、水平分量均等，考虑到各向异性导致P波水平方向速度大于垂直方向，故由各向异性造成的旅行时偏差增大。

(5) 各向异性对微地震波旅行时的影响与震源机制无关。图 2中，模型相同、震源机制不同时(每一列)，波场旅行时无明显差异。从图 3也可看出震源机制对微地震波旅行时偏差无影响。地震波旅行时仅由传播距离及速度决定，传播距离一定，速度由式(14)确定，因此震源机制对微地震波旅行时无关。

2.1.2 各向异性对振幅的影响

图 2中同一震源机制(每一行)均采用相同增益，通过横向对比可以看出，各向异性模型的微地震波振幅明显不同于各向同性模型。

首先分析各向异性对直达P波振幅的影响。图 5为不同震源(ISO、DC1、DC2)模型M_E、M_D和M₀的z分量初至P波最大振幅及相应偏差(MADP)分布图(沿y方向测线，序列号202~402)，从中可以看出：

(1) ε、δ对微地震P波振幅的影响较为显著。图 5b中，DC1震源激发，M_E和M_D的MADP最高可达60%，说明考虑各向异性影响后，微地震波振幅最高可增大0.6倍。

(2) 从爆炸源的最大绝对振幅(图 5a左)发现，M₀、M_E、M_D三组模型的最大绝对振幅均沿震源左右对称；剪切源的最大绝对振幅(图 5a中、图 5a右)相对于震源中心对称，且极性出现反转。剪切源与爆炸源最大绝对振幅分布的明显不同是由震源机制造成的。此外，图 2中不同方向微地震波的振幅不同，沿各向异性对称轴方向(z方向)振幅最强，垂直方向(x方向)振幅最弱，这是因为沿x方向微地震波传播速度最快，传播路径最远，因而能量衰减最多，振幅最弱。这与波的空间传播效应、各向异性导致的垂直和水平能量分配不一致有关。

(3) 不同各向异性参数对微地震波场的影响程度不同。从图 5b可看出，在近炮检距处(检波器序号250~350之间，对应炮检距约小于1350m)，δ对微地震波振幅影响明显强于ε；远炮检距处(检波器序号202~250和350~402之间，对应炮检距大于1350m)，ε对微地震波振幅的影响相对较强。这表明ε、δ对振幅的影响程度随炮检距发生变化。

(4) 不同震源机制各向异性对微地震波振幅的影响规律可能相似。从图 5b可看出，不同震源机制M_E、M_D相较于M₀的MADP数值不同，表明各向异性对振幅能量的影响存在差异。但是，图 5b中MADP随炮检距的变化曲线形态相似，这是因为MADP反映的是相对振幅偏差，在一定程度上减弱了震源机制的影响。

爆炸源在各向异性介质中不产生S波，故本文仅分析剪切震源各向异性对微地震直达S波振幅的影响。图 6为剪切源M_G和M₀模型的z分量直达S波最大振幅分布图(沿y方向测线)。从图 6可以看出，M_G与M₀的最大绝对振幅图形重合，偏差为零。这是因为直达S波z分量为SV震相(表 2)，V_SV不受γ影响(式(14)可以证明), 即γ对微地震直达S波z分量振幅无影响。

根据表 2可知，y方向检波器的直达S波SH震相分布在水平分量，考虑地震勘探通常使用单分量(一般为z分量)检波器，不常使用水平分量, 故本文暂不分析SH震相。关于γ对直达S波SH震相振幅的影响可作为后续研究方向。

2.2 层状VTI介质模型

设计的层状VTI介质模型如图 7所示。模型网格节点数为201×201×201，网格大小为10m×10m×10m。震源采用主频为20Hz的雷克子波，位于(1000m, 1000m, 1250m)处。模型中间层为各向异性页岩层，各向异性参数同表 1，其他参数见表 4。观测系统同图 1：在地表过震中布设两条垂直交叉测线(分别沿x和y方向，图 7)，每条测线布有201个检波器，间隔10m。基于三组震源(ISO、DC1、DC2震源，表 1)，对四组介质模型分别正演模拟了402个检波器的地震波形记录；时间采样间隔为0.5ms，采样时长为1.25s。

2.2.1 各向异性对地震波场旅行时影响

图 8为t=0.3s时不同震源机制x-z切面的微地震波场快照；图 9为不同震源机制M_E与M₀的z分量初至P波、S波旅行时偏差图；图 10为不同震源、不同各向异性模型的波形记录图。分析图 8~图 10，可以看出：

(1) 图 9中P波旅行时偏差最大为0.012s(相当于24个采样点，达子波周期(0.05s)的24%，S波旅行时偏差最大为0.015s(相当于30个采样点)，达子波周期(0.05s)的30%，表明各向异性对微地震波旅行时的影响较为显著。

(2) 同一震源机制，不同各向异性介质z分量初至P波(图 10红色箭头)、初至S波(图 10蓝色箭头)相位均不同。此外, 从图 8b可以看出，ε>0时，P波在x方向上传播最快，在z方向上传播最慢。

(3) 同一震源机制，不同各向异性介质模型微地震波z分量的初至SV波，在x、z方向上的旅行时相同。

(4) 随着炮检距的增加，地震波的传播路径变长，各向异性对微地震波旅行时的影响累积增加。

(5) 不同的震源机制，微地震波初至旅行时无明显区别(图 8、图 9)，各向异性对微地震波旅行时的影响与震源机制无关。

与均匀介质各向异性相比，层状VTI介质模型各向异性相与各向同性模型P波的旅行时偏差值变小。考虑到震源位置相同，检波器布设在地表，均匀介质(图 1)和层状VTI模型(图 7)两者的差异主要体现在建立各向异性模型的方式(均匀介质为全空间各向异性，VTI介质仅中间页岩层存在各向异性)不同。震源位于模型中心，因此两种模型微地震波的差异主要由第一层(图 7)介质参数(速度、密度、各向异性)的不同而引起。相较于均匀介质，层状VTI介质模型的总体各向异性程度减弱，因此VTI模型P波的旅行时偏差(子波周期的24%)小于均匀介质(子波周期的30%)。

2.2.2 各向异性对地震波场振幅影响

首先考虑ε、δ对微地震直达P波振幅的影响。图 11为不同震源M_E、M_D和M₀模型z分量初至P波的最大绝对振幅值及相应偏差图(沿y方向测线，序列号202~402)。通过分析可得到以下认识：

(1) 各向异性对微地震波振幅有明显影响。从图 11b可看到，DC1震源的MADP最高可达30%。

(2) 爆炸源的最大绝对振幅(图 11a)沿震源左右对称，随炮检距的增加而逐渐减小，在震源处最大；剪切源的最大绝对振幅(图 11a中、11a左)出现极性反转，沿震源中心对称。剪切源与爆炸源最大绝对振幅图形的不同，本文认为这与震源机制有关。

(3) 不同震源，M_E和M_D的MADP变化规律相似(图 11b)，表明各向异性对振幅能量的影响规律相似，这也是由于由式(13)计算的MADP反映相对偏差在一定程度上减弱了震源机制的影响。

(4) 观察同一震源模型M_E和M_D的MADP(例如图 11b左)变化可知，ε、δ对振幅的影响程度随炮检距增加而发生变化。在近炮检距处(检波器序号260~340之间，对应炮检距小于1310m)，δ对微地震波振幅的影响明显强于ε；远炮检距处(检波器序号202~260和340~402之间，对应炮检距大于1300m)，ε对微地震波振幅的影响较强。

(5) 与均匀模型相对比，层状VTI介质模型MADP数值普遍较小，不超过30%。这是因为层状VTI介质模型各向异性仅存在于中间页岩层处，各向异性程度总体减弱，故MADP数值降低。

(6) 相比均匀模型，层状VTI模型ε、δ对振幅影响程度发生反转的位置不同，对应的炮检移由约1350m处减小到约1310m。这是因为检波器位于地表，接收到的两个模型微地震波信号特征的差异主要由第一层介质不同导致。

同前面均匀介质模型相似，DC1、DC2震源直达S波z分量(图 12)M_G与M₀得到的最大绝对振幅图形重合，偏差为零，说明微地震直达SV波(表 3)振幅不受γ影响。

3 结论与建议

本文针对页岩储层各向异性对被动源微地震波场特征的影响展开研究。首先基于三维各向异性介质弹性波动方程和交错网格有限差分算法，正演模拟了被动源的微地震波场。通过对比不同震源、不同程度各向异性介质模型的理论波形记录，计算其旅行时、振幅相对偏差，定量分析了均匀及层状VTI模型中各向异性对微地震波场特征的影响。

根据实际微地震事件的记录特点，本文主要评估了垂直z分量微地震波受各向异性的影响。数值算例的结果表明，对于所采用的地面微地震监测系统及均匀(或平层状VTI)介质模型，得到如下结论：(1)各向异性的微地震记录与各向同性介质的波场记录相比，旅行时和振幅相对偏差均较为显著(旅行时偏差可达子波周期的30%，振幅相对偏差则高达60%)；(2)各向异性对微地震波传播速度具有方向性的影响。加入各向异性后(ε>0)，直达P波在垂直各向异性对称轴的方向上传播最快，平行对称轴方向上传播最慢，直达SV波在垂直、平行各向异性对称轴方向上传播速度相同；(3)炮检距越大，各向异性对波形的影响越显著；(4)不同各向异性参数对波场影响的程度不同，近炮检距处(炮检距约小于1310m)，δ对微地震波振幅影响明显强于ε，远炮检距处(炮检距约大于1310m)，ε对微地震波振幅的影响较强；γ对直达S波z分量(SV震相)的振幅无影响；(5)微地震波场特征随各向异性的变化规律与震源机制无关。

各向异性对地震记录的影响，由地层产状、各向异性对称轴方向及地震波传播路径等因素决定。进行页岩气开采和利用微地震信息进行储层评价时，需要充分考虑各向异性的影响。忽略各向异性，将影响微地震震源机制的正确判断，这直接影响储层评价工作的进行。关于各向异性对被动源微地震波场影响的研究仍存在很多问题，本文建立的模型较为简单，具有倾斜对称轴的TTI模型可能更符合实际情况。此外，各向异性对直达S波SH震相的影响也需进一步研究。本文为定量评估页岩各向异性的影响提供了一种研究方法和思路，在实际微地震传播特征研究中应具有良好应用前景。