人工地震目录的评估及其在青藏高原东北缘的应用

赵文涛，罗纲,2*，靳锡波，孙云强

1 武汉大学测绘学院，武汉 430079 2 武汉大学地球空间环境与大地测量教育部重点实验室，武汉 430079 3 福建农林大学交通与土木工程学院，福州 350002

0 引言

完备的地震目录是进行地震预测预报研究的重要基础.现代地震观测的历史仅有一百年左右，远小于多数大地震的复发周期；而古地震目录存在数据不完整及不精确等问题.因此，现有的地震目录难以用于挖掘大地震的规律和特点(Kagan and Jackson，1991；Goes，1996；李红等，2015).一种较为有效的解决以上问题的方法是使用长时间尺度的完备的人工地震目录(Robinson et al.，2011).

人工地震目录是基于对地震发生机制的认识，受已有的地震、地形变等观测数据的约束，可以通过多种方法模拟得到.Burridge和Knopoff(1967)提出了弹簧-滑块模型(BK模型).朱元清和石耀霖(1991)采用了串联多个耦合非线性单元的动力学模型.李丽等(1997)运用了弹簧-滑块-阻尼器组合单元非线性动力学模型.Bak和Tang(1989)给出了二维细胞自动机模型.朱守彪等(2006)改进了三维细胞自动机模型.我们也研发了可模拟三维断层系统人工地震目录的有限元数值模型(Luo and Liu，2018；孙云强和罗纲，2018；孙云强等，2019；Sun et al.，2020；Gao et al.，2020).

基于数值模拟得到的人工地震目录已被广泛应用，但很少有研究将人工地震目录与实际观测得到的地震目录进行对比分析.这导致人工地震目录的准确度或与观测的古地震数据的匹配度方面的评价方法较为缺乏.因此，本文开发了平均绝对误差法和余弦相似度法来尝试解决人工地震目录的准确度评定问题，并将其应用于我们前期模拟得到的青藏高原东北缘断层系统的地震目录(孙云强和罗纲，2018；孙云强等，2019)(图1)；通过将模拟得到的6个算例的人工地震目录与过去近9千年的区域古地震数据进行匹配对比(图2)，获取了两者之间的匹配度，给出了人工地震目录的准确度信息；最后，使用匹配较好的人工地震目录计算了大地震迁移概率.

1 有限元模型

本文使用的有限元模型是我们前期开发的青藏高原东北缘复杂断层系统三维动力学模型(孙云强和罗纲，2018；孙云强等，2019；Sun et al.，2020).该模型包含了青藏高原东北缘断裂系统中八条主要断裂，分别是海原断裂、香山天景山断裂、烟筒山断裂、牛首山断裂、贺兰山断裂、罗山断裂、黄河断裂和云雾山小官山断裂(图1).模型大小为500 km×500 km×100 km.该模型在深度方向上分为20 km厚的孕震的弹塑性上地壳和80 km厚的黏弹性的中下地壳上地幔.该两层的杨氏模量分别为8.25×1010Pa和1.1×1011Pa，泊松比均为0.25.

图1 青藏高原东北缘的断层、地震活动及有限元模型(修改自孙云强和罗纲，2018)(a)青藏高原东北缘的断层与地震活动.图中的黑色虚线矩形框是有限元模型的范围，红色圆圈代表震级在5级以上的历史地震，地震数据来自国家地震科学数据中心(https:∥data.earthquake.cn/)；震源机制解数据来自GCMT(http:∥www.globalcmt.org/CMTsearch.html).(b)三维黏弹塑性有限元模型的网格和边界条件.(c)有限元模型的断层系统.Fig.1 Faults and seismicity in northeastern Tibetan Plateau，and the finite-element model (Modified from Sun and Luo，2018)(a)Faults and seismicity in northeastern Tibetan Plateau.Black dashed rectangle is the domain of the finite-element model.Red circles are the locations of historic earthquakes with magnitude greater than 5 (from National Earthquake Data Center，https:∥data.earthquake.cn/).Focal mechanism solutions are from GCMT (http:∥www.globalcmt.org/CMTsearch.html).(b)Mesh and boundary conditions of the three dimensional visco-elasto-plastic finite-element model.(c)Fault system in the finite-element model.

我们的模型边界条件为：模型顶面是自由表面；底面法向位移固定，水平方向自由；侧面法向(水平方向)的边界速度由震间GPS速度场插值得到，而其剪切方向(垂直方向)是自由的.模型在重力和GPS速度边界条件下进行加载，模拟得到了近10万年的长期地震目录，其中模型前5万年的模拟计算是应力演化进入稳态加载状态的过程，不参与本文的匹配分析(孙云强等，2019).在模拟过程中，模型用应变软化的塑性单元模拟断层，使用Drucker-Prager塑性屈服准则判断地震的发生.当断层单元的应力积累达到屈服极限时，降低该断层单元的内聚力，从而产生同震滑动及相应的变形与应力；当断层单元在新的较低的内聚力下达到平衡时，地震结束；此时，将该断层单元的内聚力恢复到初始值，模型在断层单元上又开始累积应力，进入到下一次地震的震间加载阶段.上述屈服的断层单元可以是一个或多个.此过程可以重复，因此，形成地震循环及地震在断层系统中各个断层上的迁移.关于数值模型的控制方程、地震活动模拟方法等的详细描述，请参见我们的前期研究(Luo and Liu，2010，2012，2018；孙云强和罗纲，2018；孙云强等，2019；Sun et al.，2020；Gao et al.，2020).

2 数据

2.1 古地震数据

本文使用的古地震数据来自闵伟等(2000)的研究.他们研究了青藏高原东北缘的海原断裂、中卫-同心断裂(即香山天景山断裂的中东部分(俞岗等，2013))、罗山东麓断裂和贺兰山东麓断裂等主要断裂的古地震活动.这四条主要断裂在过去近9千年共发生了24次大地震(图2)(闵伟等，2000).其中，海原断裂发生大地震的次数最多，为10次；中卫-同心断裂、罗山东麓断裂和贺兰山东麓断裂发生大地震的次数较少，分别为5次、4次和5次.

图2 青藏高原东北缘的古地震数据(数据来自闵伟等，2000)(a)—(d)中的红色、绿色、蓝色和紫色的矩形分别表示海原断裂(HY fault)、香山天景山断裂(XT fault)、罗山断裂(LS fault)和贺兰山断裂(HLS fault)上的古地震.矩形的长度为包含时间误差的古地震发生时间范围.注意：在海原断裂的古地震序列中(图2a)，深红色部分表示两个古地震发生时间范围部分重叠了.Fig.2 Paleoseismic data in northeastern Tibetan Plateau (Data from Min et al.，2000)The red，green，blue and purple rectangles represent earthquakes on Haiyuan fault (a)，Xiangshan-Tianjingshan fault (b)，Luoshan fault (c)，and Helanshan fault (d).The length of every rectangle shows the occurrence time including errors.Please note that in the paleoseismic sequence on Haiyuan fault (Fig.2a)，the dark red part shows that the time spans of two events are partially overlapped.

2.2 人工合成地震目录数据

本文使用了我们前期数值模拟研究中的6个算例所产生的人工地震目录，并选取了稳定加载状态下(模型时间5万年之后)的目录数据进行分析(表1)(孙云强等，2019).

表1 有限元模型的6个算例(6个人工地震目录)的参数对比Table 1 Parameter differences in six cases (six synthetic seismic catalogs)of the finite-element model

算例1产生的人工地震目录显示：在四个主要断裂中，海原断裂和香山天景山断裂的地震活动性明显强于罗山断裂和贺兰山断裂(图3).其余5个算例也具有相同的地震活动特点.为了减少古地震序列与人工地震目录震级差异的影响，并更好地将两者进行匹配分析，本文综合考虑了断裂带的地震地质研究结果和数值模型结果，设置了震级阈值对人工地震目录中的地震进行过滤(Burchfiel et al.，1991；闵伟等，2000；Lin et al.，2015；孙云强等，2019).其中，海原断裂和香山天景山断裂的震级阈值选为7级，罗山断裂和贺兰山断裂的阈值为6级.

3 匹配对比方法

3.1 单断裂匹配

单断裂匹配是将不同断裂的人工地震目录数据分别取出，然后与同一断裂的古地震序列数据进行匹配.由于古地震序列和人工地震目录的地震数量不一致，不利于直接匹配对比这两种数据，因此本文需采用地震数目固定(与待匹配断裂的古地震数目一致)的滑动窗口方法对人工地震目录进行子序列划分，从而得到所有子序列.以海原断裂为例，海原断裂有10个古地震，因此本文先从相应的人工地震目录(图3b)的第1个地震开始，依次选取10个地震作为第1个子序列；然后再从该人工地震目录的第2个地震开始，同样选取10个地震作为第2个子序列；重复上述过程，就能得到所有子序列.

图3 算例1的人工合成地震目录(数据来自孙云强等，2019)(a)模拟的青藏高原东北缘断层系统上的地震活动;(b)—(e)分别为模拟的海原断裂、香山天景山断裂、罗山断裂、贺兰山断裂上的地震活动.横轴表示时间，纵轴是地震矩震级.不同颜色代表不同断裂上的地震.Fig.3 The synthetic seismic catalog from case 1 (Data from Sun et al.，2019)(a)The modeled seismicity on fault system in northeastern Tibetan Plateau；(b)—(e)show the modeled seismicity on Haiyuan fault，Xiangshan-Tianjingshan fault，Luoshan fault and Helanshan fault.Horizontal axis is time and vertical axis is seismic moment magnitude.Different colors denote earthquakes on different faults.

由于子序列的地震时间是模型时间，与古地震时间在数值大小上有较大的差异，因此本文假定两种序列的第1个地震时间相同，采用了公式(1)对子序列的地震时间进行平移处理：

ts′i,j=tsi,j-(tsi,1-to1),

(1)

其中，ts′i,j为平移后的第i个子序列中第j个地震的时间，tsi,j为平移前的第i个子序列中第j个地震的时间，to1为古地震序列中第1个地震的时间.平移后，本文使用了两种不同的匹配方法，计算子序列与对应古地震序列的匹配度.

3.1.1 平均绝对误差

平均绝对误差是一种反映两个序列实际误差大小的指标，在机器学习等领域中被广泛使用(Qi et al.，2020).该指标计算了两个长度相等序列的偏差的绝对值的平均，能避免计算过程中出现的误差相互抵消的问题(茆诗松，2004).平均绝对误差越接近于0，则说明两个序列越相近(贾俊平等，2009).计算公式如下：

(2)

dj=|ts′i,j-toj|,

(3)

其中，toj为古地震序列中第j个地震的时间，n为古地震数量，dj是两种数据的地震时间差的绝对值，Di为平均绝对误差结果.

本文也计算了由古地震时间误差所导致的平均绝对误差变化范围，其上下界由以下公式得到：

(4)

(5)

(6)

(7)

3.1.2 余弦相似度

余弦相似度是一种常见的相似性度量指标，在文本匹配、序列模式挖掘和信号比较等方面有广泛的应用(张振亚等，2005；杨毅明，2012；廖清科，2015).该指标将两个序列视为向量空间中的两个向量，通过计算两者之间的夹角余弦值衡量两个序列的相似度.余弦相似度的取值范围为[-1,1]，越接近1，说明两个序列越相近(廖清科，2015).在使用余弦相似度法计算时，为了避免地震时间本身的大小对计算结果的影响，本文将两种数据的地震时间序列转换为时间间隔序列.古地震序列的时间间隔计算公式如下：

dtoi=toi+1-toi.

(8)

子序列的时间间隔计算方法与古地震序列不同，还包含了古地震时间的影响，计算公式如下：

(9)

将两种数据的时间间隔序列代入(10)式即可计算得到余弦相似度：

(10)

考虑到古地震数据时间存在误差，本文采用一种遍历的方式计算余弦相似度的上下界，从而得到其可能的范围.遍历方法如下：首先将一个断裂的每一个古地震时间范围(如[toi-σi,toi+σi])五等分，每个古地震得到的五个时间点作为该地震时间所有可能的取值.若一个古地震的时间在文献记录中没有误差或误差不清晰，则认为该地震时间唯一.之后根据每一个古地震所有可能的时间取值，替换古地震序列中的地震时间并按时间顺序重新排序，即可得到多个不同的古地震时间序列.基于上述结果，在每一次匹配时，将所有可能的古地震时间序列与对应断裂的子序列匹配计算，计算结果中的最大值和最小值就分别为此次计算的上下界.以海原断裂为例，该断裂有10个古地震，其中有2个地震的时间误差为0，因此存在58种可能的古地震时间序列.每次与子序列匹配计算时，都需要进行58次计算，计算结果中的最大值和最小值分别为此次计算的上下界.

3.2 多断裂匹配

以单断裂匹配结果为基础，本文也进行了多断裂匹配，包括双断裂综合匹配(海原断裂和香山天景山断裂)和四断裂综合匹配.这两种匹配策略都是以海原断裂的子序列为参照，根据其他断裂与海原断裂第一个古地震的时间差的范围进行匹配(未使用公式(1)对多断裂匹配中的子序列地震时间进行平移处理).

(11)

四断裂综合匹配的方式与双断裂综合匹配类似，区别是考虑的断层数目不同.四断裂综合匹配通过(12)式进行计算：

(12)

4 序列匹配结果及大震迁移概率计算

4.1 参考算例(算例1)的匹配结果

以算例1为例，本文分别对平均绝对误差法和余弦相似度法的匹配结果进行了分析与评价(图4和图5).对于平均绝对误差法的匹配结果，我们以其下界(公式(5)所得结果)小于或等于0.1 ka作为该方法匹配结果较好的标准.单断裂匹配结果表明，各断裂的人工地震目录都存在与对应的古地震序列匹配较好的时间点，例如海原断裂和香山天景山断裂在55 ka左右；罗山断裂在54 ka左右；贺兰山断裂在58 ka左右.但不同断裂匹配较好的时间点存在差异，即在一个断裂匹配较好的时间点，其他断裂的人工地震目录与古地震数据的匹配未必较好.因此，需要用条件更严格的多断裂匹配方法来约束匹配结果.对于双断裂综合匹配，海原断裂和香山天景山断裂存在共同匹配较好的时间点，例如在55 ka左右.对于四断裂综合匹配，图中只有少数匹配较好的时间点，例如在56 ka左右.

图4 算例1的平均绝对误差匹配结果(a)—(f)分别为算例1中，海原断裂、香山天景山断裂、罗山断裂、贺兰山断裂、双断裂和四断裂的平均绝对误差匹配结果.(a)—(d)的横坐标为子序列的第一个地震时间，(e)—(f)的横坐标为海原断裂子序列的第一个地震时间，(a)—(f)的纵坐标为平均绝对误差.图中每一个时间点对应一对上下顶点，其中上顶点由公式(4)计算得到，下顶点由公式(5)计算得到，中间的区域代表可能的平均绝对误差取值.Fig.4 The MAE (Mean Absolute Error)matching results from case 1(a)—(f)show MAE matching results of Haiyuan fault，Xiangshan-Tianjingshan fault，Luoshan fault，Helanshan fault，two faults and four faults in case 1.In (a)—(d)，the horizontal axis is the time of the first event in the subseries.While in (e)—(f)，the horizontal axis is the time of the first event in the subseries of the Haiyuan fault.The vertical axis is MAE.Each time point in the figure corresponds to a pair of upper and lower vertices.The upper vertice is calculated by equation (4).The lower vertice is calculated by equation (5).The region between them includes possible MAE values.

图5 算例1的余弦相似度匹配结果(a)—(f)分别为算例1中，海原断裂、香山天景山断裂、罗山断裂、贺兰山断裂、双断裂和四断裂的余弦相似度匹配结果.(a)—(f)的横坐标与图4相同，纵坐标为余弦相似度.图中每一个时间点对应两个顶点，分别为计算结果中的上下界，中间的区域则是可能的余弦相似度取值.Fig.5 The cosine similarity matching results from case 1(a)—(f)show cosine similarity matching results of Haiyuan fault，Xiangshan-Tianjingshan fault，Luoshan fault，Helanshan fault，two faults and four faults in case 1.In (a)—(f)，the horizontal axis is the same as that in Fig.4.The vertical axis is cosine similarity.Each time point in the figure corresponds to two vertices，which are the upper and lower bounds in the calculation.The region between them includes possible cosine similarity values.

对于余弦相似度法的匹配结果，我们以其上界大于或等于0.4作为该方法匹配结果较好的标准.算例1的单断裂匹配、双断裂匹配和四断裂综合匹配的结果均显示在53 ka和56 ka左右匹配较好.将两种方法匹配较好的时间点进行对比(图4和图5)，两者有部分相同之处.例如，它们的四断裂综合匹配均显示在56 ka左右有较好的匹配结果.

图4与图5表明平均绝对误差法和余弦相似度法的匹配结果存在部分差异.本文认为其原因是两种方法的匹配方式不同.平均绝对误差法是通过计算两个序列在数值上的偏差来进行匹配，而余弦相似度法是通过计算两个序列对应的向量在方向上的差异来进行匹配.对于与古地震序列相近的子序列，采用这两种方法匹配计算时，都能有较高的匹配度.所以本文认为平均绝对误差法和余弦相似度法能互相补充，两种方法共同匹配较好的时间点对应的子序列与古地震序列更相近.

在青藏高原东北缘地区，海原断裂和香山天景山断裂是两条最重要的断裂；它们的滑动速度最快、地震活动最强(Burchfiel et al.，1991；Lin et al.，2015).因此，为了突出这两条断裂带的重要性，本文在之后进一步的分析中，使用了两种匹配方法共有的双断裂(海原断裂和香山天景山断裂)综合匹配较好的时间点.

图6展示了其余5个算例的平均绝对误差法和余弦相似度法的双断裂综合匹配结果.表2展示了6个算例的双断裂综合匹配较好的时间点.结果表明并非每个算例都存在双断裂综合匹配较好的时间点，如算例3、5和6(表2)，因此这些算例的人工地震目录不能与古地震序列很好地匹配.算例1、2和4的人工地震目录与古地震序列更相近(表2)，所以本文仅对算例1、2和4做进一步的分析.

图6 算例2—6的双断裂匹配结果(a1)—(a5)分别为算例2—6的平均绝对误差法的双断裂匹配结果.(b1)—(b5)分别为算例2—6的余弦相似度法的双断裂匹配结果.Fig.6 The matching results of two faults from case 2 to case 6(a1)—(a5)are the MAE matching results of two faults from case 2 to case 6.(b1)—(b5)are the cosine similarity matching results of two faults from case 2 to case 6.

4.2 地震迁移概率计算

基于4.1节中的结果，本文计算了算例1、2和4(表1)中，海原断裂或香山天景山断裂发生大地震后，下一次大地震发生在区域四条主要断裂上的概率，即大地震从海原断裂或香山天景山断裂迁移到区域某断裂的概率(我们称为某断裂的地震迁移概率).本文的地震迁移概率定义与前人研究中的定义相同(孙云强等，2019)，均表示上一次大地震在某条断裂上发生后，下一次大地震在区域每条断裂发生的个数与下一次大地震在这四条断裂上发生的总数的比值.

由于三个算例都存在多个共同匹配较好的时间点，难以全部用于计算分析，故本文在每个算例中仅使用一个匹配最好的时间点(表2).为了与古地震序列的时间跨度保持一致，本文是以选用的时间点为起始，在对应的人工地震目录中向后截取9 ka作为之后计算分析的时间段.

表2 六个算例的匹配结果及用于地震迁移概率计算的时间段Table 2 Matching results of the six cases and time spans for the calculation of the probability of earthquake migration

4.2.1 大地震在海原断裂发生后迁移到四条主要断裂的概率

综上所述，古地震数据的计算结果显示，大地震迁移到海原断裂的概率最大，其次是香山天景山断裂；算例1、2和4的计算结果比较接近，均显示大地震迁移到海原断裂的概率最大，其次是香山天景山断裂(图7).

图7 海原断裂发生大地震后，大地震迁移到四条主要断裂的概率横坐标从左到右分别为海原断裂(HYF)、香山天景山断裂(XTF)、罗山断裂(LSF)、贺兰山断裂(HLSF).纵坐标为地震迁移概率.图中蓝色、红色、绿色以及橙色柱状图分别代表算例1、2和4以及古地震数据的结果.Fig.7 Probability of the next big earthquake on the four major faults after one big earthquake occurring on the Haiyuan faultHYF：Haiyuan fault，XTF：Xiangshan-Tianjingshan fault，LSF：Luoshan fault，and HLSF：Helanshan fault.Vertical axis is the probability of earthquake migration.The columns in blue，red，green，and orange show the results from case 1，case 2，case 4，and the paleoseismic data，respectively.

4.2.2 大地震在香山天景山断裂发生后迁移到四条主要断裂的概率

图8 香山天景山断裂发生大地震后，大地震迁移到四条主要断裂的概率横坐标和纵坐标与图7相同.图中蓝色、红色、绿色以及橙色柱状图分别代表算例1、2和4以及古地震数据的结果.Fig.8 Probability of the next big earthquake on the four major faults after one big earthquake occurring on the Xiangshan-Tianjingshan faultThe horizontal axis and the vertical axis are the same as those in Fig.7.The columns in blue，red，green，and orange show the results from case 1，case 2，case 4，and the paleoseismic data，respectively.

综上所述，古地震数据显示大地震迁移到罗山断裂和贺兰山断裂的概率最大；算例1、2和4的计算结果存在着较大差别(图8).算例1显示大地震迁移到贺兰山断裂的概率最大，其次是罗山断裂.算例2显示大地震迁移到海原断裂的概率最大，其次是香山天景山断裂和罗山断裂.算例4显示大地震迁移到罗山断裂的概率最大，其次是贺兰山断裂和香山天景山断裂(图8).

5 讨论

人工合成地震目录能在一定程度上弥补当前地震目录记录时间短、不完备的缺点，但其与观测得到的地震数据的匹配程度却很少在研究中被考虑.因此，本文开发并使用了平均绝对误差法和余弦相似度法，对青藏高原东北缘地区四条断裂的古地震数据(闵伟等，2000)及从该地区地震活动数值模拟研究得到的长期人工地震目录结果数据(孙云强和罗纲，2018；孙云强等，2019)，进行时间上的匹配对比，分析并给出了匹配度较高的人工地震目录.根据匹配较好的人工地震目录，本文还计算分析了大地震在海原断裂或香山天景山断裂发生后，迁移到区域四条主要断裂上的概率.

当上一次大地震发生在香山天景山断裂时，三个算例的地震迁移概率结果有比较大的差别，其中只有算例1和4的地震迁移概率结果与古地震数据的计算结果较为接近，但优于孙云强等(2019)的研究结果.因此，此情形下的地震迁移概率对岩石圈黏度值比较敏感.这与孙云强等(2019)的结论有所不同.从目前我们的数值模拟研究结果来看，流变结构与黏度值会影响地震活动，但具体细节还需进一步的探究.

本文使用的人工合成地震目录并不能完全与古地震数据相匹配.我们认为可能存在以下两个主要原因.其一是古地震数据在时间上存在误差.其二是数值模型与地震地质研究(闵伟等，2000)所给定的断裂带长度不同.例如，地震地质研究给定的中卫-同心断裂仅表示数值模型中香山天景山断裂的中东部分(俞岗等，2013).如果我们能够将数值模型断裂单元与古地震破裂位置相对应，那么我们就能够从空间上更好地评估人工地震目录与古地震序列的匹配程度.这些问题可以在未来研究中得到更深入的调查与探索.

6 结论

本文开发并使用了平均绝对误差法和余弦相似度法，对青藏高原东北缘数值模拟产生的人工地震目录进行评估，得到了与古地震数据匹配度较高的人工地震目录.基于此目录，本文也计算了海原断裂及香山天景山断裂发生大地震后，大地震迁移到区域四条主要断裂的概率.我们得出如下结论.

(1)平均绝对误差法和余弦相似度法，在人工地震目录与古地震数据的匹配评估上都具有可用性.它们还可以相互约束，共同搜索得到与古地震数据匹配度较高的人工地震目录.

(2)匹配度较高的人工地震目录计算得到的地震迁移概率优于未进行匹配的人工地震目录计算结果.基于经过匹配的人工地震目录计算得到的地震迁移概率显示：当大地震在海原断裂上发生后，下一次区域的大地震在海原断裂上发生的概率最大，约为47%，其次是香山天景山断裂，约为23%～27%，均与古地震数据计算结果接近.

(3)数值模拟得到的人工地震目录，经过古地震数据匹配筛选后，才能更好地用于地震迁移概率计算及地震危险性分析.

致谢感谢北京大学蔡永恩教授的建议与支持及两位审稿专家的建设性意见.