三角形网络和判据图法在地质灾害气象预警中的应用

王正维，陈迎辉，楚国华

（1.河北省张家口地质环境监测站，河北张家口 075000；2.石家庄经济学院，河北石家庄 050031）

0 引言

自2003年6月中国气象局和国土资源部合作开展全国地质灾害气象预警工作以来，各省市陆续开展了此项工作，在很大程度上减少了人员伤亡和财产损失，为防灾减灾工作开辟了一条新的途径，已成为保护广大人民群众生命财产安全的重要手段。但地质灾害气象预警工作的理论基础和技术方法仍在进一步研究完善过程中。在理论研究方面刘传正等对我国地质灾害的预警方法开展了深入系统的研究，提出了全面考虑地质灾害潜势度、危险度和危害度的“三度”概念［1、2］。目前国内开展此项工作的主要技术方法有考虑前期过程雨量的综合判据图法、网格剖分概率量化综合预警模型法、考虑降雨特征量的预警指标法以及基于GIS和神经网络的地质灾害预警系统等［3］。由于地质灾害本身的复杂性，预警理论方法还处于研究探索阶段，各种方法应根据实际应用效果加以调整。本文以河北省张家口市为例，以改进的综合判据图法为理论依据，利用Delaunay三角网和边界限定自动剖分三角网进行地质灾害预警的自动计算、判别，并开发了预警软件。

1 张家口市地质灾害气象预警工作现状

张家口市位于河北省西北部，面积3.7×104km2，既有山区、盆地，又有高原、丘陵，地貌及地质构造复杂，地质环境脆弱，地质灾害具有种类多、分布范围广、突发性强、危害性大等特点，是地质灾害多发区。截止2008年底，已发现各类地质灾害及隐患708处，主要威胁着村庄、居民、农田、矿山、水库、学校、交通干线、乡间公路和水利设施等，在很大程度上影响着经济建设的快速发展。因此，在张家口市开展地质灾害气象预警工作意义重大。

从2005年6月开始，张家口市国土资源局和张家口市气象局联合开展了汛期地质灾害气象预警工作。市气象台提供雨量实测和预测资料，张家口市地质环境监测站根据气象台提供的实测预测雨量结合相关地质资料进行综合研究分析，在与各有关方面充分会商后做出地质灾害预警预报。根据掌握和收集的大量资料，将全市划分了16个地质灾害预警区。根据已经发生的地质灾害单体或群体前期累计过程雨量，参考相邻区域和地质环境条件相似区域的预警判据，确定了各个预警区的地质灾害发生的临界过程雨量判据。在2005～2006年的地质灾害气象预警工作中，采用各雨量站点预警日前期过程雨量累加，与雨量站点所在地质灾害预警区预警判据对比，分析判断后得出地质灾害预警级别。由于采用的网络数据传输、计算机软件辅助计算、预警级别人工判别方法，速度慢、效率低，计算和判别过程误差大、自动化程度低，因此，需要开发一套预警全过程自动计算判别软件，实现预警过程的自动化，以提高计算速度、准确度，消除人工判别误差，提升地质灾害气象预警工作效率。

2 预警理论方法

张家口市地质灾害方面研究程度低，各方面资料较少，因此，实现按网格剖分的地质灾害潜势度预警方法比较困难，所以，考虑当日雨量与累积雨量的综合判别表方法作为预警的基础理论依据［4］。根据以往的研究成果，将张家口市划分为16个地质灾害预警区，进而确定了各个预警区地质灾害发生的临界雨量判据。

目前，国内从国家级到省级的地质灾害气象预警业务普遍采用前期累计雨量、有效累计雨量与最新日雨量等预警指标，张家口市地质灾害气象预警分析雨量采用预警日当日预测雨量和2、4、7、10、15日前期累计雨量进行综合预警分析。按照国土资源部和中国气象局的联合规定，地质灾害气象预警分为5级。

国内最新研究成果表明［5-7］，地质灾害发生的规模和类型与降雨量有明显的正相关，而不同降雨类型下灾害的发生情况也存在明显差异，地质灾害的发生存在明显的滞后性，一般集中在降雨后72h内发生，从理论和国内有关观测事实来说，地质灾害发生的当日激发雨量对地质灾害的发生有着重要影响，5天以内降水情况对地质灾害的影响随时间延续而明显减弱，超过5天以上的前期降雨情况对地质灾害发生的影响很小。因此，作者采用的6个时间段根据预警判据计算的预警结果对预警日综合预警结果取值影响程度各不相同，距预警日时间越长，影响越小。本文使用权重系数来表示前期6个时间段的预警计算结果对综合预警结果的影响程度，将计算取得的各预警级别各时间段的权重系数之和称为预警指数。只有当预警指数达到或超过某一设定值时，即达到预警指数临界值，该预警级别计算结果才有效，这个计算过程可以有效避免以往综合预警级别取各时间段预警结果极大值时距离预警日较远时间段的降雨情况对综合预警结果的过大影响，体现距离预警日越近的降雨情况对综合预警结果影响越大的特点。经过对相关资料的研究成果的分析，初步为6个时间段预警结果分别设定不同的权重系数，见表1：

综合预警结果取值时，按预警级别将6个时间段的预警计算结果权重进行相加，如果6个时间段预警结果有一项以上为某一预警级别时，则要在该级别的权重之和上加上基本权重，最后取各个预警级别中预警指数（即总权重）达到预警指数临界值的级别中最高级别为预警日综合预警级别。根据6个时间段预警结果权重和预警指数临界值取值可以看出，当日计算结果对综合预警级别取值起决定性作用，5日以后时段预警结果对综合预警级别取值作用很小。

表1 预警指数临界值与权重系数Table 1 Critical value of early warning index and weight coefficients

Delaunay［8］三角形网格具有算法简单、生成容易、有较强的控制性、单元具有较好的形状和较强的边界适应性，在二维空间计算方面能够表现出良好的性能，因此，作者在对张家口市地质灾害气象预警采用Delaunay三角形网络剖分方法，充分利用Delaunay三角形网络剖分方法，充分利用Delaunay三角形网络的性能优势进行等值线一成和图形计算。

以张家口市界为边界限定条件，对预警区域进行Delaunay三角形网格边界限定自动剖分，剖分取得的等边三角形网络与各预警区多边形进行区对区空间判别计算，用自动剖分三角形的形心代表整个三角形区域。通过上述计算，将预警区的数据信息转换到自动剖分三角形的形心上，可以得出每一个计算三角形的6个时间段雨量预警判据。

对张家口市域内的210个自动雨量站点进行凸域点集Delaunay三角形网格自动构网，建立预警区的雨量站点三角网。在具体计算时，从雨量数据库中取出前期过程雨量，利用三角网的网对网计算进行预警三角网空间插值，计算出每一个预警三角形的过程雨量，进而计算出各时间段预警级别和综合预警级别。

对预警日综合预警级别计算完成后，预警结果按等级进行三角形区域拓扑融合计算，确定各个预警区域级别、面积、边界后，进行计算机输出和保存［9］。

3 预警方法的实现与软件开发

3.1 雨量数据库的设计与实现

雨量数据库用于存储各自动雨量站点的实测与预测雨量，将其按日期顺序设计为固定长度的数据存储区块，以预警计算日为关键字段，实现雨量数据按日期快速随机定位、查询、存取，数据库设计为固定长度，设置雨量数据初始值，可以在读取过程雨量时判断数据存储状态。

3.2 预警区的建立

将张家口市现有16个预警区作为基本预警区（图1中实线表示），基本预警区的外部包容边界为市界，为了便于进行市界边界处的自动剖分三角网对预警区的区对区判别计算，将市界外侧一定距离范围内的区域作为扩展预警区，扩展预警区各项数据与相邻的基本预警区相同。

图1 预警区与雨量站点三角网Fig.1 Early warning area and triangular networking of rainfall stands

3.3 雨量站点三角网的建立

将全市232个自动雨量站点按凸域二维点集合的形式进行Delaunay三角形自动构网（图1中虚线表示）。利用雨量站点三角网可以实现对自动剖分三角网的三角形雨量插值。

3.4 自动剖分三角网的实现

考虑到计算机性能和实际预警精度需要，设定三角形边长为2.5km，每个三角形的面积2.7km2，以张家口市界为限定边界，建立自动剖分等边三角网，对所有自动生成的三角形进行空间多边形（三角形）对市界多边形自动检测判别，保证所得三角形与基本预警区的相关性，实际运算得到1.6万多个三角形。

图2 程序设计主要流程图Fig.2 Process of the program design

3.5 预警三角网的实现

预警三角网是在自动剖分三角网的基础上构建的，将自动剖分三角网与所有预警区多边形进行区对区相交判别分析，得到预警三角网。设计预警三角网的数据结构不存储自动剖分三角网数据项中的位置与空间关系数据，而是建立自动剖分三角形与预警三角形的相互指针链接，可以快速实现两个三角网中相对应三角形的数据交换，减少内存占用，提高程序运行速度。

3.6 软件设计的主要流程

软件开发的基本思路是以雨量数据库、雨量站点三角网和地质灾害预警三角网为基础，雨量数据存储于雨量数据库中，利用雨量站点三角网计算出各节点各时段雨量值，进一步自动生成等值线；利用雨量站点三角网对地质灾害预警三角网形心的自动插值，计算出每个地质灾害预警三角形各时段的累计雨量值，通过与预警判据对比计算出各时段预警级别，进而计算出综合预警级别，软件设计的主要流程见图2。

3.7 软件开发与操作

以Microsoft Visual C++6.0为开发工具，借助MFC类库的强大功能进行软件开发，为加快运行速度，提高程序灵活性，摆脱商品软件二次开发的使用限制，全部代码由软件自身实现，数据库格式使用自定义格式。

图3 程序主要功能模块Fig.3 Main function modules of the program

地质灾害气象预警程序由建立预警模型、导入雨量数据、预警计算和显示输出四个功能模块组成（图3），建立预警模型和导入雨量数据是地质灾害预警的基础，预警计算模块实现了该程序的核心功能显示保存模块实现了预警计算结果的对外输出保存功能。

图4 2008年7月29日实测15日降雨等值线Fig.4 Contours of 15-day rainfall measured on July 29，2008

软件操作主要分为五类：基础数据库操作、雨量数据库操作、预警计算、预警结果的显示存储、辅助工具选项。

基础数据库操作主要包括雨量数据库的建立、预警区数据文件的建立、边界限定自动剖分三角网的生成、预警三角网的构建、雨量站点三角网的构建。通过菜单操作，控制参数从文本文件输入，只要输入相关数据文件名称，便可完成操作。雨量数据库操作主要包括成批导入实测雨量、成批导入预测雨量、导入实测雨量和导入预测雨量，这此操作将文本文件的数据导入雨量数据库中。

图5 2008年8月8日预测10日降雨等值线Fig.5 Contours of 10-day rainfall estimated on August 8，2008

预警计算部分可以实现指定预警日期的地质灾害自动预警，还可以实现历史日期的预警结果显示。

显示编辑和预警/等值线部分主要是对预警计算结果的显示输出，具有放大、缩小、原始尺寸显示、设置缩放系数、存储图像、显示图例等功能，显示内容可以选择6个时间段的实测或预测雨量等值线、预警结果，可以同时选择显示基础界线、雨量站点三角网、预警三角网。显示的内容均可以保存为JPEG格式的图像文件。

所有的参数可以通过配置文件修改，软件设计可以适应不同地区、不同条件的地质灾害预警工作。

4 软件应用

使用该软件系统进行了张家口市2007年和2008年各136天的汛期地质灾害气象预警，软件运行正常，预警结果与人工计算预测结果基本相同，共预警9次，其中气象部门实测降雨量与预测降雨量基本相同的6次均成功预警，避免了人员伤亡和经济损失。在自动剖分三角形边长2.5km、三角形数量1.6万多个的条件下，使用P4 3.0GHz双核处理器、512M内存的计算机进行计算，预警计算时占用内存43M，计算时间1分05秒。实践证明，该软件可以满足市级地质灾害预警工作的时间要求和精度要求（图4～图7）。

图6 2008年9月2日预警等级计算结果Fig.6 Rand of forewarning by calculation on September 2，2008

图7 2008年9月8日预警等级计算结果Fig.7 Rank of forewarning by calculation on September 8，2008

5 结论

本文在地质灾害气象预警的综合判据图法的基础上，对地质灾害预警方法进行了改进，对不同时段按综合判据图法计算的预警等级给予不同的权重系数，并引入预警指数和预警指数临界值进行综合预警等级计算，体现了距离预警日越近的降雨情况对综合预警结果影响越大的特点。

使用Delaunay三角网理论方法进行地质灾害预警区边界限定自动剖分、构建雨量站点三角网，利用Delaunay三角网的特性进行雨量自动插值、等值线追踪、预警结果计算和输出，开发了预警软件，采用与数据存储结构相适应的优化算法，提高了程序运行速度，实现了各项基础数据库建立、雨量数据库存取、预警计算全过程自动化，参数使用配置文件读取，软件模型应用灵活、适应性好，在基础资料充分的情况下、软件可以扩展构建成三角形单元网格剖分概率量化综合预警模型。

软件未进行不同地区和不同条件下的测试，各时间段的预警结果权重系数、预警指数、预警指数临界值的取值方法和大小还需要在应用中进一步探讨。

受气象部门降雨量预报精度和准确度限制，预警日预测雨量只能取预报的可能降雨量值的最大值，在一定程度上降低了地质灾害气象预警精度和准确度。在得到气象部门等相关方面技术支持的条件下，软件还可以实现自动导入各雨量站点实测和预测雨量、预警结果网络传输等功能，进一步提高地质灾害气象预警工作的效率和准确度，更好的为地质灾害防治工作服务。

［1］刘传正.地质灾害预警工程体系探讨［J］.水文地质工程地质，2000，27（4）：1-4.

［2］刘传正.突发性地质灾害的监测预警问题［J］.水文地质工程地质，2001，28（4）：1-4.

［3］李晓.重庆地区的强降雨过程与地质灾害的相关分析［J］.中国地质灾害与防治学报，1995，6（3）：39-42.

［4］陈百炼，杨胜元，杨森林，等.基于GIS的地质灾害气象预警方法初探［J］.中国地质灾害与防治学报，2005，16（4）：93-96.

［5］钟荫乾.滑坡与降雨关系及其预报［J］.中国地质灾害与防治学报，1998，9（4）：81-86.

［6］潘懋，李铁锋.灾害地质学［M］.北京：北京大学出版社，2002，10-11.

［7］傅鹤林，周中，卜翠松，等.地质灾害预测预报国内外现状［J］.湘南学院学报，2006，27（2）：43-48.

［8］邬吉明，沈隆均，张景琳.Delaunay三角网格的一种快速生成法［J］.数值计算与计算机应用，2001，22（4）：267-275.

［9］章孝灿，黄智才，戴企成，等.GIS中基于拓扑结构和凸壳技术的快速TIN生成算法［J］.计算机学报，2002，25（11）：1212-1218.

［10］周培德.计算几何——算法分析与设计［M］.北京：清华大学出版社，2000，103-109.

［11］卢朝阳，吴成柯.任意多边形内带特征约束的散列数据的最优三角剖分.计算机设计与图形学学报，1997，9（4）：302-308.