一种基于距离变换和分水岭算法的地震空区自动识别方法

王萍 陈皓一 侯谨毅

摘   要：地震空区的识别、分析是目前中期地震预报的重要手段之一，而传统的人工绘制预报方法難以取得理想的效果，计算机视觉方法提供了解决问题的新途径. 鉴于此，提出一种基于图像处理的地震空区自动识别方法，输入历史地震的文本信息，通过距离变换、阈值分割、分水岭算法等计算机视觉方法进行处理，通过迭代比较和特征参数筛选有效地震空区，输出地震空区的分布图像以及相应的特征参数. 通过具体案例进行实验，研究表明：此方法可以获得内部连通、边缘清晰的地震空区图像;与专家标定相比较，此方法的召回率为81.25%，准确率为92.86%. 本文方法为地震研究工作者进行地震预测业务及相关研究提供了有力的工具.

关键词：地震空区;计算机视觉;特征参数;识别;迭代比较

中图分类号：TP391.4                          文献标志码：A

Abstract：The identification and analysis of seismic gaps is one of the important means of medium-term earthquake prediction. However， it is difficult to achieve the desired effect by the traditional artificial method. Computer vision provides a new way to solve the problem. In this paper， an automatic identification method for seismic gaps based on image processing methods is proposed. The input is the text information of the historical earthquake. It is processed by computer vision methods such as distance transformation， threshold segmentation and watershed algorithm. The effective seismic gap is screened by iterative comparison and feature parameters. The output is the distribution image of the seismic gaps and its corresponding characteristic parameters. In addition， the algorithm of this paper is tested through a certain case. The test suggests that the algorithm of this paper can clearly obtain the seismic gaps with internal connectivity and clear external contour. Compared with the expert calibration， the recall rate of this algorithm is 81.25%， and the accuracy is 92.86%. This method provides a powerful tool for seismic researchers to conduct earthquake prediction business and related research.

Key words：seismic gap;computer vision;characteristic parameters;identification;iterative comparison

自从地震空区的概念被Fedotov[1]正式提出以来，在实际地震预报工作中得到了广泛应用. Mogi[2]将地震空区分为两类，第一类地震空区是指板块边界带、活动断裂带上的最大地震破裂带的空间内，较长时间没有发生大地震的区域，未来大震将发生在此区域内，以填充此空区;第二类地震空区为前兆空区，是指未来大震发生前，小震活动减少的区域，这个区域将以断层的潜在破裂段为中心发展成为震源区[3]. 此外，还有学者提出背景空区和孕震空区等前兆空区[4-5]及地震空白区[6]等概念. 地震空区与不久前发生过大震的邻近区域相比，具有更高的应变力，因此将成为未来大地震发生的地点[7].

同时，不同的地震空区特征参数对应不同的地震等级. 例如，5～6级地震的空区长轴一般为100～300 km，围空的孔径方位角不大于120°等. 因此对地震空区特征参数的提取与分析可以作为地震等级预报的重要参考[8-10].

由于地震空区对于中长期地震预测的重要性 [11-14]，学者们进行了深入的研究. 王炜等[15]使用空区内外频次比和应变释放等来帮助识别地震空区;韩渭宾等[16]提出了一种基于R-t图和D-t图圈定地震空区的方法;陆远忠等[17]研究了在强震发生前地震活动模式的三维数值模型;Pei等[18]通过 Pg成像发现在空区的上地壳存在低速带;Wang等[19]利用P波和S波的趋势变化证明了地震空区与低速带有关;吕坚等[20]通过系统性地对《中国震例》中震级大于5的246个震例进行梳理，研究了震级与空区的关系;Hardebeck等[21]提出了一种基于震源机制确定区域平均应力的方法;Yang等[22]将该方法应用于地震空区，通过波形拟合反演出汶川-芦山地区地震空区周围的471个围空震源机制解;梁春涛等[23]也从应力的角度分析了地震空区的形成; Yin等[24]对地震空区的复发周期在华北盆地进行了探究;方进等[25]基于三角位错机制自动剖分断层，提取断层参数，反演了2015年尼泊尔MW7.8地震的发震断层单元及滑动空区分布.

上述研究成果在地震空区的分析和我国的地震预测预报中发挥了重要的作用，但由于不同学者的研究重点和认定标准各不相同，且我国对地震空区的识别圈定大多采用人工或半人工方法，存在以下缺点：

1）地震点数量多，筛选过程繁琐，人工效率低.

2）人工绘制空区具有较大的主观性和随机性，

空区大小、形状缺乏一致性.

3）人工寻找空区容易造成空区漏报，影响结果准确性.

4）人工提取空区特征过程复杂，测量过程容易产生误差，影响震级预测结果.

针对现有技术的不足和地震空区的特点，本文从计算机视觉的角度出发，提出一种地震空区的自动检测方法，同时可以获得地震空区的特征参数. 本文使用天津地震局记录的8 849条地震记录进行测试，测试结果与人工标定结果进行对比，得到了专家的肯定. 本文方法全自动性，避免了人工操作的随机性和不确定性，提升了地震空区的绘制效率和准确度.

1   地震空区自动检测方法

提出的地震空区自动检测方法的流程图如图1所示. 该算法输入历史地震文本数据，先将其转换为地震点分布的二值图像;使用两次距离变换和阈值分割、分水岭算法从地震点分布图像中分割出地震空区潜在目标区域，在分割过程中使用等间隔的不同阈值进行处理，通过基于不同阈值产生的图像上的目标区域的位置关系设置迭代规则，迭代比较以消除无效目标区域，整合为一张潜在地震空区分布图;对于剩余的潜在空区，利用基于轮廓的检测算法计算其空间特征，进一步消除不符合标准的地震空区，得到真正的地震空区分布图，并获得其特征参数. 1.1   生成地震点分布图

算法输入关于历史地震记录的文本信息，每条记录包括地震点的地理坐标（经纬度值）和震级. 算法的第一步是得到真实历史地震点的分布图像，先选定研究的时间范围和区域范围，提取在规定范围内的历史地震记录，通过坐标转换公式将地震点地理坐标文本信息转换为地震点分布图上的坐标，即：

式中：Plo和Pla分别为文本记录的某一地震点的经度值与纬度值;Plo min和Pla max分别为所选定范围内的所有地震点地理信息中的经度最小值和纬度最大值; Red为最大纬度值所属纬度圈的半径长度;Rec为最大纬度值所属纬度圈到地球球心的球半径长度.

坐标转换完成后，地震点经纬度文本信息转换成为地震点分布图上的坐标，将地震点分布图中地震点坐标处的灰度值设置为1，其余区域设置为0，得到地震点分布的二值图像.

1.2   获取潜在地震空区分布图像

为了获得潜在地震空区分布图像，首先要提取含有地震空区的潜在目标区域. 对地震点分布的二值图像依次进行距离变换、阈值分割、距离变换、阈值分段、分水岭算法（简称DSDSW过程）的处理.

进行第一步距离变换时，假设地震点分布二值图上白色地震点的集合为B，其他黑色背景点的集合为F. 对于F中的每个点（x，y），它与B的最小距离为：

通过式（3）进行距离变换后，地震点分布二值图像被转换为距离图像（图2（a）），图2（a）中越亮的地方表示距地震点越远，反之越近.

对距离图像用阈值T进行分割，得到表示地震空区位置的局部最大值区域（图2（b）），图2（b）中的白色区域为地震空区潜在目标区域所处的基本位置. 但此刻局部最大值区域的形状、大小均不规则，不能表示完整的地震空区. 为了解决这一问题，第二次使用距离变换公式（3）对图2（b）进行处理，其中白色点为前景点B，黑色点为背景点F;之后使用阈值T对其进行分段，将距离变换后距离超过T的点全部置为黑色，结果如图2（c）所示. 图2（c）中的每个类圆形封闭区域代表完整的地震空区潜在目标区域.

通过分水岭算法[26-27]进行图像分割，得到边界清晰的地震空区潜在目标区域，获得地震空区的潜在目标区域分布图.

为保证潜在地震空区检测的全面性和准确性，在上述DSDSW过程中以相等间隔增加分割过程中的阈值，获得经不同阈值处理的N （N≥3） 张地震空区潜在目标区域分布图. 假设经最小阈值处理得到的潜在目标区域分布图像为P1，经最大阈值处理得到的潜在目标区域图像为PN，则所获得的潜在目标区域分布图像的集合可以表示为{Pi}，i = 1，…，N.

对{Pi}（i < N）进行迭代处理，以消除无效的地震空区潜在目标区域. 迭代规则基于相邻阈值处理得到的图像Pi与Pi + 1上地震空区潜在目标区域的位置关系而设立. 假设Ra与Rb分别为图像Pi与Pi + 1（i

1）Ra和Rb所处位置重叠，且重叠区域的面积大于Ra面积的70% ;

2）Ra与Pi + 1中多个潜在目标区域重叠，且重叠区域的总面积大于Ra面积的70% .

根据上述规则，按所选阈值的递增顺序对{Pi }，i =1，…，N中的所有图像进行迭代比較. 在每一步迭代完成后，将比较过的两张图像上保留的潜在目标区域合并成一张图像，所有迭代过程完成后获得潜在地震空区的分布图像.

1.3   地震空区特征参数的提取与地震空区的识别

获得潜在地震空区分布图像后，还需要筛选地震空区的特征参数，以获得真正的地震空区分布图. 本文选取的地震空区特征参数为地震空区长轴长度L和孔径方位角度数α[28].

L和α的计算方法如图3所示，图3中黑色轮廓线为圈定的地震空区，实心小黑点为围空地震点，具体计算方法如下：

L的计算方法. 通过轮廓检测算法对潜在地震空区分布图像中的潜在空区求最小外包矩形（图3中的外侧矩形框），读取潜在空区对应的最小外包矩形的边长，选择较大值作为潜在地震空区长轴长度值L.

α的计算方法. 通过轮廓检测算法对潜在地震空区分布图像中的潜在空区求最大内接矩形（图3中的内侧矩形框），将内接矩形中心点作为空区中心点（图3中的实心小灰点）. 沿潜在空区轮廓线进行顺序遍历，获得围空地震点的坐标，计算每两个相邻的围空地震点与中心点形成的孔径方位角度数α.

在所有潜在地震空区的特征参数计算完成后，进一步对地震空区进行筛选. 一个有效的地震空区应满足：（1） L > 100;（2） 对任一α，α<120°. 在潜在地震空区分布图上，同时满足这两个条件的潜在地震空区判定为有效的地震空区. 此时，得到最终的地震空区分布图像.

2   案例分析

使用案例来演示本文算法的效果，同时也详细展示了算法每个步骤的执行结果，最后与专家标定结果进行比较.

本案例的研究数据由天津地震局提供，所研究的时间范围在2008年1月1日零点至2017年12月31日24点内、地理范围在中国大陆108°～125° E与30°～43° N内，包含每次地震的地理位置和震级，总计8 849条.

根据地震空区自动检测方法，首先对地震文本数据进行坐标变换. 将地震点的经纬度坐标通过坐标变换公式（3）转换为平面坐标X和Y，获得地震点的分布图像，其中最小经度值为108° E，最大纬度值为43° N，因此地震点分布图的坐标原点为（108，43）.

对地震点分布图像进行DSDSW处理，在进行阈值分割时，按照阈值初值为25，公差为5，最大值为80的等差数列，选择12个阈值，获得对应12张地震空区潜在目标区域分布图像（图4）. 图4中每个分图下方的数字表示所选阈值.

迭代比较图4中具有相邻阈值的图像以去除无效的地震空区潜在目标区域.

使用实例进一步说明该比较过程（图5）. 图5（a）和5（b）分别为阈值为65和70时所获得的地震空区潜在目标区域分布图像，其中潜在目标区域个数分别为7和6. 表1列出了图5（a）和5（b）中每个潜在目标区域的面积.

将图5（a）和5（b）进行比较，根据1.2节中的潜在目标区域比较规则，图5（a）中的无效潜在地震空区目标区域将被移除，表2给出了比较过程中每对重叠潜在目标区域的参数. 如根据潜在目标区域比较规则1，图5（a）中的潜在目标区域1、3、4、5、6、7分别与图5（b）中的潜在目标区域1、2、3、4、5、6重叠，且重叠区域面积大于图5（a）中潜在目标区域面积的70%，因此移除图5（a）中的潜在目标区域1、3、4、5、6、7;而图5（a）中的潜在目标区域2与图5（b）的背景相交，因此图5（a）中的潜在目标区域2被保留，将图5（a）中保留的潜在目标区域2与图5（b）进行合并，得到比较结果（图5（c））.

通过上述方法将图4中的12张地震空区潜在目标区域分布图像按照阈值由小到大的顺序进行迭代比较，得到潜在地震空区分布图像（图6），其中潜在地震空区为109个.

此时图6还不是最终结果，还将通过特征参数进一步过滤这些地震空区. 利用1.3节中的方法，计算图6中的109个潜在地震空区的特征参数. 其中14个地震空区符合1.3节中的筛选标准，表3列出了这14个地震空区的特征参数计算结果，其余95个无效的潜在地震空区被移除，得到最终的地震空区分布图如图7所示. 图7中分布在地震空区周围的实心小黑点为形成地震空区的地震点，文字为重要城市名称的汉语拼音，汉语拼音左下角白点为城市所在位置.

本文结果（图7）中的1～13号地震空区与专家标定结果（图8）中的1～13号吻合，图7中14号地震空区未出现在专家标记结果中，召回率为81.25%，准确率为92.86%. 漏报出现在图8的14号、15号和16号位置，其中16号地震空区是图7的11号空区的组成部分，11号空区呈“┐”形，在专家标定时被分为图8的11号和16号两部分;误报出现在图7的14号地震空区位置，该位置并未被专家标定为空区，14号地震空区作为本文识别结果中最小的地震空区，其形状并不规则，但符合特征参数的筛选条件.

综合漏报和误报的分析结果，地震空区一般呈圆形或椭圆形，而图7中11号、13号和14号空区为不规则形状，可能空区还未完全形成或分别由两个空区组成，今后还会由于新的地震而分离，在未来的研究中可加入有关形状的筛选规则.

本文结果经天津市地震局专家认证，所得地震空区分布结果有极大参考价值. 其中11号地震空区在2018年2月9日，河南南阳市淅川县（32.83°N，111.56°E）发生了4.3级地震，但发生区域靠近空区南部边缘，并非在空区中央，将有可能进一步形成地震空区.

3   结   论

本文提出了一种基于计算机视觉的方法解决地震空区自动识别的问题. 输入地震点地理文本信息，通过距离变换、阈值分割、分水岭算法、迭代比较、轮廓检测等图像处理方法，输出地震空区的分布图和相应的特征参数. 本文实验结果与专家标定相比较，召回率为81.25%，准确率为92.86%. 实验结果经专家认证准确有效.

与现有人工或半人工繪制技术相比，地震空区的自动识别避免了人工绘制的随机性和主观性，保证了空区绘制的全面性、规范性和一致性，提高了检测效率与预测准确性，为地震研究工作者进行地震预测业务及相关研究提供了有力的工具. 同时，提取的地震空区特征参数可以用于地震震级的判断，便于对地震预测进行科学、细化地分析.

参考文献

[1]    FEDOTOV S A. Regularities of the distribution of strong earthquakes in Kamchatka，the Kuril Island and north-eastern Japan [J]. Trudy Insitute Physical Earth Academic Science，1965（36）：66—93.

[2]    MOGI K. Two kinds of seismic gaps [J]. Pure and Applied Geophysics，1979，117（6）：1172—1186.

[3]    陸远忠，宋俊高. 地震空区与地震预报[M]. 北京：地震出版社. 1989：1—55.LU Y Z，SONG J G. Seismic gaps and earthquake prediction [M]. Beijing：Seismological Press，1989：1—55. （In Chinese）

[4]    陆远忠，沈建文，宋俊高. 地震空区与“逼近地震”[J]. 地震学报，1982（4）：3—12.LU Y Z，SHEN J W，SONG J G. Seismic gaps and “imminent earthquakes” [J]. Acta Seismologica Sinica，1982（4）：3—12. （In Chinese）

[5]    陆远忠，王炜，沈建文. 地震空区与“逼近地震”续[J]. 地震学报，1983，5（2）：129—144.LU Y Z，WANG W，SHEN J W. Seismic gaps and imminent earthquakes （contitued） [J]. Acta Seismologica Sinica，1983，5（2）：129—144. （In Chinese）

[6]    梅世蓉. 地震科学研究论文选集[M]. 北京：地震出版社. 1993：15—30.MEI S R. Selection of earthquake research paper [M]. Beijing：Seismological Press，1993：15—30. （In Chinese）

[7]    SYKES L R. Aftershock zones of great earthquakes，seis-micity gaps，and earthquake prediction for Alaska and the Aleutians[J]. Journal of Geophysical Research，1971，76（32）：8021—8041.

[8]    张敏. 2000年6月6日甘肃景泰5.9级地震前中短期震兆特征研究[J]. 高原地震，2002（4）：25—30.ZHANG M. Study on characteristics of the mid-short term precursors before Jingtai earthquake with Ms5.9 on Jun.6 2000 in Gansu province [J].  Plateau Earthquake Research，2002（4）：25—30. （In Chinese）

[9]    张国民， 傅征祥， 桂燮泰. 地震预报引论[M]. 北京：
科学出版社，2001：
12—64.ZHANG G M，FU Z X，GUI X T. Introduction to earthquake prediction[M]. Beijing：
Science Press， 2001：12—64. （In Chinese）

[10]  曲延军，王海涛，邬成栋，等. 中国大陆地震空区统计特征分析[J]. 地震学报，2010，32（5）：544—556.QU Y J，WANG H T，WU C D，et al. Analysis on statistic characteristic of seismic gap in Chinese mainland [J]. Acta Seismologica Sinica，2010，32（5）：544—556. （In Chinese）

[11]  GENTILI S，GIOVAMBATTISTA R D，PERESAN A. Seismic quiescence preceding the 2016 central Italy earthquakes [J]. Physics of the Earth & Planetary Interiors，2017，272：27—33.

[12]  WYSS M，WIEMER S. Two current seismic quiescences within 40 km of Tokyo [J]. Geophysical Journal of the Royal As-tronomical Society，2018，128（2）：459—473.

[13]  ROMAN D C，RODGERS M，GEIRSSON H，et al. Assessing the likelihood and magnitude of volcanic explosions based on seismic quiescence [J]. Earth & Planetary Science Letters，2016，450：20—28.

[14]  TRAITANGWONG P，PAILOPLEE S. Precursory seismic quiescence along the sagaing fault zone，Central Myanmar - application of the region-time-length algorithm [J]. Geosciences Journal，2017，21（4）：1—10.

[15]  王炜，宋俊高，郑兆苾，等. 安徽省及邻区地震预报方案的研究[J]. 地震，1989（4）：1-8.WANG W，SONG J G，ZHENG Z B，et al. Research on earthquake prediction program in and around Anhui province [J]. Earthquake，1989（4）：1—8. （In Chinese）

[16]  韩渭宾，江道崇，杨虹，等. 借助 R-t 图、D-t 图、结合震中分布图在四川地区圈定空区的研究报告[J]. 四川地震，1989（4）：85—97.HAN W B，JIANG D C，YANG H，et al. Research report on identification of seismic gaps in Sichuan province based on R-t map，D-t map and epicenter distribution map [J]. Earthquake Research in Sichuan，1989（4）：85—97. （In Chinese）

[17]  陆远忠，叶金铎，蒋淳，等. 中国强震前兆地震活动图像机理的三维数值模拟研究[J]. 地球物理学报，2007，50（2）：177—186.LU Y Z，YE J D，JIANG C，et al. 3D numerical simulation on the mechanism of precursory seismicity pattern before strong earthquake in China [J]. Chinese Journal of Geophysics，2007，50（2）：177—186. （In Chinese）

[18]  PEI S，ZHANG H，SU J，et al. Ductile gap between the Wenchuan and lushan earthquakes revealed from the two-dimensional pg Seismic tomography [J]. Scientific Reports，2014，4：6489.

[19]  WANG X S，LU J，XIE Z J，et al. Focal mechanisms and tectonic stress field in the north-south seismic belt of China [J]. Chinese Journal of Geophysics，2015，58（11）：4149—4162.

[20]  呂坚，宋美琴，周龙泉，等. 基于《中国震例》的地震空区和地震条带统计特征[J]. 地震，2016，36（4）：22—34.L? J，SONG M Q，ZHOU L Q，et al. Statistical characteristics of seismic gaps and seismic belts based on earthquake cases in China[J]. Earthquake，2016，36（4）：22—34. （In Chinese）

[21]  HARDEBECK J L，MICHAEL A J. Stress orientations at intermediate angles to the san andreas fault，California [J]. Journal of Geophysical Research Solid Earth，2004，109（B11）：B11303.

[22]  YANG Y，LIANG C，LI Z，et al. Stress distribution near the seismic gap between Wenchuan and Lushan earthquakes [J]. Pure and Applied Geophysics，2017，174（6）：2257—2267.

[23]  梁春涛，黄焱羚，王朝亮，等. 汶川和芦山地震之间地震空区综合研究进展[J]. 地球物理学报，2018，61（5）：350—364.LIANG C T，HUANG Y L，WANG C L，et al. Progress in the studies of the seismic gap between the 2008 Wenchuan and 2013 Lushan earthquakes [J]. Chinese Journal of Geophysics，2018，61（5）：350—364. （In Chinese）

[24]  YIN A，YU X J，SHEN Z K，等. 华北盆地可能存在的一个地震空区及其高地震危险性[J]. 世界地震译丛，2018（1）：28—34.YIN A，YU X J，SHEN Z K，et al. A possible seismic gap and high earthquake hazard in the north China basin[J]. Translated World Seismology，2018（1）：
28—34. （In Chinese）

[25]  方进，许才军，温扬茂. 基于方差分量估计的2015年尼泊尔MW7.8地震同震滑动分布[J]. 地球物理学报，2019，62（3）：113—129.FANG J，XU C J，WEN Y M. Coseismic slip distribution of 2015 Gorkha （Nepal） Mw7.8 earthquake determined using the Helmert variance component estimation[J]. Chinese Journal of Geophysics，2019，62（3）：113—129. （In Chinese）

[26]  徐天芝，张贵仓，贾园. 基于形态学梯度的分水岭彩色图像分割[J]. 计算机工程与应用，2016，52（11）：200—203.XU T Z，ZHANG G C，JIA Y  Color image segmentation based on morphology gradients and watershed algorithm [J]. Computer Engineering & Applications，2016，52（11）：200—203. （In Chinese）

[27]  杨高波，张兆扬. 基于小波域分水岭变换的视频对象提取算法[J]. 湖南大学学报（自然科学版），2005，32（5）：108—113.YANG G B，ZHANG Z Y. An automatic video object extraction algorithm based on watershed in the wavelet domain [J]. Journal of Hunan University （Natural Sciences），2005，32（5）：108—113. （In Chinese）

[28]  JIAN L，MEIQIN S，LONGQUAN Z，et al. Statistical characteris-tics of seismic gaps and seismic belts based on earthquake cases in China [J]. Earthquake，2016，36（4）：22—34.