0 引言

地震发生后，相关部门需要及时向政府和决策部门提供震后破坏情况，地震烈度是判定地震破坏影响的重要参数。以往，为了快速判定地震影响，地震部门通常根据地震参数结合本地或者本区域地震烈度衰减关系，大致估算地震影响场，或者通过震后现场调查获取地震烈度。然而，地震烈度衰减关系大多根据历史震例得出，并不能反映某次地震的真实情况；现场调查则可能需要几小时到几天，严重影响震后应急的时效性。后来，随着各地强震台站的建设，又有研究者通过利用插值网格点周围的几个台站记录，结合地震烈度衰减关系得到网格点的烈度值（张红才，2008；金星等，2010）。以上均为在测震台网密度较小情况下采取的措施。随着国家“十五”项目及国家地震烈度速报与预警工程项目的开展，各地布设的地震监测站网密度不断提高。随着地震预警项目的开展，北京已建成199个地震监测台站，且台间距小于10 km。利用这种高密度地震台网，能在震后比较真实地反映地震动情况，也为直接利用仪器烈度值生成烈度等值线提供了可能（王士成等，2017）。目前，北京市地震局地震烈度速报系统已能自动生成各个仪器点位的烈度值（白立新等，2019）。然而，这些仪器站点在空间上呈离散分布状态，仅记录各站点的烈度实测值，并未产出传统所需烈度影响场。

地震发生后，可从地震台站直接获取地震动数据记录，利用烈度等值线直观了解地震对整个区域的影响。根据地震烈度影响场，可以在后续应急和快速评估中给出不同破坏程度的影响范围，进而为政府决策提供可靠依据。因此，迫切需要将地震烈度速报系统产出的离散烈度点生成线状或面状的影响范围。

将离散点生成线状或者面状的等值线（面）图，需要在现有实测值基础上进行空间插值。空间插值分析算法是一种将离散点的测量数据转换成连续表面的算法。地统计插值是基于空间自相关性的，由观测数据产生具有统计关系的曲面，代表插值法是Kriging（克里金）插值法（李海涛等，2019）。本文根据北京市高密度地震监测台网生成的实测烈度值，开发1套基于Arcpy的程序，利用ArcGIS平台提供的Geostatistical Analyst模块进行克里金插值，生成了系列发生在北京及周边地区地震对北京产生影响的地震烈度等值线分布图。

1 克里金插值原理

采用克里金差值法，对未被地震台站覆盖区域进行插值计算，从而生成等值线图。

克里金插值是一种广泛应用于空间插值的地学统计方法，是在空间相关范围分析基础上，用相关范围的采样点来估计待定点属性值。与其它插值方法相比，克里金法的显著特点是使误差的方差最小，因此也称最优无偏估计。借助克里金插值，可将仪器烈度散点值及其相关距离拟合出变差函数，同时在插值范围划分网格，得到等间距的网格点。通过求解克里金方程组求得各网格点权系数，从而得到各网格点的烈度估计值，将网格上估计值相等的点用平滑曲线连接，就可得到烈度等值线图（胡杨，2017），具体流程见图 1。

2 基于ArcGIS的Python二次开发

ArcGIS软件整合了GIS与数据库、软件工程、人工智能、网络技术及其他多方面的计算机主流技术，具有强大的地图和分析软件功能，是地理信息系统（GIS）软件、位置智能和地图绘制领域的全球市场领导者。ArcGIS在不同时期提供了不同的定制开发方式。使用Python语言进行ArcGIS的定制开发，是近年ESRI力推的开发方式。其基于ArcGIS平台的强大功能，使用Python语言进行改造和定制，可满足数据处理和分析需要，从而提升ArcGIS数据的处理效率，更好地实现ArcGIS内部任务自动化（Zandbergen，2014；Pimpler，2017）。

Python是一种解释型、面向对象、动态数据类型的高级程序设计语言，具有丰富和强大的库，使开发人员能够借助函数库实现某些复杂功能（埃里克·马瑟斯，2020）。ArcGIS在10.0以上版本引入一个Python站点包Arcpy，不少开发者基于该站点包，根据需要进行二次开发。如：檀斌等（2022）根据不同烈度衰减模型，利用Arcpy的managemen模块和analyst模块，快速生成地震烈度圈；薛蕾等（2019）使用Python语言编写了地震烈度速报软件，可以实现自动触发、分析计算作图和推送功能。周崴等（2020）、夏苏琼等（2022）采用Python语言，结合ArcGIS系统自带的Arcpy站点包，通过调用特定模块、函数和类，利用多种数据源，实现了专题图的自动化批量制作。

3 基于Arcpy生成地震烈度等值线

北京地区地震监测台站台间距虽然接近10 km，但是基于现有实测仪器烈度值进行烈度等值线绘制仍有一定差距，有必要在现有结果基础上进行插值，以获得理想的等值线图。

3.1 数据准备

北京地震台现有1套地震烈度速报计算系统，在震后根据实时地震记录，能够生成各台站实测仪器烈度值及强震报告等。各台站仪器烈度值包含在一个intensity.txt文件中，该.txt文件包含台站名、台站经纬度信息及震后记录的峰值速度、峰值加速度和烈度值等。示例信息如下：

3.2 Arcpy Geostatistical Analyst模块

使用Arcpy的Geostatistical Analyst模块进行台站仪器烈度值插值。该模块可提供一组功能全面的工具，以创建并显示、分析和了解空间现象的表面。其通过存储于点要素图层或栅格图层或多边形质心的测量值创建连续表面或地图。在使用该模块进行地统计分析时，需提供以下参数，具体见表 1。

在利用Arcpy的Geostatistical Analyst模块生成地统计图层之前，需要提供一个地统计模型源，是一个地统计图层或地统计模型（XML）。本项目模型源在ArcGIS桌面平台ArcMap利用Geostatistical Analyst扩展模块的Kriging插值方法生成。

3.3 地震烈度等值线生成流程

基于Arcpy的地震烈度等值线生成流程见图 2。

以台站仪器烈度文件intensity.txt生成的shp文件为输入数据集（in_datasets），以ArcMap生成的克里金插值图层kriging.ly为地统计模型源，通过Arcpy.ga.GACreateGeostatisticalLayer模块生成地统计图层，并利用Arcpy.ga.GALayerToContour模块将地统计图层转换成等值线，结合绘图边界shp文件、震中shp文件、天地图底图，利用Arcpy.mapping模块实现自动化制图，输出地震烈度等值线专题图。

4 主要过程程序代码 4.1 生成地统计图层

使用地统计分析模块前，检查模块许可，具体代码如下

try:

if Arcpy.CheckExtension(“GeoStats”) == “Available”:

Arcpy.CheckOutExtension(“GeoStats”) else:

# raise a custom exception

raise LicenseError

使用现有地统计图层创建新的地统计图层，具体代码如下

inLayer = self.moduleDir + r“shpData” + os.sep + “Kriging.lyr”#要分析的地统计图

inData = “intensity F1 = intensity”#要分析的地统计对象

outLayer = “intensity2” # 输出的地统计图层

Arcpy.GACreateGeostatisticalLayer_ga(inLayer, inData, outLayer) # 模型源、对象、由该工具生成的地统计图层

4.2 等值线生成

以下代码主要设置等值线的类型、平滑度、间隔等相关参数。

contour_type = “Filled_contour”#表示地统计图层的等值线类型，用面表示地统计图层。self.contourFC = self.workspacePath+os.sep + “contour”

contour_quality = “Presentation”# 确定等值线制图表达的平滑度

classification_type = “Manual”# 指定如何计算等值线间隔

classes_count = “”# 指定输出要素类的类数，如果将classification_type设置为Manual，则该参数将不适用

classes_breaks = [1, 1.49, 2.49, 3.49, 4.49, 5.49, 6.49, 7.49, 8.49, 9.49, 10.49, 11.49]

# 将classification_type设置为Manual时的中断值列表。

# Check out the ArcGIS Geostatistical Analyst extension license Arcpy.CheckOutExtension (“GeoStats”)# Execute GALayerToContour # 将地统计图层导出为等值线要素类

Arcpy.GALayerToContour_ga(outLayer, contour_type, self.contourFC, contour_quality, classification_type, classes_count, classes_breaks)

5 震例测试

收集2019—2022年京津冀范围内对北京地区影响烈度达Ⅱ度以上地震，得到8个案例（表 2），利用本程序生成地震烈度等值线（面），绘制仪器烈度等值线推测图，图中未绘出Ⅰ度区，是因为烈度为Ⅰ度时影响较小，程序不生成该烈度影响范围。文中列出北京顺义区双丰街道M_L 2.3、北京朝阳区M_L 2.7、河北唐山市古冶M_S 5.1、北京门头沟区妙峰山镇M_L 3.6地震的仪器烈度等值线图，见图 3—图 6。以上几个震例震级偏小，达不到破坏地震的程度，生成的烈度等值线数量有限，反映的内容不丰富，但是图中的Ⅱ度、Ⅲ度等值线仍较好地包络了相同烈度值台站范围，能够比较直观地反映地震的影响。

6 结束语

Python语言结合ArcGIS提供的地理信息处理工具，可以高效实现地理数据分析和地图自动化等功能。基于ArcGIS平台，利用Arcpy站点包，通过Python语言对现有仪器实测地震烈度值执行插值计算，可以快速实现烈度等值线（面）的绘制，绘制的等值线图能够较好地反映地震影响范围。