CN108919347A - 基于vmd的地震信号随机噪声压制方法 - Google Patents

Abstract

本发明公开了一种基于vmd的地震信号随机噪声压制方法，该方法包括以下步骤：选择一个时间窗口将原始含噪二维地震信号作傅里叶变换，变换至f‑x域；对每一个频率切片数据进行模态分解；将分解得到的BIMF分量组合得到滤波后的信号；将信号作傅里叶逆变换，变换回t‑x域；重复对下一个时间窗口进行同上操作，地震记录全部处理完毕后，即得到最终的二维地震去噪结果。同时，本发明中将一维变分模态分解拓展至二维复数形式，可用于三维地震去噪处理。基于VMD的地震信号随机噪声压制方法具有优异的噪声压制、幅值保持性能的同时，还具备较高的计算效率，可满足高维大尺度地震数据的处理要求。

Description

基于vmd的地震信号随机噪声压制方法

技术领域

本发明涉及地球物理勘探的地震资料处理技术领域，尤其涉及地震信号中随机噪声压制和有效能量保持处理技术。

技术背景

地震信号随机噪声衰减是地震数据处理领域的热点和难点问题。从观测数据中有效去除随机噪声干扰，提高信号信噪比和分辨率，是正反演计算和地质解释的前提。目前，地震信号随机噪声压制方法主要有：基于滤波理论的方法、基于小波域变换方法、基于矩阵理论的变换方法和基于信号分解理论方法等。近年来，基于信号分解的经验模态分解(Empirical Mode Decomposition，EMD)方法，由于在处理非平稳及非线性数据上具有很高的信噪比，已成为信号分解领域研究的热点。

完备集合经验模态分解(Complete Ensemble EMD，CEEMD)避免了EMD方法中模态混叠的问题，有效改善集合经验模态分解EEMD方法中加入随机白噪声对原始信号造成污染的缺点，但CEEMD方法中依然采用的是递归迭代式筛选分解过程，非平稳地震信号极值点插值和包络计算过程依然耗时过长，在处理多维度和多尺度地震数据时存在一定限制。

为解决这些问题，进一步提升信号分解的精度，本发明提出基于变分模态分解(Variational Modal Decomposition，VMD)的地震信号随机噪声压制处理方法。通过将模态分解过程转变为变分模型的最优解迭代求解过程，以每个分量的估计聚集带宽之和最小为约束，双重对偶上升使得每个分量的频率中心及带宽在频域内进行自适应分离。与EMD的递归“筛选”模式不同，VMD将信号分解转化非递归、变分模型泛函极值求解问题，由变分模型优化目标函数求解过程可知，其实质是多个维纳滤波组的自适应拓展，具有优良的噪声鲁棒性，对于非平稳的时序信号分解更为适用。

发明内容

本发明的目的在于提供一种基于VMD的地震信号随机噪声压制方法，以解决地震勘探任务中采集的地震数据被随机噪声干扰的问题，该方法具有优异的噪声压制、幅值保持性能的同时，还具备较高的计算效率，可满足高维大尺度地震数据的处理要求。

本发明的技术方案为：一种基于vmd的地震信号随机噪声压制方法，其技术方案步骤包括：

S1：地震数据的时域信号转换为频域信号：选择一个时间窗口将原始含噪二维地震信号d(x,t)作傅里叶变换，变换至f-x域；

S2：信号的频段分解：对每一个频段数据进行变分模态分解VMD(VariationalMode Decomposition)，生成带限固有模态函数(Band-Limited Intrinsic ModeFunction)分量，简称BIMF分量；

S3：信号的组合：将VMD分解得到的BIMF分量进行组合，得到新的滤波信号；

S4：地震数据的频域信号转换为时域信号：将S3步骤得的滤波信号作傅里叶逆变换，变换回t-x域；

S5：迭代过程：重复对下一个时间窗口进行同上操作，地震记录全部处理完毕后，即得到最终的二维地震去噪结果。

所述S2步骤的VMD分解过程可转换至变分泛函最优求解过程，其分解包括：

通过在频率域中估计固有模态函数BIMF(Band-Limited Intrinsic ModeFunction)分量频率带宽目标函数，其数学表达式如下：

其中，K为预设分解尺度个数，t是时间变量，u_k即为VMD分解后具备带限性质的BIMF分量，ω_k为对应模态的频率中心，是关于时间t的导数，δ(t)为Dirac(狄利克雷)冲击函数，*为卷积符号，f为原始的频率域实值信号，是平方L2范数。的意义是通过Hilbert变换，将每个模态函数u_k变为解析信号，使实值信号u_k转变为复值，以获得u_k的单边频谱。式(1)中的变分问题使得每个BIMF分量的频谱带宽在其中心频率ω_k附近，且要求BIMF的带宽具备稀疏性。

通过将VMD方法的变分优化框架拓展至复数空间，因此专利要求2中的等式改写为如下形式：

采用增广Lagrange函数来计算其最优解。通过引入二次惩罚因子α和拉格朗日乘法算子λ，将约束性变分问题转换为如下的非约束性变分形式：

其中，二次惩罚因子α是控制数据保真度的均衡参数，用于平衡变分正则项和二次约束项，在含噪声情形时可保证信号重构精度；Lagrange乘子λ(x)可以保证模型约束条件的严格性；然后，通过采用交替方向乘子算法ADMM(Alternate Direction Method ofMultipliers)求解式，具体步骤如下：

MM_1：初始化λ¹，n＝0；

MM_2：n＝n+1，执行主循环；

MM_3：For k＝1:K-1，执行第一个内循环更新u_k：

MM_4：k＝K，结束第一个内循环；

:MM_5：For k＝1:K-1，执行第二个内循环更新ω_k：

:MM_6：k＝K，结束第二个内循环；

MM_7：对于所有ω_k＞0，双重对偶上升，更新λ：

其中，τ表示噪声容限参数；在噪声压制任务中(而不是信号的重构)，将更新参数τ＝0，以得到更好的去噪效果。

MM_8：给定判定精度ε＞0，重复步骤2)-7)，直至满足迭代停止条件：

结束迭代，即得到K个带限BIMF分量。

ADMM的求解过程包含VMD的模态更新和频率中心更新。其中，ω_k频率中心的更新由对应模态的能量谱重心得到，u_k模态更新对应于1/αω²的Wiener滤波器结构，α为白噪声方差，1/ω²表示信号的能量谱为低通形式。参数α控制着Wiener滤波器的宽度，本发明称为保真度均衡参数。增大α值，Wiener滤波器宽度变窄，可以滤除更多的噪声，但也使得其包含更少的真实峰信息，同时算法趋于发散不收敛的几率增加，反之亦然。

BIMF分量的个数K值的取值是VMD算法的关键问题，不同的模态个数会对分解结果产生影响，从而影响最终解的评价。本发明通过计算BIMF分量瞬时频率均值变化的极大值确定最佳的模态分解个数。

本发明的优点：VMD方法的分解过程可转换至变分泛函最优求解过程，与EMD类方法的迭代筛选模式不同，信号分解过程被转移到变分框架，以每个带限BIMF分量的估计带宽之和最小为约束，通过增广Lagrange目标函数将变分问题由约束性变为非约束性，采用交替方向乘子ADMM算法寻求变分泛函的最优解达到信号自适应分解的目的。ADMM中频率中心及带宽交替更新对偶上升，使得两者同时达到最优趋势，一次性得到所有BIMF分量，具有更高的时间效率，每个模态分量在频谱上均具有带限特性，实现信号频带的自适应剖分。BIMF分量的更新过程中具备维纳滤波特性，因而VMD可看作是维纳滤波的多重化和自适应阶推广，即分解出的各组分能量已进行了类似维纳滤波操作，该过程具有坚实的理论基础。一维VMD可拓展至频率域二维复数的广义形式。基于VMD的地震信号随机噪声压制处理方法具有优异的噪声压制、幅值保持性能的同时，还具备较高的计算效率，可满足高维大尺度地震数据的处理要求。

附图说明

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

附图1是一种基于VMD的地震信号随机噪声压制处理方法实现流程图；

附图2是各组分信号及VMD分解BIMF分量；

附图3是附图2合成信号的频谱分布与VMD分解后的3个BIMF分量信号的频谱图(双对数坐标)；

附图4为合成地震记录Sigmoid模型；

附图5为Sigmoid模型的含噪记录及其采用各方法去噪结果；

附图6a为实际含噪地震资料采用Curvelet 3D方法去噪结果；

附图6b为实际含噪地震资料采用LM3D方法去噪结果；

附图6c为实际含噪地震资料采用BM4D方法去噪结果；

附图6d为实际含噪地震资料采用VMD 2D方法去噪结果；

附图7a为实际含噪地震资料Curvelet 3D去噪结果与噪声剖面的局部相似度图；

附图7b为实际含噪地震资料NLM3D去噪结果与噪声剖面的局部相似度图；

附图7c为实际含噪地震资料BM4D去噪结果与噪声剖面的局部相似度图；

附图7d为实际含噪地震资料VMD 2D去噪结果与噪声剖面的局部相似度图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明基于VMD的二维地震信号去噪技术方案的具体步骤如下：

步骤1：选择一个时间窗口将原始含噪二维地震信号d(x,t)作傅里叶变换，变换至f-x域；

步骤2：对每一个频段数据进行变分模态分解。首先，估计频率域中固有模态函数分量频率带宽目标函数；然后，将VMD方法的变分优化框架拓展至复数空间，采用增广Lagrange函数将约束性变分问题转换为非约束性变分形式；最后，采用交替方向乘子算法(ADMM)对来计算其最优解，得到K个带限BIMF分量；

步骤3：将VMD分解得到的BIMF分量进行组合，生成滤波后的信号。

步骤4：将信号作傅里叶逆变换，变换回t-x域。

步骤5：重复对下一个时间窗口进行同上操作，地震记录全部处理完毕后，即得到最终的二维地震去噪结果。

在二维地震信号去噪方法基础上，本发明进一步提出三维地震信号去噪方法。通过二维VMD拓展到复值，拓展后的目标优化函数为：

其中，▽为矢量梯度算子，为二维平面的波数向量。与一维VMD求解类似，通过引入二次惩罚和拉格朗日乘数(增广Lagrange)来重建约束变分框架，由ADMM进行优化求解，其二维VMD的三维地震随机噪声压制算法具体步骤如下：

D_1：选择一个时间窗口将原始含噪地震信号d(m,h,t)作傅里叶变换，变换至f-m-h域；

D_2：对每一个频率切片数据进行二维复数VMD分解；

:D_3：将VMD分解得到的BIMF分量组合得到滤波后的信号；

D_4：将信号作傅里叶逆变换，变换回t-m-h域；

:D_5：重复对下一个时间窗口进行同上操作；

D_6：地震记录全部处理完毕后，即得到最终的三维地震去噪结果。

附图2为合成信号s(t)经VMD分解得到的3个BIMF分量信号，观察图中(a-c)，各个BIMF分量信号与原始组分信号几乎一致，在BIMF分量局部放大图(a'-c')中，仅在信号两端点处出现微小误差。

附图3为原始合成信号的频谱分布与VMD分解后的3个BIMF分量信号的频谱图(双对数坐标)，图中每个重构BIMF模态分量的频谱曲线都有一个最高峰值，即对应其中心频率，与原始信号的期望中心频率2、24和288Hz高度一致，因此，VMD分解得到的3个BIMF分量最后均成功捕捉到了这三个中心频率。

附图4为一个含有多倾斜地层、一个不整合面、一个断层和多个正弦状地层组成的理论合成地震记录Sigmoid模型，含256道，每道256个时间采样。VMD方法中带限BIMF分量K＝4，保真度均衡参数α＝2000，频率中心ω_k初始化采用匹配追踪方法，时间轴窗口长度取64个采样点，空间轴窗口长度取64道，时间轴和空间轴的滑动步长重叠度均设置为50％。

附图5(a-f)显示了地震数据增加噪声后各方法去噪结果，附图5(a’-f’)显示了附图5(a-f)对应的噪声剖面间局部相似度图。附图5(a)为加入20％高斯白噪声后数据，附图5(a’)为纯净信号与噪声间的局部相似度图。由图可知，五种方法均达到一定的去噪效果，但附图5(b)Curvelet方法去噪结果发生一定的变形，在正弦状同相轴周围发生明显伪影现象，断层趋势基本消失，剖面的面貌比较模糊，而附图5(b’)可见Curvelet方法明显有效能量泄露；附图5(c)的NLM去噪结果中仍保留有肉眼可视的噪声信息，对应局部相似度图表现为线性和双曲线有效能量均现泄露，对原图像的结构信息保护不够，这与非局部均值(NonLocal Mean，NLM)算法中在计算图像块相似性时只考虑块的平移特性有关；CEEMD方法相对Curvelet和NLM方法在SNR数值上有些许提升，但附图5(d)去噪结果中仍保留部分噪声信息，断层部分尚可辨识；对应局部相似度图，在正弦状同相轴区域有些许能量损失；附图5(e)BM3D方法的去噪效果较好，但也可见部分低能量有效信号被当做噪声被去除；附图5(f)VMD的去噪结果中信噪比SNR提升均最为明显，优于其余四种方法，合成模型中的倾斜、正弦地层以及断层位置都得到了很好的保留，噪声得到有效压制，视觉表现上原始纯净模型最为接近；对比局部相似度图可知VMD方法有效能量泄露程度最轻微，保真度性能最好。

图6分别显示了四种方法去噪结果，图7显示了对应去噪结果与噪声剖面局部相似度图。对比各图可知，Curvelet 3D方法去噪结果过于平滑，同相轴边缘细节信息被涂抹丢失，导致较难辨别，其对应局部相似度表现较多的有效能量损失；LM3D方法去噪结果中可见局部块状细节消失，可见部分噪声依然分布于去噪结果中，其对应局部相似度出现较多的局部高异常区域；BM4D的去噪表现有所改善，剖面视觉效果表现较好，局部相似度图表现较少的有效信息损失；本文VMD 2D方法的去噪结果表现最优，整个剖面同相轴细节、断裂走向表现清晰，局部相似度图呈低能均匀分布，有效同相轴信息丢失最少。

Claims (3)

1.一种基于VMD的地震信号随机噪声压制方法，其步骤包括：

S1：选择一个时间窗口将原始含噪二维地震信号d(x,t)作傅里叶变换，变换至f-x域；

S2：对每一个频段数据进行变分模态分解，生成带限固有模态函数(BIMF)分量；

S3：将VMD分解得到的BIMF分量进行组合，生成滤波后的信号；

S4：将信号作傅里叶逆变换，变换回t-x域；

S5：重复对下一个时间窗口进行同上操作，地震记录全部处理完毕后，即得到最终的二维地震去噪结果。

2.根据权利要求1所述的基于VMD的地震信号随机噪声压制方法，其特征在于：所述S2步骤的VMD分解过程包括：

首先通过在频率域中估计固有模态函数BIMF分量频率带宽目标函数，其数学表达式如下：

其中，K为预设分解尺度个数，t是时间变量，u_k即为VMD分解后具备带限性质的BIMF分量，ω_k为对应模态的频率中心，是关于时间t的导数，δ(t)为Dirac(狄利克雷)冲击函数，*为卷积符号，f为原始的频率域实值信号，是平方L2范数；的意义是通过Hilbert变换，将每个模态函数u_k变为解析信号，使实值信号u_k转变为复值，以获得u_k的单边频谱；式(1)中的变分问题使得每个BIMF分量的频谱带宽在其中心频率ω_k附近，且要求BIMF的带宽具备稀疏性；

然后，通过将VMD方法的变分优化框架拓展至复数空间，因此专利要求2中的等式改写为如下形式：

采用增广Lagrange函数来计算其最优解；通过引入二次惩罚因子α和拉格朗日乘法算子λ，将约束性变分问题转换为如下的非约束性变分形式：

其中，二次惩罚因子α是控制数据保真度的均衡参数，用于平衡变分正则项和二次约束项，在含噪声情形时可保证信号重构精度；Lagrange乘子λ(x)可以保证模型约束条件的严格性；

最后，通过采用交替方向乘子算法ADMM求解式，具体步骤如下：

MM_1：初始化n＝0；

MM_2：n＝n+1，执行主循环；

MM_3：For k＝1:K-1，执行第一个内循环更新u_k：

MM_4：k＝K，结束第一个内循环；

:MM_5：For k＝1:K-1，执行第二个内循环更新ω_k：

:MM_6：k＝K，结束第二个内循环；

MM_7：对于所有ω_k＞0，双重对偶上升，更新λ：

其中，τ表示噪声容限参数；在噪声压制任务中，将更新参数τ＝0，以得到更好的去噪效果；

MM_8：给定判定精度ε＞0，重复步骤2)-7)，直至满足迭代停止条件：

结束迭代，即得到K个带限BIMF分量；

ADMM的求解过程包含VMD的模态更新和频率中心更新；其中，ω_k频率中心的更新由对应模态的能量谱重心得到，u_k模态更新对应于1/αω²的Wiener滤波器结构，α为白噪声方差，1/ω²表示信号的能量谱为低通形式；参数α控制着Wiener滤波器的宽度，本发明称为保真度均衡参数；增大α值，Wiener滤波器宽度变窄，可以滤除更多的噪声，但也使得其包含更少的真实峰信息，同时算法趋于发散不收敛的几率增加，反之亦然；

BIMF分量的个数K值的取值是VMD算法的关键问题，不同的模态个数会对分解结果产生影响，从而影响最终解的评价；本发明通过计算BIMF分量瞬时频率均值变化的极大值确定最佳的模态分解个数。

3.根据权利要求1所述的所述的基于VMD的地震信号随机噪声压制方法，其特征在于：将二维VMD运用于三维地震数据的随机噪声压制处理；通过二维VMD拓展到复值，拓展后的目标优化函数为：

其中，为矢量梯度算子，为二维平面的波数向量；与一维VMD求解类似，通过引入二次惩罚和拉格朗日乘数(增广Lagrange)来重建约束变分框架，由ADMM进行优化求解，其二维VMD的三维地震随机噪声压制算法具体步骤如下：

D_1：选择一个时间窗口将原始含噪地震信号d(m,h,t)作傅里叶变换，变换至f-m-h域；

D_2：对每一个频率切片数据进行二维复数VMD分解；

:D_3：将VMD分解得到的BIMF分量组合得到滤波后的信号；

D_4：将信号作傅里叶逆变换，变换回t-m-h域；

:D_5：重复对下一个时间窗口进行同上操作；

D_6：地震记录全部处理完毕后，即得到最终的三维地震去噪结果。