搞信号处理的人应该都有这种体会:实测数据永远没有教科书里那么干净。机械设备加速度传感器采回来的振动信号、心电监护仪上的生理信号、桥梁结构的动态响应信号,哪一个不是叠了一堆环境噪声、工频干扰和随机脉冲?传统的FIR/IIR滤波器应对平稳信号还行,一旦碰上非线性、非平稳的实测信号,就明显力不从心——该滤的噪声没滤干净,该留的特征被削了尖。所以这些年基于自适应分解的降噪思路越来越流行,其中ICEEMDAN(改进完全自适应噪声集合经验模态分解)配合小波降噪重构,再加上排列熵(PE)做评估,属于实战里打磨出来非常实用的一个组合。这篇文章就把这套方法从原理到MATLAB实现完整拆一遍,包含参数怎么设、代码怎么写、实际用起来有哪些坑,给正在做振动分析、故障诊断或者生理信号处理的朋友一个可参考的落地版本。
1. 方法原理与整体设计思路
1.1 从EMD到ICEEMDAN:每一代都在解决什么问题
ICEEMDAN的全称是Improved Complete Ensemble Empirical Mode Decomposition with Adaptive Noise,也就是“改进的完全自适应噪声集合经验模态分解”。它的整个家族谱系值得先理清楚,因为你只有知道前面几代方法卡在哪,才能真正理解ICEEMDAN改了什么、为什么这样改。
第一代是经典EMD(经验模态分解)。它的核心思想是把一个复杂信号自适应地分解成若干个本征模态函数(IMF),从高频到低频逐层剥离。优点是完全数据驱动,不需要预设基函数,对非线性非平稳信号尤其友好。但痛点也显而易见:模态混叠。所谓模态混叠,就是本来应该属于不同时间尺度的成分被搅进了同一个IMF里,比如一段信号里同时存在连续振荡和一个间歇性脉冲,EMD会把这些成分混在一起,后续分析直接失真。
第二代EEMD(集合经验模态分解)的思路是“以噪治噪”:往原始信号里多次加入白噪声,分别做EMD,然后把各次结果平均。白噪声的加入会打散信号中的间歇成分,让不同尺度的特征被分到不同的IMF里。这个方法有效缓解了模态混叠,但引入了新问题——白噪声加入后,分解结果里会残留噪声成分,重构后的信号不再是完全干净的,而且计算量成倍上升。
第三代CEEMDAN(完全自适应噪声集合经验模态分解)做了一次重要改进:每一层模态分解时都注入自适应白噪声,并且是在“残差”里加,而不是像EEMD那样直接在原始信号上加。这一改,重构误差大幅降低,计算效率也有提升。但CEEMDAN依然有隐患:某些IMF里还是可能残留少量噪声分量,而且分解过程中噪声是从早期阶段一直传递到后期阶段的。
第四代就是ICEEMDAN。它的关键改变在于噪声添加方式——不再把噪声直接加到残差信号上,而是生成一种特殊的噪声序列,并让算法先在信号中“估算”出噪声分量的水平,再决定每一层该去除多少。这种方式有效抑制了残余噪声的泄漏,也让模态分解结果更干净、更稳定。翻译成大白话就是:前几代算法是“把噪声撒进去再想办法去掉”,ICEEMDAN是“先搞清楚信号里有多少噪声,再精准剔除”。所以,论分解精度和重构保真度,ICEEMDAN在目前这个家族里是最好用的。
1.2 排列熵(PE)为什么能当噪声评估的裁判
把信号分解成IMF之后,接下来的核心问题就变成:哪些IMF是有效信号,哪些IMF是噪声主导?总不能用眼睛一个个看频谱、拍脑门做决定。工程上需要的是一个量化指标,能够自动给出一个判定阈值。排列熵(Permutation Entropy,PE)在这个场景里非常合适。
排列熵的基本逻辑并不复杂:把一个时间序列按特定延迟嵌入维度重构,得到一组组“局部排列模式”,然后统计这些模式出现的概率,算一个香农熵。信号的规律性越强、确定性越高,排列模式就越集中,熵值越低;信号越杂乱无章、越接近随机噪声,各种排列模式都冒出来,熵值就越高。实测下来,白噪声和各种随机噪声的PE值普遍在0.8以上,纯净的周期信号或调幅信号的PE值通常在0.2到0.5之间。这种巨大的数值差异,就是因为噪声的“无秩序”和有效信号的“有秩序”在排列模式分布上表现得极为分明。
相比其他熵指标,PE有几个实打实的优势。第一,计算量小。样本熵和近似熵需要对每个向量计算距离,复杂度高,信号一长就跑得很痛苦;PE只做排序和统计,速度快得多,实测上万点数据也不会觉得卡。第二,参数少,只有嵌入维数m和时间延迟τ两个参数,调参非常省心。第三,抗干扰能力强,对信号的幅值变化不敏感,也就是说噪声大也不影响它的相对判断能力。在工程现场,这意味着你可以很稳定地用同一个阈值去处理不同幅值等级的数据。
1.3 方法组合的完整处理链路
这套方法把ICEEMDAN、PE和小波降噪串成了一条完整流水线,整体流程非常清晰:
第一步,对原始含噪信号执行ICEEMDAN分解,得到一组从高频到低频排列的IMF分量和一个残差项。这一步的关键是确定噪声标准差倍数和集合次数,直接影响分解效果。
第二步,对每个IMF计算排列熵PE值。这一步相当于给每个分量做了一个“噪声体检”,根据PE值与预设阈值的大小关系,把分量划分为信号主导分量和噪声主导分量两类。
第三步,对噪声主导的IMF分量执行小波阈值降噪。小波降噪本身也是一套成熟方法,这里不是对全频段瞎处理,而是只针对已经被判定为噪声主导的分量做温处理,信号主导的分量原样保留。这样既去了噪声,又不伤有效信号的特征尖峰和瞬态成分。
第四步,将未处理的分量与降噪后的分量一起重构,得到最终的去噪信号,然后就可以送入后续的特征提取、故障诊断等环节。
这种组合方式的巧妙之处在于取长补短。ICEEMDAN解决的是非线性非平稳信号的分解问题,但它自己不会判断分解结果里哪些组分是噪声;PE解决的是判断问题,但光有判断没有处理手段;小波降噪解决的是处理问题,但如果你不经筛选地对所有高频分量一刀切,就会连有效的高频瞬态特征一起干掉。三个阶段各管一段,按顺序接力,整个系统就完整了。
2. 排列熵(PE)的评估机制与参数选择
2.1 PE到底怎么算:数学过程通俗版
排列熵的原理说起来数学味很浓,但落地理解其实很直观。给你一个时间序列x(1), x(2), ..., x(N),要算它的PE,核心就四步。
第一步,相空间重构。选定嵌入维数m和时间延迟τ,把序列切成若干个长度为m的子序列。以m=3、τ=1为例,第一个子序列是[x(1), x(2), x(3)],第二个是[x(2), x(3), x(4)],依此类推。每个子序列里有3个数值。
第二步,排列映射。看每个子序列内部的数值大小顺序。比如子序列[4, 2, 6],从小到大排列是2、4、6,对应的序号分别是2、1、3,那这个子序列的排列模式就是(2, 1, 3)。所有子序列都按照这套规则映射成一个排列模式。
第三步,统计概率。对所有的排列模式进行统计,计算每种模式出现的概率。当m=3时,理论上最多有3! = 6种排列模式。
第四步,算香农熵并归一化。把各模式概率代入香农熵公式H = -∑p·log(p),再除以log(m!),得到归一化排列熵,数值范围在0到1之间。这个归一化处理非常关键,因为它让不同长度的时间序列之间有了可比性。
整个过程用MATLAB写出来也就二十来行,后面我会给完整代码。理解这一步的关键在于:白噪声这种完全随机的序列,各种排列模式出现的概率几乎均等,所以熵值接近1;正弦波这种规则信号,排列模式高度集中,所以熵值很低。
2.2 嵌入维数和时间延迟怎么选
PE的两个参数——嵌入维数m和时间延迟τ,虽然看着简单,实际上对计算结果影响很大,这里直接说结论。
嵌入维数m的取值范围一般是3到7。m太小,比如m=2,排列模式只有2种,区分能力太差,信号稍微复杂一点就分不开;m太大,比如m>10,排列模式数量暴增到m!种,要准确估计每种模式的概率就需要非常长的信号,样本量不够的话,统计偏差会很大。常规推荐m=5或m=6,对大多数工程信号都能取得稳定且敏感的结果。如果信号本身周期性强、特征明显,m=4也够用;如果是处理高频噪声严重的信号,可以适当提高到m=7,但不要更高。
时间延迟τ的选择相对灵活,常用取值为1到3。τ=1是最常见的选择,适合采样率适中、信号时间尺度比较短的场景;如果信号的采样率很高,相邻采样点之间变化很小,可以取τ=2或τ=3,相当于从时间轴上隔几个点再看变化,避免因为过度采样导致排列模式被“稀释”。有个实用经验是:先取τ=1跑一遍,如果结果显示噪声IMF和信号IMF的PE值拉不开差距,再尝试增大τ,直到两类分量的区分度最大。
还有一个容易被忽略的点:嵌入维数m和信号长度N之间要满足约N >> m!的条件。比如m=6时,排列模式理论上有720种,信号长度如果只有500点,统计出来的概率就不够可靠。所以信号较短时,建议把m降低到4或5,别硬套大维数。
2.3 PE阈值如何设置:给一个可复用的标定思路
PE阈值的选择直接决定哪些IMF被处理、哪些被保留,是整个方法里最需要经验的部分。我不建议直接照抄某个文献里的固定数字,因为不同信号的幅值特性、采样率、噪声类型都不同,最好是现场标定一次阈值。
标定思路其实很简单。第一步,从原始信号里截取一段明显没有有效信号、只有噪声的片段,比如机械停机阶段的传感器输出、心电信号里已知的噪声段,用这段噪声跑一遍和主信号相同的ICEEMDAN分解和PE计算,得到噪声主导分量的PE参考范围。第二步,对包含明显特征成分的信号段做同样的分解和PE计算,得到信号主导分量的PE参考范围。第三步,取两组数值的中间区域作为分界阈值。
从我的实测经验来看,大多数场景下PE阈值落在0.55到0.75之间。纯噪声分量的PE值通常在0.85以上,而有效信号分量的PE值通常在0.4以下,中间留出的空档其实挺大。所以实际操作中,可以先设0.6试跑,然后观察重构信号的质量——如果重构后信号里还有明显的高频毛刺,说明阈值设高了,把噪声分量当成了信号分量;如果重构后信号变得过于平滑、包络特征被削平,说明阈值设低了,把有效分量误伤了。根据结果往反方向调,一般调两三轮就能找到适合自己的值。
3. MATLAB实现与参数配置全解
3.1 运行环境与依赖工具
在动手写代码之前,先说清楚跑这套流程需要什么环境。
MATLAB版本建议R2019b以上,因为新版本对内置信号处理工具箱里的函数优化得更好,读代码、调试代码的体验也更好。需要安装Signal Processing Toolbox,因为后续涉及小波变换的函数wavedec、wrcoef、wthresh等都在这个工具箱里,少一个都跑不起来。如果版本太老,这部分函数缺失,就得自己手写小波变换,工作量会大很多。
至于ICEEMDAN的代码,它不在MATLAB官方工具箱里,需要用论文作者公开提供的算法包。自己的代码目录里把ICEEMDAN相关函数放进去,通过addpath引入路径即可。为了不涉及版权问题,我下面给出自己封装的调用框架和核心流程说明,你拿到任何一版标准实现都可以直接套用这个框架。
3.2 ICEEMDAN分解的MATLAB代码实现
ICEEMDAN的完整实现代码很长,这里给出的是封装调用框架,重点在于理解输入的每个参数含义:
% ===================== ICEEMDAN 分解调用示例 ===================== % 输入: % x 原始含噪信号(列向量,长度N) % Nstd 噪声标准差相对输入信号的比值,常用0.2 % NE 集合平均次数,常用100~500,信号越复杂取越大 % MaxIter 最大筛分迭代次数,常用200~500 % 输出: % IMFs N*m矩阵,m为分解得到的IMF个数,每列是一个IMF分量 % R 残差项,代表信号的整体趋势 % % 这里假设你已下载标准ICEEMDAN函数包,并添加到了MATLAB路径中 x = x(:); % 统一为列向量 Nstd = 0.2; NE = 200; MaxIter = 500; [IMFs, R] = iceemdan(x, Nstd, NE, MaxIter); % 分解完成,查看各IMF的波形 figure; for k = 1:size(IMFs, 2) subplot(size(IMFs, 2) + 1, 1, k); plot(IMFs(:, k)); title(['IMF ', num2str(k)]); end subplot(size(IMFs, 2) + 1, 1, size(IMFs, 2) + 1); plot(R); title('残差R');重点说参数。Nstd是噪声标准差倍数,取值一般为0.1到0.3。取值太小,额外引入的噪声不足以驱动模态分离,分解效果趋近于原始EMD;取值太大,会引入过于明显的噪声痕迹,虽然之后可以靠小波降噪兜底,但整体重构误差会变大。NE是集合次数,它决定计算精度与耗时的平衡——NE太小,统计稳定性差,可能出现分解结果不稳定;NE太大,计算时间成倍增长。我平时处理几万点数据用200次就够,十几万点数据降到100次也能跑,再大就建议先切段处理。
3.3 PE计算与IMF筛选的代码实现
排列熵的计算代码不长,这里给出完整的可运行版本。这个函数我自己一直在用,也调过几次,稳定性和计算速度都经过了验证。
function pe = calc_permutation_entropy(x, m, tau) % 计算归一化排列熵 % 输入: % x 一维信号(列向量) % m 嵌入维数,常用4~6 % tau 时间延迟,常用1~3 % 输出: % pe 归一化排列熵,范围0~1 x = x(:); N = length(x); n = N - (m - 1) * tau; % 重构子序列个数 if n < 100 error('信号长度太短,排列熵统计结果不可靠,请减小m或tau'); end % 第一步:相空间重构,构建时间延迟矩阵 seqMat = zeros(n, m); for i = 1:n seqMat(i, :) = x(i : tau : i + (m - 1) * tau)'; end % 第二步:对每个子序列排序,得到排列模式编码 [~, order] = sort(seqMat, 2); % 每行升序排序,记录原位置 % order的每一行就是该子序列的排列模式,例如[1 3 2] % 第三步:统计唯一排列模式的频数并计算概率 [~, ~, ic] = unique(order, 'rows'); p = accumarray(ic, 1) / n; % 第四步:计算香农熵并归一化 h = -sum(p .* log(p)); pe = h / log(factorial(m)); % 除以最大熵 log(m!) end调用这个函数来筛选IMF,原来需要对每个IMF都算PE,再根据阈值判断。批量筛选的套路见下面这段代码:
% 批量计算各IMF的排列熵 m = 5; tau = 1; pe_values = zeros(size(IMFs, 2), 1); for k = 1:size(IMFs, 2) pe_values(k) = calc_permutation_entropy(IMFs(:, k), m, tau); end % 设置PE阈值,高于阈值判定为噪声主导分量 pe_threshold = 0.60; noise_imf_idx = find(pe_values > pe_threshold); signal_imf_idx = find(pe_values <= pe_threshold); disp('各IMF的PE值:'); disp(pe_values); disp(['判定为噪声主导的IMF编号:', num2str(noise_imf_idx')]); disp(['判定为信号主导的IMF编号:', num2str(signal_imf_idx')]);这里再强调一下m的取值。上面对每个IMF单独算PE时,需要注意该IMF的有效长度可能比原始信号短一些,如果m=7而IMF长度只有2000点,7!=5040种排列模式,2000个样本点根本喂不饱统计需求。所以我的建议是:批量计算PE时,先用信号长度反推允许的最大m。经验公式就是信号点数至少是m!的3到5倍。比如IMF长度5000点,那m=6(720种模式)有条件用;长度只有800点,老老实实用m=4或m=5。
3.4 小波降噪与重构的代码实现
将噪声主导的IMF筛选出来之后,需要对它们执行小波阈值降噪。这里给一个完整的函数,可以直接套用。核心逻辑是:对每一个待处理IMF,先做小波分解,再用软阈值处理细节系数,最后重构。
function imf_denoised = wavelet_denoise(imf, wname, level, scale) % 小波软阈值降噪 % 输入: % imf 待降噪分量(列向量) % wname 小波基名称,如'sym8'、'db4' % level 小波分解层数 % scale 阈值缩放系数,0.6~1.0,默认0.8 % 输出: % imf_denoised 降噪后的分量 imf = imf(:); N = length(imf); % 第一层小波分解 [C, L] = wavedec(imf, level, wname); % 提取各层细节系数 detail_coeffs = cell(level, 1); for j = 1:level detail_coeffs{j} = detcoef(C, L, j); end % 估计噪声标准差:用第一层细节系数的中位数 sigma = median(abs(detail_coeffs{1})) / 0.6745; % 通用阈值:sigma * sqrt(2 * log(N)) thr = sigma * sqrt(2 * log(N)) * scale; % 对每层细节系数做软阈值处理 for j = 1:level detail_coeffs{j} = wthresh(detail_coeffs{j}, 's', thr); end % 重构信号 % 先把阈值处理后的细节系数重新放回小波系数结构C中 C_denoised = C; app_coeffs = appcoef(C, L, wname, level); % 保留近似系数不变 C_denoised = wthcoef('d', C, L, 1:level, thr); % 直接对所有细节层做软阈值 % 或者采用逐层阈值替代方式 % 这里改用逐层手动替换方式,更灵活 recon_start = 0; for j = 1:level idx_start = sum(L(1:end-j)) + 1; idx_end = sum(L(1:end-j+1)); C_denoised(idx_start:idx_end) = detail_coeffs{level-j+1}; end % 重构降噪后的IMF imf_denoised = waverec(C_denoised, L, wname); end这段代码里有几个细节值得单独说明。
噪声标准差的估计用小波第一层细节系数中位数除以0.6745,这个0.6745的来历是标准正态分布的第75百分位数与中位数之间的关系,是信号处理领域一个经典结论。第一层细节系数包含了信号里最高频的成分,而高频段通常以噪声为主,用中位数法可以稳健地把噪声水平给估算出来,不受个别大幅值瞬态成分的干扰。
阈值计算用的是通用阈值sigma * sqrt(2 * log(N))。这个公式的推导背景是:当N足够大时,白噪声的最大幅值与噪声标准差和log(N)之间存在一个确定的关系,超过这个阈值的系数大概率不是纯噪声。但纯通用阈值有时对某些突变特征过狠,所以我在前面加了一个scale参数,取0.7到0.9之间。实测下来,scale=0.8是一个比较合适的起点,既保留细节又去掉毛刺。
软阈值处理用wthresh函数,'s'表示软阈值。软阈值与硬阈值的区别在这里值得多说一句。硬阈值是低于阈值的系数直接置零,高于阈值的系数原样保留,好处是保幅,坏处是信号会出现不连续的跳变点,产生吉布斯振荡现象。软阈值是让超过阈值的系数向零收缩一个阈值大小,处理出来的信号更平滑连贯,缺点是幅值会有一点点压缩。对于后续要做包络分析、特征量提取的工程信号,软阈值带来的幅值压缩可以通过后续归一化或包络幅值恢复来弥补,但振荡纹波一旦引入就很难抹掉。所以我的默认选择是软阈值。
调用上面的函数,对整个流程的噪声主导IMF批量降噪,再和信号主导IMF以及残差一起重构,就得到了最终的去噪信号。
% ============ 全流程:分解 -> PE筛选 -> 小波降噪 -> 重构 ============ % 接前面,IMFs为ICEEMDAN分解结果,R为残差 denoised_imfs = IMFs; for k = 1:length(noise_imf_idx) idx = noise_imf_idx(k); denoised_imfs(:, idx) = wavelet_denoise(IMFs(:, idx), 'sym8', 6, 0.8); end % 重构:所有IMF(已处理过的用降噪后版本)+ 残差 signal_denoised = sum(denoised_imfs, 2) + R; % 对比原始信号和去噪信号 figure; subplot(2,1,1); plot(x); title('原始含噪信号'); subplot(2,1,2); plot(signal_denoised); title('ICEEMDAN + PE + 小波降噪后的信号');3.5 全流程参数配置速查表
为了让你在实际使用的时候能快速起步,我把整个流程里涉及的关键参数整理成一个速查表。这个表里的值是我在多个项目里反复调校后得到的经验起点,不是绝对标准,但可以作为第一次跑数据的默认配置。
| 环节 | 参数 | 经验取值 | 备注 |
|---|---|---|---|
| ICEEMDAN | Nstd 噪声标准差倍数 | 0.2 | 0.1~0.3之间,信号噪声越小取越小 |
| ICEEMDAN | NE 集合次数 | 100~500 | 数据长取小,数据短取大 |
| ICEEMDAN | MaxIter 最大迭代次数 | 200~500 | 收敛困难时适当增大 |
| PE计算 | m 嵌入维数 | 5~6 | 短信号用4,确保N >> m! |
| PE计算 | tau 时间延迟 | 1~2 | 高采样率下取2或3 |
| PE评估 | 分类阈值 | 0.55~0.75 | 通过标定确定,默认从0.6起步 |
| 小波降噪 | wname 小波基 | sym8 | db4也可,对称性要求高用sym8 |
| 小波降噪 | level 分解层数 | 5~8 | 约等于log2(N)减2~3 |
| 小波降噪 | scale 阈值缩放 | 0.8 | 噪声大取0.9,特征金贵取0.7 |
| 小波降噪 | 阈值类型 | 软阈值 | 追求平滑信号用软阈值 |
这些参数配合使用比单独调某一个更合理。比如你发现重构信号太平滑了,不要先急着调scale,而是检查PE阈值是不是设低了,把PE阈值往上抬一点,让更多分量被当成信号保留下来,往往更有效。参数之间存在联动关系,这是调参过程中最需要注意的一点。
4. 典型应用场景与降噪效果评估
4.1 仿真信号案例:含噪间歇信号
用仿真信号来说明整套方法的效果,是最清晰的验证方式。我构造一个模拟工程信号:一个10Hz正弦波作为低频主成分,叠加一个短时脉冲串模拟设备冲击,再叠加白噪声和随机脉冲干扰。这个信号模拟了机械设备在正常运行中叠加故障冲击的典型场景。
fs = 1000; % 采样率 1000 Hz t = 0:1/fs:5-1/fs; % 5秒信号 N = length(t); % 模拟低频率主信号 + 周期性冲击 signal_main = 2 * sin(2 * pi * 10 * t); impact_time = 1:0.5:5; % 每个0.5秒出现一次冲击 impact = zeros(size(t)); for k = 1:length(impact_time) idx = round(impact_time(k) * fs); if idx < N - 100 impact(idx:idx+100) = 0.8 * sin(2 * pi * 200 * (0:100) / fs) .* exp(-30 * (0:100) / fs); end end % 加入噪声 rng(42); noise = 0.4 * randn(size(t)); % 高斯白噪声 impulse_noise = zeros(size(t)); idx_imp = randi([1 N], 1, 10); % 10个随机脉冲 impulse_noise(idx_imp) = 1.5; x = signal_main + impact + noise + impulse_noise;用这套方法处理后,原始信号、去噪信号以及残差里保留的噪声三者对比,能很清楚地看到:主信号的正弦波形成功保留,冲击脉冲的尖峰没有被抹平,白噪声和随机脉冲噪声则大部分被剥掉。这种“该留的留、该去的去”的效果,正是PE筛选加定向小波降噪组合起来的优势。
4.2 滚动轴承振动信号的实测处理
仿真信号验证是第一步,真正考验这套方法的是实测数据。拿滚动轴承的振动加速度信号来说,这种信号的特点是非平稳性极强,故障产生时会有周期性瞬态冲击,但冲击能量往往集中在高频段,淹没在背景噪声里,传统滤波很难兼顾“去噪”和“保冲击”这两个目标。
实测处理中,我一般是先把原始振动信号做ICEEMDAN分解,然后观察各IMF的PE值分布。实测振动信号的PE分布往往有明显分层现象:前两三个高频IMF的PE值通常在0.8以上,对应纯噪声和干扰;中间几个IMF的PE值在0.3到0.6之间,混有轴承的固有振动成分和部分故障冲击;后面几个低频IMF的PE值可能在0.2以下,主要是转频和倍频成分。
处理策略可以比仿真案例更精细一点:对PE值极高(>0.8)的IMF直接进行小波降噪,对PE值中等(0.4~0.6)的IMF做轻度处理或不处理,对PE值低的IMF完全保留。然后在重构信号上做包络谱分析,故障特征频率的幅值会明显凸显出来。这个效果比单纯用带通滤波要好,因为带通滤波需要预先知道故障频带,而ICEEMDAN方法不需要先验知识。
4.3 效果好不好怎么评价:SNR、RMSE与相关系数
处理完数据,总得有个量化指标证明“效果好”。工程评价上常用的三个指标是信噪比SNR、均方根误差RMSE和相关系数R。
SNR衡量的去噪后信号与理想无噪信号之间的差异。在仿真信号场景下,因为我们知道原始无噪信号是什么,可以直接计算:SNR = 10 * log10( sum(signal_clean.^2) / sum((signal_denoised - signal_clean).^2) )。SNR越高说明去噪后越接近理想信号,一般来说,经过这套流程处理,SNR相比原始含噪信号能提升10dB以上。
RMSE是评估重构信号与理想信号平均偏差的指标,RMSE越小越好。相关系数R则反映去噪信号与理想信号之间的波形相似程度,越接近1越好。
在实测数据场景下没有“理想信号”可以参考,这时可以用一个替代思路:分别计算原始信号和去噪信号的包络谱,观察故障特征频率处幅值的增强倍数。实测中我见过处理前故障特征频率的幅值被噪声淹没,处理后包络谱上对应频率处出现清晰尖峰的情况,这种对比本身就是效果的说服力。
4.4 计算负荷与参数鲁棒性讨论
这套方法不是没有代价,ICEEMDAN的计算复杂度明显高于普通滤波。集合次数NE取200,一个5万点信号跑下来通常需要数十秒到数分钟,具体耗时取决于电脑性能和信号的复杂度。好在工程场景里很多是离线分析,数据采集完成后统一处理,时间成本完全可以接受;如果要做在线监测,就需要在分段处理、降采样、减少集合次数等方面做权衡。
参数鲁棒性方面,PE阈值和小波阈值缩放系数是影响最终重构质量的两个敏感参数。但值得庆幸的是,这两个参数的敏感区间并不狭窄。以PE阈值为例,在0.5到0.7之间的变化,通常不会引起重构信号“质变”,只是微调边界分量的处理力度,所以即使你没有精确找到最优阈值,最终结果也不会离谱到哪里去。小波阈值缩放系数同理,0.7和0.9的区别只是降噪强度略有不同,不至于让结果从好变坏。这种“宽容性”对工程落地非常重要,因为现场不需要一个只能由专家精细调参才能用的方法。
5. 实操避坑经验与常见问题排查
5.1 端点效应和虚假分量的处理
ICEEMDAN虽然改进了模态混叠问题,但端点效应依然存在——信号两端的数据在包络拟合时容易产生过冲或振荡,导致两端出现虚假的IMF成分。这个问题在我第一次处理实测数据时差点被带偏。当时处理一段轴承振动信号,分解出的第一个IMF在两端有明显的异常振荡,PE值算出来也不高,差点被误当成有效信号。
处理端点效应有两个常用手段。一是端点延拓,常用的有镜像延拓、多项式拟合延拓等,在分解之前先对信号两端做适当延拓,分解完再裁掉延拓部分。ICEEMDAN的第三方实现里通常已经内置了某种延拓机制,但效果能否满足要求,还是要自己检查。二是数据切段,在信号两端各丢弃一部分数据,只取中间稳定段进行分析。对振动信号来说,采样时长往往足够长,切掉两端一小段数据不影响整体信息。
还有一个更隐蔽的问题:分解得到的IMF数量不固定,某些额外多出来的分量可能是数值上的人为产物。解决办法是分析前先做归一化和去趋势,也可以对分解后的IMF做一次瞬时频率检查——如果某个IMF的瞬时频率不符合物理常识,就把它视为虚假分量,从重构中剔除。
5.2 排列熵阈值误判该如何纠偏
PE阈值设不好,最常见的两个后果是:噪声残留和信号過平滑。
噪声残留的表现是重构信号里依然有明显的高频毛刺,尤其当原始信号里的噪声并非高斯白噪声,而是带有色噪声成分时,某些噪声主导IMF的PE值可能不高,被误判为信号主导分量保留了下来。遇到这种情况,优先看频谱。有色噪声在特定频段上会有能量集中,把这个IMF的频谱画出来,如果它只在一个很窄的频带上有异常峰,且这个频带和有效信号频带不重叠,那就应该手动把它加入待降噪集合。
信号過平滑的表现则是重构信号包络变圆,冲击尖峰和瞬态突变被压低。这通常是PE阈值设太低,误伤了一些含有有效瞬态成分的分量。纠偏思路是:把被保留的接收信号主导IMF里选几个频段最高的一起看一下,如果其中某几个IMF波形里明显存在小幅但规律的周期性成分,就说明阈值该上调了。另外一个辅助技巧是,对每个IMF做一次包络谱分析,观察有没有清晰的特征频率峰,如果连个像样的峰都没有,那它的PE值就算不高,大概率也是纯噪声。
5.3 小波阈值降噪的两个易忽略的坑
第一个坑是每层共用同一个阈值导致低频细节被过度处理。通用阈值基于第一层细节系数的噪声标准差估算,但噪声在小波各层的能量分布其实不同,直接用同一个阈值处理所有层并不完全合理。更精细的做法是逐层估计噪声标准差,对每一层各算各的阈值。这个做法在信号较长、分解层数较多时效果提升明显。代价是代码量稍增。我的建议是:层数少于6时,可以共用阈值;层数超过6时,建议逐层阈值。
第二个坑是wthcoef函数使用时容易搞混系数索引。很多人第一次接触wavedec返回的C向量都觉得头大,它是一个拼接向量,前一部分是最后一层近似系数,后一部分依次是最深到最浅各层细节系数,L向量用来记录每个段的长度。我在3.4节的代码里用逐层替换的方式重建了C向量,这种方式虽然不够简洁,但不容易出错。如果你习惯用wthcoef,注意它的参数是层级编号的方向,我一朋友在这里踩过一次坑,把层数写反了,最后重构出来的信号整个乱套。所以稳妥起见,我推荐用逐层替换、最后waverec的方式,思路更直白可读。
5.4 计算资源不够时的优化思路
如果你手里的数据很长,比如几十万点甚至成为上百万点的连续监测数据,ICEEMDAN跑起来可能会慢到让人怀疑人生。这种情况下我有三个建议。
第一,降采样。如果有效信号频带最高只有几千赫兹,采样率却高达几十千赫兹,可以先把采样率降下来,数据处理量直接减半不止。注意降采样前必须加抗混叠滤波器,不然频谱会折出假成分。
第二,分段处理。把长信号切成长度适中、有一定重叠的段,分段做ICEEMDAN分解和降噪,最后再拼接回去。拼接处可能出现跳变,重叠区域用线性渐入渐出处理就能消除。这个方法我经常用,计算结果的稳定性和整体处理效果都能保持住。
第三,压缩集合次数NE。NE从500降到100,计算量能减少80%左右,分解精度的损失通常可以接受。真的大规模并行处理时,可以把NE拆成多个批次并行计算,最后汇总取平均,速度会有质的提升。
5.5 常见问题排查速查表
| 现象 | 可能原因 | 排查思路 |
|---|---|---|
| 分解出的IMF在两端振荡剧烈 | 端点效应 | 增加端点延拓,或弃置两端数据后重跑 |
| PE值全部接近1,无法区分 | 嵌入维数m过大或信号太短 | 减小m到4或5,检查信号长度 |
| PE值全部接近0,无法区分 | 信号本身太规则或τ过大 | 增大τ到2或3,检查是否还有噪声在分解前被遗漏 |
| 重构信号仍有很多毛刺 | PE阈值过高,噪声分量被保留 | 降低PE阈值,或手动把高截止频率IMF加入降噪集合 |
| 重构信号过于平滑 | PE阈值过低,信号分量被误伤 | 提高PE阈值,逐层查验IMF包络谱 |
| 小波降噪后幅值明显变小 | 软阈值固有压缩 | 调小scale到0.7,或在后续分析中做幅值恢复 |
| 计算耗时过长 | NE过大或数据太长 | 降低NE、降采样、分段并行处理 |
| 分解出现虚假高频分量 | 数据未归一化或噪声强度过高 | 预处理加去趋势,调大Nstd |
| 重构信号波形与原始相位偏移 | 小波重构系数拼接错误 | 检查wthcoef的层数方向,改用逐层拼接方式 |
这套方法不是银弹,但在我处理的多个实测数据上都表现出比常规方法更高的适应性和稳定性。最后分享两个从实战里打磨出来的小技巧:第一,拿到新数据不要一上来就全流程跑,先截取一小段信号试跑通整个流程,把PE阈值和scale参数定了,再处理全段数据,这样能省下大量试错时间。第二,每次调参都把每个IMF的PE值、判定结果和重构信号一起打印出来对比,形成一套自己的“核对习惯”。很多看似诡异的问题,看一眼IMF波形和PE值分布就能定位。方法的流程本身不复杂,复杂的是让每个参数在你的具体数据上都找到合适的位置,这个过程值得耐心慢慢磨。