地震b值来源于古登堡-里克特关系式logN=a-bM(Gutenberg and Richter, 1941)，式中M代表震级，N代表震级大于震级M的地震个数，地震b值直接反映了大小地震之间的比例关系.从实验室得到岩石破裂应力与b值的关系(Scholz, 1968)开始，大量观测结果(Stefan and Benoit, 1996; Schorlemmer et al., 2005; Murru et al., 2007)都充分表明了b值与差应力成负相关关系.另外，地震b值与地震活动规律也关系密切(Mogi, 1967; Smith, 1981; Frohlich and Davis, 1993; Öncel et al., 2001)即大地震发生的区域常对应低b值区域(Schorlemmer et al., 2004; Nakaya, 2006; Nuannin et al., 2012)，这些结果都可以用b值与差应力成负相关的规律来解释.因此在一定意义上把b值当作地壳的“应力计”来使用是可行的，可以帮助分析在实际工作中测量较为困难的应力状态.

2008年5月12日在龙门山地区发生了震级达8级的汶川地震，必然会对龙门山地区的应力状态产生重大影响，且5年后又发生M_S7.0级的芦山地震，进一步影响龙门山地区的应力状态，本文期望通过研究该地区b值的变化趋势研究应力的演化过程.前人对龙门山地区的b值也做过很多研究(Zhao and Wu, 2008; 王俊等，2009；田文，2013；Zhang and Zhou, 2016)，结果差别较大.我们认为b值的计算方法、方法中参数的选取、b值在深度上的差异、地震定位精度等都可能影响到b值的准确性，因此本文首先对b值的计算方法和参数选择进行深入讨论，然后采用绝对定位和相对定位相结合的方法对龙门山地区的地震进行了重新定位，在此基础上进行较为精确的b值计算并进一步探讨b值所揭示的应力演化过程及其构造意义.

b值的计算方法主要有最小二乘法和最大似然法两种，最小二乘法通过直线拟合M-logN的关系，得到直线的斜率继而得到b值；最大似然法是由地震的概率密度函数得来，其公式为这里M表示震级的平均值，M_c表示起算震级(Aki, 1965).

首先利用合成完全精确的地震目录来评估两种算法的优缺点及其影响因素.取b=1，最大震级为6级，即利用关系式logN=-M+6，也即M=6-logN，N从1到10000分别取值，从而计算得到10000条精确且完整的地震震级目录.利用该完整地震目录，计算不同起算震级条件下的b值，发现随着起算震级的增大，最大似然法逐渐偏离真实值，而最小二乘法对精确数据而言，直线拟合斜率精确为-1，即b=1.这是因为最大似然法是由概率密度函数近似得来，潜在的对样本数量要求较高，当地震个数过少时，对最大似然方法的影响较大.如图 1所示，当地震个数为100个时(M_c=4.0) 其误差为3.3%，而地震个数为10个时(M_c=5.0)，其误差甚至达26.2%，可见最大似然法对地震个数非常敏感.特别是在做b值的空间扫描时，格子内的地震数目常在20~50个之间，个数较少，这种情况最大似然法在格子内计算的精度会受到很大影响.

图 1

Fig. 1

图 1 完全精确数据时不同起算震级下最小二乘和最大似然结果比较 圆点线为最大似然的结果，方格点线为最小二乘的结果. Fig. 1 The comparison between least square and maximum likelihood method under different threshold magnitude for accurate synthetic data The circles show the result of maximum likelihood method and the squares are the result of least square method.

然后同样利用合成数据，分别对目录的完整性，震级测定误差，震级四舍五入误差等分别对两种算法进行了测试.发现最大似然法具有很好的抗误差能力，但对目录的完整性比较敏感，因为最大似然方法作为概率方法，对地震个数特别敏感，小地震个数多，如果缺失会造成较大影响；而最小二乘方法对大地震和小地震的权重比较接近，因而对目录的完整性没有最大似然法敏感，但对大地震的震级误差较敏感，并且一些特大地震可能已经不遵守G-R定律.针对最小二乘法的此问题，本文提出最小地震数约束的方法，即要求震级间隔内最少的地震个数大于某个给定值(本文取15)，如图 2为大地震存在较大震级误差且都偏大的极端情况，采用最小地震数约束的方法仍然可以获得较为可靠的计算结果.在震级加入5%的随机误差的一般情况下，新方法也具有更好的结果.本文采用改进的最小二乘方法计算b值，同时深度剖面中把最大似然方法的结果作为参考.

图 2

Fig. 2

图 2 传统最小二乘法(a)和最小地震数约束的最小二乘法(b)结果比较 两图是在大震级存在较大误差的特殊情况下的比较，改进的方法能够大幅度提高结果的可靠性. Fig. 2 Comparison between the traditional least square method (a) and this improved method (b) constrained by minimum earthquake number The two graphs show comparison for the special case of large magnitude with large error. The new method significantly improved the reliability of result.

由于我们希望计算b值的深度变化，因此对地震震源深度的要求较高，本文首先利用VELEST震源位置和速度模型联合反演的方法(Kissling et al., 1994)同时得到地震的绝对位置和最小一维速度模型，然后再将二者用于相对定位的双差重定位HYPODD算法中，从而获得较为精确的地震位置.前人从结构的角度也对龙门山地区地震进行了重新定位(黄媛等，2008；朱艾斓等，2008；陈九辉等，2009)，但时间序列较短，本文对1992年到2014年的全部2级以上地震进行了重新定位.图 3为本文获得的最小一维速度模型和其他两种模型的比较，其中使用反演得到的最小一维速度模型来进行重定位其得到的结果残差值为0.163 s，相较于赵珠模型(赵珠等，1997)重定位残差值为0.173 s和人工地震(嘉世旭等，2014)重定位残差值0.218 s更小.

图 3

Fig. 3

图 3 三种速度模型的比较 灰色实线为本文联合反演的结果，黑色实线为赵珠模型 (赵珠等，1997)，点线为人工地震模型(嘉世旭等，2014). Fig. 3 Comparison of three 1D velocity models The gray real line is our model obtained by joint inversion, black real line is the Zhao′s model (Zhao et al., 1997) and the dotted line is the model obtained from active source (Jia et al., 2014).

定位前后震源水平位置没有太大的变化，而震源深度在重定位前后变化范围达10 km左右(图 4).可见重定位的工作是很有必要的.

图 4

Fig. 4

图 4 重定位前后震源深度的变化 Fig. 4 Changes of focal depth before and after relocation

使用重定位后的结果进行b值的计算，双差定位法不能把全部的地震都定出来，这里定出约70%的地震，剩余约30%的地震位置由VELEST绝对定位方法给出.

分析汶川地震前后各时间段内的b值拟合情况，发现汶川地震前，起算震级取2级时，地震目录较完整，且起算震级取2.5最小二乘的结果与起算震级取2最小二乘的结果差别很小，范围在0.01~0.08；汶川地震后，2008年5月的数据发现明显的不完整，起算震级取到4级时，地震目录才完整，这可能是因为地震后台站遭到破坏及大量余震数据重叠而很难识别造成的.所以在计算b值随时间的变化时，2008年5月起算震级取为4级，而计算b值的空间变化时，去掉了2008年5月不完整的地震数据；对于2008年5月以后各时间段内的地震数据进行分析发现当起算震级取2时，地震目录较完整，且与起算震级取2.5最小二乘的结果相差很小.

我们把所有地震分别向垂直于断层带和沿断层带方向投影，对深度剖面进行b值扫描，扫描格子大小分别为20 km×10 km，50 km×10 km，步长为0.5 km.同时要求格子内地震数最少为20个，且最大震级和最小震级的差大于1.5级.图 5，图 6分别为垂直断层带，汶川地震前后b值随深度的变化，因为投影采用全部地震，地震个数较多，最小二乘和最大似然结果比较接近.

图 5

Fig. 5

图 5 垂直于断层带方向汶川地震前b值随深度的变化 (a)最小二乘结果；(b)最大似然结果.横坐标为相对龙门山断裂带中央(104.5°E, 32°N)的距离，纵坐标表示深度，黑色五角星为主震位置.冷色调为低b值，暖色调为高b值. Fig. 5 Variation of b-value in depth for all earthquakes projected to perpendicular profile to the Longmenshan fault zone before the Wenchuan Earthquake The horizontal coordinate is the distance in kilometer from the center (104.5°E, 32°N) of the Longmenshan fault zone and the vertical coordinate is the depth in kilometer, the black star represents the location of the main shock. The cold color represents low b-value and the warm color indicates high b-value. Figure (a) and (b) show the result of least square and maximum likelihood, respectively.

图 6

Fig. 6

图 6 垂直于断层带方向汶川地震后b值随深度的变化 图中坐标、颜色均与图 5相同. Fig. 6 Variation of b-value same as figure 5 but for the earthquakes after the Wenchuan Earthquake The two coordinates and colors are same as figure 5.

我们可以看到b值随深度的变化情况具有明显的成层性，由于汶川地震后的数据量较汶川地震前要多很多，因而其图像反映更加明显，具体表现为以大约15 km为界，15 km以上的区域为高b值区，15 km以下的区域为低b值区.这一现象反映了b值可以揭示脆性韧性层的差异.在15 km以上岩石呈脆性，表现为b值较高；15 km以下的岩石呈韧性，表现为b值较低，于是在b值随深度变化的图上可以看到明显的分层结构.这一结果与实验室得到的围压与b值之间的结论一致，即随着岩石围压的增大，b值逐渐降低，这是因为随着围压的增大岩石逐渐趋于韧性的结果(Amitrano，2003).

从图 5，6可以看出，b值随深度变化很大，同时大部分地震都发生在15 km以上，因此在计算b值的时空变化时只考虑15 km以上的区域，如果采用全部地震来计算b值，则难以分清b值的变化来源于深部还是浅部.

将汶川地震前后15 km以上的地震分为几个时间段分别计算b值，b值随时间的变化趋势如图 7所示，在汶川地震前逐渐降低，并在汶川地震后达到最低，然后逐渐恢复，直到芦山地震发生，芦山地震后有略有降低.Zhang和Zhou(2016)也发现了类似的变化趋势.这种变化趋势可能反映了龙门山断裂带应力积累、释放的演化过程.

图 7

Fig. 7

图 7 15 km以上地震b值随时间的变化 Fig. 7 Temporal change of average b-value with error bar for earthquakes above 15 km

根据b值与差应力成反比的实验室观测(Scholz, 1968)，汶川地震前b值出现下降，可能显示应力不断积累的过程，汶川地震后即2008年5月的b值明显降低(同震变化)，可能反映局部的差应力增大，因为很多余震都在震后短时间发生，并且大多发生在应力增加区(解朝娣等，2010)由该区域的地震计算得到的b值也自然反映应力增加区的应力状态.2008年6月之后b值缓慢恢复，反映差应力逐渐下降并低于震前的水平.2013年4月20日发生芦山地震后，b值又出现一个小幅下降.全部数据计算得到的b值时间变化反映的是整个区域平均的变化趋势，进一步计算不同时间段的b值空间分布才能看到汶川地震和芦山地震震源区b值和应力演化及两强震间的相互作用.

对龙门山地区(102.5°E—106°E，30.5°N—33°N)进行b值的空间扫描，格子大小为0.5°×0.5°，步长为0.05°，要求格子内地震数目大于20个，最大最小震级差大于1.5级，可以得到龙门山地区b值的空间分布图，图 8a所示为汶川地震前(1992年6月到2008年5月)b值的空间分布图，可以看到红框内的区域是明显的b值低异常区，反映了该区域的高差应力状态，这一区域的位置正好对应于映秀附近的彭灌杂岩区，可见彭灌杂岩区为应力高度集中的区域，这可能是震后最大位错发生于该区域的根本原因.如果将龙门山断裂带分成南北两个部分分别做其汶川地震前的b值，进行比较可以看到汶川地震前北段b值为1.45，南段为1.12，北段b值明显高于南段的b值，反映了北段的差应力水平要低于南段的差应力水平，这与龙门山断裂带北段以走滑断层为主，南段以逆冲断层为主的构造特征相一致，通常逆断层差应力水平高而b值低(Schorlemmer and Wiemer, 2005).

图 8

Fig. 8

图 8 b值空间分布图随时间的变化 (a) 1992年6月到2008年4月；(b) 2008年6月到2008年12月；(c) 2011年到2013年4月19日；(d) 2013年4月20日到2014年；图中黑色线条为断层，深紫色线条为地表破裂，大小品红色五角星分别为汶川地震震中和芦山地震震中，空心圆点为震级大于3级的地震，图a中红色方框区域为彭灌杂岩区大致的范围，冷色调为低b值，暖色调为高b值. Fig. 8 Spatial distribution map of b-value in four time periods Figure a, b, c, d represent b-value in four time periods, June 1992 to April 2008 before Wenchuan earthquake, June 2008 to December 2008 after Wenchuan earthquake, Jan 2011 to April 19, 2013 before Lushan earthquake, April 20, 2013 to Dec 2014 after Lushan earthquake, respectively. The black and deep purple lines indicate active faults and surface rupture, respectively. The big and small magenta stars are the epicenter of Wenchuan earthquake and the Lushan earthquake. The circles are earthquakes for magnitude largerthan 3.0, and the red box in Figure (a) indicates the Pengguan complex rock. The color is same as figure 5.

不同时间段内b值空间分布图的变化也能反映4.2节中同震变化和震后断层愈合的过程，如图 8所示为4个时间段(1992年6月到2008年4月、2008年6月到2008年12月、2011年到2013年4月19日、2013年4月20日到2014年)内b值空间分布图的变化.由于2008年5月的地震目录在起算震级取4时才完整，地震个数过少，无法计算出b值空间分布图，但根据2008年5月总体的b值可知这一时间段内b值相较其他时间段要偏低.

从图 8中可以看出，汶川地震后映秀和北川两处位错最大的区域b值先升高，而其之间的区域b值略有降低，而后至芦山地震前整个地震带的b值都达到一个较高的水平.可以看出映秀北川地区即震后破裂最大的地区b值上升速度最快，反映这里在震后应力调整最快.芦山地震后在芦山震源附近区域b值有所降低，可能与汶川地震同样在震后短时间内出现局部差应力增加的原因一致，但几乎在整个汶川地震破裂区都出现了b值降低的趋势，这意味着芦山地震的发生促进了汶川地震破裂区应力的增加，说明两地震之间存在直接的联系，即两地震之间的空段可能是软弱区.因为如果该空段非常坚硬，则南部芦山地震的发生不可能对很远的北部地区的应力情况产生影响，只有该空区是软弱区，当青藏高原作为一个整体向东推挤的过程中，南部芦山地震发生后，青藏高原向东运动会增加，才会造成北部的汶川破裂区挤压增加，产生应力增加，b值降低的现象.汶川地震和芦山地震之间的空段根据速度结构(Pei et al., 2014)的结果已经发现为韧性变形区，本文对b值时空变化研究，进一步证实了该结果，并且发现芦山地震起到了促进汶川破裂区的愈合的作用.

本文主要对龙门山断裂带的地震b值进行了研究，通过震源位置和速度模型联合反演的方法得到了最优一维速度模型和绝对位置，并用此速度模型对龙门山断裂带的地震进行了双差重定位的相对约束，然后根据重定位结果精确计算了b值的时空分布，获得了汶川地震前后不同时间段内b值大小的变化趋势及b值空间分布规律，并进一步探讨了b值所揭示的龙门山断裂带的结构构造特征、应力环境变化与断层愈合过程.

主要结果与结论有：

(1) b值在深度上的分层性明显，15 km以上b值较高，15 km以下b值较低，这主要反映上部岩石为脆性，下部岩石为韧性的特征；汶川地震前15 km以上断裂带北段b值远高于南段，反映以走滑断层为主的北段差应力低于以逆冲断层为主的南段.

(2) 汶川地震前的b值分布揭示了区域应力分布，在映秀附近的彭灌杂岩区呈现低b值异常，意味着该区域存在较大的差应力，这可能是震后最大位错发生在该区域的根本原因.

(3) 汶川地震前后区域平均b值的时间演化，揭示了龙门山断裂带的应力演化过程，汶川地震破裂区在震前b值出现下降，而后逐渐增加并超过震前水平，且在芦山地震之后该区又发生明显b值降低，这揭示了汶川地震断层的同震变化和愈合过程.分时段b值空间分布结果显示，在芦山地震之后，汶川地震破裂区b值下降，即应力增加，这意味着汶川芦山地震之间的空段存在韧性变形，且芦山地震的发生促进了汶川地震破裂区的愈合.