0 引言

地震发生后，政府和地震部门需要快速了解地震影响范围和可能造成的灾害损失，以正确部署应急响应工作。在救灾人员到达现场前的盲场快速评估阶段，基于地震三要素，利用震前构建的数据库和模型预估地震影响范围和灾害损失程度，是下一步抗震救灾工作的重要依据。

在应急指挥技术系统中，对地震进行震害评估的第1步是快速绘制地震影响场。目前，灾情快速评估系统通常采用距震中最近活动断层方位作为影响场的长轴方向，全国分几个大区构建地震烈度衰减模型，基于震级生成地震影响场。然而，在某些震例中，影响场方向可能与极震区长轴方位相差较大，如云南鲁甸6.5级、新疆精河6.6级地震，这导致极震区得不到快速有效的救援。因此，迅速而准确地绘制地震影响场的方向，是科学有效地进行地震应急指挥决策部署的前提条件。

地震影响场的构建和共享主要基于GIS平台（周文等，2012；孙哲等，2018；陈雅慧等，2018），常用的有C、Java等编程语言，ArcMap、MapInfo等软件。这些软件运行速度快，性能稳定，但价格不菲，开发难度大，周期长，一般委托专业公司进行开发。然而，地震影响场专业性较强，模型较多，与断层、强震等紧密联系，随着研究的深入需时常调整算法，相应的程序代码也要更新，这使委托开发具有一定的难度。目前，有些省地震局已开始尝试自主研发地震影响场并对其进行技术改进（李铂等，2012；周文等，2012；杨彦明等，2017；孙哲等，2018）。以Python为例的开源软件可免费使用，编程速度快，工作量小，运行速度基本满足需求，而且技术算法的更新方便快捷。Python语言以“多面手”而闻名，它可以是面向对象编程语言，也可以是过程式脚本语言，甚至还可以是函数式编程语言；其原生数据类型是对象，并且所有的Python库都是模块化的，遵循了基本对象的结构和行为（Lawhead，2015）。Python标准库和网络库提供很多功能强大的模块用于地理空间分析，如OGR通用矢量库、Shapely通用几何库、Fiona地理处理包等。本文基于Python语言实现V度及以上地震影响场自动生成，并剔除海域和境外部分。此外，尝试使用早期余震分布结合震源机制解评估极震区长轴方位的偏移，以实现地震影响场的多时段动态修正。

1 地震烈度衰减模型

地震影响场的绘制是一项多学科交叉的基础性研究工作，也是实际应用的过程。首先，要根据每个区域的地震地质特点和历史地震等震线资料，建立适合该区域的地震烈度衰减关系；然后，采取一定方法判别影响场长轴方位；最后，以震中位置为中心点，基于衰减关系和震级绘制地震影响场。我国地震烈度衰减关系的统计分析中，通常假设地震震源为点源模型，地震烈度衰减采用椭圆模型，衰减公式为

$ I = a + bM + C\lg (R + {R_0}) $

(1)

其中，I为烈度；M为震级；R为震中距（单位：km）；a、b、C、R₀为回归系数，长、短轴各有一套不同的系数。地震烈度随震中距的增大而减小，而地震烈度衰减关系具有较强的地域性，各个地区地壳介质的吸收、几何扩散及散射情况的不同，使得地震波的能量表现出不同的规律性衰减。

此外，大地震时，地表破裂尺度的存在会导致经典的点源模型与真实情况有出入，如汶川8.0级地震的等震线并非椭圆，烈度衰减呈现出长矩形衰减模式（李志强等，2008）。因此，在处理强震时需考虑地表破裂尺度的影响再对地震影响场进行修正，使其逼近实际等震线（李铂等，2012）。

“十五”期间建设完成的地震应急指挥技术系统中，采用距震中最近活动断层方向作为影响场的方向，全国分几个大区构建地震烈度衰减模型，将地震参数输入模型生成六度及以上地震影响场，且未剔除海域和国境线以外的部分。本研究中，收集整理了全国地震烈度衰减关系，并进行归纳和细化，发震后根据震中所在区域读取相应的衰减关系，该方法较使用统一衰减模型方法适用性更强。

2 基于余震、震源机制解动态修正影响场

地震影响场作为灾情评估和辅助决策的依据，其正确与否可能会影响应急指挥调度和决策部署。然而，在某些震例中，影响场快速评估结果与实际等震线调查结果有出入，可能导致极震区得不到快速有效的救援。因此，需充分考虑震源本身的特性，参考震源机制解、余震序列分布等对影响场进行动态修正。

震源机制解（地震断层面解）反映了震源本身的特性，可通过波形分析获得矩张量并以震源球的形式将其展现出来。震源机制解包含断层走向、倾角和滑动角等3个参数，1次地震通常给出2组节面。对于走滑型断层引起的地震，震源机制解的2组节面约呈90°夹角，可利用余震分布并结合震源所在位置的地质构造特征综合判定影响场方向。对于倾滑型为主的地震，震源机制解给出的2个节面方向几乎位于1条直线上，可直接利用其进行影响场长轴方向的修正。此外，有研究者认为震源机制解对地震烈度衰减关系有显著的影响，可利用震源机制解对地震影响场的长、短轴进行初步修正（郑韵等，2016）。

大地震后可能伴随着一系列余震，Kisslinger等（1991）、Richter（1995）认为主震后的早期余震一般发生在主破裂面上，蒋海昆等（2007）认为余震分布尺度粗略地与主震破裂尺度相当。因此，对于余震丰富的地震，可以利用余震序列对地震影响场进行动态修正，以使其逼近真实烈度分布。国内研究者已开展了一些基于余震评估宏观震中和极震区的研究，可根据地震特征进行影响场动态修正，如利用余震序列（白仙富等，2011；王伟锞等，2011；杨天青等，2015）、灾情信息（杨天青等，2016）等资料。本研究中，通过读取震源机制解和余震分布的优势方位，动态修正了影响场长轴方向。

3 程序编制与实现 3.1 总体设计思路

2004年开始运行的地震速报信息共享服务系统（Earthquake Instant Messenger，EQIM）可在震后快速提供地震速报信息（杨陈等，2009）。读取EQIM消息，获取震中经纬度、震级和震源深度等信息，将其作为绘制地震影响场的输入参数。调用数据库中的地震烈度衰减经验公式和活动断裂分布情况，计算衰减椭圆方向及各烈度区的长、短半轴，可快速绘制地震影响场。

随着时间的推移，获取地震的震源机制解和余震序列分布情况，读取它们的方向可以进行影响场长轴方向的动态修正。程序主要包括以下几方面功能：①读取速报信息；②生成震中文件和地震影响场文件；③输出消息；④动态修正并输出消息。程序流程如图 1所示。

3.2 Python模块介绍

本研究中使用的Python模块主要有：fiona、shapely.geometry、math、logging、os、time等。Fiona采用Python中内置的数据结构表示矢量数据，要素以GeoJSON格式表示，使用Python内置的字典（dict）结构组织；1个图层包含在1个集合中（Collection），如：{'type': 'Point', 'coordinates': (100.07, 25.89)}。使用时，可以对该集合进行遍历，得到所需要素。shapely.geometry用于操作和分析几何对象，主要使用到其中几个函数，shape()用于返回1个与对象坐标相同的独立几何结构；mapping()将对象转为GeoJSON格式，为消息传递做准备；LineString()和Polygon()分别生成线和面对象。此外，该程序还使用了其他模块，math模块用于数学运算；logging标准库用于日志记录；os模块用来处理文件和目录；time用于返回时间戳和显示运行时间。

3.3 程序设计实现

（1）设置工作目录。为每个真实地震或测试地震新建1个文件夹，以地震ID命名，用于存放生成的震中文件和影响场文件。代码示例

output_work_path = r"f:\202002result"

output_work_path = os.path.join(output_work_path, str(dict_eqim['Event_id']))

if not os.path.exists(output_work_path):

    os.mkdir(output_work_path)

若同一地震生成不同批次的影响场，则在ID后附加“01”、“02”等。此外，为便于查阅，增加日志功能和运算时间显示功能。

（2）生成震中文件。接收地震消息后，首先要判断震中是否位于国内，如果在国内，则提取地点信息，并生成地震震中文件。震中文件一般采用WGS-84坐标系（EPSG:4326），属性表包含地震ID、发震时间、震源深度、震级、地点和经纬度。示例

schema = {'geometry': 'Point',

       'properties': {'ID': 'str', 'TIME': 'str', 'DEPTH': 'float', 'GRADE_MS_': 'float',

                  'LOCATION': 'str', 'POINT_X': 'float', 'POINT_Y': 'float'}}

epi_layer = fiona.open(epicenter_filename, mode='w', driver='ESRI Shapefile',

                   schema=schema, crs='EPSG:4326', encoding='utf-8')

feature = {'geometry': mapping(point),

'properties': {'ID': dict_eqim['Event_id'], 'TIME': dict_eqim['Event_id'], 'DEPTH': dict_eqim['depth'], 'GRADE_MS_': dict_eqim['magnitude'], 'LOCATION': dict_eqim['location'], 'POINT_X': dict_eqim['lon'], 'POINT_Y': dict_eqim['lat']}}

epi_layer.write(feature)

epi_layer.close()

（3）生成影响场。根据震中经纬度判定震中所在区域，读取相应的衰减关系和活动断层情况。以震中为椭圆中心点，距离最近的活动断层方向为椭圆的旋转角，从Ⅴ度烈度开始生成衰减椭圆的点。由于震中经纬度的单位是度，衰减公式的单位是千米，所以先将震中经纬度进行兰勃特方位等积（Lambert azimuthal equal area）投影，然后再进行运算。为了加快运行速度，减少不必要的运算，可根据长轴判定需要多少个点绘制椭圆。如R＜10 km时仅需计算36个点；R＞150 km时则需计算180个点，以使生成的线光滑。生成的点放置于列表中，用shapely.geometry.LineString()命令生成线集合rings()。绘制时从最高烈度开始，第1个烈度椭圆即极震区，第2个烈度椭圆减掉极震区是第2个影响场，以此类推直到Ⅴ度，然后用Polygon()命令生成面。新建影响场文件，同时打开国界文件，判定影响场是否与中国国境接壤，若有延伸到海域和境外的部分，则作空间交割，最后生成的影响场只含国内部分。示例

schema = {'geometry': 'MultiPolygon', 'properties': {'Intensity': 'int', 'Area': 'float'}}

intensity_layer = fiona.open(intensity_file_name, mode='w', driver='ESRI Shapefile', schema=schema, crs='EPSG:4326', encoding='utf-8')

bound_layer = fiona.open(boundryfile, 'r', driver='ESRI Shapefile')

boundary_geom = shape(bound_layer[0]['geometry'])

for index, ring in enumerate(rings):

    calc_polygon = Polygon(ring)

    if index＞0:

        calc_polygon = Polygon(rings[index]).difference(Polygon(rings[index - 1]))

if calc_polygon.intersects(boundary_geom) is True:

      calc_polygon = calc_polygon.intersection(boundary_geom)

以交割后的面文件为投影区域，投影后求出各个烈度圈的面积并保存在属性表中，将影响场文件存储于服务器中。

（4）文件输出和消息传递。生成的震中文件和影响场文件均为shp格式，为了进行下一步的人员伤亡和灾害损失评估计算，需要把影响场shp文件转换成GeoJSON，以消息的方式输出。示例

features = [feat for feat in fiona.open(intensity_file_name, 'r', encoding='utf-8', )]

intensity_geojson = {"type": "FeatureCollection",

                        "crs": {"type": "name", "properties": {"name": "urn:ogc:def:crs:OGC:1.3:CRS84"}}, "features": features}

3.4 需要注意的问题

目前的地理信息软件库在处理经度180°附近或纬度90°附近的计算时可能发生错误。由于接收的真实地震消息包含国外的地震，例如发生于新西兰的地震，所以本研究在生成经度＞175°的地震影响场时，先将震中挪到175°进行计算，衰减椭圆计算完毕后再移回震中位置绘制地震影响场。

4 应用实例 4.1 影响场产出过程

影响场产出过程如图 2所示。地震触发后，程序启动generate_intensity主模块，设置工作目录，确定震中文件名和影响场文件名，然后启动generate_epicenter子模块和generate_CreatShp子模块，分别产出震中shp文件和影响场shp文件，最后启动CreatShp_geojson子模块，生成影响场GeoJSON文件并以消息方式输出。

4.2 影响场产出实例

以2020年3月23日新疆阿克苏地区拜城县5.0级地震和【测试】云南省怒江傈僳族自治州泸水市6.0级地震为例对程序进行测试，并基于QGIS产出影响场图件（图 3，图 4）。

5 结束语

地震发生后，地震影响场的快速绘制是灾害评估的基础。获取地震三要素后，依据震级、活动断裂分布情况和震中所在区域的地震烈度衰减关系，可快速绘制第1时段地震影响场。随着时间的推移，陆续获取震源机制解和余震序列，则可综合考虑这些结果动态修正地震影响场，产出多时段应急图件，为地震应急处置提供参考依据。现今，影响场从算法到技术都已经成熟，然而委托专业公司开发收费高，开发周期长，修改算法复杂。基于开源软件生成和修正影响场完全免费，有效缓解了经费上的压力。更重要的是，代码量并不大，地震工作者可以自行完成程序编制，并且可时常结合业务对模型和算法进行调整。该方法充分反映了相关领域的最新进展，在实际地震应急工作中具有较高的实用价值。