基于Mapping Toolbox的地震图件产出

欧阳小波，汤曜伟

（广东省地震局，广东 广州 510000）

0 引言

Matlab 是由MathWorks 公司推出的一款强大的数值计算软件[1]，是数值计算、符号运算和图形处理等多种功能的强有力实现工具，近几年来Matlab 这个强大的科学计算软件包已经得到业界的广泛认可，并已深入到了各行各业的众多学科、在各大公司、科研机构、大学校园得到了广泛应用。其具有以下优点：（1）代码简洁直观，以矩阵为运算基础，内含有丰富的函数库，在计算方面避免使用者花费大量时间去编写不必要的小程序；（2）数值计算方面，科学计算功能强大，涵盖了从解微分方程，数字信号处理FFT，到最新的最优化技术诸如神经网络、遗传算法等；（3）数据可视化方面，Matlab提供了多种基本图件处理函数，如plot、scatter、contour 等，可以很容易的绘制出高质量的工程图谱；（4）Matlab 提供了强大的编程接口，支持Matlab 与其他编程语言进行数据交换。在地理信息数据可视化应用方面，Matlab 推出了专用于地理信息图件绘制的工具箱Mapping toolbox，在Matlab 官方的Mapping toolbox user guide 中是这样介绍Mapping toolbox 的，“地图工具箱包含了用于在Matlab 环境下创建地图、显示、分析或操作地理空间数据的一系列扩展函数库和图形用户接口（GUI），用户能够轻松的综合多种类型的不同地理数据源并根据它们的空间位置关系完美的展示出来。总的来说，Mapping toolbox 函数库提供了以下主要功能：（1）地理空间数据的导入和接入；（2）地理图像与数据的网格化处理；（3）地图在不同地理坐标的投射；（4）地图的展示和互动函数（如抓取要素、点击位置）；（5）支持向量及栅格数据的地理信息计算等。由此可见Mapping toolbox 在与地理信息有关的数据的可视化功能方面十分的强大。

在防震减灾工作领域，无论是地震监测预报还是灾害防御、应急救援。数据的可视化工作有着不可替代的作用，在地震监测预报科研方面，震中分布图、震源深度分布图、结合研究地区的地质构造状况，地震研究人员可以分析与研究地震与对应地区位置地质构造之间的关系，地震热点地区地震的时空迁移变化趋势，探查并了解其中蕴含的规律，为今后的地震预测研究积累经验以及区域震害防御提供指导意见。而在灾害防御、应急救援方面，地震烈度区划图又为抗震设防提供重要依据。鉴于Matlab 自带的强大的数据科学计算能力以及Mapping toolbox 提供的多种地图绘制功能，本文尝试以Mapping toolbox 作为出图工具，针对几种常见地震图件，根据其产出规范要求，利用Matlab 绘制图件，并给出对应实现代码。作为Matlab 这种强大的第四代计算机语言在地震图件产出中的研究探索。

1 Mapping toolbox 绘制地震图件流程

从地震图件的角度来看，有文献将地震图件大体分成了两类，第一类是反映某个或者某几个地震专题要素的空间分布状态，例如历史地震分布图、台站分布图、重点目标及危险源分布图等。第二类是以地震专题统计数据为基础的统计地图，用于反映专题信息随着地理位置变化而分布的情况。第一类地震图件一般采用散点图的表现形式，第二类图件一般采用等高线、分级色彩图等表现形式[2]。对于上面的二类地震图件，不管如何变化，在Maltab 的Mapping toolbox 绘制过程中，均需要遵循以下几个流程来绘制图件。步骤1：基础绘图数据导入；步骤2：地理底图的导入；步骤3：使用各种绘图函数（点、线、面）等绘制图形；步骤4：图形其他要素添加以及标注等。下面就这两类图件按照这四个步骤绘图中出现的技术难点加以介绍。

1.1 地理底图导入

Mapping toolbox 对目前主流的地理信息数据诸如shapefile、GeoTIFF 和KML 都有着对应的数据导入和写函数，另外，在手头资源不充足的情况下，Mapping toolbox 也提供了网络地图服务器数据导入功能。本文以美国环境系统研究所公司（ESRI）开发的一种空间数据开放格式shapefile 为例，说明Mapping toolbox 地理底图导入操作，shaperead 函数的基本操作是Map=shaperead（filename），直接读入名为filename 的shapefile 文件，返回值Map 是一个包含大量数据的地理数据结构体。shaperead 函数第二种读入形式为Map=shaperead（filename，Name，Value，...）Names 是地理数据结构的其他属性部分，Value 是Name 的值，通过指定Name 和Value 值，可以有选择性导入某区域地图或某个地理要素到matlab 运算空间。然后再通过函数mapshow（Map） 就可以将导入到空间中的地理底图显示出来。

1.2 基础绘图数据导入转换

前述说过，Matlab 的Mapping toolbox 能够轻松的综合多种类型的不同地理数据源并根据它们的空间位置关系完美的展示出来，在Mapping toolbox 中，向量数据与栅格数据是其主要处理对象，对于非地理行业内规范的地理空间数据源如我们地震行业的地震编目数据，按Mapping toolbox 中对于地理空间数据的定义来说，地震编目数据可以称之为一种地理空间数据，对如不同的非规范化的数据源，Mapping toolbox 提供了一种地理数据结构体的结构来统一创建和显示非规范的数据，这也是一种向量数据，在这种Mapping toolbox 规定的地理数据结构体中，有几个成员要素是必须的，第一个是表征几何形状的geometry，其值分别是点Point，多点MultiPoint，线Line 或者多边形Polygon。第二个是经纬Lat 和Lon，最后是补充属性特征如地名、街道名等，对应编目数据的属性就是震级大小。

1.3 重要的绘制函数

对于第一类地震图件，使用mapshow 与geoshow 函数能够很好的完成创建好的地理数据结构体的绘制，matlab 也提供了plotm、linem等函数完成点、线的绘制，对于等高线可以使用contourm 函数来完成该类线条的绘制工作。使用contourm 函数需要的是网格化数据，对于给定的经纬度范围，需要用到网格创建函数meshgrid，其调用格式如下：［X，Y］=meshgrid（x，y）其涵义为通过数据重复在一维数组X，Y 的每一个交叉点上创建网格点，当X 和Y 都是长度为n 的一维数组时，产生的X，Y 就是x，y 网格化之后的二维矩阵，在实际的地震烈度图件绘制时候，输入与x，y 对应的经纬度数组序列，即可以产生经纬度网格平面。另外在生成网格之后，如果生成的网格较密，那么就需要对原始绘图数据进行插值处理，填补原始数据是空白的网格点，特别是对于第二类地震图件，再利用mesh，surf，contour 等绘图时，当利用meshgrid 函数完成经纬度网格化之后，由于观测数据是散乱分布的，想要获得平滑的数据插值，就需要进行插值运算。其调用格式为：Z=girddata（x，y，z，X，Y，method） x，y，z 是原始地震数据点，X，Y 是由x，Y 利用meshgrid 网格化得到的网格坐标，method 指定了插值方法，Matlab 提供了linear（线性），cubic（立方体），nearest（最近邻点）和v4（Matlab4 网格点）四种插值算法，其中返回值Z 是插值结果，用cubic 与v4 得到的插值曲面连续光滑，而linear 和nearest 则不连续，默认情况下Matlab 会采用linear 线性插值。

1.4 图件的修饰

Matlab 提供了各类图件修饰的函数legend、colorbar 等来给图件进行要素标识，在震中分布图中可以用legend 函数标识不同震级大小的地震。clabel、colorbar 颜色渐变条可以用于等震线的绘制。

2 实例绘制

2.1 绘制震中分布图

地震震中分布图是地震图件中一种重要的基础地震图件，用于绘制地震震中分布图的基础数据来源主要是各省级台网产出的地震目录。目前其最新的行业标准规范是中华人民共和国地震行业标准（DB/T66—2016），已于2017 年05 年01 日开始实施，在其规范性附录中对CSF 地震目录和地震观测报告交换格式有明确详细的范例[3]。如图1 所示。

图1 CSF 地震目录标准格式范例Fig.1 The example of CSF seismic catalog standard format

按照前述的绘图流程，根据地震目录格式规范和Mapping toolbox 的要求，我们要进行绘图工作，首先要将规范的地震目录格式转换成符合Mapping toolbox 要求的地理数据结构体。地震震中分布图都是采用散点图的形式，其结构体geometry 值应选择point，建立一个SeismicInformation 的结构体存储地震目录中的关键数据，现将代码展示如下：

［a，b］=size（catalog）；%导入之后目录长度

［SeismicInformation（1：a）.Geometry］ = deal（′Point′）；%创建了一个SeismicInformation 结构体，其Geometry 成员值为点Point

［SeismicInformation（1：a）.lat］=catalog（a，（27：33））；% 导入目录中的经纬度

［SeismicInformation（1：a）.lon］=catalog（a，（35：42））；

［SeismicInformation（1：a）.Mag］=catalog（a，（47：50）；%结构体补充属性震级

这样一个完整的用于绘制震中分布图的地理数据结构体就创建完成。

然后再导入地理底图，在这里使用shaperead函数导入整个中国区域地理底图，并制定显示区域（此处显示广东省）。

ChinaMap=shaperead（′china.shp′）；%读入整个中国地图

mapshow（ChinaMap）；%软件初始化显示整个中国地图

axis off 消除坐标轴背景

size=length（ChinaMap）；%q 确定中国地图结构体数组长度

Province=’guangdong’；

cla；

for i=1：size %搜素整个地图结构体数组

if strcmp（ChinaMap（i）. Name，Province）==1

mapshow（ChinaMap（i））；%仅显示整个中国地图中为guangdong 省的区域

hold on；

end

end

再按震级大小分别显示SeismicInformation结构体

for i=1：length（SeismicInformation）

if SeismicInformation（i）. Mag＜=1

mapshow（SeismicInformation，′Marker′，′o′，′MarkerSize′，4，′MarkerFaceColor′，′c′，′MarkerEdgeColor′，′k′）；

........

end

end

本文利用广东台网2018 年6 月至12 月数据绘制震中分布图如图2 所示。

图2 地震震中分布图例子Fig.2 The example of epicenter distribution

2.2 绘制地震烈度图

地震烈度图是描述地震引起地面震动及其影响程度的图件，一般用等震线描述，等震线又是指不同地震烈度和地震震动强度的分界线，地震烈度的具体数据一般按照地震烈度表由现场宏观考察得到或者由强震计地面峰值加速度值得到，由于笔者未有相关地震烈度数据来源，利用地震发生之后各台站地震计速度峰值数据代替作图，为显示图件效果。前述地理底图的代码已展示，此处主要展示等值线代码。

lat= Stationlocatio（：，1）；%台站纬度

longt= Stationlocatio（：，2）；%台站经度

Value=Stationlocatio（：，3）；%地震峰值加速度值

xi=linspace（min（lat），max（lat），50）；

yi=linspace（min（longt），max（longt），50）；

［xi，yi］=meshgrid（xi，yi）；%按50 个点网格化所有经纬度

zi=griddata（lat，longt，Value，xi，yi，′linear′）；%网格化之后线性插值

［c，h］=contour（yi，xi，zi）；%绘制等震线

Clabel；%在等震线上标注数值

Colorbar；%加颜色渐变条

作图其效果如图3 所示。

图3 地震烈度等震线示例图Fig.3 The isoseismal line of the seismic intensity

3 结语

本文通过以上的两类不同类型地震图件的专题绘制研究，阐述了利用Matlab 绘制图件的四个基本步骤，展示了Matlab 的Mapping toobox强大绘图功能以及Matlab 极为简洁的编程语言,经过实例验证，Matlab 能够绘制出满足要求的具有较好表现效果的地震图件，更好的为防震减灾工作服务。最后，需要说明的是在日常的地震图件绘图中，地震图件的质量很大程度上依赖于自备的地理底图库，即地理信息库的丰满程度，因此对于绘图者来说，建立和保存自身常用的制图底图模板非常重要。