一种基于变分原理的地震数据信号去噪方法与流程

本发明涉及石油勘探领域，尤其涉及一种基于变分原理的地震数据信号去噪方法。

背景技术：

地震信号测量的方法不是很多且各有优缺点，目前检索到的估计地震信号的非线性差异测量技术有以下几种：

1、互相关方法

该方法是利用随机信号的平稳性，计算两道信号之间的互相关，通过提取互相关的最大值确定信号的时移量。该方法只适应与平稳信号，计算量大。

2、多项式拟合方法

该方法假设地震信号在局部邻域道内的时差符合一个待定的多项式，通过拟合该多项式提取不同的时间差。该方法需要较大的计算量来估计未知参数，并且假定了不同道之间的时差关系是低阶平滑的，无法处理强变化的时移情况。

3、非线性beamforming方法

该方法利用局部的抛物二阶时差近似去搜寻局部的相关事件，同样地，该方法需要估计未知的非线性系数，计算量巨大，并且基于二阶近似的方法无法适用非规则的信号时移差异。

4、常规的dtw方法。

传统的dtw方法需要通过累积和回溯两步寻找到最短路径，可以提取非线性的时移量，但是计算量大，并且易受噪音影响，陷入局部极值中，使得结果不稳定。

这些方法各有优缺点，但其仍然无法达到人们对地震信号测量操作的效率及稳定性等方面的需求。

技术实现要素：

本发明目的是针对上述问题，提供一种基于变分原理的地震数据信号去噪方法。

为了实现上述目的，本发明的技术方案是：

一种基于变分原理的地震数据信号去噪方法，包括以下步骤：

s1、在地震信号处理模型中输入测试信号f和参考信号g并选取采样点；

s2、将测试信号和参考信号进行归一化处理，构造成为差异矩阵；

s3、将差异矩阵转换为慢度矩阵；

s4、基于费马原理和惠更斯-菲涅尔原理，采用最小路径射线追踪方法在慢度矩阵中提取全局最小旅行时及对应的最短路径；

s5、将获得的最短路径中每个采样点的旅行时时移量施加在测试信号上，获得去噪后的地震数据信号。

进一步的，所述步骤s1中的测试信号f和参考信号g中有效信号的个数要求相同，其两侧端点也相同并且要求测试信号f和参考信号g均具有次序性。

进一步的，所述步骤s2中差异矩阵的构造公式为：

d(x,y)＝[f(x)-g(y)]²,x＝1,...,mt；y＝1,...,nt

其中，mt是测试信号的采样点数，nt是参考信号的采样点数。

进一步的，所述步骤s3中差异矩阵转换为慢度矩阵的转换公式为：

其中，vmin和vmax分别为最小和最大速度值，dmax为差异矩阵的最大值。

进一步的，所述步骤s4中在慢度矩阵中提取全局最小旅行时及对应的最短路径的方法包括以下步骤：

s41、将慢度矩阵分为矩形网格，在矩形网格单元顶点处建立节点，并确定起始节点和终止节点；

s42、将起始节点设为初始震源节点s0，围绕s0搜索所有可能的临近计算节点

s43、计算从当前震源节点s0到计算节点的旅行时，从中选择最小的旅行时记为对应的计算节点记为g1，则：

s44、把当前计算节点g1设为新的震源节点s1，重新搜索当前震源节点的计算节点并计算对应的旅行时，最小旅行时记为对应的计算节点记为g2；

s45、从初始震源节点s0到当前选择的最新计算节点的总旅行时为对应的路径为

s46、重复上述过程，直到最终获得的计算节点为终止节点，获得全局最小旅行时及对应的路径：

与现有技术相比，本发明具有的优点和积极效果是：

本发明中采用基于变分原理的dtw方法测量两段地震信号的差异，有着其他方法不具备的优势，其具体优势和特点表现在以下几个方面：

第一、方法效果的可靠性。本方法通过形成的慢度场进行射线追踪，提取最小旅行时和最短路径，可以有效地避免局部极值的产生，并可有效地压制噪音的干扰，求取的最短路径更稳健。另外，同常规的dtw方法相比，该方法计算效率更高。

第二、理论上进行了提升。将dtw问题转换为两个核心问题，一个是差异矩阵的构建，一个是寻找最优路径问题。差异矩阵的构建可以选择信息的属性进行提取，并选择其他抽象距离来定义；最优路径求取在变分原理的框架下求取，提供了多种方法来求取信号的扰动时差，为非线性时差测量提供新了求解思路。

第三、具有非常广泛的应用价值。本方法处理两段信号的差异，在地震勘探中有诸多这样需求，比如叠前地震数据叠加去噪、插值；成像剖面处理；成像道集后处理；参数反演中模拟数据与实际数据的差异测量等等方面。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为计算节点搜索方式示意图；

图2为本发明的方法流程图；

图3为待处理地震信号的示意图；其中左侧a信号为参考信号，右侧b信号为待处理信号；

图4a为构建的差异矩阵示意图；

图4b为构建的慢度矩阵示意图；

图5a为前两次获得的最短路径示意图；

图5b为最后获得的全局最短路径示意图；

图6为处理后的地震信号示意图；其中左侧a信号为参考信号，右侧b信号为处理后的信号；

图7a为原始的共炮道集；

图7b为处理后的炮道集。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

在地震信号处理中，往往需要对两道信号进行波形匹配，或测量两段信号的差异。一种有效的方法是采用动态时间规划方法测量两段信号的差异，该方法是一个典型的优化问题，其问题在于用满足一定条件的差异矩阵描述测试信号和参考信号时间差异关系，求解这两个信号匹配时，累计距离最小时所对应的最短路径。常规的dtw方法需要通过累积和回溯两步寻找到最短路径，可以提取非线性的时移量，但是计算量大，并且容易陷入局部极值中，不稳定。本发明提出基于变分原理求解dtw，将dtw技术主要分为差异矩阵形成和求取最短路径两个方面。差异矩阵的形成采用信号之间的欧式距离来实现；最短路径的求取则基于变分原理，采用最小路径射线追踪方法提取全局最小旅行时和最短路径。从理论和实现两个方面进行了修改，不仅提高了计算效率，增强了最短路径求取的稳定性，而且为非线性信号测量打开了新的视角。

本发明目的是提供一种基于变分原理的dtw技术测量两段信号的差异的方法，并将该方法应用到叠前地震数据去噪，优化叠加中。

本发明具体包括以下步骤：

1、差异矩阵形成：

首先，差异矩阵的构造前提要求信号满足以下条件：

a：一致性条件：测试信号f(t)和参考信号g(t)两段信号中包含相同的有效信息，即有效信号的个数是一致的；

b：单调性条件：有效的地震信号具有次序性；

c：边界条件：地震信号两侧端点一致。

在满足上述要求的假设上，差异矩阵可有下式构造：

d(x,y)＝[f(x)-g(y)]²,x＝1,...,mt；y＝1,...,nt

其中，mt和nt分别是两段信号的采样点数。

2、慢度矩阵的构建：

将差异矩阵转换为慢度矩阵：

其中，vmin和vmax为给定的最小和最大速度值，dmax为差异矩阵的最大值。

3、构造变分问题：

将求取两段信号的非线性时移量，转换为求取慢度场中最短旅行时和最短路径问题，把慢度场视为二维图像，对应的变分问题为：

其中，α(x,y)为加权系数，l0和l1是边界值。

4、提取最短路径：

基于费马原理和惠更斯-菲涅尔原理，采用最小路径射线追踪方法提取全局最小旅行时和路径。主要流程如下：

(1)节点设置:

将慢度矩阵分为矩形网格,在矩形网格单元顶点处建立节点，并确定起始节点和终止节点。

(2)节点搜寻方式：

根据信号满足单调性条件，只需在第一象限中，搜寻从当前二次源节点si到相邻节点gj的所有可能节点，搜索方式如图1所示。

(3)旅行时计算：

初始节点设为初始震源节点s0，围绕s0搜索可能的临近计算节点计算从当前震源节点s0到计算节点的加权旅行时。从中选择最小的旅行时记为为对应的节点为g1。

把当前节点g1作为新的二次源s1，重新搜索当前震源的计算节点并计算对应的旅行时，最小旅行时记为对应的计算节点记为g2。从初始震源到当前选择的最新计算节点的总的旅行时为对应的路径为重复上述过程，直到最终获得的计算节点到达终止节点，获得全局最小旅行时及对应的路径：

若获得的计算节点到达差异矩阵边界时，其节点非终点节点。此时沿着差距矩阵边界搜索，直到到达终点节点为止。

(4)更新：

为避免局部极值，交换初始节点和终止节点，重新计算最小旅行时及其对应的路径，根据信号的边界条件，若计算的路径为全局最短路径，需满足以下条件：

若不满足，选择相同路径的边界点作为新的起始和终止点，初始化参数，重新追踪，直到获得最终的全局最小路径。

实现流程如图2所示；

①输入地震数据：在地震信号处理模型中输入两道时移地震数据。

②差异矩阵构造：将数据归一化处理后，形成差异矩阵。

③慢度矩阵构建：将差异矩阵转换为慢度矩阵。

④提取最短路径：基于费马原理和惠更斯-菲涅尔原理，采用最小路径射线追踪方法提取全局最小旅行时和路径。

⑤输出：将获得的最短路径对应的每个采样点的旅行时时移量施加在测试信号上，获得更新后的信号。

下面以两段理论数据为例说明本发明的具体实现流程。然后以xx油田xx地区二维叠前地震数据为目标靶区，应用本方法对该资料进行优化叠加，进行去噪，以验证本方法的效果。该实际地震资料时间长度6000ms，时间采样间隔为4ms，采样点数1500。采用上述方法对该资料进行处理。

处理过程与效果分析：

实施例1：

1：输入两段地震数据(如图3所示)，其中左侧a为参考信号，右侧b为测试信号。

2：构造差异矩阵(如图4a所示),并将其转换为慢度矩阵(如图4b所示)。

3：提取最短路径，(如图5a、图5b所示)，其中从起始节点到终止节点计算的路径与从终止节点到起始节点获得的路径的曲线不一致，提取相同路径中的边界(图5a中a-f六个点)，重新进行最短路径射线追踪，最后获得最短路径如图5b所示。

4：将获得的最短路径对应的每个采样点的旅行时时移量施加在测试信号上，获得更新后的信号(如图6所示)。

分析可以发现，无论是向后时移的地震子波还是向前时移的地震子波，都可以正确的校正到正确的位置上，该实验例证了本方法提出的有效性。

实施例2：

将本方法应用到实际叠前地震资料优化叠加去噪过程中。

1：首先将炮道集划分为几个小区域，在每个小区域内，通过叠加所有道获得参考道。

2：在每个小区域内，将所有地震道分别与参考道计算最短路径，并获得时移扰动量，获得更新后的地震道。其估计过程同实验1流程。

3：最后在每个区域内将更新后的叠加，从而实现叠前地震数据去噪。

其中图7a为处理前的炮道集，图7b为处理后的炮道集。可以看到，原始炮集中含有较多的相干噪音，深层的同相轴被噪音掩盖，无法识别。经过处理之后，噪音得到有效的去除，无论是浅层还是深层，同相轴连续性更好。