地壳应力场的时空分布特征,特别是地壳深部应力状态,是解决地震发生机理、区域地壳稳定性、震源物理过程等地球动力学有关科学问题的基础。地质调查与大地测量结果表明,华北平原现今不存在大规模的拉张变形,区内主要断裂带的滑移速率低(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)构建了一组可调整的阻尼参数,引入平滑约束来抑制相邻网格应力模型的差异,提出了区域尺度的应力张量阻尼反演方法。本文基于京津冀地区及邻区的中小地震震源机制解,将该方法应用于对该区精细地壳应力场特征的研究。

本文参考Lin和Okubo(2016)的做法,对京津冀地区中小地震震源机制解的破裂类型比例进行定量估计,取X_f=0.33,将[-1, 1]分成3个相等的部分,震源机制解也相应分为走滑断层、正断层和逆断层(图1)。对上述地震进行统计后发现,京津冀地区有1 315个走滑型震源机制解,占整体数量的60.1%,628个正断型地震,占28.7%,244个逆断型地震,占11.2%,该区震源错动类型主要以走滑型和正断型地震为主。不同类型的震源机制解空间分布(图1)显示,走滑型地震在整个区域弥散分布,而正断型和逆断型地震基本上沿着张渤地震带和华北平原地震带分布,没有形成明显的单一类型的丛集区。

图1 京津冀地区中小地震震源机制解类型分布
Fig.1 The spatial distribution of focal mechanism solution types in the Beijing-Tianjin-Hebei region

将震源机制解按节面走向、倾角、滑动角和P轴方位角绘制成直方图进行统计(图2)。由于绝大多数地震震级太小,无法区分真正的断层面和辅助断层面,故将2 个节面等同看待,参数合并统计。节面走向大体呈各个方位均匀分布,在NNE—SSW和NWW—SEE向稍显优势(图2a); 节面倾角以接近直立的高倾角最多,表明以走向滑动为主,几乎没有近水平的节面(图2b); 节面滑动角以走向滑动为主,左旋、右旋走滑数量相仿,正倾滑明显多于逆倾滑(图2c)。P轴方位角在NEE—EW和SWW—EW向存在显著优势取向,可以看作是京津冀地区最大水平压应力方向(图2d)。上述震源机制解的统计特征与前人研究结果(李钦祖,1980; 许忠淮等,1983; 武敏捷等,2012; 刘静等,2016; 刘丽等,2017)基本一致。

采用1°×1°的步长,将京津冀地区(36.0°～43.0°N,113°～120°E)划分网格,使用MSATSI软件包反演该地区精细的地壳应力场。该方法通过在反演过程中加上一个合理的阻尼,使反演的

图2 震源机制解节面参数统计直方图
Fig.2 Histograms are shown for the nodal plane paremeters of focal mechanism solution

相邻区域的应力张量变化最小,消除人为划分反演区域导致的应力偏转假象,更好地凸显整体应力场特征(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=(σ₂-σ₁)/(σ₃-σ₁)(4)

图3 应力场反演模型长度—数据拟合残差的折中曲线
Fig.3 Trade-off curve between model length and data fitting misfit in the stress field inversion modle

式中:σ₁,σ₂,σ₃分别表示最大、中间、最小主压应力。

从表1和图4可以看出,最优主压应力轴σ₁的优势方向为NEE—EW向,且有大约半数的网格

表1 京津冀地区应力场反演结果(只列出有结果的网格)
Tab.1 Results of stress field inverted in the Beijing-Tianjin-Hebei region(only give the grids with results)

倾角接近直立,表明应力状态为正断层类型,而网格编号05,09,23的最优主压应力轴σ₁方向为NNW—NWW向,倾角都接近直立,表明最大主应力轴垂直,这与前人用中强震震源机制解和综合震源机制解研究华北地区地壳应力场的结果有所差异(李钦祖,1980,1982; 魏光兴等,1982; 许忠淮等,1983),主要表现为最优主应力轴的倾角直立或水平,这可能与本文所用的震源机制解多为小震、易受局部构造运动控制、且具有随机性特征有关,反映了整体背景应力场下的非均匀性和复杂性; 最优主张应力轴σ₃的优势方向为SSE—NWW向,且倾角均接近水平。

从图4中最大主压应力σ₁轴的指示方向可以看出,在39°N以北的首都圈地区最大主压应力轴自西向东存在由NE向近EW向的偏转,这与前人对首都圈构造应力场的研究结果基本一致(张红艳等,2009; 武敏捷等,2012; 黄骥超,万永革,2015; 樊文杰等,2019),而且本文的分区更加精细,应力场的变化连续可靠。由最大、中间、最小主压应力轴最优解推断的张渤地震带陆地段的应力状态表现为,自西向东存在一个正断层—走滑断层—正断层的转换过程,这与张红艳等(2009)、

图4 1°×1°网格划分的应力场反演结果
Fig.4 Results of stress field inversion on 1°×1° grid

武敏捷等(2012)、樊文杰等(2019)以走滑型为主的应力场特征存在一定差异; 同时,在黄骥超和万永革(2015)对首都圈的构造应力场反演结果中,北京—张家口—大同地区和唐山及邻区都表现出少量正断层性质。39°N以南地区现今构造应力场保持稳定,最优主压应力轴σ₁为NEE—SWW向,大部分网格应力状态显示走滑型。

相对应力大小R值是中间主应力σ₂轴的值,是更接近于最大主压应力轴σ₁还是最小主压应力轴σ₃的大小的度量。当R值趋于中间值0.5时,表示3个应力主轴方位可完全分辨; 而当R值趋于取值端点0或者1时,表示压缩或拉张应力状态完全不能分辨,只能确定它在垂直于拉张或压缩应力主轴的平面内。一般认为当R值小于0.5时,获得的应力状态偏拉张性质,反之则为偏压缩性质(Guiraud et al,1989)。从表1和图5可以看出,所有网格的R值均小于0.5,表明京津冀地区

图5 1°×1°网格划分的相对应力大小R值(网格中心矩形的颜色表示这个网格的R值,中括号内的数值范围表示95%置信度下R值的范围)
Fig.5 The relative stress magnitude R on 1°×1° grid(the color-filled squares in each grid represent the corresponding R values,the numbers in the brackets give the 95% confidence intervals of the solutions)

的应力状态呈现偏拉张性质。同时从图4可知,最大主应力轴σ₁的不确定范围最大,最大主压应力轴σ₁和中间应力轴σ₂表现的压应力状态是一致的,两轴无法区分(Guiraud et al,1989; 万永革等,2011); 而最小主压应力轴σ₃的不确定度变化范围相对稳定,我们可以推断京津冀地区现今地壳应力场处于一个相对统一的NNW—SSE向的拉张作用控制下。

基于京津冀地区及邻区2 187个中小地震震源机制解,通过计算标量断层类型值将其划分为3种震源错动类型,采用MSATSI软件包反演该地区1°×1°网格的精细地壳应力场,研究结果表明:研究区震源机制类型以走滑断层和正断层为主,占比约90%; 不同类型的震源机制解没有形成明显的单一类型的丛集区,反映出小地震的发生除了受背景应力场的控制,还受局部构造条件控制且具有一定的随机性。这为我们利用小地震的震源机制解来推断构造应力主轴的方向提供了更大的优越性,震源机制解的一对节面解中既含有构造应力场的信息,也可能含有局部介质不均匀结构的影响(许忠淮等,1983)。

地壳应力场反演结果显示,所有网格的相对应力大小R值均小于0.5,表明京津冀地区应力状态偏拉张性质,而且最小主压应力轴σ₃的不确定度变化范围相对稳定,表明现今京津冀地区地壳应力场处于一个相对统一的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)。该结果与该区使用中小地震震源机制解反演的构造应力场最优解基本一致,同时也和前人使用中强地震震源机制解和综合震源机制解得到的华北地区一致性良好的应力场吻合,表明中小地震震源机制的结果也可以较好地用于区域构造应力场的研究。京津冀地区是中国大陆较早实现数字化地震观测网络的地区,目前积累了丰富的中小地震波形资料,开展中小地震震源机制解及构造应力场的动态演化研究可以为该区地震动力学研究提供具有物理意义的信息,更好地服务于该地区的防震减灾综合研究。