建设地震监测台网的主要科学目的之一就是一旦地震发生后,第一时间开展地震速报工作。当前基于实时地震监测台网建立地震预警系统和地震应急控制系统在国际上已形成潮流。日本、美国、墨西哥、土耳其、意大利、瑞典、罗马尼亚、中国台湾等地震多发国家和地区发展了多套预警系统,如日本铁路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图单元的重新剖分,又不失精度,是地震预警系统中需要解决的问题。

在实际地震监测台网运行中,地震台站通常会因为中断或其他原因未记录到波形信号,如果

地震恰发生在中断台站所在的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)。

图6 中断台站形成的空区图以及四种不同情况下三角形形成示意图
(a)空区图;(b)中断台站处于网内时,形成封闭区域情形1;(c)中断台站处于网内时,形成封闭区域情形2;
(d)中断台站处于网缘时,形成半封闭区域情形1;(e)中断台站处于网缘时,形成半封闭区域情形2
Fig.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边同侧,即T₅与T₆同号,T₇与T₈同号时,连接AC形成新三角形△AOC、△ABC(见图6b、图6d),反之不能形成新的三角形,

如图6c、6e所示。

{T₅=AB×AC=AB.y*AC.x-AB.x*AC.y

T₆=AB×AO=AB.y*AO.x-AB.x*AO.y

T₇=CB×CA=CB.y*CA.x-CB.x*CA.y

T₈=CB×CO=CB.y*CO.x-CB.x*CO.y.(3)

(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

为测试算法改进后的实际效果,我们进行了计算时间的测试,这也是在地震预警系统中必须解决好的一个科学问题。本文程序编写采用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三角网格优化算法,提高了计算速度,并且考虑了实际台网运行过程中,地震台站中断时,V图重新生成的情况。其突出的优点是在不改变全局V图的前提下,局部区域采用三角网格生长法,完成对中断台站V图单元的重新剖分,即使在网缘台站中断的情况下,也能有效地约束凸壳的生成。通过计算时间的对比,可以看出在同等效果的基础上,计算用时要优于除去中断台站,重新剖分Delaunay三角网格。

V图在地震监测预报中具有潜在应用前景,如根据触发台站所在V图单元的集中性判定是否有地震事件发生,根据V图单元面积,搜索得到区域内相关地理信息及应急资源信息,同时结合其它台站记录信息,迅速判定震中位置,为地震早期预警系统提供有效的技术支持。

本文在撰写过程中得到国家科技支撑课题组(2012BAK19B04)和地震科技星火计划项目组(XH12029)成员的帮助,在此表示衷心感谢。