一種基于距離變換和分水嶺算法的地震空區自動識別方法

王萍 陳皓一 侯謹毅

摘? ?要：地震空區的識別、分析是目前中期地震預報的重要手段之一，而傳統的人工繪制預報方法難以取得理想的效果，計算機視覺方法提供了解決問題的新途徑. 鑒于此，提出一種基于圖像處理的地震空區自動識別方法，輸入歷史地震的文本信息，通過距離變換、閾值分割、分水嶺算法等計算機視覺方法進行處理，通過迭代比較和特征參數篩選有效地震空區，輸出地震空區的分布圖像以及相應的特征參數. 通過具體案例進行實驗，研究表明：此方法可以獲得內部連通、邊緣清晰的地震空區圖像;與專家標定相比較，此方法的召回率為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.