一种基于POVMD和FDTW的时变转速下行星齿轮箱故障诊断方法与流程

一种基于povmd和fdtw的时变转速下行星齿轮箱故障诊断方法
技术领域
1.本发明属于旋转机械故障诊断与振动信号分析技术领域，涉及一种旋转机械故障诊断方法，具体涉及一种在时变转速下，无转速测量装置时，基于参数优化变分模态分解(parameter optimizationvariationalmode decomposition，povmd)和快速动态时间扭曲(fast dynamic time warping,fdtw)方法的行星齿轮箱故障诊断方法，可用于行星齿轮箱中齿轮和轴承的故障诊断。


背景技术：

2.行星齿轮箱具有体积小，传动比大、承载力强、运行平稳和工作效率高等优点，已被广泛应用在直升机、风力发电、重型卡车和舰船等大型复杂机械设备中。然而，行星齿轮箱作为旋转机械的主要传动部件，长期工作在冲击、交变载荷、工况大幅度变化的情况下，其关键零部件(太阳轮、行星轮、齿圈和滚动轴承等)，极易发生故障，从而引起巨大的经济损失和恶性事故。因此，对行星齿轮箱的故障诊断，及时发现故障信息，具有重大的经济与安全意义。
3.在变转速下，当行星齿轮箱内的齿轮或滚动轴承出现局部损伤时，从箱体上采集到的振动响应信号表现为非平稳、非线性和耦合调制等特征。传统的信号分析方法如快速傅里叶变换，分析此类信号，因为速度的波动，会出现频谱模糊的现象，不能准确得到故障信息。为了消除由于时变速度下频谱模糊的现象，可以采用阶次跟踪的方法。但是传统的阶次跟踪的方法需要额外的装置测量轴转速的信息，增加了成本，并且在实际情况下，额外的测量轴转速装置无法安装。为了克服此类问题，无转速测量的阶次跟踪的方法被提出来，基于时频分布的时频脊线提取来估计轴转速信息，但是受限于时频分辨率的影响，轴转速信息很难提取准确。因此，提出了自适应信号方法对于齿轮箱故障信号的轴转速信息的估计。经验模式分解(empirical mode decomposition,emd)方法作为较早使用的自适应信号分解方法，可以根据信号自身的局部尺度特征自适应地分解为一系列固有模态函数之和，从而揭示信号的内部本质。虽然emd方法具有很多优势，但存在模态混叠、端点效应以及理论基础缺乏等问题。为了弥补emd方法的不足，出现了集合经验模式分解(ensemble empirical mode decomposition,eemd)和噪声辅助分析的总体局部平均分解(ensemble local mean decomposition,elmd)等改进方法，但eemd和elmd中添加白噪声的两个关键参数(噪声的振幅和集成次数)选取机理尚未明确，同时添加噪声后在信号重构过程中存在残余噪声污染和运算量大等问题。不同于上述的分析方法，变分模态分解(variational mode decomposition，vmd)是一种新颖的自适应信号分解算法，用于非平稳信号的分解。因为具有完备的数学理论基础，而备受学者们的青睐。但是，vmd算法在应用时，模态个数k和惩罚因子α对vmd的分解结果有较大的影响，其中k数目过多，造成中心频率重复出现，过少则分解不完全；α值过大，分量的频带带宽变窄，过小则频带带宽变宽，影响模态分量识别，进而影响轴转速的估计。
4.此外，为了进一步提升无转速测量装置的阶次跟踪技术，需要从时域中考虑消除由转速变化带来的频谱模糊的现象。快速动态时间扭曲(fast dynamic time warping，fdtw),该算法通过扭曲由上述轴转速估计后的归一化轴振动信号和假定恒定转速下的轴振动信号的路径，由该路径扭曲原始振动信号，再重采样，消除由转速带来频谱模糊的现象。


技术实现要素：

5.鉴于上述行星齿轮箱时变转速下带来的频谱模糊的现象，导致故障信息无法准确得到问题和基于时频图像时频分辨率不高，导致转速信息无法准确得到的问题，本发明提供一种基于参数优化的变分模态分解和快速动态时间扭曲(povmd
‑
fdtw)的故障诊断方法，克服了上述的缺点。
6.为实现上述目的，本发明采用如下技术方案：
7.一种基于povmd
‑
fdtw的时变转速下行星齿轮箱故障诊断方法，其特征在于：包括如下步骤：
8.步骤s1：将加速度传感器安装在行星齿轮箱壳体的正上方采集故障振动信号x(t)(t为时间序列)；
9.步骤s2：优化对行星齿轮箱故障信号x(t)进行vmd分解的模态个数k和惩罚因子α；
10.步骤s3：基于步骤s2中的参数，对行星齿轮箱故障信号x(t)进行povmd分解；
11.步骤s4：计算分解后每个本征模态函数(imf)的瞬时频率，选择啮合频率f
m
附近的imf；
12.步骤s5：基于步骤s4所选择的imf，进行希尔伯特变换，并归一化所选择的imf为x
im
f,得到瞬时相位φ
inst
(t)信息，对时间t进行求导，再除以齿圈齿数z
r
和传动比i得到太阳轮绝对转频信息；
13.步骤s6：基于步骤s5中的太阳轮轴转速信息，建立时变转速下轴振动信号x
sh
(t)；
14.步骤s7：建立定转速下的轴参考信号y(t)。
15.步骤s8：对步骤s6中获取的时变转速下轴振动信号x
sh
(t)和步骤s7获取的定转速下的轴参考信号y(t)实施fdtw算法，实现它们之间的最佳对准，获得最佳扭曲路径w，能够消除由转速变化带来的两个信号之间的相位差；
16.步骤s9：基于步骤s8获取的最佳扭曲路径w对行星齿轮箱故障信号x(t)进行扭曲，得到扭曲信号x1(k)，并重采样，使得x1(k)恢复到原来信号的长度，得到重采后的信号z(n)；
17.步骤s10：对于步骤s9获得的重采样信号z(n)，进行快速傅里叶变换，得到阶次谱，判断故障信息；
18.进一步的，步骤s2中优化vmd分解的模态个数k和惩罚因子α过程为：首先确定最佳的模态个数，递减模态个数，使得在最窄带宽的情况下，可以分解出没有频率混叠的模态信号；为了确定每个模态信号的最佳带宽，首先递减一次模态个数，然后递减惩罚因子值，使模态间不发生带宽混叠的情况，在计算到带宽混叠的情况时回退一步惩罚因子的值，得到最佳的模态个数k和惩罚因子α。
19.进一步的，步骤s3中基于步骤s2中的参数，对行星齿轮箱故障信号x(t)进行povmd分解具体过程为：对于第k个模式，通过希尔伯特变换计算相关的分析信号，以获得单边频谱。通过指数将每个模式的频谱移动到基带，同时也将估计中心频率调整到模式所对应功率谱的重心位置。通过输入信号的高斯平滑，也就是梯度的二范数的平方，来度量带宽。因此，约束变分问题由下式给出：
[0020][0021]
式中，f(t)是给定的信号，u
k
是第k个模式，ω
k
是中心频率，k是模式的数量，是函数对时间t的偏导数，δ(t)是狄拉克分布，*是卷积运算符。
[0022]
引入二次惩罚项α和拉格朗日乘法算子λ，使得问题转换为非约束问题。因为二次惩罚项在有限权重下的良好收敛性和λ对约束的严格执行。从而给出了包含λ和二次惩罚项α的最优解的方程：
[0023][0024]
使用交替方向乘法器优化算法求解上式：
[0025][0026][0027]
进一步的，步骤s4计算分解后每个本征模态函数(imf)的瞬时频率的具体过程为：对每个imf进行希尔伯特变换，求出每个imf的瞬时相位信息，并对时间求导，得出每个imf的瞬时频率信息，找到啮合频率附近的imf,进而后续的分析。
[0028]
进一步的，步骤s5得到太阳轮轴转速的具体过程为：基于步骤s4所选择的imf归一化的信号为如下式：
[0029][0030][0031]
式中，是x
imf
(t)的希尔伯特变换。瞬时相位和太阳轮轴转速信息计算如下式：
[0032]
[0033][0034]
式中，z
r
为齿圈的齿数，i为行星齿轮箱的传动比。
[0035]
进一步的，步骤s6建立时变转速下轴振动信号的x
sh
(t)的具体过程为：
[0036][0037]
初始相位θ由0到2π的微小增量，找到和x
sh
(t)的最小均方差的相位角估计的。
[0038]
进一步的，步骤s7建立定转速下的轴振动信号y(t)的具体过程为：
[0039]
y(t)＝cos(2πf
s
t β)
[0040]
定转速下的速度f
s
由步骤s6中x
sh
(t)的峰值个数p来确定，f
s
的取值范围为到其中t
zong
为采样总时间，初始相位信息β由0到2π的微小增量，找到y(t)和x
sh
(t)的最小均方差的相位角估计的。
[0041]
进一步的，步骤s8采用fdtw算法实现的时变转速下轴振动信号x
sh
(t)和转速下的轴参考信号y(t)的最佳对准路径的具体过程为：x
sh
(t)和y(t)的长度分别为m和n:
[0042]
x
sh
(t)＝x1,x2,...,x
i
,...,x
m
[0043]
y(t)＝y1,y2,...,y
i
,...,y
n
[0044]
扭曲路径，w:
[0045]
w＝w1,w2,...,w
k
,...,w
k
[0046]
k是弯曲的长度，弯曲路径的第k个元素为w
k
＝(i,j)。其中i是来自x
sh
(t)的时间序列索引，j是来自y(t)的时间序列索引，如果弯曲路径为(i,j)，则代表时间序列x
sh
(t)的第i个元素和时间序列y(t)的第j个元素对齐。
[0047]
用欧式距离寻找一条最短的路径，获得x
sh
(t)和y(t)两个时间序列间的最短距离即计算公式如下式：
[0048][0049]
d(x
sh
(i),y(j))＝(x
sh
(i)
‑
y(j))2[0050]
式中，d表示两个点之间的距离，通过计算两个时间序列的不同点之间的距离得到一个m
×
n的距离矩阵：
[0051][0052]
选择一条最小的扭曲代价：
[0053][0054]
则迭代过程中的扭曲代价距离如下：
[0055][0056]
初始值为d(2,2)＝d(x
sh
(1),y(1)) d(x
sh
(2),y(2))
[0057]
则得到m
×
n大小的矩阵d,从矩阵d找出一条距离最短的一条扭曲路径为w。
[0058]
扭曲路径必须满足以下条件：
[0059]
(1)w1＝(1,1),扭曲路径必须在每个时间序列的开始处开始；
[0060]
(2)w
k
＝(n,m),扭曲路径必须在两个时间序列的末尾结束；
[0061]
而快速动态时间扭曲是在迭代过程的不断缩小迭代范围，减少迭代的复杂度。
[0062]
进一步的，步骤s9运用最佳扭曲路径w对行星齿轮箱故障信号x(n)进行扭曲，并重新采样使得扭曲信号恢复到原始长度的具体过程为：从步骤s8获得的最佳扭曲路径w获得路径w(i)，用此路径扭曲行星齿轮箱故障信号x(n)，得到扭曲信号x1(k)，重采样后的信号z(n)。重采样到原始信号长度的原理是，从扭曲路径w中的w(j)开始到扭曲路径w的长度k,如果w
k
‑1中的j”和j不相等，则z(n)＝x1(k)，否则取两者中的最大的那个。
[0063]
进一步的，步骤s10的具体过程为：对于重采样信号z(n)，进行快速傅里叶变化，得到阶次谱，从而判断故障信息。
[0064]
本发明与现有技术的优势在于具有以下特点：
[0065]
(1)本发明不仅发挥了povmd算法自适应参数优化分解信号得到单分量信号并进而估计轴转频的优势，而且还利用fdtw算法从时域中消除由转速波动带来的频谱模糊的现象，解决了传统方法等角度重采样的插值误差、额外需要测量转速装置和需要精确转速信息的问题。
[0066]
(2)本发明的povmd算法中参数优化克服了模态个数k数目过多，造成中心频率重复出现，过少则分解不完全的问题和惩罚因子α值过大，分量的频带带宽变窄，过小则频带带宽变宽，影响模态分量识别的问题。
[0067]
(3)本发明的fdtw算法与之前无转速测量装置的阶次跟踪技术想法不同，通过扭曲所测量的时变转速下的振动信号将原始的时间刻度投影到新的时间刻度上，其中轴的时变速度朝着恒定的参考速度挤压，因此在这个变换的时间尺度上，其频率随速度变化的模糊频谱也被挤压成各自的单独峰值。
[0068]
(4)本发明除了测量时变转速下的齿轮箱的振动信号外，所需要的轴转速信息是近似的，不需要像传统的阶次分析所需转速信息那么精确。
附图说明
[0069]
图1是本发明诊断方法流程图；
[0070]
图2是本发明实施例中时变转速下行星齿轮箱故障试验台示意图；
[0071]
图3是本发明实施列中太阳轮断齿故障振动信号的时域波形及其傅里叶频谱；
[0072]
图4是本发明实施列中太阳轮故障信号的povmd分解；
[0073]
图5是本发明实施列中太阳轮故障信号的povmd分解每个imf的瞬时频率；
[0074]
图6是本发明实施列中太阳轮轴转速估计的结果；
[0075]
图7是本发明实施列中太阳轮轴时变转速下振动信号的结果；
[0076]
图8是本发明实施列中参考信号示意图；
[0077]
图9是本发明实施列中fdtw算法实施后扭曲路径图；
[0078]
图10是本发明实施列中fdtw算法实施扭曲后并重采样的信号图；
[0079]
图11是本发明实施列中最后诊断的阶次谱。
[0080]
图中：1驱动电机；2变频器；3扭矩传感器；4轴承座；5采集卡；6行星齿轮箱；7张力控制器；8磁粉加载器。
具体实施方式
[0081]
以下结合附图和实施例对本发明的具体实施过程做进一步说明。
[0082]
本发明提供一种基于povmd
‑
fdtw的时变转速下行星齿轮箱故障诊断方法，流程图如图1所示，其主要步骤如下：
[0083]
步骤s1：将加速度传感器安装在行星齿轮箱壳体的正上方采集故障振动信号x(t)(t为时间序列)；
[0084]
步骤s2：优化对行星齿轮箱故障信号x(t)进行vmd分解的模态个数k和惩罚因子α；
[0085]
步骤s3：基于步骤s2中的参数，对行星齿轮箱故障信号x(t)进行povmd分解；
[0086]
步骤s4：计算分解后每个本征模态函数(imf)的瞬时频率，选择啮合频率f
m
附近的imf；
[0087]
步骤s5：基于步骤s4所选择的imf，进行希尔伯特变换，并归一化所选择的imf为得到瞬时相位φ
inst
(t)信息，对时间t进行求导，再除以齿圈齿数z
r
和传动比i得到太阳轮绝对转频信息；
[0088]
步骤s6：基于步骤s5中的太阳轮轴转速信息，建立时变转速下轴振动信号x
sh
(t)；
[0089]
步骤s7：建立定转速下的轴参考信号y(t)。
[0090]
步骤s8：对步骤s6中获取的时变转速下轴振动信号x
sh
(t)和步骤s7获取的定转速下的轴参考信号y(t)实施fdtw算法，实现它们之间的最佳对准，获得最佳扭曲路径w，能够消除由转速变化带来的两个信号之间的相位差；
[0091]
步骤s9：基于步骤s8获取的最佳扭曲路径w对行星齿轮箱故障信号x(n)进行扭曲，得到扭曲信号x1(k)，并重采样，使得x1(k)恢复到原来信号的长度，得到重采后的信号z(n)；
[0092]
步骤s10：对于步骤s9获得的重采样信号z(n)，进行快速傅里叶变换，得到阶次谱，判断故障信息；
[0093]
实施例1：
[0094]
搭建行星齿轮箱故障诊断实验台，如图2所示。该实验台主要是由1(驱动电机)、6(行星齿轮箱)、8(磁粉加载器)等构成。其中行星齿轮箱的结构参数见表1。通过线切割技术在太阳轮的某个轮齿上切除一部分作为故障。实验时将故障太阳轮安装在行星齿轮箱内，进行实验数据采集。加速度传感器安装在行星齿轮箱外壳的垂直、水平与轴向的测点上。电机输入转速5
‑
8hz线性增加，数据采样频率设置为12800hz，采样时长为2s。在该实验条件下，可以计算得到行星齿轮箱太阳轮局部故障相对应阶次如表2所示，其中m和k为正整数。
[0095]
表1行星齿轮箱中齿轮参数(单位/个)
[0096][0097]
表2太阳轮局部故障相应特征阶次
[0098][0099]
本发明基于povmd
‑
fdtw的时变转速下行星齿轮箱故障诊断方法的流程如图1所示，具体包括以下步骤：
[0100]
步骤一：利用加速度传感器采集时变转速下所述故障诊断实验台中行星齿轮箱太阳轮故障振动信号。图3为垂直方向加速度传感器获取到的太阳轮断齿故障信号的时域波形及其傅里叶频谱。由图3可知，振动信号的时域波形中难以观察到有规律的周期性冲击特征，频谱中太阳轮故障特征频率由于转速的变化，频谱出现模糊的现象，噪声与其他无关干扰频率的幅值较为突出。因此，通过传统的时域和频域分析方法无法从原始故障信号中提取表征太阳轮健康状态的特征信息。
[0101]
步骤二：将优化vmd的两个参数模态个数k和惩罚因子α，参数优化的具体步骤为:首先确定最佳的模态个数，递减模态个数，使得在最窄带宽的情况下，可以分解出没有频率混叠的模态信号；为了确定每个模态信号的最佳带宽，首先递减一次模态个数，然后递减惩罚因子值，使模态间不发生带宽混叠的情况，在计算到带宽混叠的情况时回退一步惩罚因子的值，得到最佳的模态个数k和惩罚因子α。本发明在这个实例中分解的模态个数k＝10，惩罚因子α＝3300。
[0102]
步骤三：基于步骤二的参数信息，将步骤一中采集到的故障振动信号进行povmd分解，如图4所示，具体步骤如下：给出了包含λ和二次惩罚项α的最优解的方程：
[0103][0104]
使用交替方向乘法器优化算法求解上式：
[0105][0106][0107]
步骤四：对步骤三每个imf进行希尔伯特变换，求出每个imf的瞬时相位信息，并对时间求导，得出每个imf的瞬时频率信息，如图5所示，找到啮合频率附近的imf,进而后续的分析。
[0108]
步骤五：基于步骤四所选择的imf得到太阳轮轴转速，如图6所示，其具体过程为：imf归一化的信号为如下式：
[0109][0110][0111]
式中，是x
imf
(t)的希尔伯特变换。瞬时相位和太阳轮轴转速信息计算如下式：
[0112][0113][0114]
式中，z
r
为齿圈的齿数，i为行星齿轮箱的传动比。
[0115]
步骤六：建立时变转速下轴振动信号的x
sh
(t)，如图7所示，其具体过程为：
[0116][0117]
初始相位θ由0到2π的微小增量，找到和x
sh
(t)的最小均方差的相位角估计的。
[0118]
步骤七:基于步骤六建立定转速下的轴振动信号y(t)，如图8所示，其具体过程为：y(t)＝cos(2πf
s
t β)
[0119]
定转速下的速度f
s
由步骤六中x
sh
(t)的峰值个数p来确定，f
s
的取值范围为到初始相位信息β由0到2π的微小增量，找到y(t)和x
sh
(t)的最小均方差的相位角估计的。
[0120]
步骤八:采用fdtw算法实现的时变转速下轴振动信号x
sh
(t)和转速下的轴参考信号y(t)的最佳扭曲路径w，如图9所示，其具体过程为：x
sh
(t)和y(t)的长度分别为m和n:
[0121]
x
sh
(t)＝x1,x2,...,x
i
,...,x
m
[0122]
y(t)＝y1,y2,...,y
i
,...,y
n
[0123]
扭曲路径，w:
[0124]
w＝w1,w2,...,w
k
,...,w
k
[0125]
k是弯曲的长度，弯曲路径的第k个元素为w
k
＝(i,j)。其中i是来自x
sh
(t)的时间序列索引，j是来自y(t)的时间序列索引，如果弯曲路径为(i,j)，则代表时间序列x
sh
(t)的第i个元素和时间序列的第j个元素对齐。
[0126]
用欧式距离寻找一条最短的路径，获得x
sh
(t)和y(t)两个时间序列间的最短距离即计算公式如下式：
[0127][0128]
d(x
sh
(i),y(j))＝(x
sh
(i)
‑
y(j))2[0129]
式中，d表示两个点之间的距离，通过计算两个时间序列的不同点之间的距离得到
一个m
×
n的距离矩阵：
[0130][0131]
选择一条最小的扭曲代价：
[0132][0133]
则迭代过程中的扭曲代价距离如下：
[0134][0135]
初始值为d(2,2)＝d(x
sh
(1),y(1)) d(x
sh
(2),y(2))
[0136]
则得到m
×
n大小的矩阵d,从矩阵d找出一条距离最短的一条扭曲路径为w。
[0137]
步骤九：基于步骤八获得的最佳扭曲路径w对行星齿轮箱故障信号x(n)进行扭曲，并重新采样使得扭曲信号恢复到原始长度，如图10所示，其具体过程为：从步骤八获得的最佳扭曲路径w获得路径w(i)，用此路径扭曲行星齿轮箱故障信号x(n)，得到扭曲信号x1(k)，重采样后的信号z(n)。重采样到原始信号长度的原理是，从扭曲路径w中的w(j)开始到扭曲路径w的长度k,如果w
k
‑1中的j”和j不相等，则z(n)＝x1(k)，否则取两者中的最大的那个。
[0138]
对步骤九中获得的重采样信号z(n)进行快速傅里叶变换分析，得到阶次谱，如图11所示。由图11可知，在阶次谱有很明显的太阳轮的绝对转频阶次o,啮合阶次o
m
，边频阶次o
m
±
mo
s
±
ko,m和k为正整数(m＝1,2,3；k＝1,2)。因此上述分析表明了太阳轮出现了故障，这与实验设置相一致。从而验证了本发明的方法的有效性和可行性。
[0139]
综上所述，本发明基于povmd
‑
fdtw的时变转速下行星齿轮箱故障诊断方法，能够成功地提取微弱的故障特征信息，实现行星齿轮箱早期故障的准确识别与诊断。fdtw通过扭曲所测量的时变转速下的振动信号将原始的时间刻度投影到新的时间刻度上，其中轴的时变速度朝着恒定的参考速度挤压，因此在这个变换的时间尺度上，其频率随速度变化的模糊频谱也被挤压成各自的单独峰值。测量时变转速下的齿轮箱的振动信号外，所需要的轴转速信息是近似的，不需要像传统的阶次分析所需转速信息那么精确。
[0140]
以上所述仅为本发明的较好实施例，凡依本发明申请专利范围所做的均等变化修饰，皆应属本发明的涵盖范围。