一种参数自适应分解结构的轴承故障特征增强方法

1.本发明涉及一种参数自适应分解结构的轴承故障特征增强方法，属于故障诊断技术与信号处理分析技术领域。


背景技术：

2.轴承作为旋转机械的支撑部件，其健康程度直接影响旋转机械的运行精度和寿命。轴承发生故障时，滚动体与滚道在故障位置处的接触刚度发生变化，对应的瞬时角速度(ias)将发生规律性变化，则ias信号包含丰富的轴承故障信息。因此，基于ias信号的轴承状态监测和故障诊断是故障诊断领域的热点之一。然而，编码器安装误差的抑制对滚动轴承故障特征提取是不可忽视的，其干扰分量的抑制对于轴承故障特征提取尤为重要。
3.一方面，由于编码器内径与轴装配误差，编码器安装误差在工程应用中无法避免；编码器安装误差的能量幅值与转速成正相关，即编码器安装误差对轴承信号的调制作用随着转速上升而显著提升；不同程度的偏心误差和倾斜误差产生的总体误差是不尽相同，即对轴承故障分量的干扰程度是不同的。另一方面，早期轴承引起的ias变化相对较弱，常常淹没于测量噪声和编码器安装误差中，导致其故障辨识较为困难。
4.综上，编码器安装误差的消除对于有效揭示轴承故障特征至关重要。


技术实现要素：

5.由于编码器安装误差在工程应用中无法避免，以及编码器安装误差对滚动轴承分量的调制作用随着转速的上升而增大，进而导致轴承特征可能无法有效辨识。为解决该问题，本发明提出一种参数自适应分解结构(pds)的轴承故障特征增强方法，该方法在savitzky-golay滤波器的基础上，基于参数自适应分解结构，结合iidf指标，自适应获得优化的savitzky-golay滤波器的优化参数，以有效消除编码器安装误差，进而增强轴承故障特征。
6.如图1所示，本发明基于参数自适应分解结构的轴承故障特征增强方法如下：
7.步骤1：获得包含轴承故障信息的瞬时角速度信号；
8.通过picoscope信号采集系统获取光学编码器的包含轴承故障信息的瞬时角位移和对应的时间，采用向前差分法计算获得瞬时角速度iasi信号，其计算式如下
[0009][0010]
式中iasi表示第i时刻的瞬时角速度，i＝1,2,3,
…
，δφ＝2π/n，δti＝t
i+1-ti；n表示编码器的光栅格数。
[0011]
步骤2：通过pds获得优化滤波长度m
op
[0012]
2-1、假设参数的取值范围为[a
l
→
p
,b
l
→
p
]，分解等级l＝1,2,3,...，p表示该等级内参数范围的个数，p＝1,2,3,...。特殊的，如果l＝1，则p＝1。滤波长度m
l
(d
l
)＝a
l
→
p
,a
l
→
p
+s
l
,...,floor((b
l
→
p-a
l
→
p
)/s
l
)，等级l范围中的离散序列个数d
l
＝1,2,...,∑
p floor
((b
l
→
p-a
l
→
p
)/s
l
)，floor(
·
)表示向下取整操作，m
l
(d
l
)的增量为s
l
，s
l
与s
l+1
比值是ε，ε》1；
[0013]
2-2、基于固定滤波阶次p和m
l
(d
l
)＝a
l
→
p
,a
l
→
p
+s
l
,...,floor((b
l
→
p-a
l
→
p
)/s
l
)，离散序列为d
l
，采用savitzky-golay滤波器对原始信号iasi进行滤波，计算式为
[0014][0015]
式中表示滤波后信号，n＝(m
l
(d
l
)-1)/2+1,(m
l
(d
l
)-1)/2+2,(m
l
(d
l
)-1)/2+3,...,length(iasi)，length(
·
)表示取数据长度操作，滤波系数可表示为
[0016][0017]
式中d
t
＝[q,g1,...,g
p
]
t
，g
p
＝-((m
l
(d
l
)-1)/2)
p
,...,(-1)
p
,0,1
p
,...,((m
l
(d
l
)-1)/2)
p
，q＝-(m
l
(d
l
)-1)/2,...,-1,1,1,...,(m
l
(d
l
)-1)/2，[
·
]
t
表示转置运算。
[0018]
2-3、原始信号与滤波信号的差值为剩余信号可表示为
[0019][0020]
进一步地，对原始信号进行包络谱分析，其计算式为
[0021][0022]
式中f{
·
}和f-1
{
·
}分别表示傅里叶变换和傅里叶反变换，h(n)可表示为
[0023][0024]
式中数据长度进一步地，包络解析信号可表示为
[0025][0026]
2-4、采用iidf指标表征包络解析信号中包含轴承故障信息的丰富性，计算式为
[0027][0028]
式中γ1＝f
reb-f
t-2fv，γ2＝f
v-2f
t
，γ3＝f
reb-2f
t
，γ4＝(h1+1)f
reb-f
t-2fv，γ5＝h1f
reb-f
t-(h2+1)fv，γ6＝h1f
reb
+f
t-h2fv，γ7＝h1f
reb-f
t
，f
t
＝0.02f
reb
。f
reb
表示轴
承故障特征阶次，fv表示轴承故障特征相关的调制频率，h1和h2分别表示轴承故障特征阶次和调制分量特征阶次。χ1可表示为
[0029][0030]
式中mad＝1.4826
·
κ(|χ
2-κ(χ2)|)，κ(
·
)表示取中值，|
·
|表示取绝对值，通过公式(8)可获得包络解析信号中包含轴承故障信息丰富性的指标
[0031]
2-5、如果floor(s
l+1
/2)《＝1满足，则输出优化滤波长度参数m
op
，可表示为
[0032][0033]
式中argmax{
·
}表示返回最大时对应的优化参数m
op
；
[0034]
如果floor(s
l+1
/2)《＝1不满足，进一步缩小滤波长度m
l+1
(d
l+1
)的取值范围和精细化步长s
l+1
，一种自适应阈值被提出，其可表示为
[0035][0036]
式中0.5《μ《1表示幅值系数，max{
·
}为取最大值操作，min{
·
}为取最小值操作。
[0037]
2-6、分解等级l+1时的参数范围[a
l+1
→
p
,b
l+1
→
p
]，其可表示为
[0038][0039]
式中γ[
·
]表示获取大于阈值th
l
时对应的参数范围[a
l+1
→
p
,b
l+1
→
p
]。进一步地，采用参数增量s
l+1
划分确定的参数范围[a
l+1
→
p
,b
l+1
→
p
]，其中s
l+1
＝s
l
/ε，m
l+1
(d
l+1
)＝a
l+1
→
p
,a
l+1
→
p
+s
l+1
,...,floor((b
l+1
→
p-a
l+1
)/s
l+1
)；返回步骤2-2，进行迭代操作，直至满足公式(10)，并输出优化滤波长度参数m
op
。
[0040]
步骤3：基于确定的优化参数m
op
和固定滤波阶次p，采用savitzky-golay滤波器对原始信号iasi进行滤波，通过公式(4)获得最优剩余信号r
op
。
[0041]
步骤4：通过公式(7)对r
op
进行包络谱分析，以揭示轴承故障特征。
[0042]
本发明的有益效果是：
[0043]
(1)本发明提出一种自适应参数划分结构(pds)，解决了savitzky-golay滤波器依赖于经验确定优化参数的缺陷；
[0044]
(2)本发明提出自适应参数划分结构，其具有兼顾低计算成本和高精度的优势，在工业大数据的工况下，自适应参数划分结构的优势更加明显；
[0045]
(3)本文发明通过合理设置ε、μ和s
l
可提高算法计算效率，并可增强pds的适应性；
[0046]
(4)本发明基于自适应参数的savitzky-golay滤波器消除了编码器安装误差，增强了轴承故障特征。
附图说明
[0047]
图1为本发明方法流程示意图；
[0048]
图2为参数自适应结构(pds)示意图；
[0049]
图3为实施例1中仿真信号的时域波形，其中图(a)为编码器安装误差波形，图(b)为具有噪声的轴承故障波形，图(c)为叠加信号，图(d)是图(c)的包络谱；
[0050]
图4为实施例1中不同理论获得的结果，其中图(a)为cpw处理后的波形，图(b)为图(a)的包络谱，图(c)为快速谱相干算法处理结果，图(d)为本文所提方法确定最优滤波参数，图(e)为所提方法的剩余信号，图(f)为图(e)的包络谱；
[0051]
图5为实施例2中实验台；
[0052]
图6为实施例2中采集的ias信号图(a)和对应的包络谱图(b)；
[0053]
图7为实施例2中不同理论获得的结果，其中图(a)为cpw处理后的波形，图(b)为图(a)的包络谱，图(c)为快速谱相干算法处理结果，图(d)为本文所提方法确定最优滤波参数，图(e)为所提方法的剩余信号，图(f)为图(e)的包络谱。
具体实施方式
[0054]
下面结合本发明实施例中的附图，对本文发明实例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员没有做出创造性劳动前提下所获得的所有其他实施例，都属于本文发明保护的范围，本实施例中如无特殊说明的方法，均为常规方法。
[0055]
在下面的描述中阐述了很多具体细节以便充分理解本发明，但是本发明还可以采用其他不同于在此描述的其他方式来实施，本领域技术人员可以在不违背本发明内涵的情况下做类似推广，因此本发明不受下面公开的具体实施例的限制。
[0056]
实施例1：本实例描述的是本发明方法用于仿真瞬时角速度轴承外圈故障特征的提取，如图3所示，具体过程包括：
[0057]
仿真分析采用的轴承外圈故障信号模型计算式为：
[0058][0059][0060]
式中w(θ)表示平均角速度，wo(θ)表示编码器安装误差，ξ表示阻尼系数，fn表示固有频率，ψ＝θ-jθ-τj，角度序列θ＝2π/n,4π/n,6π/n...，n(θ)表示编码器测量噪声，ρ＝δr/r是几何中心和旋转中心的偏心距离的比值，r为编码器孔直径，δr是几何中心和旋转中心的偏心距离，β表示编码器旋轴与转轴之间的倾斜角度，初始角度θe∈[φe,2kπ+φe]，θ
t
∈[φ
t
,2kπ+φ
t
]，c表示轴承故障冲击幅值。在本仿真信号中，n＝5000，c＝0.0002，θ＝n/3.56，f
reb
＝3.56
×
，《τj》＝0，max{τj}＝0.0071rad，fn＝50
×
，
×
表示特征阶次，ξ＝0.03，w(θ)＝5rad/s，ρ＝0.01，β＝0.03，信噪比snr＝-1，fv＝0，h1＝3，h2＝0，a1→1＝5，b1→1＝n，s1＝10，p＝3，ε＝2，μ＝0.5，φe＝0，φ
t
＝0。
[0061]
步骤1：仿真编码器安装误差波形如图3(a)所示，轴承外圈故障和测量噪声叠加波形如图3(b)所示，图3(a)和3(b)的叠加如图3(c)，对应的包络谱如图3(d)所示。可以发现，编码器安装误差谱线占主导地位，轴承故障特征谱线无法有效辨识。因此，抑制编码器安装误差对于揭示轴承故障特征是至关重要的。
[0062]
步骤2：获得最优的参数m
op
；
[0063]
2-1、基于m1(d1)＝[a1→1,b1→1]，初始参数增量s1＝10，滤波阶次p＝3，滤波系数
如图2中level＝1所示，计算式为
[0064][0065]
式中d
t
＝[q,g1,...,g
p
]
t
，g
p
＝-((m
l
(d
l
)-1)/2)
p
,...,(-1)
p
,0,1
p
,...,((m
l
(d
l
)-1)/2)
p
，q＝-(m
l
(d
l
)-1)/2,...,-1,1,1,...,(m
l
(d
l
)-1)/2，[
·
]
t
表示转置运算，d1＝1,2,...,length((a1→
1-b1→1)/s1)。进一步地，基于滤波系数采用savitzky-golay滤波器对原始信号iasi信号进行滤波，其计算式为
[0066][0067]
进一步地，对应的剩余信号可表示为
[0068][0069]
进一步地，对剩余信号进行包络谱分析，计算式为
[0070][0071]
2-2、采用iidf指标评估包络解析信号中包含轴承故障信息的丰富性，其可表示为
[0072][0073]
式中fv＝0，h1＝3，h2＝0，f
reb
＝3.56
×
，f
t
＝0.02
·freb
＝0.0712，γ1＝f
reb-f
t-2fv＝3.4888，γ3＝f
reb-2f
t
＝3.546，γ7＝h1f
reb-f
t
＝h1f
reb-0.0712，χ1的计算式为
[0074][0075]
式中mad＝1.4826
·
κ(|χ
2-κ(χ2)|)，κ(
·
)表示中值操作，|
·
|表示绝对值操作；
[0076]
2-3、基于s2＝s1/ε＝5，则floor(s2/2)《＝1不满足，基于图2所示的pds结构，进一步缩小滤波长度m1(d1)的取值范围和精细化m1(d1)，μ＝0.5和iidf
1d1
用于获取自适应阈值方法可表示为
[0077][0078]
2-4、获得分解等级level＝2时参数范围[a
s2
→
p
,a
e2
→
p
]，其可表示为
[0079][0080]
2-5、采用参数增量s2＝s1/ε划分确定的参数范围[a2→
p
,b2→
p
]，其中ε＝2，s2＝s1/2＝5，m2(d2)＝a2→
p
,a2→
p
+s2,...,floor((b2→
p-a2→
p
)/s2)。
[0081]
2-6、基于图2所示pds结构，如果floor(s
l+1
/2)《＝1满足，获得最大时对应的优化参数m
op
，其计算式为
[0082][0083]
如图4(d)所示，获取的优化参数m
op
＝2471。
[0084]
步骤3：基于获取的优化参数m
op
＝2471和滤波阶次p＝3，采用savitzky-golay滤波器对原始信号进行滤波，其剩余信号如图4(e)所示，可以发现编码器安装误差被有效抑制。
[0085]
步骤4：对步骤3获得的剩余信号进行包络谱分析，如图4(f)所示，可以发现，轴承故障冲击分量被有效增强。
[0086]
步骤5：为进一步验证本文所提方法的有效性，采用传统cpw算法对原始信号进行分析，处理后的信号和对应的频谱如图4(a)和4(b)所示，此外，采用快速谱相干算法对原始信号分析，循环频率α
max
＝30hz，窗宽nw＝512，结果如图4(c)所示，可见与轴承故障相关的特征谱线均无法有效辨识。
[0087]
实施例2：本实施例描述的是本发明方法用于实际滚动轴承外圈故障特征的提取
[0088]
本实施例中采用轴承测试试验台，如图5所示，在实验台上安装一个resatron光学编码器，编码器线数n＝5000，并采用106采样率的picoscope高速采集装置获取对应的角度信息和时间信息。本试验台的轴承类型是nu206e(nb＝13,eb＝9.525,e
p
＝46,α＝0)，为模拟轴承外圈故障，采用线切割方式在外圈上加工宽度大约0.5mm，深度约为0.5mm的槽；由以下计算式得到轴承外圈故障特征频率f
reb
为5.15
×
。
[0089][0090]
步骤1：获得的iasi轴承外圈故障波形如图6(a)所示；
[0091]
步骤2：直接对原始iasi信号进行包络阶次谱分析，如图6(b)所示。可以发现，编码器安装误差的谱线占主导地位，然而轴承故障特征谱线无法有效辨识。因此，编码器安装误差的干扰需要被抑制以增强轴承故障特征。
[0092]
步骤3：通过pds结构确定优化滤波参数m
op
[0093]
3-1、基于m1(d1)＝[a1→1,b1→1]＝[5,5000]，滤波阶次p＝3，m1(d1)的增量s1＝20，获得滤波系数其计算式为
[0094][0095]
式中d
t
＝[q,g1,...,g
p
]
t
，g
p
＝-((m
l
(d
l
)-1)/2)
p
,...,(-1)
p
,0,1
p
,...,((m
l
(d
l
)-1)/2)
p
，q＝-(m
l
(d
l
)-1)/2,...,-1,1,1,...,(m
l
(d
l
)-1)/2，[
·
]
t
表示转置运算，d1＝1,2,...,length((a1→
1-b1→1)/s1)。
[0096]
3-2、基于滤波系数采用savitzky-golay滤波器对原始iasi信号进行滤波，计算式为
[0097][0098]
进一步地，对应的剩余信号可表示为
[0099][0100]
进一步地，对剩余信号进行包络谱分析，计算式为
[0101][0102]
进一步地，对剩余信号进行包络谱分析。
[0103]
3-3、采用iidf指标评估包络解析信号中包含轴承故障信息的丰富性，其可表示为
[0104][0105]
式中fv＝0，h1＝3，h2＝0，f
reb
＝5.15
×
，f
t
＝0.02
·freb
＝0.103，γ1＝f
reb-f
t-2fv＝5.047，γ3＝f
reb-2f
t
＝4.944，γ7＝h1f
reb-f
t
＝h1f
reb-0.103，χ1的计算式为。
[0106][0107]
式中mad＝1.4826
·
κ(|χ
2-κ(χ2)|)，κ(
·
)表示中值操作，|
·
|表示绝对值操作。
[0108]
3-5、基于ε＝2，s2＝s1/ε＝10，则floor(s2/2)《＝1不满足，基于图2所示的pds结构，进一步缩小滤波长度m1的取值范围和精细化m1，μ＝0.5和iidf
1d1
用于获取自适应阈值
[0109][0110]
3-6、获得分解等级l＝2时的参数范围[a2→
p
,b2→
p
]，如图2中level＝2所示，其可表示为
[0111][0112]
3-7、采用参数增量s2＝s1/ε划分确定的参数范围[a2→
p
,b2→
p
]，如图2中level＝2，其中ε＝2，s2＝s1/2＝10，m2＝a2→
p
,a2→
p
+s2,...,floor((b2→
p-a2→
p
)/s2)。
[0113]
2-6、基于图2的pds结构所示，如果floor(s
l+1
/2)《＝1满足，获得最大时对应的优化参数m
op
，其计算式为
[0114]
[0115]
如图7(d)所示，确定的优化参数m
op
＝147；
[0116]
步骤3：基于p＝3和m
op
＝147，采用savitzky-golay滤波器对原始信号进行滤波，其剩余信号如图7(e)所示。可以发现编码器安装误差被有效抑制；
[0117]
步骤4：对步骤3获得的剩余信号进行包络谱分析，如图7(f)所示，可以发现，轴承故障冲击分量被有效增强；
[0118]
步骤5：为进一步验证本文所提方法的有效性，采用传统cpw算法对原始信号进行分析，处理后的信号和对应的频谱如图7(a)和7(b)所示，此外，采用快速谱相干算法对原始信号分析，循环频率α
max
＝30hz，窗宽nw＝512，结果如图7(c)所示，可见，可以发现，与轴承故障相关的特征谱线均无法有效辨识。
[0119]
本文中应用了具体的实例对本发明的原理及实施方式进行了参数，以上实例的说明只是用于帮助理解本发明及核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变指出，综上，本说明书内容不应理解为本发明的限制。