基于IVMD的单通道盲源分离方法及其应用*

汤 杰 , 陈 剑,2，杨 斌

(1.合肥工业大学 噪声振动研究所，合肥 230009；2.安徽省汽车NVH技术研究中心，合肥 230009)

0 引言

盲源分离[1](Blind Source Separation)自提出以来在机械故障诊断领域得到了广泛的应用[2]。但是，传统盲源分离方法如ICA算法，JADE算法要求传感器数目必须大于或等于拟分离分量的数目；实际中，在没有振源先验知识的情况下，很难事先确定安装传感器个数；同时，受到安装空间和其他客观因素限制，也无法安装足够多数量的传感器，使得传统盲源分离方法的应用受到限制。因此，研究欠定盲源分离甚至单通道混合信号的盲源分离问题，对于机械故障诊断具有重要意义。

近年来，为解决单通道观测信号盲源分离的问题，不少学者提出了很多有效的方法，如文献[3-4]将经验模态分解(EMD)局域均值分解(LMD)与传统盲源分离方法相结合，实现机械故障信号的单通道盲源分离，上述欠定盲源分离方法往往忽略了实际测量信号非平稳性的问题，加之EMD,LMD在分解过程中模态混叠、过包络、欠包络、端点效应等缺陷，使得分离结果难如人意。为了克服实际测量信号非平稳性对盲源分离结果的影响，文献[5]提出基于时频分析的盲源分离方法。文献[6]将小波分解和基于时频分析的盲源分离方法相结合，用于滚动轴承振动信号单通道盲源分离中，但是小波基的选择较为困难，使得该方法在应用中存在着局限性。为了避免EMD、LMD分解过程中存在的缺点和小波分解过程中小波基选择困难的问题，作为一种自适应信号分解新方法—变分模态分解[8](VMD)被引入到欠定盲源分离中,文献[9]将VMD与FastICA算法相结合，实现单通道滚动轴承复合故障的分离与诊断；但是，对于该方法，VMD分解的层数对分解结果有较大的影响，同时该方法也忽略了实际信号的非平稳性。

基于上述研究，本文提出基于改进变分模态分解(Improved Variational Mode Decomposition, IVMD)和时频分析(Time-Frequency Analysis,TFA)的单通道盲源分离方法。该方法针对VMD分解过程中的分解层数选取问题，提出一种根据谱相关系数确定分解层数的IVMD方法，然后利用IVMD将单通道观测信号分解，实现信号的升维，并根据信号源的数目重构观测信号；对重构后的观测信号进行基于时频分析的盲源分离(Blind Source Separation Based on Time-Frequency Analysis,TFA-BSS)，得到源信号的估计信号。

1 基本理论

1.1 变分模态分解

变分模态分解[8]是以经典维纳滤波、希尔伯特变换和频率混合这三个概念为基础的变分问题求解方法，通过搜寻约束变分模型最优解来实现信号自适应分解。

(1)

(2)

(3)

利用Parseval/Plancherel傅里叶等距变换到频域后用ω-ωk代替ω并将其转化为非负频率区间积分的形式，则优化问题的最终解为：

(4)

(5)

1.2 基于时频分析的盲源分离

TFA-BSS[5]主要包括如下两部分：

(1)对观测信号进行白化处理。

假设x(t)=(x1(t),x2(t),…,xm(t))T为m维观测信号，s(t)=(s1(t),s2(t),…,sn(t))T为n维源信号，则两者之间的关系可以表示成：x(t)=As(t)+n(t)，式中A为混合矩阵，n(t)为加性白噪声，且m≥n。观测信号x(t)的自相关矩阵定义如下：

Rxx=E[x(t)x(t)*]

(6)

其中，上标*表示复共轭，对Rxx进行特征值分解得到其特征值和相应的特征向量，Rxx前n个最大的特征值和相应的特征向量分别用λ1,λ2,…,λn和h1,h2,…,hn表示。根据Rxx剩余m-n个最小特征值估计噪声的方差σ2和白化矩阵W。

(7)

白化矩阵为：

(8)

其中，上标H表示复共轭转置。经过噪声补偿和白化后的观测信号可用下式表示：

z(t)=Wx(t)=WAs(t)=Us(t)

(9)

(2)矩阵联合对角化，对上式左右两端进行时频变换(时频变换采用平滑伪Wigner-Ville)，得到源信号与白化观测信号时频分布矩阵的关系：

DZZ(t,f)=WADss(t,f)AHWH=UDss(t,f)UH

(10)

将源信号空间时频分布矩阵Dss(t,f)的对角元素标记为自项，表示各源信号自身的时频分布；而非对角元素标记为互项，表示两个源信号之间的时频分布。按一定规则选择L个自项点，并对自项位置处的时频分布矩阵进行联合对角化，估计酉阵U；根据上述估计的白化矩阵W和酉矩阵U，可以得到源信号的估计信号为:s(t)=UHWx(t)。

2 基于IVMD的单通道盲源分离

2.1 改进变分模态分解

由VMD算法原理可知，其处理信号时需要预先设定BLIMFs分量的个数K，然而受不同设备不同工况的限制，K值通常难以准确设定。关于VMD分解K值设定的问题，文献[10]以BLIMFs分量与原信号的互相关系数为适应度函数，采用粒子群算法优化VMD的参数，以选择合适的分解层数。文献[11]为减少人为主观选择VMD参数存在的弊端，提出了基于包络谱特征因子的参数自动搜索策略。文献[12]以BLIMFs分量的局部峭度最大化为目标，采用人工鱼群算法优化变分模态分解的参数。但是上述方法计算量较大，在实际中不易应用。鉴于VMD是在频域中实现非递归、变分模式的分解，笔者提出以BLIMFs分量与原信号间谱相关系数的最小值作为原算法迭代停止条件。谱相关系数[13]可定量评价两信号频谱的相似程度,谱相关系数的定义如下式：

(11)

式中,|X(k)|、|Y(K)|分别表示两信号的傅里叶变换的模，n是频域离散值序列号。根据VMD的带通滤波特性，当原始信号出现过分解时，在分解出的K个BLIMFs中，会存在与原信号谱相关系数较小的模态分量，IVMD具体实现过程如下：

(1) 初始化K=1，α=2000，τ=0.95；

(2)K=K+1,执行外层循环；

(4)n=n+1，开始内层循环；

(7) 由上述得到K个BLIMFs分量，计算各BLIMFs分量与原始信号的谱相关系数：ρ1,ρ2,ρi,…,ρK;ρi表示第i个模态分量与原信号间的谱相关系数；

(8) 在上述求得的谱相关系数中找出最小值ρmin=min{ρ1,ρ2,ρi,…,ρK}，判断ρmin是否小于给定阈值;若小于，结束外层循环，确定分解层数为K-1层；否则，重复步骤(2)～(8)。

2.2 算法流程

根据章节2.1所述的IVMD原理对单通道观测信号进行分解，得到K个BLIMFs。然后将单通道观测信号与各模态分量及其残余项一起组成虚拟多维观测信号，实现信号的升维，以满足盲源分离对观测矩阵正定的要求。将上述虚拟多维观测信号的自相关矩阵进行奇异值分解，采用贝叶斯信息准则得到信号源数目的估计(参考文献[14]),根据相关系数最大原则选择若干BLIMFs分量与单通道观测信号及其分解残余项共同组成新的观测矩阵；以上欠定盲分离的流程如图1。

图1 信号盲分离流程

3 仿真分析

本节利用仿真信号探讨结合VMD与TFA-BSS方法的有效性。轴承发生故障时，因为受力不均匀会引发振动信号的幅值调制甚至频率调制。据此，仿真源信号表达式如公式(12)，设置采样频率为6000Hz，采样时间为0.5s，为了更清晰的表示信号波形，图中仅表示出0.1s的波形。源信号的波形和频谱如图2和图3所示，其中f1=100Hz,f2=150Hz,f3=180Hz,f4=500Hz,f5=1000Hz。

(12)

(a)源信号1 波形

(b)源信号2波形

(c)源信号3波形

(a)源信号1 波形

(b)源信号2波形

(c)源信号3波形

实际应用中，传感器采集到的振动数据大多是由多个源信号混合而成。为了模拟这一现象，随机选取一个3×3的混合矩阵A，

根据章节1.3所述，采用线性混合的方法得到三个混合信号。同时，为了模拟真实数据中的噪声，向三个混合信号添加随机噪声，使得最终混合信号x(t)的信噪比达到20dB。混合信号的时域波形如图4所示。

(a)观测信号1 波形

(b)观测信号2波形

(c)观测信号3波形

假设实际应用中，仅仅监测到了一路观测信号，选取图4的观测信号1作为仿真实验的单通道观测信号x1(t)，对其进行分析、验证。根据上文所述IVMD方法将单通道观测信号分解，随着分解层数的增加，各模态分量与原信号间的谱相关系数如表1。

表1 各模态分量与原始信号的谱相关系数

从表1可以看出当分解层数为5时，第5个分量和原信号的谱相关系数小于设定阈值0.1，可以认为出现了过分解。因此，K值设为4，将观测信号x1(t)分解为4个有限带宽固有模态函数(限于文章篇幅，此处未将VMD处理的分量在图中表示出)，将原始单通道信号与各模态分量及其残余项一起组成虚拟6维观测信X=(x,c1,c2,...,c4,r)T。通过计算虚拟多维观测信号的特征值，根据贝叶斯信息准则判断信号源个数为3，然后根据表1的谱相关系数，选取BLIMF1和BLIMF2与单通道观测信号x1(t)及其分解残余量r组成四维观测信号，经过TFA-BSS算法分离后得到4个估计信号(图中仅表示出与源信号1、2、3相对应的估计信号；另一估计信号为噪声，未在图中表示。)如图5所示。

(a)估计信号1 波形

(b)估计信号2波形

(c)估计信号2波形

采用本文所述基于IVMD的单通道盲源分离方法，结果如图5；从图5明显可以看出s3对应于估计信号1，s2对应于估计信号2，s1对应于估计信号3，并且估计信号2和估计信号3产生了一定的模态混叠。为了定量比较分离性能，采用相关系数对分离结果进行评价；同时为了对比研究，采用文献[4]和文献[5]所述方法，对上述仿真信号进行分离，三种方法分离结果评价指标见表2。

表2 不同算法分离结果相关系数对比

对比表2中数据，本方法分离信号s1、s3的相关系数明显优于另外两种；s2虽说出现了一定的失真但相对于另外两种方法，相关系数有所提高。

4 实验分析

4.1 数据来源

为了验证本方法在实际故障诊断中的有效性，以美国凯斯西储大学轴承数据中心 (Case Western Reserve University (CWRU) Bearing Data Center Website)发布的故障数据为研究对象，开展验证性分析。轴承型号为: SKF6025-2RS。传感器分别安装在电机驱动端、电机风扇端、试验台基座上。本文分析数据来源于驱动端加速度传感器。使用电火花加工技术分别在轴承内圈、外圈上布置单点故障，电机转速为1730r/min, 故障直径为0.5334mm，采样频率为12kHz,根据电机转速计算外圈故障特征频率为103.36Hz，内圈故障特征频率为156.13Hz。滚动轴承的有关参数如表3所示。

表3 滚动轴承参数

4.2 数据分析

试验采集的轴承外圈故障信号波形与内圈故障信号波形如图 6a 和图 6b所示。

(a)内圈故障信号

(b)外圈故障信号

将上述内圈故障信号和外圈故障信号通过一个随机1×2矩阵A混合，得到单通道观测信号，信号时域波形如图7所示。

A=[0.73 0.24]

图7 单通道复合故障观测信号波形

为了对比分析，笔者首先采用文献[15]所述的频谱自相关分析方法对信号进行特征提取，得到的结果如图8所示。可以看出，内圈故障特征频率及其倍频处有很明显的峰值，外圈故障特征被转频和内圈故障所淹没，无法发现。

图8 单通道复合故障观测信号频谱自相关分析

采用本文所述方法，单通道复合故障观测信号随着分解层数的增加最小谱相关系数的变化见图9所示；当分解层数K=10时出现小于0.1的谱相关系数，表明已过分解，故分解层数K取9。对重组虚拟观测矩阵，采用基于贝叶斯信息准则的源数估计方法估计出源数为5，然后利用TFA-BSS方法对重构观测信号进行盲分离，得到分离信号时域波形见图10(此处仅表示出与故障有关的两个分离信号)，分离信号的频谱自相关分析如图11所示。

图9 最小谱相关系数的变化

(a)分离信号1时域图

(b) 分离信号2时域图

(a)分离信号1频谱自相关分析

(b)分离信号2频谱自相关分析

从图11a可以明显看到主峰值出现在104Hz处，这与轴承外圈故障频率理论值103.36Hz非常接近(实际值和理论值之间的差异可能与电机转速不稳定、频率分辨率设置、轴承打滑等因素有关),同时外圈故障频率的二倍频也被凸显出来。相对于图8，外圈故障辨识度有了很大的提高。图11b中156.7Hz处出现的峰值，与内圈故障频率156.13Hz较接近，并且在其二倍频处也出现明显峰值，这与图8的结果一致。综上分析可知，原故障信号中不仅包含内圈故障成分，同时也包含外圈故障成分，外圈故障成分在原信号的频谱自相关分析中没有表现出来是由于被转频和内圈故障所淹没。通过盲源分离的方法将两故障信号分配在不同的通道中，便于故障的准确辨识。

5 结论

(1) 以原信号和各BLIMFs分量间的最小谱相关系数为准则对VMD的迭代停止条件进行了改进，实现了分量个数的自适应确定。

(2) 针对单观测通道盲源分离的问题，将改进后的变分模态分解与基于时频分析的盲源分离相结合，实现单通道观测信号的盲源分离。通过与CWD-TFA-BSS，EMD-TFA-BSS方法的对比，验证了本方法的优势。

(3) 通过对实际滚动轴承混合故障信号的分析，证明本方法对滚动轴承复合故障诊断更加准确。