青藏高原在西南受印度板块的挤压、在北部受西伯利亚刚性块体的抵阻而形成大量褶皱和断裂带^[1], 各个主断裂带将青藏高原东南部划分成多个部分, 其地表形变并非一致。青藏高原东南部岩石圈变形异常强烈, 地震活动性非常活跃, 地质构造特征和应力分布复杂多变, 因此也是地震危险性评估的难点地区。

震级-频度经验关系公式^[2]为

$\log _{10} N=a-b M$,

(1)

式中：N指震级M大于等于M_c的地震累计次数, M_c为最小完整性震级; 地震b值可反映区域内应力的分布与变化。已有的岩石力学实验、数值模拟实验和实际数据分析等都验证了b值与区域应力存在负相关关系^[3-6]。根据二者的负相关关系, 利用地震b值可估计区域内应力大小。许多研究发现大地震前后在破裂区附近的b值变化显著, 反映区域内震前期应力累积和同震应力释放过程^[7-15]。地震b值能指示区域地震活动及应力大小分布情况, 因此, 对地震活动性分析和地震危险性评价具有重要意义。

青藏高原东南部近年来发生多次强震活动, 地震数据记录比较完备, 为该应力场研究提供了丰富的地震资料。前人通过地震目录研究b值获得大量研究结果。刘雁冰和裴顺平^[7]通过研究汶川地震区发现, 在汶川地震和芦山地震后b值出现下降, 然后逐渐升高, 在映秀附近地区震前期存在低b值异常; 刘艳辉等^[10]则将汶川地震和芦山地震前后b值的变化具体化, 提出b值-时间曲线在震前震后存在水平—负增长—正增长—负增长—水平的变化, 认为汶川地震前川北存在低b值的异常区; 史海霞等^[11]发现汶川地震前b值在较长时间内存在趋势性降低的现象, 汶川地震前半年快速显著下降; 钱晓东和秦嘉政^[12]根据小江断裂带周边地震活动和各区域b值, 认为小江断裂带东川段、石屏—建水段和宜良—嵩明段应力高; 易桂喜等^[13]认为川滇东部b值空间分布差异大, 并指出小江断裂带、安宁河—则木河断裂带和鲜水河断裂带上应力较大的区域; 刘静伟等^[14]分析研究鲜水河断裂带中段存在低异常b值区, 认为该区域具有较高应力水平; 朱艾斓等^[15]根据b值分布指出鲜水河断裂带、安宁河断裂带和则木河断裂带上可能存在的应力较大的闭锁段和应力水平低的蠕滑段位置。

为进一步研究青藏高原东南部地震活动特点、b值及应力时空分布特征, 本文以青藏高原东南部为研究区域, 采用美国地震学研究联合会(Incorporated Research Institutions for Seismology, IRIS)提供的1970—2019年的地震目录, 对研究区域地震活动分布进行分析, 利用地震b值估计青藏高原东南部的应力分布及其变化特征, 为地震危险性评估和抗震减灾等工作提供参考。

1 青藏高原东南部构造背景

青藏东南部是中国大陆地震活动最频繁的地区之一。该地区地质构造和深部结构复杂, 活动断裂带主要包括鲜水河—小江断裂系统与龙门山断裂带, 鲜水河—小江断裂系统自北向南发育为西北—东南向转化成近南北向, 龙门山断裂带发育为东北—西南向, 以上断裂带将青藏东南分为川滇块体、华南块体和巴颜喀拉块体3个主要的构造块体(图 1)。由于青藏高原物质东流与高原的隆起, 青藏东南地区地壳运动绕喜马拉雅东构造结呈顺时针旋转^[16-18], 华南块体运动相对稳定, 川滇块体地壳运动方向由近WE旋转转为NS方向。

鲜水河—小江断裂带是大型左旋走滑断裂系统, 主要包含鲜水河断裂带、安宁河断裂带、则木河断裂带、大凉山断裂带和小江断裂带。鲜水河—小江断裂带地震活跃, 是青藏高原活动最剧烈的断裂带之一, 20世纪以来发生多次7级以上强震^[19-20]。近年来该断裂带发生2次较大地震, 2014年11月22日M_W 6.1康定地震发生于鲜水河断裂带东南段, 2014年8月3日M_W 6.2级鲁甸地震发生于大凉山断裂带西侧(图 1)。

龙门山断裂带是巴颜喀拉块体与稳定的华南块体的分界、川西震区与相对刚硬的四川盆地的分界, 走向NE-SW方向(图 1)。虽然龙门山断裂带现今的活动速率较小, 但该断裂带孕育了2008年5月12日M_W 7.9汶川地震和2013年4月20日M_W 6.6芦山地震(图 1)。汶川地震在龙门山断裂系统的中央断裂和前山断裂上分别形成了长度250 km和70 km的地表破裂, 最大垂直和右旋水平错动分别达6.5 m和4.8 m, 沿破裂带平均错动达2~3 m ^[21-23]。

2 b值计算方法与地震数据 2.1 b值计算方法

目前, b值计算方法主要有2种, 即最小二乘法和最大似然法。前者是对震级-频度关系log₁₀N-M分布进行最小二乘线性回归估算b值。而后者最大似然法计算结果相对稳定可靠且更常用, 该方法由Aki^[24]最先提出：

$b=\frac{1}{\ln 10\left(\bar{M}-M_{\mathrm{C}}\right)}$,

(2)

其中：M是平均震级, M_C为最小完整性震级, 指对于某较小震级M_C, 震级M≥M_C的所有地震都能被检测和记录。为提高b值的计算精度和稳定性, 本文采用Utsu^[25]对公式(2)的修正公式

$b=\frac{1}{\ln 10\left(\bar{M}-\left(M_{\mathrm{C}}-\Delta M / 2\right)\right)}$,

(3)

其中: ΔM为震级划分的最小间隔, 本文选取ΔM=0.1。式(3)将M_C替换为M_C -ΔM/2, 因为M_C -ΔM/2更接近真实的最小完整性震级。本文采用最大曲率法确定M_C。

b值的不确定度常用δb来评估：

$\delta b=b / \sqrt{N}$

(4)

此处N为样本量。

2.2 地震数据

本文的地震数据资料于IRIS的网站(https://www.iris.edu/hq/)下载。该网站提供了全球范围内1970年至今连续的地震资料, 数据资源时间跨度广、涵盖震级范围大, 相对全面丰富。选用IRIS网站公布的地震目录, 首先是考虑到该网站地震数据的可获取性和可共享性, 也便于其他研究者验证我们的结果。IRIS的地震资料包含的信息有：地震发生时间、震源经纬度、深度、震级和震级类型(m_b、M_L、M_S和M_W)、地震发生城市和数据来源等。为确保计算结果的准确性, 将地震资料中的震级类型统一化。下载的地震数据中绝大多数使用的震级类型为m_b, 为减少震级换算带来的误差, 故将各类震级统一换算转化为m_b类型的震级, 采用的是李莹甄等^[26]的经验公式

$m_{\mathrm{b}}=0.71 M_{\mathrm{L}}-0.009 M_{\mathrm{L}}^2+1.86$

(5)

$m_{\mathrm{b}}=0.56 M_{\mathrm{S}}+2.57$,

(6)

$m_{\mathrm{b}}=0.73 M_{\mathrm{W}}+1.37$

(7)

为调查青藏高原东南部的地震资料的完整性及地震分布特征, 统计了该区域1970—2019年地震数据的震级-时间和震级-深度分布情况(图 2)。从图 2可以看出, 该区域1970—1980年4级以下地震记录缺失比较严重; 1980—1995年3级以下记录有一定缺失; 1995年1月—2015年8月, 2~3级以上地震记录完整性最好; 而在2015年之后, IRIS未记录4级以下的小地震(图 2(a))。青藏东南的地震主要分布在0~100 km深度范围, 5.5级以上的地震分布于深度50 km以内, 7级以上的大地震分布于深度40 km甚至更浅的深度(图 2(b))。

3 b值的空间分布

青藏高原东南部的震源深度主要集中于20 km以上的脆性上地壳^[27-28], 参考王椿镛等^[29]的地壳速度结构的研究, 并结合青藏高原东南部地震目录资料的震源深度特征(图 2(b)), 深度大于40 km的地震数目较少, 为充分探究该区域b值和应力的空间分布特征, 将地震事件按震源深度分为3层：0~20、20~40和大于40 km。将不同深度层的地震事件按经纬度投影到水平剖面, 计算每个网格点的b值和M_C, 可得到整个水平剖面上的b值和M_C分布(图 3)。

青藏高原东南部的水平方向b值显示(图 3(a))：区域西南部b值比东北部的b值略高, 龙门山断裂带附近b值较低, 推测龙门山断裂带附近所受应力相对较大; 安宁河断裂带北部b值较高, 推测此处附近所受应力相对较小; 则木河断裂带、大凉山断裂带和小江断裂带交汇处及其东侧b值较低, 推测此处所受应力较大。另外, 沿龙门山断裂带走向可识别出汶川地震与芦山地震之间的地震空区, 该空区的b值比其周围的b值高; 安宁河—则木河断裂带及其西侧的b值高于大凉山断裂带附近及其东侧的b值, 推测安宁河—则木河断裂带附近的应力小于大凉山断裂带周边应力。

将深层20~40 km(图 3(b))与浅层0~20 km(图 3(a))结果对比, 在龙门山断裂带附近地区两层b值都比较低; 四川盆地和安宁河—则木河断裂带西侧的区域在深层的b值比浅层小得多, 说明应力水平随深度的增加而增大。

在青藏高原东南部, 震源深度＞40 km的地震数目较少(图 2(b)), 搜索半径R在最大值内无法搜索出足够的地震数目N, 所以未能计算深度＞40 km的水平剖面上的b值。

图 3(c)和3(d)分别显示了在0~20 km和20~40 km剖面上的M_C分布。浅层的M_C大部分在3.7~3.9, 局部地区(小江断裂带西侧)在3.9~4.0;深层的M_C在3.8~4.0。表明浅层的地震事件完整性比深层完整性好。

由于使用的地震数据可能存在震源深度定位的误差, 因此我们也测试了不同分层深度范围, 调查不同的分层深度对b值分布结果的影响(图 4和图 5)。在原分层深度范围(0~20 km和20~40 km, 图 3)的基础上, 计算了不分层(所有深度, 图 4)及其他测试深度范围的b值分布, 包括0~16和16~40 km、0~18和18~40 km(图 5)。

不分层的结果(图 4(a))与0~20 km层的结果(图 3(a))较为接近, 是由于该区域的震源深度主要集中于20 km以上的脆性上地壳^[27-28]。浅层(0~16/0~18/0~20 km)的结果比较接近, 均可反映青藏高原上地壳脆性层中的相对高应力(b值相对较低)和相对低应力(b值相对较高)区域(图 5(a)~5(c))。对比深层(16~40/18~40/20~40 km)的结果(图 5(d)~5(f)), 随着数据深度增加, b值有降低的趋势。这与前人关于围压由浅至深逐渐增大、岩石由脆性逐渐转化为韧性会得到低b值的实验结果^[30]相一致。该结果也与图 3的结果一致。

4 b值的时间变化

近些年来青藏高原东南部发生4次较大的地震：2008年汶川地震、2013年芦山地震与2014年的鲁甸地震和康定地震(图 1, 图 6(a))。通过b值-时间变化(图 6(b))与震级-时间关系(图 6(a))对照, 可以发现大地震发生时该区域b值明显降低且处于较低水平。汶川地震前的2002—2008年, 该区域的b值有持续降低的趋势, 应与区域内的构造应力一直不断累积有关; 2008—2010年, 汶川地震发生后区域b值呈升高趋势; 2010—2015年区域b值持续下降, 与2013年芦山地震和2014年的鲁甸地震、康定地震的孕育有一定的相关性; 而更早的1970—1995年以及震后的2015年后, b值的变化明显、波动剧烈(图 6(b)), 可能与区域应力水平的变化有关, 但期间的地震记录缺失较严重, 导致计算的b值可靠性与准确性有所降低。1995—2015年期间的地震数据资料比较完备, 其δb明显小于其他时间段(图 6(b)), 选取该时间段内的地震数据估计研究区域震级-频度关系的b值(图 6(c))。该时间段青藏高原东南部b值估计在1.04左右, 最小完整性震级M_C为3.6。

分别对2008汶川地震、2013芦山地震、2014鲁甸地震和康定地震震中(图 1)附近的地震活动区域进行b值时间扫描(图 7)。考虑到以上地震发生的间隔, 为分别调查它们对震中及邻区的影响并保证有足够的地震事件数量, 以及便于跟前人研究结果^[10]对比讨论, 以一年为单位进行b值时间扫描。在2008年, 汶川、芦山和康定的b值都出现明显的下降, 汶川的b值降低尤为显著, 距汶川较近的芦山比距汶川略远的康定b值降低多, 而距汶川更远的鲁甸的b值在2008年波动较小, 甚至出现略微的上升趋势(图 7)。在2013年, 震中芦山和近处的汶川b值明显降低, 而较远的康定和鲁甸的b值受芦山地震的影响较小。2014年的鲁甸地震和康定地震的震中b值都在该年明显下降, 汶川和芦山的b值受这2次地震的影响较小, 一方面是由于汶川和芦山距鲁甸较远, 另一方面这2次地震震级在6级左右, 相对前2次的震级较小, 故影响也较小。综上所述, 震级越大对b值(应力)的影响范围越广, 中小震对b值(应力)的影响范围比较有限; 距震中越近, b值的降低越明显, 而距震中较远的区域甚至会出现b值升高的现象。马鸿庆^[31]和周翠英等^[32]在研究华北和山东强震前后地震活动性特征时也发现了类似震中的近场b值降低而远场b值升高的现象。

以汶川地震为例, 调查了大地震发生前后一段时间内地震活动及b值随时间变化情况。利用震源深度0~40 km的地震数据(图 8), 对汶川地震前后的地震活动和b值变化进行调查和分析(图 9)。在对汶川地震前后时间段划分时, 根据汶川地震前后地震活动性的变化, 在保证b值的稳定性以及满足网格点上地震数目N≥50的情况下, 尽可能详细地了解龙门山断裂带上b值分布的时间变化, 因而将2008年5月12日前后分为T1~T8时间段。

T1时间段, 汶川地震发生前3~5年, 鲜水河断裂带南端和安宁河断裂带北端以及小江断裂带南北端附近的地震活动略多(图 8(a)); 龙门山断裂带附近地震活动相对平静, 前人研究也发现汶川地震及其他大地震前存在“地震平静现象”^{[11, 33-39]}。由于没有足够数目的地震活动, 无法准确计算龙门山断裂带及其附近在该时间段的b值。

T2时间段, 龙门山断裂带附近的地震活动性有增加的趋势(图 8(b)); 但地震活动的数目仍不足, 故无法计算此段时间龙门山断裂带及其附近的b值分布。

T3时间段, 即震后20 d左右, 龙门山断裂带附近的地震活动急剧增加(图 8(c)), 汶川地震后发生大量的余震; 由于震后20 d的时间较短, 只有龙门山断裂带附近发生了大量的地震(图 8(c))。b值空间分布结果(图 9(a))只显示了龙门山断裂带附近的b值, 在0.7~1.2, 震中附近b值较低。

T4时间段, 龙门山断裂带附近地震活动仍保持较活跃的程度(图 8(d)), 但与上一阶段相比频率有所降低, 地震主要发生在龙门山断裂带的中部和东北部; 断裂带附近b值在0.9~1.4, T4时间段的b值(图 9(b))较T3的b值(图 9(a))整体有略微升高, 尤其是龙门山断裂带中段较T3时间段升高明显, 说明余震活动也释放了应力, 使龙门山断裂带的应力水平降低。

T5时间段, 地震活动较活跃(图 8(e)); 龙门山断裂带附近在该月b值(图 9(c))较上月明显升高, 在震中附近尤为明显, 推测此时龙门山断裂带应力水平仍在下降, 断裂带中部b值高而东北端b值低。

T6时间段, 龙门山断裂带附近地震活动(图 8(f))比上月降低; 此时震中附近100 km左右的范围内b值达到最高, 龙门山断裂带中段b值最高(图 9(d))。

T7时间段, 震后4~7个月龙门山断裂带附近的地震活动仍然明显高于其他地区(图 8(g)), 龙门山断裂带b值开始回落(图 9(e)), 断层面开始闭锁, 构造加载导致应力增加。

T8时间段, 地震活动聚集于龙门山断裂带中部至东北部(图 8(h)), 该裂带附近的地震活动仍然比青藏高原东南部的其他地区频繁; 龙门山断裂带附近b值继续降低, 应力累积。

T3~T8时间段, 龙门山断裂带的b值经历了先升高后降低的过程。b值升高反映汶川地震及大量余震使震前累积的应力得以释放, 在震后3个月左右b值达到最高即所受应力最小(图 9(d)), 此后b值开始回落(图 9(e)和9(f)), 应力开始逐渐累积增大。

5 讨论 5.1 b值空间分布与构造特征

龙门山断裂带附近b值较低(图 3~图 5), 推测是因为青藏高原隆升, 龙门山断裂带西北侧的巴颜喀拉块体自北西向南东推覆和挤压于华南块体之上, 受华南块体内四川盆地高强度块体的坚强抵阻^[40-41], 应力在逆冲的龙门山断裂带附近累积, 断层面受较强烈的挤压^[42-46]; 另外, 龙门山山脉在四川盆地的西北边缘挤压隆生, 高程载荷在龙门山断裂带附近施加额外的挤压, 造成该区域所受应力较大, b值较低。图 3(a)暗示汶川地震和芦山地震之间的地震空区应力较小。

大凉山断裂带附近b值低于安宁河—则木河断裂带(图 3, 图 4), 推测大凉山断裂带的所受应力更大。青藏东南部自西向东由NW-SE逆冲应力场变为NS拉张应力场的转变中, 大凉山断裂带承担着重要的作用, 分担着更多华南块体和川滇块体的相对运动, 侧面说明了大凉山断裂带在鲜水河—小江断裂带中段“截弯取直”的趋势和结果^[47-48]。

对比不同深度的水平剖面b值结果(图 3, 图 5)可以发现, 深层的b值比浅层的b值低, 反映出深层的应力水平大于浅层, 与实验室发现的b值与围压负相关的现象^[30], 以及随着围压的增大岩石的韧性增强的实验结果^{[30, 49-50]}相吻合。

5.2 b值时间变化与地震孕育

汶川地震前5~6年, 青藏高原东南部的b值一直存在降低趋势(图 6), 直至2008年该地区b值迅速而显著下降(图 6(b), 图 9), 与刘艳辉等^[10]估算的结果趋势接近。汶川地震及大量的余震发生后, 龙门山断裂带附近的b值在2~3个月内迅速显著升高, 达到峰值后逐渐降低(图 9)。由于震前处于应力累积的状态, 应力增大, b值降低; 大地震和余震的发生使累积的应力得以释放, 应力减小, b值升高; 震后由于断层闭锁和构造加载导致应力累积, 因而b值开始回落; 至下一次大地震发生前, 应力累积程度达最高、b值也再次达到低谷。其他区域(例如日本和中国台湾)的大地震发生前后b值变化也有类似的降低—升高—降低的过程^[51-53]。地震的发生伴随着应力的累积与释放, 应力变化会造成b值的波动, 大地震可以使b值显著升高^{[7-8, 54-55]}。

2013年芦山地震、2014年鲁甸地震和康定地震发生前, b值也存在下降趋势, 但较汶川地震的b值变化幅度小, 这是因为震级小释放的应力少, b值波动较小。对比汶川地震和芦山地震前b值的变化(图 6), 可以发现地震震级越大, 震前b值负增长持续时间越久, 地震发生后b值降低越显著。

5.3 b值时空扫描的不确定性

在本研究中, 图 6(b)1995—2015年的b值不确定度δb(公式(4))大致在0.05以内; b值空间分布的网格点上(图 3~图 5, 图 9)的不确定度δb大致不超过0.15。

1995年以前和2015年以后青藏高原东南部的地震记录不够完整, 早期(1995年以前)的记录在各个震级都有一定的缺失, 晚期(2015年以后)的记录则是4级以下的小地震缺失严重。青藏高原比地震频繁的板块交界处的地震台网密度低, 检测到小地震的能力相对略差。IRIS网站关于青藏高原东南部的地震记录的震级类型(M_W、m_b、M_L、M_S等)不统一, 计算b值前需要通过经验公式进行换算, 震级换算带来的误差也影响到b值计算的准确度。

本文没有对震源位置进行重新定位, 原因有两个方面：其一, 本文研究区域为青藏高原东南部, 空间尺度较大, IRIS网站的数据相对精度较高; 其二, 参考前人在该区域的工作^{[13, 56-58]}, 并结合本文研究区域的范围尺度, 故未对震源位置进行重新定位, 一般相对小尺度问题的相关b值估计是需要考虑震源位置重定位的^{[7, 15, 59]}。

6 结论

通过对青藏高原东南部地震活动的分析及b值的时空扫描研究, 得到主要结论如下：

1) 在鲜水河—小江断裂系统的中部, 大凉山断裂带附近的b值比安宁河—则木河断裂带低, 推测大凉山断裂带承担较多的应力, 是未来地震风险较高的断裂带。

2) 研究区浅层(0~20 km)比深层(20~40 km)的b值高, 表明浅层区域比深层应力小, 与浅层区域围压低、岩石倾向脆性破裂, 而深层区域围压高、岩石倾向韧性变形的特征一致的。

3) 汶川地震发生前后, 震源及周边区域的b值存在降低—升高—降低的过程, 揭示了区域内应力的累积—释放—累积的演化过程。震级越大, b值波动越大, 震前b值降低趋势持续越久。

4) 龙门山断裂带及其附近地区的b值较低。巴颜喀拉块体自北西向南东推覆于华南块体之上, 应力在逆冲的龙门山断裂带附近累积, 断层面受较强烈的挤压; 另外, 青藏东部的高程载荷也对龙门山断裂带施加了额外的挤压, 因此造成该区域应力较大。此外, 汶川地震和芦山地震之间的地震空区应力较小。