基于开源GIS的活动断层矢量数据融合方案初探

2017-12-02中国地震局地质研究所中国地震局活动构造与火山重点实验室吴熙彦

电子世界 2017年22期

中国地震局地质研究所中国地震局活动构造与火山重点实验室吴熙彦

中国地震局地质研究所中国地震局活动构造与火山重点实验室吴熙彦

活动断层是确定区域内最大发震震级和破坏强度的重要资料。随着计算机技术的发展，在国家支持下，地震活动断层调查相关项目也获得了大量活动断层调查空间数据。将这些数据进行综合，提取关键信息，形成一套格式、结构统一的数据，有利于提高效率和进一步的数据分析。

矢量数据融合；活动断层；数据；pyshp

一、全国活动断层数据融合意义

中国位于世界两大地震带——环太平洋地震带与欧亚地震带之间，地震频发，造成的灾害损失巨大。然而，地震短临预报仍是世界难题，短期内无法攻克。目前，国际和国内防震减灾的一项主要措施是通过确定某区域可能发生地震时造成的破坏强度，设定抗震设防烈度，规定建构筑物要达到的抗震要求，减少地震灾害造成的损失。活动断层则是确定区域内最大发震震级和破坏强度的重要资料。随着计算机技术的发展，在国家支持下，地震活动断层调查相关项目也获得了大量活动断层调查空间数据。这些数据获取的时间跨度大，来源不一。因此，数据的格式、结构和制作标准等并不一致。使用时需要调用和查看不同的数据集，调用不同的图层，效率较低。同类数据分层分包存储，也不利于空间查询、统计等分析研究。将这些数据进行综合，提取关键信息，形成一套格式、结构统一的数据，有利于提高效率和进一步的数据分析。

二、全国活动断层数据融合流程设计

不同时期、不同来源的活动断层GIS数据融合的内容包括空间图形和属性值[1]。将数据融合的思路分解为图形的融合和属性值的融合，以先图形后属性的方式进行流程设计。数据融合的步骤如下。

（1）转换文件格式：原有数据采用了不同的文件格式。数据处理的第一步，是将所有数据统一转换为shape file文件。Shapefile文件是被广泛接受的一种GIS数据文件。常用的GIS软件都能将GIS数据转换为这种格式。

（2）转换空间坐标系：不同来源的GIS数据可能在不同的空间坐标系下，应在处理之前转换到统一的坐标系下。目前，我国最新的坐标系统为2000国家大地坐标系，它是国家基础地理信息数据（基础地理底图）所采用的坐标系，由国务院批准自2008年7月1日启用。为推广2000国家大地坐标系的使用，国家测绘地理信息局公布了现有测绘成果转换到该坐标系的技术指南。依据技术指南，可以将现存数据都转换到2000国家大地坐标系。数据转换到统一坐标系后，可实现多来源空间数据叠加，进行编辑和修改。

（3）分区处理数据：空间数据叠加后，可以看到部分区域内只有一层图层，每个对象只有一种图形表达，而部分区域有一层以上的重叠数据。对于只有唯一空间数据的区域，不需进行空间操作，直接保留。对于重叠区域，须进行图形合并。因此，在图形合并前，需将将数据集切块，切出重叠区的部分待下一步处理。

（4）重叠区图形数据合并：地图合并技术一般包括匹配实体和将匹配的同名实体合并两个过程[2][3]。实体自动匹配主要依靠几何相似度、拓扑、数据模型和属性等。同名实体自动合并主要根据预设权重和待合并图形几何形态生成新的图形。然而，以线要素表达的活动断层线的几何形态并非地质对象的真实形态，即地壳断裂带的真实几何形态，而是地质学家对地质对象的一种解释和抽象。同一实体的不同解释和抽象可能使几何形态相似程度较低。因此，活动断层的几何形态的，并不具备很好地自动匹配和自动合并的条件。因此，采用以下三步合并数据。1）选取重叠区某一图层为最终保留层（称基准图层），用人机交互的方式修改该图层；2）为每一图层增加“同名实体编号”字段。所有同名的实体赋予唯一的一个编号填入该字段；3）利用“同名实体编号”连接其他数据属性表到基准图层，将同名实体的属性全部合并入基准图层。最终重叠区仅保留基准图层一层数据。

（5）确定属性表结构：合并重叠区和非重叠区的空间数据，形成一层新的图层。新数据的属性表是原有数据所有属性的集合。新属性表存在重复信息，需要重新设计与规划。根据属性内容，重新设计属性表结构，在不损失信息条件下精简表结构。

（6）属性融合：根据重新设计的表结构，为属性表新增字段，准备填入信息融合后的属性值。属性处理有3种情况。1）属性值直接赋给新表中对应属性字段；2）多个信息重复的属性值，合并到新表中的一个属性字段；3）多个信息不同的属性值，合并到新表中的一个属性字段。根据新形成的数据和新设计的表结构，设计属性迁移和合并规则。

（7）删除冗余字段：将上一步形成的新数据中的冗余字段删除，获得最终成果数据。

图1 流程图

三、方案可行性分析

1.pyshp介绍

方案的第1-5步都较容易实现。方案实现的关键是第6步属性融合。笔者采用了pyshp作为实现第6步的试验工具。pyshp全称Python Shapefile Library，是基于python语言的开源函数库，提供了读写ESRI shapefile格式文件的方法。pyshp函数库具有编译工具简单、免安装、代码开源、使用简单的特点，并且提供了简单易懂的使用说明和示例代码，在ESRI shapefile编辑方面便利，非常适合地理数据融合工具。pyshp提供三个基础类:Reader、Writer、Editor。Reader 提供了ESRI shapefile的文件头信息、几何图形和属性表的读取方法；Writer提供了ESRI shapefile的文件头信息、几何图形和属性表的写入方法；Editor是Writer子类，继承了Writer的所有方法，主要用于已有文件的编辑。pyshp可以在任何pyshon平台直接使用import命令导入模块文件调用函数。本文使用的是pyshp 1.2.3。

2.关键步骤的实现

这里仅讨论关键的两种属性处理。

（1）多个信息重复的属性值，合并到新表中的一个属性字段。例如，断层活动年代的编码方式不同，一部分代码值为{a1,a2,a3,...ai}，另一部分代码值为{b1,b2,b3,...ai}。实际上，a1和b1代表相同属性。该属性的新代码值为{c1,c2,c3,...ci}。此时，设计一个代码映射表，将代码一一映射到新代码集。核心python代码如下：

iSetField=self.getEsfFieldIndex(esf, SetField)#新字段

iReferField=self.getEsfFieldIndex(esf, ReferField)#旧代码值所在字段nRecord=0

for xRecord in esf.records:

#如果新代码字段是空的，就将旧代码值的对应新代码填入

if xRecord[iSetField] is None or xRecord[iSetField]==0:

reSValue=str(u(xRecord[iReferField]))

xRecord[iSetField]=int(AgeDic.get(reSValue,0))#从代码映射表获得新代码

nRecord+=1

return nRecord

（2）多个信息不同的属性值，合并到新表中的一个属性字段。例如，原数据有一个表示断层是“出露”还是“隐伏”的字段，还有另一个表示“正断”还是“逆断”的字段。两种信息可以合并，成为新的断层性质字段。新编代码的前几位代表前一个字段值，后几位代表后一个字段值。核心python代码如下：

for xRecord in esf.records:

Sfirvalue=str(u(xRecord[iFirField])).strip() #获得第一个字段值