CN112284727B - 一种基于卷积极大极小凹罚算法的旋转机械故障诊断方法 - Google Patents

Abstract

本发明公开一种基于卷积极大极小凹罚算法的旋转机械故障诊断方法，包括步骤：步骤S1、采集旋转机械振动加速度响应信号及转速信号，并记录可能存在的故障特征信息；步骤S2、从采集到的数据中截取一段振动信号，以移不变K‑SVD方法从中获取模式；步骤S3、以获取到的模式作为输入，通过卷积极大极小凹罚算法求解稀疏系数，将稀疏系数与模式卷积，得到重构故障特征；步骤S4、分析提取的故障特征信号的时域特征与解调谱特征，确定故障类型，完成故障诊断。本发明提出以卷积极大极小凹罚作为罚函数，不仅改善了现有方法存在的幅值低估问题，而且将算法求解过程从时域转换到频域中进行，极大提高了运算效率。

Description

一种基于卷积极大极小凹罚算法的旋转机械故障诊断方法

技术领域

本发明属于旋转机械故障诊断领域，更具体地，涉及一种基于卷积极大极小凹罚算法的旋转机械故障诊断方法。

背景技术

旋转机械服役在复杂多变的环境中，若出现故障将导致严重的后果，因此对旋转机械进行状态监测就显得至关重要。旋转机械的振动响应信号中包含大量的噪声，如何从含噪振动信号中准确而快速地提取出故障特征是机械故障诊断中的一个难题。

谱峭度(Kurtogram)、经验模态分解(Empirical Mode Decomposition,EMD)及卷积稀疏表示(Convolutional Sparse Coding,CSC)等信号处理方法被应用于机械故障诊断中，但仍然方法存在着不足。如谱峭度方法可能导致共振带被错误分割，降低共振带内故障特征的信噪比(丁康,黄志东,林慧斌.一种谱峭度和Morlet小波的滚动轴承微弱故障诊断方法[J].振动工程学报,2014,27(01):128-135.)。经验模态分解方法易造成模态混叠，进而影响特征提取的精度(胡爱军,孙敬敬,向玲.经验模态分解中的模态混叠问题[J].振动.测试与诊断,2011,31(04):429-434+532-533.)。旋转机械局部型故障响应信号表现出卷积特性，可表示为一冲击模式与稀疏系数的卷积，现有的卷积稀疏表示方法能有效地提取出冲击特征(施莹,庄哲,林建辉.基于卷积稀疏表示及等距映射的轴承故障诊断[J].振动.测试与诊断,2019,39(05):1081-1088+1138.)，但此类方法普遍采用l₁正则项作为目标函数的正则项，会系统性地低估冲击特征的幅值，造成无法恢复的误差。

发明内容

本发明的目的是针对现有卷积稀疏表示方法会低估冲击特征幅值的缺陷，提供一种基于卷积极大极小凹罚算法的旋转机械局部型故障诊断方法。该方法以卷积极大极小凹罚作为罚函数，改善了幅值低估问题，此外还将稀疏系数求解过程从时域转换到频域进行，极大提高了运算效率，能准确而高效地提取出故障特征。

本发明至少通过如下技术方案之一实现。

一种基于卷积极大极小凹罚算法的旋转机械故障诊断方法，包括以下步骤：

步骤S1、采集旋转机械振动加速度响应信号及转速信号，并记录可能存在的故障特征信息；

步骤S2、从采集到的数据中截取一段振动信号，以移不变K-SVD方法从中获取模式，模式为一波形呈现冲击衰减特性的向量；

步骤S3、以获取到的模式作为输入，通过卷积极大极小凹罚算法求解稀疏系数，将稀疏系数与模式卷积，得到重构故障特征；

步骤S4、分析提取的故障特征信号的时域特征与解调谱特征，确定故障类型，完成故障诊断。

优选的，所述步骤S1中具体包括：

S11、将三向加速度传感器安装于设备旋转轴承座处，并连接好数据采集仪与便携式计算机；

S12、设置数据采集仪的采样频率为f_s，采集并保存对应的振动加速度时域信号，以光电脉冲编码器采集输出轴的转速信号，根据设备结构分析并记录可能出现的齿轮以及轴承的故障类型，记录由运行工况决定的故障特征频率。

优选的，所述步骤S2具体包括：

S21、从采集到的振动加速度信号中截取一段信号x∈R^N×1，R^N×1表示N维实数列向量；

S22、设置模式的长度为l个采样点，以移不变K-SVD方法从信号x中获取模式d∈R^{l ×1}，包含故障信息的模式应呈现冲击衰减波形；

移不变K-SVD(Shift Invariant K-SVD,SI K-SVD)方法在获取模式时可直接对一维信号作处理而无需分段，其求解过程是一个交替更新的过程，包括稀疏系数更新和模式更新两个阶段，在稀疏系数更新阶段以匹配追踪(Matching Pursuit,MP)算法进行求解，其目标函数为：

式中y是使目标函数取最小值的稀疏系数，||y||₀是稀疏系数的稀疏度，s为稀疏系数y的稀疏度阈值，τ为稀疏系数中非0时刻，y_τ是时刻τ对应的稀疏系数值，T_τ∈R^N×l是用于将模式d时移到τ时刻处的时移算子；

在模式更新阶段，模式将根据其支撑集相继进行更新，支撑集表示为σ＝{τ|y_τ≠0}，σ是所有非0稀疏系数中对应时刻的集合，令残量

则通过奇异值分解(Singular Value Decomposition,SVD)方法来更新模式d，目标函数为

其中算子adj(T_τ)是时移算子T_τ的伴随矩阵，能够从残量r中提取位于时刻τ与模式相同长度的信号片段，移不变K-SVD方法以匹配追踪与奇异值分解在迭代过程中交替更新稀疏系数与模式，当停止迭代时，输出模式d。

优选的，所述步骤S3中具体包括：

S31、判断模式d的波形是否呈现冲击衰减特性，若是，则执行后续步骤，若否，则表明设备尚未出现局部故障；

S32、将模式d作为输入，设正则化参数为λ，通过卷积极大极小凹罚算法(Convolutional Generalized Minimax Concave,CGMC)求解如下目标函数：

其中，y为稀疏系数，z是辅助变量，

为向量l₂范数的平方，||₁||为向量的l₁范数，γ是用于控制目标函数凸性的参数，

是使是使卷积极大极小凹罚函数f(y,z)取最大值的辅助变量，

是使函数

取最小值的稀疏系数，求解上式所得稀疏系数

表示冲击发生的时刻与幅度；

S33、将模式d与稀疏系数

卷积，得到重构故障特征信号

其中*为卷积运算符。

优选的，求解上式所得稀疏系数

具体步骤如下：

(1)输入信号x、模式d、正则化参数λ、常量μ、参数γ以及最大迭代次数I_max；

(2)初始化辅助变量s⁽⁰⁾和v⁽⁰⁾为零向量，令迭代次数i＝0，模式d补零为并归一化为d₁∈R^N×1；

(3)频域更新辅助变量w⁽ⁱ⁾

式中

表示向量的傅里叶变换，

为向量傅里叶变换的共轭，均可由FFT(FastFourier Transform)快速计算，

表示向量的傅里叶逆变换，可通过IFFT(InverseFast Fourier Transform)快速计算，频域更新w⁽ⁱ⁾的计算量以FFT和IFFT的计算量为主，所需的时间复杂度及空间复杂度都较低；

(4)频域更新辅助变量u⁽ⁱ⁾

(5)更新辅助变量s⁽ⁱ⁺¹⁾

s⁽ⁱ⁺¹⁾＝T_λμ(w⁽ⁱ⁾)

其中T_λμ为软阈值函数，λμ为其阈值，

为向量点乘运算，绝对值函数|·|、符号函数sign(·)和max(·)函数均是对向量的每一个分量做处理；

(6)更新辅助变量v⁽ⁱ⁺¹⁾

v⁽ⁱ⁺¹⁾＝T_λμ(u⁽ⁱ⁾)

(7)令i＝i+1，若i>I_max，则进行后续步骤(8)，若i≤I_max，则继续重复步骤(3)-步骤(7)；

(8)输出稀疏系数

上述求解步骤中涉及到的正则化参数λ设置为λ＝αδ，其中α取3.3～10，δ为信号x中噪声的估计标准差，由δ＝MAD(x)/0.6745估计得到，其中MAD(x)为信号x的绝对中位差，MAD(x)＝median(|x-median(x)|)，median(·)为向量的中位数，常量

其中常量

关键参数γ取0.8，最大迭代次数I_max为1000次。

优选的，所述步骤S4中，通过分析提取的故障特征信号的时域特征与解调谱特征，即可确定故障类型，完成故障诊断。

本发明与现有技术相比，具有如下优点和有益效果：

(1)该方法以卷积极大极小凹罚作为惩罚项，能有效改善现有卷积稀疏表示方法存在的幅值低估问题；

(2)该方法将卷积极大极小凹罚算法的求解过程由时域转换到频域进行，所需的时间复杂度较低，能极大提高运算效率；

(3)该方法适用性广，旋转机械局部型故障的振动响应信号都具有移不变特性，故该方法适用于对所有旋转机械的局部型故障进行诊断。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍。附图构成本申请的一部分，但仅是作为体现发明概念的非限制性示例，并非用于做出任何限制。

图1是本发明基于卷积极大极小凹罚算法的旋转机械故障诊断方法的流程图；

图2是本实施例中所使用的滚动轴承信号；

图3是本实施例中得到的冲击模式；

图4是本实施例中由卷积极大极小凹罚算法重构得到的冲击特征信号；

图5是本实施例中冲击特征信号的解调谱。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本实施例提供的一种基于卷积极大极小凹罚算法的旋转机械故障诊断方法，如图1所示，以旋转机械中发生局部故障的滚动轴承为研究对象，对本发明进行进一步说明。实验所采用的滚动轴承型号为N205M，其节径为38mm、滚动体直径为6.5mm、滚动体个数为13个，接触角为0。本实施例包括以下步骤：

S1、采集旋转机械振动加速度响应信号及转速信号，并记录可能存在的故障特征信息，如齿轮的局部故障特征频率以及轴承的内圈、外圈与滚动体故障的特征频率，具体包括如下步骤：

S11、将三向加速度传感器安装于设备旋转轴承座处，并连接好数据采集仪与便携式计算机；

S12、设置数据采集仪的采样频率为f_s＝100kHz，采集并保存对应的振动加速度时域信号，以光电脉冲编码器采集输出轴的转速信号，输出轴转速为500rpm，根据设备结构分析并记录可能出现的轴承外圈故障，记录由运行工况决定的外圈故障特征频率为f_o＝44.9Hz。

S2、从采集到的数据中截取一段振动信号，以移不变K-SVD方法从中获取模式，具体步骤如下：

S21、如图2所示，从采集到的振动信号中截取时长为1s的振动信号x，其点数为N＝100000点；

S22、设置模式d的长度为l＝550个点，以移不变K-SVD方法从信号x中获取的模式d∈R^l×1，包含故障信息的模式应呈现冲击衰减波形，R^l×1表示l维实数列向量；

移不变K-SVD(Shift Invariant K-SVD,SI K-SVD)方法在获取模式时可直接对一维信号作处理而无需分段，其求解过程是一个交替更新的过程，包括稀疏系数更新和模式更新两个阶段，在稀疏系数更新阶段以匹配追踪(Matching Pursuit,MP)算法进行求解，其目标函数为：

式中y是使目标函数取最小值的稀疏系数，||y||₀是稀疏系数的稀疏度，s为稀疏系数y的稀疏度阈值，τ为稀疏系数中非0时刻，y_τ是时刻τ对应的稀疏系数值，T_τ∈R^N×l是用于将模式d时移到τ时刻处的时移算子；

在模式更新阶段，模式将根据其支撑集相继进行更新，支撑集表示为σ＝{τ|y_τ≠0}，σ是所有非0稀疏系数中对应时刻的集合，令残量

则通过奇异值分解(Singular Value Decomposition,SVD)方法来更新模式d，目标函数为

其中算子adj(T_τ)是时移算子T_τ的伴随矩阵，能够从残量r中提取位于时刻τ与模式相同长度的信号片段，移不变K-SVD方法以匹配追踪与奇异值分解在迭代过程中交替更新稀疏系数与模式，当迭代次数达到100时停止迭代，输出模式d。

S3、以获取到的模式作为输入，通过卷积极大极小凹罚算法求解稀疏系数，将稀疏系数与模式卷积，得到重构故障特征，使用的卷积极大极小凹罚算法兼具稀疏系数求解的准确性与计算的高效性，具体包括如下步骤：

S31、图3所示模式d的波形很好地反映了冲击衰减特性，表明其中包含故障信息，因此可执行后续步骤；

S32、将模式d作为输入，设置正则化参数λ＝4.7δ，其中噪声标准差的估计值δ＝0.2698，以卷积极大极小凹罚算法求解得到稀疏系数

卷积极大极小凹罚算法求解如下目标函数：

其中，y为稀疏系数，z是辅助变量，

为向量l₂范数的平方，||·||₁为向量的l₁范数，γ是用于控制目标函数凸性的参数，

是使是使卷积极大极小凹罚函数f(y,z)取最大值的辅助变量，

是使函数

取最小值的稀疏系数，求解上式所得稀疏系数

表示冲击发生的时刻与幅度，

由冲击特征的移不变特性，上式中模式d与稀疏系数y的卷积d*y也可表示为移不变字典D∈R^N×N与稀疏系数y的乘积，即d*y＝Dy，将模式d补零为并归一化为原子d₁∈R^N×1，字典D的第一列为d₁，D是由d₁时移而成的循环矩阵，根据循环矩阵的性质，等式右侧Dy通过下式转换到频域进行快速计算

式中

表示向量的傅里叶变换，

为向量傅里叶变换的共轭，均可由FFT(FastFourier Transform)快速计算，

向量的傅里叶逆变换，可通过IFFT(Inverse FastFourier Transform)快速计算，上式的计算由时域转换到频域后，计算量从等式左侧的O(N²)，降低到右侧的O(3Nlog₂N)，O(.)表明计算量与括号中的值处于同一数量级，因此计算效率极大提高。D^T∈R^N×N也是一循环矩阵，其第一列为

D^TD∈R^N×N同样为循环矩阵，且其第一列为

综合以上分析可知，根据循环矩阵的性质，移不变字典D、移不变字典的转置D^T以及矩阵D^TD与任意N维列向量的乘积均可以转换到频域中进行快速计算，计算量均以O(3Nlog₂N)为主。

卷积极大极小凹罚算法的目标函数具体求解步骤如下：具体步骤如下：

(1)输入信号x、模式d、正则化参数λ、常量μ、参数γ以及最大迭代次数I_max；

(2)初始化辅助变量s⁽⁰⁾和v⁽⁰⁾为零向量，令迭代次数i＝0，模式d补零为并归一化为d₁∈R^N×1；

(3)频域更新辅助变量w⁽ⁱ⁾

式中

表示向量的傅里叶变换，

为向量傅里叶变换的共轭，均可由FFT(FastFourier Transform)快速计算，

表示向量的傅里叶逆变换，可通过IFFT(InverseFast Fourier Transform)快速计算，频域更新w⁽ⁱ⁾的计算量以FFT和IFFT的计算量为主，所需的时间复杂度及空间复杂度都较低；

(4)频域更新辅助变量u⁽ⁱ⁾

(5)更新辅助变量s⁽ⁱ⁺¹⁾

s⁽ⁱ⁺¹⁾＝T_λμ(w⁽ⁱ⁾)

其中T_λμ为软阈值函数，λμ为其阈值，

为向量点乘运算，绝对值函数|·|、符号函数sign(·)和max(·)函数均是对向量的每一个分量做处理；

(6)更新辅助变量v⁽ⁱ⁺¹⁾

v⁽ⁱ⁺¹⁾＝T_λμ(u⁽ⁱ⁾)

(7)令i＝i+1，若i>I_max，则进行后续步骤(8)，若i≤I_max，则继续重复步骤(3)-步骤(7)；

(8)输出稀疏系数

上述求解步骤中涉及到的正则化参数λ设置为λ＝αδ，其中α取3.3～10，δ为信号x中噪声的估计标准差，由δ＝MAD(x)/0.6745估计得到，其中MAD(x)为信号x的绝对中位差，MAD(x)＝median(|x-median(x)|)，median(·)为向量的中位数，常量

其中常量

关键参数γ取0.8，最大迭代次数I_max为1000次。

S33、将所得模式d与所得稀疏系数

卷积，得到如图4所示的重构故障特征信号

S4、分析提取的故障特征信号的时域特征与解调谱特征，确定故障类型，完成故障诊断，具体包括如下步骤：

S4、分析图4可知，所得重构故障特征信号

时域中可以清晰地分辨出冲击发生时刻及间隔，其解调谱如图5所示，解调谱中外圈故障特征频率f_o的前6倍频非常突出，由时域和解调谱特征可断定该滚动轴承外圈存在局部型故障。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (5)

1.基于卷积极大极小凹罚算法的旋转机械故障诊断方法，其特征在于，包括以下步骤：

步骤S1、采集旋转机械振动加速度响应信号及转速信号，并记录可能存在的故障特征信息；

步骤S2、从采集到的数据中截取一段振动信号，以移不变K-SVD方法从中获取模式，模式为一波形呈现冲击衰减特性的向量；

步骤S3、以获取到的模式作为输入，通过卷积极大极小凹罚算法求解稀疏系数，将稀疏系数与模式卷积，得到重构故障特征；求解稀疏系数

具体步骤如下：

(1)输入信号x、模式d、正则化参数λ、常量μ、参数γ以及最大迭代次数I_max；

(2)初始化辅助变量s⁽⁰⁾和v⁽⁰⁾为零向量，令迭代次数i＝0，模式d补零为并归一化为d₁∈R^N×1；

(3)频域更新辅助变量w⁽ⁱ⁾

式中

表示向量的傅里叶变换，

为向量傅里叶变换的共轭，均可由FFT(Fast FourierTransform)快速计算，

表示向量的傅里叶逆变换，可通过IFFT(Inverse FastFourier Transform)快速计算，频域更新w⁽ⁱ⁾的计算量以FFT和IFFT的计算量为主，所需的时间复杂度及空间复杂度都较低；

(4)频域更新辅助变量u⁽ⁱ⁾

(5)更新辅助变量s⁽ⁱ⁺¹⁾

s⁽ⁱ⁺¹⁾＝T_λμ(w⁽ⁱ⁾)

其中T_λμ为软阈值函数，λμ为其阈值，

为向量点乘运算，绝对值函数|·|、符号函数sign(·)和max(·)函数均是对向量的每一个分量做处理；

(6)更新辅助变量v⁽ⁱ⁺¹⁾

v⁽ⁱ⁺¹⁾＝T_λμ(u⁽ⁱ⁾)

(7)令i＝i+1，若i>I_max，则进行后续步骤(8)，若i≤I_max，则继续重复步骤(3)-步骤(7)；

(8)输出稀疏系数

上述求解步骤中涉及到的正则化参数λ设置为λ＝αδ，其中α取3.3～10，δ为信号x中噪声的估计标准差，由δ＝MAD(x)/0.6745估计得到，MAD(x)为信号x的绝对中位差，MAD(x)＝median(|x-median(x)|)，median(·)为向量的中位数，常量

其中常量

关键参数γ取0.8，最大迭代次数I_max为1000次；

S33、将模式d与稀疏系数

卷积，得到重构故障特征信号

其中*为卷积运算符；

步骤S4、分析提取的故障特征信号的时域特征与解调谱特征，确定故障类型，完成故障诊断。

2.根据权利要求1所述的基于卷积极大极小凹罚算法的旋转机械故障诊断方法，其特征在于，所述步骤S1中具体包括：

S11、将三向加速度传感器安装于设备旋转轴承座处，并连接好数据采集仪与便携式计算机；

S12、设置数据采集仪的采样频率为f_s，采集并保存对应的振动加速度时域信号，以光电脉冲编码器采集输出轴的转速信号，根据设备结构分析并记录可能出现的齿轮以及轴承的故障类型，记录由运行工况决定的故障特征频率。

3.根据权利要求2所述的基于卷积极大极小凹罚算法的旋转机械故障诊断方法，其特征在于，所述步骤S2具体包括：

S21、从采集到的振动加速度信号中截取一段信号x∈R^N×1，R^N×1表示N维实数列向量；

S22、设置模式的长度为l个采样点，以移不变K-SVD方法从信号x中获取模式d∈R^l×1，包含故障信息的模式应呈现冲击衰减波形；

移不变K-SVD方法在获取模式时可直接对一维信号作处理而无需分段，其求解过程是一个交替更新的过程，包括稀疏系数更新和模式更新两个阶段，在稀疏系数更新阶段以匹配追踪算法进行求解，其目标函数为：

式中y是使目标函数取最小值的稀疏系数，||y||₀是稀疏系数的稀疏度，s为稀疏系数y的稀疏度阈值，τ为稀疏系数中非0时刻，y_τ是时刻τ对应的稀疏系数值，T_τ∈R^N×l是用于将模式d时移到τ时刻处的时移算子；

在模式更新阶段，模式将根据其支撑集相继进行更新，支撑集表示为σ＝{τ|y_τ≠0}，σ是所有非0稀疏系数中对应时刻的集合，令残量

则通过奇异值分解方法来更新模式d，目标函数为

其中算子adj(T_τ)是时移算子T_τ的伴随矩阵，能够从残量r中提取位于时刻τ与模式相同长度的信号片段，移不变K-SVD方法以匹配追踪与奇异值分解在迭代过程中交替更新稀疏系数与模式，当停止迭代时，输出模式d。

4.根据权利要求3所述的基于卷积极大极小凹罚算法的旋转机械故障诊断方法，其特征在于，所述步骤S3中具体包括：

S31、判断模式d的波形是否呈现冲击衰减特性，若是，则执行后续步骤，若否，则表明设备尚未出现局部故障；

S32、将模式d作为输入，设正则化参数为λ，通过卷积极大极小凹罚算法求解如下目标函数：

其中，y为稀疏系数，z是辅助变量，

为向量l₂范数的平方，||·||₁为向量的l₁范数，γ是用于控制目标函数凸性的参数，

是使卷积极大极小凹罚函数f(y,z)取最大值的辅助变量，

是使函数

取最小值的稀疏系数，求解上式所得稀疏系数

表示冲击发生的时刻与幅度。

5.根据权利要求4所述的基于卷积极大极小凹罚算法的旋转机械故障诊断方法，其特征在于，所述步骤S4中，通过分析提取的故障特征信号的时域特征与解调谱特征，即可确定故障类型，完成故障诊断。

Images