CN103499443B

CN103499443B - A kind of gear distress is without key phase angular domain average computation order analysis method

Info

Publication number: CN103499443B
Application number: CN201310416504.7A
Authority: CN
Inventors: 訾艳阳; 何正嘉; 曹宏瑞; 万志国; 冯超亮; 何水龙
Original assignee: Xian Jiaotong University
Current assignee: Xian Jiaotong University
Priority date: 2013-09-12
Filing date: 2013-09-12
Publication date: 2016-01-20
Anticipated expiration: 2033-09-12
Also published as: CN103499443A

Abstract

The invention provides a method for analyzing the average calculation order of the gear failure keyless phase angle domain. Firstly, the collected vibration acceleration signal is subjected to low-pass phase-preserving filtering, and then the time-frequency distribution is calculated by smoothing pseudo-Wigner-Ville distribution, and then the time-frequency distribution is calculated by Viterbi The optimal path search algorithm estimates the instantaneous frequency of the gearbox shaft, and then uses the key phase signal estimation model to integrate the instantaneous frequency point by point to obtain the estimated key phase signal, and finally combines the equal angle resampling and angle domain averaging technology to calculate the vibration acceleration signal Order analysis to obtain order spectra based on instantaneous frequency estimates. This order spectrum can fully reflect the characteristic information of gearbox faults. The invention combines the advantages of smooth pseudo Wigner-Ville distribution, Viterbi optimal path search algorithm, angle domain average technology and calculation order analysis, and can effectively diagnose the fault of the gearbox under the variable speed operating condition.

Description

An order analysis method for gear failure keyless phase angle domain average calculation

技术领域technical field

本发明属于齿轮故障分析方法领域，具体涉及一种齿轮故障无键相角域平均计算阶次分析方法。The invention belongs to the field of gear fault analysis methods, and in particular relates to a gear fault keyless phase angle field average calculation order analysis method.

背景技术Background technique

传统的齿轮箱振动测试以及对测试信号的分析、测试是在恒定转速、恒定工况下进行的，而实际工作中的齿轮箱往往是变转速、变工况负载下运行，例如风电发电机组的齿轮箱，直升机齿轮箱，汽车减速器，起重机等。以起重机为例，起重机负责吊运钢水包，工作时齿轮箱经历升速提升、降速倾倒的过程。目前，在齿轮箱变转速运行工况条件下，所用的计算阶次分析方法中，转轴的键相信号是进行齿轮箱故障计算阶次分析所必不可少的。然而，对于一些箱体相对封闭(旋转机械不可拆开)或者不可在轴端安装传感器齿轮箱，很难提取到键相信号。特别是有些齿轮减速机构在设计时考虑不周或者装配时结构紧凑，无法安装键相脉冲传感器，很难在变转速工况下实现等周期采样，导致在实际工作过程中对齿轮箱振动测试以及对测试信号的分析、测试较难完成。因此，研究有效的方法，针对变转速工况下的齿轮箱并且转速信息无法采集的情况下监测齿轮箱的运行状态，安排合理的预定备件、停机更换时间，把设备维修安排在设备故障事故发生之前，最大限度的减少突发停机事件，充分发挥设备效能，延长设备使用寿命，对国民生产的安全性、可靠性具有重要意义。The traditional gearbox vibration test and the analysis and test of the test signal are carried out under constant speed and constant working conditions, but the gearboxes in actual work are often operated under variable speed and variable working conditions, such as wind power generators. Gearboxes, helicopter gearboxes, automotive reducers, cranes, etc. Take a crane as an example. The crane is responsible for lifting and transporting the ladle. When working, the gearbox undergoes a process of speed-up, speed-down and dumping. At present, in the calculation order analysis method used under the variable speed operating condition of the gearbox, the key phase signal of the rotating shaft is indispensable for the calculation order analysis of the gearbox fault. However, for some cases that are relatively closed (the rotating machine cannot be disassembled) or the sensor gearbox cannot be installed at the shaft end, it is difficult to extract the key phase signal. In particular, some gear reduction mechanisms are poorly designed or assembled in a compact structure, and cannot be installed with a key phase pulse sensor. It is difficult to achieve equal periodic sampling under variable speed conditions, resulting in the vibration test of the gearbox during actual work. It is difficult to complete the analysis and test of the test signal. Therefore, effective methods should be studied to monitor the operation status of the gearbox under the condition of variable speed and the speed information cannot be collected, arrange reasonable scheduled spare parts, shutdown replacement time, and arrange equipment maintenance when equipment failure accidents occur Previously, it was of great significance to the safety and reliability of national production to minimize sudden downtime events, give full play to equipment performance, and extend equipment service life.

发明内容Contents of the invention

本发明的目的在于提供一种齿轮故障无键相角域平均计算阶次分析方法，该方法能够解决在键相信号不可采集的情况下进行变转速工况下故障诊断的难点问题。The purpose of the present invention is to provide a gear fault keyless phase angle domain average calculation order analysis method, which can solve the difficult problem of fault diagnosis under the condition of variable speed when the key phase signal cannot be collected.

为了实现上述目的，本发明采用的技术方案是：In order to achieve the above object, the technical scheme adopted in the present invention is:

一种齿轮故障无键相角域平均计算阶次分析方法，包括以下步骤：A gear fault keyless phase angle domain average calculation order analysis method, comprising the following steps:

1)瞬时频率估计，具体包括以下三个步骤：1) Instantaneous frequency estimation, specifically including the following three steps:

1-1)采集振动加速度信号，对采集的振动加速度信号进行低通保相滤波，得到反映转子转动的低频信号x(t)；1-1) Collect vibration acceleration signals, perform low-pass phase-preserving filtering on the collected vibration acceleration signals, and obtain low-frequency signals x(t) reflecting rotor rotation;

1-2)通过平滑伪Wigner-Ville分布计算得到的低频信号x(t)的时频分布，得到时频分布图；1-2) The time-frequency distribution of the low-frequency signal x(t) calculated by smoothing the pseudo-Wigner-Ville distribution to obtain a time-frequency distribution diagram;

1-3)在得到的时频分布图上利用Viterbi最优路径搜索算法估计齿轮箱转轴的瞬时频率(InstantaneousFrequency,IF)；。1-3) Use the Viterbi optimal path search algorithm to estimate the instantaneous frequency (Instantaneous Frequency, IF) of the gearbox shaft on the obtained time-frequency distribution map;

2)基于瞬时频率估计的阶次谱，具体包括以下四个步骤：2) Order spectrum based on instantaneous frequency estimation, which specifically includes the following four steps:

2-1)利用键相信号估计模型对步骤1-3)得到的瞬时频率进行逐点积分，得到估计键相信号；2-1) using the key-phase signal estimation model to carry out point-by-point integration of the instantaneous frequency obtained in step 1-3), to obtain an estimated key-phase signal;

2-2)利用得到的估计键相信号对振动加速度信号进行等角度重采样，得到经过等角度重采样后的信号；2-2) Using the obtained estimated key phase signal to carry out equiangular resampling to the vibration acceleration signal, to obtain the signal after equiangular resampling;

2-3)对经过等角度重采样后的信号进行角域平均处理，得到角域平均处理后的信号；2-3) performing angle-domain average processing on the signal after equiangular resampling, to obtain a signal after angle-domain average processing;

2-4)对角域平均处理后的信号进行快速傅里叶变换(FFT)，得到能够反映齿轮箱故障特征信息的基于瞬时频率估计的阶次谱，对齿轮故障进行分析。2-4) Fast Fourier transform (FFT) is performed on the signal after angular domain averaging processing to obtain an order spectrum based on instantaneous frequency estimation that can reflect the characteristic information of the gearbox fault, and analyze the gear fault.

所述的步骤1-2)的具体操作为：The concrete operation of described step 1-2) is:

步骤1-1)得到的低频信号x(t)的平滑伪Wigner-Ville分布SPW(t,ω)由式(1)所示：The smooth pseudo-Wigner-Ville distribution SPW(t,ω) of the low-frequency signal x(t) obtained in step 1-1) is expressed by formula (1):

$SPW SPW ((t t,, ω ω)) = = \frac{11}{22 π π} {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} x x ((t t + + \frac{11}{22} τ τ)) {x x}^{* *} ((t t - - \frac{11}{22} τ τ)) h h ((τ τ)) g g ((u u)) {e e}^{- - iτω iτω} dτdu dτdu - - - - - - ((11))$

由式(1)即得到时频分布图；式(1)中*表示共轭，t表示时间，ω表示频率，h(τ)和g(u)是窗函数，τ,u分别表示时间轴与频率轴的平移因子。The time-frequency distribution diagram can be obtained from formula (1); in formula (1), * means conjugate, t means time, ω means frequency, h(τ) and g(u) are window functions, τ, u represent time axis respectively The translation factor with the frequency axis.

所述的窗函数包括Gauss窗和Hamming窗。The window functions include Gauss window and Hamming window.

所述的步骤1-3)的具体操作为：The concrete operation of described step 1-3) is:

取x(t)的平滑伪Wigner-Ville分布上的时间序列[n₁,n₂]，时间点n∈[n₁,n₂]，每一时间点n上的频率值为SPW_n，把SPW_n按降序排列后的数组为记为SPW′_n；Take the time series [n ₁ ,n ₂ ] on the smooth pseudo-Wigner-Ville distribution of x(t), the time point n∈[n ₁ ,n ₂ ], the frequency value of each time point n is SPW _n , put The array after SPW _n is arranged in descending order is denoted as SPW′ _n ;

按从大到小的顺序从0开始依次按增量1分别重新赋值，赋值后的时频分布序列记为f(n,SPW′_n)；In order from large to small, start from 0 and then re-assign by increment of 1, and the time-frequency distribution sequence after assignment is recorded as f(n, SPW′ _n );

瞬时频率即为使式(2)最小的路径：The instantaneous frequency is the path that minimizes formula (2):

$\overset{^^}{ω ω} ((n no)) = = arg arg min min [[{Σ Σ}_{n no = = {n no}_{11}}^{{n no}_{22}} g g (({ω ω}_{n no},, {ω ω}_{n no + + 11})) + + {Σ Σ}_{n no = = {n no}_{11}}^{{n no}_{22}} f f ((n no,, {SPW SPW}_{n no}^{' '}))]] - - - - - - ((22))$

其中，g(ω_n,ω_n+1)是为了满足频率变化的连续性而对距离较大的频率变化加上的惩罚因子。Among them, g(ω _n ,ω _n+1 ) is a penalty factor added to frequency changes with large distances in order to satisfy the continuity of frequency changes.

所述的惩罚因子g(ω_n,ω_n+1)如式(3)所示：The penalty factor g(ω _n ,ω _n+1 ) is shown in formula (3):

$g g (({ω ω}_{n no},, {ω ω}_{n no + + 11})) = = \{\begin{matrix} 00,, & | | {ω ω}_{n no} - - {ω ω}_{n no + + 11} | | \leq \leq δ δ \\ c c ((| | {ω ω}_{n no} - - {ω ω}_{n no + + 11} | | - - δ δ)),, & | | {ω ω}_{n no} - - {ω ω}_{n no + + 11} | | > > δ δ \end{matrix} - - - - - - ((33))$

其中δ为2个连续点的瞬时频率变化的最大期望值。where δ is the maximum expected value of the instantaneous frequency change of 2 consecutive points.

所述的步骤2-1)的具体操作为：The concrete operation of described step 2-1) is:

一根转轴处于变转速工况下工作，转动瞬时频率为f(t)；经过时间t_n，转轴转过n圈，转轴转过角度如式(4)所示：A rotating shaft works under the condition of variable speed, and the instantaneous frequency of rotation is f(t); after time t _n , the rotating shaft rotates n times, and the angle of rotation of the rotating shaft is shown in formula (4):

${&Integral; &Integral;}_{00}^{{t t}_{n no}} 22 πf πf ((t t)) dt dt = = 22 πn πn - - - - - - ((44))$

对f(t)逐点积分，记录积分值为整数的t_n，在t_n位置插入一个阶跃的键相信号，即是由瞬时频率得到的估计键相信号。Integrate f(t) point by point, record the integral value t _{n as an integer, and insert a step bond-phase signal at the position of t n} _, which is the estimated bond-phase signal obtained from the instantaneous frequency.

所述的步骤2-2)的具体操作为：The concrete operation of described step 2-2) is:

根据式(5)得到等角度θ重采样时的时间节点t，进行等角度重采样，得到经过等角度重采样后的信号；According to the formula (5), the time node t when the equiangular θ resampling is obtained, and the equiangular resampling is performed to obtain the signal after the equiangular resampling;

$t t = = \frac{11}{{22 b b}_{22}} [[\sqrt{{44 b b}_{22} (({θ θ - - b b}_{00})) + + {b b}_{11}^{22}} - - {b b}_{11}]] - - - - - - ((55))$

式(5)中b₀、b₁、b₂根据式(6)得到；In formula (5), b ₀ , b ₁ , b ₂ can be obtained according to formula (6);

其中 in

所述的步骤2-3)的具体操作为：The concrete operation of described step 2-3) is:

将经过等角度重采样后的信号序列记为由角域平均处理的公式式(7)得到角域平均处理后的信号序列y_n，Denote the signal sequence after equiangular resampling as The signal sequence y _n after the angle domain average processing is obtained by the formula (7) of the angle domain average processing,

${y the y}_{n no} = = \frac{11}{P P} {Σ Σ}_{P P = = 00}^{P P - - 11} {\overset{~ ~}{y the y}}_{n no + + PN PN} - - - - - - ((77))$

其中，P为将经过等角度重采样后的信号平均分成的段数，N为各段相同的采样点数。Wherein, P is the number of segments that the equiangularly resampled signal is divided into on average, and N is the same number of sampling points in each segment.

相对于现有技术，本发明的有益效果为：Compared with the prior art, the beneficial effects of the present invention are:

本专利提出的齿轮故障无键相角域平均计算阶次分析方法方法，能够从振动加速度信号中提取转轴的转速信息，有效解决无键相信号(键相信号不可采集)时齿轮箱在变转速运行工况下故障难以诊断的问题，即能够在无键相信号时及变转速工况下完成对齿轮箱故障的诊断。该方法综合了平滑伪Wigner-Ville分布、Viterbi最优路径搜索算法、键相信号估计、角域平均处理以及计算阶次分析的优点，具有下列区别于传统方法的显著优势：The gear failure keyless phase angle domain average calculation order analysis method proposed in this patent can extract the rotational speed information of the rotating shaft from the vibration acceleration signal, and effectively solve the problem of the gear box changing the speed when the keyless phase signal (the key phase signal cannot be collected) The problem that the fault is difficult to diagnose under the operating condition, that is, the fault diagnosis of the gearbox can be completed when there is no key phase signal and under the condition of variable speed. This method combines the advantages of smooth pseudo-Wigner-Ville distribution, Viterbi optimal path search algorithm, bond phase signal estimation, angle domain average processing and calculation order analysis, and has the following significant advantages different from traditional methods:

1)在瞬时频率估计方面，本发明结合了平滑伪Wigner-Ville分布良好的时频聚焦性，以及Viterbi最优路径搜索算法遵循自然频率变化连续性原则的特点，能够有效准确的估计齿轮箱转轴的瞬时频率。1) In terms of instantaneous frequency estimation, the present invention combines the good time-frequency focus of the smooth pseudo-Wigner-Ville distribution, and the characteristics of the Viterbi optimal path search algorithm following the principle of continuity of natural frequency changes, which can effectively and accurately estimate the gearbox shaft the instantaneous frequency.

2)在基于瞬时频率估计的阶次谱计算方面，利用键相信号估计模型对瞬时频率进行逐点积分得到精确的估计键相信号，然后结合等角度重采样技术及角域平均处理技术对振动加速度信号进行计算阶次分析，有效地抵御了噪声及非相关信号的干扰，得到能够反映齿轮箱故障特征信息的基于瞬时频率估计的阶次谱。2) In terms of order spectrum calculation based on instantaneous frequency estimation, the instantaneous frequency is integrated point-by-point using the bond-phase signal estimation model to obtain an accurate estimate of the bond-phase signal, and then the vibration The calculation order analysis of the acceleration signal can effectively resist the interference of noise and non-correlated signals, and obtain the order spectrum based on instantaneous frequency estimation that can reflect the fault characteristic information of the gearbox.

3)整体上，本发明提供的方法综合了平滑伪Wigner-Ville分布、Viterbi最优路径搜索算法、键相信号估计、角域平均处理技术以及计算阶次分析的优点，利用本发明可以有效地对变转速工况下无法采集键相信号的齿轮箱进行故障诊断，提高故障诊断的准确率。3) On the whole, the method provided by the present invention combines the advantages of smooth pseudo Wigner-Ville distribution, Viterbi optimal path search algorithm, key phase signal estimation, angle domain average processing technology and calculation order analysis, and the present invention can effectively Carry out fault diagnosis for gearboxes that cannot collect key phase signals under variable speed conditions to improve the accuracy of fault diagnosis.

附图说明Description of drawings

图1(a)是Viterbi最优路径搜索算法前向搜索的示意图；Figure 1(a) is a schematic diagram of the forward search of the Viterbi optimal path search algorithm;

图1(b)是Viterbi最优路径搜索算法后向搜索的示意图；Figure 1(b) is a schematic diagram of the backward search of the Viterbi optimal path search algorithm;

图1(c)是Viterbi最优路径搜索算法的局部最优路径的示意图；Fig. 1 (c) is the schematic diagram of the locally optimal path of Viterbi optimal path search algorithm;

图1(d)是Viterbi最优路径搜索算法的全局最优路径的示意图；Figure 1(d) is a schematic diagram of the global optimal path of the Viterbi optimal path search algorithm;

图2是轴的转动模型的示意图；Fig. 2 is the schematic diagram of the rotational model of shaft;

图3(a)是仿真振动加速信号图；Fig. 3 (a) is the simulated vibration acceleration signal diagram;

图3(b)是对仿真振动加速信号的插值重采样图；Fig. 3 (b) is the interpolation resampling diagram to the simulated vibration acceleration signal;

图3(c)是仿真振动加速信号的估计键相信号图；Fig. 3 (c) is the estimation bond phase signal diagram of simulation vibration acceleration signal;

图3(d)是对估计键相信号进行等角度重采样后的信号图；Figure 3(d) is a signal diagram after equiangular resampling of the estimated bond phase signal;

图4是实验测得的振动加速度信号图；Fig. 4 is the vibration acceleration signal diagram that the experiment measures;

图5是对所得振动加速信号低通滤波后的信号图；Fig. 5 is the signal figure after the low-pass filtering of gained vibration acceleration signal;

图6是低通滤波后信号的SPWVD时频分布图，其中a为1X转频，b为2X转频；Figure 6 is the SPWVD time-frequency distribution diagram of the low-pass filtered signal, where a is 1X frequency conversion, and b is 2X frequency conversion;

图7是利用Viterbi最优路径搜索算法搜索得到的瞬时频率图；Fig. 7 is the instantaneous frequency figure that utilizes Viterbi optimal path search algorithm to search and obtain;

图8是利用键相信号估计模型得到的估计键相信号图；Fig. 8 is the estimated key-phase signal diagram obtained by using the key-phase signal estimation model;

图9是等角度重采样后的信号图；Fig. 9 is a signal diagram after equiangular resampling;

图10是基于瞬时频率估计的阶次谱图；Fig. 10 is the order spectrogram based on instantaneous frequency estimation;

图11是行星齿轮箱中心轮齿面均匀点蚀图。Figure 11 is a uniform pitting diagram of the tooth surface of the center gear of the planetary gearbox.

具体实施方式detailed description

下面结合附图对本发明的内容作进一步详细说明。The content of the present invention will be described in further detail below in conjunction with the accompanying drawings.

本发明公开了一种齿轮故障无键相角域平均计算阶次分析方法(Angle-domain-averageAndNon-keyphasorCalculatedOrderTracking,AN-COT)，具体按以下步骤实施：The invention discloses a gear fault keyless phase angle domain average calculation order analysis method (Angle-domain-averageAndNon-keyphasor Calculated Order Tracking, AN-COT), which is specifically implemented according to the following steps:

1)瞬时频率估计。对采集的振动加速度信号进行低通保相滤波，得到反映转子转动的低频信号。然后通过平滑伪Wigner-Ville分布计算低频信号的时频分布，在得到的频分布图上利用Viterbi最优路径搜索算法估计齿轮箱转轴的瞬时频率(InstantaneousFrequency,IF)。具体包括以下3个步骤：1) Instantaneous frequency estimation. A low-pass phase-preserving filter is performed on the collected vibration acceleration signal to obtain a low-frequency signal reflecting the rotor rotation. Then the time-frequency distribution of the low-frequency signal is calculated by smoothing the pseudo-Wigner-Ville distribution, and the instantaneous frequency (Instantaneous Frequency, IF) of the gearbox shaft is estimated by using the Viterbi optimal path search algorithm on the obtained frequency distribution map. Specifically, it includes the following three steps:

1-1)采集振动加速度信号，对采集的振动加速度信号进行低通保相滤波，提取反映旋转轴对信号首先进行包含输入轴转频信号在内的低通保相位滤波，去除信号中其他阶次的干扰，得到以输入轴振动信号及其低倍阶次为主要分量的振动信号(反映转子转动的低频信号)。低通滤波器的上限截止频率根据机器运转过程中的最高运转速度设定。同时保相位滤波十分重要，否则滤波后如果有相位的错动会引起IF估计与振动信号在时间轴上的不同步。1-1) Collect vibration acceleration signals, perform low-pass phase-preserving filtering on the collected vibration acceleration signals, and extract and reflect the rotating shaft pair signal. First, perform low-pass phase-preserving filtering including the input shaft frequency signal to remove other order The vibration signal (low frequency signal reflecting the rotation of the rotor) is obtained with the input shaft vibration signal and its low order as the main component. The upper limit cut-off frequency of the low-pass filter is set according to the highest operating speed during machine operation. At the same time, phase-preserving filtering is very important, otherwise, if there is a phase shift after filtering, it will cause the IF estimation and the vibration signal to be out of sync on the time axis.

1-2)时频分析：通过平滑伪Wigner-Ville分布计算得到的低频信号的时频分布，得到时频分布图。1-2) Time-frequency analysis: the time-frequency distribution of the low-frequency signal calculated by smoothing the pseudo-Wigner-Ville distribution to obtain the time-frequency distribution diagram.

低频信号x(t)的Wigner-Ville分布定义为，The Wigner-Ville distribution of a low-frequency signal x(t) is defined as,

$WVD WVD ((t t,, ω ω)) = = \frac{11}{22 π π} {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} {x x}^{* *} ((t t - - τ τ / / 22)) x x ((t t + + τ τ / / 22)) {e e}^{- - iτω iτω} dτ dτ$

Wigner-Ville具有众多的优点，特别是时频聚集性，但会引入交叉项。为了克服交叉项的不足，可以在时域或时域、频域同时加窗以克服交叉项，这就是伪Winger分布(PWVD)和平滑伪Winger分布(SPWVD)。在时域加窗后的Wigner-Ville分布称为伪Wigner-Ville分布，定义如下：Wigner-Ville has many advantages, especially time-frequency aggregation, but it will introduce cross terms. In order to overcome the deficiency of the cross term, windowing can be performed in the time domain or the time domain and the frequency domain simultaneously to overcome the cross term, which is pseudo-Winger distribution (PWVD) and smooth pseudo-Winger distribution (SPWVD). The Wigner-Ville distribution after windowing in the time domain is called the pseudo-Wigner-Ville distribution, which is defined as follows:

$PW PW ((t t,, ω ω)) = = \frac{11}{22 π π} {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} h h ((τ τ)) {x x}^{* *} ((t t - - \frac{11}{22} τ τ)) x x ((t t + + \frac{11}{22} τ τ)) {e e}^{- - iτω iτω} dτ dτ$

其中，公式中*表示共轭。h(t)是窗函数。常用的窗函数包括Gauss窗、Hamming窗等。Wherein, * in the formula means conjugation. h(t) is the window function. Commonly used window functions include Gauss window, Hamming window, etc.

在时域、频域同时加窗则可以在时间和频率两个方向同时抑制交叉项。这种两个方向加窗处理的Wigner-Ville分布称为平滑伪Wigner-Ville分布(SPWVD)，x(t)的平滑伪Wigner-Ville分布定义如下：Simultaneous windowing in the time domain and frequency domain can suppress cross-terms in both time and frequency directions. The Wigner-Ville distribution of windowing in two directions is called smooth pseudo-Wigner-Ville distribution (SPWVD), and the smooth pseudo-Wigner-Ville distribution of x(t) is defined as follows:

其中，t表示时间，ω表示频率，g(u)是用来在频率轴方向做平滑的窗函数，τ,u分别表示时间轴与频率轴的平移因子。由式(1)即得到时频分布图。Among them, t represents time, ω represents frequency, g(u) is a window function used for smoothing in the direction of the frequency axis, τ, u represent the translation factors of the time axis and the frequency axis, respectively. The time-frequency distribution diagram can be obtained by formula (1).

1-3)在得到的时频分布图上利用Viterbi最优路径搜索算法估计齿轮箱转轴的瞬时频率(InstantaneousFrequency,IF)。1-3) Use the Viterbi optimal path search algorithm to estimate the instantaneous frequency (Instantaneous Frequency, IF) of the gearbox shaft on the obtained time-frequency distribution map.

在1-2)步求得时频分布图上，采用Viterbi算法搜索谱峰最优时频路径。谱峰检测只是依次搜索时频分布上的峰值，实际的IF作为一个连续变化的物理量，应该考虑两个相邻时间点之间的变化不能过大，Viterbi算法正好考虑了这点。在时频分布上搜索IF时提出如下两个假设，On the time-frequency distribution map obtained in step 1-2), the Viterbi algorithm is used to search for the optimal time-frequency path of the spectral peak. Spectrum peak detection only searches for peaks in the time-frequency distribution in turn. As the actual IF is a continuously changing physical quantity, it should be considered that the change between two adjacent time points should not be too large. The Viterbi algorithm just takes this into account. The following two assumptions are made when searching for IF on the time-frequency distribution,

(1)若在某个时间点上，时频分布的最大值不是信号的IF，则信号的IF值是SPWVD第n(n较小)个最大值的概率很大；(1) If at a certain time point, the maximum value of the time-frequency distribution is not the IF of the signal, then the probability that the IF value of the signal is the nth (smaller n) maximum value of SPWVD is very high;

(2)2个连续时间点的IF变化不是非常大。(2) The IF changes at two consecutive time points are not very large.

仅看第一条假设，实质上是谱峰检测法。Viterbi算法的关键之处于是第二条假设，此假设反映的是旋转机械变工况运行工程中自然频率变化的连续性，不可能在极短的时间内出现很大的频率跳跃。Just looking at the first hypothesis, it is essentially a peak detection method. The key point of the Viterbi algorithm is the second assumption. This assumption reflects the continuity of natural frequency changes in rotating machinery operating under variable conditions, and it is impossible to have a large frequency jump in a very short time.

Viterbi算法归结于一个极值问题，即取SPWVD分布上的时间序列[n₁,n₂]之间的时频分布来分析，时间点n∈[n₁,n₂]，每一时间点n上的频率值为SPW_n，即The Viterbi algorithm is attributed to an extreme value problem, which is to analyze the time-frequency distribution between the time series [n ₁ , n ₂ ] on the SPWVD distribution. The time point n∈[n ₁ , n ₂ ], each time point n The frequency value on SPW _n , that is

SPW_n＝[SPW(n，ω₁)，SPW(n，ω₂)，…，SPW(n，ω_M)]SPW _n = [SPW(n, ω ₁ ), SPW(n, ω ₂ ), . . . , SPW(n, ω _M )]

M为包含的频率点数。M is the number of frequency points included.

把SPW_n按降序排列后的数组定义为SPW′_n，则f(n,SPW′_n)定义为赋值后的时频分布序列，The array of SPW _n arranged in descending order is defined as SPW′ _n , then f(n,SPW′ _n ) is defined as the time-frequency distribution sequence after assignment,

f(n,SPW′_n)＝j-1(j＝0,1,…,M)f(n,SPW′ _n )=j-1(j=0,1,…,M)

上式的意义就是把SPWVD分布时间点n上的所有频率点处的值按从大到小的顺序从0开始依次赋值按增量1分别赋值。The meaning of the above formula is to assign values at all frequency points at the SPWVD distribution time point n in order from large to small, starting from 0 and assigning values in increments of 1 respectively.

瞬时频率即为使式(2)最小的路径，The instantaneous frequency is the path that minimizes formula (2),

式中，f(n,SPW′_n)是根据假设(1)进行构造，g(ω_n,ω_n+1)根据假设(2)构造。g(ω_n,ω_n+1)是为了满足频率变化的连续性而对距离较大的频率变化加上的惩罚因子，可以保证IF路径满足以上假设(2)。g(ω_n,ω_n+1)的定义如下：In the formula, f(n,SPW′ _n ) is constructed according to assumption (1), and g(ω _n ,ω _n+1 ) is constructed according to assumption (2). g(ω _n ,ω _n+1 ) is a penalty factor added to the frequency change with a large distance in order to satisfy the continuity of the frequency change, which can ensure that the IF path satisfies the above assumption (2). g(ω _n ,ω _n+1 ) is defined as follows:

式中，δ的最优选择为2个连续点IF变化的最大期望值，即对于较小的IF变化(|x-y|≤δ)的代价函数为0,在实验中选取较小的δ(如δ＝3)可以得到很好的效果。In the formula, the optimal choice of δ is the maximum expected value of the IF change at two consecutive points, that is, the cost function for a small IF change (|x-y|≤δ) is 0, and a smaller δ is selected in the experiment (such as δ =3) can get good results.

以δ＝1，c＝5为例说明Viterbi最优路径搜索算法的实现过程，取某局部时频面进行分析，时频面上有3个时间点，8个频率点，Viterbi最优路径搜索算法的具体步骤为：Taking δ=1, c=5 as an example to illustrate the implementation process of the Viterbi optimal path search algorithm, take a local time-frequency plane for analysis, there are 3 time points and 8 frequency points on the time-frequency plane, Viterbi optimal path search The specific steps of the algorithm are:

(1)根据式f(n,SPW′_n)＝j-1(j＝0,1,…,M)计算每个点(n,ω_i)对应的f(n,ω_i)函数值，确定每个时间点n的所有频率与其前一时间点n-1的所有频率的局部最优路径并记录，选取f(n,SPW′_n)+g(ω_n,ω_n-1)最小的路径作为两时间点之间的局部最优路径(向左箭头线段)，如图1(a)所示；(1) Calculate the f(n,ω _i ) function value corresponding to each point (n,ω _i ) according to the formula f(n,SPW′ _n )=j-1(j=0,1,…,M), Determine and record the local optimal paths of all frequencies at each time point n and all frequencies at the previous time point n-1, and select the minimum f(n,SPW′ _n )+g(ω _n ,ω _n-1 ) The path is used as a local optimal path between two time points (leftward arrow line segment), as shown in Figure 1(a);

(2)同理，确定每个时间n的所有频率点与其后一时间点n+1的所有频率点的局部最优路径并记录(向右箭头线段)，如图1(b)所示；(2) In the same way, determine the local optimal path of all frequency points at each time n and all frequency points at the next time point n+1 and record (right arrow segment), as shown in Figure 1(b);

(3)根据式(3)，以两个时间点之间的局部最优路径为基准，首先沿着由步骤(1)向前搜索确定的局部最优路径计算所有路径的函数值，记录函数值最小的路径如图1(c)中向左的箭头所示；然后按照同样的方法确定向后搜索的局部最优路径，如图1(c)中向右的箭头所示；(3) According to formula (3), based on the local optimal path between two time points, first calculate the function values of all paths along the local optimal path determined by the forward search in step (1), and record the function The path with the smallest value is shown by the left arrow in Figure 1(c); then follow the same method to determine the local optimal path for backward search, as shown by the right arrow in Figure 1(c);

(4)比较取图1(c)中向前、向后搜索路径函数值的大小，选择最小函数值对应的路径为最优路径，如图1(d)所示。(4) Compare the function values of the forward and backward search paths in Figure 1(c), and select the path corresponding to the minimum function value as the optimal path, as shown in Figure 1(d).

此外，由于时频分布的边缘效应，IF估计的边缘值存在很大误差，可以去掉时频分布图上左右两边的一些数据。同时针对实际工况中齿轮箱振动信号中多分量成分的特点，需要在上述Viterbi算法搜索出一个IF的基础上，估计出其他的IF。假设现已估计出是对应信号最强的信号。将范围内的值置零，形成新的时频表示，ξ为一设定的常量。对新的时频表示重复以上Viterbi搜索算法步骤，直至将信号中各分量的IF估计都估计出来。In addition, due to the edge effect of the time-frequency distribution, there is a large error in the edge value of the IF estimate, and some data on the left and right sides of the time-frequency distribution graph can be removed. At the same time, in view of the characteristics of multi-component components in the vibration signal of the gearbox in actual working conditions, it is necessary to estimate other IFs on the basis of one IF searched by the Viterbi algorithm above. Assumptions have been estimated is the signal corresponding to the strongest signal. Will The values in the range are set to zero to form a new time-frequency representation, and ξ is a set constant. Repeat the above Viterbi search algorithm steps for the new time-frequency representation until the IF estimates of each component in the signal are estimated.

2)基于瞬时频率估计的阶次谱。利用键相信号估计模型对瞬时频率进行逐点积分得到精确的估计键相信号。然后结合角域平均技术对振动加速度信号进行计算阶次分析，得到基于瞬时频率估计的阶次谱。具体包括以下四个步骤：2) Order spectrum based on instantaneous frequency estimation. Using the key-phase signal estimation model, the instantaneous frequency is integrated point by point to obtain an accurate estimated key-phase signal. Then, combined with the angular domain averaging technique, the calculation order analysis of the vibration acceleration signal is carried out, and the order spectrum based on the instantaneous frequency estimation is obtained. Specifically, it includes the following four steps:

2-1)利用键相信号估计模型对步骤1-3)得到的瞬时频率进行逐点积分，得到估计键相信号。在得到IF的基础上，为了应用计算阶次分析，需要用IF估计键相信号。假设一根转轴处于变转速工况下工作，转速为ω(t)，如图2所示。ω(t)和转动瞬时频率f(t)的关系为ω(t)＝2πf(t)。取图2中某时间段转过的角度为Δθ,Δθ可表示为：2-1) Using the bond-phase signal estimation model to perform point-by-point integration on the instantaneous frequency obtained in step 1-3), to obtain an estimated bond-phase signal. On the basis of obtaining the IF, in order to apply the computational order analysis, it is necessary to use the IF to estimate the bond-phase signal. Assume that a rotating shaft works under the condition of variable speed, and the speed is ω(t), as shown in Figure 2. The relationship between ω(t) and the rotational instantaneous frequency f(t) is ω(t)=2πf(t). Take the angle turned in a certain period of time in Figure 2 as Δθ, and Δθ can be expressed as:

Δθ＝ω(t)dt＝2πf(t)dtΔθ=ω(t)dt=2πf(t)dt

现假定从a点作为时间的测量起点t，每一小段Δθ累加，则经过时间t₁，转子又转到a点，转子转过的角度为2π。由此得到如下方程：Assuming that point a is used as the starting point t of time measurement, and each small segment of Δθ is accumulated, then after time t ₁ , the rotor turns to point a again, and the angle the rotor rotates is 2π. This leads to the following equation:

${&Integral; &Integral;}_{00}^{{t t}_{11}} 22 πf πf ((t t)) dt dt = = 22 π π$

如果经过t_n，转轴转过n圈，转轴转过角度如式(4)所示：If after t _n , the shaft rotates n times, the rotation angle is shown in formula (4):

式(4)中，t_n是需要求解的未知量。对f(t)逐点积分，记录积分值为整数的t_n，在t_n位置插入一个阶跃的键相信号，即是由IF得到的估计键相信号。In formula (4), t _n is the unknown quantity to be solved. Integrate f(t) point by point, record the integral value t _{n as an integer, and insert a step bond-phase signal at the position of t n} _, which is the estimated bond-phase signal obtained by IF.

2-2)利用得到的估计键相信号对振动加速度信号进行等角度重采样，得到经过等角度重采样后的信号。2-2) Using the obtained estimated key phase signal to perform equiangular resampling on the vibration acceleration signal to obtain the equiangularly resampled signal.

利用转轴上的键槽标记或者反光片每一转内只能采集到一个脉冲，等角度重采样需要做的工作就是利用采集的一系列的标记每一转2π弧度的脉冲，利用相应的计算方法插值分成N段，每一段的弧度为为计算等角度重采样时刻，假设转轴在每一个小的时间段内转轴是作匀角加速度运行的。Using the keyway marks on the rotating shaft or the reflectors can only collect one pulse per revolution, the work that needs to be done for equiangular resampling is to use a series of collected pulses that mark 2π radians per revolution, and use the corresponding calculation method to interpolate Divided into N segments, the arc of each segment is In order to calculate the equiangular resampling time, it is assumed that the rotating shaft runs at a uniform angular acceleration in each small time period.

转轴转角θ的匀角速度转动方程可以描述为：The uniform angular velocity rotation equation of the shaft rotation angle θ can be described as:

θ(t)＝b₀+b₁t+b₂t² θ(t)=b ₀ +b ₁ t+b ₂ t ²

式中：b₀、b₁、b₂为未知参数，t是时间点。In the formula: b ₀ , b ₁ , b ₂ are unknown parameters, and t is the time point.

取转速脉冲传感器采集的三个连续脉冲的时间点t₁、t₂、t₃，以t₁作为角度增量的起始刻度，则t₁和t₂、t₁和t₃时间点之间的角度增量分别为可列出如下的方程：Take time points t ₁ , t ₂ , t ₃ of the three continuous pulses collected by the speed pulse sensor, and take t ₁ as the initial scale of the angle increment, then the time points between t ₁ and t ₂ , t ₁ and t ₃ The angle increments of The following equations can be listed:

由方程可解出b₀、b₁、b₂，如式(6)所示：b ₀ , b ₁ , b ₂ can be solved from the equation, as shown in formula (6):

同时，由方程θ(t)＝b₀+b₁t+b₂t²解出t：At the same time, t is solved from the equation θ(t)=b ₀ +b ₁ t+b ₂ t ² :

将式(6)代入式(5)即可求出等角度θ重采样时的时间节点t，然后进行等角度重采样，得到经过等角度重采样后的信号。图3所示给出了一个等角度重采样的示意图。在图中，在每两个采集到的键相信号中第每隔个插入一个重采样点(即每转采样8个点)，这些插入的时间点从时间轴上来看，不是等时间间隔的，但角域上看相隔同等间隔的 Substituting Equation (6) into Equation (5) can find the time node t when resampling at equal angle θ, and then perform equal angle resampling to obtain the signal after equal angle resampling. Figure 3 shows a schematic diagram of equiangular resampling. In the figure, in every two key phase signals collected every other Each inserts a resampling point (that is, sampling 8 points per revolution). These inserted time points are not at equal time intervals from the time axis, but at the same interval in the angular domain.

2-3)对经过等角度重采样后的信号进行角域平均处理，得到角域平均处理后的信号。2-3) Perform angle-domain average processing on the equiangularly resampled signal to obtain a signal after angle-domain average processing.

将经过等角度重采样后的信号序列定义为与时域同步平均思想相同，角域同步平均的算法公式为式(7)，由式(7)得到角域平均处理后的信号序列y_n。The signal sequence after equiangular resampling is defined as The same idea as the time domain synchronous averaging, the algorithm formula of the angular domain synchronous averaging is formula (7), and the signal sequence y _n after the angular domain averaging processing is obtained from the formula (7).

其中，in,

P：将经过等角度重采样后的信号平均分成的段数；P: The number of segments that divide the signal after equiangular resampling into average;

N：各段相同的采样点数。N: The same number of sampling points in each segment.

对于旋转机械或者往复机械运行过程中由角域重采样得到的有周期性的机械信号序列y(n)，假设y(n)由周期信号f(n)和随机噪声s(n)组成，即：For the periodic mechanical signal sequence y(n) obtained by angular domain resampling during the operation of rotating machinery or reciprocating machinery, it is assumed that y(n) is composed of periodic signal f(n) and random noise s(n), that is :

x(n)＝f(n)+s(n)x(n)=f(n)+s(n)

以f(n)的周期T去对信号x(n)进行截取，一共截得M段，然后将其相加，由于随机噪声的不相关性，可得到：The signal x(n) is intercepted with the period T of f(n), and a total of M segments are intercepted, and then added together, due to the uncorrelation of random noise, it can be obtained:

$x x (({n no}_{i i})) = = Mf mf (({n no}_{i i})) + + \sqrt{M m} s the s (({n no}_{i i}))$

再对x(n_i)求平均，得到平均后的的输出信号y(n_i)：Then average x(n _i ) to get the averaged output signal y(n _i ):

$y the y (({n no}_{i i})) = = f f (({n no}_{i i})) + + \frac{11}{\sqrt{M m}} s the s (({n no}_{i i}))$

此时随机噪声是原输入信号x(t)中的随机噪声的因而信噪比提高了倍。At this time, the random noise is the random noise in the original input signal x(t) Therefore, the signal-to-noise ratio is improved times.

因此，角域同步采样的滤波器同时域同步采样的滤波器特性相似，能够有效地消除噪声和非相关信号。Therefore, the filters with synchronous sampling in the angle domain have similar characteristics to the filters with synchronous sampling in the angle domain, which can effectively eliminate noise and non-correlated signals.

2-4)对角域平均处理后的信号进行快速傅里叶变换(FFT)，即得到能够反映齿轮箱故障特征信息的基于瞬时频率估计的阶次谱。2-4) Fast Fourier transform (FFT) is performed on the signal processed by the angle-domain average, that is, an order spectrum based on instantaneous frequency estimation that can reflect the fault characteristic information of the gearbox is obtained.

实施案例Implementation case

用某行星减速器试验台所测振动信号验证该发明的有效性。该试验台包括电动机，减速器，发电机，控制台四部分。行星减速器是一个二级精密减速器，减速比为32，行星减速器的结构参数如下表1所示。The effectiveness of the invention is verified by the vibration signal measured by a planetary reducer test bench. The test bench includes four parts: motor, reducer, generator and console. The planetary reducer is a two-stage precision reducer with a reduction ratio of 32. The structural parameters of the planetary reducer are shown in Table 1 below.

表1星减速器结构参数Table 1 Structural parameters of star reducer

根据表1的参数以及行星轮啮合频率计算公式计算第一级传动太阳轮和行星轮的啮合频率。行星轮啮合频率计算公式如下：Calculate the meshing frequency of the first stage transmission sun gear and planetary gear according to the parameters in Table 1 and the formula for calculating the meshing frequency of the planetary gear. The formula for calculating the meshing frequency of the planetary gear is as follows:

f_啮合＝(n₁-n_H)×z₁ f _meshing =(n ₁ -n _H )×z ₁

式中，n₁：太阳轮转速；n_H：行星架转速；z₁：太阳轮齿数。In the formula, n ₁ : speed of sun gear; n _H : speed of planet carrier; z ₁ : number of teeth of sun gear.

由上式得到f_啮合＝10.5n₁。10.5是一个重要的参数，对照前文所述的阶次概念，J_n＝10.5是行星轮出现故障时的阶次参数。From the above formula, f _meshing = 10.5n ₁ . 10.5 is an important parameter. Compared with the order concept mentioned above, J _n =10.5 is the order parameter when the planetary gear fails.

手动调节控制箱，使输入电机的速度为一个减速过程，在时域内分别采集振动信号如图4所示。采样频率为10240Hz，满足采样频率设置要求。Manually adjust the control box so that the speed of the input motor is a deceleration process, and the vibration signals are collected in the time domain as shown in Figure 4. The sampling frequency is 10240Hz, meeting the sampling frequency setting requirements.

根据特征阶次设定最大分析阶次为40，即在轴每转一周的时间内等角度地采样80个点,等角度采样的采样角度为2π/80。According to the characteristic order, the maximum analysis order is set to 40, that is, 80 points are sampled equiangularly during each revolution of the shaft, and the sampling angle of equiangular sampling is 2π/80.

对信号进行截止频率为50Hz的零相滤波，滤波后的信号如图5所示。The zero-phase filter with a cutoff frequency of 50 Hz is performed on the signal, and the filtered signal is shown in Figure 5.

利用SPWVD求取信号的时频分布。取窗函数为Hamming窗，窗口长度为信号长度的1/10，得到时频分布后，再在时频分布两端去掉1s时长的时频分布，以减小端点效应。得到的时频分布如下图6所示。可以看到时频图上明显有输入轴转频的1倍，2倍的时频分布。Use SPWVD to obtain the time-frequency distribution of the signal. The window function is taken as the Hamming window, and the window length is 1/10 of the signal length. After the time-frequency distribution is obtained, the time-frequency distribution with a duration of 1 second is removed at both ends of the time-frequency distribution to reduce the endpoint effect. The resulting time-frequency distribution is shown in Figure 6 below. It can be seen that on the time-frequency diagram, there is obviously a time-frequency distribution of 1 times and 2 times the rotational frequency of the input shaft.

利用Viterbi算法，在时频图上分别搜索出1倍转频的瞬时转频估计如图7中粗实线所示，2倍转频的瞬时转频估计(提取的2倍转频的IF的最优路径除以2，则得到1倍转频的IF估计)如图7虚线所示，真实瞬时频率曲线如图7细实线所示。从图中可以看出瞬时转频估计出的结果与真实瞬时频率曲线能够非常好的吻合。Using the Viterbi algorithm, the instantaneous frequency conversion estimation of 1 times the frequency conversion is searched on the time-frequency diagram respectively, as shown in the thick solid line in Fig. 7, the instantaneous frequency conversion estimation of 2 times the frequency conversion (the extracted IF Divide the optimal path by 2 to obtain the IF estimate of 1 times the conversion frequency) as shown in the dotted line in Figure 7, and the real instantaneous frequency curve is shown in the thin solid line in Figure 7. It can be seen from the figure that the estimated result of the instantaneous conversion frequency is in good agreement with the real instantaneous frequency curve.

首先用1倍转频IF估计给合键相估计算法，即可以得到估计的键相信号，如图8所示。Firstly, the 1-fold frequency IF is used to estimate the bond-phase estimation algorithm, that is, the estimated bond-phase signal can be obtained, as shown in Figure 8.

利用估计的键相信号，对原信号进行插值重采样可以得到等角度的重采样信号，如图9所示。Using the estimated key-phase signal, the original signal can be interpolated and resampled to obtain an equiangular resampled signal, as shown in Fig. 9 .

对角域重采样后的信号进行角域平均处理，并进行FFT变换。，得到如图10所示的阶次谱。在谱图上可以看到明显的与行星齿轮箱第一级传动中的阶次10.5有关的谱峰，如图J_n标记所示，同时两旁伴有明显的基频特征边频带，即边频带为转速阶次1。在2倍啮合频率2J_n处，间隔为基频的变频带高而窄，属于典型的齿轮均匀点蚀故障特征。打开此行星齿轮箱，发现第一级传动中心轮的每个齿面上都有图11所示的齿面均匀点蚀现象(椭圆圈中)。Angle-domain average processing is performed on the resampled signal in the angle domain, and FFT transformation is carried out. , and the order spectrum shown in Figure 10 is obtained. In the spectrogram, it can be seen that there are obvious spectral peaks related to the order 10.5 in the first-stage transmission of the planetary gearbox, as shown in Figure J _n , and there are obvious fundamental frequency characteristic sidebands on both sides, that is, sidebands It is speed order 1. At 2 times the meshing frequency 2J _n , the frequency conversion band whose interval is the fundamental frequency is high and narrow, which is a typical characteristic of uniform pitting faults of gears. Open this planetary gearbox, and find that each tooth surface of the first-stage transmission center wheel has uniform pitting on the tooth surface shown in Figure 11 (in the elliptical circle).

本专利提出的AN-COT算法，在无健相信号下，准确诊断了齿轮箱在变转速工况下的条件下的齿面均匀点蚀故障特征，提高了故障诊断的准确性，为变转速条件下的齿轮箱故障诊断提供了有效可靠的保证。The AN-COT algorithm proposed in this patent accurately diagnoses the uniform pitting fault characteristics of the tooth surface of the gearbox under the condition of variable speed without a healthy phase signal, which improves the accuracy of fault diagnosis and provides a new solution for variable speed. The fault diagnosis of gearboxes under these conditions provides an effective and reliable guarantee.

Claims

1. A gear failure keyless phase angle domain average calculation order analysis method, is characterized in that, comprises the following steps:

1) Instantaneous frequency estimation:

1-1) Collect vibration acceleration signals, perform low-pass phase-preserving filtering on the collected vibration acceleration signals, and obtain low-frequency signals x(t) reflecting rotor rotation;

1-2) The time-frequency distribution of the low-frequency signal x(t) calculated by smoothing the pseudo-Wigner-Ville distribution to obtain a time-frequency distribution diagram;

1-3) Utilize the Viterbi optimal path search algorithm to estimate the instantaneous frequency of the gearbox shaft on the obtained time-frequency distribution diagram;

2) Order spectrum based on instantaneous frequency estimation:

2-1) using the key-phase signal estimation model to carry out point-by-point integration of the instantaneous frequency obtained in step 1-3), to obtain an estimated key-phase signal;

2-2) Using the obtained estimated key phase signal to carry out equiangular resampling to the vibration acceleration signal, to obtain the signal after equiangular resampling;

2-3) performing angle-domain average processing on the signal after equiangular resampling, to obtain a signal after angle-domain average processing;

2-4) Fast Fourier transform is performed on the signal after angular domain average processing to obtain an order spectrum based on instantaneous frequency estimation that can reflect the characteristic information of the gearbox fault, and analyze the gear fault.

2. gear failure according to claim 1 no key phase angle domain average calculation order analysis method, it is characterized in that: the concrete operation of described step 1-2) is:

The smooth pseudo-Wigner-Ville distribution SPW(t,ω) of the low-frequency signal x(t) obtained in step 1-1) is expressed by formula (1):

S S P P W W ((t t,, ω ω)) = = \frac{11}{22 π π} {&Integral; &Integral;}_{- - ∞ ∞}^{+ + ∞ ∞} {&Integral; &Integral;}_{- - ∞ ∞}^{+ + ∞ ∞} x x ((t t + + \frac{11}{22} τ τ)) {x x}^{* *} ((t t - - \frac{11}{22} τ τ)) h h ((τ τ)) g g ((u u)) {e e}^{- - i i τ τ ω ω} d d τ τ d d u u - - - - - - ((11))

The time-frequency distribution diagram can be obtained from formula (1); in formula (1), * means conjugate, t means time, ω means frequency, h(τ) and g(u) are window functions, τ, u represent time axis respectively The translation factor with the frequency axis.

3. The order analysis method for gear fault keyless phase angle domain average calculation according to claim 2, characterized in that: said window function includes a Gauss window and a Hamming window.

4. according to any one of claim 1-3, the gear fault keyless phase angle domain average calculation order analysis method is characterized in that: the concrete operation of described step 1-3) is:

Take the time series [n ₁ ,n ₂ ] on the smooth pseudo-Wigner-Ville distribution of x(t), the time point n∈[n ₁ ,n ₂ ], the frequency value of each time point n is SPW _n , put The array after SPW _n is arranged in descending order is denoted as SPW′ _n ;

In order from large to small, start from 0 and then re-assign by increment of 1, and the time-frequency distribution sequence after assignment is recorded as f(n, SPW′ _n );

The instantaneous frequency is the path that minimizes formula (2):

\overset{^^}{ω ω} ((n no)) = = arg arg min min [[{Σ Σ}_{n no = = {n no}_{11}}^{{n no}_{22}} g g (({ω ω}_{n no},, {ω ω}_{n no + + 11})) + + {Σ Σ}_{n no = = {n no}_{11}}^{{n no}_{22}} f f ((n no,, {SPW SPW}_{n no}^{' '}))]] - - - - - - ((22))

Among them, g(ω _n ,ω _n+1 ) is a penalty factor added to frequency changes with large distances in order to satisfy the continuity of frequency changes.

5. The gear failure keyless phase angle domain average calculation order analysis method according to claim 4, characterized in that: the penalty factor g(ω _n ,ω _n+1 ) is as shown in formula (3):

g g (({ω ω}_{n no},, {ω ω}_{n no + + 11})) = = \{\begin{matrix} 00,, & | | {ω ω}_{n no} - - {ω ω}_{n no + + 11} | | \leq \leq δ δ \\ c c ((| | {ω ω}_{n no} - - {ω ω}_{n no + + 11} | | - - δ δ)),, & | | {ω ω}_{n no} - - {ω ω}_{n no + + 11} | | > > δ δ \end{matrix} - - - - - - ((33))

where δ is the maximum expected value of the instantaneous frequency change of 2 consecutive points.

6. the gear fault keyless phase angle domain average calculation order analysis method according to claim 4, is characterized in that: the concrete operation of described step 2-1) is:

A rotating shaft works under the condition of variable speed, and the instantaneous frequency of rotation is f(t); after time t _n , the rotating shaft rotates n times, and the angle of rotation of the rotating shaft is shown in formula (4):

{&Integral; &Integral;}_{00}^{{t t}_{n no}} 22 π π f f ((t t)) d d t t = = 22 π π n no - - - - - - ((44))

Integrate f(t) point by point, record the integral value t _{n as an integer, and insert a step bond-phase signal at the position of t n} _, which is the estimated bond-phase signal obtained from the instantaneous frequency.

7. gear failure keyless phase angle domain average calculation order analysis method according to claim 4, is characterized in that: the concrete operation of described step 2-2) is:

According to the formula (5), the time node t when the equiangular θ resampling is obtained, and the equiangular resampling is performed to obtain the signal after the equiangular resampling;

t t = = \frac{11}{22 {b b}_{22}} [[\sqrt{44 {b b}_{22} ((θ θ - - {b b}_{00})) + + {b b}_{11}^{22}} - - {b b}_{11}]] - - - - - - ((55))

In formula (5), b ₀ , b ₁ , b ₂ can be obtained according to formula (6);

in

8. gear failure keyless phase angle domain average calculation order analysis method according to claim 4, is characterized in that: the concrete operation of described step 2-3) is:

Denote the signal sequence after equiangular resampling as i=1,2,3,...,N, the signal sequence y _n after the angle domain average processing is obtained by the formula (7) of the angle domain average processing,

{y the y}_{n no} = = \frac{11}{P P} {Σ Σ}_{P P = = 00}^{P P - - 11} {\overset{~ ~}{y the y}}_{n no + + P P N N} - - - - - - ((77))

Wherein, P is the number of segments into which the equiangularly resampled signal is averagely divided, and N is the same number of sampling points in each segment.