引言

地壳中的应力是地震和其它构造运动过程的直接动力来源． 区域地壳应力场及其时空变化的研究为深入了解地震机理、 构造负载及地震应力的相互作用提供了可靠的信息．

岩石圈内真实构造应力的计算和测量都是非常困难的，但可以利用其它的一些地球物理观测资料反演出构造应力的某些特征，如应力方向、 应力相对大小等. 这一方面国内外早已开展了大量的研究，并取得了丰硕的成果(Zoback，1992； 谢富仁等，2004)． 20世纪80—90年代，许忠淮等利用小震初动符号采用区域综合震源机制解的方法，对华北和中国大陆其它地区的构造应力方向进行了大量的研究，取得了很多具有开创意义的成果(许忠淮， 1983，1989； 许忠淮，1985)； 利用野外观测的第四纪断层活动擦痕数据反演构造应力张量，也可以得到区域构造应力场的基本特征(Nemcok，Lishle，1995； 谢富仁， 1993，2001； 荆振杰等，2008； 张红艳等，2009). 此外，还可以使用跨断层形变观测资料(谢富仁等，2008)以及钻井的井孔崩落特征或水压致裂的情况对构造应力场进行研究(许忠淮，吴少武，1997； 许忠淮等，1999). 这些研究所得到的构造应力场的方向已经为近年来GPS观测的结果所证实(袁金荣等，1999； 刘峡等，2006)． 值得一提的是，1986—1990年，在国际岩石圈计划下开展了“世界应力图”项目的合作研究，综合各种资料完成了全球范围内应力图的绘制，取得了具有里程碑意义的成果(Zoback，1992； 许忠淮，1994)．

使用大量的地震震源机制解资料可以推断区域应力场的特征，对此地震学家已经发展了许多经典的方法. 例如，Michael(1984)，1987)的LSIB(linear stress inversion with bootstrapping)方法，Gephart和Forsyth(1984)的FMSI(focal mechanisms stress inversion)方法，Angelier(2002)基于滑动剪切应力分量的反演方法，以及崔效锋和谢富仁(1999)逐次收敛法等． 这些方法在实际研究中都有着广泛的应用(Petit et al，1996; Plenefisch，Bonjer，1997; Lund，Slunga，1999; Hardebeck，Hauksson，2001;Bressan et al，2003;Kastrup et al，2004;Boness，Zoback，2006;Carena，Moder，2009)．

中小地震的震源机制参数通常有很大的不确定性，但是大量的离散分布的数据可以准确地约束应力张量的方向(Hardebeck，Michael，2006)． 近年来，使用震源机制解推断构造应力特征的方法又有了新的发展(Loohuis，van Eck，1996;Abers，Gephart，2001)，一方面基于大量的资料进行小区域精细构造应力场的研究(Hardebeck，2006；Fay et al，2008)，另一方面通过增加阻尼约束、 应用概率信息准则等数学手段来更好地约束反演结果(Hardebeck，Michael，2006；Arnold，Townend，2007；Terakawa，Matsu ’ura，2008).

自20世纪80年代的“七五”、 “八五”攻关中集中计算小震震源机制以来，利用初动符号或振幅比资料求解区域中小地震的双力偶模型参数的工作逐渐成为国内分析预报部门的日常工作，积累了大量的小震震源机制解的资料． 尤其是近年来，随着数字观测台网的加密，区域小震震源机制的结果日见丰富，其质量也更加可靠(胡新亮等，2004). 也由此对不同区域开展了构造应力场的研究(吴建平等，2004； 康英等，2008； 李瑞莎等，2008； 刘泽民等，2011； 李莹甄等，2011)． 我们利用山东省地震预报研究中心所掌握的资料，使用Hardebeck和Michael(2006)的区域应力张量阻尼反演方法对山东地区的背景应力场进行研究，并分析探讨区域应力场与区域内的深大断裂 郯庐断裂带的关系．

1.   数据与资料

1.1   震源机制解数据

前人对山东地区中小地震的单震震源机制已经开展了大量的工作(魏光兴，李秉锋，1980； 魏光兴等，1982； 周翠英，魏光兴，1987； 魏光兴等，1993； 董旭光等，1999； 周翠英， 2003，2005； 山长仑等，2007； 郑建常等，2010； 郑建常，陈运泰，2012)． 自2006年以来，随着数字化测震台网的投入使用，山东省地震预报研究中心使用中国地震局数字地震实验室推广的Snoke的FOCMEC程序(Snoke et al，1984; 刘杰等，2004)以及MapSIS软件中的初动符号格点尝试法程序进行例行的震源机制求解，基本上每次省内M_L≥2.0地震后都会给出相应的震源机制解． 通过搜集汇总这些资料，我们总计得到了山东地区1970—2012年5月共387次中小地震的震源机制解(图 1a)． 采用Frohlich(2001)的震源机制量化分类和评估方法(图 1b)，可以看出这些小震的震源机制类型主要是走滑和正断的断层作用，合计约占总数的76%．

图  1  本文使用的小震资料空间分布(a)及震源机制解分类(b)图

Figure  1.  Locations of 387M_L≥2.0 earthquakes used in this paper (a) and classification of focal mechanisms on the octant equal-area projection according to Frohlich (2001) (b)

下载: 全尺寸图片 幻灯片

1.2   山东地区整体应力场特征

图 2给出了本文所使用的震源机制解的P轴和T轴在震源球上的投影分布. 可以看出其P轴方位的优势分布为ENE-WSW向，T轴方位的优势分布为NNW-SSE向． 将P，T轴方位归一化到[0，180°]之间，统计研究区内所有震源机制解的P，T轴方位的优势方向(图 3)，结果得到山东地区的平均P，T轴取向的均值及标准差分别为P=86.9°±38.6°，T=173.2°±38.5°．

图  2  本文使用的小震震源机制解的P轴(a)、T轴(b)下半球极射投影图

Figure  2.  Lower hemisphere projection ofPaxes (a) andTaxes (b) used in this study

下载: 全尺寸图片 幻灯片

图  3  本文使用的小震震源机制解的归一化P轴方位(a)与T轴方位(b)统计直方图

Figure  3.  Histogram of normalized azimuths ofPaxes (a) andTaxes (b) of the focal mechanism solutions used in this study

下载: 全尺寸图片 幻灯片

使用ZMAP软件包中的FMSI程序反演山东地区整体的应力场，得到3个应力主轴的方向(图 4). 结果表明，山东地区的最大主应力方位为82°，最小主应力方位为172°，与前面统计结果相吻合，应力相对大小的平均值R=(σ₁-σ₂)/(σ₁-σ₃)=0.67．图 4给出了3个应力主轴方向的95%置信区间． 六省市震源机制小组(1981)对山东及周边应力场的研究表明，苏鲁皖豫地区处于北东东向水平压应力和北北西向水平张应力的现代构造应力场作用之中，其地震应力场的优势分布为主压应力轴约76°±28°，主张应力轴约344°±27°. Wan(2010)基于中国地壳应力数据库和美国哈佛大学矩心矩张量目录给出了中国大陆现代构造应力场. 其中山东地区的最大主应力轴方位为70°—90°，最小主应力轴方位为160°—194°. 本文结果与上述结果基本一致．

图  4  山东地区主应力轴及其置信区间的下半球等面积投影

Figure  4.  Lower-hemisphere equal-area projections of principle stress axesσ₁， σ₂andσ₃， and their confidence regions

下载: 全尺寸图片 幻灯片

2.   理论与方法

地球内部任一点的构造应力状态可以由6个独立的参数完全定义，表示为6个分量的3×3对称矩阵或者3个特征向量与地理坐标的夹角以及对应的3个特征值． 在理想情况下，钻孔应力测量可以推断6个应力张量的参数，但在通常的地质学和地震学研究中，一般只能给出其中的4个参数，即3个主应力的方向和中等主应力相对大小的参数．

受地下介质结构、 地质构造条件及断层展布情况等因素的影响，一定区域的应力空间分布是不均匀的(如王勤彩等(2009)对汶川地震余震区应力状态分段特征的研究)． 当需要考察应力状态的空间分布特征时，通常的方法是将区域分成很多小区，对每个小区的震源机制拟合应力张量(Petit et al，1996 Hardebeck，Hauksson，2001;Bressan et al，2003)． 分区方式不同和分区大小的变化都可能对反演的应力张量结果产生影响(Maury，Cornet，2011)． 并且由于每个小区都是独立反演，在某些情况下，计算结果显示出来的应力偏转等特征有可能相当的大． 相邻两个区内主应力方向的不同或许反映了应力场真实的变化，但也可能是震源机制数据的误差或反演的约束较差造成的假象． 受实际地震资料的限制，某些小区的震源机制可能对一组应力张量都拟合得很好，因而对这样的一个小区就有多个应力模型可以满足数据，如何增加约束来选择更合适的应力模型就成为必要的问题(Hardebeck，Michael，2006)．

由于地壳中的应力分布存在连续性，因此，为了解决反演的应力模式依赖于分区的问题，并且为了更准确地得到区域范围内的应力变化的特点，Hardebeck和Michael(2006)参考其它地球物理领域中的阻尼最小二乘反演方法，提出了区域尺度的应力张量阻尼反演方法． 他们构建一组可调整的阻尼参数模型，引入平滑约束来抑制相邻单元之间应力模式的差异，然后使用最小二乘法可以得到稳定解．

2.1   阻尼区域尺度应力反演

构造应力张量有6个独立分量，对于一组地震的震源机制，设应力张量不存在各向同性分量时，σ₃₃=-(σ₁₁+σ₂₂)，则应力张量只有5个分量

Michael(1987)的反演方法通过最小化每个断层面上滑动向量和解得的剪切应力向量的差异，求解方程

式中，d 为由震源机制解求得的滑动矢量

式中，s_ij是第i次地震的单位滑动矢量的第j个分量(j=1，2，3)． 矩阵 G 由每个震源机制的断层法线矢量组成．

对于区域应力场，设研究范围分为I×J个二维空间网格． 第i行第j列的格点上的应力张量由向量 m_ij表示，各点上的观测数据由向量 d_ij表示，数据核矩阵记为 G_ij，则新的模型向量

新的数据向量

新的数据核矩阵

则对于区域应力场求解问题，方程(2)可写为

对此方程增加“平滑”约束，最小化相邻各点的应力张量的差异，引入阻尼矩阵

式中，I是5×5的单位矩阵．

阻尼反演最小化两个值的权重和为数据错配、 模型长度． 模型长度表示为向量的二阶范数，包含了每对相邻格点的每个应力张量的差异

数据与模型理论值的错配写为

阻尼最小二乘问题的解，需要同时最小化等式(9)和式(10)，即

式中，e是标量阻尼参数，控制最小化问题中数据错配和模型长度的相对权重． 由于 G，D 都非常稀疏，大括号内的矩阵也是稀疏矩阵，故使用共轭梯度方法求解等式(11) ．

2.2   “真实”最大水平主压应力方向的确定

记最大水平主应力为σ_Hmax. 识别最大水平应力的常用方法是简单地采用最大近水平主应力的水平投影方向，但该情况只有当应力张量中的其中一个主应力轴严格垂直时是适用的． 应力的Anderson模型可以很好地解释不同构造条件下观测到的应力状态. 该情况下，σ_Hmax与最大主应力轴平行，除非该轴相对于地表水平面是垂直的(此时，σ_Hmax平行于介质应力轴)． 但是，一般情况下由地震数据估计得到的应力张量的结果中，几乎没有应力轴是严格垂直于地表水平面的． 这时，σ_Hmax不能简单地采用最大近水平应力的水平投影，二者之间在方向上可能存在最大数十度的偏差(Lund，Townend，2007)． Lund和Townend(2007)从数学上推导给出了由应力张量严格计算σ_Hmax方向的公式和方法，可以得到“真实的”最大水平主应力方向．

已知构造应力张量 S，由此可以得到3个主应力σ1≥σ2≥σ3及其对应的主应力矢量 s1，s2，s3，在Aki和Richards(1980)定义的地理坐标系下(北、 东、 下为正)

作用在我们感兴趣的垂直平面上的法向正应力S_n与该法向应力的方位角α之间存在如下关系(Lund，Townend， 2007，式5)：

式中，α为法向正应力S_n的方位角． 一般情况下，S_n在区间[0，π)具有一个最大值和一个最小值，使得上式为零，分别对应了水平最大主应力和最小主应力．

3.   计算结果与分析

3.1   阻尼系数的确定

在对公式(11)进行实际反演计算过程中，阻尼参数e的取值对反演结果起着重要的作用．由于该参数控制着理论值与观测数据之间错配值和模型长度(即模型复杂程度)的相对权重，简化模型，则错配值升高，反演误差增大； 反之提高错配的相对权重，则反演误差减小，而模型逐渐变得复杂化，甚至失去了阻尼约束的意义． 因此通常是根据如图 5所示的折中曲线来选择阻尼参数． 折中曲线的拐点在e≈0.5附近，意味着低于该相对权重值，提高模型复杂程度对反演误差的改善几乎没有作用; 而提高相对权重时，随着模型的简化，反演误差会急剧增加． 因此，我们在反演过程中选择相对权重系数e≈0.5.

图  5  模型长度与数据拟合误差之间的折中曲线图虚曲线旁所标数字是阻尼参数e的取值(见式(11))

Figure  5.  Trade-off between model length and misfit Numbers along the dashed line are values of damping parameterse

下载: 全尺寸图片 幻灯片

3.2   山东地区背景应力场

对山东地区(33.5°—38.5°N，114°—124.5°E)划分0.5°×0.5°的网格，使用前面引述的Hardebeck和Michael(2006)的区域应力场阻尼反演方法，计算每个网格的最佳拟合构造应力张量; 然后使用Lund和Townend(2007)给出的公式计算最大水平主应力的方向． 为了平滑地震空白的区域，没有数据的格点也包括在反演中，但这些点上的结果没有显示． 最终反演结果见图 6．

图  6  山东地区背景应力场反演结果(0.5°网格， 阻尼系数e=0.5)(a) 最大水平主压应力方向； (b) 最大水平主张应力方向； (c) 最大水平主压应力方向误差

Figure  6.  Inverted stress results for Shandong region with damping parametere=0.5, on a 0.5°×0.5° grid (a) The maximum horizontal principal compressive stress orientation; (b) The maximum horizontal principal extensive stress orientation; (c) Errors in estimated principal stress orientations

下载: 全尺寸图片 幻灯片

图 6a给出的从鲁西到黄海海域地区最大水平主应力的方向变化，直观地反映了山东地区背景应力场的作用方式． 从图 6中可以看出，山东地区的最大主应力方向整体上为近东西方向，与前面使用ZMAP得到的结果一致，也与前人(许忠淮， 1983，1989； 王辉等，2004； 张国民等，2004)对华北地区研究得到的结果吻合； 并且最大主应力方向在空间上存在平稳的过渡变化，存在自西向东逐渐的偏转，由ENE向逐渐转化为近EW向，这一现象与汪素云等(1996)以及王绳祖和张流(2002)的研究结果一致．

对比前人对山东地区不同区域的研究结果(魏光兴等，1982； 周翠英，魏光兴，1987； 魏光兴等，1993； 周翠英等，2005)，本文结果与其也非常一致． 周翠英等(2005)根据区域综合机制解研究认为，以鲁中隆起为主的泰沂蒙山区及其周边地区主压应力轴约为60°—80°，主张应力轴约为330°—350°，本文统计结果显示该区的最大主压应力方向变化范围为67.2°—93.5°； 沂沭带以东的胶东隆起和南黄海北部凹陷地区主压应力轴优势方位约为80°—100°，主张应力轴约为350°—360°，本文的统计结果显示该区的最大主压应力方向变化范围为80.0°—112.8°，平均最大主应力方位为89.1°±7.3°． 沂沭带地区最大主压应力轴方位为72°±21°，主张应力轴方位为342°±20°(周翠英，魏光兴，1987)，我们得到这一地区的最大主应力方向为74.0°±5.5°，也与之非常一致．

图 6a显示，鲁豫交界地区主压应力方向差异变化较大，聊考断裂带中段的主压应力方向为WNW，而在其南侧的鲁西南及河南地区主压应力方向显示为ENE. 机制小组(1981)研究也曾发现，豫北地区地震应力主轴方向离散，并且认为造成应力主轴取向摆动较大、 错动类型多变的原因，可能是该区断裂性质、 运动方式复杂． 我们结合地质构造背景分析认为，除此以外，这一现象可能还与该地区的受力状态有关． 鲁西地区在地质构造属于渤海盆地的凹陷区包含了内黄隆起、 东濮凹陷、 开封凹陷、 鲁西隆起等构造单元，地质条件复杂，构造活动明显； 作为整体，鲁西块体一方面受到秦岭—大别构造带的北东向挤压作用，另一方面由于鄂尔多斯地块的旋转平移运动，对华北平原地块也产生了ESE方向的推挤作用．

南黄海北部海域的应力方向偏转较为明显. 其原因为一方面该区震源机制解的个数相对较少且处于研究区的边缘，反演误差较大(图 6c)； 另一方面可能反映了该区真实的复杂应力场(见Wan(2010)文中图 1)．

3.3   山东地区区域应力场与郯庐断裂带的关系探讨

由图 6可以发现，山东地区的构造应力场以郯庐断裂带为界，其东侧胶东半岛及南黄海北部海域地区最大主压应力方向整体呈近EW或ESE方向，而西侧沂沭断裂带及鲁西地区最大主压应力方向整体呈ENE方向． 郯庐断裂带山东段两侧最大主压应力方向的不同，反映了鲁东-黄海地块与鲁西块体受力状态的差异以及运动方式的不同．

作为超壳断裂，郯庐断裂带的断层运动性质直接反映了区域的构造应力状态． 研究发现，现今郯庐带的运动以右旋走滑为主，具有压扭性特征(魏光兴等，1993)． 本文得到的郯庐断裂带山东段两侧最大主压应力方向的不同，恰好与上述研究揭示的郯庐带的运动性质非常吻合． 应力场方向主要受两个因素的影响： 一个是主要的地质单元的构造性质； 另一是与断层的复杂性(Hardebeck，Hauksson，2001)，应力场与断层几何形态之间彼此影响、 相互作用(Townend，Zoback，2004)． 因此我们认为，山东地区郯庐带两侧应力场状态的差异，一方面反映了断裂带两侧地质单元的构造性质的不同，另一方面也与郯庐带的断层形态及其发展演化有关．

4.   讨论与结论

根据1970—2012年5月共387次中小地震的震源机制解资料，使用Hardebeck和Michael(2006)提出的区域应力场反演方法，结合共轭梯度法求解了应力张量的阻尼最小二乘最优化问题，并使用Lund和Townend(2007)推导的利用构造应力张量计算最大水平主应力方向的公式，研究了山东地区的背景构造应力场状态．

计算结果显示，山东地区整体呈现出相对统一的区域应力场. 其最大主应力方向为82°，最小主应力方向为172°； 最大水平主压应力方向主要是ENE—近EW方向，并且在空间上存在平稳的过渡变化，从西往东由ENE方向逐渐偏转为近EW方向. 这一结果与前人(许忠淮， 1983，1989； 汪素云等，1996； 王绳祖，张流，2002； 王辉等，2004； 张国民等，2004)研究结果一致; 一些小的构造单元的结果也与其他研究者(魏光兴等，1982； 周翠英，魏光兴，1987； 魏光兴等，1993； 周翠英等，2005)的结论有较好的吻合． 此外，郯庐断裂带两侧地区最大主压应力方向存在明显不同，反映了不同地质构造块体的不同受力状态．

由于山东地区地震活动水平不高且分布不均，除近几年由初动符号和振幅比得到的机制解资料相对丰富外，早期的机制解个数非常稀疏． 受原始资料的限制，我们对山东地区划分了0.5°×0.5°的网格，部分网格没有小震震源机制资料，个别子区仅有一两个震源机制，在相邻网格的约束下，这些网格也可以得到结果，但结果的可靠性肯定较差． 因此，要了解山东地区更加详细的应力场方面的信息，我们现有的数据质量是远远不够的，这有待未来进一步地研究．

审稿专家对本文的修改提出了中肯的意见和建议； 山东省地震局周翠英研究员提供了大量的原始数据，中国地震局地震预测研究所邵志刚博士提出了富有建设性的建议并给予了相关指导，张浪平博士给予了重要的帮助； 部分图件使用Stefan Wiemer和Max Wyss等开发的ZMAP 6.0生成； 本研究还得益于全国大形势分析预报青年小组成员之间的友好而融洽的讨论. 作者谨在此一并致谢．