0 引言

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

学者们在EMD方法基础上，进一步提出两种经典的改进方法，即集合经验模态分解(Ensemble EMD，EEMD)^[18]和完备集合经验模态分解(Complete EEMD，CEEMD)^[19]。其中，EEMD避免了EMD方法中的模态混叠，而CEEMD改善了EEMD方法中因辅助噪声造成的噪声污染。但是，CEEMD方法中依然采用递归迭代式筛选分解过程，非平稳地震信号极值点插值和包络计算过程耗时过长，在处理多维度和多尺度地震数据时存在一定限制。为了解决这些问题，以进一步提升信号分解的精度，本文提出全新的非线性、非递归的完全自适应时序序列分解算法——基于变分模态分解(Variational Modal Decomposition，VMD)的地震信号分解方法，并应用于二维和三维地震数据随机噪声压制中。

1 基本原理

VMD的基本思路是：在EMD的固有模态函数(Intrinsic Mode Function，IMF)分量基础上，引入更严格的频率域带宽先验约束(Band-Limited Intrinsic Mode Function，BIMF)信号，又称带限分量，解决EMD中模态混叠和高频率信号丢失问题。IMF被定义为在信号极值点个数和过零点个数相等或相差一个时，任一时刻点由局部极值点所确定包络均值为零的信号。

VMD基本原理是：在IMF信号分量基础上，通过引入调幅—调频(AM-FM)信号

$ u_{k}(t)=A_{k}(t) \cos \varphi_{k}(t) $

(1)

式中：A_k(t)和φ_k(t)分别为u_k(t)的瞬时幅值和瞬时相位。φ_k(t)为非减函数，即瞬时角频率

$ \omega_{k}(t)=\frac{\mathrm{d} \varphi_{k}(t)}{\mathrm{d} t} \geqslant 0 $

(2)

与φ_k(t)相比，A_k(t)和ω_k(t)的变化较为缓慢，即在短时范围内可认为是一个幅值和频率不变的谐波信号。将满足上述条件的BIMF赋予具体的物理意义，即VMD将输入信号自适应分解成有限个BIMF分量之和，分解过程满足每个带限BIMF分量的估计聚集带宽之和为最小。与传统EMD方法不同，在获取BIMF分量时，VMD方法的递归迭代筛选模式将信号的分解过程转换至变分问题的求解，在寻找变分模型最优解的过程中完成信号的自适应分解。

2 地震随机噪声压制 2.1 地震信号的域表示

线性地震数据可表示为

$ d(t, x)=w\left(t-\frac{x}{c}\right) $

(3)

式中：t和x为时间和炮检距轴；w为地震子波函数，如Ricker子波；c为波在介质中的传播速度。将同相轴数据沿t(时间轴)作傅里叶变换，将t-x域转换到f-x域，得到原始数据的f-x频谱

$ D(f, x)=\hat{w}(f) \mathrm{e}^{\mathrm{i} 2 \pi f \frac{x}{c}} $

(4)

将式(4)进行离散化，即可得到D(f-x)的离散表示

$ D_{f}(m) \equiv D(f, m \Delta x) \quad m=1,2, \cdots, M $

(5)

式中：Δx表示沿x轴方向的道间距；m为地震采样道数，M为总道数。统计相邻两道可表示为

$ D_{f}(m) \equiv \mathrm{e}^{\mathrm{i} 2 \pi f \frac{\Delta x}{c}} D_{f}(m-1) \quad m=1,2, \cdots, M $

(6)

式(6)是一阶递归微分方程，对应每一个频率切片D_f(m)信号，由一个复调和谐波函数(调幅—调频)组成。因此，通过在t-x域中的无噪声线性数据中找到递归系数e^i2πfΔx/c，使其在f-x域中完全可预测。类似地，t-x域中N个线性事件的叠加相当于f-x域中的N个复数谐波的叠加。

2.2 VMD框架求解

基于上述地震数据域分析，在频率域中估计BIMF分量频率带宽目标函数^[20]

$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_{\left\{ {{u_k}} \right\},\left\{ {{\omega _k}} \right\}} \left\{ {\sum\limits_{k = 1}^K {\left\| {{\partial _t}\left\{ {\left[ {\delta (t) + \frac{{\rm{i}}}{{{\rm{ \mathsf{ π} }}t}}} \right]*{u_k}(t){{\rm{e}}^{ - {\rm{i}}{\omega _k}t}}} \right\}} \right\|_2^2} } \right\},}\\ {且满足\;\sum\limits_{k = 1}^K {{u_k}} = F} \end{array} $

(7)

式中：min(·)表示函数最小化；ω_k为对应模态的瞬时频率；K为预设分解尺度个数；u_k即为VMD分解后具备带限性质的BIMF分量；δ(·)为Dirac冲击函数；*为卷积符号；F为原始的频率域实值信号。$\left[\delta(t)+\frac{\mathrm{i}}{\pi t}\right] * u_{k}(t) $的意义是通过Hilbert变换将每个模态函数u_k变为解析信号，使实值信号u_k转变为复值以获得u_k的单边频谱。式(7)中的变分问题使每个BIMF分量频带在ω_k附近且要求带宽具备稀疏性。

通过将VMD方法的变分优化框架(式(7))拓展至复数空间，而复数空间正弦信号的傅里叶频谱是单边的，因此式(7)改写为

$ \mathop {\min }\limits_{\left\{ {{u_k}} \right\},\left\{ {{\omega _k}} \right\}} \left\{ {\sum\limits_{k = 1}^K {\left\| {{\partial _x}\left[ {{u_k}(x)} \right]{{\rm{e}}^{ - {\rm{i}}{\omega _k}x}}} \right\|_2^2} } \right\},且满足\;\sum\limits_{k = 1}^K {{u_k}} = F $

(8)

为求解式(8)中变分问题，采用增广Lagrange函数^[14]计算其最优解。通过引入二次惩罚因子α和拉格朗日乘法算子λ，将约束性变分问题转换为非约束性变分形式

$ \begin{array}{*{20}{c}} {L\left( {\left\{ {{u_k}} \right\},\left\{ {{\omega _k}} \right\},\lambda } \right): = \alpha \sum\limits_{k = 1}^K {\left\| {{\partial _x}\left[ {{u_k}(x)} \right]{{\rm{e}}^{ - {\rm{i}}{\omega _k}x}}} \right\|_2^2} + }\\ {\left\| {F - \sum\limits_{k = 1}^K {{u_k}} (x)} \right\|_2^2 + \left\langle {\lambda ,F - \sum\limits_{k = 1}^K {{u_k}} (x)} \right\rangle } \end{array} $

(9)

式中：L(·)表示目标函数；“:”是一个函数的标记；〈·, ·〉表示內积；二次惩罚因子α是控制数据保真度的均衡参数，用于平衡变分正则项和二次约束项，在含噪声情形时可保证信号重构精度；λ可以保证模型约束条件的严格性。

采用交替方向乘子算法(Alternate Direction Method of Multipliers，ADMM)^[22]求解式(9)，具体步骤如下：

(1) 初始化u_k⁽¹⁾、ω_k⁽¹⁾、λ⁽¹⁾，n=0；

(2) n=n+1，执行主循环；

(3) For k=1:K-1，执行第一个内循环更新u_k

$ u_k^{(n + 1)} \leftarrow \arg \mathop {\min }\limits_{{u_k}} L\left( {\left\{ {u_{i < k}^{(n + 1)}} \right\},\left\{ {u_{i \ge k}^{(n)}} \right\},\left\{ {\omega _i^{(n)}} \right\},{\lambda ^{(n)}}} \right) $

(10)

式中arg min表示函数L(·)达到最小时u_kⁿ⁺¹的值。

(4) k=K，结束第一个内循环；

(5) For k=1:K-1，执行第二个内循环更新ω_k

$ \omega _k^{(n + 1)} \leftarrow \arg \mathop {\min }\limits_{{\omega _k}} L\left( {\left\{ {u_i^{(n + 1)}} \right\},\left\{ {\omega _{i < k}^{(n + 1)}} \right\},\left\{ {\omega _{i \ge k}^{(n)}} \right\},{\lambda ^{(n)}}} \right) $

(11)

(6) k=K，结束第二个内循环；

(7) 对于所有ω_k＞0, k∈[1, K]，双重对偶上升，更新λ

$ \lambda^{(n+1)} \leftarrow \lambda^{(n)}+\tau\left(f-\sum\limits_{k=1}^{K} u_{k}^{(n+1)}\right) $

(12)

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

(8) 给定判定精度ε＞0, k∈[1, K]，重复步骤(2)~(7)，直至满足迭代停止条件

$ \sum\limits_{k=1}^{K} \frac{\left\|u_{k}^{(n+1)}-u_{k}^{(n)}\right\|_{2}^{2}}{\left\|u_{k}^{(n)}\right\|_{2}^{2}}<\varepsilon $

(13)

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

从上述步骤可知，ADMM的求解过程包含VMD的模态和频率中心更新。其中，ω_k的更新由对应模态的能量谱重新得到，u_k模态更新对应于1/βω²的Wiener滤波器结构(β为白噪声方差，1/ω²表示信号的能量谱为低通形式)。参数α控制着Wiener滤波器的宽度，本文称为保真度均衡参数。增大α值，Wiener滤波器宽度变窄，可以滤除更多的噪声，但也使得BIMF分量包含更少的真实信号峰值信息，同时算法趋于发散不收敛的几率增加，反之亦然。式(7)的增广Lagrange函数形式与Tikhonov正则一致，具有保真项和正则项。其中，VMD方法中分解的带限BIMF分量具有明显的稀疏特性，因而正则项可用稀疏先验。α与正则参数类似，用于控制保真项和正则项的权重，使BIMF模态分量的更新公式(式(10)~式(12))具备Wiener滤波特性，这正是VMD方法具有较强抗噪性能的关键因素。

2.3 基于VMD的二维地震信号去噪方法

根据上述理论，基于VMD方法的二维地震信号随机噪声压制处理算法流程如下：

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

(2) 对每一个频率切片数据进行复数VMD分解；

(3) 将VMD分解得到的BIMF分量组合得到滤波后的信号；

(4) 将信号作傅里叶逆变换回t-x域；

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

2.4 基于VMD的三维地震信号去噪方法

在三维地震数据中，平面波可表示为

$ d(t,\mathit{\boldsymbol{x}}) = w\left( {t - \frac{{\langle \mathit{\boldsymbol{x}},\mathit{\boldsymbol{n}}\rangle }}{c}} \right) $

(14)

式中：x=(m, h)为中心点和炮检距坐标；n为波的传播单位方向。三维地震模型的频率切片D_f是由沿n方向的一个f-m-h域复谐波模式组成。将d(t, x)沿时间轴上作傅里叶变换得到f-m-h谱，每个频率切片数据即为D_f。相应地t-m-h域中p个平面波的叠加与f-m-h域中p个复谐波叠加等效。本文将二维VMD拓展到复值，拓展后的目标优化函数为

$ \mathop {\min }\limits_{\left\{ {{u_k}} \right\},\left\{ {{v_k}} \right\}} \left\{ {\sum\limits_{k = 1}^K {\left\| {\nabla \left[ {{u_k}(\mathit{\boldsymbol{x}}){{\rm{e}}^{ - {\rm{i}}\left( {{\mathit{\boldsymbol{v}}_k},\mathit{\boldsymbol{x}}} \right)}}} \right]} \right\|_2^2} } \right\},且满足\;\sum\limits_{k = 1}^K {{u_k}} = F $

(15)

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

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

(2) 对每一个频率切片数据进行二维复数VMD分解；

(3) 组合VMD分解的BIMF分量得到滤波后的信号；

(4) 将信号作傅里叶逆变换回t-m-h域；

(5) 对下一个时间窗口重复以上操作；

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

3 算例 3.1 VMD参数优选

VMD变分模型目标函数是非凸的，算法的收敛性与参数的设置有密切联系，有三个重要参数，即带限BIMF分量个数K、ω_k和α。其中，K依据BIMF分量瞬时频率均值的极大值点确定；ω_k采用匹配追踪算法(Matching Pursuits，MP)^[23]确定；α与原始信号的噪声水平有关，用来控制保真项和正则项的权重，本文采用“L”曲线法则确定α。为验证VMD方法的有效性，本文将三个不同中心频率余弦信号累加构成合成信号s(t)

$ \left\{ \begin{array}{l} {s_i}(t) = \cos \left( {2{\rm{ \mathsf{ π} }}{f_i}t} \right)\;\;\;\;{f_i} = 2,24,288(i = 1,2,3)\\ s(t) = \sum\limits_{i = 1}^3 {{s_i}} (t) \end{array} \right. $

(16)

图 1a、图 1c、图 1e显示了合成信号s(t)经VMD分解得到的3个BIMF分量信号，图 1b、图 1d、图 1f分别为上述相应的BIMF分量信号局部放大结果。从分解结果可以看出，各个BIMF分量信号与原始组分信号几乎一致，仅在信号两端点处出现微小误差。

图 2为原始合成信号的频谱分布与VMD分解后的3个BIMF分量信号的频谱图(双对数坐标)，每个重构BIMF模态分量的频谱都有一个最高峰值，对应其中心频率，与原始信号的期望中心频率2、24、288Hz高度一致(非等间距横坐标)，VMD分解得到的3个BIMF分量均能成功捕捉对应中心频率。

图 3为分解各模态信号的频率中心随迭代次数的变化曲线(半对数坐标)。由图可见，图 3a迭代更新过程较慢，在150次后才找到最优的中心频率，而图 3b经3次迭代即可获得最佳频率中心分布。

图 4是一组由四个线性事件组成的合成地震数据、加噪数据及其对应的f-x频谱图、f-k频谱图。对图 4b中第40个采样点切片含噪数据作VMD分解，结果如图 5所示。设置参数α=2000，K=4，以匹配追踪算法初始化ω_k分解得到四个近似谐波信号，各个BIMF分量的频谱与原始信号的频谱峰值高度一致。

3.2 二维地震去噪实例

本文采用信噪比(SNR)、局部相似度(Local Similarity，LS)^[24]、结构相似度(Structural Simila-rity Index，SSIM)^[25]评价算法的有效性。局部相似度定义为

$ \mathrm{LS}=\sqrt{c_{1}^{\mathrm{T}} c_{2}} $

(17)

式中：T表示转置；c₁、c₂通过最小二乘问题最小化得到

$ {c_1} = \arg \mathop {\min }\limits_{{c_1}} \left\| {\mathit{\boldsymbol{a}} - \mathit{\boldsymbol{B}}{c_1}} \right\|_2^2 $

(18)

$ {c_2} = \arg \mathop {\min }\limits_{{c_2}} \left\| {\mathit{\boldsymbol{b}} - \mathit{\boldsymbol{A}}{c_2}} \right\|_2^2 $

(19)

式中：a和b为矩阵；A、B分别为a、b构成的对角算子。LS值变化范围为[0, 1]，LS值越大表示去噪结果与噪声剖面的局部相似程度越高，有效能量泄漏越严重，算法的幅值保持能力越差。

结构相似度定义为

$ \operatorname{SSIM}(x, y)=\frac{\left(2 \mu_{x} \mu_{y}+C_{1}\right)\left(2 \sigma_{x y}+C_{2}\right)}{\left(\mu_{x}^{2}+\mu_{y}^{2}+C_{1}\right)\left(\sigma_{x}^{2}+\sigma_{y}^{2}+C_{2}\right)} $

(20)

式中：μ_x、μ_y分别表示图像x、y的均值；σ_x、σ_y分别表示图像x、y的标准差；σ_x、σ_y分别表示图像x、y的方差；σ_xy代表图像x和y的协方差；C₁、C₂为常数，C₁=(k₁L)²，C₂=(k₂L)²，L为像素动态范围，一般k₁=0.01、k₂=0.03、L=255；SSIM值变化范围为[0, 1]，值越大表示两个图像的结构相似程度越高，去噪性能越强。

图 6表示一个含有多层倾斜地层、一个不整合面、一条断层和多套正弦起伏形态地层组成的理论合成地震记录Sigmoid模型^[26]，含256道，每道有256个时间采样。为验证VMD方法的有效性，本文将VMD与Curvelet(曲波)、NLM(非局部均值)、BM 3D(三维块匹配)、CEEMD四种方法进行对比。

在本实验中，VMD方法中K=4，α=2000，ω_k初始化采用匹配追踪方法，时间轴窗口长度取64个采样点，空间轴窗口长度取64道，时间轴和空间轴的滑动步长重叠度均设置为50%；CEEMD方法集总次数I取100；Curvelet、NLM和BM3D三种方法采用工具包中的默认参数。算法运行平台参数为MATLAB 2017a Parallel Pool×8、Win 10×64、i7-7700K CPU 4.2GHz。

图 7a~图 7f左图为地震数据加噪后用不同方法去噪的结果，右图为相应的噪声剖面局部相似度。由图可知，五种方法均达到一定的去噪效果，但Curvelet方法去噪结果有一定的变形，在正弦状起伏的强反射同相轴周围发生明显伪影现象，断层现象不明显，弱反射同相轴连续性差，对应的局部相似度图表明有效能量泄漏明显(图 7b)；NLM去噪结果中仍保留有肉眼可视的噪声信息，对应的局部相似度图表明线性和双曲线有效能量均出现泄漏，对原图像的结构信息保护不够，这与NLM算法中在计算图像块相似性时只考虑块的平移特性有关(图 7c)；相对Curvelet和NLM方法，CEEMD方法在SNR数值上有些提升，但去噪结果中仍保留部分噪声信息，断层部分尚可辨识，对应的局部相似度图在正弦状起伏的强反射同相轴有些许能量损失(图 7d)；BM 3D方法的去噪效果较好，但也可见部分低能量有效信号被当做噪声被去除(图 7e)；VMD的去噪结果中SNR提升最为明显，优于其余四种方法，合成模型中的倾斜、正弦状起伏地层以及断层位置都得到了很好的保留，噪声得到有效压制，视觉表现上与原始无噪模型最为接近。对比局部相似度图可知，VMD方法有效能量泄漏程度最轻微，保真度性能最好(图 7f)。

图 8给出了五种方法的第175道信号去噪结果对比，可以看出VMD去噪结果与无噪信号最为接近。因此，VMD方法在去除强噪声的同时，可最大程度的减少有效地震能量的损失。

Curvelet、NLM、BM 3D、CEEMD和VMD五种方法处理耗时分别为0.724、11.291、181.039、2.207和2.513s，Curvelet去噪处理效率最高，BM3D处理效率最低。

3.3 三维地震去噪实例

处理效率是考虑算法是否适用的重要指标，在三维甚至五维地震资料的随机噪声压制处理中，尤其需考虑计算耗时。本文将Curvelet 3D、NLM 3D、BM 4D和二维VMD(VMD 2D)四种方法运用于三维实际含噪地震数据处理。取某地随机噪声分布较多的三维叠后实际资料，数据尺度为221×271×752(依次表示x、y和时间三个方向的采样点数)，采样间隔分别为20m和4ms。取实际数据x、y方向的第100个切片，如图 9所示。

图 10为四种方法去噪结果，图 11为相应的噪声剖面局部相似度图，去噪结果与噪声剖面的局部相似程度越高，表明有效能量泄漏越严重。对比各图可知，Curvelet 3D方法去噪结果过于平滑，同相轴边缘细节信息模糊，较难辨别，对应局部相似度表现较多的有效能量损失；NLM 3D方法去噪结果中局部块状细节消失，可见部分噪声依然分布于去噪结果中，对应局部相似度出现较多的局部高异常区域；BM 4D方法的去噪效果有所改善，剖面视觉效果较好，局部相似度表现为有效信息损失少；本文二维VMD 2D方法的去噪结果最优，剖面上同相轴细节、断裂走向清晰，局部相似度呈低能均匀分布，有效同相轴信息丢失最少。

Sigmoid模型和实际资料的去噪结果数值统计如表 1所示，可见Curvelet方法的SNR数值提升最小，但时间效率最高；CEEMD方法处理耗时最长；与CEEMD递归式模态分解方法不同，VMD在变分框架约束下一次性求解得到所有带限BIMF分量，在SNR和SSIM指标上均达最优，去噪结果与原始无噪信号在细节上最为接近。由于在匹配追踪优化下快速定位到频率中心，VMD方法时间效率高于NLM、CEEMD和BM 3D方法。三维实际地震资料的去噪结果LS同样验证其有效性。因此，VMD方法具备优异的噪声压制、幅值保持性能的同时，还有具有较高的计算效率，可以满足实时处理要求。

4 结束语

本文在EMD信号分解理论基础上，提出了在频率域内复数VMD方法。通过将信号分解过程转移到变分框架内，引入增广Lagrange函数迭代求解各分量的频率中心及带宽参数，生成窄带约束模态分量，使VMD方法在噪声压制和有效能量保持间寻求最优平衡。实验结果表明，在二维地震数据的随机噪声压制处理方面，与CEEMD、Curvelet、NLM和BM 3D方法相比，VMD方法在SNR、SSIM和视觉效果方面均保持最优趋势，处理后同相轴细节变得清晰连续，微弱信号被有效保留；在三维地震数据的随机噪声压制处理方面，与Curvelet 3D、NLM 3D和BM 4D方法相比，VMD方法具备优异的噪声压制性能且计算效率高。