基于改进S变换和ICA的相关源分离方法*

韦成龙 ， 周以齐 ， 李 瑞 ， 于 刚

(1.山东大学机械工程学院 济南，250061) (2. 济南大学自动化与电气工程学院 济南，250061)

引 言

机械系统运行时各个部件同时工作，由于受到不同振动源的激励，测试采集的振动信号往往是多源混合的结果。如何从复杂的观测信号中快速准确地分离出源信号，对于振源识别和贡献量估计具有重要意义。独立成分分析作为盲源分离的常用方法，在源未知的情况下，仅利用观测信号和先验知识分离和提取信号源的过程，广泛应用于噪声控制、振源识别和故障诊断等领域[1-2]。

ICA以源信号相互统计独立作为假设条件，基于非高斯性最大化原理分离观测变量，而在实际工程应用中，严格的统计独立往往过于理想化。尤其对于动力机械来说，振源众多且传递路径复杂，发动机作为主要的动力源，振动信号具有明显的冲击和时变特征，且各部件级联耦合，观测信号中独立成分和相关成分混合在一起，难以获得理想的分离效果[3-4]。当源信号存在统计相关成分时，分离获得的估计源与真实源之间将相差一个与源向量协方差有关的因子向量[5]。

在解决相关性源信号分离的问题上，国内外学者开展了相关研究，主要方法有子带分解[6-7]、核典型相关[8]和稀疏成分分析[9-10]等。文献[7]提出基于小波包分解的相关源分离方法，利用互信息作为独立性测度重构观测信号，获得分离矩阵，但分解性能受先验知识和小波分解层数的影响。文献[8]采用核典型相关分析的方法将非线性信号映射到核特征空间中，将非线性问题转化为线性，解决非线性相关信号的分离，而核化过程并不能有效消除相关项。文献[9]利用稀疏信号延混合矩阵列方向向量的线性聚类特征，能够实现源数欠定条件下的估计，当信源的稀疏度不够时，分离误差较大。

笔者借鉴子带分解和稀疏分量分析的思想，以相关源中部分成分统计独立为前提假设，提出了一种基于改进S变换和ICA的相关源分离方法。在S变换[11-12]中引入参数可调的尺度因子，将时域观测信号扩展到时频域，利用混合观测信号在时频域中实部和虚部的方向向量，检测并剔除源信号中的相关成分，避免了独立子带方法中信号分解层数的确定和互信息在同频成分检测上的失效问题，通过重构观测信号以满足盲源分离的统计独立条件，进而采用负熵最大化的独立成分分析获得分离矩阵，实现相关源的有效分离。

1 分离模型

盲源分离是在混合系统未知的情况下仅利用观测信号恢复源信号的过程。假设正定或超定条件下，即n个源信号瞬时混合被m个传感器接收(n≤m)，将观测信号中的多变量数据简化为线性混合问题，其数学模型可表示为

x(t)=As(t)+n(t)

(1)

其中：A∈Rm×n为列满秩混合矩阵；x=(x1,x2,…,xm)T为观测变量；s=(s1,s2,…,sn)T为n维未知的信号源；n(t)为附加噪声。

观测变量x(t)作为信号源s(t)的线性加权，估计源信号y(t)就转化为寻找分离矩阵W的过程，即

y(t)=Wx(t)

(2)

当激励源存在同频分量时，信号之间统计相关，无法采用以独立源作为假设前提的标准ICA方法直接分离出来。假设相关源由独立成分sD(t)和非独立成分sID(t)两部分组成

s(t)=sD(t)+sID(t)

(3)

观测信号可看作独立成分和相关成分的线性叠加，在瞬时混合模型下，各成分的混合系数相同，分离矩阵W可通过部分独立成分的混合信息获得。

yD(t)=Wxrec(t)=WAsD(t)

(4)

其中：xrec(t)为独立成分重构的观测信号；yD(t)为剔除相关项后的估计源。

笔者通过时频预处理方法，从混合变量中剔除源信号的相关成分，对剩余信号进行重构，从而满足统计独立的条件。对重构后的信号进行盲分离，估计出分离矩阵W，由式(2)恢复相关源信号y(t)。

2 算法原理

2.1 改进S变换

S变换作为一种短时傅里叶变换的可变窗扩展和连续小波的相位修正，能够有效克服窗宽固定的缺点，时频变换后各分量的相位与原始信号保持直接联系，对工程测试中非平稳信号具有较好的时频特性[13-15]。标准S变换是在短时傅里叶变换的基础上引入时宽与频率倒数成正比关系的高斯窗，定义为

其中：ST(τ,f)为信号x(t)的S变换；g(τ-t,f)为窗函数；τ为平移因子，用于控制高斯窗在时间轴t上的位置；f为信号频率。

高斯窗函数中标准差被改造为关于频率f的倒数，然而标准S变换的窗宽变化趋势固定，时频谱唯一。为了提高其适应性，修正窗函数模式固定的缺陷，引入尺度因子λ,α,β，得到改进S变换的表达式为

(7)

其逆变换(inverse modified S-transform，简称 IMST)表示为

(8)

图1为引入尺度因子后，窗函数半高时宽随频率的变化曲线。当λ=1,α=1,β=0时，可以看作标准S变换。当0<α<1时，时窗宽度的减小速率降低，高频部分频率分辨率提高，而α>1时，通常用于捕捉高频瞬态信号。当0<λ<1时，频率分辨率提高，窗宽呈非线性变化，趋于平缓，λ>1时，则相反。β通常取正值，用于提高低频区域的时间分辨率。λ,α,β均具有控制窗口宽度和变化率的作用，可针对实际需求优化三者参数，以获得最优时频分辨率。

图1 窗函数半高时宽随频率变化曲线Fig.1 The curve of temporal full width at half maximum versus frequency

2.2 相关成分检测

观测信号x(t)经MST得到一个关于时间t和频率f的二维复时频矩阵x(t,f),aj=[a1j,a2j,…,amj]T为混合矩阵A的第j列向量，则瞬时混合模型在时频域中的表达式为

(9)

仅存在两个激励源s1,s2时，式(9)可写作

x(t,f)=a1s1(t,f)+a2s2(t,f)

(10)

满足统计独立的条件下，混合矩阵中的每个时频点仅为其中一个独立源的响应。此时，假设在时间tp和频率fq处，s1(tp,fq)≠0，s2(tp,fq)=0，式(10)可表示为

x(tp,fq)=a1s1(tp,fq)

(11)

其实部和虚部变换分别为

Re{x(tp,fq)}=a1Re{s1(tp,fq)}

(12)

Im{x(tp,fq)}=a1Im{s1(tp,fq)}

(13)

由式(12)，(13)可知，x(tp,fq)实部和虚部的方向向量与混合矩阵的列向量方向一致。当两激励源在时间tp和频率fq处频率相同，即存在相关成分时，源信号s1(tp,fq)≠0，s2(tp,fq)≠0，对应观测信号在时频域中实部和虚部变换为

Re{x(tp,fq)}=a1Re{s1(tp,fq)}+a2Re{s2(tp,fq)}

(14)

Im{x(tp,fq)}=a1Im{s1(tp,fq)}+a2Im{s2(tp,fq)}

(15)

此时，只有在式(16)成立，相关成分的相位相差0°或180°的情况下，x(tp,fq)实部和虚部的列向量方向一致。

(16)

实际激励源中，相关成分的相位往往是随机的，同相位的概率相当低，所对应的观测信号在复时频矩阵中的实部和虚部向量会存在一个角度差Δθ。将其扩展为m×n维混合信号，时频域中实部和虚部向量之间的夹角可表示为

(17)

式(17)表明，当cosΔθ=1时，实部和虚部向量方向一致，观测信号中对应成分独立，而cosΔθ<1时，对应成分相关。考虑到工程测试中测点位置和噪声的影响，很难严格满足条件，通过设置阈值ε来降低其影响，ε取接近于1的值，当cosΔθ<ε时，检测为相关成分并将其剔除。阈值的确定需要根据实测数据来选择，至少保留10%以上的数据用来保证分离矩阵的估计。

图2 相关源分离算法流程图Fig.2 Correlated source separation flow chart

2.3 相关源估计

基于MST和ICA的相关源分离算法流程如图2所示。分离出相关成分后，利用IMST重构剩余信号，从而满足ICA对信号源统计独立的要求，然后以最大化负熵作为独立性判据，采用快速固定点算法[1](fast fixed-point independent component analysis，简称 FastICA)获得分离矩阵。负熵作为一种非高斯性的度量，较峭度具有更好的鲁棒性和统计特性，以非二次函数近似的形式估计负熵[16]

(18)

其中：y为标准化输入变量；υ为具有零均值和单位方差的高斯随机变量；G为非二次函数。

采用渐进方差的形式，选取最优非二次函数G(y)

(19)

其中：1≤α≤2时，取G(y)=log2cosh(y)；α<1，即超高斯程度很高时，取G(y)=-exp(-y2/2)。

通过极大化E{G(WTz)}获得负熵的最大近似值，根据牛顿迭代更新分离矩阵W，如式(20)所示，每轮迭代对W进行正交化和标准化，直到收敛为止。

W+=E{zG(WTz)}-E{zG′(WTz)}W

(20)

(21)

3 仿真分析

为了验证提出方法的有效性，构造4个仿真信号模拟机械振源。s1，s2分别为20 Hz和120 Hz谐波信号，模拟中低频简谐振动，s3为调幅信号，模拟幅值调制特征；s4为指数衰减信号，模拟周期性冲击。在s1～s4中分别加入相位随机的40 Hz谐波成分sd。设采样频率为1 024 Hz，采样时间为1 s，信号函数为

(22)

图3为加入40 Hz同频相关成分后源信号的时域波形。采用皮尔森相关系数衡量信号之间相似程度，假设两个随机变量x和y，相关系数ρ定义为

图3 源信号时域波形Fig.3 Waveform of source signals

(23)

ρ的绝对值越接近于1，说明两变量之间波形的相似度越高，相关性也越强。

计算仿真源之间相关系数如表1所示。加入同频成分后，源信号呈现不同程度的相关性。随机生成(0,1)均匀分布的4×4混合矩阵A，以线性叠加的方式混合相关源得到4个观测信号x1～x4，其时域波形如图4所示。

表1 源信号相关系数Tab.1 Source signal correlation coefficient

图4 观测信号时域波形Fig.4 Waveform of observed signals

为使待分离信号满足统计独立条件，去除相关成分的干扰，首先利用MST提取时频系数，取尺度参数α=0.3，β=5，λ=1，计算观测信号在各时频点上实部和虚部向量的夹角Δθ，设置阈值ε=0.9，将cosΔθ<ε成分置零，再经IMST重构剩余信号。图5,6分别为观测信号x1和重构信号xrec1的MST时频谱。对比看出，经本方法处理后，观测信号中的40 Hz相关成分得到有效去除。然后采用ICA分离重构信号，利用分离矩阵得到估计源y1～y4。

图5 观测信号x1时频谱Fig.5 Time-frequency spectrum of observed signal x1

图6 重构信号xrec1时频谱Fig.6 Time-frequency spectrum of reconstruction signal xxec1

图7 分离信号独立分量时域波形Fig.7 Waveform of separated independent components

图8 估计源时域波形Fig.8 Waveform of estimated source signals

图7为分离重构信号得到独立分量yD1～yD4的时域波形。可以看出，分离源中20 Hz和120 Hz的中低频谐波成分，调幅成分以及瞬态冲击成分被有效分离出来，可用于相关源信号中独立特征的提取。图8为本研究方法得到的估计源y1～y4的时域波形。对比图3，波形上与加入40 Hz相关成分后的仿真源基本一致，并计算分离信号与仿真源之间的相关系数，如表2所示。y1，y2，y3，y4分别与仿真源s4，s2，s3，s1的相关系数为1，说明在无噪条件下，几乎可以实现相关源信号的完全分离与估计。

表2 分离信号相关系数Tab.2 Separated signal correlation coefficient

相关系数仅用来衡量两两变量间的相似程度，为了进一步评价算法的分离性能，采用Amari性能指数[18](performance index，简称 PI)来度量混合矩阵A和分离矩阵W之间的差异，定义为

(24)

其中：gij作为全局矩阵G=WA的第(i,j)个元素，PI越接近于零，说明算法分离性能越好，当超过0.2时，可认为分离失效。

在x1～x4中加入信噪比为5～50 dB的高斯噪声，对比本分离方法与直接采用FastICA分离时的PI值，如图9所示。可以看出，随着信噪比的降低，算法的分解性能变差，而本研究方法的PI指数明显小于直接使用FastICA的分解结果，在信噪比20 dB以上的有噪条件下，对相关源混合信号仍能够取得较好的分离效果。

图9 Amari性能指数对比Fig.9 Comparison of Amari performance index

4 实例分析

工程机械中观测变量通常为多激励源的综合响应，以某型挖掘机为测试对象，发动机作为主要动力源，振动信号经4个减振装置传递到车架上并混合，将每个减振点看作一个激励，利用车架振动信号估计振源。测试采用8通道NI数据采集系统和压电式加速度振动传感器，在空挡最大油门工况下，同时测量四点激励与上车架的振动信号。发动机转速为2 280 r/min，采样频率为5 120 Hz,实验样车和上车架传感器测点布置如图10(a)和图10(b)所示。

图10 挖掘机测试图Fig.10 Experiment of excavator

由于车架振动来源于发动机和减振器的共同作用，激励信号之间存在不同程度的相关性，四点减振处振动信号的相关系数最大为0.36。利用本研究方法分离车架振动信号，取尺度参数α=0.3，β=56，λ=0.8，设置阈值ε=0.9，去除相关项后重构独立混合信号，最后利用ICA获得分离矩阵W,估计出振源信号，确定各个振源对车架振动的贡献量。

图11 上车架振动信号Fig.11 Upper frame vibration signals

图12 分离信号Fig.12 Separated signals

图13 分离信号时频谱Fig.13 Time-frequency spectrum of separated signals

图11，图12分别为车架振动信号x1～x4和本研究方法分离得到的估计源y1～y4的时域波形和频谱图。图13为分离信号的MST时频谱。可以看出，4个分离信号的能量主要在1 000 Hz以下，以76 Hz和380 Hz频率成分为主要峰值，对应发动机发火基频和10阶转频，y1中600 Hz附近周期性冲击成分也较为明显，应为缸体内气体燃烧产生的振荡冲击引起。

分离信号与激励信号之间的相关系数如表3所示。估计源y1～y4分别与左后减振、左前减振、右前减振、右后减振处的振动信号相关系数达到0.8以上，为显著相关，说明分离出来信号与4个激励点处振动信号的波形基本一致。本研究方法能够从车架混合信号中成功分离出激励源，对于工程测试信号，依然能够取得较好的分离效果。

利用分离矩阵与估计信号定量计算各个振源对车架振动的贡献比，如表4所示。可以看出，右侧前后减振处激励对车架振动贡献较大，其中，右前减振作为第1主振源，贡献比为40.4%。结合车架振动以及分离信号的谱分析结果推测，右侧车架的380 Hz峰值频率成分应主要由右前减振处振动激励引起，可通过针对性的调节右前减振装置、优化车架结构，来达到改善车架振动的目的。

表3 分离信号相关系数Tab.3 Separated signal correlation coefficient

表4 振源贡献百分比Tab.4 Percentage of vibration source contribution

5 结束语

针对传统盲源分离方法对于统计独立条件的限制，提出了一种基于改进S变换和独立成分分析的相关源分离方法，用来解决工程应用中相关性机械振源在混合观测信号中的分离和估计问题。该方法在线性瞬时混合盲源分离模型的基础上，针对独立成分和相关成分共存的信号源，以相关成分相位随机为假设，利用其在时频域中的方向向量，有效剔除混合信号中的相关分量，避免了相关成分对分离结果的影响。在时频化处理中，以S变换为基础引入窗宽可调的尺度因子，改善时频分辨率固定的缺陷，提高了复杂工程测试信号的适应性。通过仿真和实测分析验证了提出方法在相关源分离中的有效性，分离性能优于未对相关项进行预处理的盲源分离结果，通过阈值的合理设置，能够从混合信号准确分离出相关源信号，得到各个振源对混合信号的贡献比。该方法适用于正定、超定条件下含有交叉频率成分的相关源分离，在信源识别和故障诊断等领域具有一定的工程应用价值。