基于变分模态分解参数优化的地震随机噪声去除方法

2019-10-16李钟晓

北京化工大学学报(自然科学版) 2019年5期

徐智唐刚* 刘伟李钟晓

(1.北京化工大学机电工程学院，北京 100029； 2.青岛大学电子信息学院，青岛 266071)

引言

由于人工、环境和地质条件等因素的影响，采集的地震数据常含有大量随机噪声，这会严重影响地震数据的后续处理和解释，因此有效压制随机噪声对提高地震数据的信噪比具有重要意义。

小波分析因其多分辨的特点，在随机噪声压制方面具有较好的效果[1-3]。陈晓玉[4]使用小波阈值去噪降低了地震数据中的随机噪声。崔少华等[5-6]研究了分解层数、重构个数、去噪阈值以及阈值处理方式对地震数据去噪效果的影响。但受基函数固定、多分辨率恒定等限制，小波分析缺乏自适应性。

经验模态分解(empirical mode decomposition, EMD)可以自适应压制地震数据中的噪声，提升信噪比[7-8]，但EMD方法存在模态混叠和端点效应等不足，导致其分解结果不稳定、不唯一。集合经验模态分解[9]和完备集合经验模态分解[10]能在一定程度上缓解模态混叠问题，但同时也极大增加了计算量，且端点效应问题仍然存在。为解决这一问题，何元等[11]和康佳星等[12]使用变分模态分解(variational mode decomposition, VMD)[13]进行地震数据处理，结果都表明VMD方法抗噪性更好，运算效率更高，且去噪效果仅依赖于分解参数的选取。然而，目前VMD参数大都依赖经验选取，结果具有较大随机性和偶然性。本文提出使用估计信噪比作为评价参数，自适应地优化VMD的参数分解和有效模态选取过程，以改善地震数据去噪效果。

1 理论方法

1.1 变分模态分解

变分模态分解是将信号f(t)分解为有限个具有限定带宽的模态函数uk(t)，每个模态函数围绕在该模态的中心频率ωk附近，变分问题的解为最小带宽之和。

1.1.1变分模态分解模型构建

首先对每个模态分量进行Hilbert变换。

(1)

式(1)中，t为时间，uk(t)为第k个模态分量。

其次对每个模态预设一个中心频率，将各模态得到的解析信号频谱通过移频方式移至基带上。

(2)

式(2)中，ωk为第k个模态预设的中心频率。

最后计算解调信号梯度的范数，估计出各频带的带宽。

(3)

1.1.2变分模态分解模型求解

首先引入二次惩罚因子α和拉格朗日乘法算子λ(t)，将约束变分问题转化为无约束变分问题，其表达式为

(t)]e-jωkt2+f(t)-∑kuk(t)2+λ(t),f(t)-

∑kuk(t)

(4)

然后利用乘子交替方向算法(ADMM)求取无约束变分问题的鞍点uk和ωk，即为该变分问题的解。uk和ωk的具体更新步骤如下。

1) 更新uk将式(4)转化为等价最小化问题

(5)

利用二范数下的Parseval/Plancherel傅里叶等距，通过式(6)在傅里叶域求解式(5)

(6)

使用ω-ωk代替式(6)中ω，结果如式(7)

(7)

利用重构保真项中实数信号的埃尔米特对称性将式(7)写成非负频率上的半空间积分

(8)

利用式(9)求解式(8)的二次优化问题

(9)

2) 更新ωk由于重构保真项中不含有ωk，所以式(4)可以简化为

(10)

将式(10)转换到傅里叶域

(11)

对公式(11)求解如下

(12)

最后迭代拉格朗日算子λ

(13)

1.1.3VMD算法步骤

综合1.1.1和1.1.2两节可得VMD算法的步骤如下：

(2)根据公式(9)和公式(12)迭代计算uk和ωk；

(3)根据公式(13)迭代计算λ；

1.2 敏感参数优化

1.2.1敏感参数

分数高斯噪声(fractional Gaussian noise，FGN)是由分数布朗运动(fractional Brown motion，FBM)表示的时间离散序列。本文采用FGN的数值模拟实验对VMD等效滤波特性进行分析[14]。设xH(m)(m=…，-1,0,1，…)为FGN序列，自相关函数pH(l)=E{xH(m)xH(m+l)}满足式(14)[15]

(14)

式中，H为影响FGN特性的Hurst指数，σ为FGN的标准差。0

使用式(14)对VMD方法进行实验可知，VMD分解等价于一系列的带通滤波器，带通滤波器主要由分解参数K进行控制，K越大，带通滤波器的个数越多，带宽越小。由VMD实验结果还可以得出，VMD地震数据去噪严重依赖于模态分量个数K的大小。由于K值主要根据经验设定，具有偶然性和随机性；并且与经验模态分解类似，分解后的模态中有若干个模态可分为有效模态(包含有效信息)和噪声模态(包含噪声信息)，说明对有效模态的选取也会影响去噪效果。故K的选择及有效模态的选取成为VMD地震数据降噪处理的关键步骤。

1.2.2估计信噪比

通常采用信噪比来定量评价地震数据的去噪效果，但实际地震数据由于没有真实信号可作为参考而无法直接计算其信噪比，所以需要对地震数据进行信噪比估计。刘洋等[16]比较研究了叠加法、时域奇异值分解法和频域奇异值分解法对地震数据信噪比的估计效果，结果表明频域奇异值分解法在相邻道的地震数据间同时存在振幅变化和相位变化时都可以很好地估计信噪比，对于不同地震数据的估计能力优于其他两种方法。故本文选取频域奇异值分解法作为信噪比估计的方法。

在同时考虑振幅变化和相位变化时，地震信号可以表示为

(15)

(16)

任意两道的频域互相关系数为

(17)

(18)

对矩阵R进行奇异值分解后可以写成R=UΛVT，其中U是RRT的特征值对应的特征向量矩阵，Λ=diag(σ1,σ2,…,σN)是RRT特征值的非负平方根按递减顺序组成的对角矩阵；V是RTR的特征值对应的特征向量矩阵。可以证明，矩阵R的奇异值为

(19)

(20)

(21)

故信号的能量可由式(22)计算得出

(22)

估计信噪比可由式(23)计算得出

(23)

式中，Es为信号能量，En为噪声能量。

1.2.3基于估计信噪比的VMD敏感参数优化

本文建立了基于估计信噪比的变分模态分解敏感参数优化方法，具体步骤如下。

(1)设置参数。设对照分解模态数目K0=0，对照含噪信号估计信噪比为xSNR0。

(2)确定K的范围。VMD的分量个数至少应为2，故K的下限取为2；EMD分解时产生模态混叠且存在残余分量，其上限值应为EMD分解模态分量个数KEMD。故K的取值范围为[2,KEMD]。

(3)在取值范围内选取K值。使用VMD方法分解地震数据得到K个模态，使用频域奇异值分解方法求取模态估计信噪比xSNRm。xSNRm>0时该模态中的有效信号能量强于噪声能量，故将该模态作为有效模态；xSNRm<0时该模态中的噪声能量强于有效信号能量，故将该模态作为噪声模态。将有效模态叠加得到重构信号，使用频域奇异值分解方法求取重构信号估计信噪比xSNRr，并与xSNR0进行比较。若xSNRr>xSNR0，则xSNR0=xSNRr，K0=K；否则K0和xSNR0均保持不变。

(4)遍历取值范围内的K值后迭代更新K0。

(5)使用K0对地震数据进行VMD分解重构，重构结果即为基于估计信噪比的变分模态分解参数优化方法的去噪结果。

基于VMD参数优化的地震数据去噪流程如图1所示。

图1 VMD参数优化流程图Fig.1 VMD parameter optimization flow chart

1.3 去噪质量评价指标

为了对比各种方法对模型数据去噪处理的效果，以信噪比、互相关系数和均方根误差[17]作为地震去噪质量的评价指标。

信噪比(signal- noise ratio，SNR)是真实信号与随机噪声能量比值的函数(下文称为实际信噪比)，其计算公式为

(24)

式(24)中，Ps和Pn分别代表信号和噪声的能量，f(i)为原始信号，(i)为重构信号，N′为信号长度。xSNR值越大越好。

均方根误差(root-mean-square error，RMSE)代表重构信号同原始信号之间的偏差，其计算公式为

(25)

xRMSE越小越好。

互相关系数(mutual relationship，MR)表示原始信号与重构信号线性相关的程度，其计算公式为

(26)

式(25)中，σ(f(t))为原始信号标准差，σ((t)为重构信号标准差，cov (f(t),(t))为原始信号和重构信号的协方差。xMR越接近1，相关程度越高，去噪效果越好。

2 数值实验

2.1 估计信噪比

使用合成的二维地震剖面(图2)对信噪比估计方法进行实验验证。分别加入强度为[1 dBW,10 dBW]的随机噪声，得到估计信噪比与实际信噪比如表1所示。

图2 合成的二维地震剖面Fig.2 Synthesized 2D seismic profile intensity

噪声强度/dBW实际信噪比/dB估计信噪比/dB15.84445.539524.84624.748133.74193.872342.83602.951751.84792.003360.84290.97377-0.1575-0.03358-1.1568-1.13919-2.1484-2.219110-3.1591-3.3638

由表1可以看出，估计信噪比随噪声强度的增加而减小，能有效反映地震数据的噪声程度；并且估计信噪比与实际信噪比误差很小，可以将估计信噪比作为有效模态选取和参数优化的评价参数。

2.2 仿真模型数据

2.2.1二维剖面

为了验证本文方法的有效性，首先进行了仿真模型数据实验。该模型数据包含1 000道，每道1 000个采样点，采样间隔2 ms。原始数据如图3(a)所示，加入随机噪声后数据如图3(b)所示。

图3 二维数据各方法去噪结果对比Fig.3 Comparison of the results of methods for obtaining 2D data

EMD分解模态数目为12个，故参数搜索范围为[2,12]。使用本文方法进行搜索筛选得到最优分解模态数K=4，再使用K=4对地震数据进行VMD分解重构，结果如图3(c)所示，其差剖面如图3(d)所示。比较图3(c)、(b)可以看出，处理后噪声明显减小，有效信号得以凸显；由图3(d)看出，差剖面中仅包含噪声，不含同相轴信息，说明本文方法可以在有效压制地震信号噪声的同时完整保留有效地质信息。

为了验证本文方法的优越性，同时使用小波阈值去噪方法及EMD去噪方法进行对比实验。小波阈值去噪使用symmlet 6小波进行信号分解，分解层数为3[6],利用“3σ”方法选取阈值后通过硬阈值方法处理小波系数[5]，再对处理后小波系数进行重构，结果如图3(e)所示，其差剖面如图3(f)所示。EMD去噪通过对信号的EMD分解选取包含有效信息的模态进行重构[8]，得到结果如图3(g)所示，其差剖面如图3(h)所示。可以看出，小波分析和EMD均有一定的去噪效果，但仍有部分噪声无法去除。

综上可得，参数优化的VMD相较小波变换及EMD方法能够在进一步压制噪声的同时完整保留有效信息，提高地震数据的信噪比。

图4 单道仿真模型数据各方法处理结果对比Fig.4 Comparison of experimental single channel simulation model data processing results

为避免主观判断去噪效果可能出现的偏差，通过计算信噪比、互相关系数和均方根误差对去噪效果进行量化分析，结果如表2所示。由表2可以得出，与含噪地震剖面相比，EMD方法信噪比提升7.886 0 db，小波方法信噪比提升8.518 7 db，而本文方法信噪比提升16.674 1 db；EMD方法均方根误差降低了1.885 4，占总误差的59.66%，小波方法均方根误差降低了1.974 9，占总误差的62.50%，而本文方法均方根误差降低了2.696 6，占总误差的85.34%；EMD方法互相关系数提升了0.293 3，占含噪数据互相关系数误差的68.40%，小波方法互相关系数提升了0.309 1，占含噪数据互相关系数误差的72.08%，而本文方法相关系数提升了0.407 1，占含噪数据互相关系数误差的94.94%。评价指标定量分析结果表明，小波分析、EMD和参数优化的VMD方法均能压制地震信号中的噪声，但参数优化的VMD方法压制噪声效果最好。

表2 二维数据处理结果的评价指标

2.2.2单道数据

从2.2.1节仿真模型数据中抽取第100道数据对本文方法进行验证。该道数据包含1 000个采样点，采样间隔2 ms。原始信号如图4(a)所示，加入随机噪声后数据如图4(b)所示。利用参数优化的VMD方法去噪后的单道数据如图4(c)所示，去除的噪声如图4(d)所示。由图4(a)、(b)对比可以看出，含噪信号噪声很大，有效信号被噪声淹没；由图4(c)、(b)对比可知，经参数优化VMD方法处理后的干扰噪声得到有效压制，子波信号被较好地恢复，从而显著增强了信号的信噪比和分辨率。

利用小波变换和EMD进行对比实验，小波去噪的单道处理结果和去除的噪声分别如图4(e)、(f)所示，EMD去噪的单道处理结果和去除的噪声分别如图4(g)、(h)所示。由图4各方法处理结果对比可以看出，参数优化的VMD方法比小波分析和EMD方法的子波恢复能力更强，压制噪声的效果更明显。

图5 实际地震数据处理结果Fig.5 Actual seismic data processing results

由表3可以得出，EMD方法信噪比提升8.494 5 db，小波阈值去噪方法信噪比提升9.273 0 db，而本文方法信噪比提升16.988 1 db；EMD方法均方根误差降低了1.967 4，占总误差的61.25%，小波阈值去噪方法均方根误差降低了2.074 1，占总误差的64.57%，而本文方法均方根误差降低了2.743 9，占总误差的85.43%；EMD方法互相关系数提升了0.291 7，占含噪数据互相关系数误差的71.16%，小波阈值去噪方法相关系数提升了0.312 7，占含噪数据互相关系数误差的76.29%，而本文方法相关系数提升了0.39，占含噪数据互相关系数误差的95.15%。以上评价指标计算结果表明参数优化VMD的去噪效果明显优于小波去噪和EMD去噪。

表3 单道模型数据处理结果的评价参数

2.3 实际数据

本文所用实际地震数据包含800道，每道数据包含800个采样点，采样点时间间隔4 ms，原始实际地震数据如图5(a)所示。

设EMD分解的K值为10，利用本文方法在[2,10]范围内对K进行筛选搜索，求得最优K值为5，使用K=5对地震数据进行VMD分解重构如图5(b)所示，差剖面如图5(c)所示。对比图5(a)、(b)可以看出，实际地震数据中的噪声被去除，有效信号得以保留，去噪效果明显，显著提高了信号的信噪比。

通过小波阈值去噪方法处理的结果如图5(d)所示，差剖面如图5(e)所示；通过EMD方法处理的结果如图5(f)所示，差剖面如图5(g)所示。

由图5(b)、(d)和(f)对比可以看出，小波分析和EMD分解有一定的去噪效果，但是在0～200道处噪声干扰大，噪声压制效果不佳；参数优化的VMD方法的整体噪声压制效果明显优于上述两种方法。通过图5(c)可以看出，参数优化的VMD方法去除的噪声中不含有效地质信息，同样证明本文方法可以在保留有效地震信息的同时显著压制噪声，是一种较优的地震信号噪声去除方法。