地震灾害风险普查中的危险性计算与数据处理

邵 磊，胡 刚

（1. 湖南省地震局, 长沙 410004；2. 湖南省震灾风险防治中心, 长沙 410001；3. 中国地震局地球物理研究所, 北京 100081）

0 引言

地震危险性分析是地震灾害风险普查工作中的重要环节，可为地震灾害风险评估与区划提供危险性输入[1]。根据地震危险性图编制与技术规范的要求，采用《中国地震动参数区划图：GB18306—2015》使用的潜源划分方案、地震动衰减关系、地震区带划分方案及地震活动性等参数，进行地震危险性概率计算[2]，并根据1∶100 万宏观场地类别进行场地调整、建立基岩和场地峰值加速度矢量数据库[3]。此次工作任务有国家级、省级、市级以及区县级，按照要求需建立国家级30″、省级30″或6″的标准格网；市级和区县级因地区而异，有些要求3″、甚至1″格网。格网精度的不同往往对应不同级别的数据量，在计算与后续数据处理上可能会产生一系列问题，如计算点经纬度坐标是否为格网中心位置、地震危险性概率分析是选择逐点计算还是通过插值实现、矢量数据shp 文件的制作、矢量数据间拓扑检查是否满足要求等问题。与常规的地震安全性评价中地震危险性分析工作相比，本次风险普查工作的数据计算与处理工作量巨大，往往都是十万、百万甚至千万级别，且需要按照要求建立矢量数据库，需要熟悉相关软件和一定的编程处理数据经验才可得到满足要求的结果。本文没有涉及地震危险性概率分析方法的原理，而是根据笔者在从事风险普查地震危险性计算中积累的经验、做法以及可能出现的问题，运用Matlab、Python 代码结合ArcGIS 等软件平台，主要从数据处理角度进行梳理，可省时、高效的完成该部分工作。在后续地震灾害风险普查成为常态化工作中，可为继续从事该方向的科研工作者提供一些研究思路，部分做法在地震安全性评价、地震区划工作中也有一定参考价值。

1 标准格网的建立

根据地震灾害风险普查工作的要求，一般是以某个省、市、区县的普查行政边界为工作目标区（以下简称目标区），生成所要求精度的标准格网，以各格网中心经纬度位置为计算点，进行地震危险性计算与分析。推荐使用Matlab 程序编写较短的代码即可实现标准格网的自动生成与计算点经纬度坐标的输出，主要使用函数shaperead 读取目标区边界数据矢量特征和属性，结合函数inpolygon 加以判断控制点是否在目标区内，同时使用函数shapewrite 可输出所需目标区格网矢量数据shp 文件[4]。使用程序代码操作十分便捷、运行速度快、可重复性操作效率高。为满足质检要求，需按照风险普查所要求的某个省份格网范围开始起算，程序在循环上设置满足所要求格网大小的一定倍数即可快速判断。考虑到与边界相交的格网中心点在边界内或外的问题是否满足数据质检规则，此处建议按照比目标区行政边界稍大的矩形区域或对目标区行政边界生成一定距离的缓冲区范围所有计算点，进行整体计算与数据处理，最后剪切到目标区范围所需的上传至普查系统平台的矢量文件。在生成标准格网的同时输出格网中心点坐标文件，坐标经纬度小数点后建议不小于10 位，便于后续处理满足质检不少于8 位且不四舍五入的数据要求。主要使用的matlab 代码如下。

NDP = shaperead(‘矢量边界shp 文件’);

BoundaryLon=[NDP.X];BoundaryLat=[NDP.Y];

inpolygon(Longitude,Latitude,BoundaryLon,Bound aryLat);

P(i).Geometry = 'Point';

P(i).X = Longitude;

P(i).Y = Latitude;

P(i).xcenter = Longitude;

P(i).ycenter = Latitude;

shapewrite(P,‘输出shp 文件’).

2 地震危险性计算与数据处理

2.1 地震危险性概率计算

利用ArcGIS 软件添加格网中心点坐标或程序代码直接输出的计算点坐标数据，使用《中国地震动参数区划图：GB18306—2015》的潜源等数据资料[5]，进行4 个超越概率（50年超越概率63%、10%、2%和年超越概率10-4）的基岩地震动峰值加速度的计算。由于计算点数量较大，需考虑计算效率与计算机存储问题，可从几个方面考虑：①由于较远的潜源对基岩峰值加速度的贡献很小或可忽略不计，故不要使用过多潜源参与计算，一般选择目标区周边300～400 km 范围内的潜源即可满足要求；②尽量使用较高配置的计算机，尤其是硬盘与内存空间充足，能满足大量数据的一次性输出与存储需求；③若计算点过多，超过百万级别，为了保持计算与输出数据的稳定性，可考虑分批计算，最后对计算结果数据进行合并。计算完成后，使用简短的Matlab 或Python 程序代码提取出不同超越概率水平的基岩峰值加速度等数据，建议输出为逗号分隔符的csv 格式，便于后续生成矢量数据文件。主要使用的matlab 代码如下（以50年63%的峰值加速度为例）。

fidin = fopen('反应谱数据文件','r'); i=0;

while ～feof(fidin)

i=i+1;line1 = fgetl(fidin1);

PGA50a63=str2num(line1(12:20));

fprintf(fidout,'%10d,%8.3f ',[i PGA50a63]);

end

2.2 场地调整与地震危险性分级

根据地震危险性图编制技术规范（FXPC/DZ P-01），地震危险性概率分析计算得到基岩峰值加速度，对应Ⅰ1类场地峰值加速度，使用全国1∶100 万宏观场地类别，按照表1 所给值采取分段线性插值确定调整系数，进行调整转换为相应场地的峰值加速度值[3]。

表1 场地地震动峰值加速度调整系数Table 1 Adjustment coefficient of peak acceleration of site ground motion

根据年超越概率10-4的地震动峰值加速度（ax），将场地地震危险性分为4 级，其中1 级的危险程度最高，4 级的危险程度最低。分级原则为：1 级为ax≥760 gal、2 级为380 gal≤ax＜760 gal、3 级为190 gal≤ax＜380 gal、4 级为ax＜190 gal。值得注意的是，此次地震灾害风险普查全国30″格网的地震危险性结果，分级原则在此基础上有所调整。

地震动场地调整大致有以下3 种处理方法。

1）若完全基于ArcGIS 软件平台实现操作，则首先将基岩地震危险性计算结果转换为矢量数据shp 点文件[6]，或直接存储在gdb 数据库中会提高运行速度。新建场地类别与各个超越概率场地峰值加速度、地震危险性等级等字段，使用空间连接工具赋值场地类别，编写Python 代码实现根据场地类别和线性插值获取调整系数的函数，运用字段计算器对各个概率场地加速度值进行计算，再根据年超越概率10-4的地震动峰值加速度对地震危险性分级。当计算数据量达百万级别时，在矢量数据shp 文件中处理字段有时会不易操作且耗时较长，后续对加速度值进行归档用于风险评估与区划的输入还需进行一系列处理。

2）完全使用Matlab 程序来处理，读取宏观各个场地类别数据的边界，使用函数inpolygon 逐点判断场地类别，编写子函数线性插值获得场地调整系数对峰值加速度进行调整、并判断地震危险性等级。同时为了满足后续绘图分析以及地震风险评估所需地震危险性输入，可根据场地峰值加速度按照一定间隔赋值等值线值和进行烈度归档，输出所需数据文件，利用函数shapewrite 可生成矢量数据shp 文件。

3）由于不同场地类别都是由多个面对象组成，若在程序代码架构方面不是特别完善，特别是处理沿海岛屿较多的地区，循环中嵌套循环，在判别场地类别存在一定的困难，且耗时较长。此时使用ArcGIS 软件自带的空间连接工具赋值该字段属性效率更高，再提取场地类别并与峰值加速度数据进行合并处理，后续使用程序处理数据并生成矢量数据shp 文件。

2.3 矢量文件

按照地震灾害风险普查工作的要求，地震危险性分析部分需要生成3 个矢量数据shp 文件即地震动PGA 图（矢量点）、地震危险性等级图（矢量点）、地震危险性等级图（矢量面），在全国综合风险普查系统中通过质检才可上传成功，且各shp 文件对字段与顺序都有要求。

若前述数据操作均基于ArcGIS 软件，则可直接选择相应的字段数据输出，低版本（如ArcGIS 10.2）在选取字段时不支持调整顺序，可在其他软件（如Access）中调整表字段顺序，或者直接使用ArcGIS10.6 及以上的较高版本。若使用Matlab 程序代码，利用函数shapewrite 则可快速生成所要求的数据；带有后缀为shx、shp、dbf 三个文件是矢量数据必不可少的，加上后缀为prj 的投影文件即可，投影文件实际上是和其他数据同名的文件，程序操作时将所需投影文件内容写入即可（如CGCS2000）[6]，故在程序中只需要输出这4 个文件；输出数据字段名可自行设置，建议为不超过8 位字符的英文名且符合shp 数据字段命名规则。需要注意的是矢量数据shp 文件存储数据的上限为2 GB，超过则无法保存，需要存储在gdb 数据库中操作。为了顺利生成矢量文件，需预先构建所需字段的结构体，否则当数据量超过十万级别时生成shp 数据较为困难，所使用的matlab 代码示例如下（其中Geometry、X、Y 为shp 文件必须，longitude、latitude、risk 为输出字段，可根据需要增减）。

STR_P='struct(''Geometry'',values,''X'',values,''Y'',v alues,''longitude'',values,''latitude'',values,''risk'',values)';

values = cell(Grid_Points,1); P = eval(STR_P).

3 批量绘图与生成文档

3.1 批量绘图

根据地震危险性分析结果，需要基于ArcGIS软件平台编制峰值加速度和地震危险性等级等图件且制图需满足地震灾害普查相关标准的要求。若要求绘制某省各个区县相应的图件，则可按照省级行政边界整体进行地震危险性分析与计算处理，在ArcMap 中按照各区县范围实现自动出图。具体操作步骤如下：加载地震危险性分析结果、行政边界、水系、公路、铁路等矢量图层，启动数据驱动并进行设置，以区县矢量数据为索引图层；勾选区县名称字段，在数据框中选择剪裁至当前数据驱动页面范围，此时可以根据制图需求排除不需剪裁的图层，切换到布局视图；设置动态文本标题，插入边框、风险普查logo、指北针、比例尺等，调整到预期出图效果，保存该地图文档mxd 文件；最后选择pdf 格式批量导出地图。4 个不同超越概率水平的基岩和场地地震危险性结果以及地震危险性等级单独保存为不同的地图文档mxd 文件，逐个批量导出图件。

为提高批量出图效率需要注意以下几点：①不同的地图文档mxd 中需设置统一的图例；②水系、道路等数据量较大的图层，剪裁至比目标区稍大范围即可；③制图所用矢量图层均导入gdb 数据库中会提高效率；④区县数量较多时，在ArcMap 中批量导出较高分辨率的图件，耗时较长且可能无法导出，尤其是低版本的ArcGIS，此时可在ArcGIS Pro 中导入mxd 文件再批量导出地图；⑤若对ArcGIS Pro 较为熟悉，可按照上述类似步骤直接绘图并批量导出，较为高效；⑥设置地图导出至一个pdf 文件后，可批量转换为图片，输出顺序与数据驱动索引图层字段内容一致，此时可用简洁的代码新建以各区县名称命名的文件夹，将对应的地图图件重命名并导入。

3.2 批量生成word 文档

对于地震灾害风险普查工作而言，同一地区不同区县的地震危险性分析报告，内容上相对来说较为固定，主要包括任务概述、地震危险性分析原理、潜源划分、地震活动性参数的确定、衰减关系、地震危险性计算与场地调整等几个方面。可以设置一个通用的word 文档模板，利用Matlab 或Python等程序，构建代码启动Word 调用功能（Word.Application），批量修改报告名称、替换区县名称、增加数据统计表格等内容。同时调用上述各区县文件夹中的地震危险性图件，设置剪裁图件空白部分的参数、调整大小等处理后，插入相应的地震危险性word 报告中去，并添加各图件标题。运行程序之前，需在Microsoft Word 中将当前文件夹设置为受信任位置，可避免反复出现信任提醒而导致中断循环的操作。主要使用的matlab 代码如下。

copyfile('word 模板文件名','所需word 文件名');

try

Word = actxGetRunningServer('Word.Application');

catch

Word = actxserver('Word.Application');

end

Word.Visible = 1;

Selection = Word.Selection;

Selection.InlineShapes.AddPicture('图片路径+名称','True','True').

4 示例分析

以湖南省地震危险性计算与数据处理工作为例，需要生成上传至风险普查系统中的3 个地震危险性矢量数据shp 文件，以及湖南省14 个地市州与122 个区县地震危险性图件与报告，下面从计算过程与处理平台、效率等方面进行简要分析。采用30″标准格网（约28 万个计算点），计算出50年超越概率63%、10%、2%以及年超越概率10-4基岩峰值加速度并进行场地调整；批量绘制4 个超越概率的场地地震危险性图和地震危险性等级图（每个市或区县均5 张图件、共计680 张图），批量生成每个市或区县的地震危险性分析word 报告；所使用的计算机配置为Windows10 操作系统、内存容量32 GB、硬盘容量2 TB、inteli7 处理器。

表2 为湖南省地震危险性分析工作的具体操作步骤、处理平台与处理效率，其中地震危险性分析计算使用SEC R2019 软件，其他数据处理与分析步骤主要使用Matlab、Python、ArcGIS 等软件。可以看出，利用程序编写代码处理数据效率很高、灵活性强且易于操作，生成矢量数据shp 文件也比在ArcGIS 软件中加载数据再导出更加高效；利用ArcGIS 数据驱动批量出图对这种可重复性工作值得推荐；最后利用程序代码实现批量生成14 个市和122 个区县的word 报告，并将共计680 张图件批量插入相应的word 中相应位置并增加图名，极其高效地完成了工作。

表2 湖南省地震危险性分析与数据处理一览表Table 2 List of seismic risk analysis and data processing in Hunan province

5 讨论与结论

本文没有涉及地震危险性概率分析方法原理及相关内容，但在实际工作中需要深入分析，比如潜源参数、地震带参数取值[7]、衰减关系适用性、不同超越概率峰值加速度间比例关系等。在大量数据计算和处理之前，可先在目标区内选取一定数量的控制点进行地震危险性计算，并对不同超越概率水平的峰值加速度与潜源贡献等结果进行分析，综合判定计算结果的合理性。若需要计算反应谱，可根据数据计算量情况对不同周期点单独或合并计算。部分地区还需要特殊处理，如湖南省南部跨越长江中游地震带与华南沿海地震带，分别对应中强地震区与东部活跃区的衰减关系[8]，可分别计算并取峰值加速度的较大值做为最终结果。

地震危险性等级图（矢量面）可根据计算结果，使用ArcGIS 软件的克里金插值和等值线工具结合实现，也结合其他专业软件（如Surfer）提供的不同算法生成等值线。对于危险性等级不复杂的情况，使用ArcGIS 软件的剪切面工具手工操作，可自动满足拓扑检查的要求，也可结合ArcGIS 二次开发生成矢量面，此问题需要进一步研究解决。

需要具备一定的编程基础才能解决可能出现的各种问题并高效地完成工作。Matlab 和Python程序在数据处理领域便于使用，且ArcGIS 软件内嵌Python 程序易于二次开发。如对于湖南省地震危险性分析工作，需要在全省数据基础上得到14 个地市州与122 个区县的地震动PGA 和地震危险性等级的矢量数据shp 文件，可以使用python 代码调用剪裁工具arcpy.Clip_analysis ('输入shp','剪裁面shp','输出shp')循环实现[9]。另外，大量绘图工作均基于ArcGIS 软件平台，当危险性矢量点数据量很大时操作不便，可将点转为栅格进行处理。

地震灾害风险普查工作中的危险性概率计算数据量往往很大，为了避免存储空间不足、内存溢出等问题，使用高配置的计算机是十分必要的，在计算与后续数据处理中均需考虑是否分批进行。笔者在前述配置计算机上做过测试，单次成功计算并处理约500 万个计算点的数据，且利用Matlab 程序代码生成此级别的矢量点数据shp 文件只需几分钟，更多的数据只要内存不溢出均可处理。需要注意的是单个矢量数据shp 文件的大小上限为2 GB[10]，超出则需要在gdb 数据库中存储和处理。另外，在较粗格网地震危险性计算结果基础上，使用克里金插值等方式获取较细格网数据，可大大减少地震危险性计算的时间，但需要对结果的合理性进行论证。

文中所述数据处理、批量出图与生成word 文档等方法，在地震安全性评价工作也可参考。如线状工程的沿线区划工作，可参考标准格网建立所使用的方法，外加线状工程一定距离缓冲区加以约束，可快速生成区划范围所需间隔的计算点。在区域性地震安全性评价工作中，大量的地震动时程与规准谱标定图件，亦可利用程序构建简洁的代码调用Word.Application、设置图片大小等参数批量快速的导入word 报告。

文中所述来源于笔者实际工作中的探索，并写了相应的Matlab 和Python 程序代码，由于水平有限可能存在不足，在后续工作中与感兴趣科研工作者交流分享。

致谢 非常感谢地震灾害风险普查专题培训工作中专家们的悉心指导及工作中同事们的帮助，也感谢审稿专家提出的宝贵建议。