引言

地壳应力状态分析在地学研究领域中占有重要地位，绘制地壳应力图并准确表达应力场分布特征一直是科学家们关心的问题。目前国内外专家多采用世界应力图计划（The World Stress Map，WSM）提供的可视化工具CASMO和CASMI实现地应力图的绘制（Roth等，2001；Tingay等，2005；胡幸平，2018），其中CASMO为在线绘图网页，CASMI为单机版人机交互绘图软件。此外，国内相关科研机构还基于商业GIS软件（如ArcGIS、MapInfo等）开发了单机版系统，用于绘制地应力图，如基于GIS的地震分析预报系统（陆远忠等，2000）、中国大陆地壳应力环境基础数据库查询分析系统（崔效锋等，2004）、基于WebGIS的中国大陆地壳应力环境系统（郝天平等，2012）等。近年来，上述工具软件在地应力研究中发挥了重要辅助作用。

然而，上述软件均由用户通过人机交互控件输入或选择参数，无法添加自定义要素精修图形细节，如果用户额外添加图层，需借助ArcGIS、Illustrator、Photoshop等商业软件编辑。上述软件在降低绘制应力图难度的同时失去了绘图灵活性，间接增加了用户绘图的复杂度。人机交互式绘图难以实现批量化图件处理，影响震后地应力场应急产品的快速产出，不利于跟进应急管理部对地震应急产品快速产出的要求。为解决以上问题，本文编写易于理解、支持代码改写且利于批量化处理的程序，为用户绘制应力图提供技术性帮助。

GMT作为脚本式绘图工具，同时也是CASMO、CASMI等软件的底层工具，具有开源免费、可精确控制绘图细节、易批量制图的优点，目前已在大地测量及地震成果中得到广泛应用（占伟等，2010；张永奇等，2013；Fischer等，2014；刘珠妹等，2018；Seredkina等，2018）。然而，GMT脚本不易理解和学习，且其内置符号不满足国际通用地应力符号（图 1）绘制需要。因此，本文研究GMT自定义符号设计方法与传参技术，设计利用GMT软件进行地应力图绘制的整体流程，最终实现地应力数据矢量和应力状态信息的准确表达，极大地简化了用户批量绘制应力图的工作。本文研究成果已在由中国地震台网中心主导、湖北省地震局开发的《震后趋势判定技术系统》中得到应用，实现了震情触发下的地应力图智能成图和产品自动发布，以帮助科研人员快速分析地震周边构造应力场，从而提高震后趋势会商效率。

图 1  全球应力分布图（2016）

Figure 1.  World stress map(2016)

下载: 全尺寸图片 幻灯片

1.   GMT软件概况

GMT是由Wessel等（1998）在哥伦比亚读研期间共同开发的一款开源免费、兼容性强的通用制图软件，具有强大的绘图和数据处理功能。开发者遵循Unix模块化设计，将一个多功能的复杂程序划分成若干个简单、功能单一的程序模块，有利于程序优化、分工和调试。

与目前常用的绘图软件ArcGIS、MapGIS、MapSIS相比，GMT通过纯命令式操作，不仅内存占有率低、执行效率高，且方便用户批量化产出图件。GMT支持30多种投影方式的地图绘制，默认输出格式为PS（PostScript）文件，优点是可供跨平台、高质量的图形打印。

GMT常用于绘制二维网格文件，除内置颜色及符号供用户绘图直接使用外，GMT还支持自定义颜色及自定义符号，目的是提高绘图的灵活度。其中颜色由色标文件CPT（Color Palette Table）控制，自定义符号由符号文件决定。相比创建CPT文件，符号文件的自定义更复杂，要求用户使用GMT中特有的宏命令设计，这是本文绘制应力符号的关键部分。

时至今日，GMT软件仍在持续维护更新，用户可在Unix、macOS、Windows系统下安装配置后使用GMT脚本进行绘图和数据处理。

2.   应力图实现思路

国际上通用的应力图绘制样式如图 1所示（Heidbach等，2018），该图显示WSM2016数据库发布的全球A-C质量应力分布，图中表达了7类应力符号，分别对应震源机制解、钻孔崩落、钻孔诱发张裂隙、钻孔槽、套芯解除、水压致裂、断层活动7个地应力指标。符号长度和颜色分别对应数据质量评价结果（质量越好则符号越大）和数据指示的应力机制（即应力状态，红色为正断层，绿色为走滑断层，蓝色为逆断层，黑色为不确定性断层）。

地应力图中应力符号较复杂，需考虑方位指向性，GMT自带的部分标准几何符号（图 2）无法直接表示应力符号，因此本文利用GMT宏语言环境下的自定义符号技术设计7个独立的符号文件filename.def，用于生成上述7个应力指标符号。设计过程包括单个符号文件的设计、自定义符号的传参技术、自定义符号的表达。

图 2  GMT中的内置符号

Figure 2.  Custom symbols in GMT

下载: 全尺寸图片 幻灯片

2.1   单个符号文件（filename.def）的设计

自定义符号文件（后缀为.def）实质上是在宏语言的环境中创建的。符号文件中存储多行宏命令，可将多个独立要素重组为复杂的自定义符号。每行命令格式为“参数代码”，表示绘制点、线、面等独立要素，其中“参数”定义了要素的坐标、尺寸、角度等信息（如圆的参数为x y size），“代码”表示绘制要素的类型（如圆的代码为c）。“参数”和“代码”之间以空格或逗号分隔，如应力指标震源机制解符号由直线和圆重组而成（见图 3（a）），其符号文件focal_mec.def的命令结构为：

图 3  应力指标符号

Figure 3.  Stress indicator symbols

下载: 全尺寸图片 幻灯片

0 0 1 y-W2p

0 0 0.2c-Gblack

其中第1行表示以（0，0）为中心坐标，画长度为1、粗细为2个像素的竖直线（代码为y）；第2行表示以（0，0）为圆心坐标，画直径为0.2、填充色为黑色的圆（代码为c）。

GMT提供各类宏命令用于自定义符号的绘制，且支持各种复杂符号样式的组合绘制。自定义符号语言中可用的部分绘图命令如表 1所示，利用以上原理，本文设计的7个应力指标符号如图 3（b）所示。

表 1  自定义符号语言中可用的部分绘图命令

Table 1.  Some drawing commands available in symbolic languages

名称	代码	目的	参数
旋转	R	旋转坐标系	${\rm{ \mathsf{ α} [a]}}$
圆形	c	画圆	$ \rm{x}，\rm{y}，\rm{size}$
菱形	d	画菱形	$ \rm{x}，\rm{y}，\rm{size}$
六边形	h	画六边形	$ \rm{x}，\rm{y}，\rm{size}$
倒三角形	i	画倒三角形	$ \rm{x}，\rm{y}，\rm{size}$
文字	l	添加文字	$ \rm{x}，\rm{y}，\rm{string}$
绘制矩形	r	绘制矩形	$ \rm{x}，\rm{y}，\rm{windth}，\rm{heigth}$
正方形	s	绘制正方形	$ \rm{x}，\rm{y}，\rm{size}$
正三角形	t	绘制正三角形	$ \rm{x}，\rm{y}，\rm{size}$
十字架	x	绘制一个十字架	$ \rm{x}，\rm{y}，\rm{size}$
x轴上的线段	-	在x轴上绘制线段	$ \rm{x}，\rm{y}，\rm{size}$
y轴上的线段	y	在y轴上绘制线段	$ \rm{x}，\rm{y}，\rm{size}$

下载: 导出CSV 

| 显示表格

2.2   自定义符号传参技术

根据前文所述的符号文件设计原理绘制出的符号较单一，不具备任何角度信息。而方位指向性恰恰是应力符号的本质特征，用于反映研究区域水平最大主应力的优势分布。因此本文需使用自定义符号传参技术，即引入“符号变量”将水平最大主应力方位角赋予符号自身。

“符号变量”在宏语言环境下由指令“N”定义，其在filename.def文件中的命令格式为“N：n[types]”，其中n为符号变量个数；types为指令N定义的符号变量类型，包括$ \rm{a}、\rm{l}、\rm{r}、\rm{s}、\rm{o}$，分别表示地理方位角、符号长度、旋转角度、字符串类型、其他类型。仍以应力指标震源机制解符号为例，在focal_mec.def中宏命令前添加以下语句：

N：1 a

＄1 R

-90 R

“N：1 a”表示该文件定义了1个符号变量，类型为方位角；“＄1 R”表示引用第1个变量值旋转坐标系；“-90 R”表示将默认的笛卡尔坐标系逆时针旋转90°到地理坐标系。

2.3   应力符号的表达

创建好的符号文件被放入GMT脚本所在的文件夹中，由模块psxy调用绘制符号，其调用格式为-Sk filename/size，其中filename为符号文件名，size表示符号大小，符号颜色直接由-G选项决定，红色对应正断层，绿色对应走滑断层，蓝色对应逆断层。如应力指标震源机制解符号的方向分别为45°、100°、120°，应力数据质量分别为C、B、A，应力状态分别为正断层、逆断层、走滑断层。应力方向、数据质量及应力状态在表达效果上如图 4所示。

图 4  具有方向、质量、应力状态信息的应力指标震源机制解符号

Figure 4.  Symbols of seismic mechanism with direction, data quality, stress state information

下载: 全尺寸图片 幻灯片

2.4   GMT脚本绘制应力图总流程

以龙门山应力区为例，在Windows平台使用GMT5.4.5绘制应力图的总流程如下（其中数据源采用世界应力图计划发布的WSM2016数据库）：

（1）添加地理数据

rem绘制地形底图

gmt psbasemap-R101.2/106/29.5/33.5-JM5i-Ba-K > stressmap.ps

gmt grdimage-R-J earth_relief_ 30s.grd-Cglobe-I-K-O > > stressmap.ps

rem绘制省界

gmt psxy CN-border-La.dat-J-R-W0.5p-O-K > > stressmap.ps

rem添加震中所在的省会城市

gmt psxy CN-capitals.dat-J-R-Sa0.3c-W0.2p, black-Gred-K-O > > stressmap.ps

gmt pstext CN-capitals.dat-J-R-F+f14p, 39, red+j-Dj0.15c/0.15c -K-O > > stressmap.ps

rem绘制市名市界

gmt psxy %CITY%-R%R%-J%J%-W0.25p, 100-O-K > > stressmap.ps

gmt pstext %CITYNAME%-R%R%-J%J%-D0/0.2-F+f7, 39, gray50 -O-K > > stressmap.ps

为直观表达应力区地理信息，本例添加了地形数据、省市边界数据等。绘制地形图采用30s高精度网格数据earth_relief_30s.grd，地形图颜色由色标文件globe控制。

（2）添加应力符号

rem绘制应力符号

gawk-F, "{if(＄3~/FM./ & & (＄6==\"C\") & & (＄4 == \"TF\" || ＄4 == \"TS\")) print ＄2, ＄1, ＄5}" stress_data_longmenshan.csv | gmt psxy-R-J-Skfocal_mec/1-W0.5p, blue-Gblue-K-O > > stressmap.ps

rem绘制断层线

gmt psxy XuFault.gmt-R-J-Sqd3c:+Lh+f11, STHeiti-Regular--GB-EUC-H-W0.5p, red-O-K > > stressmap.ps

stress_data_longmenshan.csv为龙门山断裂带区域应力数据表，数据格式如表 2所示，表 2中共有6列数据，分别为纬度、经度、应力指标类型、应力状态、水平最大主应力方位、数据质量。

表 2  应力数据

Table 2.  Stress data

LAT	LON	TYPE	REGIME	S1AZ
32.500	104.3	FMS	TS	108	C
31.500	104.0	FMS	TF	116	C
31.155	103.7	BO	TF	121	C
…	…	…	…	…	…

下载: 导出CSV 

| 显示表格

gawk是在Windows下使用的数据处理工具，在此处用于筛选应力指标类型为震源机制解（代码为FMA、FMF、FMS）、应力状态为逆断层（代码为TF、TS）、应力数据质量为C的资料，并传递给GMT绘图模块psxy，由选项-Sk调用符号文件绘制应力符号。应力符号方向、大小、颜色分别对应水平最大主应力方位、数据质量和应力状态，应力符号的样式与图 1所示的国际标准保持一致。

另外，本文为增强龙门山区域应力状态可视度，添加了断层线（以红色标记）包括重要断层的名称（如文县断裂），其中断层数据由中国地震局地质研究所活动断层探测数据汇交与共享管理中心提供（徐锡伟等，2016）。

3.   GMT应力图绘制的结果与分析

本文通过GMT自定义符号技术解决了用户制作应力场分布图时绘制应力符号的难题，实现了在Windows10系统下使用GMT5.4.5绘制符合标准的应力场分布图（图 5）。相比国内其他学者（张红艳等，2013）通过传统方式绘制的应力图（图 6）和WSM可视化工具CASMO（http://www.world-stress-map.org/casmo/）在线绘制的应力图（图 7、图 8），本文使用GMT绘制应力图的方法在绘图技巧和信息表达效果上均有所提高。

图 5  基于GMT软件绘制的龙门山应力区应力场分布

Figure 5.  Distribution of stress field in Longmen mountain stress zone based on GMT software

下载: 全尺寸图片 幻灯片

图 6  传统方式绘制的龙门山应力区应力场分布

Figure 6.  Distribution of stress field in Longmen mountain stress zone by traditional methods

下载: 全尺寸图片 幻灯片

图 7  基于CASMO可视化工具绘制的龙门山应力区应力场分布

Figure 7.  Distribution of stress field in Longmen mountain stress zone based on CASMO visualization tool

下载: 全尺寸图片 幻灯片

图 8  CASMO在线绘制地应力图首页

Figure 8.  The website homepage of CASMO used drawing stress map online

下载: 全尺寸图片 幻灯片

3.1   绘图灵活

在相同的应力数据背景下，与CASMO在线绘制网站相比，GMT软件在绘图方式上具有更强的灵活性。

（1）CASMO界面信息主要包括应力数据的筛选、绘图区域设置、投影方式、图例位置及输出格式（图 8）。优点是模板固定、出图快，缺点是无法灵活添加其他自定义数据。虽然CASMO底层依然调用了GMT，但将其封装为可交互式网页，降低了GMT脚本式绘图的灵活度，图 5所示的地形、断层等图层信息的直观表达在CASMO模板下无法实现。

（2）由于CASMO模板已固定，图形涉及的具体细节不易控制，如图 7中坐标轴刻度间隔、底图颜色等常规地图属性设置过于简单。而GMT是纯命令行，所有绘图操作均由命令实现，不仅执行效率高，且能精确控制图形显示，如线条属性、坐标轴设置、画笔属性等。

3.2   信息丰富

相比国内传统方式，本文使用GMT绘制的应力图在信息表达上更直观，具体体现在以下方面：

（1）比较图 5、图 6中的应力分布情况可知，图 5数据分布较密集，经统计达231条记录，其中震源机制解数据211个，比例高达91%，可有效反映地壳深部应力信息。而图 6中显示的数据共86个，其中震源机制解为31个（张红艳等，2013）。出现此现象的原因可能是图 6数据是从中国大陆地壳应力环境基础数据库（谢富仁等，2004）中筛选出来的，而该基础数据库中震源机制解的地震发震时间为1920—2005年，距离芦山地震有近8年的历史地震震源机制解缺失。

（2）相比图 6中大小一致的应力符号，图 5中的符号有大小之分，符号越大代表应力数据的质量越好。应力符号的规格设计可帮助用户直观区分不同指标类型的应力数据质量。

（3）图 6中的应力符号仅具有方向性，而图 5中的应力符号除含有方向、大小的意义外，其颜色的分布还可直观表达研究区域应力状态，如图中红色表示正断层，绿色表示走滑断层，蓝色表示逆断层，黑色表示不确定性，4种断层类型是由构造应力方位的空间取向定义的（Zoback，1992）。龙门山应力区覆盖范围内的断层类型统计如表 3所示，逆断层和走滑断层占比较大，表明该断裂带发生的地震以逆冲型、走滑型为主。

表 3  龙门山应力区断层类型统计

Table 3.  Statistic of fault types in Longmen mountain stress area

项目	正断层	逆断层	走滑断层	不确定型
断层数目（条）	12	113	105	1
占总断层数目百分比（%）	5%	49%	45%	1%

下载: 导出CSV 

| 显示表格

4.   地应力图的自动产出

为简化应急产品绘制程序，避免由人工编辑与整理带来的信息错漏、时间效率低等问题，亟须实现地震科学产品的自动化产出与发布。为此本文在自定义符号技术研究的基础上，使用C#语言与GMT脚本联合编程，实现了地震区域构造应力场应急产品的智能化产出。目前该技术已应用于由中国地震台网中心主导、湖北省地震局开发的《震后趋势判定技术系统》中。

此系统可由震情自动触发或人工输入地震参数手动触发，自动产出地震应急产品。除本文提到的地应力图外，还包括震中分布图、震中附近历史震源机制、震中附近地震序列类型、历史地震活动统计等震情会商产品。最终产品在一定权限控制下以Word版和PPT版的形式发布，推送到服务器及微信客户端上（图 9）。

图 9  微信客户端界面

Figure 9.  WeChat client interface

下载: 全尺寸图片 幻灯片

为了更好地与已有研究结果进行对比，本文以2013年4月20日芦山7.0级地震为例，手动输入参数触发系统。该系统自动生成的芦山地震周边构造应力场图件如图 10所示，除区域应力场分布图外，还包括优势方位玫瑰花样图及主应力轴空间取向图。由优势方位玫瑰花样图可知，水平最大主应力分布大体符合已有研究（张红艳等，2013）计算的龙门山应力区最大主应力方向。结合应力分布图和主应力轴空间取向分析，该断裂带为大型逆冲断层，以近EW向挤压、近NS向拉张为主要特征，与已有研究结论一致（崔进业，2018；张红艳等，2013；张红艳，2015）。

图 10  芦山地震区域构造应力场图

Figure 10.  Structural stress field map of Lushan earthquake area

下载: 全尺寸图片 幻灯片

震后趋势判定技术系统已于2018年10月开始试运行，至今仍在不断维护与更新中，在地震趋势会商研判会上发挥了重要作用。

5.   结论

使用GMT绘制地应力图解决了用户依赖商业软件单机模式生成地应力图产品带来的滞后性、封闭性、复杂性等问题，本文研究成果可总结为以下几点：

（1）利用GMT自定义符号传参技术，本文创建了7个地学界通用应力指标符号标准的符号文件，既保证了脚本绘图的灵活性，又实现了地应力数据可视化表达的丰富性。

（2）通过调研国内外地应力分布图的符号模式和制图要素规范，本文设计了完整的地应力场GMT绘图流程，用于地震应急产品的生产。使用人员可在完全不了解GMT命令的情况下替换图层数据，生成制图规范、要素完整的地应力场应急产品。

（3）本文将地应力场图的绘制成果集成到《震后趋势判定技术系统》中，在一定程度上降低了人工处理图件信息的出错率，提高了震中区域构造应力场图等地震应急产品的智能处理程度，已应用于地震趋势会商工作中。