引言

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

小波分析因其多分辨的特点，在随机噪声压制方面具有较好的效果^[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)分解为有限个具有限定带宽的模态函数u_k(t)，每个模态函数围绕在该模态的中心频率ω_k附近，变分问题的解为最小带宽之和。

1.1.1 变分模态分解模型构建

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

$ \left[ {\delta (t) + \frac{{\rm{j}}}{{{\rm{ \mathsf{ π} }}t}}} \right] * {u_k}(t) $

(1)

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

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

$ \left[ {\left[ {\delta (t) + \frac{{\rm{j}}}{{{\rm{ \mathsf{ π} }}t}}} \right] * {u_k}(t)} \right]{{\rm{e}}^{ - {\rm{j}}{\omega _k}t}} $

(2)

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

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

$ \left\{ {\begin{array}{*{20}{l}} {\mathop {\min }\limits_{\left\{ {{u_k}} \right\},\left\{ {{\omega _k}} \right\}} \left\{ {\sum\limits_k {{{\left\| {{\partial _t}\left[ {\left[ {\delta (t) + \frac{{\rm{j}}}{{{\rm{ \mathsf{ π} }}t}}} \right] * {u_k}(t)} \right]{{\rm{e}}^{ - {\rm{j}}{\omega _k}t}}} \right\|}^2}} } \right\}}\\ {{\rm{s}}{\rm{.t}}{\rm{. }}\sum\limits_k {{u_k}} = f(t)} \end{array}} \right. $

(3)

式(3)即为所构造的变分模型。式中，‖·‖²为二范数，{u_k}等同于{u₁, …, u_K}，{ω_k}等同于{ω₁, …, ω_K}，$\sum\limits_k {} $等同于$\sum\limits_{k = 1}^k {} $。

1.1.2 变分模态分解模型求解

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

$ \begin{array}{l} L\left( {\left\{ {{u_k}} \right\},\left\{ {{\omega _k}} \right\},\lambda } \right) = \alpha \sum\limits_k {\left\| {{\partial _t}\left[ {\left[ {\delta (t) + \frac{j}{{\rm{ \mathsf{ π} } t}}} \right]{u_k}} \right.} \right.} \\ {\left. {\left. {(t)} \right]{{\rm{e}}^{ - {\rm{j}}{\omega _k}t}}} \right\|^2} + {\left\| {f(t) - \sum\limits_k {{u_k}} (t)} \right\|^2} + \left\langle {\lambda (t),f(t) - } \right.\\ \left. {\sum\limits_k {{u_k}} (t)} \right\rangle \end{array} $

(4)

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

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

$ \begin{array}{l} u_k^{n + 1} = \mathop {\arg \min }\limits_{{u_k}} \left\{ {\alpha \left\| {{\partial _t}\left[ {\left[ {\delta (t) + \frac{{\rm{j}}}{{\rm{ \mathsf{ π} } t}}} \right]{u_k}(t)} \right]} \right.} \right.\\ \left. {{{\left. {{{\rm{e}}^{ - {\rm{j}}{\omega _k}t}}} \right\|}^2} + {{\left\| {f(t) - \sum\limits_i {{u_i}} (t) + \frac{{\lambda (t)}}{2}} \right\|}^2}} \right\} \end{array} $

(5)

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

$ \begin{array}{l} \hat u_k^{n + 1} = \mathop {\arg \min }\limits_{{u_k},{{\hat u}_k}} \left\{ {\alpha \left\| {{\rm{j}}\omega \left[ {\left( {1 + {\rm sgn} \left( {\omega + {\omega _k}} \right)} \right){{\hat u}_k}} \right.} \right.} \right.\\ {\left. {\left( {\omega + {\omega _k}} \right)} \right\|^2} + {\left\| {\hat f(\omega ) - \sum\limits_i {{{\hat u}_i}} (\omega ) + \frac{{\hat \lambda (\omega )}}{2}} \right\|^2} \end{array} $

(6)

使用ω－ω_k代替式(6)中ω，结果如式(7)

$ \begin{array}{l} \hat u_k^{n + 1} = \mathop {\arg \min }\limits_{{u_k},{{\hat u}_k}} \left\{ {\alpha \left\| {{\rm{j}}\left( {\omega - {\omega _k}} \right)\left[ {(1 + {\rm sgn} (\omega )){{\hat u}_k}} \right.} \right.} \right.\\ \left. {{{\left. {(\omega )} \right\|}^2} + {{\left\| {\hat f(\omega ) - \sum\limits_i {{{\hat u}_i}} (\omega ) + \frac{{\hat \lambda (\omega )}}{2}} \right\|}^2}} \right\} \end{array} $

(7)

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

$ \begin{array}{l} \hat u_k^{n + 1} = \mathop {\arg \min }\limits_{{u_k},{{\hat u}_k}} \left\{ {\int_0^\infty 4 \alpha {{\left( {\omega - {\omega _k}} \right)}^2}{{\left| {{{\hat u}_k}(\omega )} \right|}^2} + 2} \right.\\ \left. {{{\left| {\hat f(\omega ) - \sum\limits_i {{{\hat u}_i}} (\omega ) + \frac{{\hat \lambda (\omega )}}{2}} \right|}^2}{\rm{d}}\omega } \right\} \end{array} $

(8)

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

$ \hat u_k^{n + 1}(\omega ) = \frac{{\hat f(\omega ) - \sum\limits_{i \ne k} {{{\hat u}_i}} (\omega ) + \frac{{\hat \lambda (\omega )}}{2}}}{{1 + 2\alpha {{\left( {\omega - {\omega _k}} \right)}^2}}} $

(9)

$\hat u_k^{n + 1}$即为u_k的更新结果。

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

$ \begin{array}{l} \omega _k^{n + 1} = \mathop {\arg \min }\limits_{{\omega _k}} \left\{ {\left\| {{\partial _t}\left[ {\left[ {\delta (t) + \frac{j}{{\rm{ \mathsf{ π} } t}}} \right]{u_k}(t)} \right]} \right.} \right.\\ \left. {{{\left. {{{\rm{e}}^{ - {\rm{j}}{\omega _k}t}}} \right\|}^2}} \right\} \end{array} $

(10)

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

$ \omega _k^{n + 1} = \mathop {\arg \min }\limits_{{\omega _k}} \left\{ {\int_0^\infty {{{\left( {\omega - {\omega _k}} \right)}^2}} {{\left| {{{\hat u}_k}(\omega )} \right|}^2}{\rm{d}}\omega } \right\} $

(11)

对公式(11)求解如下

$ \omega _k^{n + 1} = \frac{{\int_0^\infty \omega {{\hat u}_k}{{(\omega )}^2}{\rm{d}}\omega }}{{\int_0^\infty {{{\hat u}_k}} {{(\omega )}^2}{\rm{d}}\omega }} $

(12)

最后迭代拉格朗日算子λ

$ {\hat \lambda ^{n + 1}}(\omega ) = {\hat \lambda ^n}(\omega ) + \tau \left[ {\hat f(\omega ) - \sum\limits_k {\hat u_k^{n + 1}} (\omega )} \right] $

(13)

w_kⁿ⁺¹即为w_k的更新结果。

1.1.3 VMD算法步骤

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

(1) 初始化{$\hat u_k^{ 1}$}、{ω_k¹}、{${\hat \lambda ^1}$}和n；

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

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

(4) 给定误差x_tol，若$\sum\limits_k {\frac{{{{\left\| {\hat u_k^{n + 1} - \hat u_k^n} \right\|}^2}}}{{{{\left\| {\hat u_k^n} \right\|}^2}}}} < {x_{{\rm{tol}}}}$，则停止迭代，否则重复步骤(2)、(3)。

1.2 敏感参数优化 1.2.1 敏感参数

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

$ {p_H}(l) = \frac{{{\sigma ^2}}}{2}\left( {{{\left| {l - 1} \right|}^{2H}} - 2{{\left| l \right|}^{2H}} + {{\left| {l + 1} \right|}^{2H}}} \right) $

(14)

式中，H为影响FGN特性的Hurst指数，σ为FGN的标准差。0 < H < 0.5时FGN序列为负相关；0.5 < H < 1时FGN序列为正相关；H=0.5时FGN序列不相关，且等同于高斯噪声。

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

1.2.2 估计信噪比

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

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

$ {\mathit{\boldsymbol{\tilde s}}_l} = {\alpha _l}{{\rm{e}}^{{\rm{j}}{\theta _l}}}\mathit{\boldsymbol{\tilde \omega }} $

(15)

式(15)中，${{\mathit{\boldsymbol{\tilde s}}}_l}$为地震信号s_l的频域列向量，${\mathit{\boldsymbol{\tilde \omega }}}$为子波列向量ω的频域列向量，α_l为第l道振幅系数，θ_l为第l道相位角。故频域内各道数据向量可写为

$ {\mathit{\boldsymbol{\tilde d}}_l} = {\mathit{\boldsymbol{\tilde s}}_l} + {\mathit{\boldsymbol{\tilde n}}_l} = {\alpha _l}{{\rm{e}}^{{\rm{j}}{\theta _l}}}\mathit{\boldsymbol{\tilde \omega }} + {\mathit{\boldsymbol{\tilde n}}_l} $

(16)

式(16)中，${{\mathit{\boldsymbol{\tilde d}}}_l}$为含噪地震信号d_l的频域列向量；${{\mathit{\boldsymbol{\tilde n}}}_l}$为噪声n_l的频域列向量。

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

$ \begin{array}{l} {\mathit{\boldsymbol{R}}_{{{\mathit{\boldsymbol{\tilde d}}}_i}{{\mathit{\boldsymbol{\tilde d}}}_l}}} = {{\mathit{\boldsymbol{\tilde d}}}_i} * {{\mathit{\boldsymbol{\tilde d}}}_l} = {\alpha _i}{\alpha _l}{{\rm{e}}^{{\rm{j}}\left( {{\theta _l} - {\theta _i}} \right)}}\mathit{\boldsymbol{\tilde \omega }} * \mathit{\boldsymbol{\tilde \omega }} + {\alpha _i}{{\rm{e}}^{ - {\rm{j}}{\theta _i}}}\mathit{\boldsymbol{\tilde \omega }}*{{\mathit{\boldsymbol{\tilde n}}}_l} + \\ {\alpha _l}{{\rm{e}}^{{\rm{j}}\theta }}\mathit{\boldsymbol{\tilde \omega }} * {{\mathit{\boldsymbol{\tilde n}}}_i} + {{\mathit{\boldsymbol{\tilde n}}}_i} * {{\mathit{\boldsymbol{\tilde n}}}_l} \end{array} $

(17)

假设在频域内噪声与信号互不相关，各道之间噪声互不相干，且各道噪声方差均为σ²。记‖ω‖²= $\mathit{\boldsymbol{\tilde \omega }} * \mathit{\boldsymbol{\tilde \omega }}$，则

$ {\mathit{\boldsymbol{R}}_{{{\mathit{\boldsymbol{\tilde d}}}_i},{{\mathit{\boldsymbol{\tilde d}}}_l}}} = {\alpha _i}{\alpha _l}{{\rm{e}}^{{\rm{j}}\left( {{\theta _l} - {\theta _i}} \right)}}{\left\| \mathit{\boldsymbol{\omega }} \right\|^2} + {\delta _{il}}{\sigma ^2} $

(18)

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

$ \sigma _1^2 = \sum\limits_{l = 1}^N {\alpha _l^2} {\left\| \mathit{\boldsymbol{\omega }} \right\|^2} + {\sigma ^2} = \sum\limits_{l = 1}^N {{{\left\| {{{\mathit{\boldsymbol{\tilde s}}}_l}} \right\|}^2}} + {\sigma ^2} $

(19)

$ \sigma _2^2 = \sigma _3^2 = \cdots = \sigma _N^2 = {\sigma ^2} $

(20)

式中，N为奇异值个数。由式(19)、(20)可以看出，信号能量集中在σ₁中，噪声的能量平均分布在所有奇异值中。通过σ₂至σ_N的噪声方差σ²可以求得噪声平均能量E_n为

$ {\bar E_{\rm{n}}} = {\sigma ^2} = \frac{1}{{N - 1}}\sum\limits_{l = 2}^N {\sigma _l^2} $

(21)

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

$ {E_{\rm{s}}} = \sum\limits_{l = 1}^N {{{\left\| {{{\mathit{\boldsymbol{\tilde s}}}_l}} \right\|}^2}} = \sigma _1^2 - {\bar E_{\rm{n}}} = \sigma _1^2 - \frac{1}{{N - 1}}\sum\limits_{l = 2}^N {\sigma _l^2} $

(22)

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

$ {x_{{\rm{ESNR}}}} = \frac{{{E_s}}}{{{E_{\rm{n}}}}} = \frac{{{E_{\rm{s}}}}}{{N{{\bar E}_{\rm{n}}}}} = \frac{{\sigma _1^2 - \frac{1}{{N - 1}}\sum\limits_{l = 2}^N {\sigma _l^2} }}{{\frac{N}{{N - 1}}\sum\limits_{l = 2}^N {\sigma _l^2} }} $

(23)

式中，E_s为信号能量，E_n为噪声能量。

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

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

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

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

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

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

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

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

1.3 去噪质量评价指标

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

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

$ {x_{{\rm{SNR}}}} = 10\lg \left( {\frac{{{P_{\rm{s}}}}}{{{P_{\rm{n}}}}}} \right) = 10{\rm{lg}}\left( {\frac{{\frac{1}{{{N^\prime }}}\sum\limits_{i = 1}^{{N^\prime }} {{f^2}} (i)}}{{\frac{1}{{{N^\prime }}}\sum\limits_{i = 1}^{{N^\prime }} {{{\left[ {f\left( i \right) - \hat f\left( i \right)} \right]}^2}} }}} \right) $

(24)

式(24)中，P_s和P_n分别代表信号和噪声的能量，f(i)为原始信号，$\hat f\left( i \right)$为重构信号，N′为信号长度。x_SNR值越大越好。

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

$ {x_{{\rm{RMSE}}}} = \sqrt {\frac{1}{{N'}}\sum\limits_{i = 1}^{N'} {{{\left[ {f\left( i \right) - \hat f\left( i \right)} \right]}^2}} } $

(25)

x_RMSE越小越好。

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

$ {x_{{\rm{MR}}}} = \frac{{ {\rm cov} (f(t),\hat f(t))}}{{\sigma (f(t))\sigma (\hat f(t))}} $

(26)

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

2 数值实验 2.1 估计信噪比

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

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

2.2 仿真模型数据 2.2.1 二维剖面

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

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方法能够在进一步压制噪声的同时完整保留有效信息，提高地震数据的信噪比。

为避免主观判断去噪效果可能出现的偏差，通过计算信噪比、互相关系数和均方根误差对去噪效果进行量化分析，结果如表 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.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方法的子波恢复能力更强，压制噪声的效果更明显。

由表 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去噪。

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方法去除的噪声中不含有效地质信息，同样证明本文方法可以在保留有效地震信息的同时显著压制噪声，是一种较优的地震信号噪声去除方法。

3 结论

VMD分解可以有效去除地震数据中的噪声, 相较小波分析和EMD方法，VMD压制噪声效果更明显。利用仿真模型数据和实际数据对多种方法进行实验验证的结果表明，基于估计信噪比的参数优化的VMD方法去噪效果良好，能够有效克服人为选取参数的主观随机性，提升了VMD分解压制噪声效果的稳定性。

致谢: 在本文完成过程中，哈尔滨工业大学马坚伟教授和于四伟博士提供了有益讨论和帮助，东方地球物理公司提供了实际地震数据，在此一并表示感谢。