0 引言

盲源分离(Blind Source Separation,BSS)是指在源信号与混合通道参数均未知的条件下,仅通过传感器观测信号来估计源信号和未知通道参数的一种新兴信号处理新方法^[1]。其广泛应用于通信、语音、图像分析和机械设备故障诊断等领域。目前,在机械设备故障诊断盲分离中,往往假设观测信号的数目不少于机械设备源信号的数目,而这个假设在机械设备故障诊断实践中并不成立^[2]。

因此,在机械设备故障诊断中,恰定信号盲源分离的方法并不能满足实际需求,需要研究一种能够适用于观察信号少于源信号的欠定信号盲分离方法。

为了解决欠定信号盲分离问题,大多数研究人员通过一定的信号处理方法将欠定信号升维为恰定信号,再通过恰定信号盲分离方法估计源信号。孟宗等^[3]提出一种融合小波分解与时频分析的单通道盲源分离方法,并在轴承故障诊断实验中验证了算法的有效性,但是母小波选择的困难使得该方法具有一定的局限性。为了使欠定信号能够自适应地分解为不同的分量,Mijovic等^[4]提出了基于经验模态分解(Empirical Mode Decomposition,EMD)和独立成分分析(Independent Component Analysis,ICA)的单通道信号欠定盲源分离方法。但是EMD算法存在模态混叠和端点效应^[5],并且算法的物理意义不明确^[6]。针对EMD存在的模态混叠现象,Wu等^[7]提出总体经验模态分解(Ensemble EMD,EEMD)方法。孟宗等^[8]提出利用EEMD将信号从一维转换为多维,将欠定的盲源分离问题转化为恰定的盲源分离问题。但是,EEMD通过加噪处理的有效性难以保障,白噪声的污染效应会进入所有模态而不仅仅局限于高频部分。因此,加噪处理的副作用很大,上百次的加噪处理很可能会使信号面目全非^[9]。

王金良等^[9-10]提出一种新的数据分析方法,即极点对称模态分解(Extreme-point Symmetric Mode Decomposition,ESMD)方法。此方法是EMD的新发展,具有较强的自适应性和更强的局部特性,适合于非平稳、非线性信号的处理,可将原始信号中不同尺度的振荡或趋势分量逐级分离出来^[11]。传统盲源分离方法通常将观察信号置于时域中来进行分析,对于具有非平稳、非线性的机械故障信号,有必要借助时频分析(Time-Frequency Analysis,TFA)的优点,将盲源分离拓展到时频域^[3]。本文结合ESMD和基于时频分析的盲源分离方法的优点,提出利用ESMD实现观察信号的升维,将欠定条件下的盲源分离转换成恰定条件下的盲源分离问题;然后,根据机械设备故障信号非平稳的特点,采用平滑伪Wigner-Ville分布将新观察信号拓展到时频域估计源信号$\hat{s}$(t)。

1 极点对称模态分解

极点对称模态分解是经验模态分解的新发展,其采用直接插值法分解各个模态,相对于EMD的外包络线插值来说,其内部均线有着较小的幅值,能更好地降低由插值带来的幅值不确定性,在一定程度上抑制了端点效应^[9];同时,ESMD借用“最小二乘法”思想来优化最后剩余模态,使其成为整个信号序列分解过程中的“自适应全局均线”,以此来确定分解过程中的最优筛选次数,也在一定程度上抑制了EMD方法的“模态混叠”(或称“频率交叉”)^[11-12]。ESMD能将信号分解为数个表征信号特征时间尺度的本征模函数(Intrinsic Mode Function,IMF)和一条最优自适应全局均线(Adaptive Global Mean,AGM),其分解过程^[9]如下:

1) 找出数据Y(t)的所有极值点(极大值点和极小值点),并将它们依次记为E_i(i=1,2,…,n);

2) 对相邻极点用线段连接,并将线段中点依次记为F_i(i=1,2,…,n-1);

3) 通过一定方式补充左、右边界中点F₀和F_n;

4) 利用所获取的n+1个中点构造p条插值线L₁,L₂,…,L_p(p≤1) ,并计算它们的均值曲线L^*=(L₁+L₂+…+L_p)/ p;

5) 对Y-L^*重复步骤1) ~4) 直到|L^*|≤ε(ε是预先设定的容许误差)或筛选次数达到了预先设定的最大值K,此时分解出第一个模态M₁(t);

6) 对Y-M₁重复步骤1) ~5) 依次获得M₂(t),M₃(t),…,直到最后余量R(t)只剩一定数量的极点;

7) 让最大筛选次数K于整数区间[K_min,K_max]内变化并重复步骤1) ~6) 得到一系列分解结果,进而计算方差比率σ/σ₀,并画出其随K的变化图,其中σ和σ₀分别是Y(t)-R(t)的相对标准差和原始数据Y(t)的标准差;

8) 于区间[K_min,K_max]中选出对应最小方差比率σ/σ₀(意味着R(t)为数据的最优拟合曲线)的最大筛选次数K₀,据此重复步骤1) ~6) ,输出分解结果。

经过分解,原时间序列Y(t)可用式(1) 表示:

$Y(t)=\sum\limits_{i=1}^{n}{{{M}_{i}}(t)+R(t)}$

(1)

其中：Y(t)为时间序列；M_i(t)为第i个模态分量；R(t)为最优自适应全局均线；n为模态分量个数。

在机械故障诊断中,由于大多数设备所处环境复杂,故障信息常常掩埋在背景噪声之中。为了分析故障信息,经常需要先通过滤波来抑制噪声,再采用ESMD方法将信号从高频到低频依次分离出来,从而实现信号的升维。

2 源信号个数估计及最优观察信号选取

本文现以盲源分离的瞬时混合模型为研究对象。设s(t)=(s₁(t),s₂(t),…,s_N(t))^T,x(t)=(x₁(t),x₂(t),…,x_M(t))^T,其中：s(t)表示N维未知的相互独立源信号；x(t)表示M 维观测信号,且M≤N。令A为未知混合矩阵,则瞬时混合模型为：

$x\left( t \right)=As\left( t \right)$

(2)

盲源分离就是根据观测信号x(t)找到分离矩阵B,从而恢复源信号s(t)=Bx(t)。

2.1 源信号个数估计

在对已经进行预处理(滤波)的机械故障观察信号x(t)进行欠定盲源分离时,为了能够准确估计源信号,需要对源信号的个数进行估计。对观察信号x(t)进行ESMD,得到模态分量M_i(t)(i=1,2,…,n)和最优自适应全局均线R(t)。观察信号x(t)中的各个源信号以及噪声信号均包含在各个分解分量中,则利用上述子带信号组成的多维观测信号x_imf(t)=(M₁(t),M₂(t),…,M_n(t),R(t))^T预估源信号个数。

多维观测信号x_imf(t)的相关矩阵为:

${{R}_{x}}=E\left[ {{x}_{imf}}\left( t \right){{x}_{imf}}^{H}\left( t \right) \right]$

(3)

其中：x_imf^H(t)表示x_imf(t)的复数共轭。

当噪声是白噪声时,式(3) 可以改写为：

${{R}_{x}}=E\left[ s\left( t \right){{s}^{H}}\left( t \right) \right]+{{\sigma }^{2}}{{I}_{M-N}}$

(4)

其中：M为信号x_imf(t)的维数；N为R_x 主特征值个数；I_M-N为M-N维单位矩阵；σ²为噪声功率。

对相关矩阵R_x进行奇异值分解

${{R}_{x}}={{V}_{s}}{{\Lambda }_{s}}V_{S}^{T}+{{V}_{z}}{{\Lambda }_{z}}V_{z}^{T}$

(5)

其中：Λ_s=diag(λ₁,λ₂,…,λ_N),且λ₁≤ λ₂≤…≤λ_N;Λ_z是M-N个噪声特征值,且Λ_z=diag(λ_N+1,λ_N+2,…,λ_M)=σ²I_M-N。

理论上R_x的M-N个噪声特征值等于噪声功率σ²。因此,在假定精确估计协方差矩阵和噪声功率相对较小的前提下,可以通过判断R_x最小特征值的个数估计噪声子空间维数,从而确定源信号个数。但是在实际情况下,噪声特征值不可能完全相等,并且主特征值和噪声特征值之间的阈值选择较困难,往往不能准确预估源信号的个数。对此,可以采用贝叶斯信息准则(Bayesian Information Criterion,BIC)预估源信号与噪声子空间的维数。基于贝叶斯模型,Minka^[13]提出MIBS(Minaka Bayesian)选择模型,其能够有效估计源信号与噪声子空间维数。MIBS可用贝叶斯信息准则近似为：

$BIC(k)=(\prod\limits_{j=1}^{k}{{{\lambda }_{j}}{{)}^{-N/2}}}{{\tilde{\sigma }}_{k}}^{-N(l-k)/2}{{N}^{{{d}_{k}}}}$

(6)

其中：${{\tilde{\sigma }}^{2}}=(\sum\limits_{j=k+1}^{l}{{{\lambda }_{j}}})/(l-k)$；d_k=k(k-2l-1) /4,1≤k≤l,l为非零特征值个数。

通过寻找使代价函数BIC(k)最大的k=m,m即为估计的源信号个数。

2.2 最优观察信号选取

在估计出源信号个数之后,将混合信号x(t)和多维信号x_imf(t)组成新的多维观测信号x_{new_observe}(t)=(x(t),x_imf(t)^T)^T。需要从新的多维观测信号x_{new_observe}(t)中选取相同数目的最优观察信号,从而使分离达到最优效果。从新的多维观测信号x_{new_observe}(t)的结构可看出,其由两部分组成:混合信号x(t)和经过ESMD得到的各个分量{M₁(t),M₂(t),…,M_n(t),R(t)}。混合信号x(t)作为源信号的线性混叠,其包含了各个源信号信息,因此其应该得到保留。在此,需要从M₁(t),M₂(t),…,M_n(t),R(t)中选择m-1个最优观测信号。

为了找到最优观测信号,需要找到与源信号关系比较大的观测信号。对此,可以根据M₁(t),M₂(t),…,M_n(t),R(t)与包含所有源信号信息的x(t)的相关系数来选择最优观察信号。选择相关系数最大的前m-1个模态分量与x(t)构成多通道观察信号x_mul(t),再采用基于时频分析的盲分离方法估计源信号。

3 基于时频分析的盲源分离方法

机械设备故障信号具有非平稳特征,而时频分析提供了时间域和频率域的联合分布信息,清楚地描述了信号频率随时间变化的关系,成为处理非平稳信号的有力工具。引入平滑伪Wigner-Ville分布,将盲源分离方法由传统的时域分析拓展到时频域,能够有效刻画机械设备故障信号中的各个故障特征分量,从而更加准确地估计源信号。

对于上述处理后得到的多通道观察信号x_mul(t),其基于时频分析的盲源分离算法步骤^{[3, 14]}如下:

1) 计算观察信号的协方差矩阵R_xx(τ)。

2) 计算协方差矩阵R_xx(τ)按降序排列的特征值λ₁,λ₂,…,λ_m和相应的特征矢量h₁,h₂,…,h_m,并估计噪声方差${{\hat{\sigma }}^{2}}$和白化矩$\hat{W}$：

${{\hat{\sigma }}^{2}}=\frac{1}{m-N}\sum\limits_{i=N+1}^{m}{{{\lambda }_{i}}}$

(7)

$\begin{align} & \hat{W}=[{{({{\lambda }_{1}}-{{{\hat{\sigma }}}^{2}})}^{-1/2}}{{h}_{1}},{{({{\lambda }_{2}}-{{{\hat{\sigma }}}^{2}})}^{-1/2}}{{h}_{2}},\cdots \\ & ,{{({{\lambda }_{N}}-{{{\hat{\sigma }}}^{2}})}^{-1/2}}{{h}_{N}}{{]}^{H}} \\ \end{align}$

(8)

对式(2) 左乘一个白化矩阵$\hat{W}$,得到白化后的观测信号：

$z(t)=\hat{W}{{x}_{mul}}(t)=\hat{W}As(t)=Us(t)$

(9)

由式(9) 可以看出,白化信号是源信号的一个“酉矩阵混合”。

3) 同时对式(2) 左右两端进行时频变换,得到源信号和观测信号的时频分布矩阵关系：

${{D}_{xx}}(t,f)=A{{D}_{ss}}(t,f){{A}^{H}}$

(10)

其中：D_xx(t,f)为观测信号的时频分布；D_ss(t,f)为源信号的时频分布。将D_xx(t,f)左右分别乘以$\hat{W}$,得到白化的时频分布矩阵：

${{D}_{zz}}(t,f)=\hat{W}A{{D}_{ss}}(t,f){{A}^{H}}{{\hat{W}}^{H}}$

(11)

4) 按一定规则选择K个自项,并对自项位置处的时频分布矩阵$\{D_{zz}^{\varphi }({{t}_{i}},{{f}_{i}}),i=1,2,\cdots ,K\}$进行联合对角化,估计酉矩阵$\hat{U}$。

5) 根据上述估计的白化矩阵$\hat{W}$和酉矩阵$\hat{U}$,可以得到源信号的估计信号为：

$\hat{s}(t)={{\hat{U}}^{H}}\hat{W}{{x}_{mul}}(t)$

(12)

4 仿真信号分析

为了验证本文算法的有效性,使用仿真数据模拟机械设备故障进行分析。仿真源信号如式(13) ,采样点个数为512,采样频率为1000Hz。

$\left\{ \begin{matrix} {{s}_{1}}(t)=2sin(2\pi {{f}_{1}}t+10) \\ {{s}_{2}}(t)=(1+sin(2\pi {{f}_{2}}t))sin(2\pi {{f}_{3}}t) \\ {{s}_{3}}(t)=2cos(2\pi {{f}_{4}}t-10) \\ \end{matrix} \right.$

(13)

其中: f₁=20Hz, f₂=10Hz, f₃=100Hz, f₄=50Hz。源信号时域波形图和频谱图如图 1所示。

在机械设备故障诊断中,传感器采集到振动数据大多是由多个源信号混合而成。为了模拟这一现象,随机选取一个3×3的混合矩阵A,根据式(2) 得到三个混合信号。同时,为了模拟真实数据中的噪声,向三个混合信号添加随机噪声,使得最终混合信号x(t)的信噪比达到20dB。混合矩阵A如式(14) 所示,混合信号的时域波形图和频谱图如图 2所示。

$A=\left[ \begin{matrix} -0.5081 & 0.5110 & 0.5285 \\ -0.2767 & 0.1712 & -0.9426 \\ 0.4015 & -0.4222 & -0.5161 \\ \end{matrix} \right]$

(14)

选择第一个混合信号作为本文的观察信号x₁(t),为了能够成功分离出各个源信号,需要估计源信号的个数,然后选取最优观察信号,将欠定盲分离问题转换为恰定盲分离问题。对观察信号x₁(t)进行ESMD,得到M₁(t)、M₂(t)、M₃(t)、M₄(t)、M₅(t)等5个模态和一个最优自适应全局均线R(t),各个分量如图 3所示。

将ESMD各个分量组成多维观测信号x_imf(t)=(M₁(t),M₂(t),…,M_n(t),R(t))^T,计算其相关矩阵R_x。再对相关矩阵R_x进行奇异值分解,特征值Λ=diag(λ₁,λ₂,…,λ₆)由大到小排列,见表 1。通过特征值个数判断源信号个数为3个,进一步由BIC计算获得源信号个数为3,同时前三个特征值的累积占比达到96.21%。

为了选取最优观察信号,各个模态与混合信号x₁(t)的相关系数见表 2。从表 2中可知M₁(t)与M₂(t)相关系数最大,则新的观察信号为x_{new_observe}(t)=(x₁(t),M₁(t),M₂(t))^T,新观察信号的时域波形图及时频分布如图 4~5所示。

根据获得的观察信号,采用基于时频分析的盲源分离方法估计源信号,估计源信号时域图和频域图如图 6所示。图 7是采用EMD混合信号x₁(t)而得到的观察信号,基于EMD和时频分析的盲源分离算法(BSS algorithm based on EMD and TFA,EMD-TFA-BSS)的结果的时域波形图和频谱图如图 8所示。对比图 4和图 7可以发现,图 4和图 7后两个观察信号都存在一定的模态混叠和幅值畸变。但是,基于ESMD方法得到的观察信号要优于基于EMD方法得到的观察信号,这主要是由于ESMD方法采用直接插值法分解各个模态,相对于EMD的外包络线插值来说,其内部均线有着较小的幅值,能更好地降低由插值带来的幅值不确定性;同时,ESMD也在一定程度上抑制了模态混叠。通过对比图 6和图 8的估计源信号发现,两种方法得到的估计源信号在幅值上都与源信号存在差异,同时估计出的源信号顺序也都发生了变化,这恰好说明了盲源分离的“两个不确定性”。图 6和图 8的时域波形图显示,基于ESMD方法的估计源信号明显优于基于EMD方法的估计源信号;从频域图也可以看出,基于EMD方法的估计源信号在主频附近也存在着幅值较大的其他频率信号。

为了定量比较两种算法,采用相关性评价标准对上述两种算法进行评价。两信号相关系数定义^[8]如下:

$\rho ({{\hat{s}}_{i}},{{s}_{j}})=\frac{|cov({{{\hat{s}}}_{i}},{{s}_{j}})|}{\sqrt{cov({{{\hat{s}}}_{i}},{{{\hat{s}}}_{i}})cov({{s}_{j}},{{s}_{j}})}}$

(15)

其中：cov为协方差；${{\hat{s}}_{i}}$表示第i个估计源信号；s_j表示与${{\hat{s}}_{i}}$对应的源信号；0≤ρ(${{\hat{s}}_{i}}$,s_j)≤1,ρ(${{\hat{s}}_{i}}$,s_j)越接近1,说明估计源信号和源信号越相似。两种盲分离算法的相关系数如表 3所示。表 3中相关系数按照源信号的顺序排列。与EMD-TFA-BSS方法相比,基于ESMD和时频分析的盲源分离算法(ESMD-TFA-BSS)的估计信号与仿真信号的相关系数分别提高了12.35%、0.80%、13.00%,ESMD-TFA-BSS方法分离精度明显优于EMD-TFA-BSS方法。

5 实验分析

为了验证本文算法在工程实践中的有效性,以美国西储大学滚动轴承点蚀故障数据为研究对象。在滚动轴承故障实验中,滚动轴承型号为6205-2RS JEM SKF深沟球轴承,该轴承滚珠个数为9个,滚动体直径为7.9mm,轴承节径为39mm。电机转速为1797r/min,点蚀故障直径为0.1778mm,振动信号采样频率为12kHz。

轴承在发生内圈故障或外圈故障时,滚动体每次经过故障点都会产生一个周期性的冲击,对于内圈故障,其故障频率可以用式(16) 表示：

${{f}_{inner}}=\frac{r}{60}\times \frac{m}{2}\times (1+\frac{d}{D}cos(\alpha ))$

(16)

对于外圈故障,其故障频率可以用式(17) 表示

${{f}_{out}}=\frac{r}{60}\times \frac{m}{2}\times (1-\frac{d}{D}cos(\alpha ))$

(17)

其中：r为电机转速；m为滚动体个数；d为滚动体直径；D为轴承节径；α为接触角,且α=0°。

按照式(16) 和式(17) 计算可得,轴承内圈故障频率理论值为162.1Hz,轴承外圈故障频率理论值为107.5Hz。将内圈故障数据和外圈故障数据经过一个2×2的混合矩阵A混合,得到两路混合观测信号。选择其中的一路作为观测信号,观测信号时域波形及包络频谱如图 9所示。

为了分离出源信号,首先采用ESMD方法分解单通道观测信号,利用BIC估计源信号个数,同时采用相关系数法选出最优观测信号,完成观测信号的升维。在将欠定盲分离问题转换成恰定盲分离问题后,采用基于时频分析的盲源分离方法,从新的观测信号中估计源信号,估计源信号的时域波形以及包络频谱如图 10所示。由于分离信号成分比较复杂,单从分离信号时域波形中并不能确定轴承的故障类型,需要借助分离信号的包络频谱图进行分析。从图 10第一个包络频谱图中可以明显看到主峰值出现在164.1Hz处,

与轴承内圈故障频率理论值162.1Hz非常接近,2Hz的差别可能是由于电机转速无法恒定在1797r/min,而是在1797r/min附近波动引起的^[3]。第二个包络频谱图分别在105.5Hz、210.9Hz、539.1Hz处出现了峰值,最大幅值出现在105.5Hz处,这与轴承外圈故障频率理论值107.5Hz较接近。显然,210.9Hz和539.1Hz分别对应着轴承外圈故障二倍频和五倍频。因此,可以初步判断轴承发生了内圈故障和外圈故障。

为了对比,采用EMD-TFA-BSS方法对上述单通道观测信号进行盲源分离,分离信号的包络频谱图如图 11所示。从包络频谱图中也可以得出轴承发生了内圈故障和外圈故障的结论,但是可以发现EMD-TFA-BSS方法得到的分离信号1的包络频谱图在0~100Hz、500~700Hz均出现较大的幅值,在一定程度上会干扰对故障类型的判断。出现这种情况有可能是由于单通道观测信号在采用EMD方法分解时,产生了较ESMD方法严重的“模态混叠”。

6 结语

针对实际故障诊断中,现场传感器个数往往小于源信号的个数,传统盲源分离方法通常难以实现源信号有效分离的问题,本文提出采用ESMD方法分解观测信号,实现单通道观察信号的升维,将欠定盲分离问题转换成比较容易解决的恰定信号盲分离问题。考虑到实际机械设备故障中振动信号的幅值和频率都随时间变化,具有非平稳特性,本文结合时频分析在非平稳信号分析中的优点提出在时频域下进行恰定信号盲源分离,仿真和实验研究结果证明了该方法的有效性和可行性。

本文主要考虑了观测信号中各个源信号占比大致相同的情况,对于观测信号中某些源信号占比较小的情况将是下一步研究的重点。