趙 瑞，劉希強

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

0 引言

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

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

在實時地震監(jiān)測臺網(wǎng)中，當(dāng)在臺站中斷記錄時段發(fā)生地震，利用V圖進(jìn)行地震定位必須要調(diào)整V圖布局，若考慮所有正常臺站重新計算V圖則復(fù)雜度會明顯加大，如何能在不改變?nèi)諺圖的前提下，快速完成對中斷臺站V圖單元的重新剖分，又不失精度，是地震預(yù)警系統(tǒng)中需要解決的問題。

1 V圖生成算法

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

1.1 初始凸殼的生成

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

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

1.2 點集預(yù)處理

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

圖2 點集預(yù)處理(a)M點位于R1區(qū);(b)N點位于R1區(qū)中①區(qū);(c)N點位于R1區(qū)中②區(qū);(d)N點位于R1區(qū)中③區(qū);(e)M點位于R1區(qū);(f)M點位于R2區(qū);(g)M點位于R3區(qū);(h)M點位于R4區(qū)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 逐點插入算法

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

1.4 Delaunay三角剖分優(yōu)化算法

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

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

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

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

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

1.5 V圖生成算法

存儲每個臺站所涉及三角形的個數(shù)及編號，按一定順序?qū)θ切瓮饨訄A心進(jìn)行排序連接，即得到此臺站的V圖單元。同理，依次遍歷所有臺站，得到研究區(qū)域的V圖見圖5。

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

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

在實際地震監(jiān)測臺網(wǎng)運行中，地震臺站通常會因為中斷或其他原因未記錄到波形信號，如果地震恰發(fā)生在中斷臺站所在的V圖單元內(nèi)，則根據(jù)首先觸發(fā)臺站判定地震發(fā)生的區(qū)域時會產(chǎn)生較大的震中判定誤差，從而影響地震定位的精度。所以在地震發(fā)生時需要實時判定正常運行臺站的數(shù)量，快速準(zhǔn)確地對中斷臺站的V圖單元進(jìn)行調(diào)整和歸并。

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

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

圖6 中斷臺站形成的空區(qū)圖以及四種不同情況下三角形形成示意圖(a)空區(qū)圖;(b)中斷臺站處于網(wǎng)內(nèi)時，形成封閉區(qū)域情形1;(c)中斷臺站處于網(wǎng)內(nèi)時，形成封閉區(qū)域情形2;(d)中斷臺站處于網(wǎng)緣時，形成半封閉區(qū)域情形1;(e)中斷臺站處于網(wǎng)緣時，形成半封閉區(qū)域情形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所示，中斷點設(shè)為點O，選取空區(qū)中任意兩條相交邊AB、BC，當(dāng)AO、AC在AB邊同側(cè)且CA、CO在CB邊同側(cè)，即T5與T6同號，T7與 T8同號時，連接 AC形成新三角形△AOC、△ABC(見圖6b、6d)，反之不能形成新的三角形，如圖6c、6e所示。

(2)當(dāng)空區(qū)內(nèi)所有三角形形成后，在空區(qū)內(nèi)對其進(jìn)行優(yōu)化，使之成為Delaunay三角形。針對圖6a中存在5個中斷臺站的情形，應(yīng)用上述思路，生成的新的Delaunay三角剖分圖見圖7a，3個空區(qū)新生成的三角形用不同顏色標(biāo)識出。

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

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

3 計算時間測試

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

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

表1 Delaunay三角網(wǎng)格剖分算法改進(jìn)前后對比表Tab.1 The comparison of delaunay triangulation algorithm before and after improvement

第二種是模擬山東地區(qū)存在5個中斷臺站，除去中斷點，其余74個臺站點重新剖分生成Delaunay三角網(wǎng)格與本文提出的局部采用改進(jìn)后的三角網(wǎng)格生長法所用的時間進(jìn)行對比分析。前者生成的Delaunay三角網(wǎng)格如圖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圖生成算法優(yōu)于傳統(tǒng)的逐點插入法，在實際臺網(wǎng)中，臺站數(shù)量較少的情況下，基本上是瞬間完成的，并且地震臺站中斷時，局部采用改進(jìn)后的三角網(wǎng)格生長法與除去中斷點、重新剖分Delaunay三角網(wǎng)格得到的結(jié)果是一致的，但在計算用時上優(yōu)于后者。

4 結(jié)語

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

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

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

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

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

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

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

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

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

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

殷海濤，劉希強，李杰，等.2012.現(xiàn)今地震預(yù)警技術(shù)及其在國內(nèi)發(fā)展?fàn)顩r的探討［J］.中國地震，28(1):1－9.

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

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.