刘珠妹 李盛乐

1 中国地震局地震研究所（地震大地测量重点实验室），武汉市洪山侧路40号，430071

2015－04－25发生的尼泊尔Mw7.8（M8.1）地震是尼泊尔81a来遭遇的最强烈地震。此次地震强余震不断，共发生Mw5 以上余震17 次，且地震破裂由西向东迁移，对我国西藏等地区造成一定的影响。

大地震发生后，对后续强余震发生时间的研究是震后趋势判定的工作重点之一。谷继成等［1］对我国1931～1977年28个地震序列进行研究，提出地震序列中相邻余震之间“等待时间”规律，并以此判断地震类型，以及进行大震后强余震发生时间的短临预测。但在实际工作中，该方法更多地被用于地震序列类型的判定［2－4］，而较少应用于强余震发生时间的预测［5－8］。本文首先从理论上对“等待时间”强余震预测法的残差结构改变问题进行分析，在此基础上给出了基于加权最小二乘拟合的模型改进方法，并将其应用于尼泊尔地震的强余震预测。

1 问题提出与模型改进

1.1 “等待时间”强余震预测模型

“等待时间”即一个地震序列中相邻两组强余震之间的时间间隔。谷继成［1］发现，主余型地震序列中第i与第i－1个强余震的等待时间Δti与第i个强余震的发震时刻ti满足以下关系：

根据已知的地震余震序列，采用最小二乘法即可求得式（1）的线性拟合参数β0、β1，进而通过最近一次强余震的发震时间ti预测后续余震的发震时间。

1.2 基于加权最小二乘的强余震预测模型

“等待时间法”对观测数据采用了对数变换，使其满足线性模型假设。但同时，对数变换也改变了模型残差的方差性质。在最小二乘求解线性方程时，改变了最小二乘法的准则函数［9］。

设Δt为相邻强余震的等待时间，t为强余震的发震时间，则“等待时间－发震时间”模型的普通最小二乘准则函数为：

而对数变换后的最小二乘准则为：

从而

忽略常数项，对数变化后的最小二乘准则函数变为：

将原始“等待时间－发震时间”数据作对数变换相当于对原始数据进行加权，且权值为。加权的结果导致等待时间Δt值较小时，权值较大；随着能量的衰减，相邻两次强余震等待时间增大时，权值较小。因此利用普通最小二乘法求解拟合参数时，最近的余震样本距主震发震时间较远，权值较低，对下一次强震时间预测控制力不足，从而影响预测的精度。

为抵消上述加权结果，在对数变换后线性拟合时，同样可对数据加权：

令权值wi＝Δti，即可在一定程度上弥补双对数变化对余震预测误差的影响，从而提高预测的精度。利用加权改进的“等待时间法”进行强余震预测的流程见图1。

2 尼泊尔地震强余震预测

2.1 余震目录准备

谷继承［11］建议，对于7 级以上强震，一般筛选低于主震2.5～3级的余震作为已知样本进行后续预测，确定合适的序列震级下限。本文分别针对尼泊尔地震Mb4.5、Mb5以上余震目录进行“等待时间”分析。数据采用美国ANSS 地震目录，时 间 截 至05－12 07：05（UTC 时）发 生 的Mw7.3强余震。

图1 “等待时间法”强余震模型预测流程Fig.1 Flow of waiting time prediction for strong aftershocks

图2（a）中余震下限为Mb4.5时，序列在“发震时间－等待时间”双对数坐标轴中的线性拟合决定系数R2为0.68，样本均方误差MSE为0.26，余震序列分布较分散，不适宜进行后续余震时间预测；而图2（b）中，5级以上余震序列线性拟合决定系数R2为0.95，MSE为0.06，线性关系更为显著，序列分布更集中。因此，选择Mb5以上尼泊尔余震目录进行“等待时间”强余震预测实验。

表1 给出了尼泊尔Mw7.8 地震序列中Mw5.0以上地震目录，并计算每个余震样本距离主震的发震时间ti、等待时间Δti和各自的对数值。其中序号13的余震与后续余震的间隔时间（0.28h）小于前一个强余震等待时间（7.88h）的1／20，文中将其视为同一次能量释放，取较大震级的余震参与分析［11］。

修正目录前后，样本的线性关系发生明显变化。图3（a）中，未筛选目录拟合系数R2为0.532，修正后R2提高到0.873，达到0.01的显著水平（图3（b））。

图2 不同余震下限样本的线性关系Fig.2 The linear relationship of samples with different minimum magnitudes

2.2 实验分析

选取整个序列中通过线性回归R检验及F检验，表现为显著线性关系的N个样本进行多组强余震预测实验。首先对样本分别进行普通最小二乘拟合和以Δti为权值的加权最小二乘拟合。将线性拟合参数β0、β1代入式（1），得到拟合线性方程。再穷举Δt＝i（i＝1，2，3，…，n），tN＋1＝tN＋ΔtN－1／20＋i，作出穷举线（图4）。穷举线与拟合线的交点即为预测的第N＋1个强余震，交点x坐标值经乘幂变换后即为强余震的发震时间。

表1 尼泊尔Mw7.8地震5级以上余震序列统计表Tab.1 Mw7.8Nepal earthquake sequences with Mw≥5.0

图3 筛选目录前后样本的线性关系Fig.3 The linear relationship of samples before and after sequence sieving

图4 穷举法求解余震发震时间（样本数N＝10）Fig.4 Solving equation by exhaustion methods（sample number＝10）

表2中“最大样本时间”指参与建模的最后余震距主震的发震时间。除N＝11／13的两组实验外，其他5组经过加权改正的余震预测精度均有不同程度的提高。总体平均绝对残差由原模型的0.173降为0.16，模型预测精度得到改善。预测的置信区间跨度也由原模型的0.590 缩小为0.125，说明预测结果的不确定性降低，预测实用性有了一定的提高。

3 分析与讨论

为进一步验证加权改正模型的适用性，本文另选取几组典型的7级以上主余型地震的强余震样本进行模型改正前后的强余震预测对比实验。表3给出不同震例分别采用普通最小二乘和加权最小二乘得到的后续余震预测值、预测残差以及真实发震时间与实际误差。

由表3 及图5 可见，经过改进的“等待时间法”应用于多个震例后，模型预测精度均有不同程度的提高。精度提高最大的乌恰震例，原模型预测误差为0.323，模型改正后预测误差减小为0.014。芦山余震预测误差，也由0.112 减小为0.071。6 个震例的平均绝对残差由原模型的0.24降为0.045，平均降幅高达71%。

虽然经过模型改正后，余震预测的误差有了不同程度的降低，但由对数模型转化为真实时间后，实际预测误差却长达数10h，甚至上百h。由表3可见，于田、通海等震例改正模型的预测误差均小于0.1，而转化为实际时间后，后续余震的预测发震时刻和实际发震时刻相差长达4～5d。

表2 原模型和加权改正的“等待时间”余震预测结果统计Tab.2 Statistical results of aftershock predictionwith two methods

表3 典型震例的“等待时间”法强余震预测分析结果统计Tab.3 Statistical results of aftershock prediction among typical earthquake instances

图5 多震例的模型改正前后预测误差比较Fig.5 Error comparation of predicting model before and after correction

图6将样本距主震时间分组为“1d以内”、“3 d以内”、“7d以内”、“15d 以内”以及“15d 以上”，分析模型误差、实际误差与样本天数的关系。由图中可知，改进后的预测方法虽然模型误差控制在了较小的范围内（＜0.1），但随着发震时间的推移，3d以后实际时间误差可达数天，已不满足应急所需。由于“等待时间法”样本数要达到5个以上才具有稳定的拟合结果，对于短时间内发生较多强余震的主余型震例，利用本文方法进行强余震预测有效性较高；反之，对强余震发震较稀疏、时间跨度较大的震例而言，利用“等待时间法”进行余震预测的效率则不甚明显。

图6 模型误差和实际误差与样本时间的关系Fig.6 Relationship among model error、true error and sample time

［1］谷继成，谢小碧，赵莉.强余震的时间分布特征及其理论解释［J］.地球物理学报，1979（1）：32－46（Gu Jicheng，Xie Xiaobi，Zhao Li.On Temporal Distribution of Large Aftershocks of the Sequence of a Major Earthquake and Preliminary Theoretical Explanation［J］.Chinese Journal of Geophysics，1979（1）：32－46）

［2］蒋海昆，黎明晓，吴琼，等.汶川8.0级地震序列及相关问题讨论［J］.地震地质，2008（3）：746－758（Jiang Haikun，Li Mingxiao，Wu Qiong，et al.Features of the May 12 M8.0 Wenchuan Earthquake Sequence and Discussion on Relevant Problems［J］.Seismology and Geology，2008（3）：746－758）

［3］华爱军，刁守中，王红卫.强余震“等待时间”判别方案预报效能的研究［J］.地壳形变与地震，1996（2）：47－52（Hua Aijun，Diao Shouzhong，Wang Hongwei.Study on Prediction Ability of Method for Distinguishing the Waiting Time of Strong Aftershock［J］.Crustal Deformation and Earthquake，1996，16（2）：47－52）

［4］孟令媛，周龙泉，张小涛，等.2014年新疆于田Ms7.3地震序列和震源特征初步分析［J］.中国地震，2014（2）：159－167（Meng Lingyuan，Zhou Longquan，Zhang Xiaotao，et al.Study the Characteristics of the Yutian，Xinjiang Ms 7.3 Earthquake，February 12，2014［J］.Earthquake Research in China，2014（2）：159－167）

［5］曲延军.新疆部分地震的余震序列特征及强余震预报［J］.内陆地震，1990（1）：87－96（Qu Yanjun.The Characteristics of Aftershock Sequence and the Prediction of Strong Aftershock in Some Parts of Xinjiang［J］.Inland Earthquake，1990（1）：87－96）

［6］秦乃岗.广东及邻区地震频度衰减系数和余震等待时间［J］.华南地震，1990，10（4）：42－49（Qin Naigang.The Method of Earthquake Frequency Attenuation Coefficient and Aftershock Duration Time Applied in Guangdong［J］.South China Seismological Journal，1990，10（4）：42－49）

［7］赵小艳，韩立波，徐甫坤，等.2014年云南鲁甸6.5级地震序列跟踪分析研究［J］.地震研究，2014（4）：508－514（Zhao Xiaoyan，Han Libo，Xu Fukun，et al.Research on Tracking Analysis for Ludian Ms6.5 Earthquake Sequence in Yunnan in 2014［J］.Journal of Seismological Research，2014（4）：508－514）

［8］平建军，刘荣环，贾炯，等.地震序列较强余震灰色及最小二乘拟合预测方法的应用研究［J］.华北地震科学，2005（1）：6－13（Ping Jianjun，Liu Ronghuan，Jia Jiong，et al.The Application Study of the Predict Method of the Gray and the Method of Least Square in the Earthquake Array Strong Aftershock［J］.North China Earthquake Sciences，2005（1）：6－13）

［9］耿鸿江.对数变换的问题及其改进［J］.水文，1996（2）：13－16（Geng Hongjiang.Problems in the Logarithmic Transformation and Improvements［J］.Journal of China Hydrology，1996（2）：13－16）

［10］杨安华，彭清娥，贺昌政，等.双对数模型对模型模拟误差的放缩问题探讨［J］.数学的实践与认识，2006（8）：135－139（Yang Anhua，Peng Qinge，He Changzheng，et al.The Discussing of the Problem of Enlarging or Reducing Simulation Error by Double Logarithmic Model［J］.Mathematics in Practice and Theory，2006（8）：135－139）

［11］谷继成.地震类型的判别与强余震的预测［J］.地震，1987（6）：1－9（Gu Jicheng.A Method of Distinguishing Earthquake Type and of Predicting Large Aftershocks［J］.Earthquake，1987（6）：1－9）