CN110542557A - 图像积分驱动的机床大数据周期性故障特征快速分析方法 - Google Patents

Abstract

本发明涉及一种图像积分驱动的机床大数据周期性故障特征快速分析方法，对振动信号进行小波包分解实现多种耦合成分的分离，对每一个小波包重构序列的时域包络进行三次样条插值以提高采样密度，对插值后的图像进行二值化转换并进行形态学腐蚀及膨胀处理以消除图像上的非光滑点；选取图像的所有列和相连的部分列进行图像的逐步展开，根据取值非零像素点在空间上的关联性识别周期性特征产生的若干闭合连通区域，计算各区域的加权中心值并平均可以智能地识别周期性特征的出现频率，能够有效避免随机噪声和偶发性冲击的干扰，对振动信号的故障特征快速准确的识别。

Description

图像积分驱动的机床大数据周期性故障特征快速分析方法

技术领域

本发明涉机床故障诊断分析领域，具体涉及一种图像积分驱动的机床大数据周期性故障特征快速分析方法。

背景技术

高精度数控机床的动力由伺服电机作用于机械传动链并最终传递到主轴对工件材料进行去除，常期运行于非平稳工况下，传动链中的齿轮、轴承、转子等不可避免地发生疲劳从而产生局部损伤故障。为了在不停机拆解机床的前提下对传动链的健康状态进行评估，常采用振动测试及振动分析的方法，然而由旋转机械局部损伤故障所产生的微弱特征常淹没于其它成分的干扰中，因此必须通过信号处理方法挖掘其中的异常非平稳特征。

目前常通过时间尺度分析方法对信号进行分解，但分解后小波包信号的周期性冲击特征需要进行进一步分析，如“数控铣床电主轴故障诊断及仿真分析[J]，朱冬，王玲，殷国富，文华.工具技术,2019(1):125-128”；还有一种常见的方法是对小波包时域序列进行希尔伯特包络解调，在包络解调域上识别故障频率及其倍频成分，从而确认故障的发生，如“共振解调与小波降噪在电机故障诊断中的应用[J],张雄希，刘振兴，电机与控制学报,2010(6):66-70”。但由于机床振动测试所产生的数据量巨大，如果对每一条记录的数据都进行人工检验，将耗费大量的人力成本，为了减少分析对人工的依赖并提高特征分析的智能化程度，可以引入统计学中的峭度指标对小波包信号进行冲击性评估，然而这种统计过程中忽略了波形的几何特征，且容易将随机噪声干扰及偶发性冲击作为故障敏感特征进行错误识别。

发明内容

针对现有技术的缺陷，本发明的目的在于提供一种图像积分驱动的机床大数据周期性故障特征快速分析方法，模拟了周期性特征的人工肉眼识别方法，对图像从上往下逐步展开搜索，通过几何特征的机器视觉读取进行特征自动定位识别，可以在无人工干预的前提下对特征进行快速识别，识别准确且效率高，具有较好的工程应用推广价值。

为了达到上述目的，本发明所采用的技术方案是：

图像积分驱动的机床大数据周期性故障特征快速分析方法，包括以下步骤：

步骤1)在机床的非旋转部位安装振动加速度传感器采集振动信号，这些部位包括主轴的静止部分，伺服电机及变速齿轮箱的轴承座；振动信号序列 {x(n)|n＝1,2,…,N}的采样频率为f_s且采样长度为N，根据监测点附近机械传动链结构计算零部件(轴承、齿轮等)可能出现故障的特征频率，该故障特征频率的集合记为f_c＝{f_c(q)|且1≤q≤M}，其中表示整数集合，M为故障特征频率里总个数。

步骤2)采用Daubechies小波基函数对振动信号进行j层小波包分解并且进行单枝重构，产生2^j个小波包子空间：

{wp_j,k(n)|j＝1,2,…,J；k＝1,2,…,2^j；n＝1,2,…,N}，其中小波包子空间的顺序k按照理论频率通带的能量重心由小到大进行排列。

步骤3)对各重构小波包子空间{wp_j,k(n)}进行希尔伯特包络解调，获得时域上的包络序列Env{wp_j,k(n)}。

步骤4)采用三次样条插值对时域包络序列Env{wp_j,k(n)}进行插值，使插值包络序列Envc{wp_j,k(n)|n＝1,2,…,4N-3}的采样频率提升为4f_s，令 {wp_j,k(n)}的上界M_j,k表示为

其中算子max{·}表示取序列的最大值。

步骤5)对插值包络序列Envc{wp_j,k(n)}进行图像二值化，得到图像 I_j,k(u,v)的维度为(1000)×(4N-3),即行索引u(且1≤u≤1000)，列索引 v(且1≤v≤1000)。依次对I_j,k(u,v)中各列的像素进行初始化，初始化的方法为：对于图像中各像素按照如下表达式赋值:

步骤6)为了得到较为光滑的几何特征，对二值化图像进行5次形态学腐蚀及5次形态学膨胀处理，得到

步骤7)对步骤6)得到的图像进行9次分割，第i次分割选后的图像子集定义为：

即选取图像中第1000-100i+1到第1000行的所有像素。采用边界搜索算法对图像分割子集内灰度值为0的像素单联通区域(该区域的边界及内部所有像素的灰度值为0，而区域边界以外的相邻像素点灰度值为255) 进行识别，得到m_i个闭合区域,表示为当信号的数据点数较多时将造成很大的计算量，为此最多搜索10个单连通区域，对其它的区域停止搜索。计算从中识别的像素单通区域横坐标的加权重心：

计算小波包子空间{wp_j,k(n)}中单连通区域横坐标加权重心的平均间隔：

再判断其对应的频率是否与{f_c}中的某个特征频率f_c(q)接近。

这里的容许误差ε＝3Hz。对于图像分割子集如果某个满足上述条件，则将其原重构信号wp_j,k(n)及其对应的故障特征频率加入特征备选集合 E中。如果对于小波包子空间的所有内的单连通区域都不存在满足条件(1)，则说明wp_j,k(n)中不存在故障特征，继续搜索其它小波包子空间。

步骤9)若集合E非空，计算E中各小波包子空间的时域峭度值，选出峭度值最大的成分对应的小波包子空间的波形及特征频率进行观察，确认故障周期性冲击故障特征。如果特征备选集合E为空集，则说明信号中不存在故障特征。

进一步地，所述的区域边界搜索方法，具体包括以下步骤：

a)初始化图像子集的索引i＝1；

b)针对图像子集初始化单连通区域的索引初始化行方向的索引u＝1000，列方向的索引v＝0；

c)则针对图像子集中第个单连通区域，对输入的二值图像子集中的像素p(u,v)为起点进行从上往下，从左往右的搜索。当搜索到一个取值为0的像素点，表示为更新即在保持列索引v的值不变前提下对行索引的值逐次进行自减(u←u-1)搜索灰度值为0的像素，若该列不存在像素为0的点，则对列索引进行自加(v←v+1)。重复上述步骤直到发现灰度值为0的像素，然后进入步骤e)；若不存在这样的像素，则对图像子集索引i进行自加(i←i+1),并返回步骤b)。

d)在p(u,v)所在的行向下搜索所有灰度值连续为0的像素点，该连续灰度值为0像素集合的最下方像素记为满足下式的约束条件：

将上述所有满足要求的像素点加入集合

e)对于步骤d)中的取值为0集合的各像素逐一检索其右边相邻的像素。如果对于某个行索引满足则停止搜索并以该像素为起点在该行同时向上向下两个方向搜索灰度值连续取为0的其它像素点，更新及的数值，更新条件为：

如果上述更新条件无法满足，即该列像素中所有像素灰度值都为255，则对v进行自加(v←v+1)，对进行自加若j≤1000且则返回步骤c)；否则结束计算并将通过上述步骤得到图像子集中总数为m_i各灰度值为0的像素单连通区域表示为

有益效果：

1、本发明对振动信号进行小波包分解实现多种耦合成分的分离，对每一个小波包重构序列的时域包络进行三次样条插值以提高采样密度，对插值后的图像进行二值化转换并进行形态学腐蚀及膨胀处理以消除图像上的非光滑点；选取图像的所有列和相连的部分列进行图像的逐步展开，根据取值非零像素点在空间上的关联性识别周期性特征产生的若干闭合连通区域，计算各区域的加权中心值并平均可以智能地识别周期性特征的出现频率，如果与机床的故障特征频率相接近则可以确定周期性特征的存在性，再结合备选特征的时域峭度值可以找到最优特征，能够有效避免随机噪声和偶发性冲击的干扰，对振动信号的故障特征快速准确的识别。

2、本发明模拟了周期性冲击故障特征的人工肉眼识别方法，对图像从上往下逐步展开搜索，通过几何特征的机器视觉读取进行特征自动定位识别，本发明所提出的方法可以在无人工干预的前提下对特征进行快速识别，适用于机床传动系统的轴承、齿轮及蜗轮蜗杆等机械零部件的机械故障诊断，具有较好的工程应用推广价值。

附图说明

图1振动信号的小波包子空间包络曲线；

图2对细化后包络曲线进行二值转换得到的二维灰度图像；

图3从图像子集中搜索到的灰度值为0的所有像素单连通区域。

具体实施方式

下面结合附图和具体实施例对本发明做进一步详细的说明。

本发明在高速铣削工艺系统的非加工区域安装振动加速度传感器获取动态信号，信号采集过程中需要使用抗混叠滤波处理。

本发明图像积分驱动的机床大数据周期性故障特征快速分析方法具体包括以下步骤

步骤1).在机床的非旋转部位安装振动加速度传感器采集振动信号，这些部位包括主轴的静止部分，伺服电机及变速齿轮箱的轴承座；振动信号序列{x(n)|n＝1,2,…,N}的采样频率为f_s且采样长度为N，根据监测点附近机械传动链结构计算零部件(轴承、齿轮等)可能出现故障的特征频率，该故障特征频率的集合记为f_c＝{f_c(q)|且1≤q≤M}，其中表示整数集合，M为故障特征频率里总个数。

步骤2).采用Daubechies小波基函数对振动信号进行j层小波包分解并且进行单枝重构，产生2^j个小波包子空间：

{wp_j,k(n)|j＝1,2,…,J；k＝1,2,…,2^j；n＝1,2,…,N}，其中小波包子空间的顺序k按照理论频率通带的能量重心由小到大进行排列。

步骤3).对各重构小波包子空间{wp_j,k(n)}进行希尔伯特包络解调，获得时域上的包络序列Env{wp_j,k(n)}，如图1所示。

步骤4).采用三次样条插值对时域包络序列Env{wp_j,k(n)}进行插值，使插值包络序列Envc{wp_j,k(n)|n＝1,2,…,4N-3}的采样频率提升为4f_s，令 {wp_j,k(n)}的上界M_j,k表示为：

其中算子max{·}表示取序列的最大值。

步骤5).对插值包络序列Envc{wp_j,k(n)}进行图像二值化，得到图像 I_j,k(u,v)的维度为(1000)×(4N-3),即行索引u(且1≤u≤1000)，列索引 v(且1≤v≤1000)。依次对I_j,k(u,v)中各列的像素进行初始化，初始化的方法为：对于图像中各像素按照如下表达式赋值:

依照该步骤得到的二值图像如图2所示。

步骤6).为了得到较为光滑的几何特征，对二值化图像进行5次形态学腐蚀及5次形态学膨胀处理，得到

步骤7).对于输入的小波包子空间wp_j,k(n)，按照步骤6)的图像进行9 次分割，第i次分割选后的图像子集表示为：

即选取图像中第1000-100i+1到第1000行的所有像素。依照该步骤得到的二值图像子集如图3所示。采用边界搜索算法对图像分割子集内灰度值为0的像素单联通区域(该区域的边界及内部所有像素的灰度值为0，而区域边界以外的相邻像素点灰度值为255)进行识别，得到m_i个闭合区域,表示为由于故障特征体现为一系列冲击波包，在图像上体现为一系列单连通区域，识别了这些区域，就可以计算冲击的平均间隔。当信号的数据点数较多时,搜索所有的像素单连通区域计算量很大。因此当实际存在的像素单连通区域数量超过10个时，仅搜索前10个像素单连通区域。像素单连通区域的边界搜索算法步骤如下：

a)初始化图像分割子集的索引i＝1；

b)针对图像子集初始化单连通区域的索引初始化行方向的索引u＝1000，列方向的索引v＝0；

c)则针对图像子集中第个单连通区域，对输入的二值图像子集中的像素p(u,v)为起点进行从上往下，从左往右的搜索。当搜索到一个取值为0的像素点，表示为更新即在保持列索引v的值不变前提下对行索引的值逐次进行自减(u←u-1)搜索灰度值为0的像素，若该列不存在像素为0的点，则对列索引进行自加(v←v+1)。重复上述步骤直到发现灰度值为0的像素，然后进入步骤e)；若不存在这样的像素，则对图像子集索引i进行自加(i←i+1),并返回步骤b)。

d)在p(u,v)所在的行向下搜索所有灰度值连续为0的像素点，该连续灰度值为0像素集合的最下方像素记为满足下式的约束条件：

将所有满足上述要求的像素点加入集合

e)对于步骤d)中的取值为0集合的各像素逐一检索其右边相邻的像素。如果对于某个行索引满足则停止搜索并以该像素为起点在该行同时向上向下两个方向搜索灰度值连续取为0的其它像素点，更新及的数值，更新条件为：

如果上述更新条件无法满足，即该列像素中所有像素灰度值都为255，则对v进行自加(v←v+1)，对进行自加若j≤1000且则返回步骤c)；否则结束计算，并将通过上述步骤得到图像子集中总数为m_i的灰度值为0的像素单连通区域记为

然后计算所识别的各像素单连通区域的加权重心：

步骤8).计算中识别的像素单连通区域横坐标重心的平均间隔：

再计算时间间隔对应的频率判断是否与某个特征频率f_c(q)接近：

式(1)中容许误差ε＝3Hz。对于图像分割子集当某个满足上述判定条件，则将小波包子空间重构信号wp_j,k(n)及其对应的故障特征频率加入特征备选集合E中,并结束对图像子集索引的迭代检索。如果对于输入小波包子空间的所有内的单连通区域都不存在满足条件(1)，则说明该小波包子空间中不存在故障特征，继续搜索其它小波包子空间。

步骤9)若集合E非空，计算E中各小波包子空间的时域峭度值，选出峭度值最大的成分对应的小波包子空间的波形及特征频率进行观察，确认故障周期性冲击故障特征。如果特征备选集合E为空集，则说明信号中不存在故障特征。

下面结合一个工程实施案例对本发明的方法流程及附图进行进一步说明。

该工程实施案例中的振动信号采集自某故障模拟实验台，该试验台由直流电机驱动，负载由两个完全相同的质量块模拟。采用人工摩擦源模拟转子系统的动静碰摩故障，故障特征频率集合为f_c＝{36.67Hz}。

采用电涡流传感器获取竖直方向的振动位移信号x(n)，信号的采样频率 f_s＝2000Hz，信号的采样长度N＝1024；对信号x(n)进行j＝3层小波包分解，产生8个小波包子空间序列；搜索各小波包子空间中是否存在周期性冲击故障特征。对于某个wp_j,k(n)，其中j＝3，k＝1,2,…8。首先采样希尔伯特包络解调方法得到其包络曲线Env{wp_j,k(n)}，如图1所示；然后采用三次样条差值对 Env{wp_j,k(n)}进行采样细化，得到维度为长度为(4096-3)＝4093的包络解调曲线Envc{wp_j,k(n)}；按照Envc{wp_j,k(n)}中最大幅值的1.2倍(1.2M_j,k)对该曲线进行离散化，并生成二值图像，图像的维度为1000×4093。图像的横坐标为细化后的采样点索引；纵坐标对值域[0,1.2M_j.k]的区间进行1000等分，当p(u,v)的数值不大于Envc(u,v)时像素的灰度值取为0，否则取为255。wp_3,3(n)的二值化灰度图像如图2所示。取该二值图像的分割子集I_j,k,i并识别图像中灰度值为0 的像素单连通区域。图3为wp_3,3中i＝5时的图像分割子集，该图像子集中实际存在18个冲击特征对应的单连通区域；实际计算中为兼顾算法的有效性及效率，仅计算前10个单连通区域的横坐标能量重心平均间隔，结果为0.0273s，对应的频率为36.63Hz；该结果与系统的故障特征频率36.67Hz十分接近，因此将该小波包子空间加入集合E中；对8个小波包子空间按照上述方法分别进行特征提取，最终可以发现集合E非空，因此确认结构中存在机械故障。

本发明模拟了周期性特征的人工肉眼识别方法，对图像从上往下逐步展开，通过几何特征的机器视觉读取进行特征自动定位识别，所提出的方法可以在无人工干预的前提下对特征进行快速识别，具有较好的工程应用推广价值。

Claims (2)

1.图像积分驱动的机床大数据周期性故障特征快速分析方法，其特征在于，包括以下步骤：

步骤1)在机床的非旋转部位安装振动加速度传感器采集振动信号，振动信号序列{x(n)|n＝1,2,…,N}的采样频率为f_s且采样长度为N，根据监测点附近机械传动链结构计算零部件可能出现故障的特征频率，该故障特征频率的集合记为其中表示整数集合，M为故障特征频率里总个数；

步骤2)采用Daubechies小波基函数对振动信号进行j层小波包分解并且进行单枝重构，产生2^j个小波包子空间：{wp_j,k(n)|j＝1,2,…,J；k＝1,2,…,2^j；n＝1,2,…,N}，其中小波包子空间的顺序k按照理论频率通带的能量重心由小到大进行排列；

步骤3)对各重构小波包序列{wp_j,k(n)}进行希尔伯特包络解调，获得时域上的包络序列Env{wp_j,k(n)}；

步骤4)采用三次样条插值对时域包络序列Env{wp_j,k(n)}进行插值，使插值包络序列Envc{wp_j,k(n)|n＝1,2,…,4N-3}的采样频率提升为4f_s，令{wp_j,k(n)}的上界M_j,k表示为：

其中算子max{·}表示取序列的最大值；

步骤5)对插值包络序列Envc{wp_j,k(n)}进行图像二值化，得到图像I_j,k(u,v)的维度为(1000)×(4N-3)，即行索引列索引依次对I_j,k(u,v)中各列的像素进行初始化，初始化的方法为：对于图像中各像素按照如下表达式赋值：

步骤6)为了得到较为光滑的几何特征，对二值化图像进行5次形态学腐蚀及5次形态学膨胀处理，得到

步骤7)对步骤6)得到的图像进行9次分割，第i次分割选后的图像子集定义为

即选取图像中第1000-100i+1到第1000行的所有像素；采用边界搜索算法对图像分割子集内灰度值为0的像素单联通区域进行识别，单联通区域的边界及内部所有像素的灰度值为0，而区域边界以外的相邻像素点灰度值为255，得到m_i个闭合区域，表示为当信号的数据点数较多时将造成很大的计算量，为此最多搜索10个单连通区域，对其它的区域停止搜索，计算从中识别的像素单通区域横坐标的加权重心：

步骤8)计算小波包子空间{wp_j,k(n)}中单连通区域横坐标加权重心的平均间隔：

再判断其对应的频率是否与{f_c}中的某个特征频率f_c(q)接近；

此处容许误差ε＝3Hz，对于图像分割子集如果某个满足上述条件，则将其原重构信号wp_j,k(n)及其对应的故障特征频率加入特征备选集合E中；如果对于小波包子空间的所有内的单连通区域都不存在满足条件(1)，则说明wp_j,k(n)中不存在故障特征，继续搜索其它小波包子空间；

步骤9)若集合E非空，计算E中各小波包子空间的时域峭度值，选出峭度值最大的成分对应的小波包子空间的波形及特征频率进行观察，确认故障周期性冲击故障特征；如果特征备选集合E为空集，则说明信号中不存在故障特征。

2.如权利要求1所述的图像积分驱动的机床大数据周期性故障特征快速分析方法，其特征在于，所述的区域边界搜索方法，具体包括以下步骤：

a)初始化图像子集的索引i＝1；

b)针对图像子集初始化单连通区域的索引初始化行方向的索引u＝1000，列方向的索引v＝0；

c)则针对图像子集中第个单连通区域，对输入的二值图像子集中的像素p(u,v)为起点进行从上往下，从左往右的搜索；当搜索到一个取值为0的像素点，表示为更新即在保持列索引v的值不变前提下对行索引的值逐次进行自减u←u-1搜索灰度值为0的像素，若该列不存在像素为0的点，则对列索引进行自加v←v+1；重复上述步骤直到发现灰度值为0的像素，然后进入步骤e)；若不存在这样的像素，则对图像子集索引i进行自加i←i+1，并返回步骤b)；

d)在p(u,v)所在的行向下搜索所有灰度值连续为0的像素点，该连续灰度值为0像素集合的最下方像素记为满足下式的约束条件：

将上述所有满足要求的像素点加入集合

e)对于步骤d)中的取值为0集合的各像素逐一检索其右边相邻的像素，如果对于某个行索引满足则停止搜索并以该像素为起点在该行同时向上向下两个方向搜索灰度值连续取为0的其它像素点，更新及的数值，更新条件为：

如果上述更新条件无法满足，即该列像素中所有像素灰度值都为255，则对v进行自加v←v+1，对进行自加若j≤1000且则返回步骤c)；否则结束计算并将通过上述步骤得到图像子集中总数为m_i各灰度值为0的像素单连通区域表示为