一种历史地震图波形自动跟踪方法

摘要

一种历史地震图波形自动跟踪方法，属于油气地球物理勘探技术领域。其特征包括：(1)利用微分几何曲线提取算法从地震图中提取出地震波曲线点，以及地震波曲线的局部方向；(2)依据已经提取出来的曲线局部方向将提取出来的曲线点，连成多条曲线段；(3)最后依据地震波曲线的自身特性，如地震波曲线在其拐点处呈现出的对称特性以及地震波曲线围绕其中轴线做往返运动的震动特性，以回溯的思路对地震波曲线段进行跟踪。本发明的效果和益处是提高了地震波波形自动跟踪的精度，无需人工干预，能够处理有交叉线和断点的地震图。

说明书

技术领域

本发明属于油气地球物理勘探技术领域，涉及到对历史地震图的数字化， 特别涉及到数字化过程中波形自动跟踪的方法。

背景技术

在数字化地震观测系统建成之前的百年间，人们通过使用模拟笔的方法对 地震图进行记录，留下了许多宝贵的无可复制的地震资料。这些历史地震图一 般使用纸张、胶带等进行记录，具有明显的保存时间限制，为了方便使用现代 地震学方法对其进行分析同时也为了延长这些资料的保存时间，我们需要对其 进行数字化。

近年来，越来越多的学者专家投入到历史地震图的数字化工作中。其中如 何高效并且自动的跟踪地震波成为历史地震图数字化的关键步骤。目前已有的 跟踪方法主要分为两大类。第一类跟踪算方法属于半自动范畴，这些方法需要 大量的人机交互过程，过程很繁琐，比较耗时，如2003年Bromirki等设计的 SeisDig系统，以及2005年Pintore等设计的Teseo2系统。第二类跟踪方法可以 实现对地震图波形的全自动跟踪，然而这些方法一般依赖于地震图图像的局部 特征，对于交叉线和断点较多的地震图，他们不能给出很好的结果，如Church 等人在2013年提出的基于形态学的波形自动跟踪算法，Mao等人在2014年提 出的基于逐点跟踪策略的波形自动跟踪算法。

发明内容

本发明的目的是提供了一种历史地震图波形自动跟踪方法，该方法将利用 地震波曲线的自身特性，如地震波曲线在其拐点处呈现出的对称特性以及地震 波曲线围绕其中轴线做往返运动的震动特性，采用回溯的思路完成对波形的自 动跟踪。该方法提高了地震波波形自动跟踪的精度，无需人工干预，能够处理 有交叉线和断点的地震图。

本发明的技术方案是：

一、技术方案的原理

(1)技术方案的基本概念。

历史地震图：以滚筒的形式由模拟笔记录到纸张、胶带上的原始地震资料。

目标地震图：历史地震图经过扫描仪得到的光栅图像，是本发明的处理对 象。

跟踪：从目标地震图中，提取并且识别出不同的地震波曲线。

地震波曲线：目标地震图中，地震波波形曲线往往具有一定的线宽，本方 案中，地震波曲线指的是地震波波形曲线的中心线，并且是本方案的跟踪对象。

曲线局部走向：即地震波曲线的切线方向。

拐点：即地震波曲线走向发生变化的点。

地震波中轴线：地震波曲线由模拟笔记录完成，模拟笔在记录过程中沿其 平衡位置上下摆动，本方案称这个平衡位置为地震波中轴线。

(2)技术方案的具体原理。

本发明方法由三部分组成，首先依据微分几何曲线提取算法从目标地震图 中提取出地震波曲线点以及地震波曲线的局部走向；然后结合地震波的局部走 向将提取出来的地震波曲线点连接成多条曲线段；最后依据地震波曲线的自身 特性，如地震波曲线在其拐点处呈现出的对称特性以及地震波曲线围绕其中轴 线做往返运动的震动特性，以回溯的思路对地震波曲线段进行跟踪。下面将对 上述每一步的技术原理进行详细介绍,并在介绍地震波跟踪算法之前，将地震波 的一些性质与其计算方法也给与相应的介绍。

1)地震波曲线点以及曲线局部走向的提取

对于目标地震图I，沿地震波波形曲线的切线方向做灰度剖面，该灰度剖面 的极值点即为地震波曲线点。由极值点的性质可知，对于地震波曲线点来说， 它沿切线方向的一阶方向导数为零二阶方向导数的绝对值较大。本发明取I的二 阶方向导数最大的方向为地震波波形曲线的切线方向，记上述切线方向为 其中为该方向与X轴的夹角。I沿任意方向(cosθ,sinθ)的二 阶方向导数为：I_θθ＝I_xxcos²θ+2I_xycosθsinθ+I_yysin²θ。令可以 求得I的二阶方向导数最大的方向,该方向即为法线方向，如公式(1)所示：

$\hat{\theta }=\arctan \left(\frac{I_{\mathrm{xx}}-I_{\mathrm{yy}}+\sqrt{{(I_{\mathrm{xx}}-I_{\mathrm{yy}})}^2+4I_{\mathrm{xy}}^2}}{2I_{\mathrm{xy}}}\right)---$

其中I_xx,I_xy,I_yx,I_yy为I的各阶偏导。一旦地震波波形曲线的法线方向被计算出 来，与该方向垂直的方向即为地震波曲线的局部走向。

由于目标地震图是离散的，其最小计量单位为像素，在本发明中认为像素 为正方形，对于某一像素点(x,y)，它的定义范围为直 接计算上述灰度剖面的极值点，仅能得到像素级别的精度，为了得到亚像素级 别的精度，我们需要对每个像素点进行插值。图像处理中泰勒展开式是获取亚 像素的一种重要手段。对像素点(x,y)沿方向做泰勒展开：

$\left(\begin{array}{c}\tilde{I}(x+t·\cos \hat{\theta },y+t·\sin \hat{\theta })\approx I(x,y)+(t·\cos \hat{\theta },t·\sin \hat{\theta })\left(\begin{array}{c}{I}_x\\ I_y\end{array}\right)\\ +\frac{1}{2}(t·\cos \hat{\theta },t·\sin \hat{\theta })\left(\begin{array}{cc}{I}_{\mathrm{xx}}& I_{\mathrm{xy}}\\ I_{\mathrm{xy}}& I_{\mathrm{yy}}\end{array}\right)\left(\begin{array}{c}t·\cos \hat{\theta }\\ t·\sin \hat{\theta }\end{array}\right)\end{array}\right)---\left(2\right)$

其中为估计的亚像素点，t为亚像素点与展开点(x,y) 的距离，由于在一个像素内部做展开故t<0.5。这样亚像素级别的极值点即地 震波曲线点可以通过公式(3)、(4)求得。

$\frac{\partial \tilde{I}(x+t·\cos \hat{\theta },y+t·\sin \hat{\theta })}{\partial t}=0---\left(3\right)$

$\left|\frac{{\partial }^2\tilde{I}(x+t·\cos \hat{\theta },y+t·\sin \hat{\theta })}{\partial t^2}\right|>{\gamma }_1---\left(4\right)$

其中γ₁为用户设置的阈值。

2)地震波曲线点的连接策略

记已提取的地震波曲线点的集合为C，从中选取二阶导数较大的点作为连 线的初始种子点。对于每一个初始种子点，本发明将按照图3所示的策略进行 连线。从初始种子点开始，沿当前地震波曲线方向搜索初始种子点的三个邻域 像素内是否存在地震波曲线点，例如初始种子点为(x,y)，当前地震波曲线方向 为θ且θ∈[-22.5°,22.5°]，本发明将搜索(x+1,y-1)，(x+1,y),(x+1,y+1)这 三个邻域像素内是否存在地震波曲线点，如果不存在，本次连线结束，重新选 择新的初始种子点重复连线过程；如果存在，标记这些曲线点为候选种子点， 按公式(5)计算这些候选种子点与初始种子点的代价函数值，删除代价大于阈 值的候选种子点，记该阈值为γ₂；删除后如果一个候选种子点都没剩下，声明 本次连线结束，重新选择新的初始种子点重复连线过程；删除操作后如果还有 多个候选种子点，取代价最小的候选种子点作为最优连线点，而后更新初始种 子点为已经找到的最优连线点；重复上述策略。本过程是一个迭代过程，当搜 索不到候选种子点时，该迭代过程结束。

$\mathrm{Cost}=|{\alpha }_1-{\alpha }_2|,\mathrm{Cost}\in [0,\frac{\pi }{2}]---\left(5\right)$

式中Cost表示代价函数，α₁表示地震波曲线在初始种子点的方向，α₂表示地震 波曲线在候选种子点的方向。

3)地震波曲线的自身性质以及计算方法

在地震波曲线点提取完成后，本方案对提取出来的地震波曲线点进行了初 步连接，将其连成一段段小的曲线段。为了完成地震图波形的跟踪，需要设计 算法对这些小的曲线段进行跟踪，即那些属于同一条地震波曲线的曲线段被划 分为一类。由于地震波曲线的自身性质在跟踪算法中极为重要，本发明在介绍 跟踪算法之前，首先介绍地震波曲线的一些重要的性质。本发明通过对大量地 震图的研究得出地震波曲线的以下三点重要性质。

性质1：对于一笔画下来的相邻的两个地震波曲线段来说，这两条曲线段曲 线走向相差很小。

曲线段走向可按如下策略计算：对于任一条曲线线段L，对其首尾两端M个 曲线点的方向进行加权平均计算其首尾两端的曲线走向，如公式(6)所示，取 L的最左端为其首部，最右端为其尾部。

${\beta }_1=\frac{1}{M}\underset{k=1}{\overset{M}{\Sigma }}d\left(L\left(k\right)\right),{\beta }_2=\frac{1}{M}\underset{k=N-M}{\overset{N}{\Sigma }}d\left(L\left(k\right)\right)---\left(6\right)$

其中β₁，β₂为曲线首部和尾部走向，d(L(k))表示种子线段L在第k个点的局 部方向，N表示种子线段L包含的曲线点数。

性质2：地震波曲线由模拟笔的摆动绘制而成，存在大量的拐点。对于这 些拐点来说，位于拐点处两端的地震波曲线段具有一定的对称性和走向相反的 特性。

关于曲线走向是否相反的判断，通过公式(6)计算曲线的走向并作出判断。 关于两条曲线对称性的判断，本发明采用如下策略：记这两条曲线为L₁，L₂。 将这两条曲线分别等分为N份，计算切分后的每一小段曲线的斜率。而后按公 式(7)判断两条曲线的对称性。公式(7)中K(i)为曲线L₁切分后的第i条曲 线段的斜率，k(i)为曲线L₂切分后的第i条曲线段的斜率，γ₃是用来判断对称 性的阈值，由人工设定。

$\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}|K\left(i\right)-k\left(i\right)|<{\gamma }_3---\left(7\right)$

性质3：每条地震波曲线都有一条属于自己的中轴线，地震波曲线围绕其做 往返运动，并且不同的地震波曲线具有不同的中轴线。

本发明采用IRLS迭代加权最小二乘法对地震波曲线的中轴线进行估计。记 当前地震波曲线为Γ，关于Γ的中轴线具体估计步骤如下：

Step1：对Γ的纵坐标进行平均求和，给出中轴线的初始估计，记为L_inital。

Step2:对Γ的纵坐标进行加权求和，Γ中每个点的权值按公式(8)计算， 记该次估计的中心线为L_next。

Step3:判断|L_next-L_inital|<1是否成立，如果成立转向Step5，否则转入Step4。

Step4：更新L_inital为L_next，转入Step2.

Step5:中轴线估计结束。

${\omega }_i=\frac{1}{\max (\sigma ,|y_i-L_{\mathrm{inital}}\left|\right)}---\left(8\right)$

公式(8)中，ω_i为Γ中第i个点纵坐标的权值，y_i为Γ中第i个点纵坐标， σ是一个很小的值，主要是为了防止除数为零的情况发生，取0.001。

3)基于回溯的跟踪方法

回溯算法也叫做试探法，它是一种系统的搜索问题解的方法。回溯的基本 思想为：从一条路线出发开始跟踪，能进则进，不能进则退回来，换一条路再 试试。本发明将回溯的思路引入到地震图波形的跟踪上。

将整幅地震图看作是一个有向图G；记初步连接地震波曲线点过程中得到 的地震波曲线段为G的有向边。则地震图波形的跟踪问题，可以看做是从有向 图G中按照地震波波形的特点寻找最优路径的问题。对于有向图中最优路径的 选取可以采用回溯的算法给予实现，其基本思路为：从G中的一条有向边出发 往前走，碰到交叉点，首先判断该点是否为拐点。如果不是拐点选择与当前有 向边走向相差最小的一条有向边作为最合适的路径并往前走；如果是拐点执行 回溯的思路，在拐点处对跟踪路径进行分叉，沿拐点处的每条有向边都分出一 条新的路径，记这些路径为Ψ，按照IRLS迭代加权最小二乘法对每条路径的曲 线的中轴线进行估计并计算这些估计的中轴线与当前跟踪的地震波中轴线的差 值，取差值最小的路径最为最合适的路径并往前走。整个跟踪过程是个迭代过 程，每次碰到交叉点，就按照上述规则进行路径选取，当在G中搜索不到新的 有向边时，迭代过程结束。

上述跟踪过程中需要用到当前跟踪的地震波中轴线。我们利用已经跟踪到 的曲线段，通过IRLS迭代加权最小二乘法对该中轴线进行估计。

二、技术方案的步骤

S1获取目标地震图的灰度图片记为I，I的像素灰度级为0到255；

S2将I与高斯模糊核进行卷积，用于压制I中噪声，通过sobel算子计算I 的各阶偏导数I_x,I_y,I_xx,I_xy,I_yy。

S3依据公式(1)计算地震波波形曲线的切线方向，取与该方向垂直的方向 为地震波曲线的局部走向。设置公式(4)中的阈值，依据公式(3)、(4)提取 出地震波曲线点。

S4记步骤S2中提取的曲线点的集合为Φ，对于Φ中每一点沿公式(1)给 出的方向计算其二阶方向导数，对这些二阶方向导数进行归一化处理，取二阶 方向导数大于0.2的曲线点作为连线的起始种子点。对于每一个起始种子点，按 图3所示的流程连接成曲线段,流程图中的阈值设定为0.12*π。

S5将整幅地震图看作是一个有向图G；记步骤S3中得到的所有地震波曲 线段为G的有向边，按照公式(6)计算G中所有有向边的首尾走向。从G中选 取用于跟踪地震波曲线的初始跟踪路径，初始跟踪路径按照如下方案选取：将 有向图G中所有的有向边按照它们的起始点的横坐标大小从小到大进行排序， 取排在前10％的边作为初始跟踪路径。

S5.1从步骤S5中取出第一条初始跟踪路径。

S5.2在初始跟踪路径尾部节点处搜索有向边作为候选的跟踪路径，该搜索 区域为半径为2.5倍的地震波波形曲线宽度的圆形区域，如果一条候选的跟踪路 径都没找到跳到步骤S5.8；如果找到了候选路径，进入步骤S5.3。

S5.3利用本次跟踪过程已经跟踪过的路径，结合IRLS迭代加权最小二乘 法对当前跟踪的地震波曲线的中轴线进行估计。

S5.4按照公式(7)计算所有候选跟踪路径与初始路径之间的相似性，计算 过程中N取4。如果这些候选跟踪路径中存在某条路径和初始路径相似并且这 两条路径走向相反，则认为本次跟踪过程中碰到了拐点，并转入步骤S5.6，否 则转入步骤S5.5。

S5.5从步骤S5.2中所有候选跟踪路径中，选取首部走向和初始跟踪路径尾 部走向差值最小的路径，作为后续跟踪路径，并转入步骤S5.7。

S5.6执行回溯策略，在拐点处进行路径分叉，拐点处的每条候选跟踪路径 都作为后续跟踪路径分出多条新的路径，记这些路径为Ψ，采用IRLS迭代加权 最小二乘法对每条路径的中轴线进行估计并计算这些中轴线与当前跟踪的地震 波中轴线的差值，取差值最小的路径作为本次跟踪过程中的后续跟踪路径，舍 弃其他路径。

S5.7更新初始跟踪线段为找到的后续跟踪路径，转入步骤S5.2.

S5.8本次跟踪过程结束，保存本次跟踪过程中的路径，存入计算机。

S5.9按顺序取步骤S5中未执行跟踪过程的起始种子线段，重复步骤 S5.2-5.8。

本发明的和效果和益处是，本发明将地震图的局部特征和地震波波形的宏 观特性相结合并借鉴回溯的思想，提高了地震波波形自动跟踪的精度，无需大 量的人工干预，可以较好的处理交叉线断点较多的地震图。

附图说明

附图1是经过扫描仪得到的目标地震示意图。

附图2是历史地震图波形自动跟踪算法流程图。

附图3是地震波曲线点连接成线的流程图。

具体实施方式

以下结合技术方案和附图详细叙述本发明的具体实施方式。

实施例：

本例以冰岛地震图的一部分为实验目标，该地震图由大连地震台记录于 2000年7月21号。实验平台CPU为I3-2120,运行环境matlab 2012a。

步骤1、获取实例地震图的灰度图片记为I，I的像素灰度级为0到255

步骤2、将I与σ＝1.4的高斯模糊核与进行卷积，用于压制I中噪声，采用 sobel算子求取卷积后的图像的各阶偏导数，I_x,I_y,I_xx,I_xy,I_yy。

步骤3、依据公式(1)计算地震波波形曲线的切线方向，取与该方向垂直 的方向为地震波曲线的局部走向。设置公式(4)中的阈值γ₁＝1.2，依据公式 (3)、(4)提取出地震波曲线点，对每一个曲线点沿公式(1)给出的方向向计 算其二阶方向导数，从这些曲线点中删除二阶方向导数小于0的点，余下的点 作为最终提取出来的地震波曲线点。

步骤4、对步骤3中地震波曲线点的二阶方向导数进行归一化，取二阶方向 导数大于0.2的曲线点作为连线的起始种子点。对于每一个起始种子点，按图3 所示的流程连接成曲线段，连线过程的γ₃设定为0.12π。

步骤5、将整幅地震图看作是一个有向图G；记步骤S3中得到的地震波曲 线段为G的有向边，按照公式(6)计算G中所有有向边的首尾走向，公式(6) 中M取8。从G中选取用于跟踪地震波曲线的初始跟踪路径，初始跟踪路径按 照如下方案选取：将有向图G中所有的有向边按照它们的起始点的横坐标大小 从小到大进行排序，取排在前10％的有向边作为初始跟踪路径。

步骤5.1、从步骤5中的取出第一条初始路径。

步骤5.2、在初始路径尾部节点处搜索有向边作为候选的跟踪路径，该搜索 区域为半径为15个像素的圆形区域，如果一条候选的跟踪路径都没找到跳到步 骤5.8；如果找到了候选跟踪路径，进入步骤5.3。

步骤5.3利用本次跟踪过程已经跟踪过的路径，结合IRLS迭代加权最小二 乘法对当前跟踪的地震波曲线的中轴线进行估计。

步骤5.4按照公式(7)计算所有候选跟踪路径与初始路径之间的相似性， 计算过程中N取4，γ₂取0.1。如果这些候选跟踪路径中存在某条路径和初始路 径的相似并且这两条路径走向相反，则认为本次跟踪过程中遇到了拐点，转入 步骤5.6，否则转入步骤5.5。

步骤5.5从步骤5.2中所有候选跟踪路径中，选取首部走向和初始跟踪路径 尾部走向差值最小的路径，作为后续跟踪路径，并转入步骤5.7。

步骤5.6执行回溯策略，在拐点处进行路径分叉，拐点处的每条候选跟踪 路径都作为后续跟踪路径分出多条新的路径，记这些路径为Ψ，按照IRLS迭代 加权最小二乘法对每条路径的中轴线进行估计并计算这些估计的中轴线与当前 跟踪的地震波中轴线的差值，取差值最小的路径作为本次跟踪过程中的后续跟 踪路径，舍弃其他路径。

步骤5.7更新初始跟踪线段为找到的后续跟踪路径，转入步骤5.2.

步骤5.8声明本次跟踪过程结束，保存本次跟踪过程中的路径，存入计算 机。

步骤5.9按顺序取步骤5中未执行跟踪过程的跟踪路径，重复步骤5.2-5.9。