循环双谱及在周期平稳类故障中的应用

摘要

论述了未知循环频率的周期平稳信号循环双谱的估计方法；研究了在循环双谱的循环累积量计算中涉及变量的简化存储方法，提出了该变量矩阵是对称阵，通过算法可以三角阵的元素来表述，克服了循环双谱传统估计方法计算量较大的缺陷，提高了运算效率；提出了循环双谱对调相故障信号的分析能力，以及对加性噪声的处理能力，仿真并验证了该方法的有效性，并将其应用于旋转机械状态分析。

正文

数学上通常把均值与自相关函数呈现周期性或几乎为周期性的信号称为周期性平稳随机信号。由于旋转机械的旋转周期性以及与外界的作用, 其振动信号可以看成是周期过程与随机过程的组合, 表现为周期平稳性 ^[1,2]。而基于信号的高阶循环统计量的循环双谱能够有效地利用高阶循环统计量抑制噪声。对于齿轮箱的调制故障信号, 非平稳过程循环双谱比双谱分析更为有效, 但循环双谱估计的实时性差, 所需的数据量大, 短数据将影响数据估计精度 ^[3]。故需对循环双谱实现的时效进行分析。

笔者对循环双谱的实现方法进行了探索, 研究了未知循环频率的周期平稳信号循环双谱的估计方法。由于循环双谱的计算量比较大, 计算时间长, 故研究了能够减少计算时间的方法和循环双谱在调相信号中的应用, 以及对含加性噪声的信号处理, 并应用循环双谱来分析旋转机械状态。

《1 循环双谱理论与实现》

1 循环双谱理论与实现

《1.1基本实现方法》

1.1基本实现方法

根据循环平稳定义, 信号的一阶、二阶、三阶循环矩分别表示为 ^[4]:

$\begin{array}{l} Μ_{1 x}^{α} = 〈 x (t) e^{- j 2 π α t} 〉_{t} (1) \\ Μ_{2 x}^{α} (τ) = 〈 x (t) x (t + τ) e^{- j 2 π α t} 〉_{t} (2) \\ Μ_{3 x}^{α} (τ_{1}, τ_{2}) = 〈 x (t) x (t + τ_{1}) x (t + τ_{2}) e^{- j 2 π α t} 〉_{t} (3) \end{array}$ $\begin{array}{l} Μ_{1 x}^{α} = 〈 x (t) e^{- j 2 π α t} 〉_{t} (1) \\ Μ_{2 x}^{α} (τ) = 〈 x (t) x (t + τ) e^{- j 2 π α t} 〉_{t} (2) \\ Μ_{3 x}^{α} (τ_{1}, τ_{2}) = 〈 x (t) x (t + τ_{1}) x (t + τ_{2}) e^{- j 2 π α t} 〉_{t} (3) \end{array}$

式中α表示循环频率, 〈·〉_t表示时间平均, x (t) 是采样信号, τ_i表示采样的延时。

循环累积函数可用较低阶循环矩函数表示 ^[5]:

$\begin{array}{l} C_{k x}^{α} (τ_{1} ‚ \dots, τ_{k - 1}) = \sum_{U_{p = 1}^{q} Ι_{p} = Ι} [(- 1)^{q - 1} (q - 1)! \\ \sum_{α + \dots α_{q} = β} \prod_{p = 1}^{q} Μ_{n_{p} x}^{α_{p}} (τ Ι_{p})] (4) \end{array}$ $\begin{array}{l} C_{k x}^{α} (τ_{1} ‚ \dots, τ_{k - 1}) = \sum_{U_{p = 1}^{q} Ι_{p} = Ι} [(- 1)^{q - 1} (q - 1)! \\ \sum_{α + \dots α_{q} = β} \prod_{p = 1}^{q} Μ_{n_{p} x}^{α_{p}} (τ Ι_{p})] (4) \end{array}$

式中 $\sum_{U_{p = 1}^{q} Ι_{p} = Ι}$ $\sum_{U_{p = 1}^{q} Ι_{p} = Ι}$ 表示在指数符集I={0, 1, …, k-1}的所有无交连的非空集合分割 (1≤q≤k) 内的求和, n_p=|I_p|表示分割I_p中的元素个数。

信号的三阶循环累积量可以通过信号的一阶、二阶、三阶循环矩求得:

$\begin{array}{l} C_{3 x}^{α} (τ_{1}, τ_{2}) = Μ_{3 x}^{α} (τ_{1}, τ_{2}) + 2 (Μ_{1 x}^{α})^{3} - \\ Μ_{1 x}^{α} [Μ_{2 x}^{α} (τ_{1}) + Μ_{2 x}^{α} (τ_{2}) + Μ_{2 x}^{α} (τ_{2} - τ_{2})] (5) \end{array}$ $\begin{array}{l} C_{3 x}^{α} (τ_{1}, τ_{2}) = Μ_{3 x}^{α} (τ_{1}, τ_{2}) + 2 (Μ_{1 x}^{α})^{3} - \\ Μ_{1 x}^{α} [Μ_{2 x}^{α} (τ_{1}) + Μ_{2 x}^{α} (τ_{2}) + Μ_{2 x}^{α} (τ_{2} - τ_{2})] (5) \end{array}$

三阶累积量具有优良性质:对任何平稳的高斯或非高斯噪声以及非平稳的高斯噪声信号, 三阶累积量的值均为零, 从而使三阶累积量具有消除信号中噪声的功能, 是分析被污染信号特征的较好方法。三阶循环平稳的应用主要是基于三阶循环累积量的循环双谱表示, 循环双谱定义为

$S_{3 x}^{α} (ω_{1} ‚ ω_{2}) = \sum_{τ_{1} = - \infty}^{\infty} \sum_{τ_{2} = - \infty}^{\infty} C_{3 x}^{α} (τ_{1}, τ_{2}) e^{- j (τ_{1} ω_{1} + τ_{2} ω_{2})} (6)$ $S_{3 x}^{α} (ω_{1} ‚ ω_{2}) = \sum_{τ_{1} = - \infty}^{\infty} \sum_{τ_{2} = - \infty}^{\infty} C_{3 x}^{α} (τ_{1}, τ_{2}) e^{- j (τ_{1} ω_{1} + τ_{2} ω_{2})} (6)$

循环双谱表示循环频率下循环累积量的谱对双频率轴的分布是三维的。一般的循环双谱只能表示一个循环频率下的信号特征, 即信号对某一循环频率的循环双谱。这种表示方法在信号特征频率已知的情况下是有效的, 但在信号特征频率不明确的情况下, 不能检测信号特征 ^[6,7]。循环频率区间内计算信号的特征表示方法既能有效抑制信号中的噪声, 又能对信号中的特征具有检测能力, 对于信号x (t) , 这个过程包括两步:

步1 根据式 (5) , 计算信号x (t) 在循环频率α处的三阶循环累积量, 即计算循环频率集α∈{α_j}内所有循环累积量之和。如α未知, 可以计算一段区间内以一定步长为循环频率的累积量之和:

$C_{3 x}^{α} (τ_{1}, τ_{2}) = \sum_{α \in {α_{k}}} C_{x}^{α_{k}} (τ_{1}, τ_{2}) (7)$ $C_{3 x}^{α} (τ_{1}, τ_{2}) = \sum_{α \in {α_{k}}} C_{x}^{α_{k}} (τ_{1}, τ_{2}) (7)$

步2 根据式 (6) 计算C $_{3 x}^{{α}}$ $_{3 x}^{{α}}$ (τ₁, τ₂) 的二维傅立叶变换得到在频率集α∈{α_k}上的循环双谱。

上述方法虽然不能表示循环双谱对循环频率的分布情况, 但能表示循环双谱特征对双频率的分布, 并能在一个三维空间内表示信号在多个循环频率中的特征, 增加了该方法检测信号特征的能力。

《1.2简化算法》

1.2简化算法

为了减少计算, 提高实现的效率, 准备如下:

对于一个零均值的循环平稳过程{x (t) }, 有

$C_{3 x}^{α} (τ_{1}, τ_{2}) = Μ_{3 x}^{α} (τ_{1}, τ_{2}) (8)$ $C_{3 x}^{α} (τ_{1}, τ_{2}) = Μ_{3 x}^{α} (τ_{1}, τ_{2}) (8)$

这样不但可以省去计算一阶和二阶循环矩, 而且还将在后处理中减少大量的复数乘法运算, 直接计算三阶循环矩即可获得三阶循环累积量。

计算三阶循环矩时, 研究三阶滞后积的问题:

$L_{3 x} = x (t) x (t + τ_{1}) x (t + τ_{2}) (9)$ $L_{3 x} = x (t) x (t + τ_{1}) x (t + τ_{2}) (9)$

x (t) 为N点采样数据, 采样频率为f_s, 则t分别为0, 1, …, (N-1) 采样时刻的数据;τ₁, τ₂表示采样的延时, 指延时采样周期的个数。

计算中的数据要重复使用, 故在实现中, 采用保存变量, 避免重复计算。采用如下方法, 既可以减少变量存储, 又可以减少计算量:

L_3x用三维数组来表示, 则为N′×N×N 维, 计t为第一维, 则τ₁, τ₂为第二、第三维, 第一维固定, 第二、第三维组成:

$L^{'}_{k x} = [\begin{matrix} L^{'}_{k x} (1, 1) & L^{'}_{k x} (1, 2) & \dots & L^{'}_{k x} (1, Ν) \\ L^{'}_{k x} (2, 1) & L^{'}_{k x} (2, 2) & \dots & L^{'}_{k x} (2, Ν) \\ ⋮ & ⋮ & ⋮ \\ L^{'}_{k x} (Ν, 1) & L^{'}_{k x} (Ν, 2) & \dots & L^{'}_{k x} (Ν, Ν) \end{matrix}] 。$ $L^{'}_{k x} = [\begin{matrix} L^{'}_{k x} (1, 1) & L^{'}_{k x} (1, 2) & \dots & L^{'}_{k x} (1, Ν) \\ L^{'}_{k x} (2, 1) & L^{'}_{k x} (2, 2) & \dots & L^{'}_{k x} (2, Ν) \\ ⋮ & ⋮ & ⋮ \\ L^{'}_{k x} (Ν, 1) & L^{'}_{k x} (Ν, 2) & \dots & L^{'}_{k x} (Ν, Ν) \end{matrix}] 。$

由式 (9) 知, L′_kx (τ₁, τ₂) 与τ₁, τ₂次序无关, 可见是对称阵, 可以简化成一维的数组, 所以只需要下对角阵的数组元素和 (1+N) N/2个单元即可。原来的纪录数为N², 当N较大时, 即可减少大量的计算及占用的内存。

计算时按照L″_kx (id) =L′_kx (τ₁, τ₂) , 存放时按照τ₁≥τ₂, L″_kx (id) =L′_kx (τ₁, τ₂) , 其中id= (1+τ₁) τ₁/2+ (τ₂-τ₁) 。计算中读取时, 要判断 $i d = {\begin{cases} (1 + τ_{1}) τ_{1} / 2 + (τ_{2} - τ_{1}) τ_{1} \geq τ_{2} \\ (1 + τ_{2}) τ_{2} / 2 + (τ_{1} - τ_{2}) τ_{1} < τ_{2} \end{cases}$ $i d = {\begin{cases} (1 + τ_{1}) τ_{1} / 2 + (τ_{2} - τ_{1}) τ_{1} \geq τ_{2} \\ (1 + τ_{2}) τ_{2} / 2 + (τ_{1} - τ_{2}) τ_{1} < τ_{2} \end{cases}$ , 取出L′_kx (τ₁, τ₂) =L″_kx (id) , 其中id= (1+τ₁) τ₁/2+ (τ₂-τ₁) 。

类似上述方法, L_3x (t, τ₁, τ₂) 存放时按照τ₁≥τ₂, L′_3x (t, id) = L_3x (t, τ₁, τ₂) , 同上, id= (1+τ₁) τ₁/2+ (τ₂-τ₁) 。计算中读取时, 要判断 $i d = {\begin{cases} (1 + τ_{1}) τ_{1} / 2 + (τ_{2} - τ_{1}) τ_{1} \geq τ_{2} \\ (1 + τ_{2}) τ_{2} / 2 + (τ_{1} - τ_{2}) τ_{1} < τ_{2} \end{cases}$ $i d = {\begin{cases} (1 + τ_{1}) τ_{1} / 2 + (τ_{2} - τ_{1}) τ_{1} \geq τ_{2} \\ (1 + τ_{2}) τ_{2} / 2 + (τ_{1} - τ_{2}) τ_{1} < τ_{2} \end{cases}$ , 取出L_3x (t, τ₁, τ₂) =L′_3x (t, id) 。

若L_3x为N′×N×N 维, 则可以减少的复数乘法运算的次数为原来的

$[(1 + Ν) Ν Ν^{'} / 2] / (n^{'} \times Ν \times Ν) = (1 + Ν) / 2 Ν (10)$ $[(1 + Ν) Ν Ν^{'} / 2] / (n^{'} \times Ν \times Ν) = (1 + Ν) / 2 Ν (10)$

然后再计算C_3x (τ₁, τ₂) ;以N=300为例计算三阶循环累积量, 改进后的计算时间大约为原来的6.1/9.5;再运用二维的Fourier变换, 即可完成循环双谱估计。

《2 循环双谱仿真》

2 循环双谱仿真

当齿轮出现裂纹后, 啮合刚度降低, 在传递动力的过程中, 形成有规律的相位变动 (相位调制) ^[8], 忽略幅值变化时的振动信号, 数学表达式为

$x (t) = \sum_{m = 1}^{Μ} x_{m} \cos (2 π m f_{z} t + b_{m} (t) + φ_{m}) (11)$ $x (t) = \sum_{m = 1}^{Μ} x_{m} \cos (2 π m f_{z} t + b_{m} (t) + φ_{m}) (11)$

其中b_m (t) 为相位调制函数:

$b_{m} (t) = \sum_{i = 1}^{Ν} B_{m, j} \cos (2 π i f_{c} t + β_{m, i}) (12)$ $b_{m} (t) = \sum_{i = 1}^{Ν} B_{m, j} \cos (2 π i f_{c} t + β_{m, i}) (12)$

设采集到的信号为:

$\begin{array}{l} y (t) = x (t) + C s i n (2 π f_{a} t) (13) \\ y_{1} (t) = y (t) + a (14) \end{array}$ $\begin{array}{l} y (t) = x (t) + C s i n (2 π f_{a} t) (13) \\ y_{1} (t) = y (t) + a (14) \end{array}$

式 (11) 至式 (14) 中, M=1; N=1; φ_m=0;β_{m, i}=0; f_z=17; f_c=5; f_a=41; C=1; 以及σ (a) =0.125; E{a}=0, a为白噪声。

以调相信号与正弦信号的混合信号为对象进行仿真实验。信号y (t) , y₁ (t) 的循环双谱见图1。

由图1a看到, 有一频率为17 Hz的突出谱线, 还可以得到 (17±5) Hz的频率, 以及41 Hz的频率成分, 这些频率成分还存在频率的相互作用。

由图1b同样可以看到如同图1a中的频率成分及其相互作用, 循环双谱具有一定的降噪功能。

《图1》

图1 仿真信号y, y₁循环双谱图 Fig.1 The cyclic bispectrum of the simulation testing signal (y and y₁)

《3 循环双谱应用》

3 循环双谱应用

齿轮箱振动信号为周期平稳性 ^[9]。以某化工厂的R05齿轮箱为对象, 用笔者提出的方法进行分析。该齿轮箱由5根轴组成, 内部有4对齿轮啮合, 其上装有4个振动加速度传感器, 其中测点3位于Ⅴ轴的齿轮箱壳体垂直方向;由于种种原因, Ⅳ轴与Ⅴ轴及其齿轮对经常发生故障 ^[10]。

驱动电机的转速为1 017 r/min, 经过4级减速, 输出的Ⅴ轴转速约为5 r/min。该测点信号数据的采样长度为1 024点, 采样频率为704 Hz, 采样时做过高通滤波, 滤除3 Hz以下频率成分。图2分别绘出了该采样信号的循环双谱图。该齿轮箱的主要特征频率如表1所示。

《图2》

图2 采样数据的循环双谱分析 Fig.2 The analysis of industrial vibration data by cyclic bispectrum

表1 R05齿轮箱主要特征频率表 Table 1 The characteristic frequency of the gearbox R05 Hz

《表1》

项目	齿数	1×	2×	3×	4×	5×	6×	7×	8×
电机工频		16.95	33.90	50.85	67.80	84.75	101.70	118.65	135.60
电机和Ⅱ轴啮合	15/68	254.25	508.50	726.75	1 017.00	1 525.50	1 779.75	2 034.00
Ⅱ和Ⅲ轴啮合	16/63	59.82	119.65	179.47	239.29	299.12	358.94	418.76	478.59
Ⅲ和Ⅳ轴啮合	15/54	14.24	28.49	42.73	56.97	71.22	85.46	99.71	113.95
Ⅳ和Ⅴ轴啮合	14/48	3.69	7.39	11.08	14.77	18.46	22.16	25.85	29.54

由图2看出该采样的频率成分包含啮合频率3.7 Hz的2×, 3×, 5×, 6×, 7×, 8×以及电机工频19.7 Hz的2倍频;主要频率成分为啮合频率的2×, 3×以及电机工频的2倍频。由此可以看出, 利用循环双谱从复杂信号看出系统的特征频率, 从而能够更好地分析旋转机械设备的运行状态。

《4 结论》

4 结论

针对循环双谱的计算量大, 计算时间长的特点, 提出了循环双谱的简化实现方法, 并给出了该方法在仿真信号和某齿轮箱振动信号上的应用, 得出下列结论:

1) 提出未知循环频率的周期平稳信号循环双谱的估计方法, 即计算一段区间内以一定步长为循环频率的累积量之和来代替循环频率集内所有循环累积量之和。

2) 循环累积量的计算中涉及三阶滞后积问题的简化存储, 滞后积矩阵是对称阵, 通过算法可以用下三角阵的元素来表述, 这样既可以减少存储空间的占用, 又可以减少计算时间。

3) 分析了仿真故障信号的特性, 指出循环双谱对调相信号具有一定的处理能力和降噪能力;通过现场采样数据进行分析, 表明循环双谱对旋转机械的状态分析具有一定的能力。

展示更多