基于OpenCV的模拟地震记录矢量化软件

史鹏飞+张翔+丰继林+王茂发+韩定良

摘 要：随着地震观测技术与手段的不断更新，数字化地震观测逐渐取代传统的模拟地震观测。为了保存历史數据，如何有效提取模拟地震记录中的珍贵信息是地震工作者的研究重点。该项目提出了基于OpenCV的模拟地震记录矢量化算法，使用C#开发相应软件。该软件利用有震波形区域提取方式可以快速、准确地将模拟地震记录矢量化。

关键词：模拟地震记录；数字化；OpenCV；有震波形区域提取；波形追踪

0 引言

世界上第一张地震图要追溯到1889年，英国人米尔恩和尤因利用安置在德国波茨坦的现代地震仪记录下了发生在日本的一次地震。此后70年的时间里，人类观测的地震信息都是以模拟地震记录的方式存在。模拟地震记录蕴含着大量宝贵的地震信息，具有不可估量的价值。直到19世纪70年代，美国率先建成了由十个数字化台站组成的数字化台网。随后数字地震记录技术逐渐取代模拟地震记录技术，一直发展沿用至今。这造成了地震数据应用研究中，模拟地震信息与数字地震信息之间有一条明显的数据断裂。如今大多数模拟地震记录存放在仓库里经受着自然的侵蚀，蕴藏着宝贵信息的模拟地震记录无法被计算机处理。若能有效提取这些模拟地震信息，将会推动地震研究工作快速发展。

1 研究现状

目前为止，对于模拟地震记录的研究主要分为两个阶段：存储与矢量化。

存储是将模拟地震记录中的信息由易损坏的纸质信息转化为方便存储与备份的数字图片信息。2012年开始，哈佛大学投入大量精力进行相关模拟地震记录电子存档工作，该项计划被命名为哈佛大学地震波形归档计划。2015年，河北地震局蒋宏毅等人研发了一套模拟地震记录数字化存储管理系统，实现了模拟地震记录的数字化存储与查询。但数字存储只是将纸质模拟地震记录扫描为数字图片，没有进行矢量化，计算机仍无法识别和处理模拟地震记录，也无法挖掘记录中的宝贵信息‘5]。

许多科研工作者进行模拟地震记录矢量化的研究：Teves-Costa（ 1999）使用CADcore和AutoCAD开展模拟地震记录的矢量化工作；Pintore（2005）基于人工神经网络提出了Teseo方法，可以手动和自动矢量化模拟地震记录；徐涛等（2014）利用Matlab GUI开发了一个交互的模拟地震矢量化程序；王茂发等（2016）提出了CSF算法，并开发了相应的程序。虽然矢量化研究工作取得了显著成绩，但许多矢量化工作只是将一小部分含有复杂波形的图片截取出来单独进行矢量化研究，脱离了整张图片使得后期波形拼接与时间拼接变得难以实现。 鉴于跨平台计算机视觉库OpenCV在人脸识别、车牌识别、动态手势检测等方面的准确性与高效性，因此，本文基于OpenCV设计了新的模拟地震记录矢量化算法，并使用程序设计语言C#开发软件。其中包含分块提取波动区域、自动提取复杂波等特有的功能模块，能够快速、准确地处理整张模拟地震记录。

2 软件矢量化过程

2.1 流程

本软件实现模拟地震记录矢量化的过程共分为8个步骤（图1所示）：加载地震图纸、图纸二值化、搜索起点、有震波形区域提取、扫描波动区域、判断波形类型、拼接和反演。下面给出每一个步骤的具体描述。

2.2加载地震图纸

将整张图纸加载至矢量化软件中是模拟地震记录矢量化最基本的工作。早期的模拟监测是将一段时间（一般为1天）内3个方向（东西、南北、上下）的地震波动记录在一张图纸上。观察整张图（如图2）纸可知，波形密集，包含的信息量巨大，该项目使用C#常用控件PictureBox实现了软件的图片加载功能，可以加载、处理整张模拟地震记录。

2.3 图纸二值化

许多因素如光线、图纸的新旧程度等，会对模拟地震记录扫描而成的图片产生影响，所以必须对图纸进行二值化处理，它是在图片预处理阶段减少图片噪声非常有效的方式。一张灰度的光栅图片类似一个二维的数组，每一个元素在数组中的值在0-255之间。利用OpenCV的Threshold函数对图片进行二值化处理：设置一个阈值（目前实验获得经验数值为155），将大于阈值的像素点置为255，小于阈值的置为0，如图3所示：

2.4搜索起点

为了确保数据的一致性，软件在提取模拟地震记录中的信息时，始终在同一个坐标系下处理。以图片的左上角为原点，图片的上边缘为x轴正值方向，图片的左边缘为y轴正值方向。取原点附近的一个点c（x，0）为起点，沿着垂直方向向下搜索，每次步长为一个像素点。搜素算法具体描述如下：

（I）水平方向

水平方向x值保持不变（公式1），为了避免图片左侧空白边缘对搜索结果的影响，一般x取值稍大于0（实验取值为10）。 X_i=x（i=l，2…m）

（1）

（2）垂直方向

垂直方向y值由0不断变大，每次以一个步长即一个像素向下搜索，并判定像素值的变化。当遇到相邻两个像素点值跳变时记录这个像素点坐标到数组中（x，Y_n1），然后继续向下搜索，再次遇到跳变时继续记录（x，Y_n2）…。

则第一条波的起点坐标为公式（2）：

直到搜索y值到图纸的下边缘，数组中记录的点（x，Y_n1），（x，Y_n2），……，（x，_nm）。（点的编号为1，2…m）这样每一条波形起点的信息都可以确定。第h（l≤h≤m/2）条波的起点坐标信息为（公式3）：

2.5有震波形区域提取

模拟地震记录中最重要的数据是记录地震发生时的复杂波，矢量化工作的难点是矢量化这震波形区域的波形。通过观察发现这些复杂波只是整张图纸很小的一个局部，所以采取有震波形区域提取的方式：先找到复杂波所在的第一个区域，在这个波动区域内进行矢量化操作，接着寻找复杂波所在的第二个区域进行矢量化操作，直至完成所有波动区域的波形矢量化（图4）。所有的矢量化工作基于图纸的局部却没有脱离整体图纸，确保了波形信息与时间信息拼接操作的可行性。如果在整张模拟地震图纸上进行矢量化复杂波的操作，这样既增加难度，也降低效率和准确率。endprint

2.6扫描波动区域

2.6.1 判断波形类别

首先，使用搜索起點算法将波动区域的所有波形起点找到。然后，按照Y值由小变大的顺序以每条波的起点水平扫描这条波。扫描算法具体描述如下：

（1）选取在此区域内第一条波的起点（X_c1，y_c2）。

（2）水平方向，X_C1不变。垂直方向，以一个像素点为步长，先沿着v轴正方向搜索至（ X_C1，Y_C1+△y） （Ay最大为波形宽度的平均值，若在此区域内出现像素跳变点则记录其坐标（ X_c1.y_D1），若不出现跳变点则记录点（ X_C1，y_D1）=（X_C1，Y_C1+△y）；后沿着y轴负方向搜索至（X_C1，Y_C1-△y），若在此区域内出现像素跳变点则记录其坐标（X_C1，y_D2），若不出现跳变点则（X_C1，y_D2）=（X_C1，Y_C1-△y）。然后记录在（X_C1，Y_C1±△y）范围内中点坐标同时记录该点的像素值。

（3）水平方向：X_C2=X_C1+1。

（4）重复（2-3）操作直至波动区域边界。

使用扫描算法扫描平滑波时，由于Ay的限制，可以基本消除与复杂波交叉的影响，追踪到平滑波的中线，所以记录下来平滑波中线的像素基本都为黑点。但由于复杂波是有振幅的，而△_y值最大为线宽，所以扫描复杂波时许多中点的像素点为白色（图5）。因此，当扫描一条波形记录下的白色像素点多于其他波形时，则判定该波为复杂波同时进行标记。

2.6.2 追踪平滑波

判断波形类型后，根据标记信息进行平滑波的跟踪。利用改进的扫描算法追踪平滑波：扫描时，若在区域（y+△y，y-△y）}H现像素跳变点，记录此范围内的中点并将此范围的所有黑色像素点置为白色像素点，反之则不变。这样既保证记录平滑坡中线坐标的正确性，也不会破坏复杂波的波形，追踪效果如图5所示：

2.6.3 去除噪声

平滑波追踪完成后，复杂波追踪变得更加简单。但仪器记录和保存时各种因素可能会导致图片存在干扰信息，为了避免这些噪声对追踪复杂波的影响，必须再次进行去噪处理。首先，利用OpenCV对图片进行SmoothMedian处理，消除椒烟噪声与斑点噪声。其次，利用FindContours函数提取所有图形的轮廓，记录每一个轮廓的id，计算每一个轮廓的周长。最后，通过比较保留周长最长的轮廓，其余全部去除——这样就保存了复杂波信息而去除了噪声，图6是保留的复杂波轮廓。

2.6.4追踪复杂波

如同上面一样，可以计算出复杂波起始处的中点坐标，然后以此点为起点追踪复杂波。在追踪过程中，值得注意的是，在波峰和波谷处需要记录此处的v值上边缘值与下边缘值，这样可以正确反映地震的能量释放情况。追踪复杂波的算法与追踪平滑波的算法相似，不同之处在于：没有△y的限制，当中间点的y值大于相邻两点的v值时（波谷，v轴正方向向下），记录其下边缘的v值；当中间点的v值小于相邻两点的v值日寸（波峰），记录其上边缘的v值；搜素到区域边界，图7中红线为波形的中线，蓝点为提取的关键点。

2.7 拼接和反演

在前期工作中已经实现了拼接与反演，这里只做简单描述。每条波形上的点都以坐标的形式存储在数组中并且标有唯一的id，首先根据点的坐标信息进行波形的拼接，然后根据图纸上标注的时间信息进行时间的拼接，最后将波形信息反演，如图8所示：

3 结束语

该项目设计了基于OpenCV的模拟地震记录矢量化软件，分8个步骤矢量化模拟地震记录的栅格图片，通过给出的测试结果可以看出，反演出来的波形较好的符合原始图纸中的波形。下一步，将进一步提升算法效率和改善软件，期望开发出自动校验模块以确保矢量化的质量。endprint