基于稀疏特征的连轧机故障信号分离方法*

2019-12-31严保康周凤星宁博文李维刚

振动、测试与诊断 2019年6期

严保康，周凤星，宁博文，李维刚

(武汉科技大学信息科学与工程学院武汉，430081)

1 问题的引出

在冶金轧制过程中，连轧机机组起着关键的作用，直接决定产品的质量精度和性能指标。连轧机机组在长期连续、高强度工作中，经常会出现轧辊表层裂纹、轧辊表层局部剥落、轧辊支撑轴承失效等故障，进而影响产品质量精度，甚至导致设备停产及重大人员伤亡等事故，因此，连轧机的状态监测和故障诊断对保证生产稳定性、提高产品质量、降低生产成本具有重要意义[1-3]。图1所示为某大型钢厂的精轧连轧机示意图，该机组共有7台精轧机组成(F1～F7)，每台精轧机核心部件主要包括工作辊和支撑辊，其中工作辊直接与坯料接触使其变形，支撑辊主要是增加工作辊的刚度，防止工作辊轧制时断辊及过大的弹跳，在工作辊和支撑辊两端均有轴承座进行支撑。每台轧机的主要故障点在工作辊表层、支撑辊表层以及8个轴承座，而通常每台轧机只布置4个振动传感器，分别布置在上、下支撑辊和上、下工作辊的轴承座上，在轧机出现复合故障或连轧机组其他轧机出现故障时，由于振动在材料中的传递，使得所采集的振动数据均为不同部件振动信号强耦合的混合数据，难以直接根据振动数据进行故障分析[4-5]。因此，连轧机故障信号的盲分离对故障的精准定位和诊断具有重大意义。

图1 精轧连轧机示意图

盲源分离方法仅根据测量信号可得到混合模型及各源信号，在理论上为多源振动分析提供了强有力的技术手段。基于信号高阶统计特性的独立分量分析方法在信号处理和特征提取中有着广泛的应用，但是在工程应用中，由于强噪声的干扰及各振源交叉频率成分的存在，导致独立成分分析算法(independent component correlation algorithm，简称ICA)无法满足独立性假设而失效。

李纪永等[6]针对线性时频相交混合信号，提出白化-均匀加权非正交联合对角化方法估计混合矩阵，能分离出断齿与磨损复合故障源。Bofill等[7]针对具有6个信号源稀疏信号、2个传感器的盲分离问题，提出了势函数方法，但是参数的设定缺乏理论指导，难以适用到更多的传感器情形。何昭水等[8]克服了Bofill方法的局限性，建立了K-聚类与主成分分析相结合的方法，可以方便地应用于3个及以上传感器的情况，且实现简单、混叠矩阵估计精度高。于刚等[9]针对相关机械源盲分离问题，提出了一种基于信号稀疏特性的盲源分离方法，使得相关成分被有效剔除，从而准确地估计出混合矩阵。Sadhu等[10]通过小波变换推断出全局结构模态信息，再通过主成分分析，达到盲源分离的目的。

笔者针对连轧机故障的强耦合特征和盲源分离过程中噪声对分离效果的影响问题，提出了一种基于稀疏特征的连轧机故障信号分离方法，该方法结合稀疏分解方法对复杂非平稳信号有精确表示能力，对混合信号进行稀疏表示，提取出混合信号中的微弱故障特征信息，再通过稀疏原子的相似性对原子进行聚类分析，确定混合信号的稀疏矩阵和盲源个数，实现连轧机混合故障信号的盲分离。

2 基于稀疏分量的盲分离算法

2.1 盲源分离的基本概念

盲源分离(blind source separation, 简称BSS)是指在源信号和传输信道参数未知的情况下，仅通过一系列传感器测量的信号分离出各源信号的一种信号处理方法。BSS的数学模型可以描述为

x(t)=As(t) (t=1,2,…,T)

(1)

其中：x(t)=(x1(t),x2(t),…,xm(t))T为通过传感器测量的m维向量；s(t)=(s1(t),s2(t),…,sn(t))T为n个未知源构成的n维向量；A为一个m×n的混叠矩阵。

盲分离的目的就是在相差一个排列矩阵和尺度矩阵的意义上使得分离信号与源信号的波形保持一致。在分离过程中，混叠矩阵的先验知识是未知的，需要在源信号统计独立的前提下进行估算[7]。

2.2 基于稀疏分量的盲分离基本原理

稀疏性信号是指信号在很小部分的时刻取值较大，而在很多时刻取值为零或者接近于零。当源信号具有较好的稀疏性时，观测信号在同一时刻出现2个或以上稀疏源信号的可能性很小[9]。对于属于同一源信号的较大幅值点正比于观测信号中的混合系数，因此，可以用信号的稀疏特征估计源信号的混合矩阵。

n个源信号混合的m个观测信号，可以展开为

(2)

其中：si(t)为第i个源信号；xj(t)为第j个观测信号；aij为第j个源信号到第i个观测信号的混合系数[8]。

假设在某一时刻点t，只有源信号si(t)取值非零，其他源信号取值很小或接近于零，则式(2)可近似地写成

(3)

进一步写成

(4)

因此，在信号具有稀疏性时，某一个源信号si(t)取值占优时的时刻可以确定一条直线，这条直线的方向为混合矩阵的第i个列向量(a1i,a2i,…,ami)T。这样，每个源信号将分别确定一条直线，n维源信号将确定n条直线。通过混合信号散点图中的直线方向，即可确定混合矩阵列向量方向，因此，通常采用“K-均值聚类”方法估计混合矩阵[9]。然而因为噪声的干扰，“K-均值聚类”方法对于直线方向的估计准确性和鲁棒性不好，无法达到对混合矩阵的良好估计。

3 基于稀疏特征的连轧机信号分离方法

连轧机设备在轴承或轧辊出现裂纹、点蚀等故障时，会产生冲击信号。由于设备的复杂性和轧制现场的恶劣性，冲击信号呈现以下特点：a.微弱的冲击信号淹没在强背景噪声下；b.冲击信号具有稀疏性；c.信号具有强耦合性，且振源位置通常超过3处。如何将耦合在一起的微弱冲击成分提取并分离出来是连轧机故障信号分离的关键。笔者提出的基于稀疏特征的连轧机信号分离方法的流程如图2所示。

图2 基于稀疏特征的连轧机信号分离方法流程图

3.1 基于时频谱的稀疏特征提取方法

广义S变换(generalized S-transform, 简称GST)实质上是加窗的傅里叶变换，窗的宽度是随着频率而变化的，频率越高，宽度越窄；频率越低，宽度越宽，从而具备类似多分辨率的特性[11-12]。GST的正反变换定义表示如下

其中：λ和p为调节参数，且λ>0,p>0。

由于冲击类故障信号中的冲击成分是振荡衰减信号，且具有稀疏性，因此在滑动窗进入冲击信号和离开冲击信号的整个区间，都会与冲击信号进行有效的内积运算，从而在时频谱上呈现出一定强度的幅值；而在其他非冲击信号区间，则不会与冲击信号进行有效的内积运算，从而在时频谱上呈现出微弱的幅值。因此，在冲击类故障信号广义S变换的时频谱的特点主要有：a.在冲击信号的振荡频率周围会呈现出能量集中区域；b.该能量集中区域会分布在每一个冲击信号的时频范围内。

由Heisenberg测不准原理可知，为了兼顾时间分辨率和频率分辨率，让信号在低频段的时窗较宽，从而获得较高的频率分辨率；而在高频段的时窗较窄，从而获得较高的时间分辨率[13-15]。

对于一个稀疏性的冲击信号，其GST频谱随着调节参数的变化而变化，但是能量主要集中在冲击出现的时频附近，且能量集中区域的形状各不相同。随着调节参数的增大，能量集中区域在时间轴上越来越窄，在频率轴上越来越宽。在与故障信号一致的调节因子下的时频谱图，其能量最为集中。

为获取稀疏分解过程中的各个最佳原子，需要将多分辨率下GST时频谱中各区域的最大值点提取出来，从而构建稀疏表示的最佳原子序列[18]。

基于时频谱分割的稀疏特征提取方法算法流程如图3所示。

图3 基于时频谱分割的稀疏特征提取方法算法流程

经过GST所获取不同调节参数下的时频谱图，尽管能量集中区域的时频尺寸是变化的，但是能量集中区域的极值点位置都是分布在冲击信号时频点附近。因此，需要对所有时频谱图做以下处理：a.二值化处理成二值时频谱；b.通过图像融合处理生成最佳二值时频谱；c.通过图像分割分离出最佳二值时频谱中各独立的连通域。

1) 二值化处理。二值化处理是将多尺度所有时频谱图通过适当的阈值选取而获得可以反映局部特征的二值化图像，时频谱的二值化处理使得数据量大幅减少，有利于对时频谱的局部特征分割。

阈值的选取对生成二值图像十分关键，阈值太大，尽管很大程度上抑制了噪声，但是也会过滤掉部分微弱的冲击成分；而阈值太小，则引入了大量噪声，同时增加后续图像分割的计算量，并降低了算法的稀疏度和精度。因此，二值化处理的阈值满足以下条件

(6)

其中：a为时频谱对应的数据矩阵；N为矩阵a的数据总个数[16-17]。

二值化处理的数学表达式为

(7)

其中：i为矩阵a的行索引值；j为矩阵a的列索引值。

2) 二值时频谱融合。根据时频谱图的能量分布性质，将同一个冲击成分的所有调节参数下的二值时频谱图进行融合，可得到能量更加集中的区域。笔者所采用的为交集融合，最终形成最佳二值时频谱图。

3) 二值时频谱分离。最佳二值时频谱中包含若干个连通域，为提取各连通域下对应时频谱的能量最大值，需将这些连通域分离出来。通过搜索二值化数组中各像素之间的连接关系，可分离出各连通域并进行顺序标注。

通过时频谱分割所获取的连通域即为各冲击信号时频能量最大值所在的区域。将各连通域映射到原始不同调节参数下的时频谱图中，即可得到各连通域在不同调节参数下的时频能量极大值序列，提取各序列下的最大值及其所对应的时频因子，即可精准定位各冲击信号出现的时频位置。根据提取的连通域下的一组最佳时频因子，结合故障过完备原子库的数学模型，即可构建一组用于稀疏表示故障信号的原子集合，最终得到故障信号的稀疏表示。

3.2 源个数的估计

传统的源个数的估计是依据数据时域散点图中的直线方向个数，常采用“K-均值聚类”方法[7-9]。该类方法对2个或3个源信号的源个数估计效果较好，而当源个数超过3个，则由于时域散点图的维数过多而难以凑效。

根据不同结构的设备其振动信号振荡衰减特性的差异性可知，不同振源产生的冲击信号的结构不同，特别是尺度因子和频率因子不同，而同一振源产生的冲击信号具有相似性[18]。因此，将所有稀疏表示提取的原子以尺度因子为横坐标、频率因子为纵坐标做出原子分布图，然后根据K-均值聚类方法对其进行分类，聚类的个数即为源个数。每个聚类的聚类中心所对应的原子可近似地替代该聚类下的其他原子，以简化后面的混合矩阵计算。

3.3 混合矩阵的估计

传统的基于稀疏分量的盲分离方法可以根据混合信号的散点图得到一系列代表列向量方向的直线，然而该方法坐标维度与观测信号数量一致，因此，仅适合欠定盲源分离，而对多路传感器的连轧机故障信号分析难以凑效。通过以下步骤可以更好地估算出混合矩阵的系数。

1) 更新稀疏表示。根据上一步的聚类分析，将每个聚类下的稀疏表示原子用该聚类的聚类中心所对应的原子替代，并重新计算各通道信号在新的原子集合下的稀疏表示。

2) 计算聚类系数。混合信号中的各通道信号均包含各振源信号，即各聚类信号，计算每个通道各聚类信号对应幅值的平均值，并按照同一顺序排列，得到一个由聚类系数组成的矩阵，该矩阵描述了各振源信号在各通道的能量衰减情况。

3) 估计混合矩阵系数。将上一步计算的矩阵每一列均进行归一化处理，即可估算出混合矩阵。

通过以上步骤估计的混合矩阵，首先，对振源和观测通道数量没有限制，可用于3个及以上观测通道的混合矩阵估计；其次，由于信号经过了稀疏特征提取，过滤了绝大部分噪声，从而使得混合矩阵的系数更加精确，克服了由散点图拟合列向量方向时的偏移问题。

4 仿真分析

仿真实验采用4个振源、2个观测通道的仿真数据。分别通过本研究方法和传统基于稀疏分量的盲源分离方法进行比较和分析。

根据式(8)模拟振源中的冲击成分

(8)

其中：λ为尺度因子；u为位移因子；f为频率因子。

4个振源冲击成分的时频因子如表1所示。

表1 冲击成分的时频因子

混合矩阵如式(9)所示

(9)

加噪后，2个观测通道的波形如图4所示，其中S1，S2分别表示混合后两个通道的信号。

图4 2个观测通道的波形

采用传统基于稀疏分量的盲源分离方法所得到的散点图及方向向量如图5所示。

图5 传统方法的散点图及方向向量

根据图5的方向向量可估算出混合矩阵为

(10)

图5和式(10)的结果表明，在信号掺杂部分噪声后，散点图的大多数点集中在噪声附件，根据K-均值聚类计算的方向向量具有一定的偏差。

采用本研究方法，首先通过基于时频谱分割的稀疏特征提取方法提取出各通道的稀疏信号，如图6所示。其中：SP1为S1稀疏表示信号；SR1为S1稀疏表示后的残差信号；SP2为S2稀疏表示信号；SR2为S2稀疏表示后的残差信号。

图6 本研究方法提取的稀疏信号

通过图6可以看出，2个通道信号中的冲击成分被准确提取出来了，而且各残差信号中基本没有冲击成分，稀疏表示的信号将噪声滤除了。

表2所示为各通道信号稀疏表示的原子时频因子信息。

表2 稀疏提取的原子时频因子

从表2可以看出，提取的时频因子信息与构成信号冲击成分的时频因子信息基本一致，由于噪声的影响，部分原子信息有微小的偏差。

图7所示为根据提取原子的结构信息产生的散点图及聚类中心。

图7 本研究方法的散点图及聚类中心

由图7可以看出，结构相似的原子主要聚集在4个点附近，根据4个聚类中心，即可完成稀疏表示的更新、聚类系数和混合矩阵的计算。计算后的混合矩阵为

(11)

比较式(9～11)不难发现，式(11)更加接近实际的混合矩阵，即抗干扰能力更强。

5 故障实例分析

通过旋转机械故障实验平台，在4个可更换轴承位置放置4路加速度传感器，并分别加上不同载荷。其中，3个轴承外圈有裂纹故障，1个轴承正常，分别采集4个测点的振动数据，并采用本方法进行分析。

图8所示为4个测点采集的振动数据，图9所示为本方法稀疏表示后的信号，图10所示为采用本方法分离出的各振源信号。

图8 4个测点采集的振动数据

图9 本研究方法稀疏表示后的信号

图10 本研究方法分离出的各振源信号

图8的4路振动信号S1，S2，S3和S4含有部分噪声。经过本研究方法的稀疏特征提取，可以看出图9中各自的稀疏表示信号SP1，SP2，SP3和SP4基本滤掉了噪声，保留了冲击成分。根据原子相似性原则进行混合矩阵的估算，最终得到图10所示的各独立振源信号C1，C2和C3。由于噪声的干扰，各源信号中冲击成分的幅值有微弱的波动。