實時地震監測臺網中Voronoi圖生成算法的改進*

趙 瑞，劉希強

(山東省地震局，山東濟南250014)

0 引言

建設地震監測臺網的主要科學目的之一就是一旦地震發生后，第一時間開展地震速報工作。當前基于實時地震監測臺網建立地震預警系統和地震應急控制系統在國際上已形成潮流。日本、美國、墨西哥、土耳其、意大利、瑞典、羅馬尼亞、中國臺灣等地震多發國家和地區發展了多套預警系統，如日本鐵路UrEDAS系統、墨西哥SAS和SASO系統、美國ElarmS和PreSEIS系統、土耳其PreSEIS和SOSEWIN系統等，并有一些成功預警且收到減災實效的先例(Gasparini et al，2007;Wu et al，2007;Allen et al，2009;殷海濤等，2012)。在地震預警技術中，Voronoi圖(以下簡稱V圖)在地震事件觸發判斷、地震定位等方面已得到初步應用(Rydelek，Pujol，2004;Horiuchi et al，2005;Satriano et al，2008;王慶民等，2012);同時在地理信息系統、計算機圖形學、氣象學、虛擬實現、機械工程、機器人等領域也有著廣泛的應用(劉金義，劉爽，2004;范會川等，2010)。

V圖最早可追溯至1908年，Georgy Voronoi首先在數學上限定了每個離散點的有效作用范圍，提出了二維 V圖(Voronoi，1908)。1934年，Boris Delaunay提出Delaunay三角網，V圖與Delaunay三角網格互為對偶圖，可通過Delaunay三角網格尋找與每個離散點相關的三角形，順序連接每個三角形的外接圓心得到多邊形，即為V圖(孫繼忠等，2010)。經過眾多專家學者的研究，Delaunay三角網的生成算法日趨成熟。目前主要有分割—歸并法、逐點插入法和三角網格生長法及每種方法的改進算法，其中前兩種方法應用較多(Aurenhammer，1991;Delaunay，1934;Okabe et al，1994)。每種算法各有其優缺點，分割—歸并算法適合離散點較多的情況，但該算法含有大量迭代運算，空間復雜度高，并且凸殼形成時，算法較復雜，時間復雜度也較高(周培德，2000)。逐點插入法算法簡單，占用內存少，空間和時間復雜度不高，但沒有很好地解決凸殼問題。三角網格生長法算法簡單，以基邊為循環條件，產生凸殼時較困難，同時易產生邊與邊相交，不利于Delaunay三角網格優化。武曉波等(1999，2000)將逐點插入法與分割—歸并法相結合，提高了計算速度，降低了算法復雜度，但同樣沒有很好解決凸殼、初始三角網的生成問題。王強等(2010)在此合成算法的基礎上，提出一種改進的Delaunay三角網生成算法，采用虛擬四邊形凸殼，較好解決了初始三角網的生成問題，但是在刪除虛擬頂點時，若想保證實際凸殼的正確性，必須提前在程序中對離散的數據點進行預處理排序，排序原理沒有給出很好的說明，同時采用這種合并算法對Delaunay三角網優化時，復雜度高于其它算法。

在實時地震監測臺網中，當在臺站中斷記錄時段發生地震，利用V圖進行地震定位必須要調整V圖布局，若考慮所有正常臺站重新計算V圖則復雜度會明顯加大，如何能在不改變全局V圖的前提下，快速完成對中斷臺站V圖單元的重新剖分，又不失精度，是地震預警系統中需要解決的問題。

1 V圖生成算法

以山東測震臺網79個臺站為例，生成V圖的主要技術思路是選擇算法較簡單的逐點插入法，采用向量關系解決節點定位(王強等，2010)，同時對點集進行預處理、Delaunay三角網生成進行改進，從而保證凸殼的客觀性。

1.1 初始凸殼的生成

采用虛擬矩形方式生成初始凸殼，由4個頂點(［Min(x)、Min(y)］、 ［Min(x)、Max(y)］、［Max(x)、Min(y)］、［Max(x)、Max(y)］)形成矩形，同時判定頂點是否屬于原有節點，如果不屬于，則將該虛擬頂點加入到離散點集進行循環計算。分兩種情形生成初始三角形:一是當4個頂點不相重復，則按圖1a生成初始三角形;二是當4個點有相重點時，則按圖1b生成初始三角形。根據節點所在位置，在圖1上顯示R1、R2、R3和R4四個區域。

圖1 初始三角形生成圖情形1(a)和情形2(b)Fig.1 Generating the initialization triangle in situation 1(a)and situation 2(b)

1.2 點集預處理

為保證凸殼的有效生成，即刪除原始虛擬凸殼頂點時不出現凹陷面，需要在逐點插入算法前，對節點進行排序。圖2a中當前節點M位于第R1區域時，連接MA，MB，MC生成3個新的三角形區域①、②、③。圖2b中當后節點N位于①區時，做線段CM的延長線至AB，交點為C1，當點N位于△C1BM內部時(圖2b中陰影區)，能夠滿足BNMC的連線為凸面，則將N點排在M點后，反之，調換M與N點順序。圖2c中當點N位于②區時，做線段BM的延長線至AC，交點為B1，當點N位于△C1BM內部時，點BMNC的連線為凸面，則將點N排在M點后，反之，調換M與N點的順序。當點N位于圖2d陰影區內時，點BMC的連線均為凸面，則點N排在點M后。總之，當點N位于圖1第R1區中如圖2e所示陰影區內時，不用調換M與N點的順序，反之需要調換順序。同理，當點N位于圖1第R2、R3和R4區域內的陰影區時，不需要調換M與N點的順序(圖2f～h)。

圖2 點集預處理(a)M點位于R1區;(b)N點位于R1區中①區;(c)N點位于R1區中②區;(d)N點位于R1區中③區;(e)M點位于R1區;(f)M點位于R2區;(g)M點位于R3區;(h)M點位于R4區Fig.2 Preprocessing of point set(a)point M locates in zone R1;(b)point N locates in zone①of zone R1;(c)point N locates in zone②of zone R1;(d)point N locates in zone③of zone R1;(e)point M locates in zone R1;(f)point M locates in zone R2;(g)point M locates in zone R3;(h)point M locates in zone R4

1.3 逐點插入算法

根據王強等(2010)向量叉乘算法來判定點與三角形的位置。逐點插入時，每當新邊、新三角形生成時，其編號權重設置為1。

1.4 Delaunay三角剖分優化算法

圖3 Delaunay三角優化前(a)和優化后(b)Fig.3 Delaunay triangles before(a)and(b)after optimization

判斷新生成的三角形是否為Delaunay三角形，即判斷是否滿足Delaunay三角形的空圓特性。首先判斷四邊形ABCD是否為凸四邊形。為避免復雜的平方、開方計算，同樣引入向量叉乘的概念。如圖3所示有兩個待優化的三角形ABD和BCD，連接AD，DC，CB，AB，得到 4個向量，如式(1)所示:

當四邊形為凸四邊形時，向量DA、CB在DC的同側，并且向量AD、BC在AB的同側，即T1與T2同號，且T3與T4同號(如式(2)所示)。確定為凸多邊形后，分別求取△ABD、△BDC的外接圓圓心，判定是否滿足空圓性。如果圓內有其他頂點，則交換四邊形的對角線。原△ABD及△BDC的權重設置為0，邊BD的權重設置為0，新生成的兩個三角形，即△ADC及△ABC的權重設置為1，邊AC的權重設置為1。

為降低Delaunay三角網格優化的復雜度，在逐點插入過程中進行優化。當所有離散點插入之后，再次對全體三角形進行優化，以期全部達到Delaunay三角剖分的準則。按照本文選取的研究區域(33.0～40.0°N，114.5～122.5°E)，將經緯度轉換成直角坐標，研究范圍為800 km×800 km。由山東地區測震臺網生成的Delaunay三角剖分圖見圖4。

圖4 由山東測震臺網生成的Delaunay三角剖分圖Fig.4 Delaunay triangulation diagram generated by Shandong Seismic Network

1.5 V圖生成算法

存儲每個臺站所涉及三角形的個數及編號，按一定順序對三角形外接圓心進行排序連接，即得到此臺站的V圖單元。同理，依次遍歷所有臺站，得到研究區域的V圖見圖5。

圖5 由山東測震臺網生成的Voronoi圖Fig.5 Voronoi diagram generated by Shandong Seismic Network

2 臺站中斷時V圖生成算法

在實際地震監測臺網運行中，地震臺站通常會因為中斷或其他原因未記錄到波形信號，如果地震恰發生在中斷臺站所在的V圖單元內，則根據首先觸發臺站判定地震發生的區域時會產生較大的震中判定誤差，從而影響地震定位的精度。所以在地震發生時需要實時判定正常運行臺站的數量，快速準確地對中斷臺站的V圖單元進行調整和歸并。

本文給出3種臺站中斷例子，以說明V圖生成算法及其結果有效性。假設有5個臺站中斷，分別是55號曲阜臺(代碼QUF，下同)、64號泰安臺(TIA)、35號濟南臺(JIN)、56號榮成臺(RCH)、31號河口臺(HEK)，5個中斷臺站中HEK臺屬于孤立點，且屬于網內的臺站中斷情形;JIN、TIA和QUF臺屬于片區中斷臺;RCH臺屬位于凸殼邊界上臺站中斷情形(圖4)。

首先判斷中斷臺站是否為孤立中斷點，從而保證因臺站中斷所形成的中斷區域邊界頂點均為運行良好的臺站。將此區域內的三角形及邊的權重設置為0，形成一個空區(圖6a)，記錄此空區邊界上的點號及邊號。當中斷臺站處于網內時，空區為封閉區域(圖6b、c)，當中斷臺站處于網緣時，只有邊界線，為半封閉區域(圖6d、e)。當多個中斷臺站相鄰時，此空區不一定為凸多邊形，有可能存在凹陷的邊界(圖6c)，這時無論采取分割—歸并法、逐點插入法或傳統的三角網格生長法均不能完全保證空區邊界點所重新形成的凸殼與空區邊界線重合，并且易存在空區新凸殼與空區邊界線相交的情況，這為小空區內Delaunay三角剖分、V圖生成增加了難度。筆者在傳統的三角網格生長法上提出了一種改進方法，使其保證空區凸殼的完整，步驟如下:

圖6 中斷臺站形成的空區圖以及四種不同情況下三角形形成示意圖(a)空區圖;(b)中斷臺站處于網內時，形成封閉區域情形1;(c)中斷臺站處于網內時，形成封閉區域情形2;(d)中斷臺站處于網緣時，形成半封閉區域情形1;(e)中斷臺站處于網緣時，形成半封閉區域情形2Fig.6 Empty areas diagram formed by the interrupt stations and the schematic of triangular generated in four different cases(a)empty area;(b)forming a enclosed area under the case 1 when the interrupt stations in the network;(c)forming a enclosed area under the case 2 when the interrupt stations in the network;(d)forming a semi-closed area under the case 1 when the interrupt stations in the network edge;(e)forming a semi-closed area under the case 2 when the interrupt stations in the network edge

(1)如圖6b～e所示，中斷點設為點O，選取空區中任意兩條相交邊AB、BC，當AO、AC在AB邊同側且CA、CO在CB邊同側，即T5與T6同號，T7與 T8同號時，連接 AC形成新三角形△AOC、△ABC(見圖6b、6d)，反之不能形成新的三角形，如圖6c、6e所示。

(2)當空區內所有三角形形成后，在空區內對其進行優化，使之成為Delaunay三角形。針對圖6a中存在5個中斷臺站的情形，應用上述思路，生成的新的Delaunay三角剖分圖見圖7a，3個空區新生成的三角形用不同顏色標識出。

(3)遍歷所有臺站，生成新V圖(圖7b)，其中由于臺站中斷而涉及V圖變化的區域用不同顏色給出。

圖7 考慮地震臺站中斷時生成的Delaunay三角剖分圖(a)和新生成的Voronoi圖(b)Fig.7 Delaunay triangulation(a)and the new Voronoi diagram(b)generated by considering the interrupted seismic stations

3 計算時間測試

為測試算法改進后的實際效果，我們進行了計算時間的測試，這也是在地震預警系統中必須解決好的一個科學問題。本文程序編寫采用Fortran語言，時間對比均采取相同的硬件條件(聯想臺式機，主頻2.83G)。下面分為兩種情況進行討論:

第一種是臺站不存在中斷的情況，采用改進后的逐點插入法與傳統的逐點插入法進行對比。分別模擬1 000個、2 000個、5 000個、10 000個節點數據，運算得到Delaunay三角網格的CPU計算時間，對比如表1。

表1 Delaunay三角網格剖分算法改進前后對比表Tab.1 The comparison of delaunay triangulation algorithm before and after improvement

第二種是模擬山東地區存在5個中斷臺站，除去中斷點，其余74個臺站點重新剖分生成Delaunay三角網格與本文提出的局部采用改進后的三角網格生長法所用的時間進行對比分析。前者生成的Delaunay三角網格如圖8所示。這兩種方法所用的時間分別是0.35和0.10 s，生成V圖的計算時間分別是0.52 s和0.15 s。

圖8 除去5個中斷臺站后新生成的Delaunay三角剖分圖Fig.8 The new generation Delaunay triangulation by removing five interrupt stations

通過兩種計算時間對比，可以看出本文提出的V圖生成算法優于傳統的逐點插入法，在實際臺網中，臺站數量較少的情況下，基本上是瞬間完成的，并且地震臺站中斷時，局部采用改進后的三角網格生長法與除去中斷點、重新剖分Delaunay三角網格得到的結果是一致的，但在計算用時上優于后者。

4 結語

本文對實時地震監測臺網V圖的逐點插入算法的細節進行了梳理，明確了點集預處理的準則，改進了Delaunay三角網格優化算法，提高了計算速度，并且考慮了實際臺網運行過程中，地震臺站中斷時，V圖重新生成的情況。其突出的優點是在不改變全局V圖的前提下，局部區域采用三角網格生長法，完成對中斷臺站V圖單元的重新剖分，即使在網緣臺站中斷的情況下，也能有效地約束凸殼的生成。通過計算時間的對比，可以看出在同等效果的基礎上，計算用時要優于除去中斷臺站，重新剖分Delaunay三角網格。

V圖在地震監測預報中具有潛在應用前景，如根據觸發臺站所在V圖單元的集中性判定是否有地震事件發生，根據V圖單元面積，搜索得到區域內相關地理信息及應急資源信息，同時結合其它臺站記錄信息，迅速判定震中位置，為地震早期預警系統提供有效的技術支持。

本文在撰寫過程中得到國家科技支撐課題組(2012BAK19B04)和地震科技星火計劃項目組(XH12029)成員的幫助，在此表示衷心感謝。

范會川，周慧娟，賈利民.2010.Voronoi圖和Delaunay三角網在鐵路應急管理中的運用［J］.中國科技信息，(8):245－247.

劉金義，劉爽.2004.Voronoi圖應用綜述［J］.工程圖學學報，25(2):125－132.

孫繼忠，胡艷，馬永強.2010.基于Delaunay三角剖分生成Voronoi圖算法［J］.計算機應用，30(1):75－77.

王強，鄭逢斌，喬保軍，等.2010.一種改進的Delaunay三角網生成算法［J］.計算機應用與軟件，27(8):138－140.

王慶民，劉希強，沈得秀.2012.Voronoi圖和雙曲線聯合方法在地震快速定位中應用［J］.西北地震學報，34(3):234－244.

武曉波，王世新，肖春生.1999.Delaunay三角網的生成算法研究［J］.測繪學報，28(1):28 －35.

武曉波，王世新，肖春生.2000.一種生成Delaunay三角網的合成算法［J］.遙感學報，4(1):32－35.

殷海濤，劉希強，李杰，等.2012.現今地震預警技術及其在國內發展狀況的探討［J］.中國地震，28(1):1－9.

周培德.2000.計算幾何算法分析與設計［M］.北京:清華大學出版社.

Allen R.M.，Gasparini P.，Kamigaichi O.，et al.2009.The Status of Earthquake Early Warning aroud the World:An Introductory Overview［J］.Seismol.Res.Lett.，80(5):682 －693.

Aurenhammer F..1991.Voronoi diagrams-a survey of a fundamental geometric data structure［J］.ACM Computing Surveys，23(3):345－405.

Delaunay B Sur la Sphere Vide.1934.Bulletin of the Academy of Sciences of the USSR［J］.Classe des Sciences Mathematiques et Naturelles，8:793 －800.

Gasparini P，Manfredi G，Zschau J.2007.Earthquake Early Warning System［M］.Germany:Springer-Verlag，249 －282.

Horiuchi S.，Negishi H.，Abe K.，et al.2005.An Automatic Processing System for Broadcasting Earthquake Alarms［J］.Bulletin of the Seismological Society of America，95(2):708 －718.

Okabe A.，Boots B.，Sugihara K..1994.Nearest neighborhood operations with generalized Voronoi diagrams:A review［J］.International Journal of Geographical Information Systems，8(1):43 －71.

Rydelek P.，Pujol J..2004.Real-Time Seismic Warning with a Two-Station Subarray［J］.Bulletin of the Seismological Society of America，94(4):1546－1550.

Satriano C.，Lomax A.，Zollo A..2008.Real-time evolutionary earthquake location for seismic early warning［J］.Bulletin of the Seismo-logical Society of America，98(3):1482 －1494.

Voronoi G..1908.Nouvelles Applications des Parameters Continus，a la Theorie des Formes Quadratiques.Deuxieme Memorte:Recherches sur les Parrallelloedres Primitifs［J］.Jounal fur die Reine und Angewandte Matnematik，134:198 －287.

Wu Y.M.，Kanamori H.，Allen R.M.，et al.2007.Experiment using the tau-c and pd method for earthquake early warning in southern California［J］.Geophys.J.Int.，170:711 －717.