1 引言

众多研究者已经对微地震监测开展了广泛的研究，通过分析处理水力压裂过程中产生的微震信号，我们可以反演得到微地震发生的震源位置、震级大小、震源机制等信息，从而得到储层中裂缝的发育情况.这些信息为最终的储层改造评价提供了有力证据(Eisner et al.，2010，2011; Kendall et al.，2011; Maxwell，2011).

Backus和Mulcahy(1976a，1976b)最早对微地震的震源机制进行了定义，他们认为地层有限空间中的应力不均衡导致了岩石的破裂，从而产生了微地震现象；并且，对于不同的应力分布情况，将会产生不同的破裂现象.对于水力压裂中的微地震震源机制反演和储层改造评价也有众多的研究者进行了详细的分析和研究(Julian et al.，1998; Eaton，2009; Baig and Urbancic，2010). 根据所观测微震数据的类型，微地震震源机制反演方法主要分为两大类：基于井中观测数据的震源机制反演方法和基于地表观测数据的震源机制反演方法.但是，井中观测数据容易受到观测方位不足，观测数据类型有限和各类观测噪声的影响，其对应的震源机制反演精度有限(Vavryčuk，2001，2007; Eaton and Forouhideh，2010，2011).为此，Song和Toksz(2011)，Song等(2014)和Grechka(2015)等人分别提出了利用微地震近场数据和介质的各向异性特征，进行单井观测时的震矩张量反演工作.另一方面，对于地表监测的低信噪比微地震数据，利用基于偏移理论的反演方法可以有效地压制各类噪声的干扰，并同时进行微地震定位和震源机制反演工作(Gharti et al.，2011; Dando et al.，2014).然而，观测数据误差所导致的震矩张量反演结果误差，会引起岩石破裂面的不同属性参数解释结果出现不同程度的偏差(Du and Warpinski，2011).Grechka和Blias(2013)在各向同性介质的假设条件下通过数值计算与实际数据分析发现地质模型的速度误差对于微震源定位的影响并不严重，但是对于地层衰减所带来的影响，仍需要做进一步的分析.为此，Eaton等(2014)利用水力压裂过程中的不同射孔记录计算了被压裂地层的Q值，从而对微地震数据进行了吸收衰减补偿，并获得了岩石破裂类型与微震源频谱特征间的关系.随着工业界对水力压裂技术的要求不断增加，Leaney等(2011)提出了基于VTI介质的震源机制反演方法.Das和Zoback(2013)通过反演较长持续时间的微地震震源机制分析了裂缝断层类型与震源机制间的相互关系.Urbancic和Baig(2013)利用不同方位角所记录到的微震信息，反演每个微震事件对应的破裂类型、震源主方向和应力释放模式，从而解释反演水力压裂所形成的裂缝网络特征.杨心超等(2015)基于网格点搜索算法计算了水力压裂过程中的P波初动震源机制解，进而反演得到了裂缝的方位参数.

尽管前人提出了不同的震源机制反演方法并对反演过程中存在的一些问题进行了分析，但地层的吸收衰减作用对震源机制反演的影响并没有得到充分的研究.并且微震信号具有主频高，能量小的特性，地层的吸收衰减必然会严重影响微震信号的波形、频谱和信噪比，从而降低震源机制的反演精度.因此，本文首先根据地层的吸收衰减特性建立了含衰减地层的微地震震源机制频率域反演方法.然后，通过数值计算，并利用反演分辨率矩阵分析了地层吸收衰减的计算误差对不同震源机制反演结果的影响.

2 震矩张量与微地震波场关系

根据Backus和Mulcahy(1976a，1976b)最早对震源机制的定义，任何震源都是由于地层中的过量应力 σ 造成的，通过对应力过量变化率进行时间与空间上的积分，可以用地震矩对震源进行定义：

(1)

并且，这一严格定义同样可以用震矩张量 M 的形式进行表示：

(2)

其中，M₀表示地震矩幅度，每一个元素M_ij表示沿着±i方向，大小相等，并且在j方向相距无限小的一对力所形成的力偶(i，j=1，2，3).由于任何地震的震矩张量必须满足净应力为零和净应力矩也为零的基本条件(Udías et al.，2014).因此，震矩张量 M 中的元素必须满足M_ij=M_ji这一条件，即，M 是一个对称矩阵，只有6个独立的元素(M₁₁、M₁₂、M₁₃、M₂₂、M₂₃、M₃₃).

利用震矩张量，我们可以对岩石破裂的几种基本形式进行定量的描述.例如：

爆炸型震源(ISO)

补偿线性偶极子震源(CLVD)

线性偶极子震源(LVD)

双偶极子震源(DC)

以上定义中，只是部分震源机制的表示形式.对于同一类型的震源，由于所受的力偶的不唯一性，其震矩张量的表示形式也有多样.如双偶极子震源还可以表示为：

但是，当实际地层中的岩石发生破裂时，其诱发的微地震震源是十分复杂的.利用正则化分解，可以将任意的震源机制分解为三种基本震源机制(ISO、 CLVD、DC震源)的组合(Julian et al.，1998; Chapman and Leaney，2012)，即

(3)

Aki和Richards(2002)将各向同性介质中，点震源的远场波场函数与震源机制的关系定义为：

(4)

其中，u_i(x，t)表示在位置 x 处，t时刻的波场，表示震源时间函数的时间偏导数，ρ表示地层密度，α与β分别是地层的纵波和横波速度，r是震源到检波点的距离，τ^P和τ^S分别表示纵横波从炮点到检波点的走时；R_i^P和S_i^P分别表示纵波和横波的辐射模式函数，利用爱因斯坦求和约定，可以简写为：

(5)

(6)

其中，向量 n =(n_i，n_j，n_k)表示震源到检波点的方向向量，δ_ij为克罗内克函数，M_jk表示震矩张量中的元素.

3 地层的吸收衰减效应

地震波在真实地层中传播的过程中，其能量不仅会发生几何扩散，还会由于地层介质颗粒间的内摩擦作用，孔隙流体的流动，黏弹性介质的黏性张弛等作用而产生吸收衰减和地震波频散现象.而对于地层的吸收衰减效应主要使用地层的品质因子Q值进行表征，其严格定义可表示为：

(7)

其中，E₀表示地震波的原始能量，ΔE表示地震波传播一个周期后所消散的能量(Tonn，1991).在常规的地震勘探中，由于地震子波的频带范围有限，我们假设地层的吸收衰减作用与地震波的频率无关，即地层的Q值在常规地震勘探频带范围内保持不变(Zhang and Ulrych，2002).因此，地层对地震波的衰减过程可以表示为：

(8)

其中，A(f)表示检波器接收到的波形频谱，S(f)表示震源的频谱特征，t表示炮点到检波点的走时.

但是，微地震信号的主频值和频带范围都远远高于常规地震勘探的地震信号特征.并且地层Q值也表现出频率依赖性，对于不同类型的介质和裂缝分布特征，地层的Q_P和Q_S在不同频带范围内具有较大的差异，且随着频率具有不同的变化特征(Pride et al.，2004; Brajanovski et al.，2006; Müller et al.，2010; Quintal et al.，2014).实际的应用中，为了分析和计算的方便，可以假设在微地震事件的有效频带范围内，地层的Q值仍保持其频率独立性(Eaton et al.，2014).因此，在反演的过程中本文仅利用微地震波形有效频带内的数据作为反演过程的输入数据.根据公式(4)—(6)和公式(8)，我们可以得到，简化后的受到地层吸收衰减作用影响的波场函数表达式：

(9)

根据(9)式，本文正演出主能量频带范围受到地层吸收衰减作用影响的微地震波场，以此模拟微地震事件在地层中的传播过程.

4 微震源机制反演与解释

在均匀各向同性介质中，根据公式(9)的定义，震源机制与地震波场函数之间的关系可以用矩阵的形式表示为：

(10)

其中，A_i^W表示检波器观测到的直达P波和S波数据，N表示高斯随机噪声，G_ij，k^W对应于P波和S波，可分别定义为：

(11)

(12)

另一方面，根据波场数据与震矩张量之间的关系可知，公式(10)是一个线性不适定方程组.为了求解该方程组，我们利用MATLAB中的最小二乘算法，即，通过求取观测数据与正演数据间的最小平方差 Γ(公式(13))，便可以反演得到微地震震源机制M_jk.

(13)

在实际的计算过程中，本文通过傅里叶变换将时间域数据转换到频率域进行计算，以提高反演的速度.

对于反演得到的震矩张量，我们将采用Hudson图示法的形式进行解释分析.这一图示法最早由Hudson等(1989)提出，对于任意的震矩张量进行正则化分解，首先计算其特征值并按照大小进行排序：

(14)

然后，定义各向同性组分M及其对应的偏差矩M′ ₁、M′ ₂和M′ ₃分别为：

(15)

(16)

根据M′ ₃的大小，可以定义张性组分k为：

(17)

剪切组分T为：

(18)

通过这一参数化运算，震矩张量的6个独立元素可以简化为两个元素(T、k)进行表示(如图 1).其中，横坐标表示剪切组分T值，纵坐标表示张性组分k值，T、k的取值范围均为-1~1.通过该图我们可以清楚直观的对水力压裂过程中发生的全部微地震震源机制进行统计性的分类描述.

不同的震源机制位于图 1中不同的T-k坐标位置处.例如，各向同性爆炸型震源张量(M _ISO)中的膨胀震源(E)和压缩震源(I)的T-k坐标分别为(0，1)和(0，-1)，而纯剪切型震源张量(M _DC)，即双偶极子震源的T-k坐标为(0，0)，位于Hudson图示的中心位置.两类错动方向相反的补偿线性偶极子震源(M _+CLVD与 M _-CLVD)的T-k坐标则分别为(-1，0)和(1，0).

5 数值计算分析

为了分析地层的吸收衰减作用对微震震源机制反演的影响，本文设计了一个均匀各向同性的地质 模型.三口垂直观测井位于距离震源点距离为150 m 的圆弧上，每口井中模拟9个三分量检波器进行观测(如图 2).在这一观测系统下，所获得的数据具有充足的观测方位角，可以避免因观测方位角不足而带来的反演误差.另一方面，该模型的各项物性参数如表 1所示，其中，纵横波速度与地层密度参考延长组地层的速度和密度特征，地层的Q值则参考前人的研究成果进行设定，在本文的所有数值计算中均采用600 Hz的雷克子波作为震源子波；考虑到实际地层Q值的复杂性，模型中对其设定了一个较大的误差范围-90%~165%.本文主要分析在不同 地层Q值误差条件下，微震源机制反演结果变化情况.

根据图 2的地质模型与观测系统，我们利用DC型震源进行数值计算，得到了直达P波和S波的三分量微地震记录(如图 3).由于DC型震源机制具有特定能量辐射样式，三口井中所观测到的地震记录具有明显的差异，这与常规地震勘探中利用假设爆炸震源所获得的地震记录有本质的不同.

在数值计算过程中，本文利用不同类型微震源正演得到的微地震数据作为观测数据，然后根据公 式(10)反演对应的微震源机制.对于不同的观测数据，反演结果的精度也不尽相同，这是根据不适定方程(10)中的系数矩阵的条件数决定的(常旭等，1999).为分析不同的Q值误差对微震源机制反演的影响，本文采用了Grechka(2010，2015)提出的 费雷谢偏导矩阵SVD分解方法，对震矩张量 M 每个元素的反演分辨率进行了详细分析.费雷谢偏导矩阵SVD分析方法首先求取观测数据关于震矩张量6个独立元素M_ij的偏导数矩阵 F，通过对矩阵 F 进行SVD分解：

(19)

可以得到分解后的奇异值矩阵 Λ和特征向量矩阵 W，对比分析这两个矩阵我们可以获得不同观测数据对于震矩张量 M 的6个独立元素的反演分辨率(如图 4).其中，图 4a和图 4b表示在地层Q值没有误差的情况下，分别利用P波和S波联合反演和单独利用P波反演的震源机制反演分辨率矩阵.

图 4中的红色坐标点表示经过SVD分解后得到的奇异值Λ_i(i=1~6)，灰度网格代表SVD分解后的特征矩阵 W，每个奇异值Λ_i对应特征矩阵 W 中的每一列特征向量；其中，灰度网格的每个方块颜色对应于数值0~1的变化，当W_ij=0时对应白色， W_ij =1时对应黑色.理论上，当奇异值矩阵中的每个元素都是非零值时，可以反演得到唯一的震矩张量矩阵 M _ij. 但是，奇异值矩阵元素 Λ _i的相对大小与特征矩阵 W 中每一列特征向量的稀疏性同样影响着对应震矩张量M_ij的反演精度.例如，在奇异值矩阵元素 Λ _i均大于零时，若部分元素 Λ _k(k⊂i)远远小于其他元素，则元素 Λ _k所对应未知参数的反演结果很容易受到输入数据误差的影响，使整个系统的反演结果不稳定.另一方面，特征矩阵 W 对于反演结果的影响主要体现在每一列特征向量的稀疏性，特征向量越稀疏则对应未知数的约束条件越好，在图 4中主要表现为每一列的灰度网格中仅有一个方块趋近于黑色，其余5个方块的颜色趋近于白色.

当获得的地层Q值没有误差，且观测方位角充足时，无论利用直达P波和S波数据联合反演，还是单独利用直达P波进行反演，均可以得到准确的反演结果.但是，不同类型的输入数据对于震矩张量的不同元素的反演敏感度是不相同的.从图 4中可以看出，无论使用何种类型数据进行反演，经过费雷谢偏导矩阵SVD分解后，其对应的奇异值Λ_i均远大于零，这验证了当Q值没有误差时可以通过当前观测系统反演得到震矩张量.但是，对比图 4a和图 4b中的特征矩阵 W，当利用P波和S波数据联合反演时，待反演震矩张量中的M₁₁、M₁₂和M₁₃所对应的特征向量具有更好的稀疏性；而只利用直达P波进行反演时，输入数据的约束效果刚好相反，即M₂₂、M₂₃和M₃₃所对应的特征向量具有更好的稀疏性.

然而，在实际的微地震反演过程中，模型的各项参数总是存在不同幅度的误差.尤其是地层的Q值计算误差，对于具有高主频的微地震信号更容易产生较大的干扰.因此，本文在数值计算的反演过程中，假定地层的品质因子受到不同程度的误差干扰，而震源子波和地层的其他物性参数均不含有误差，以此重点分析地层不同幅度的Q值衰减误差对震源机制反演结果的影响程度.

5.1 随机微震源机制反演

在本小节中，我们利用423个不同类型的随机微震源作为输入的微地震震源(如图 5)，并根据表 1中的参数所示，将地层的品质因子误差引入反演过程中.然后，采用直达P波和S波数据联合反演得到了受到不同Q值误差影响后的微震震源机制反演结果.

图 6展示了当地层Q值误差从-90%逐渐增 大到165%的过程中，随机微震源的部分反演结果.

从图中可以看出地层的品质因子误差对震源机制的反演结果产生了严重干扰.并且，当地层的品质因子具有正向误差时，随机微震源的反演结果出现了向ISO型震源汇聚的现象；而当地层的品质因子具有 负向误差时，随机微震源反演结果开始向TC型震源汇聚.由此可以说明，当我们过高地估计地层的品质因子时，最终的反演结果更偏向于爆炸型震源，而当我们低估地层的品质因子时，最终的反演结果更偏向于TC型震源.并且，对比正向Q值误差与负向Q值误差对应的反演结果可以发现，地层品 质因子的负向误差对微地震震源机制反演的影响更大.

但是，不同类型的微震源对于地层品质因子误差的抗干扰能力并不相同，而利用随机微震源，我们不能确定地层品质因子的误差对于各种类型微震源的震源机制反演结果的影响程度.因此，本文设计了第二个数值算例，用于分析特定类型的张性震源和剪切震源在地层吸收衰减效应影响下的反演结果误差.

5.2 含随机扰动的四类基本微震源反演

在实际的水力压裂地层中，微震源并不是单一类型的震源机制组成.根据公式(3)，任意的震源机制都可以利用四类基本震源机制(ISO、DC、CLVD^-和CLVD⁺型震源)进行描述.为了模拟实际地层中复杂震源机制，本小节对四类基本震源引入一定组分含量的随机扰动震源项，即：

(20)

其中，M _ST表示基础微震源，M _Per表示随机扰动微震源项，φ表示扰动震源项所占的组分系数(本文计算中，φ值为30%).然后，利用图 2的地质模型和观测系统正演出每个微震源对应的直达P波与S波地震记录，并作为反演过程的输入数据.

通过数值计算，分析地层Q值存在误差时不同数据类型微震源反演结果的影响(如图 7).其中，图 7a、7b、7c和7d分别表示利用P波和S波数据联合反演含随机扰动ISO、DC、CLVD^-和CLVD⁺类型震源的结果，图 7e、7f、7g和7h则表示单独利用P波数据进行反演的结果.图中不同颜色代表不同地层Q值误差情况下的反演结果，当地层Q值没有误差时，真实的反演结果位于图中的“*”号位置.根据图中反演结果的变化情况，可以得出以下结论：

(1) 利用直达P波和S波联合反演的结果更容易受到地层Q值误差的影响，其反演结果主要表现在张性组分的误差变化，地层Q值减小使得所有震源的反演结果向爆炸型震源偏移，地层Q值增大则产生相反的偏移效果；

(2) 当单独利用P波数据进行反演时，其反演结果整体好于直达P波和S波联合反演的结果，但地层Q值减小对震源机制反演的误差影响更大；

(3) 不同类型的震源对于地层Q值误差的敏感度不同，地层的Q值减小对于反演结果的影响要大于地层Q值增大带来的干扰.

为了定量地分析反演结果与反演误差，我们计算了四种类型随机震源机制反演结果的剪切组分(T值)与张性组分(k值)误差(如图 8).通过对比可以看出，利用不同数据进行震源机制反演所得到的T值与k值反演误差具有明显不同的变化规律.

当地层Q值的正向误差逐渐增大，即地层Q值逐渐增大时，利用直达P波和S波联合反演的震源机制T值结果受到地层Q值误差的影响很小，而其k值反演结果对于地层Q值误差的增大很敏感；然而，利用直达P波单独反演时，不同震源类型的反演误差变化很大；其中，ISO和DC型随机震源的k值反演误差较小，T值反演误差较大，而CLVD^-和CLVD⁺型随机震源的T值误差和k值误差均随地层Q值误差的增大而增大，两者对比可知，CLVD^-型随机震源的T-k值误差的增大速率更低.

另一方面，当地层Q值的负向误差逐渐增大，即地层Q值逐渐减小时，ISO型随机震源的T-k值的反演误差出现了不同的变化规律，即单独利用直达P波反演得到的误差明显小于直达P波和S波数据联合反演的误差.DC、CLVD^-和CLVD⁺型 随机震源的T-k反演误差随Q值误差增大的速率明显大于单独利用P波反演的结果；并且，对于T值，利用直达P波和S波反演效果较好，而对于k值，单独利用直达P波反演效果较好.

综合上述所有类型震源机制的反演结果和其对应的反演误差，结合各类震源机制对于地层Q值误差的敏感度，本文认为在震源机制反演的过程中，可以利用直达P波和S波数据联合反演震源的T值，单独利用直达P波反演震源的k值，从而提高反演算法的稳定性和抗干扰性.

5.3 震源机制反演分辨率分析

为了深入地探索地层Q值误差对震源机制反演的影响，本小节同样采用费雷谢偏导矩阵SVD分解方法对不同Q值条件下的震源机制反演分辨率进行了对比分析(如图 9).其中，图 9a、9c为利用直达P波和S波联合反演时，地层分别具有-60%和150%误差时的反演分辨率；图 9b、9d为利用直达P波单独反演时，地层分别具有-60%和150%误差时的反演分辨率.将图 9与图 4对比分析可知，在反演过程中不同观测数据对地层Q值误差的敏感度是不同的.

当利用直达P波和S波数据联合反演时，地层Q值负向误差逐渐增大(如图 9a)，即地层衰减效果增强使得奇异值Λ_i的分布稀疏性降低，这会影响最终解的稳定性.另外，奇异值Λ_i所对应的特征向量 W _i 也发生了改变，这使得观测数据对M₁₁、M₁₂和M₁₃的反演约束性降低，从而导致反演结果的误差增大.另一方面，随着地层Q值正向误差逐渐增大(如图 9c)，即地层衰减效果减弱时，费雷谢偏导矩阵对应的奇异值Λ_i稀疏性变化不大，而特征向量 W ₅和W ₆的独立性明显下降，这使得从观测数据中无法求解出准确的M₂₃和M₃₃.

当利用直达P波单独反演时，无论地层Q值的误差正向增大还是负向增大，费雷谢偏导矩阵对应的特征向量矩阵 W _i 没有发生明显变化，这说明观测数据对于震源机制的约束性没有改变.但是，当地层Q值的误差负向增大时，奇异值Λ_i的稀疏性降低，最终导致利用直达P波反演时的误差小于直达P波和S波联合反演的误差，且对于地层Q值的正向误差不敏感.由此可知，当观测系统的方位角足够时，应只利用P波数据进行反演工作，以此提高反演结果的抗干扰性.

6 结论

微地震震源机制的准确反演对于水力压裂的成功实施具有重要的指导意义，通过Hudson图示法可以统计整个压裂过程中微地震震源机制的主要类型与性质，从而有效地评估整个压裂过程.由于微地震信号具有主频高，能量小的特征，即使采用井下观测系统，检波器接收到的微地震信号仍然会受到地层吸收衰减效应的严重影响.因此，本文重点分析了在均匀各向同性的吸收衰减地层条件下，微地震震源机制反演问题.考虑到实际地层Q值计算结果不稳定和误差较大的特性，本文利用数值算例分析不同地层Q值误差条件下不同观测数据对于微地震震源机制反演的分辨率变化.同时，数值计算的结果表明不同类型的微地震数据在反演过程中对于地层Q值的误差具有不同的敏感度.在观测系统的方位角足够的情况下，震源机制反演应优先考虑单独使用P波数据，或者采用直达P波和S波数据联合反演震源的T值，单独利用直达P波反演震源的k值，以此提高反演结果的稳定性和准确性.

通过理论分析和数值计算，本文也为提高微地震监测技术提供了理论基础，所获得的结论也可用于实际微地震监测的评价.然而，在实际的水力压裂过程中，目标层的物性参数产生了巨大的改变，这不仅会改变微地震震源机制的特征和微地震波场的传播规律，也会对整个反演过程产生严重的影响.这也是微地震监测技术当今所要面对的重大挑战，更深入的理论研究与详尽的岩石物理实验测试仍需开展.