基于GIS和GDAL的地质灾害汛期风险预警产品信息自动化生成算法的构建与实现

张像源

（天津市地质调查研究院，天津 300191）

0 引言

我国地质灾害种类多、分布广、危害大。为进一步减轻地质灾害风险，最大限度保障人民群众生命财产安全，从2003年起，我国大陆开始开展汛期区域地质灾害气象预警工作，并形成预警产品面向公众的发布,在地质灾害防治中发挥了 “消息树”和“发令枪”的作用,取得了明显的社会经济效益[1],极大地提升了公众社会对防范地质灾害的认知。随着该项工作的不断推进，预警产品的内涵从粗到细逐步走向规范化，包括了预警范围、等级、时段和文字说明等内容[2]，且临灾发布工作要求快速、高效和精准。然而这些信息若靠人工获取并不是件易事，比如要经技术人员通过读取预警区划图并研判后，才能粗略得出空间分布范围等定性信息，这一过程耗时耗力，且难以达到定量描述的准确效果。因此自动计算预警结果并快速生成预警产品是预警系统的一个重要功能[2],围绕聚焦解决好地质灾害可能发生的地点、成灾范围等预警预报问题[3-4]，考虑如何利用信息技术自动获取预警产品描述信息的实现势在必行。然而经检索发现国内直接进行相关研究还较少，可供参考的文献不多。为了快速形成权威、科学、符合实际的这一产品，本文提出一种可高效自动分析形成初步的预警产品描述信息的技术算法，为有关信息平台的功能模块的研发、支撑专家做出更详尽的预警产品研判提供理论支撑。

1 算法原理

1.1 技术路线

以研究区行政区划图和预警区划图[2]为数据来源，应用GIS 技术将空间位置信息和属性信息无缝结合，结合数学统计和地理知识等，精确获取预警等级在行政区划所处方位、区域占比和防治措施等精细描述的预警信息，将图面内容转换为直观的文本描述，为最终预警产品的快速生成奠定基础（图1）。

图1 技术路线Fig.1 Technical route

1.2 关键数据

为了实现目标，在开始研究前，要准备必要的空间图层作为计算的基础，即将具有不同预警等级属性的矢量栅格地质灾害预警区划图层和反映行政区名称及空间范围的研究区行政区划图层作为数据源，二者要求具有相同的地图投影参数。其中，前者一般通过网格剖分后通过多因子要素叠加进行综合评价得到，这一过程在当前地质灾害空间评价预警研究中是普遍采用的方法[5-6]，但因不是本文的重点，故不再赘述。

1.3 预警描述信息各因子获取

1.3.1 预警等级集合概念模型

很显然，前述空间图层叠加结果中，不同预警级别对应的地区预警信息处于离散状态，且一般情况下，这些信息数量较为可观，为了后续分析数据高效便捷，有规律可循，需要提前分析预警等级、地区、方位等数据间的关联性，通过聚类分析建立分类簇集合信息，这一预警等级集合概念模型设计见图2。

图2 预警等级集合概念模型示意图Fig.2 Schematic diagram of the conceptual model of early warning level aggregation

图2中可以看出，预警等级集合以三级预警等级为唯一的主键，包括了这一等级下的地区列表预警信息子集合，而该子集合以地区名称为主键，包括了一对多的分布方位和面积占比列表，该列表以分布方位为唯一索引，在实际计算时，相同分布方位需要进行去重处理，所占的面积要进行求和运算。最终通过综合计算，形成该预警等级下的总体分布区域、分布方位等综合预警描述信息。

1.3.2 空间信息叠加与判别

空间关系描述是GIS 系统的基本功能之一，GIS 的技术支持的地质灾害风险区划的最终目的是划分不同灾害等级的区域，可为地质灾害预警提供依据[7-8]。通过将预警区划图层和行政区划分区图层进行空间叠加，遍历判断每个预警等级矢量栅格单元和行政区划单元的空间拓扑关系，确定该行政单元是否包括某预警级别（图3）。如果二者拓扑关系为不相离，说明该地区具备该等级，反之如果是包含、相交、穿越等非相离关系，则还需进行两两拓扑求交运算，并重新采用交集中的预警等级区参与面积计算更具科学性和精确性。

图3 空间判别示意图Fig.3 Schematic diagram of spatial discrimination

1.3.3 预警等级分布面积占比求算

一个由N个拐点（xi,yi）确定的封闭多边形的面积如式（1）计算：

式中：i——拐点序号；

N——拐点个数；

xi——第i个拐点x坐标；

yi——第i个拐点y坐标；

A——封闭多边形的面积。

相同方位的预警等级所占行政区划单元的面积占比（R）是对预警等级广泛程度的描述,算式如下：

式中：i——拐点序号；

N——拐点个数；

R——某预警等级的面积占比；

Ai——某预警等级单元的面积；

1.3.4 预警等级分布范围描述

对分布范围则用绝大部分、大部分、局部、个别四级进行空间范围的广泛程度描述，判别指标为集合中单元格之和的占比（R），定义见表1。

表1 预警等级占比描述表Table 1 Description of the proportion of early warning levels

1.3.5 预警等级单元分布方位求算

空间方向的定性描述是用若干主方向粗略地描述空间方向。而定量描述则是用方位角来量测空间目标之间的方向关系，因此方位角是空间方向描述的一个重要手段[9]。要获取一个预警等级单元在行政区划中的分布方向，实际上是通过计算该等级相对于所处区域的方位角得到（图4），结合地理知识，根据实际情况共划分出了8 个方位角区间和对应的分布关系（图5）。

图4 方位角示意图Fig.4 Schematic diagram of the azimuth angle

图5 方位角与分布方向的映射关系图Fig.5 The mapping relationship between azimuth and distribution direction

1、方位角计算

（1）形心求算

本次方位角计算要获取行政区划单元和预警等级单元的几何形心。一个由N个拐点（xi,yi）确定的封闭多边形的中心如式（3）、式（4）计算：

2．2．1 所获病例的Apgar评分构成比 所获265份病例中，重度窒息组中，生后1min Apgar评分为1分、2分和3分的分别有2例（0．8%）、3例（1．1%）和27例（10．2%）；轻度窒息组中，生后1min Apgar评分为4分、5分、6分和7分的分别为11例（4．2%）、19例（7．2%）、56例（21．1%）和120例（45．3%）；而对照组中，生后1min Apgar评分为8分和9分的分别有25例（64．1%）和2例（5．1%）。

式中：i——拐点序号；

N——拐点个数；

A——多边形的面积，由式（1）得出；

xi——第i个拐点x坐标；

yi——第i个拐点y坐标；

Cx——几何形心x坐标；

Cy——几何形心y坐标。

（2）方位角求算

当在平面上2 个点的坐标已知时，给出方位角（十进制度数）公式如式（5）：

式中：X2——平面上终点x坐标；

X1——平面上起点x坐标；

Y2——平面上终点y坐标；

Y1——平面上起点y坐标；

α——方位角。

（3）分布方向描述

当方位角得出后，便可以根据方位角与分布方向的映射关系图（图5）获取具体的分布方向。

1.4 预警描述信息提取

经过前述步骤后，将得到一个以预警等级为索引的数据集合，其中包括各个地区的预警信息子集合。将地区预警信息子集合按照预警方位的个数进行升序排列，为了简洁，可以选取前若干个方位作为主要的方位，其余的则以“等地区”代替,而该级别的防治措施则从表2中对照获取[2]。

表2 预警等级防治措施描述表Table 2 Early warning level control measures description table

1.5 形成预警描述综合信息

通过遍历预警等级集合，形成该地区综合的预警产品描述信息。最终某地区预警信息描述格式形成通用模板举例如下：风险大（Ⅱ级）主要分布于××地区的西北部、北部、东北等局部地区，请加强防范；风险较大（Ⅲ级）主要分布于××地区的西南、南部等个别地区，请注意防范。

2 应用实现

2.1 数据准备

考虑到预警信息的范围划分取决于行政区划图，不同的划分将得到不同的预警产品描述信息，故可将行政区划图看作一个相对的变量，为了验证工作的灵活性和可扩展性，本次通过GIS 软件形成了包括11 个虚拟乡镇的地区作为预警范围实验数据进行工作（图6），该地区某日的矢量栅格地质灾害预警区划图层业已形成（图7），采用3 km×3 km 为预警等级单元，预警单元总数为5 194 个。上述文件均为Shapefile 格式，其投影参数为CGCS2000。

图6 某地区地质灾害风险预警范围Fig.6 Geological hazard risk warning range in a region

图7 某日地质灾害预警区划图层Fig.7 Geological disaster early warning zoning map on a certain day

2.2 基于GDAL 信息模块的研发

GDAL 库(Geospatial Data Abstraction Library)是一个开源的用于栅格和矢量地理空间数据格式的C++转换器库，由开源地理空间基金会在MIT（麻省理工学院）风格的开源许可下发布，目前几乎所有的GIS 和RS 软件底层都使用GDAL 来读写空间数据[10]。由于GDAL库能够很好地支持包括Shp 格式在内的很多数据格式，因此本次工作基于GDAL 和C++语言，在Visual Studio 2019 平台下实现了整个研发验证过程。研发过程中，以本算法为理论基础，结合UML 理论建立了空间分析类、因子获取类、信息处理类，并形成了类的方法、属性和事件（图8），最后编译形成了预警信息生成软件工具。

图8 模块构成图Fig.8 Module composition diagram

2.3 预警产品信息结果

利用已有的实际数据图层和软件工具，耗时不足1 分钟形成了预警产品信息（表3）。从表中可以看出，该日预警级别包括了Ⅱ级、Ⅲ级，主要是以风险较大（Ⅲ级）为主，分布于秦姜董镇、王家湾镇、白堆子镇等地，并且有方位和分布范围程度描述。依据表中信息，通过组合可形成规范格式的综合预警产品描述信息，见图9。

图9 预警产品描述信息Fig.9 The description information of the warning product

表3 某地区预警产品描述信息表Table 3 Table of product description for a specific region

2.4 与传统人工分析对比验证分析

为了进一步验证算法的正确性和优越性，将风险预警范围图（图6）和预警区划图层（图7）进行叠加，通过逐一人工判断每个乡镇区域预警区划情况，可看出所有的乡镇的结果和预警区划图层（图7）的空间展现结果完全吻合，如王家湾镇预警产品信息为风险大（Ⅱ级）、分布于西北，占比不大（图10），这一目估结果与采用本算法生成的预警产品描述信息一致，其次应用计算过程耗时短暂，可见算法能完全满足汛期地质灾害风险预警之快速高效的需要，与传统人工定性分析相比，具有规范、快速、准确的特点，没有遗漏，其优越性显著。

图10 预警产品描述信息验证Fig.10 Early warning product description information verification

3 结论

文章提出了基于GIS 技术支撑下的地质灾害风险预警产品描述信息的自动化生成技术方法，并以随机的行政区划图和已有的地质灾害预警区划图为数据源，采用GDAL 库，通过研发预警信息获取模块进行了实例验证。结果表明，采用本算法形成预警产品描述信息快速高效，结果科学规范、完整全面，进一步提高了预警信息精度，不仅大大降低了预警产品信息获取的繁琐程度，而且节省了时间成本，又兼顾了产品的准确性，可为相关信息系统功能研发及专家研判提供基础理论支撑，将显著提高预警精细化程度和工作效率，适合在实际工作推广使用。