震例研究中前兆数据曲线快速矢量化程序

2015-11-27石富强

华北地震科学 2015年1期

石富强，戴勇，张博

（1.陕西省地震局，西安 710068；2.内蒙古自治区地震局，呼和浩特 010010；3.甘肃省地震局，兰州 730000）

0 引言

历史震例的回溯性研究是地震分析预测的学习过程，也是对历史震例重新认识、不断升华的过程。广大地震科技工作者在长期的监测预报中反复推敲历史震例，寻找其中的共性问题和个性问题，并做了大量的探索、挖掘和深化研究工作［1－8］。在地震预报还处在探索阶段的今天，这项工作还将持续下去。但同时也面临着另一个问题：过去很多历史震例，尽管在当时震例总结时给出了很好的曲线形态，也汇编成册（如：《中国震例》）。但由于时间久远、观测手段革新、观测仪器更新换代、观测台点更换、数据库升级换代等等问题，过去好多台站的原始观测数据已经丢失，这无疑给震例再认识带来巨大困难。如何从有限的《中国震例》图片中提取高质量的前兆数据，是震例再研究的一项重要的基础性工作。另外，提取数据的质量直接决定着异常再认识的科学性和准确性。

目前互联网上关于曲线矢量化的程序也很多，如：Origin里的Digitize.opk插件，R2V32，Engauge Digitizer，GetdataDraghDigitizer，Graph Digitizer Scout以及Image2data等等。但作者在试用时发现存在以下问题：①这些软件针对的都是横轴为实数时间序列的曲线，而震例里面的时间序列，横轴基本都是年月日（如：20080512）式的时间序列，用这些软件提取出来的数据还需要进一步转化才能方便使用；②这些程序多需要手动逐点双击拾取点，非常累人且效率较低（如Digitize.opk和R2V32），同时取出的数据点主观依赖性较大；③在自动矢量化方面这些程序仅能针对相对单一且干净的图片，震例里面的图片大多相对复杂，有多条曲线同图对比或不间断标地震事件以及“仪器检修、调零”等标注字样，用这些软件自动矢量化时还需要借助第三方软件（如：Photoshop）手动擦除标注文字及多余曲线，增加工作量且效率不高。

为此，本文作者基于Matlab编写了一种能够快速从现有图片文件中提取历史震例前兆曲线的程序，并将其命名为Quick＿Pick（程序下载地址：http：／／pan.baidu.com／s／1c0vkmze）。

1 程序基本原理

图像的亮度值直接存于图像数组中。对于RGB图像，该数组为M×N×3阶矩阵；对于灰度图像，该数组为M×N 阶矩阵。其中，M、N 表示图像像素的行、列数。矩阵中的元素代表该点的亮度值，其范围为［0，255］，分别代表黑色和白色。

1.1 图片的读取及曲线数据筛选

程序首先读入需要被矢量化的图片，并利用matlab命令“imadjust”提高图像对比度。一般图片中，曲线基本为黑色或者彩色，而背景基本为白色，因此我们在程序里设置一个阈值（如C＝254），将亮度值小于阈值C 的像素点亮度全部修改为“0”（黑色），将亮度值大于阈值C 的像素点亮度全部修改为“255”（白色）。这样可以很容易从巨大的像素矩阵中找出疑似曲线像素点的集合R，当然其中可能包含其他曲线或者地震标注等其他文字信息。为了进一步剔除其他曲线以及多余的文字标注，我们在显示的像素点集R 图像上手动拾取一个能够包围目标曲线的多边形框，利用matlab自带的“inpolygon”命令自动拾取多边形内的像素点坐标（图2d）。至此，基本完成了数据的筛选。

1.2 像素点坐标向时间轴（T 轴）和实数轴（y轴）的转换

对于横轴为实数序列的曲线来说，将1.1节拾取的像素点坐标直接做线性等比例转换便可以得到曲线矢量化后的数据。但对于横轴时间序列为“YYMMDD”格式的地震前兆观测数据，则会涉及到闰年（366天）和平年（365天）等问题。为此，我们利用matlab 自带命令“datenum”将指定日期“YYMMDD”转换为实数时间序列；然后再将横轴和纵轴的像素点坐标（i，j）分别线性等比例转化为实数时间序列坐标下的（t，y）；最后将横轴数据转到“YYMMDD”格式下的（T，y）。至此，基本完成了图像的初步矢量化工作并提取保存了曲线数据。

1.3 数据平滑处理

图1 直接提取的像素点坐标及其简单平均结果

上述过程得到的曲线形态在整体上与原始图像非常相像，但放大到细节，会出现如图1 的锯齿形态。这是由于原始图片曲线的粗度占用多个像素点造成在一个横坐标i处拾取了多个（n）纵坐标j值，这将给频谱分析等添加不可估量的数据成分。为避免这种影响，一般在i坐标下对多个j坐标求平均值得到光滑曲线［9］。但这种求平均在遇到“毛刺”、“突跳”等单个脉冲型信号时，会直接将原来的脉冲幅度折半，有失信号的真实性。在现实中，“毛刺”、“突跳”等单个脉冲型信号往往又被认为是可靠的前兆异常［10－13］。因此，有必要尽可能保留并还原这种信号。

为此，我们将前面1.1节识别出来的时间序列通过滑动平均将曲线置于“0”附近，对于给定的横轴像素点坐标i计算拾取的n个点纵坐标平均值A＝mean（y（i，1）：y（i，n））。当A＜0 时，我们认为曲线向下脉冲，取i坐标下的最小值min（y（i，1：n））为脉冲峰值；当A＞0时，我们认为曲线向上脉冲，取i坐标下的最大值max（y（i，1：n））为脉冲峰值；当A＝0时，曲线光滑没有脉冲信号，自动服从前面提到的平均光滑处理，信号值为mean（y（i，1：n））。至此，我们全部完成了数据的筛选、转换以及平滑处理工作。下面，简单介绍程序的使用流程，以便读者有直观深入的理解。

2 程序操作流程简介

假设我们在工作中遇到图2所示的前兆数据曲线。由于各种原因该测项原始数据已经丢失，为了研究和深挖掘的需要，我们必须尽可能还原原始数据，从图1中提取较为可靠的数据来替代原始数据做深入分析。首先，我们将图片截取放在程序目录下面（图3a），运行“Quick＿Pick”读入图片并按照图3b“1”、“2”、“3”顺序，单击鼠标左键选定坐标3个控制点，程序自动跳出图3c输入框；根据前面控制点“1”、“2”、“3”分别输入X 轴起始点坐标“20090801”和“20090831”以及Y 轴起始点“58.10”和“58.45”后点“确定”键，此时程序会将图片上灰度值为“0”的全部可疑信息筛选出来显示为图3d；再单击鼠标左键生成一个尽可能简单的闭合多边形区域（顺时针或逆时针均可），使所有有用信息全部落于多边形区域内，排除“降雨”等标记的多余信息（多边形闭合的时候，双击鼠标左键）。这样，我们便完成了一张无原始数据图片的矢量化工作。提取的数据将自动保存于程序目录下“某台某测项2009年8月观测值.dat”文件中（图3e），可以直接在Mapsis日常分析软件下使用。

图2 某测项2009年8月观测值图片文件

图3 程序使用流程

3 程序检验

由于是事先假定的曲线，在矢量化完成后，我们可以拿出原始数据和提取的数据做对比，将原始数据和提取的数据在同比例的坐标中绘图，如图4所示。对比发现，提取的数据在长尺度上与原始数据形态几乎完全一致。但个别突跳点（如18日晚上19日凌晨的突跳）幅度略有降低，这是由于个别像素点灰度值不够被限所致。将18日至20日数据单独显示如图5所示，可以发现：在原始曲线较光滑的时候，提取数据与原始数据完全一致；在原始曲线变化不稳的时候，略有差距。这是因为，原始数据变化不稳定、上下波动时，在图2中像素点已经靠在一起或者重叠，而程序很难正确将他们分离。但整体来说，提取的信号基本符合原始曲线的变化规律，能够替代原始曲线做近似化的深入分析。