【时频分析】变分广义非线性模态分解算法及其应用【附MATLAB代码】

文章来源微信公众号：EW Frontier

QQ交流群：192253249

摘要

近年来提出的变分信号分解方法，如自适应线性调频模式分解（ACMD）和广义色散模式分解（GDMD）在各个领域引起了广泛的关注。然而，这些方法难以同时分离啁啾分量和色散分量。本文提出了一个变分广义非线性模式分解（VGNMD）框架来解决这个问题。VGNMD首先引入自适应时频融合与聚类（ATFFC）方法，提高噪声环境下信号时频分布（TFD）的抗噪性和分辨率，准确获取各模式的TFD。然后，模式类型的判别标准，建立了模式分为啁啾模式或色散模式的基础上，其时间-频率（TF）的脊。最后，以这些TF脊作为初始瞬时频率（IF）或初始群延迟（GD），应用变分优化算法精确重构模态并细化其IF或GD。仿真例子和实际应用蝙蝠回声定位信号分析和铁路轮轨故障诊断被认为是VGNMD的有效性。结果表明，该方法可以同时准确地提取啁啾模和色散模，适用于分析具有不连续TF模式的非线性信号。

引言

机械设备的状态监测与故障诊断对于保证其安全稳定运行具有重要意义[1-5]。信号分解通过将信号分解为一系列具有明确物理意义的子信号分量，可以抑制机械振动信号中的环境噪声和干扰分量[6，7]。随着机械设备故障诊断向精细化方向发展，信号分解在微弱故障诊断和定量识别中发挥着越来越重要的作用[8]。然而，对于轨道车辆等大型机械系统，复杂的运行环境（多源激励、时变工况等）使得状态监测信号往往呈现出多模混频、复杂非平稳调制、强噪声干扰的特征[9，10]。具体地说，由于机械系统中零件的耦合效应，一个零件的振动信号可能包含来自多个零件的多个分量，甚至可能具有两个类谐波（即，瞬时频率（IF）随时间变化的线性调频模式[11]）和脉冲状分量（即，具有随频率变化的群延迟（GD）的色散模式[12]）。准确分离这些模式并提取其IF和GD是一项非常具有挑战性和实用性的工作[13，14]。

大多数现有的信号分解方法主要致力于某类分量的分离（即，啁啾模式或色散模式）。啁啾模式通常携带关于机器故障的类型和程度的重要信息[15]。许多分解方法，如经验模式分解（EMD）[16]，局部均值分解[17]，经验小波变换[18]和变分模式分解（VMD）[19]，用于分离啁啾模式。这些方法由于其独特的优点而引起了各个领域的广泛关注，例如高适应性（例如，EMD）[20]和强噪声鲁棒性（例如，VMD）[21]。然而，上述方法中的大多数限于分解在频域中分离良好的窄带模式[22]。为了有效地分析宽带非平稳信号，最近提出了基于最佳解调算法的变分非线性啁啾模式分解（VNCMD）[23]。与VMD [24]的频域滤波器组不同，VNCMD本质上是一个时频（TF）域滤波器组，其中心频率是估计的IF，可以分解具有非常接近甚至交叉模式的啁啾信号。此外，提出了自适应线性调频模式分解（ACMD）[25]和数据驱动ACMD（DD-ACMD）[26]，以增强VNCMD的适应性。注意，虽然上述变分信号分解方法在分析非平稳信号时逐渐满足高适应性和良好的噪声鲁棒性，但它们不适用于色散信号分析，因为chirp模式模型无法表征色散效应。

频散模常用于声源定位和系统识别。对于色散信号分解，色散补偿方法（DCM）可以有效地降低色散效应，并引起了人们的极大兴趣[27]。该方法将具有频率变化GD的色散波形压缩成具有恒定GD的时间脉冲，该时间脉冲可以通过矩形时间窗口提取。受DCM和ACMD的启发，开发了广义色散模式分解（GDMD）方法来处理重叠或交叉色散模式[28]。该方法将模式分解问题转化为最优色散补偿问题，通过联合估计算法求解。本质上，最佳色散补偿是VNCMD和ACMD中的最佳解调算法的频域对偶版本。由于其良好的分解能力和噪声鲁棒性，GDMD被广泛用于结构健康监测和多通道信号分解[29]。值得注意的是，这些变分信号分解方法通过变分优化算法提取信号模式，该算法需要由初始IF或初始GD驱动。具体而言，ACMD通过希尔伯特变换获得初始IF，以驱动变分优化算法在时域中提取啁啾模式。GDMD通过提取TF脊来获得初始GD，以驱动变分优化算法在频域中提取色散模式。然而，受色散模模型的限制，这些方法很难同时提取色散模和啁啾模。此外，现实生活中的信号可能具有某些不连续的分量，这些分量仅存在于参考间隔中的有限持续时间和频带内。上述方法中的大多数在整个参考区间上定义模态模型，因此它们不能提取这些不连续分量。

与上述信号分解方法不同的是，时频分析（TFA）可以在TF域同时表征啁啾特征和色散特征。然而，由于受Heisenberg-Gabor不等式的限制，传统的时频分析方法得到的时频分布往往能量集中性较差。为了提高TFD的分辨率，研究人员开发了许多基于TF重新分配的后处理方法，如同步压缩变换（SST）[30]及其改进版本[31-33]。遗憾的是，上述方法仅将TF系数沿着一个方向移动以保持信号重构能力，因此不能同时提高时间分辨率和频率分辨率。为了解决这些问题，基于“分而治之”策略提出了时频多压缩变换（TFMST）[34]，该变换不仅可以有效地表征类谐波和类脉冲分量的混合信号，而且具有信号重构的能力。此外，TF融合也是一种通过融合不同参数的TFD来提高TF分辨率的有效技术[35-37]。值得注意的是，上述方法大多集中于获得高分辨率的TFD，并且没有给予同时分离啁啾分量和色散分量的有效策略。

如前所述，同时提取非线性信号中的啁啾模式和色散模式是一项具有挑战性的工作。受TF融合技术和变分信号分解方法ACMD和GDMD的启发，提出了一种变分广义非线性模式分解（VGNMD）框架，用于自适应地对模式进行分类，并进一步精确地提取这些模式及其IF或GD。与ACMD和GDMD不同的是，VGNMD能够自适应地获得各模式的类型及其时频脊，并以此为初始参数驱动变分优化算法，从复杂的非平稳信号中同时提取啁啾模式和色散模式。本文的主要贡献可概括如下。首先，提出了一种自适应TF融合与聚类（ATFFC）方案，提高TF分辨率，消除噪声干扰，进一步获得高质量的各模式的TFD。然后，引入了一种模式类型判别准则，根据模式的TF脊将模式分为啁啾模式和色散模式。同时，这些模式的支持区间是由TF脊的非零值得到的。最后，在这些TF脊的驱动下，变分优化算法精确地提取模式及其IF或GD。此外，用这些模式和IF或GD重建高分辨率TFD。

本文的其余部分组织如下。在第二节中，我们介绍了广义非线性信号模型，并回顾了ACMD和GDMD中的变分优化算法。拟议的VGNMD框架详见第三节。在第四节中，我们提供了一些仿真和实际信号来测试VGNMD的性能。我们在第五节结束我们的论文。

VGNMD算法流程

GNMD仿真实验

本文采用EMD、VMD、ACMD和VNGMD对不同噪声污染的信号进行了分析。图11示出了通过不同方法的噪声信号的重建结果。可以看出，由于强噪声的干扰，EMD不能正确地提取信号。尽管VMD被设置为一个相对准确的初始中心频率（200 Hz），但由于窄带限制，VMD不能获得可接受的重建结果。ACMD的重建结果在2.6 ~ 4s内没有收敛到正确的结果，因为在这个时间段内的初始IF远离真实IF。相比之下，VGNMD可以很好地重建被强噪声污染的信号。此外，EMD、VMD、ACMD和VGNMD的代价时间分别为0.05 s、0.09 s、0.23 s和7.61 s。诚然，VGNMD需要更高的计算成本。

了分析噪声强度对VGNMD分解性能的影响，我们计算了不同噪声水平下四种方法重构模式的信噪比，结果如图12所示。可以看出，四种方法的重建结果均随噪声强度的增加而降低。EMD和VMD由于不适用于宽带信号，即使在噪声强度较低的情况下也无法获得高信噪比的重构结果。相反，ACMD和VGNMD都可以在弱噪声环境下很好地重建信号，并且所提出的VGNMD在强噪声环境下也表现出优异的重建性能。

建议的VGNMD被用来分解信号。图14示出了针对噪声信号的ATFFC方案的TFD和聚类结果。可以看出，TFD的浓度显著提高，并且噪声被完全去除。此外，每一个模式的TFD精确分割和模式的持续时间和有效带宽的准确获得。

GNMD MATLAB代码

%%%%%%%%% Generalized nonlinear signal %%%%%%%%%
clc;
clear;
close all
SampFreq = 1000;
​
t1 = 0:1/SampFreq:1.5;t1 = t1(1:end-1);
t2 = 1.5:1/SampFreq:3;t2 = t2(1:end-1);
t = [t1 t2];
Sig1 = cos(2*pi*(160*t+20*t.^2+3*cos(3*pi*t)));%mode1
IF1 = 160+40*t-9*pi*sin(3*pi*t);%IF1
​
Sig21 = cos(2*pi*(70*t1+20*t1.^2));
Sig22 = zeros(1,length(t2));
Sig2 = [Sig21,Sig22];%mode2
IF2 = 70+40*t1;%IF2
​
Sig31 = zeros(1,length(t1));
Sig32 = cos(2*pi*(20*t2+20*t2.^2+3*cos(3*pi*t2)));%mode3
Sig3 = [Sig31,Sig32];
IF3 = 20+40*t2-9*pi*sin(3*pi*t2);%IF3
​
T = 3; 
Nt = 3000;
Nf = floor(Nt/2)+1;
​
f41 = (0 : 1*Nf/2-1)/T;f42 = (1*Nf/2 : Nf)/T;
Ds41 = linspace(0,0,length(f41));Ds42 = 30*exp(-j*2*pi*(0.4*f42+2*cos(2*pi*f42/100)));
GD42 = 0.4-4*pi/100*sin(2*pi*f42/100);%GD4
CDs41 = [complex(Ds41) Ds42];Dss41 = real(CDs41);
iDFs41 = [CDs41,conj(fliplr(CDs41(2:ceil(Nt/2))))]; ifftSig4 = ifft(iDFs41);
Sig4 = real(ifftSig4);%mode4
​
f51 = (0 : 3*Nf/5-1)/T;f52 = (3*Nf/5 : Nf)/T;
Ds51 = linspace(0,0,length(f51));Ds52 = 30*exp(-j*2*pi*(0.8*f52+0.0005*f52.^2));
GD52 = 0.8+0.001*f52;%GD5
CDs51 = [complex(Ds51) Ds52];Dss51 = real(CDs51);
iDFs51 = [CDs51,conj(fliplr(CDs51(2:ceil(Nt/2))))]; ifftSig5 = ifft(iDFs51);
Sig5 = real(ifftSig5);%mode5
​
f61 = (0 : 7*Nf/10-1)/T;f62 = (7*Nf/10 : Nf)/T;
Ds61 = linspace(0,0,length(f61));Ds62 = 30*exp(-j*2*pi*(1.8*f62+2*cos(2*pi*f62/100)));
GD62 = 1.8-4*pi/100*sin(2*pi*f62/100);%GD6
CDs61 = [complex(Ds61) Ds62];Dss61 = real(CDs61);
iDFs61 = [CDs61,conj(fliplr(CDs61(2:ceil(Nt/2))))]; ifftSig6 = ifft(iDFs61);
Sig6 = real(ifftSig6);%mode6
​
f71 = (0 : 8*Nf/10-1)/T;f72 = (8*Nf/10 : Nf)/T;
Ds71 = linspace(0,0,length(f71));Ds72 = 30*exp(-j*2*pi*(2.2*f72+0.0005*f72.^2));
GD72 = 2.2+0.001*f72;%GD7
CDs71 = [complex(Ds71) Ds72];Dss71 = real(CDs71);
iDFs71 = [CDs71,conj(fliplr(CDs71(2:ceil(Nt/2))))]; ifftSig7 = ifft(iDFs71);
Sig7 = real(ifftSig7);%mode7
​
Sig = Sig1+Sig2+Sig3+Sig4+Sig5+Sig6+Sig7; %Generalized nonlinear signal
noise = addnoise(length(Sig),0,0.4);%noise
Sign = Sig+noise;
%% Time domain waveform
figure
set(gcf,'Position',[20 100 640 500]);
set(gcf,'Color','w');
plot(t,Sign,'b-','linewidth',1);
xlabel('Time (s)','FontSize',24,'FontName','Times New Roman');
ylabel('Amplitude','FontSize',24,'FontName','Times New Roman');
set(gca,'FontSize',24,'FontName','Times New Roman');
set(gca,'looseInset',[0.02 0.02 0.02 0.02])
axis([0 3 -6 6]);
%% Short-time Fourier transform
[Spec,f] = STFT(Sign',SampFreq,length(Sign),SampFreq/2);
figure
imagesc(t,f,Spec);  
set(gcf,'Position',[20 100 640 500]);  
axis([0 3 0 SampFreq/2]);
xlabel('Time (s)','FontSize',24,'FontName','Times New Roman');
ylabel('Frequency (Hz)','FontSize',24,'FontName','Times New Roman');
set(gca,'FontSize',24,'FontName','Times New Roman');
set(gca,'YDir','normal');set(gca,'looseInset',[0.02 0.02 0.02 0.02])
colormap('jet')
%% Variational generalized nonlinear mode decomposition (VGNMD)
tic
[Modet,Modef,Type,IIF_GD,EIF_GD,f,t] = VGNMD(Sign,SampFreq);
toc
%% IIF eIF
IIF_GD1 = IIF_GD{1};IIF1 = f(IIF_GD1(:,2)); IIF1 = curvesmooth(IIF1,1e-6);
IIF_GD2 = IIF_GD{2};IIF2 = f(IIF_GD2(:,2)); IIF2 = curvesmooth(IIF2,1e-6);
IIF_GD3 = IIF_GD{3};IGD3 = t(IIF_GD3(:,2)); IGD3 = curvesmooth(IGD3,1e-6);
IIF_GD4 = IIF_GD{4};IGD4 = t(IIF_GD4(:,2)); IGD4 = curvesmooth(IGD4,1e-6);
IIF_GD5 = IIF_GD{5};IIF5 = f(IIF_GD5(:,2)); IIF5 = curvesmooth(IIF5,1e-4);
IIF_GD6 = IIF_GD{6};IGD6 = t(IIF_GD6(:,2)); IGD6 = curvesmooth(IGD6,1e-6);
IIF_GD7 = IIF_GD{7};IGD7 = t(IIF_GD7(:,2)); IGD7 = curvesmooth(IGD7,1e-5);
​
figure
set(gcf,'Position',[20 100 640 500]);
plot(t,IF1,'b-','linewidth',1.5);
hold on
plot(t(IIF_GD2(:,1)),IIF2,'k-.','linewidth',1.5);
hold on
plot(t(IIF_GD2(:,1)),EIF_GD{2}(IIF_GD2(:,1)),'r--','linewidth',1.5);
​
hold on
plot(t1,IF2,'b-','linewidth',1.5);
hold on
plot(t(IIF_GD1(:,1)),IIF1,'k-.','linewidth',1.5);
hold on
plot(t(IIF_GD1(:,1)),EIF_GD{1}(IIF_GD1(:,1)),'r--','linewidth',1.5);
​
hold on
plot(t2,IF3,'b-','linewidth',1.5);
hold on
plot(t(IIF_GD5(:,1)),IIF5,'k-.','linewidth',1.5);
hold on
plot(t(IIF_GD5(:,1)),EIF_GD{5}(IIF_GD5(:,1)),'r--','linewidth',1.5);
​
hold on
plot(GD42,f42,'b-','linewidth',1.5);
hold on
plot(IGD3,f(IIF_GD3(:,1)),'k-.','linewidth',1.5);
hold on
plot(EIF_GD{3}(IIF_GD3(:,1)),f(IIF_GD3(:,1)),'m--','linewidth',1.5);
​
hold on
plot(GD52,f52,'b-','linewidth',1.5);
hold on
plot(IGD4,f(IIF_GD4(:,1)),'k-.','linewidth',1.5);
hold on
plot(EIF_GD{4}(IIF_GD4(:,1)),f(IIF_GD4(:,1)),'r--','linewidth',1.5);
​
hold on
plot(GD62,f62,'b-','linewidth',1.5);
hold on
plot(IGD6,f(IIF_GD6(:,1)),'k-.','linewidth',1.5);
hold on
plot(EIF_GD{6}(IIF_GD6(:,1)),f(IIF_GD6(:,1)),'r--','linewidth',1.5);
​
hold on
plot(GD72,f72,'b-','linewidth',1.5);
hold on
plot(IGD7,f(IIF_GD7(:,1)),'k-.','linewidth',1.5);
hold on
plot(EIF_GD{7}(IIF_GD7(:,1)),f(IIF_GD7(:,1)),'r--','linewidth',1.5);
​
xlabel('Time (s)','FontSize',24,'FontName','Times New Roman');
ylabel('Frequency (Hz)','FontSize',24,'FontName','Times New Roman');
set(gca,'FontSize',24,'FontName','Times New Roman');
axis([0 3 0 500]);
set(gca,'looseInset',[0.02 0.02 0.02 0.02])
legend('Ture','Intial','Estimated');legend('boxoff');
%% Estimated modes and errors.
figure
set(gcf,'Position',[20 100 1280 640]);
​
axes('position',[0.06 0.83 0.42 0.13]);
plot(t,Modet(2,:),'b-','linewidth',1);  % estimated mode
hold on
plot(t,Sig1-Modet(2,:),'k--','linewidth',1);
axis([0 3 -2 2]);
set(gca,'xtick',[0:1:3]);set(gca,'ytick',[-2 0 2]);
ylabel('m1','FontSize',24);
set(gca,'FontSize',24,'FontName','Times New Roman','linewidth',1);
​
axes('position',[0.56 0.83 0.42 0.13]);
plot(t,Modet(1,:),'b-','linewidth',1);  % estimated mode
hold on
plot(t,Sig2-Modet(1,:),'k--','linewidth',1);
axis([0 3 -2 2]);
set(gca,'xtick',[0:1:3]);set(gca,'ytick',[-2 0 2]);
ylabel('m2','FontSize',24);
set(gca,'FontSize',24,'FontName','Times New Roman','linewidth',1);
​
axes('position',[0.06 0.60 0.42 0.13]);
plot(t,Modet(5,:),'b-','linewidth',1);  % estimated mode
hold on
plot(t,Sig3-Modet(5,:),'k--','linewidth',1);
axis([0 3 -2 2]);
set(gca,'xtick',[0:1:3]);set(gca,'ytick',[-2 0 2]);
ylabel('m3','FontSize',24);
set(gca,'FontSize',24,'FontName','Times New Roman','linewidth',1);
​
axes('position',[0.56 0.60 0.42 0.13]);
plot(t,Modet(3,:),'b-','linewidth',1);  % estimated mode
hold on
plot(t,Sig4-Modet(3,:),'k--','linewidth',1);
axis([0 3 -4 4]);
set(gca,'xtick',[0:1:3]);set(gca,'ytick',[-4 0 4]);
ylabel('m4','FontSize',24);
set(gca,'FontSize',24,'FontName','Times New Roman','linewidth',1);
​
axes('position',[0.06 0.37 0.42 0.13]);
plot(t,Modet(4,:),'b-','linewidth',1);  % estimated mode
hold on
plot(t,Sig5-Modet(4,:),'k--','linewidth',1);
axis([0 3 -4 4]);
set(gca,'xtick',[0:1:3]);set(gca,'ytick',[-4 0 4]);
ylabel('m5','FontSize',24);
set(gca,'FontSize',24,'FontName','Times New Roman','linewidth',1);
​
axes('position',[0.56 0.37 0.42 0.13]);
plot(t,Modet(6,:),'b-','linewidth',1);  % estimated mode
hold on
plot(t,Sig6-Modet(6,:),'k--','linewidth',1);
axis([0 3 -4 4]);
set(gca,'xtick',[0:1:3]);set(gca,'ytick',[-4 0 4]);
ylabel('m6','FontSize',24);
set(gca,'FontSize',24,'FontName','Times New Roman','linewidth',1);
​
axes('position',[0.06 0.14 0.42 0.13]);
plot(t,Modet(7,:),'b-','linewidth',1);  % estimated mode
hold on
plot(t,Sig7-Modet(7,:),'k--','linewidth',1);
axis([0 3 -4 4]);
set(gca,'xtick',[0:1:3]);set(gca,'ytick',[-4 0 4]);
ylabel('m7','FontSize',24);
set(gca,'FontSize',24,'FontName','Times New Roman','linewidth',1);
xlabel('Time (s)','FontSize',24);
​
axes('position',[0.56 0.14 0.42 0.13]);
plot(t,Modet(1,:)+Modet(2,:)+Modet(3,:)+Modet(4,:)+Modet(5,:)+Modet(6,:)+Modet(7,:),'b-','linewidth',1);  % estimated mode
hold on
plot(t,Sig-(Modet(1,:)+Modet(2,:)+Modet(3,:)+Modet(4,:)+Modet(5,:)+Modet(6,:)+Modet(7,:)),'k--','linewidth',1);
axis([0 3 -6 6]);
set(gca,'xtick',[0:1:3]);set(gca,'ytick',[-6 0 6]);
ylabel('m','FontSize',24);
set(gca,'FontSize',24,'FontName','Times New Roman','linewidth',1);
xlabel('Time (s)','FontSize',24);
%% Reconstructed TFD
band = [0 SampFreq/2];
timerange = [0 t(end)];
eIF1 = EIF_GD{2};eIF2 = EIF_GD{1};eIF3 = EIF_GD{5};
eGD4 = EIF_GD{3};eGD5 = EIF_GD{4};eGD6 = EIF_GD{6};eGD7 = EIF_GD{7};
[ASpec1 fbin] = TFspec(eIF1,abs(hilbert(Modet(2,:))),band);
[ASpec2 fbin] = TFspec(eIF2,abs(hilbert(Modet(1,:))),band);
[ASpec3 fbin] = TFspec(eIF3,abs(hilbert(Modet(5,:))),band);
​
[ASpec4 tbin] = TFspec_GD(eGD4(1:length(f)),hilbert(Modef(3,:)),timerange);
[ASpec5 tbin] = TFspec_GD(eGD5(1:length(f)),hilbert(Modef(4,:)),timerange);
[ASpec6 tbin] = TFspec_GD(eGD6(1:length(f)),hilbert(Modef(6,:)),timerange);
[ASpec7 tbin] = TFspec_GD(eGD7(1:length(f)),hilbert(Modef(7,:)),timerange);
ASpec = ASpec1+ASpec2+ASpec3+ASpec4'+ASpec5'+ASpec6'+ASpec7';
figure
imagesc(t,fbin,abs(ASpec)); 
set(gcf,'Position',[20 100 640 500]);  
axis([0 3 0 500]);
set(gca,'xtick',[0:1:3]);set(gca,'ytick',[0:100:500]);
xlabel('Time (s)','FontSize',24,'FontName','Times New Roman');
ylabel('Frequency (Hz)','FontSize',24,'FontName','Times New Roman');
set(gca,'FontSize',24,'FontName','Times New Roman');
set(gca,'YDir','normal');set(gca,'looseInset',[0.02 0.02 0.02 0.02])
colormap('jet')
​

VGNMD 仿真结果