CN106680716A - 一种基于无位置传感器的永磁电机轴承故障诊断的方法 - Google Patents

Abstract

本发明公开了一种基于无位置传感器的永磁电机轴承故障诊断的方法。其实现步骤为：首先，从电机控制芯片中获得电压与电流，并利用无位置传感器算法获得电机的转子位置角与转速，然后，通过滑动平均滤波器提取出转速中的交流成分，并根据位置角度对其进行角域重采样，最后，对重采样速度信号进行频域分析，根据频域信息判断轴承故障是否发生。本发明计算量小，不要求额外的传感器，可以对非稳态条件下运行的电机进行故障诊断，并且比起普通的基于电流信号的故障诊断方法具有更高的精度。

Description

一种基于无位置传感器的永磁电机轴承故障诊断的方法

技术领域

本发明属于故障诊断技术领域，尤其是一种永磁电机轴承故障诊断的方法。

背景技术

永磁电机具有体积小，功率密度高的优点而在精密机床、交通运输、国防工业中得到了广泛的应用。然而，电机在长时间的运行后，可能会出现各种故障。在故障初期，若没有发现故障并进行合适的处理，故障会进一步扩大，并造成整个系统的非预期停机，并因此会带来巨大的经济损失甚至给工作人员带来生命危险，因此，对一些重要场合的电机进行故障诊断具有重大的意义。

永磁电机的故障主要分为定子匝间短路故障、转子退磁故障、轴承故障及气隙偏心故障。其中，轴承故障占到了所有故障的40％。

电机的轴承发生故障时，每当滚珠与滚道的故障点之间相互接触时，都会产生特定周期的脉冲冲击力，进而导致电机产生噪声与振动。同时，电机的负载转矩也会产生周期性的扰动，进而使得电机的电流，电压等发生变化。因此，利用加速度传感器采集振动信号，或者用电流传感器采集电流信号，并分析其频谱。

然而，加速度传感器存在附加成本，噪音较大，安装不便等缺点。而电机电流信号分析方法中的基波信号的幅值比故障分类的幅值要大很多，容易对故障信号的采集产生干扰。而且，当电机处于非平稳工况时，传统基于稳态数据的频谱分析方法无法得到有意义的结果，而基于时频分析的方法计算复杂，难以进行快速运算。

发明内容

为了克服现有技术的不足，本发明提供了一种基于无位置传感器的永磁电机轴承故障诊断的方法。该方法简单，实时性强，可以有效地实现轴承故障诊断。

本发明的目的是通过以下技术方案实现的：一种适用于永磁电机轴承故障诊断方法，包括以下步骤：

(1)通过测量逆变器直流母线电压U_dc与相电流(i_a，i_b，i_c)，在控制芯片内部得到αβ轴电压(v_α，v_β)与电流(i_α，i_β)，并通过无位置传感器算法计算得到电机的转子位置角θ(t)与转速ω(t)。

(2)利用滑动平均滤波器，分别提取出转速ω(t)中的直流成分ω_DC(t)与交流成分ω_AC(t)。

(3)对转子位置角θ(t)根据ω_DC(t)进行矫正，并对ω_AC(t)进行角域重采样，得到重采样后转速ω_RS(θ)。

(4)计算轴承故障信号特征次数k_fault，并计算ω_RS(θ)的部分频谱F(k)以及故障阈值T_fault，当频谱中k次谐波处频谱幅值超过阈值T_fault，则证明发生了轴承故障。

进一步的，所述步骤(1)中，θ(t)与ω(t)通过以下方法得到：

(2.1)根据测量得到的电机相电流与直流母线电压，在控制芯片内部得到电机的αβ轴电流i_α,i_β与电压v_α,v_β，并记i＝[i_α i_β]^T，v＝[v_α v_β]^T，上角标^T表示矩阵的转置

(2.2)构造下述观测器，求出电机的αβ轴反电势

其中，为辅助变量，G＝gI为观测器的增益矩阵，为ω_re的函数，A₂₂＝ω_reJ，均为与参数有关的矩阵，其中，为常数矩阵，R为定子电阻，L_d为d轴电感，L_q为q轴电感，可以通过测量或参数辨识等手段得到，ω_re为电机的实际电气角速度，由下文中提到速度辨识环节反馈获得。

(2.3)转子位置角θ(t)通过下式获得

其中，atan2表示值域为[0,2π)的取反正切函数。

(2.4)电机转速ω(t)通过以下步骤获得

(2.4.1)将反电势通过下式，标幺化为

(2.4.2)构建自适应观测器：

其中，为观测器的状态变量，为辨识得到的电气角速度，H＝hI(h＞0)为观测器增益，为观测器误差量，为常数矩阵。

(2.4.3)通过PI控制器，得到辨识的电机电气角速度

其中，k_p,k_i为PI控制器的系数。

(2.4.4)电机转速其中n_pp为电机极对数，而(1.2)中的ω_re使用辨识得到的电气角速度

所述步骤(2)中，ω_DC(t)与ω_AC(t)通过以下方法得到：

(3.1)构造滑动平均滤波器，求出ω_DC(t)与ω_AC(t)，步骤如下：

(3.1.1)生成长度为2L的数组S，并令其初值全部为0。L为一较大的正整数。定义变量index，作为数组下标的计数值，初值为0。定义转速直流部分ω_DC(t)，初值为0。

(3.1.2)每当步骤(2)中的无位置传感器算法一个步长T结束时，按照下式更新ω_DC(t)：

其中，ω_DC(t-LT)为当前步长计算得到的电机转速直流部分，其表示2LT时刻之前到0时刻的速度均值，其与当前时刻速度之间存在LT的滞后，ω_DC(t-(L+1)T)为上一步计算得到的电机转速直流部分，ω(t)为此时无位置传感器算法计算得到的转速，S(index)为数组S中，下标为index的值。

(3.1.3)计算转速交流部分ω_AC(t-LT)＝ω(t-LT)-ω_DC(t-LT)

(3.1.4)将当前转速ω(t)存入数组S下标为index的位置，然后将下标index加1，当index等于数组长度2L时，将index重新置0。

(3.1.5)重复(3.1.2)至(3.1.4)，直到无位置传感器停止计算。

(3.2)将转速交流部分ω_AC(t-LT)前L个数据剔除后，得到对应LT时刻至t_All-LT时刻的转速交流部分。

所述步骤(3)中，ω_RS(θ)通过以下方法计算得到：

(4.1)随着速度的变化，步骤(1)中获得的转子位置角θ(t)存在不同的相位延时，通过查表法可以通过ω_DC(t)的值，获得矫正后的转子位置角θ_adj(t)。

(4.2)为了与转速交流部分ω_AC(t)在时域上对齐，将步骤(4.1)中矫正后的转子位置角θ_adj(t)剔除前L个数据，记作θ_t(t)。

(4.3)构建等角度间隔序列θ_a(n)＝nθ_δ，其中，n∈[0,N·M)为非负整数，为给定的角度间隔，其在角度域中等效为时域中的采样时间，需满足香农采样定理；M为需要进行分析的旋转圈数，其在角度域中等效为时域中的采样时长，决定了角度域的频谱分辨率。

(4.4)对于θ_a(n)中一个特定点θ_a(k)＝kθ_δ，在序列θ(t)中找到一个时刻t_s满足下式：

θ_t(t_s)＜θ_a(k)＜θ_t(t_s+T)

对应的ω_AC(t)角域重采样值如下式所示：

(4.5)对序列θ_a(n)中的每一个点均施加(4.3)中的处理，得到角域重采样后的转速ω_RS(θ)。

所述步骤(4)中，轴承故障信号特征次数k_fault，频谱F(k)以及故障阈值T_fault通过以下方法计算得到：

(5.1)根据轴承故障位置的不同，计算故障信号特征次数：

其中，k_out,k_in,k_ball,k_cage分别为外滚道故障，内滚道故障，滚动体故障，保持架故障时故障信号特征次数，D_p,D_b,N,β分别为故障轴承的节径，滚动体直径，滚动体数量，接触角。

(5.2)为了节省计算量，针对角度域重采样数据ω_RS(θ)，如图6所示，构建一组阶数为k的Goertzel滤波器来计算频谱F(k)，其中k∈{0,θ_δ,2θ_δ,3θ_δ,4θ_δ,...,K_maxθ_δ}，K_maxθ_δ为预设的最大计算次数，应至少选择大于故障特征次数k_fault两倍的数。下式中，s_k[-2]＝s_k[-1]＝0：

s_k[n]＝x[n]+2cos(k)s_k[n-1]-s_k[n-2]

y_k[n]＝s_k[n]-e^-jks_k[n-1]

F(k)＝y_k[N·M]

(5.3)将频谱F(k)的幅值表示为常用对数的形式F_log(k)，利用以下步骤计算故障阈值T_fault：

(5.3.1)定义滑窗A_k＝[max(0,k-W),min(N·M,k+W)]，将窗口A_k内频谱值的平均值作为中心频率点k的频谱值，构建平均值频谱F_mean(k)。

(5.3.2)定义滑窗B_k＝[max(0,k-L),min(N·M,k+L)]，将窗口B_k内频谱值的中位数F_median(k)＝median{F_mean(j)}|_{max(-k,-L)≤j≤min(L,N·M-k)}作为中心频率点k的频谱值，构建中位数频谱F_median(k)。

(5.3.3)取中位数频谱F_median(k)的最大值，即为故障阈值T_fault＝max[F_median(k)]。

本发明的有益效果列举如下：

1)可以仅使用电机运行时的电流传感器与直流母线电压传感器，无需额外的传感器，降低了成本；

2)算法可以一边进行数据的采集，一边在线进行计算，而无需进行数据的存储，计算简单，易于实施；

3)可以用较小的计算量处理电机的非平稳信号；

4)能够自动计算故障阈值，判断故障是否发生。

附图说明

图1为本发明实施方案的一个结构示意图；

图2为本发明的一个实施步骤框图；

图3为转子角度与电机转速观测器的示意图；

图4为滑动平均滤波器的示意图；

图5为角域重采样的示意图；

图6为频谱计算的示意图；

图7为故障阈值计算的示意图。

具体实施方式

为了更加具体地描述本发明，下面结合附图及具体实施方式对本发明的技术方案进行详细说明。

图1给出了本发明的实验装置图。通过电流传感器与电压传感器获取电机的三相电流(i_a，i_b，i_c)与直流母线电压U_dc。对于星型接法的电机，其三相电流之和为零，所以也可只采集两相电流，并据此计算出第三相电流的值。通过控制芯片获得αβ轴电压(v_α，v_β)与电流(i_α，i_β)，并送入到本发明提出的故障诊断算法中，即可判断出电机轴承故障与否。

本发明具体的实施步骤如图2所示：

(1)根据测量得到的电机相电流与直流母线电压，在控制芯片内部得到电机的αβ轴电流i_α,i_β与电压v_α,v_β，并记i＝[i_α i_β]^T，v＝[v_α v_β]^T，上角标^T表示矩阵的转置

(2)构造下述观测器，求出电机的αβ轴反电势

其中，为辅助变量，G＝gI为观测器的增益矩阵，为ω_re的函数，A₂₂＝ω_reJ，均为与参数有关的矩阵，其中，为常数矩阵，R为定子电阻，L_d为d轴电感，L_q为q轴电感，可以通过测量或参数辨识等手段得到，ω_re为电机的实际电气角速度，由下文中提到速度辨识环节反馈获得。

(3)转子位置角θ(t)通过下式获得

其中，atan2表示值域为[0,2π)的取反正切函数。

(4)电机转速ω(t)通过以下步骤获得

(4.1)将反电势通过下式，标幺化为

(4.2)构建自适应观测器：

其中，为观测器的状态变量，为辨识得到的电气角速度，H＝hI(h＞0)为观测器增益，为观测器误差量，为常数矩阵。

(4.3)通过PI控制器，得到辨识的电机电气角速度

其中，k_p,k_i为PI控制器的系数。

(4.4)电机转速其中n_pp为电机极对数，而(1.2)中的ω_re使用辨识得到的电气角速度

(5)为了计算转速的直流部分与交流部分，首先需要生成长度为2L的数组S，并令其初值全部为0，其中L为一较大的正整数；定义变量index，作为数组下标的计数值，初值为0；定义转速直流部分ω_DC(t)，初值为0；

(5.1)如图4所示，步骤(4)每更新一次电机转速ω(t)，便按照下式更新ω_DC(t)：

其中，ω_DC(t-LT)为当前步长计算得到的电机转速直流部分，其表示2LT时刻之前到0时刻的速度均值，其与当前时刻速度之间存在LT的滞后；ω_DC(t-(L+1)T)为上一步计算得到的电机转速直流部分，S(index)为数组S中，下标为index的值，其值在更新前为0，或者保存了滞后2LT时间的电机转速ω(t-2LT)。

(5.2)计算得到的直流成分ω_DC(t-LT)总是比最新转速ω(t)滞后LT的时间，在计算交流成分ω_AC(t)时，需要从数组S中取出时间上与ω_DC(t-LT)对应的ω(t-LT)，再进行计算；

ω_AC(t-LT)＝ω(t-LT)-ω_DC(t-LT)

(5.3)将当前转速ω(t)存入数组S下标为index的位置，然后将下标index加1；当index等于数组长度2L时，将index重新置0，用新数据ω(t)逐个将旧数据ω(t-2LT)替换掉。

(6)针对(2)中的观测器，计算出不同转速下对应的转子位置角θ(t)滞后大小，然后根据转速直流分量ω_DC(t)，对θ(t)进行矫正；然后，为了与转速交流部分ω_AC(t)在时域上对齐，将矫正后的转子位置角剔除前L个数据，记作θ_t(t)。

(7)以矫正后的θ(t)为基准，对ω_AC(t)进行角域重采样，方法如图5所示：

(7.1)构建等角度间隔序列θ_a(n)＝nθ_δ，其中，n∈[0,N·M)为非负整数，为给定的角度间隔，其在角度域中等效为时域中的采样时间，需满足香农采样定理；M为需要进行分析的旋转圈数，其在角度域中等效为时域中的采样时长，决定了角度域的频谱分辨率；

(7.2)对于等角度间隔序列θ_a(n)＝nθ_δ中的任意一个点θ_a(k)＝kθ_δ，在序列θ(t)中找到相邻的两个点θ_t(t_s)与θ_t(t_s+T)，使其满足θ_t(t_s)＜θ_a(k)＜θ_t(t_s+T)；

(7.3)然后，在角度域中的两点(θ_t(t_s),ω_AC(t_s))与(θ_t(t_s+T),ω_AC(t_s+T))之间进行线性插值，求出θ_a(k)对应的重采样后转速；

(7.4)对等角度间隔序列θ_a(n)＝nθ_δ中所有点均进行步骤(7.3)的操作，得到重采样转速ω_RS(θ)；

(8)根据轴承故障位置的不同，计算故障信号特征次数：

其中，k_out,k_in,k_ball,k_cage分别为外滚道故障，内滚道故障，滚动体故障，保持架故障时故障信号特征次数，D_p,D_b,N,β分别为故障轴承的节径，滚动体直径，滚动体数量，接触角。

(9)为了节省计算量，针对角度域重采样数据ω_RS(θ)，如图6所示，构建一组阶数为k的Goertzel滤波器来计算频谱F(k)，其中k∈{0,θ_δ,2θ_δ,3θ_δ,4θ_δ,...,K_maxθ_δ}，K_maxθ_δ为预设的最大计算次数，应至少选择大于故障特征次数k_fault两倍的数。下式中，s_k[-2]＝s_k[-1]＝0：

s_k[n]＝x[n]+2cos(k)s_k[n-1]-s_k[n-2]

y_k[n]＝s_k[n]-e^-jks_k[n-1]

F(k)＝y_k[N·M]

(10)如图7的第1幅子图所示，将频谱F(k)的幅值表示为常用对数的形式F_log(k)，利用以下步骤计算故障阈值T_fault：

(10.1)定义滑窗A_k＝[max(0,k-W),min(N·M,k+W)]，将窗口A_k内频谱值的平均值作为中心频率点k的频谱值，构建平均值频谱F_mean(k),如图7的第2幅子图所示。

(10.2)定义滑窗B_k＝[max(0,k-L),min(N·M,k+L)]，将窗口B_k内频谱值的中位数F_median(k)＝median{F_mean(j)}|_{max(-k,-L)≤j≤min(L,N·M-k)}作为中心频率点k的频谱值，构建中位数频谱F_median(k)，如图7的第3幅子图所示。

(10.3)取中位数频谱F_median(k)的最大值，即为故障阈值T_fault＝max[F_median(k)]，如图7的第4幅子图所示。

(11)在各个故障特征次数k_fault处，比较F_log(k)与故障阈值T_fault的大小，当频谱中k次谐波处频谱幅值超过阈值T_fault，则证明发生了轴承故障。

(12)利用本发明提出的方法，可以仅使用电机运行时的电流传感器与直流母线电压传感器，无需额外的传感器，降低了成本。算法可以一边进行数据的采集，一边在线进行计算，而无需进行数据的存储，计算简单，易于实施。可以用较小的计算量处理电机的非平稳信号。并且能够自动计算故障阈值，判断故障是否发生。

Claims (5)

1.一种基于无位置传感器的永磁电机轴承故障诊断的方法，其特征在于，

(1)通过测量逆变器直流母线电压U_dc与相电流(i_a，i_b，i_c)，在控制芯片内部得到αβ轴电压(v_α，v_β)与电流(i_α，i_β)，并通过无位置传感器算法计算得到电机的转子位置角θ(t)与转速ω(t)；

(2)利用滑动平均滤波器，分别提取出转速ω(t)中的直流成分ω_DC(t)与交流成分ω_AC(t)；

(3)对转子位置角θ(t)根据ω_DC(t)进行矫正，并对ω_AC(t)进行角域重采样，得到重采样后转速ω_RS(θ)；

(4)计算轴承故障信号特征次数k_fault，并计算ω_RS(θ)的部分频谱F(k)以及故障阈值T_fault，当频谱中k次谐波处频谱幅值超过阈值T_fault，则证明发生了轴承故障。

2.根据权利要求1所述的方法，其特征在于，

进一步的，所述步骤(1)中，转子位置角θ(t)与转速ω(t)通过以下方法得到：

(2.1)根据测量得到的电机相电流与直流母线电压，在控制芯片内部得到电机的αβ轴电流i_α,i_β与电压v_α,v_β，并记i＝[i_α i_β]^T，v＝[v_α v_β]^T，上角标T表示矩阵的转置；

(2.2)构造下述观测器，求出电机的αβ轴反电势

$\stackrel{\·}{\mathrm{\ξ}}=({\mathrm{GA}}_{12}+A_{22})\mathrm{\ξ}+({\mathrm{GA}}_{11}-{\mathrm{GA}}_{12}G-A_{22}G)i+{\mathrm{GB}}_{1}v$

$\hat{e}=\mathrm{\ξ}-G\hat{i}$

其中，为辅助变量，G＝gI为观测器的增益矩阵，为ω_re的函数，A₂₂＝ω_reJ，均为与参数有关的矩阵，其中，为常数矩阵，R为定子电阻，L_d为d轴电感，L_q为q轴电感，可以通过测量或参数辨识等手段得到，ω_re为电机的实际电气角速度；

(2.3)转子位置角θ(t)通过下式获得

$\mathrm{\θ}\left(t\right)=atan2\left(\frac{-{\hat{e}}_{\mathrm{\α}}}{{\hat{e}}_{\mathrm{\β}}}\right)$

其中，atan2表示值域为[0,2π)的取反正切函数；

(2.4)电机转速ω(t)通过以下步骤获得：

(2.4.1)将反电势通过下式，标幺化为

${\hat{e}}_n=\frac{1}{\sqrt{{\hat{e}}_{\mathrm{\α}}^2+{\hat{e}}_{\mathrm{\β}}^2}}\left[\begin{array}{c}{\hat{e}}_{\mathrm{\α}}\\ {\hat{e}}_{\mathrm{\β}}\end{array}\right],$

(2.4.2)构建自适应观测器：

${\stackrel{\·}{\tilde{e}}}_n={\widehat{\mathrm{\ω}}}_{re}J{\tilde{e}}_n+H\mathrm{\ϵ},$

其中，为观测器的状态变量，为辨识得到的电气角速度，H＝hI(h＞0)为观测器增益，为观测器误差量，为常数矩阵；

(2.4.3)通过PI控制器，得到辨识的电机电气角速度

${\widehat{\mathrm{\ω}}}_{re}=(k_p+\frac{k_i}{s})\left({\mathrm{\ϵ}}^{T}J{\tilde{e}}_n\right)$

其中，k_p,k_i为PI控制器的系数；

(2.4.4)电机转速其中n_pp为电机极对数，而(2.2)中的ω_re使用辨识得到的电气角速度

3.根据权利要求1或2所述的方法，其特征在于，

所述步骤(2)中，ω_DC(t)与ω_AC(t)通过以下方法得到：

(3.1)构造滑动平均滤波器，求出ω_DC(t)与ω_AC(t)，步骤如下：

(3.1.1)生成长度为2L的数组S，并令其初值全部为0，L为一较大的正整数，定义变量index，作为数组下标的计数值，初值为0，定义转速直流部分ω_DC(t)，初值为0；

(3.1.2)每当步骤(2)中的无位置传感器算法一个步长T结束时，按照下式更新ω_DC(t)：

${\mathrm{\ω}}_{DC}(t-LT)={\mathrm{\ω}}_{DC}(t-(L+1\left)T\right)+\frac{\[\mathrm{\ω}\left(t\right)-S\left(index\right)\]}{2L}$

其中，ω_DC(t-LT)为当前步长计算得到的电机转速直流部分，其表示2LT时刻之前到0时刻的速度均值，其与当前时刻速度之间存在LT的滞后，ω_DC(t-(L+1)T)为上一步计算得到的电机转速直流部分，ω(t)为此时无位置传感器算法计算得到的转速，S(index)为数组S中，下标为index的值；

(3.1.3)计算转速交流部分ω_AC(t-LT)＝ω(t-LT)-ω_DC(t-LT)；

(3.1.4)将当前转速ω(t)存入数组S下标为index的位置，然后将下标index加1，当index等于数组长度2L时，将index重新置0；

(3.1.5)重复(3.1.2)至(3.1.4)，直到无位置传感器停止计算；

(3.2)将转速交流部分ω_AC(t-LT)前L个数据剔除后，得到对应LT时刻至t_All-LT时刻的转速交流部分。

4.根据权利要求1-3任一项所述的方法，其特征在于，

所述步骤(3)中，ω_RS(θ)通过以下方法计算得到：

(4.1)随着速度的变化，步骤(1)中获得的转子位置角θ(t)存在不同的相位延时，通过查表法可以通过ω_DC(t)的值，获得矫正后的转子位置角θ_adj(t)；

(4.2)为了与转速交流部分ω_AC(t)在时域上对齐，将步骤(4.1)中矫正后的转子位置角θ_adj(t)剔除前L个数据，记作θ_t(t)；

(4.3)构建等角度间隔序列θ_a(n)＝nθ_δ，其中，n∈[0,N·M)为非负整数，为给定的角度间隔，其在角度域中等效为时域中的采样时间，需满足香农采样定理；M为需要进行分析的旋转圈数，其在角度域中等效为时域中的采样时长，决定了角度域的频谱分辨率；

(4.4)对于θ_a(n)中一个特定点θ_a(k)＝kθ_δ，在序列θ(t)中找到一个时刻t_s满足下式：

θ_t(t_s)＜θ_a(k)＜θ_t(t_s+T)，

对应的ω_AC(t)角域重采样值如下式所示：

${\mathrm{\ω}}_{RS}\left(k\right)={\mathrm{\ω}}_{AC}\left(t_s\right)+\frac{{\mathrm{\ω}}_{AC}(t_s+T)-{\mathrm{\ω}}_{AC}\left(t_s\right)}{{\mathrm{\θ}}_t(t_s+T)-{\mathrm{\θ}}_t\left(t_s\right)}\[{\mathrm{\θ}}_a\left(k\right)-{\mathrm{\θ}}_t\left(t_s\right)\];$

(4.5)对序列θ_a(n)中的每一个点均施加(4.3)中的处理，得到角域重采样后的转速ω_RS(θ)。

5.根据权利要求1-4任一项所述的方法，其特征在于，

所述步骤(4)中，轴承故障信号特征次数k_fault，频谱F(k)以及故障阈值T_fault通过以下方法计算得到：

(5.1)根据轴承故障位置的不同，计算故障信号特征次数：

$k_{fault}=\left\{\begin{array}{c}{k}_{out}=\frac{N}{2}(1-\frac{D_{b}cos\mathrm{\β}}{D_p})\\ k_{in}=\frac{N}{2}(1+\frac{D_b\cos \mathrm{\β}}{D_p})\\ k_{ball}=\frac{D_p}{D_b}(1-{\left(\frac{D_b\cos \mathrm{\β}}{D_p}\right)}^2)\\ k_{cage}=\frac{1}{2}(1-\frac{D_b\cos \mathrm{\β}}{D_p})\end{array}\right.,$

其中，k_out,k_in,k_ball,k_cage分别为外滚道故障，内滚道故障，滚动体故障，保持架故障时故障信号特征次数，D_p,D_b,N,β分别为故障轴承的节径，滚动体直径，滚动体数量，接触角；

(5.2)为了节省计算量，针对角度域重采样数据ω_RS(θ)，如图6所示，构建一组阶数为k的Goertzel滤波器来计算频谱F(k)，其中k∈{0,θ_δ,2θ_δ,3θ_δ,4θ_δ,...,K_maxθ_δ}，K_maxθ_δ为预设的最大计算次数，应至少选择大于故障特征次数k_fault两倍的数，下式中，s_k[-2]＝s_k[-1]＝0：

$\begin{array}{c}{s}_k\[n\]=x\[n\]+2\cos \left(k\right)s_k\[n-1\]-s_k\[n-2\]\\ y_k\[n\]=s_k\[n\]-e^{-jk}{s}_k\[n-1\]\\ F\left(k\right)=y_k\[N\·M\]\end{array},$

(5.3)将频谱F(k)的幅值表示为常用对数的形式F_log(k)，利用以下步骤计算故障阈值T_fault：

(5.3.1)定义滑窗A_k＝[max(0,k-W),min(N·M,k+W)]，将窗口A_k内频谱值的平均值作为中心频率点k的频谱值，构建平均值频谱F_mean(k)；

(5.3.2)定义滑窗B_k＝[max(0,k-L),min(N·M,k+L)]，将窗口B_k内频谱值的中位数F_median(k)＝median{F_mean(j)}|_{max(-k,-L)≤j≤min(L,N·M-k)}作为中心频率点k的频谱值，构建中位数频谱F_median(k)；

(5.3.3)取中位数频谱F_median(k)的最大值，即为故障阈值T_fault＝max[F_median(k)]。