V.Vavryuk

引言

目前人们提出了几种根据地震震源机制确定构造应力的方法（参见 Maury et al，2013）。最普遍应用的方法为 Michael（1984）、Gephart与Forsyth（1984）、Angelier（2002）提出并由Lund与Slunga（1999）、Hardebeck与 Michael（2006）、Arnold 与Townend（2007）、Maury等（2013）及他人修改和发展的方法。这些方法通常假定：（1）构造应力在区域内均匀；（2）地震发生在方位变化的先期存在断层上；（3）滑动向量表明断层上剪应力的方向（所谓的 Wallace－Bott假 设；参 见 Wallace 1951；Bott，1959）。如果上述假定得到满足，应力反演方法就能确定应力张量的4个参数：确定主应力方向的3个角度σ1、σ2、σ3和形状比R（Gephart and Forsyth，1984）：

这些方法不能反演应力张量的另外两个参数。因此应力张量的迹通常假定为零（Michael，1984）：

并且应力张量是归一化的。

然而，由震源机制反演应力要面对一个共同的基本麻烦。为了应用Wallace－Bott假设，应力反演需要确定哪个节面为断层。这是困难的，通常需要一些额外的信息（如，地质证据或震源丛集的指示），因为震源机制固有的模糊性不能区分断层和辅助节面。不幸的是，如果断层和辅助面在应力反演中交换，结果会有偏差且不正确。由于这个原因，Lund和Slunga（1999）提出应用所谓的断层不稳定性约束，并用给定应力场中更不稳定的那个节面确定断层。采用该准则，他们修改了Gephart和Forsyth（1984）的方法并提高了其效率。本文按照同样的方式修改Michael（1984）的方法。这样就导致应力不再像Michael初始的方法那样一步计算出，而是以迭代方式计算出。迭代过程需要区分断层面，并需要确定更准确的应力场。我们用数值试验表明，迭代联合反演应力和断层取向的方法较Michael方法更稳健，并且能够得到更准确的形状比值。最后，采用克里特岛中部和捷克共和国波希米亚西部震群的实际数据验证了所提反演方法的运行效率。

1 MICHAEL方法

Michael（1984）提出的应力反演方法采用断层上的正应力σn和剪应力τ表达式：

式中δik为克罗内克符号，T为沿断层的牵引力，n为断层法向，N为断层剪应力矢量的单位向量。随后，（4）式修改为：

为了能估计（4）式右边，Michael（1984）采用了Wallace－Bott假设并把沿断层的剪应力方向N看作是剪切运动的滑动方向s。他进一步假定活动断层上的剪应力τ对所有地震具有相同的值。因为该方法不能确定绝对应力值，τ在（5）式中归一化为1。其后，（5）式写成矩阵形式为：

式中t为应力分量的向量。

A为根据断层法向n计算的3×5矩阵，

s为滑动向量的单位方向。将（8）式扩展为已知断层法向n和滑动方向s的K 个地震震源机制，我们得到应力张量的5个未知分量的3 K个线性方程组。最后，联合（2）式，采用L2范数的广义线性反演（Lay and Wallace，1995的6.4节）得到

根据以上方程，Michael方法的基本缺陷为需要知道断层的取向（Michael，1987）。如果Michael方法采用不正确的断层面取向，得到的应力张量的准确性就会降低。Michael（1987）进行了一系列的数值试验并发现，形状比很容易被扭曲。另一方面，该方法运行速度很快，并能重复运行。因此解的置信范围可以根据标准的重取样方法估计（Michael，1987）。如果震源机制的断层面取向未知，每个节面在重取样过程中有50%被选为断层的概率。

2 断层不稳定性约束

震源机制中断层面的模糊性在其他的应力反演方法中也引起了困难。例如，Gephart和Forsyth（1984）提出的计算断层面上剪应力牵引方向与断层滑动方向的方法是对两个节面均做为断层面进行试算，选择产生较小剪应力牵引方向和断层滑动方向偏离较小的节面为断层面。这种处理策略对于无噪声数据容易求解，但不正式，并且对于有噪声的实际数据效率较低（参见Lund and Slunga，1999）。Michael（1987）评述了一些其他‘断层选择算法’并在断层滑动数据集进行了检测。他也提到采用断层破裂准则的可能性。Lund和Slunga（1999）进一步详尽说明这种方法，他们分析应力场中两个节面的取向，发现两个节面之一稳定性更差，并且更容易受到剪切断层活动的影响。这种对破裂的敏感性采用莫尔―库仑破裂准则来量化（参见图1）。根据该准则（Beeler et al，2000；Scholz，2002；Zoback，2010），激活断层上的剪切牵引力τ必须超过临界值τc，τc根据内聚力C、断层摩擦系数μ、正应力σn和孔隙压力p计算，即：

式中

断层面就是库仑破裂应力Δτ具有较大值的节面。在估计（10）式的过程中，我们需要（11）式中的摩擦系数μ。内聚力C和孔隙压力p并不需要，因为我们比较的是应力的相对差别。裂纹的摩擦系数μ在实验室的岩石样品中测定，范围多数在0.6～0.8之间（Byerlee，1978）。地壳断层的摩擦系数值是相似的（Vavrycˇuk，2011），但对于一些大尺度断层，如圣安德烈斯断层，也有0.2～0.4的较低值的报道（Scholz，2002）。

图1 莫尔－库仑破裂准则。τ和σ分别为剪切和有效正应力；σ1、σ2和σ3为有效主应力。红区（原图为彩色图——译注）表示满足莫尔―库仑破裂准则的可能断层面取向。蓝点表示应力张量最优取向主断层，C表示内聚力。图的上、下半平面对应于共轭断层

图2 断层不稳定性在莫尔圆中的定义。红色实心圆（原图为彩色图——译注）标出主断层上失稳特征I＝1的牵引力。黑色实心圆标出任意取向断层上不稳定性为I的牵引力

牵引力τ和τc在（10）和（11）式中总假定为正，因此原则上，仅绘出图1中莫尔圆的上半平面就足够了（参见Jaeger et al，2007，图4.9；Zoback，2010，图5.9）。按下列方式考虑剪切牵引力τ的符号可以获得整个莫尔圆：

式中N为沿断层的剪切牵引力τ的方向，e（3）为最小主应力的方向。绘出整个莫尔圆对于区分共轭断层更方便（参见图1的左图）。

区分断层面的另一种方法是估计Vavrycˇuk等（2013，式3）提出的断层不稳定性I：

式中τc及σc为最优取向断层面的剪切牵引力和有效正应力（图2的红点），τ和σ为所分析断层面上的剪切牵引力和有效正应力（图2的黑色实心圆）。因为（13）式与绝对应力值不相关，断层不稳定性I仅由摩擦系数μ、形状比R及定义为断层面与主应力偏离的方向余弦n来估计。如果我们按下式规定应力张量：

式中正值表示压缩，我们得到

并最终得到

式中

图3 应力反演数值试验所用的数据例子。该图给出了满足莫尔—库仑破裂准则的200个无噪声震源机制数据。左右两边的图为全部和减少震源机制多样性的数据组。（a，b）莫尔圆图；（c，d）P／T轴；（e，f）对应节面。在（c）和（d）图中，P轴用红色圆圈（原图为彩色图——译注）标示，T轴用蓝色十字标示。在（c）中P／T轴组成了所谓的蝴蝶两翼，（d）中为一翼，参见Vavrycˇuk（2011）。应力轴σ1、σ2 和σ3 的方位角和倾角分别为：115°／65°，228°／10°和322°／23°。形状比R为0.7，内聚力C为0.85，孔隙压力p为0，摩擦系数μ为0.6

因子R为（1）式所定义的形状比，n为主应力方向坐标系中表示的断层法向。断层不稳定性I范围为0（最稳定断层）到1（最不稳定断层）。最不稳定断层为剪切断层活动的最优取向断层（图2中的红色实心圆），称作主断层（Vavrycˇuk，2011）。注意（16）式表达的断层不稳定性与Vavrycˇuk等（2013，式6）似乎不同，原因是那里的正应力表示拉张，而不是本文中定义的压缩。

3 迭代应力反演

Lund和Slunga（1999）将断层不稳定性约束用于Gephart和Forsyth（1984）的应力反演，提高了该方法的效率。这里我们将该约束用于 Michael（1984）的方法。因为Gephart和Forsyth的算法为非线性算法，拟合差采用网格求得，断层不稳定性约束可以方便地用于拟合每次应力的估计中。相对于Gephart和Forsyth的方法，Michael的方法为线性的，引入断层不稳定性约束会使应力反演只能迭代求解。首先，采用标准的Michael方法，不考虑约束及断层面取向的任何信息。找到主轴方向和形状比后，用这些值估计所有震源机制节面（16）式的不稳定性。断层面就是稳定性较差的节面。在第一次迭代过程中的断层面取向用在第二次迭代的Michael方法中。这种步骤持续进行直至应力收敛于最优值。

图4 主应力方向（左图）和形状比（右图）作为震源机制数目和噪声水平函数的平均误差。采用减少震源机制多样性的数据（单翼数据）进行反演。上图：Michael方法；下图：迭代反演方法。应力方向误差为真实应力轴和反演应力轴的平均值。误差在50次随机重复反演中进一步平均。应力方向的误差单位为度，形状比的误差单位为百分比

当用（16）式估计断层不稳定性时，需要摩擦系数μ的值。正如前面所提到的，断层摩擦系数范围在0.2～0.8之间，但通常是未知的。然而数值试验表明，反演对摩擦系数μ相当不敏感，因此在反演中给定摩擦系数为平均值，比如0.6即可。另外一种方法是选用不同的摩擦系数进行反演，选取能够产生反演数据的断层的总体具有最高不稳定性的摩擦系数。这种方法用于后面的合成数据测试及实际数据的应用中。

4 精确度测试

迭代应力反演的稳健性和准确性可以采用数值模拟来检验。我们给出自20到200的不同震源机制数的数值试验表现。选择的震源机制满足莫尔—库仑准则（参见图3）。而后，矩张量采用其范数的0～100%均匀噪声进行污染。这里的范数采用矩张量本征值的绝对最大值（所谓的 “谱范数”）来计算。矩张量的噪声分解到被污染的震源机制的走向、倾角和滑动角之中，采用这些震源机制反演应力场。实际断层法向和含有噪声的断层法向的偏离值范围为0°到35°。因为用于反演的震源机制根据矩张量计算，因此没有从两个节面中区分断层面的信息。采用50次加入随机噪声，并采用震源机制两种类型的子集进行重复反演。第一个子集包含投影到莫尔圆图的上半平面和下半平面（图3a）。根据 Vavrycˇuk（2011），P／T 轴组成了所谓的 “蝴蝶两翼”模式（图3c）。第二个数据子集包含了相同的震源机制数，但它们的多样性大大减少。选择的震源机制仅投影到莫尔圆的上半平面（图3b），并在P／T轴绘图中仅占了一翼（图3d）。各种不同噪声的两组数据反演的主应力方向和形状比与真值进行比较，并估计误差（参见图4，5）。与此相似，反演过程从节面中识别的断层也与真值进行比较，从而估计成功识别率（参见图6）。

为了避免不稳定性约束中摩擦系数对反演的影响，按照间隔为0.05自0.2到1.2的摩擦系数，重复多次反演。每次反演中识别出的断层的总体不稳定性在反演中可以给出。能够产生最高总体断层不稳定性的摩擦系数看作最优值。

图5 与图4标注相同，但对应的是完全震源机制多样性的数据组（两翼数据）

图6 采用迭代反演得到的断层成功识别率（左图）及摩擦系数误差（右图）随震源机制数目和噪声水平的变化。上图：减少震源机制多样性的结果。下图：完全震源机制多样性结果。这些值是由50次随机噪声实现计算的。识别成功率为正确识别断层的相对数目

数值试验表明，应力反演的准确性根据反演中所用的震源机制数目变化和数据的噪声水平而变化。Michael方法和迭代反演方法对于断层法向和滑动方向误差为20°的20个污染的震源机制数据均能得到误差小于12°的主轴方向（参见图4和图5的左图）。对于更准确的震源机制或更大的数据集，精确度会更高。两种方法的差别几乎看不到，迭代方法对于较低噪声水平的数据的结果精度略高。其理由比较简单。如果断层没有识别出，Michael反演永远不能得到正确的应力方向，即使是采用大量无噪声的震源机制数据也是如此。相对而言，迭代应力反演克服了这种缺陷，因为断层取向在迭代过程中确定。然而，Michael方法和迭代联合反演方法的明显差别表现在形状比的确定上（参见图4和图5的右图）。即使采用大量的无噪声震源机制数据，得到的形状比的平均误差差不多为Michael反演误差的20%。采用有噪声的震源机制数据误差会更大。除了反演震源机制的数目较少和噪声较大的情况，迭代方法得到的平均误差小10%。

另外，迭代反演方法能够区分震源机制中的断层面及断层的总体摩擦系数。图6（左图）表明，采用单翼平均误差为15°的噪声震源机制数据，断层识别正确率大于85%。对于双翼震源机制，成功率达90%。估计的摩擦系数的值对于相同类型的数据精度为10%～15%（图6，右图）。

5 在实际数据中的应用

5.1 里克特岛中部

迭代联合反演的效率采用Angelier（1979）和 Michael（1984，1987）已经分析的里克特岛中部的数据来证实。该数据为现场测量的38个断层取向和滑动方向，所以震源机制中的断层面和辅助面是已知的。在应力反演过程中，我们通过随机选择两个节面的断层修改输入数据，并检验迭代过程能否识别出正确的断层面。

图7给出了输入震源机制的P／T轴、节线及应力反演结果。得到的最大压缩轴几乎垂直（图7a），莫尔圆表明（图7c）输入数据覆盖了两个蝴蝶翼。从38个震源机制中，成功识别了36个断层面（图7d）。形状比为0.88，最优摩擦系数为0.75。这里形状比的定义与 Michael（1984）定义的φ＝（σ2－σ3）／（σ1－σ3）不同，按照 Michael的定义为0.12。

为了分析主轴方向和形状比的不确定性，用有噪声的震源机制重复进行应力反演。我们重复反演了噪声数据1 000次。断层取向和滑动方向的平均误差为6°。因为我们不知道输入数据的精度，因此该值看上去有些任意。因为我们的目的是比较两种反演方法的效率，更为精确地确定误差并不是关键。

表1和图8表明，不管断层面从节面中随机选择（图8a）还是在反演过程中正确区分（图8b），主应力方向均具有相似的精度。然而形状比的行为明显不同。随机选择断层面的数据得到的形状比分布相当宽广，其最大值随着反演中选择断层的正确性而明显不同。最大误差对于随机选择断层和正确选择断层分别为44%和12%。

图7 用克里特岛中部数据迭代反演的应力（Angelier，1979；Michael，1984，1987）。（a）得到的主轴方向及P／T轴。（b）节线。（c）断层在莫尔圆上的位置（蓝加号）（原图为彩色图——译注）。（d）应力反演中识别的断层。（a）中的P轴和T轴分别采用红空心圆和蓝加号表示

表1 应力反演结果

图8 用Michael方法（a）和迭代反演方法（b）得到的主应力方向的置信范围和形状比的直方图（c）。（a）和（b）中的红、绿、蓝色（原图为彩色图——译注）分别对应于σ1、σ2和σ3应力方向。（c）中的灰色和红色直方图分别给出了Michael方法和迭代方法得到的形状比直方图

5.2 波希米亚西部震群区

迭代联合应力反演进一步用捷克共和国波希米亚西部2008年的震群数据举例说明（Fischer et al，2010；Vavrycˇuk，2011；Davi and Vavrycˇuk，2012；Vavrycˇuk et al，2013；Fischer et al，2014）。由于微震被很多地方地震台记录到，在震源球上有很好的台站覆盖，高质量的数据较为适合精确提取震源机制和构造应力场。

我们选择了具有较高信噪比的所有WEBNET台站记录的167个微震，对直达P波的垂直分量进行矩张量反演。格林函数采用射线方法计算（Cˇervený，2001），完全矩张量在时间域采用广义线性反演确定（Vavrycˇuk and Kühn，2012）。震源机制的精度采用输入数据加噪声和位置偏离（震中误差为250m，深度误差为500m）的重复反演进行评估。断层法向和滑动方向的平均误差估计为6°。

震源区的构造应力根据167个震源机制分别用Michael（1984）方法和迭代联合应力反演方法计算（图9；参见表1）。误差的计算与克里特岛中部数据相似，采用反演添加噪声的震源机制数据给出。我们加入噪声计算了1 000次，得到了与震源机制相似的精度（断层法向和滑动方向的平均误差为6°）。在迭代应力反演中，最终的应力通过3～4次迭代得到。断层面覆盖了蝴蝶两翼（图9c）。应力以水平σ3轴为特征，σ1和σ2轴明显偏离了水平和垂直方向（图9a）。主应力轴的方向及其不确定性对于Michael方法和迭代应力反演方法几乎一致（图10a和b）。这两种方法的主要差别在于形状比。与合成数据检验结果一样，Michael方法得到了具有较大误差的形状比偏差值（图10c）。迭代反演既能得到精确的主应力方向，也能得到精确的形状比值（表1）。此外，迭代反演还能够区分断层面。正如期望的那样，断层面集聚在莫尔—库仑破裂准则的有效区内（图9c）。估计得到的断层取向被震源丛集分析所确认（Vavrycˇuk et al，2013）。震源沿两个断层系丛集（图9d）被认为 与 该 区 域 的 主 断 层 有 关 （Vavrycˇuk，2011，图10）。反演中发现最优的摩擦系数为0.85。

图9 与图8相同，但对应于波希米亚西部的数据。（d）中的蓝色线（原图为彩色图——译注）给出了根据地震震源丛集分析得到的波希米亚西部地区的主断层（参见Vavrycˇuk et al，2013）

6 讨论

Michael方法获取主应力方向快速而准确。即使随机选择震源机制节面作为断层面也可以得到相当合理的结果。然而，当断层面没有正确选择时，形状比的准确性大大下降。即使采用大量精确的震源机制数据，估计的形状比值也是非常粗略的。因此，比起主应力方向，形状比的值对于是否正确选择断层面更加敏感，将辅助节面代替断层面将导致较大误差。

形状比的精确值可以用联合反演应力和断层取向的迭代算法求得。在反演中，采用断层不稳定性约束，将那些更不稳定、更易于断裂的节面作为断层面。在反演中加入断层不稳定性约束会导致必须采用迭代算法，不是一步就计算出应力，而是迭代计算出应力。在初始迭代时，随机选择断层面，采用Michael方法求解。在第二次或高次迭代中，通过非稳定性约束区分断层面，然后采用Michael方法进行求解。数值实验表明，经过3～4次迭代就可以收敛。因为主应力方向即使第一次迭代也很容易获得，后续迭代主要是区分出断层面，并提高形状比精度。这就是收敛如此之快的原因。

在分析实际数据集时，获得的应力方向、形状比、断层面选择及总体摩擦系数的不确定性不像 Michael（1987）那样用重取样技术计算。这里计算的不确定性是无噪声数据和有噪声数据反演1 000次所得结果之间的最大差别。该方法由于下述原因而具有优势：首先，该方法考虑了由于对断层取向和滑动方向施加了不同的噪声水平，一些节面相对于其他节面更加不确定。其次，得到的不确定性更加实际。例如，形状比的直方图可以非常宽广（图8c和10c），并且真值不必靠近直方图的最大值。

迭代反演方法的效率用克里特岛中部和捷克共和国波希米亚西部的数据进行了检验。对于这两组数据，该方法效果很好，得到了明显比从震源机制中随机选择断层面进行反演更准确的形状比结果。在绘制莫尔圆图时，反演中识别的断层表现了较高的断层不稳定性，并且集中在莫尔—库仑破裂准则有效区内。这就清楚地表示数据与断层不稳定性模型一致，应用迭代应力反演是合适的。如果反演中识别的断层没有表现出较高的不稳定性，断层不稳定性模型可能不完全适合，迭代反演的结果将不可靠。这可能是研究区的震源机制具有较大误差，或者应力场更加复杂，或者当地构造应力环境所致。

Michael方法非常流行，已用在数个应力反演软件中，如SATSI（Hardebeck，2006；Hardebeck and Michael，2006）或其Matlab修改版本 MSATSI（Martínez－Garzón et al，2013，2014；Ickrath et al，2014）。由于迭代应力反演基于的是Michael方法，因此可以很容易地用于这些软件中，提高结果的精度。本研究反演的Matlab软件叫做STRESSINVERSE，在 网 页 上 提 供 （http：／／www.ig.cas.cz／stress－inverse，最后访问时间2014年6月27日）。