基于WOA-VMD的瞬变电磁探测信号降噪方法

戚庭野，卫会汝，冯国瑞，张新军，余传涛，赵德康，杜孙稳

(1. 太原理工大学矿业工程学院，山西太原，030024；2. 山西省绿色采矿工程技术研究中心，山西太原，030024)

随着我国矿产资源的开发利用，地下资源开采后遗留的采空区面积逐年增加。地下采空区坍塌往往会造成巨大经济损失和人员伤亡[1−2]，因此，准确探测并掌握采空区的空间分布对矿山的发展和矿区城市化进程具有重要意义[3]。目前，采空区探测手段众多，其中瞬变电磁法(transient electromagnetic method，TEM)因其勘探深度大、体积效应小、对低阻层分辨率高和对高阻层探测能力强被广泛应用[4−5]。瞬变电磁法通过对接收的二次感应电动势计算反演得到相关的地层信息。因此，接收到的感应电动势信号质量直接关系到解释工作的可靠性。然而，在野外勘探环境中，瞬变电磁信号不可避免地被周围的高压电线、仪器设备、作业车辆以及自身传感器产生的电磁噪声影响，导致很难从接收的信号中提取出真实的感应电动势的衰减特征，目标区域的地质信息得不到有效反映[5]。因此，有必要提出一种滤除信号中的噪声成分的方法，这成为采空区准确预测的前提条件。

当前研究中，研究者根据瞬变电磁衰减曲线的特点，提出了诸多降噪方法，其中，时频分析方法是有效识别和分离瞬变电磁信号噪声的主要手段之一。DAI 等[6−8]根据理论电磁响应的瞬时特性和实测噪声特点，研究了小波变换方法在抑制电磁干扰方面的优势。虽然单独使用小波变换可以显著降低噪声，但细节系数中保留噪声仍需要联合其他方法进行解决。解海军等[9]采用混合小波去噪方法克服了传统阈值和单一阈值去噪方法的缺陷，在提高瞬变电磁信号信噪比的同时保证了数据质量，但没有解决小波基选取的问题。此外，经验模态分解(EMD)方法能自适应对非平稳、非线性信号进行分解，在瞬变电磁信号去噪领域得到广泛应用。

WANG 等[10]利用EMD 方法将瞬变电磁信号自适应分解为一系列的IMF 分量，通过功率检测自适应识别噪声获得有效信号，但没有考虑EMD 存在的端点效应和模态混叠现象。朱通等[11]在EMD算法的基础上，利用多项式拟合延拓来压制信号的端点失真问题，可有效地识别电磁信号中的谐波噪声。LIU 等[12]将改进的EMD 算法即集成经验模态分解算法(EEMD)应用到瞬变电磁信号降噪领域，成功去除了信号中的运动噪声，但忽略了添加的辅助高斯白噪声的分解残留问题。

鉴于 EMD 及其改进算法的不足，DRAGOMIRETSKIY 等[13]提出一种非递归式自适应分解时间序列信号的变分模态分解(variational mode decomposition，VMD)算法。该算法具有坚实的数学理论，鲁棒性强，避免了端点效应和模态混叠问题，但分解结果受分量个数K和惩罚因子α的限制。为避免人为经验选取参数造成的时间损失和效果偏差，本文将MIRJALILI 等[14]提出的鲸鱼优化算法(whale optimization algorithm，WOA)用于VMD 参数的优化。相比传统的优化算法[15]，鲸鱼优化算法不仅需要调节的参数少，寻优精度高，收敛速度快，而且全局寻优能力强，在许多领域得到广泛应用[16]。基于此，本文提出一种基于WOA-VMD 的瞬变电磁信号噪声处理方法，通过仿真模拟和现场实测数据处理，利用该方法对瞬变电磁仿真信号和实测信号进行自适应分解，并且与基于EEMD 方法的分解结果对比，以检验本文所提出的方法的有效性和可行性。

1 降噪算法

1.1 VMD算法

VMD(variational mode decomposition)是一种将信号非递归分解成数个具有准正交性的带限本征模态函数的算法，其基本原理可以表示为一个约束变分问题的求解，构造的约束变分问题表示为

采用乘子交替方向算法交替更新得到K个IMF及其对应的中心频率。

从式(2)可以看出：K和α影响着VMD 的分解性能。若设置的K较小，则信号的多个分量可能同时包含在1 个模态中；若K较大，则会导致1 个分量包含在多个模态中，迭代得到的中心频率也会重叠。对α而言，若α很大，则带宽限制就会很窄，从而导致有用的频率成分被消除；反之，冗余频率成分将会被保留下来。因此，本文提出鲸鱼优化算法来优化确定最佳参数组合(K，α)。

1.2 WOA算法

鲸鱼优化算法是一种基于座头鲸种群迭代进化搜索的元启发式优化算法。通过模拟座头鲸独特的气泡网攻击狩猎行为，实现复杂优化问题寻优，具体优化机理包括局部开发阶段和全局探索阶段2个部分[14]。

1)在局部开发阶段，鲸鱼假设当前的猎物位置为目标位置，在收缩包围机制和螺旋更新机制之间以50%的概率不断靠近猎物获得最优解，其数学模型如下：

2)在全局探索阶段，鲸鱼以其他鲸鱼位置为参考，不断更新自身位置进行随机搜索，其数学模型可表示为

式中：A的绝对值大于1；Xrand为随机选择的鲸鱼位置。

1.3 基于WOA的VMD参数优化

本文提出采用鲸鱼优化算法实现变分模态分解算法的参数优化。该算法中，参数优化的关键是选择适应度函数。排列熵作为一种度量混沌时间序列复杂性的平均熵参数，时间序列复杂度参数计算更简单，抗干扰能力更强，具有更好地鲁棒性的特点，尤其适用于存在动态噪声和观测噪声的信号[17]。基于此，本文选择排列熵作为WOAVMD算法的适应度函数，具体构建步骤如下。

1)给定1个离散时间序列{x(i)，i=1～N}，对该序列进行相空间重构，得到如下重构矩阵：

式中：r=N-(m- 1)τ，r为重构行向量的个数；j为重构矩阵的第j行分量，j=1～r；m为嵌入维数；τ为延迟时间。

2)将重构矩阵中的r个行向量分别按照元素数值进行升序排列，得到各自的符号序列S(q)=(j1,j2,…,jm)。其中，q=1～r，r≤m！；m！为m维相空间映射的符号序列的总数；j1，j1，…，jm表示各元素在原重构分量中的索引号。

3) 计算r种不同符号序列S(q)出现的概率P1，P2，…，Pr，则时间序列{x(i)，i=1～N}的排列熵定义为

基于上述提出的适应度函数，鲸鱼优化算法和变分模态分解算法联合后的信号优化分解的流程图如图1所示，具体步骤如下。

图1 WOA算法优化VMD参数的流程图Fig.1 Flow chart of WOA algorithm to optimize VMD parameters

1)输入信号，设置VMD 算法中K和α的参数范围，初始化WOA模型中的各项参数，包括种群规模、最大迭代次数、空间维度以及初始种群个体。

2)对信号进行VMD分解，利用式(6)计算初始种群中每个个体的适应度。作为该方法的适应度函数，排列熵被用来衡量参数组合的分解效果。熵越小，则时间序列分布越有规律，表示VMD处理得到的IMF 包含更多的有效信息；反之，时间序列越接近随机分布，IMF 中噪声成分更多。因此，当排列熵最小时，对应的参数最优。

3)当|A|≤1时，选择最小排列熵对应的鲸鱼位置作为局部开发的目标值，然后根据p，选择式(3)更新鲸鱼个体的位置；当|A|>1时，随机选择一个鲸鱼的位置，根据式(4)更新鲸鱼个体的位置，保留最优适应度及对应的参数组合。

4)保留更新后的鲸鱼种群位置作为新一轮的初始种群，循环迭代，直到达到所设定的最大迭代次数为止。

5)输出最优鲸鱼个体及对应的适应度。

1.4 相关模态的选择

为了有效识别和分离经VMD分解后的IMF分量中的噪声分量，在对IMF 分量进行时频分析的基础上引入相关系数，进一步判断IMF 与原信号之间的相关程度。相关系数C的定义如下[18]：

式中：ui(t)和x(t)为分别表示VMD 分解后的第i个IMF分量和原始瞬变电磁响应电压信号；E和D表示数学上的期望和方差。

C越大，对应的IMF 分量与原信号相关性越好。本文首先根据式(7)计算各个分量与原始信号的相关系数，然后将相关系数按照分量顺序进行排列，取第1个极值点位置作为有效信号分量和噪声分量的分界点。由于VMD 分解得到1 组从低频到高频的IMF 分量，而EEMD 的频率随着IMF 阶数增加而逐渐减小，因此，对于VMD方法选取极值点前的分量将作为有效分量，EEMD 方法极值点前对应的则为噪声分量，最后将有效分量重构得到去噪后的信号。

2 仿真分析

2.1 仿真信号的建立

为了测试基于WOA-VMD 算法对瞬变电磁含噪信号的降噪性能，进行仿真试验。瞬变电磁现场探测经常受到背景噪声、工频噪声及其谐波干扰的影响。

首先，为了仿真出接近真实工况的含噪瞬变电磁信号，基于Gaver-Stehfest 算法构建电阻率为20 Ω·m 的均匀半空间模型，模拟线圈边长为30 m中心回线接收的瞬变电磁信号理想衰减曲线。

其次，根据背景噪声随机分布的特点和工频噪声及其谐波干扰为固定频率特征，信号采用高斯函数和正弦函数产生3组强度分别为−130，−147和−158 dB的噪声。

最后，将3 组噪声添加到瞬变电磁理想信号中产生信噪比(SNR)分别为−5.858，1.642 1 和5.852 6 dB 的含噪信号中，具体的信号计算模型如下：

式中：x(t)为含噪声干扰的瞬变电磁信号的电压；f(t)表示理想的瞬变电磁信号的电压；g(t)为高斯白噪声；f1为工频干扰固定频率50 Hz；f2和f3分别为二次谐波和三次谐波的固定频率，分别为100 Hz和150 Hz；A1，A2和A3分别为信号振幅的系数。

2.2 WOA-VMD去噪处理

以3 种不同信噪比的仿真信号验证WOAVMD 算法的有效性。当含噪信号的信噪比为−5.858 7 dB 时，仿真信号的时域波形及频谱图如图2 所示。从图2(a)可见：晚期波形出现剧烈震荡，信号的衰减趋势被噪声严重干扰，难以准确解译地层真实地电信息。从图2(b)可见，相应的频域变化曲线在50，100 和150 Hz 处发生异常突起，且高频区被毛刺覆盖，频谱的分布特征被破坏。

图2 瞬变电磁含噪信号时域波形及频谱Fig.2 Time-domain waveform and frequency spectrum of TEM noisy signal

对建立的含噪信号利用WOA-VMD 算法进行降噪处理，使用WOA 对VMD 的参数(K，α)进行寻优，设置种群规模为50，最大迭代次数为30，K的迭代范围为(3，15)，α迭代范围为(1 000，15 000)，最小排列熵迭代曲线如图3 所示。由图3可见：最小排列熵在迭代第9次达到最小值6.426，得到最优参数组合(K，α)=(3，14 886)。

图3 VMD参数的WOA优化结果(仿真信号)Fig.3 WOA optimization results of VMD parameters(simulation signal)

设置优化后VMD的参数，分解仿真信号得到3 个IMF 分量，如图4 所示。由图4 可见：TEM 信号被成功提取到IMF1分量中，工频及其谐波干扰和高斯白噪声则混合在其余分量中。计算IMF 分量与TEM 含噪信号之间的相关系数，结果见表1。

表1 仿真信号的IMF分量相关系数Table 1 Correlation coefficient of IMFs of simulated signal

图4 WOA-VMD分解后IMF分量的时域波形Fig.4 Time domain waveform of IMFs after WOA-VMD decomposition

根据IMF分量的相关系数在IMF2处首次取得极小值0.072 7，进一步确定IMF1为有效分量，重构后的信号如图5所示。由图5可见：降噪后的信号时域波形与原信号波形相符，隐藏于仿真信号中的噪声信息基本被剔除，TEM 信号的特征信息被很好地保留下来。频谱图中固定频率处的异常突起和高频处的毛刺也被消除，信号频域曲线变得更加平滑。但由于TEM 信号为以低频为主的宽频信号，频带的重合导致VMD分解时仿真信号发生部分能量损失。值得注意的是，从图5(a)中可以看出，能量损失对信号分解影响很小，信号重构效果比较理想。

图5 瞬变电磁含噪信号经WOA-VMD降噪后的时域波形及频谱Fig.5 Time-domain waveform and frequency spectrum of TEM noisy signal after WOA-VMD denoising

2.3 不同降噪方法对比

为了验证本文提出的WOA-VMD 算法方法在处理复杂信号的优势，将其与EEMD 方法进行对比。采用EEMD 方法对相同仿真信号进行分解，得到8 个IMF 分量如图6 所示。VMD 只有3 个分量，说明VMD方法的计算效率更高，信息聚集性能更好[19]。从图6(a)中仅能分辨出IMF1 为噪声分量，EEMD 的IMF2～IMF8 分量为无趋势信号，无法区分有效信号和噪声。结合图6(b)中IMF 分量晚期时域波形及相关系数的计算结果(见表1)中首个极值点的位置，可以确定IMF1～IMF3 为噪声分量。从仿真信号中将这些噪声信号剔除获得的降噪后信号如图7所示。从图7(b)可见：去噪后的频域曲线较光滑，噪声成分在一定程度上被去除。然而，观察图7(a)发现仍有部分噪声残留在信号中，信号晚期发生上下波动，主要原因是EEMD在基于极值点定义的局部包络函数进行反复递归分解过程中发生包络估计误差[20−21]，另外，受端点效应影响，信号早期出现解调误差，波形发生轻微凹陷。

图6 EEMD分解后IMF分量的时域波形Fig.6 Time domain waveform of IMFs after EEMD decomposition

图7 瞬变电磁含噪信号经EEMD降噪后的时域波形及频谱Fig.7 Time-domain waveform and frequency spectrum of TEM noisy signal after EEMD denoising

同样地，采用WOA-VMD 和EEMD 方法的仿真信号进行降噪，结果如图8 所示。由图8 可见：虽然2种方法都保留了无干扰信号的衰减趋势，但经EEMD 去噪后的信号不仅在早期出现了畸变，而且在晚期依旧存在波动，显然，WOA-VMD 去噪后的平滑效果优于EEMD方法。

图8 瞬变电磁含噪信号降噪结果Fig.8 Denoising results of TEM noisy signals

2.4 仿真试验结论

为了定量评价WOA-VMD对TEM信号的去噪效果，利用信噪比、均方误差和平均相对误差对晚期数据进行衡量，2种方法的降噪评价指标见表2。其中，信噪比越大，均方差越小，平均相对误差越小，去噪效果越好。对比不同噪声强度情况下VMD 和EEMD 的客观降噪指标发现，VMD 方法的总体效果均优于EEMD 的分解效果，WOAVMD 方法得到的信噪比最大，得到的均方差和平均相对误差最小。结合上述仿真结果，证实本文提出的WOA-VMD 算法在处理瞬变电磁信号方面具有优越性。

表2 WOA-VMD和EEMD的降噪评价指标Table 2 Denoising evaluation index of WOA-VMD and EEMD

3 现场应用

3.1 采空区地质条件

为验证WOA-VMD 降噪算法的实用效果，将该方法用于袁家村铁矿实测瞬变电磁信号处理。袁家村铁矿位于山西省吕梁市北部，地层由浅至深依次为新生代第四系、古近系和新近系，下部富集有残坡积铁矿砾石、奥陶系下统白云岩、寒武系砾岩铁矿、灰岩、页岩交互层以及主要的含铁岩层上太古界吕梁山群。由于矿山前期被无序开采，导致地下存在大规模的采空区，为此采用ATEM-II型瞬变电磁系统探查采空区分布。

3.2 不同方降噪分析

瞬变电磁探测中，接收到的二次衰减曲线信号受作业设备以及高压电力线电磁噪声的影响发生剧烈振荡，晚期信号特征被噪声完全掩盖。本文选取某一条测线作为降噪对象，采用WOAVMD方法和EEMD方法分别对该测线上的每个测点进行降噪处理。其中，测点3 实测信号的VMD参数的WOA 优化迭代过程如图9 所示，在迭代至第5 次时排列熵取得最小值4.958，输出最优参数组合(K，α)=(2，2 044)。设置相应参数对信号进行EEMD分解和VMD分解，结果分别如图10和图11所示。由图10 可知：IMF4～IMF8 为虚假分量，IMF1～IMF3 之间发生模态混叠，通过对分量进行局部放大判定IMF1 和IMF2 为噪声分量。由图11可见：噪声分量的中心频率在50 Hz左右，波形变化与工频噪声相似。分别剔除这些噪声分量，重构后的响应数据衰减曲线如图12 所示，信号晚期的噪声明显降低。EEMD 降噪曲线在早期由于模态混叠造成信号畸变，WOA-VMD 的降噪效果更有优势，基本还原了信号的衰减趋势。

图9 VMD参数的WOA优化结果(实测信号)Fig.9 WOA optimization results of VMD parameters(measured signal)

图10 EEMD分解后IMF分量的时−频图Fig.10 Time−frequency diagram of IMFs after EEMD decomposition

图11 WOA-VMD分解后IMF分量的时−频图Fig.11 Time−frequency diagram of IMFs after WOA-VMD decomposition

图12 瞬变电磁实测信号的降噪效果Fig.12 Denoising effect of measured TEM signal

为更直观获知本文提出方法的有效性和优越性，将测线的11 个测点实测数据分别用EEMD 算法和WOA-VMD 算法进行降噪处理，抽取实测信号和降噪信号的26 个测道的感应电压绘制成感应电动势多测道曲线，并进行对比，结果如图13 所示。由图13(a)可见：信号第2测道数据畸变严重，晚期测道数据由于受噪声干扰相互交错严重，从图中无法正确识别采空区的位置。由图13(b)可见：经EEMD 处理后的测道曲线的信号早期测道数据由于算法无法有效分离不同振荡特性的信息出现失真，引入了更多的干扰信息，严重模糊了信号原有的变化趋势。由图13(c)可见：经WOA-VMD处理后，早期突变与晚期的交叉现象都被去除，信号中的噪声信息得到有效过滤，整体剖面曲线变得更平滑。以上分析表明，WOA-VMD 算法能够有效抑制信号中的噪声，较好地恢复信号中的有效信息。

图13 感应电动势剖面曲线图Fig.13 Profile curve of induced electromotive force

3.3 实测降噪数据反演

视电阻率等值线剖面直观反映探测数据的解译结果，因此，对上述实测数据和降噪后的数据进行反演计算，生成视电阻率剖面图，如图14 所示。结合地质资料分析可知，所选取的测线赋存10号矿体，矿体最高标高为1 780 m，最低标高为1 450 m，延伸140～830 m，倾角70°～80°，钻进中测不到水位，地表亦未见泉水露头，为透水不含水层。根据瞬变电磁接收信号可知，在早期，实测数据的反演结果(见图14(a))在距离20 m 和40 m处由于噪声干扰出现高阻假异常现象，经EEMD处理后得到的视电阻率剖面图(见图14(b))中，在1 480～1 520 m 处出现多处低阻假异常区，40 m 和80 m 处的图像质量因模态混叠而无参考意义，而经WOA-VMD 处理后的剖面图(见图14(c))很好地还原了早期的层状结构。随着勘察深度增加，晚期数据受噪声的影响增大，从实测数据的剖面图可以看出高阻异常区散乱分布于1 420～1 500 m，经WOA-VMD 处理后的剖面图像清晰地反映了高阻异常区呈层状分布于1 460～1 500 m。根据处理结果，在距离60 m 处布设了钻孔位置，结果显示在1 470～1 495 m 范围内出现了掉钻现象，证明了采空区的存在，这验证了本文提出方法的有效性和准确性。

图14 视电阻率等值线剖面Fig.14 Contour profile of apparent resistivity

4 结论

1)提出一种基于鲸鱼优化算法改进变分模态分解的信号降噪方法，该方法很好地解决了瞬变电磁法探测采空区时存在噪声干扰的问题。

2)在不同信噪比条件下，改进的VMD方法能够准确地滤除仿真信号中的工频干扰和随机噪声，更好地恢复整体信号的衰减特征。特别是在定量评价分析中，经WOA-VMD方法去噪后的TEM信号的信噪比分别从−5.858 7，1.642 1和5.852 6 dB提高到−1.251 2，13.546 1和25.101 5 dB，信号质量得到较大改善。

3)WOA-VMD 在瞬变电磁法实测采空区数据降噪中的应用表明，该方法能够有效去除强干扰信号中的环境噪声，提高高阻异常区的分辨率，得到更准确的采空区位置及范围，具有良好现场适用性。