2018, Vol. 38
2. 中国地震局第二监测中心，西安市西影路 316 号，710054
在中国地震局监测预报司的支持下，基于GIS的地震分析预报系统(MapSIS)于1999年开始建设，2002年初完成研制，随后在全国地震系统广泛推广使用。MapSIS系统融合近百种分析预报方法和手段，成为全国地震系统分析预报人员重要的资料处理工具[1-2]。
MapSIS系统虽于2010年完成一次大的架构调整，但仍不能满足日益增长的用户需求和日新月异的信息技术的要求。主要表现在：1)技术体系陈旧老化；2)计算效率有待提高；3)系统臃肿不易维护；4)知识产权边界模糊。其中，知识产权问题已成为目前亟待解决的问题。原有系统依赖于MapInfo商业平台，需要用户付费安装授权的MapInfo运行环境；原有的MapInfo制图模式陈旧，用户需要更丰富的可视化效果辅助分析决策。鉴于此，中国地震局监测预报司提出，加快推进MapSIS系统的升级和改版工作，建设基于免费开源制图工具的地震分析预报系统。
1 基于GMT的MapSIS制图子系统设计GMT是功能强大的开源制图工具，遵循LGPL(GNU lesser general public license)开源协议。其自带80余种制图工具，支持30多种输出方式，可运行于Windows、Linux、Mac OS三大主流操作系统[3]，近些年在地震成果的表达中应用愈加广泛[4]。
由于GMT遵循UNIX模块化思想，以命令行方式执行制图，这给普通用户的使用带来很大难度。一个要素完整的图件的绘制，往往需要数个模块命令、几十种选项，按照一定的顺序编写脚本。此外，选项可接不同格式、不同长度的参数，这些参数对中文字符、编码格式、空格敏感，导致用户在学习和使用GMT过程中极易出错。本文采用人机交互可视化方式集成GMT工具，用户不需了解任何GMT语法即可制图；引入模板化思想，移植性高、易维护，避免个人用户修改和调试代码的繁琐工作。
为确保GMT制图工具与地震产品加工产出流程的继承性和一致性，基于可视化GMT的MapSIS制图子系统包括GMT模板制图和相应的数据处理功能，制图产出的矢量Postscript文件(以下简称PS文件)可直接解码显示，形成了从数据预处理、数据编辑、加工到产出图件一套完整的GMT制图业务流程(图 1)。
|
图 1 基于GMT制图子系统功能设计 Fig. 1 The function of the GMT mapping system |
GMT采用命令行方式绘图，其优点是制图灵活，但同时参数格式复杂、规范严格、出错率高、用户编码效率低。本系统引入模板化思想，在系统配置文件中写入统一样式的地震常规或应急产品运行模板。程序调用时，根据不同产品的固定特征，在模板库中匹配查找对应模板，进而准确高效地执行制图程序。
一个完整的制图模板包括参数设置、创建底图、叠加图层、创建文件尾、格式转换等步骤(图 2)。
|
图 2 模板制图流程 Fig. 2 Template mapping process |
1) 制图默认参数、变量设置。包括格网类型、图框类型、文字字体/字号等制图默认参数，以及制图区域和数据的变量定义。如“gmtset HEADER_FONT {HEADER_FONT} set DATA={DATA}”，其中大括号内容为接口输入参数。
2) 创建底图。可选择本地现有底图，修改文件尾，继续追加图层；也可创建文件头，生成自定义底图。创建文件头的模板如“psxy {R} {J} -T -K > {PS}”。
3) 绘制地质、地理图层。模板中以{layer}标记表示用户在人机交互界面选择的断层、板块等图层，由接口循环读入。
4) 绘制数据图层。根据读入的数值跨度，自动分级分类，绘制不同级别的空间分布图、等值线图、趋势变化图等。如“psxy {DATA} {R} {J} {OFFSET} -S{PTTYPE}{PTSIZE} -W{PTFILL} -O -K > > {PS}”。
5) 绘制图例图层。首先借助echo命令，将标题、符号、文字、分割线等多种要素写入文件，再利用pslegend命令绘制图例。
6) 创建文件尾。制图最后，需要创建文件尾，若缺少文件尾，PS文件可继续追加图层，但不可以转换为栅格格式。创建文件尾的命令为“psxy {R} {J} -T -O > > {PS}”。
7) 转换图片。利用ps2raster命令将PS文件转化为jpg、bmp、png等常用图片格式。
基于模板的GMT可视化制图技术大大降低了用户操作复杂度，并提高了绘图效率。调整图像细节时，只需修改模板即可，避免反复调参出错的情况，增加了系统的可维护性。
2.2 Postscript编码支持GMT绘制的PS文件为页面描述语言(Postscript语言)编写的矢量文件，传统方法是转化为jpg、png等常用图片格式，或依赖第三方阅读器(如GSview、sumatraPDF)打开。本系统从底层集成GhostScript接口库，可实现对PS文件的解码阅读、缩放调整及格式转换。制图后可即时查看绘制结果，实现了所见即所得。
Postscript解码包括以下几个步骤：
1) 加载GhostScript动态链接库；
2) 创建实例：gsapi_new_instance()；
3) 利用gsapi_init_with_args()初始化实例，若PS文件中含有中文字符，在初始化前需要执行gsapi_set_arg_encoding命令设置字符编码为GB2312；
4) 解译Postscript代码：利用gsapi_run_*()系列函数对PS文件中的不同信息(文字、图形、图像等)进行解析；
5) 退出进程：gsapi_exit()；
6) 删除实例：gsapi_delete_instance()。
GhostScript 6.5以上版本均支持中文显示，但GMT并未作配套修改，用户仍需手工将GhostScript的gs-CJK字符集添加入GMT的字体配置文件PS_font_info.d中[5]。本系统预先完成了中文字符的匹配工作，用户只需安装覆盖字体配置文件即可。
2.3 GMT与MapInfo环境的对接MapSIS系统最初基于MapInfo平台开发，输出数据包括MapInfo表文件(包括*.Tab、*.Dat、*.Map、*.ID)、Surfer格网数据(*.grd)、常见地震目录(如*.eqt、*.txt、*.csv)、其他自定义格式文本文件等。而GMT的输入数据则需特定排列的ASCII表格数据、netCDF4格式的nc或grd文件、特定规范的颜色表cpt文件。由于在今后一段时间内，MapSIS系统将保留MapInfo平台供授权用户使用，因此MapInfo环境与GMT制图之间的数据交换是必不可少的环节。本子系统设计的数据交换接口如图 3所示。
|
图 3 数据交换接口 Fig. 3 Data exchange interface |
MapInfo Tab等表文件内部数据结构未公开，只能通过*.Mif、*.Mid交换格式进行对外转换。其中*.Mif保存着表结构和对象的空间信息(坐标、符号样式)，*.Mid则按同样顺序保存对象的属性信息。MIF2GMT接口对输入的Mif和Mid文件进行解析，提取点(POINT)、线(PLINE、PLINE MULTIPLE)、面(REGION)数据的空间信息以及对应的属性信息，输出为GMT直接可用的ASCII文本文件。
2.3.2 格网文件数据交换接口原MapSIS采用的格网数据是Surfer格式的grd文件，与GMT支持的netCDF格式grd文件后缀相同但格式不同。Sufer2GRD接口借助GMT的“xyz2grd”命令，将Surfer格网数据转为netCDF格式。
首先解析文件头中的格网行数Nrow、列数Ncol和经纬度范围Xmin、Xmax、Ymin、Ymax，计算格网间隔：
| $ {I_x} = ({X_{{\rm{max}}}}-{X_{{\rm{min}}}})/({N_{{\rm{col}}}}-1) $ | (1) |
| $ {I_y} = ({Y_{{\rm{max}}}}-{Y_{{\rm{min}}}})/({N_{{\rm{row}}}}-1) $ | (2) |
之后调用xyz2grd命令，命令格式为“xyz2grd-RXmin/Xmax/Ymin/Ymax-IIx/Iy src.grd -H5 -N-9999 -ZBLA -Gout.grd”。其中，-R选项设置了格网数据的空间范围，-I选项设置了转化后数据的格网间隔，-H选项跳过了原文件的文件头，-N选项定义了空值的表示方法，-Z选项定义了图像的存储顺序和数据格式。
2.3.3 ASCII数据交换接口ASCII数据交换接口PREDATA是将原有的ASCII编码文件重新编辑，使其满足GMT对输入文件的列数和每列数值意义的要求。该接口有多个重载，如绘制时序图功能，PREDATA的输出数据格式为“YYYYMMDDThh:mm:ss DataValue”，即第一列为时间，第二列为数值；绘制震中分布图功能，PREDATA的输出数据格式为“Lat, Long, Mag, Date”，即数据列分别为纬度、经度、震级、发震距离当前的天数。
3 测试与结果分析 3.1 功能性测试以绘制全国水位变化空间分布图为例。使用GMT制图子系统时，首先需在“工作参数”菜单中切换“GMT地图”环境选项。之后连接相应前兆数据库，选择测项“动水位”及“静水位”，设置采样率和数据时间，从数据库中读取“差值数据”(图 4)。数据转换接口PREDATA将数据库中存储的水位日变数据按照GMT格式写入本地文件。
|
图 4 从数据库中读取数据界面 Fig. 4 Form of fetching data from database |
点击“绘分布点图”，弹出“绘制空间分布图”的人机交互式对话框(图 5(a))。设定数据绘制的符号、颜色、大小，添加地质地理要素图层，最终产出“全国水位日变空间分布图”的GMT图件(图 5(b))。
|
图 5 人机交互制图参数选择对话框及全国水位日变空间分布图制图结果 Fig. 5 Parameter selecting dialog and daily variation of national water level |
其他可产出的图件类型如表 1所示。所有产品均可通过对话框形式进行人机交互的数据选择和参数调整，调用后台模板进行GMT绘图。产出图件类型基本覆盖了地震监测预报人员的日常需求。目前基于模板的GMT可视化制图接口已在多个业务系统得到应用和检验，如地震前兆数据加工处理系统[6]、数字前兆异常分析系统[7]、地震前兆台网数据跟踪分析系统[8]、基于B/S的地震分析预报信息管理与服务系统(新版APnet)[9]。
|
表 1 产出图件类型 Tab. 1 List of all types of maps |
GMT制图子系统的性能优劣取决于制图结果的显示效率。本文测试了3种PS图件屏幕显示方案：1)执行ps2raster命令将PS文件转化为jpg格式并在系统中显示；2)调用第三方工具(本次测试采用Gsview 5.0软件)打开PS文件；3)在本系统中直接解译PS文件并显示。测试采用Intel Core i5 2.5 GHz、内存8 G的PC机，对图 5(b)进行读取并显示，测试结果见表 2。
|
|
表 2 不同方案PS图件显示效率比较 Tab. 2 The comparison among 3 methods of displaying PS map |
方案1采用GMT自带的ps2raster模块将PS文件转化为jpg等栅格形式，供系统调用显示。该方案转换效率较低，且无法查看矢量文件。而借助商业软件直接读取PS文件，虽能快速显示矢量图件，但无法解决版权问题，且不适用于系统的集成。本系统采用第3种方案，利用GhostScript接口直接解译PS文件，无需用户额外安装第三方工具，即可高效显示矢量图件。
4 结语基于GMT的MapSIS制图子系统引入了可视化的GMT人机交互制图技术，使用户在不了解GMT命令的情况下，也能够通过直观的可视化方式绘制出丰富美观的GMT图件，并可在系统中实时查看矢量格式的GMT图件，进一步提升了MapSIS系统的实用性。基于语法稳定性考虑，本系统暂时采用GMT4进行制图。下一步，我们还将升级GMT版本，并引入开源GIS，实现更全面的空间查询、空间分析功能；实施定制化桌面端，解决系统的臃肿和维护问题，为实现制图标准规范化、知识产权自主化、成果共享开放化建立新技术应用典范，为分析预报技术和方法奠定平台基础。
| [1] |
路远忠, 李胜乐, 邓志辉, 等. 基于GIS的地震分析预报系统[M]. 成都: 成都地图出版社, 2002 (Lu Yuanzhong, Li Shengle, Deng Zhihui, et al. The Analytics System of Earthquake Prediction Basing on GIS[M]. Chengdu: Chengdu Cartographic Publishing House, 2002)
(  0) |
| [2] |
杜学彬, 刘君, 崔腾发, 等. 两次近距离大震前成都台视电阻率重现性、相似性和各向异性变化[J]. 地球物理学报, 2015, 58(2): 576-588 (Du Xuebin, Liu Jun, Cui Tengfa, et al. Repeatablility, Similarity and Anisotropy Changes in Apparent Resistivity Recorded by Station Chengdu at Near Distances before Two Great Earthquakes[J]. Chinese J Geophys, 2015, 58(2): 576-588 DOI:10.6038/cjg20150220)
( 0) |
| [3] |
Wessel P, Smith W H F. GMT Version 4. 5. 15 Technical Reference and Cookbook[EB/OL]. http://gmt.soest.hawaii.edu/gmt4/gmt/html/GMT_Docs.html, 2016
( 0) |
| [4] |
张永奇, 丁晓光, 韩晓飞. GMT软件在大地测量及地震成果表达中的应用[J]. 山西建筑, 2013, 39(6): 205-207 (Zhang Yongqi, Ding Xiaoguang, Han Xiaofei. Application of GMT Software in Geodetic Surveying & Seismic Results Expression[J]. Shanxi Architecture, 2013, 39(6): 205-207)
( 0) |
| [5] |
苏鹤军, 张慧, 李晨桦, 等. GMT绘图软件汉字库配置技术应用研究[J]. 地震工程学报, 2013, 35(4): 928-935 (Su Hejun, Zhang Hui, Li Chenhua, et al. Configuration Technology of Chinese Character Library of Generic Mapping Tool Drawing in Windows System[J]. China Earthquake Engineering Journal, 2013, 35(4): 928-935)
( 0) |
| [6] |
苏融. 前兆数据处理与加工系统的设计与实现[D]. 武汉: 中国地震局地震研究所, 2014 (Su Rong. Design and Implementation of Precursor Data Processing and Product Processing System[D]. Wuhan: Institute of Seismology, CEA, 2014) http://cdmd.cnki.com.cn/Article/CDMD-85404-1015526017.htm
( 0) |
| [7] |
刘珠妹, 李盛乐, 刘坚, 等. 数字前兆资料分析处理系统的设计与研制[J]. 震灾防御技术, 2013, 8(4): 416-422 (Liu Zhumei, Li Shengle, Liu Jian, et al. Design and Development of the Digital Precursory Data Analysis and Processing System[J]. Technology for Earthquake Disaster Prevention, 2013, 8(4): 416-422 DOI:10.11899/zzfy20130409)
( 0) |
| [8] |
刘高川, 李正媛, 王建国, 等. 地震前兆台网数据跟踪分析平台设计[J]. 大地测量与地球动力学, 2016, 36(9): 841-846 (Liu Gaochuan, Li Zhengyuan, Wang Jianguo, et al. The Design of the Tracking and Analysis Platform for Earthquake Precursor Observation Data[J]. Journal of Geodesy and Geodynamics, 2016, 36(9): 841-846)
( 0) |
| [9] |
郭啟倩. 基于B/S的地震分析预报信息管理与服务系统设计与实现[D]. 武汉: 中国地震局地震研究所, 2015 (Guo Qiqian. Design and Implementation of Seismic Analysis and rediction Information Management and Service System Based on B/S[D]. Wuhan: Institute of Seismology, CEA, 2015)
( 0) |
2. The Second Monitoring and Application Center, CEA, 316 Xiying Road, Xi'an 710054, China