京津冀地区地壳应力场特征
2020-09-26王晓山冯向东赵英萍
王晓山 冯向东 赵英萍
摘要：基于标量断层类型值，对京津冀地区及邻区2 187个中小地震震源机制解进行分类，统计结果显示研究区震源机制类型以走滑断层和正断层为主，P轴优势方位为NEE—EW和SWW—EW向;采用MSATSI软件包反演该区1°×1°网格的精细地壳应力场，结果表明：最大主压应力轴最优解的优势方向为NEE—EW向，与P轴优势方位一致;所有网格的相对应力大小R值均小于0.5，表明京津冀地区应力状态偏拉张性质，而且最小主压应力轴的不确定度变化范围相对稳定，表明现今京津冀地区地壳应力场处于一个相对统一的NNW—SSE向的拉张作用控制下。39°N以北地区最大主压应力轴方位最优解显示一定角度的偏转，同时最大、中等、最小主压应力轴最优解推断的应力状态由西向东存在一个正断层—走滑断层—正断层的转换过程;而39°N以南地区的现今构造应力场保持稳定，最优主压应力轴呈NEE—SWW向，大部分网格应力状态显示走滑型。构造应力场的反演结果与活动构造、GPS主应变方向和剪切波分裂的快波偏振方向等相关研究结果基本一致。
关键词：京津冀地区;震源机制;应力张量阻尼反演;地壳应力场
中图分类号：P315.727文献标识码：A文章编号：1000-0666（2020）04-0610-10
0引言
地壳应力场的时空分布特征，特别是地壳深部应力状态，是解决地震发生机理、区域地壳稳定性、震源物理过程等地球动力学有关科学问题的基础。地质调查与大地测量结果表明，华北平原现今不存在大规模的拉张变形，区内主要断裂带的滑移速率低（Shen et al，1996;郭良迁等，2008），是一个地壳变形相对较弱的地区。在这种缓慢的变形区域中，沉积盖层或浅部岩石中的直接应力测量通常不能代表深部的应力场，而震源机制解是可用于估计缓慢变形区域中震源深度处应力场的少数信息之一。
京津冀地區是我国东部强震活动频度和强度最强的地区，同时也是新构造活动最强烈的地区之一，区内分布有张家口—渤海地震构造带（张渤地震带）和华北平原地震构造带，中强地震绝大多数都发生在这2条构造带或其交会部位。该区是国内使用地震记录开展地壳应力场较早的地区，李钦祖等（1973）利用河北红山台和沙城台记录的小地震P波初动资料确定了台站周边区域的应力场;许忠淮等（1979，1983）对上述方法进行改进，形成了“综合震源机制解法”，求解了京津唐张和华北地区的构造应力场方向。首都圈数字地震台网于2002 年开始正式运行（刘瑞丰等，2008），为地震精确定位和震源机制解反演提供了有利条件，研究人员利用P波初动极性、振幅比和波形反演求解了大量的中小震震源机制解，且其质量也更加可靠（胡新亮等，2004）。基于此，不同研究人员使用不同方法开展了包括首都圈在内的京津冀地区的地壳应力场研究（兰从欣等，2005;李瑞莎等，2008;张红艳等，2009;武敏捷等，2011，2012;胡幸平，崔效锋，2013;黄骥超，万永革，2015;刘静等，2016;刘丽等，2017;樊文杰等，2019）。其中大部分研究根据构造特征和地震分布将研究区域分成几个应力小区，每个应力小区的范围都和许忠淮等（1983）对京津冀地区的分区基本一致，所不同的是张红艳等（2009）在此基础上进行了细分，但是所得的构造应力场方向基本相同;刘静等（2016）和刘丽等（2017）对河北地区整体构造应力场的研究结果也验证了该地区处于华北地区统一地壳应力场的控制下（李钦祖，1980;魏光兴等，1982;许忠淮等，1983）。
在京津冀地区震源机制解数据日益丰富的当下，武敏捷等（2011）采用细分网格加滑动步长的方式初步讨论了华北北部地区的构造应力场，由于每个网格都是独立反演，分区方式不同和分区大小的变化都可能对反演的应力张量结果产生影响（Maury et al，2013）。相邻2个网格主应力方向的不同或许反映了应力场真实的变化，但也可能是震源机制数据的误差或反演的约束较差造成的假象。地壳中的真实应力场是连续分布的，为了解决应力反演模型依赖于分区的问题，更准确地反映研究区域的应力场空间变化特征，Hardebeck和Michael（2006）构建了一组可调整的阻尼参数，引入平滑约束来抑制相邻网格应力模型的差异，提出了区域尺度的应力张量阻尼反演方法。本文基于京津冀地区及邻区的中小地震震源机制解，将该方法应用于对该区精细地壳应力场特征的研究。
1数据和方法
1.1震源机制解数据
本文所用震源机制解共2 187个，震级区间2.0≤ML≤5.5，覆盖范围包括京津冀地区以及山西、辽宁、内蒙部分地区。数据来源分为2个部分：第一部分为2002—2008年京津冀地区ML≥2.0地震的震源机制解，使用梁尚鸿等（1984）提出的利用区域地震台网垂直向P，S最大振幅比资料测定小震震源参数的方法，该方法反演结果的可靠性在前人的对比论证中得到了证实（胡新亮等，2004;付虹等，2009），且已经在首都圈的应力场研究中得到应用（李瑞莎等，2008;张红艳等，2009;武敏捷等，2011，2012）;第二部分为2009—2019年10月京津冀地区ML≥2.0地震的震源机制解，使用格点尝试法（许忠淮等，1983）利用P波初动极性计算震源机制解，P波初动极性资料主要来自“地震编目系统”下载的震相报告和从事件波形中人工读取的P波初动。为了保证震源机制解的可靠性，参与计算的P波初动极性数据一般要求10个以上，且在震中周围均匀分布，矛盾比小于0.20。胡新亮等（2004）对首都圈数字地震台网内的地震分别使用格点尝试法和垂直向最大振幅比方法测定震源机制，二者结果基本一致，表明2种方法计算结果可以联合起来进行京津冀地区地壳应力场特征分析。
1.2震源机制解分类方法
对震源机制解分类有多种方法，如使用P，T轴参数进行三角形分类（Frohlich，1992，2001）、根据震源机制解3个应力轴倾角的世界应力图划分标准（Zoback，1992）、以3个应力轴的倾角最大值判断震源机制破裂类型（刁桂苓等，2011）等，以上都是以应力轴参数对震源机制解分类。本文采用Shearer 等（2006）提出的依据节面滑动角ri对震源机制解进行分类的方法。该方法的优势在于提供单个标量断层类型值用于表征震源错动类型，标量断层类型值fptype为：fptype=r90（1）式中：当r1 1.3应力场反演方法 目前，使用震源机制解进行应力场反演应用最广泛的2种方法为：震源机制应力反演（FMSI）的网格搜索方法（Gephart，Forsyth，1984;Gephart，1990）和自助线性应力反演（LSIB）方法（Michael，1984，1987）。Hardebeck和Hauksson（2001）使用含噪声的合成数据对这2种方法进行测试，2种方法都能精确地确定应力方向;FMSI方法通常对于高质量数据更准确，而LSIB方法对于离散度大的数据更准确;LSIB方法产生的置信区间通常较为合理，而FMSI方法的置信区间通常过大。合成数据测试结果表明，只要滿足应力场均匀和断层的滑动方向与剪应力方向一致的模型假设，由LSIB方法产生的置信区间也应适用于真实数据。Hardebeck和Michael（2006）引入阻尼最小二乘反演方法对LSIB方法加以改进，提出了区域尺度的应力张量阻尼反演方法SATSI（Spatial And Temporal Stress Inversion）算法;Martínez-Garzón等（2014）将SATSI算法移植到Matlab环境中，形成一款结合可靠经典方法的、新简化用户处理及可视化工具的应力反演Matlab软件包。由于本文所用震源机制解数据量较大，且2种不同的方法计算得到的数据离散度较大，故选用MSATSI软件包反演京津冀地区地壳应力场，以期可以准确地约束应力张量的方向（Hardebeck，Michael，2006）。 2震源机制解统计特征 本文参考Lin和Okubo（2016）的做法，对京津冀地区中小地震震源机制解的破裂类型比例进行定量估计，取Xf=0.33，将[-1，1]分成3个相等的部分，震源机制解也相应分为走滑断层、正断层和逆断层（图1）。对上述地震进行统计后发现，京津冀地区有1 315个走滑型震源机制解，占整体数量的60.1%，628个正断型地震，占28.7%，244个逆断型地震，占11.2%，该区震源错动类型主要以走滑型和正断型地震为主。不同类型的震源机制解空间分布（图1）显示，走滑型地震在整个区域弥散分布，而正断型和逆断型地震基本上沿着张渤地震带和华北平原地震带分布，没有形成明显的单一类型的丛集区。 将震源机制解按节面走向、倾角、滑动角和P轴方位角绘制成直方图进行统计（图2）。由于绝大多数地震震级太小，无法区分真正的断层面和辅助断层面，故将2 个节面等同看待，参数合并统计。节面走向大体呈各个方位均匀分布，在NNE—SSW和NWW—SEE向稍显优势（图2a）;节面倾角以接近直立的高倾角最多，表明以走向滑动为主，几乎没有近水平的节面（图2b）;节面滑动角以走向滑动为主，左旋、右旋走滑数量相仿，正倾滑明显多于逆倾滑（图2c）。P轴方位角在NEE—EW和SWW—EW向存在显著优势取向，可以看作是京津冀地区最大水平压应力方向（图2d）。上述震源机制解的统计特征与前人研究结果（李钦祖，1980;许忠淮等，1983;武敏捷等，2012;刘静等，2016;刘丽等，2017）基本一致。 3地壳应力场反演 采用1°×1°的步长，将京津冀地区（36.0°～43.0°N，113°～120°E）划分网格，使用MSATSI软件包反演该地区精细的地壳应力场。该方法通过在反演过程中加上一个合理的阻尼，使反演的相邻区域的应力张量变化最小，消除人为划分反演区域导致的应力偏转假象，更好地凸显整体应力场特征（Hardebeck，Michael，2006）。在反演之前要先确定阻尼系数e，它控制着理论值与观测数据之间拟合残差和应力反演模型长度的相对权重，对反演结果起着至关重要的作用。图3为模型长度—数据拟合残差的折中曲线，其拐点在1.2，因此最佳阻尼系数取e=1.2。反演中每个网格采用的震源机制解至少为6个，随机选取其中一个节面进行构造应力场反演。反演应力场参数的置信水平的可设定范围为68%～95%，本文设定为95%。重采样的迭代次数对于解的不确定度评价很重要，重采样次数过少导致评估结果没有统计意义，过多会导致计算量大，反演中重采样次数的默认范围是1 000～5 000次（Martínez-Garzón et al，2014），本文对所有数据进行2 000次抽样估计。经过计算得到最优状态下的3个应力主轴的方向和相对应力大小（R值）（Gephart，Forsyth，1984）为：R=（σ2-σ1）/（σ3-σ1）（4）式中：σ1，σ2，σ3分别表示最大、中间、最小主压应力。 从表1和图4可以看出，最优主压应力轴σ1的优势方向为NEE—EW向，且有大约半数的网格倾角接近直立，表明应力状态为正断层类型，而网格编号05，09，23的最优主压应力轴σ1方向为NNW—NWW向，倾角都接近直立，表明最大主应力轴垂直，这与前人用中强震震源机制解和综合震源机制解研究华北地区地壳应力场的结果有所差异（李钦祖，1980，1982;魏光兴等，1982;许忠淮等，1983），主要表现为最优主应力轴的倾角直立或水平，这可能与本文所用的震源机制解多为小震、易受局部构造运动控制、且具有随机性特征有关，反映了整体背景应力场下的非均匀性和复杂性;最优主张应力轴σ3的优势方向为SSE—NWW向，且倾角均接近水平。 从图4中最大主压应力σ1轴的指示方向可以看出，在39°N以北的首都圈地区最大主压应力轴自西向东存在由NE向近EW向的偏转，这与前人对首都圈构造应力场的研究结果基本一致（张红艳等，2009;武敏捷等，2012;黄骥超，万永革，2015;樊文杰等，2019），而且本文的分区更加精细，应力场的变化连续可靠。由最大、中间、最小主压应力轴最优解推断的张渤地震带陆地段的应力状态表现为，自西向东存在一个正断层—走滑断层—正断层的转换过程，这与张红艳等（2009）、武敏捷等（2012）、樊文杰等（2019）以走滑型为主的应力场特征存在一定差异;同时，在黄骥超和万永革（2015）对首都圈的构造应力场反演结果中，北京—张家口—大同地区和唐山及邻区都表现出少量正断层性质。39°N以南地区现今构造应力场保持稳定，最优主压应力轴σ1为NEE—SWW向，大部分网格应力状态显示走滑型。 相对应力大小R值是中间主应力σ2轴的值，是更接近于最大主压应力轴σ1还是最小主压应力轴σ3的大小的度量。当R值趋于中间值0.5时，表示3个应力主轴方位可完全分辨;而当R值趋于取值端点0或者1时，表示压缩或拉张应力状态完全不能分辨，只能确定它在垂直于拉张或压缩应力主轴的平面内。一般认为当R值小于0.5时，获得的应力状态偏拉张性质，反之则为偏压缩性质（Guiraud et al，1989）。从表1和图5可以看出，所有网格的R值均小于0.5，表明京津冀地区的应力状态呈现偏拉张性质。同时从图4可知，最大主应力轴σ1的不确定范围最大，最大主压应力轴σ1和中间应力轴σ2表现的压应力状态是一致的，两轴无法区分（Guiraud et al，1989;万永革等，2011）;而最小主压应力轴σ3的不确定度变化范围相对稳定，我们可以推断京津冀地区现今地壳应力场处于一个相对统一的NNW—SSE向的拉张作用控制下。4结论和讨论 基于京津冀地区及邻区2 187个中小地震震源机制解，通过计算标量断层类型值将其划分为3种震源错动类型，采用MSATSI软件包反演该地区1°×1°网格的精细地壳应力场，研究结果表明：研究区震源机制类型以走滑断层和正断层为主，占比约90%;不同类型的震源机制解没有形成明显的单一类型的丛集区，反映出小地震的发生除了受背景应力场的控制，还受局部构造条件控制且具有一定的随机性。这为我们利用小地震的震源机制解来推断构造应力主轴的方向提供了更大的优越性，震源机制解的一对节面解中既含有构造应力场的信息，也可能含有局部介质不均匀结构的影响（许忠淮等，1983）。 地壳应力场反演结果显示，所有网格的相对应力大小R值均小于0.5，表明京津冀地区应力状态偏拉张性质，而且最小主压应力轴σ3的不确定度变化范围相对稳定，表明现今京津冀地区地壳应力场处于一个相对统一的NNW—SSE向的拉张作用控制下。这一结果与前人的研究结果比较一致，高名修（1979）认为华北地区处于引张应力场环境，有限元模拟显示主张应力大约为主压应力2～8倍（刘峡等，2010）。利用反演结果中最大、中间、最小主压应力轴最优解推断的应力状态显示，39°N以北地区的最大主压应力轴方向存在一定角度的偏转，同时应力状态由西向东存在一个正断层—走滑断层—正断层的转换过程;而39°N以南地区的现今构造应力场保持稳定，大部分网格应力状态显示走滑型。 39°N以北地区与前人研究中的张渤带陆地段和首都圈范围基本一致，大致分为3段：西段、中段和东段。西段为张家口—大同地区，应力状态表现为正断层类型，该区为山西断陷带右旋走滑在NE端形成的一个张性构造区（徐锡伟等，2002），活动构造、遥感和垂直形变研究结果都表明山西断陷带中北段至张家口地区总体以张性正断层活动为主（孙启凯等，2018;吴玉涛等，2018;申星等，2019;彭远黔等，2019），该区的中小地震震源机制解以正断型和走滑型为主（孙小入等，2019）。中段为北京—天津地区，应力状态为走滑类型和正断层，该区历史上曾发生过1679年三河—平谷地震，使用现今小震资料推断该地震的错动类型为走滑兼少量正断（王晓山，2017），2018年发生在该区的河北永清地震的震源机制为走滑兼正断类型（王晓山等，2018;李赫等，2020），GPS观测资料亦显示该区存在高剪切变形区（陈长云，2016）。东段为唐山地区，应力状态为走滑类型和正断层，杨雅琼等（2016）对唐山地震序列的现今小震震源机制分段拟合最佳应力张量，得到的结果为正断层和走滑类型兼有。综合前人多种资料的论证表明，本文的研究结果与前人结果具有较好的一致性，但仍存在少许差异，可能由所使用的反演方法、网格划分和数据量等导致，有待进一步研究分析。 GPS主应变率结果显示：京津冀大部分地区地壳变形特征主要受NEE的主压应变控制（陈长云，2016），京津冀地区的剪切波分裂的快波偏振方向为NEE—EW向（吳晶等，2008a，b;赵博等，2011;吴鹏等，2017）。该结果与该区使用中小地震震源机制解反演的构造应力场最优解基本一致，同时也和前人使用中强地震震源机制解和综合震源机制解得到的华北地区一致性良好的应力场吻合，表明中小地震震源机制的结果也可以较好地用于区域构造应力场的研究。京津冀地区是中国大陆较早实现数字化地震观测网络的地区，目前积累了丰富的中小地震波形资料，开展中小地震震源机制解及构造应力场的动态演化研究可以为该区地震动力学研究提供具有物理意义的信息，更好地服务于该地区的防震减灾综合研究。 参考文献： 陈长云.2016.张家口—渤海断裂带分段运动变形特征分析[J].地震，36（1）：1-11. 刁桂苓，徐锡伟，陈于高，等.2011.汶川MW7.9和集集MW7.6地震前应力场转换现象及其可能的前兆意义[J].地球物理学报，54（1）：128-136，doi：10.3969/j.issn.0001-5733.2011.01.014. 樊文杰，崔效锋，胡幸平，等.2019.首都圈地区震源机制解及现今构造应力场时空变化特征研究[J].地震学报，41（1）：33-45，doi：10.11939/jass.20180064. 付虹，王绍晋，龙晓帆，等.2009.利用数字化速度波形和仿真位移波形资料求解震源机制解的比较[J].地震研究，32（1）：18-24. 高名修.1979.华北断块构造区的现代引张应力场[J].地震地质，1（2）：1-12. 郭良迁，马青，杜雪松，等.2008.华北地区断层形变与地震的关系[J].大地测量与地球动力学，28（3）：14-20. 胡新亮，刁桂苓，马瑾，等.2004.利用数字地震记录的P，S振幅比资料测定小震震源机制解的可靠性分析[J].地震地质，26（2）：347-354. Gephart J W，Forsyth D W.1984.An improved method for determining the regional stress tensor using earthquake focal mechanism data：application to the San Fernando Earthquake sequence[J].J Geophys Res，89（B11）：9305-9320. Gephart J W.1990.FMSI：A FORTRAN program for inverting fault/slikenside and earthquake focal mechanism data to obtain the regional stress tensor[J].Computers and Geosciences，16（7）：953-989. Guiraud M，Laborde O，Philip H.1989.Characterization of various types of deformation and their corresponding deviatoric stress tensors using microfault analysis[J].Tectonophysics，170（3-4）：289-316. Hardebeck J L，Hauksson E.2001.Stress orientations obtained from earthquake focal mechanisms：what are appropriate uncertainty estimates？[J].Bull Seis Soc Amer，91（2）：250-262. Hardebeck L L，Michael A J.2006.Damped regional-scale stress inversions：methodology and examples for southern california and the coalinga aftershock sequence[J].J Geophys Res，111（B11）：B11310，doi：10.1029/2005JB004144. Lin G，Okubo P G.2016.A large refined catalog of earthquake relocations and focal mechanisms for the island of Hawaii and Its seismotectonic implications[J].J Geophys Res，121（7）：5031-5048，doi：10.1002/2016JB013042. Martínez-Garzón P，Kwiatek G，Ickrath M，et al.2014.MSATSI：A MATLAB package for stress inversion combining solid classic methodology，a new simplified user-handling，and a visualization tool[J].Seismol Res Lett，85（4）：896-904. Maury J，Cornet F H，Dorbath L.2013.A review of methods for determining stress fields from earthquakes focal mechanisms;application to the sierentz 1980 seismic crisis（upper Rhine graben）[J].Bull Soc geol，184（4-5）：319-334. Michael A J.1984.Determination of stress from slip data：faults and Folds[J].J Geophys Res，89（B13）：11517-11526. Michael A J.1987.Use of focal mechanisms to determine stress：a control study[J].J Geophys Res，92（B1）：357-368. Shearer P M，Prieto G A，Hauksson E.2006.Comprehensive analysis of earthquake source spectra in Southern California[J].J Geophys Res，111（B6）：303，doi：10.1029/2005JB003979. Shen Z K，Jackson D D，Ge B X.1996.Crustal deformation across and beyond the Los Angeles Basin from geodetic measurements[J].J Geophys Res，101（B12）：27957-27980. Zoback M L.1992.First- and secondorder patterns of stress in the lithosphere：the world stress map project[J].J Geophys Res，97（B8）：11703-11728，doi：10.1029/92JB00132. Characteristics of Crustal Stress Field in Beijing-Tianjin-Hebei Region
