CN113702043A - Planetary gear box fault diagnosis method based on POVMD and FDTW under time-varying rotation speed - Google Patents

Planetary gear box fault diagnosis method based on POVMD and FDTW under time-varying rotation speed Download PDF

Info

Publication number
CN113702043A
CN113702043A CN202110911167.3A CN202110911167A CN113702043A CN 113702043 A CN113702043 A CN 113702043A CN 202110911167 A CN202110911167 A CN 202110911167A CN 113702043 A CN113702043 A CN 113702043A
Authority
CN
China
Prior art keywords
signal
time
path
speed
frequency
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202110911167.3A
Other languages
Chinese (zh)
Inventor
李宏坤
孙斌
王朝阁
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Dalian University of Technology
Original Assignee
Dalian University of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Dalian University of Technology filed Critical Dalian University of Technology
Priority to CN202110911167.3A priority Critical patent/CN113702043A/en
Publication of CN113702043A publication Critical patent/CN113702043A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M13/00—Testing of machine parts
    • G01M13/04—Bearings
    • G01M13/045—Acoustic or vibration analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis

Abstract

The invention provides a fault diagnosis method of a planetary gear box under time-varying rotation speed based on POVMD and FDTW, belonging to the technical field of fault diagnosis and vibration signal analysis of rotary machines. The invention differs from the previous order tracking technical idea without a rotating speed measuring device in that the original time scale is projected onto a new time scale by distorting the vibration signal at the measured time-varying rotating speed, wherein the time-varying speed of the shaft is squeezed towards a constant reference speed, so that on this transformed time scale, the fuzzy frequency spectrum of the frequency variation with the speed is squeezed into respective individual peak values, and the required shaft rotating speed information is not required to be so accurate, so that the gear and bearing fault information of the planetary gear box at the time-varying rotating speed can be accurately identified.

Description

Planetary gear box fault diagnosis method based on POVMD and FDTW under time-varying rotation speed
Technical Field
The invention belongs to the technical field of fault diagnosis and vibration signal analysis of rotating machinery, relates to a fault diagnosis method of rotating machinery, and particularly relates to a fault diagnosis method of a planetary gearbox based on a Parameter Optimization Variational Modal Decomposition (POVMD) and Fast Dynamic Time Warping (FDTW) method under time-varying rotating speed and without a rotating speed measuring device, which can be used for fault diagnosis of gears and bearings in the planetary gearbox.
Background
The planetary gear box has the advantages of small volume, large transmission ratio, strong bearing capacity, stable operation, high working efficiency and the like, and is widely applied to large-scale complex mechanical equipment such as helicopters, wind power generation, heavy trucks, ships and warships and the like. However, the planetary gear box is used as a main transmission component of a rotary machine, and key components (a sun gear, a planet gear, a gear ring, a rolling bearing and the like) of the planetary gear box are easy to break down when the planetary gear box works under the conditions of impact, alternating load and large working condition change for a long time, so that huge economic loss and malignant accidents are caused. Therefore, the method has great economic and safety significance for diagnosing the fault of the planetary gearbox and finding out fault information in time.
Under the condition of variable rotating speed, when gears or rolling bearings in the planetary gearbox are locally damaged, vibration response signals collected from the box body are characterized by non-smoothness, nonlinearity, coupling modulation and the like. The traditional signal analysis method such as fast Fourier transform analyzes such signals, because of the fluctuation of speed, the phenomenon of spectrum blurring can occur, and the fault information can not be accurately obtained. In order to eliminate the phenomenon of spectrum blurring due to time-varying velocity, an order tracking method may be adopted. However, the conventional order tracking method requires additional devices to measure the information of the shaft rotation speed, which increases the cost, and in practical cases, the additional devices for measuring the shaft rotation speed cannot be installed. In order to overcome the problems, a method of order tracking without rotating speed measurement is provided, and shaft rotating speed information is estimated based on time-frequency ridge extraction of time-frequency distribution, but is limited by the influence of time-frequency resolution, and is difficult to extract accurately. Thus, an adaptive signal approach is proposed to estimate shaft speed information for gearbox fault signals. An Empirical Mode Decomposition (EMD) method, which is an adaptive signal decomposition method used earlier, can adaptively decompose a signal into a sum of a series of natural modal functions according to local scale characteristics of the signal itself, thereby revealing the internal nature of the signal. Although the EMD method has many advantages, there are problems of modal aliasing, end-point effects, and lack of theoretical basis. In order to make up for the defects of the EMD method, improved methods such as Ensemble Empirical Mode Decomposition (EEMD) and Ensemble Local Mean Decomposition (ELMD) of noise-aided analysis have been developed, but the selection mechanism of two key parameters (amplitude and integration times of noise) for adding white noise in EEMD and ELMD is not clear, and there are problems of residual noise pollution and large computation amount in the signal reconstruction process after noise is added. Unlike the analysis methods described above, the Variational Mode Decomposition (VMD) is a novel adaptive signal decomposition algorithm for the decomposition of non-stationary signals. Because of having a perfect mathematical theoretical basis, the method is favored by the learners. However, when the VMD algorithm is applied, the mode number K and the penalty factor α have a large influence on the decomposition result of the VMD, where too many K cause repeated occurrence of center frequency, and too few result in incomplete decomposition; if the value of alpha is too large, the bandwidth of the frequency band of the component is narrowed, and if the value of alpha is too small, the bandwidth of the frequency band is widened, so that the modal component identification is influenced, and further the estimation of the shaft rotating speed is influenced.
In addition, in order to further improve the order tracking technology of the non-rotating speed measuring device, the phenomenon of spectrum blurring caused by rotating speed change needs to be eliminated from the time domain. Fast Dynamic Time Warping (FDTW), which eliminates the phenomenon of spectrum ambiguity caused by the rotational speed by distorting the path of the normalized shaft vibration signal estimated from the shaft rotational speed and the shaft vibration signal assumed to be at a constant rotational speed, distorting the original vibration signal by the path, and resampling.
Disclosure of Invention
In view of the problem that fault information cannot be accurately obtained due to the phenomenon of spectrum blurring caused by the time-varying rotating speed of the planetary gearbox and the problem that rotating speed information cannot be accurately obtained due to the fact that time-frequency resolution is not high based on a time-frequency image, the invention provides a fault diagnosis method based on parameter optimization variational modal decomposition and fast dynamic time warping (POVMD-FDTW), and the fault diagnosis method overcomes the defects.
In order to achieve the purpose, the invention adopts the following technical scheme:
a fault diagnosis method for a planetary gear box under a time-varying rotating speed based on POVMD-FDTW is characterized by comprising the following steps: the method comprises the following steps:
step S1: an acceleration sensor is arranged right above a planetary gearbox shell to collect fault vibration signals x (t) (t is a time sequence);
step S2: optimizing the modal number K and the penalty factor alpha for VMD decomposition of the fault signal x (t) of the planetary gearbox;
step S3: performing POVMD decomposition on the planetary gearbox fault signal x (t) based on the parameters in step S2;
step S4: calculating instantaneous frequency of each eigenmode function (IMF) after decomposition, and selecting mesh frequency fmNearby IMFs;
step S5: performing Hilbert transform based on the IMF selected in step S4, and normalizing the selected IMF to ximf, obtaining the instantaneous phase phiinst(t) information, deriving time t, dividing by number of teeth Z of ring gearrObtaining the absolute rotation frequency of the sun wheel according to the transmission ratio i
Figure BDA0003201450590000021
Information;
step S6: based on the sun gear shaft rotation speed in step S5
Figure BDA0003201450590000022
Information, establishing a shaft vibration signal x at a time-varying rotational speedsh(t)；
Step S7: a shaft reference signal y (t) at a fixed rotational speed is established.
Step S8: for the shaft vibration signal x at the time-varying rotation speed acquired in step S6sh(t) and the shaft reference signal y (t) at the constant rotation speed obtained in step S7, implementing FDTW algorithm to achieve the optimal alignment therebetween, obtaining the optimal twisted path W, and being capable of eliminating the phase difference between the two signals caused by the rotation speed variation;
step S9: distorting the planetary gearbox fault signal x (t) based on the optimal distortion path W obtained in the step S8 to obtain a distortion signal x1(k) And resampled so that x1(k) Restore to the original signal length to obtain(n) the re-mined signal z;
step S10: performing fast Fourier transform on the resampled signal z (n) obtained in the step S9 to obtain an order spectrum, and judging fault information;
further, the process of optimizing the number K of modes of VMD decomposition and the penalty factor α in step S2 is as follows: firstly, determining the optimal modal number, and decreasing the modal number in a descending manner, so that a modal signal without frequency aliasing can be decomposed under the condition of the narrowest bandwidth; in order to determine the optimal bandwidth of each modal signal, firstly, the number of the modes is decreased once, then, the penalty factor value is decreased progressively, so that the condition of bandwidth aliasing does not occur among the modes, and when the condition of bandwidth aliasing is calculated, the value of one-step penalty factor is returned, and the optimal mode number K and the penalty factor alpha are obtained.
Further, in step S3, based on the parameters in step S2, the specific process of performing the popmd decomposition on the planetary gearbox fault signal x (t) is as follows: for the kth mode, the correlated analysis signal is computed by the hilbert transform to obtain a single-sided spectrum. The spectrum of each mode is moved to the baseband by an index, while the estimated center frequency is also adjusted to the position of the center of gravity of the corresponding power spectrum of the mode. The bandwidth is measured by gaussian smoothing of the input signal, i.e. the square of the two-norm of the gradient. Therefore, the constraint variation problem is given by:
Figure BDA0003201450590000031
wherein f (t) is a given signal, ukIs the kth mode, ωkIs the center frequency, k is the number of modes,
Figure BDA0003201450590000032
is the partial derivative of the function over time t, δ (t) is the dirac distribution, and is the convolution operator.
The quadratic penalty term alpha and the Lagrangian multiplier lambda are introduced, so that the problem is converted into an unconstrained problem. Because of the good convergence of the quadratic penalty term under finite weights and the strict enforcement of the λ constraint. Giving an equation containing the optimal solution for λ and the quadratic penalty term α:
Figure BDA0003201450590000033
solving the above equation using an alternating direction multiplier optimization algorithm:
Figure BDA0003201450590000041
Figure BDA0003201450590000042
further, the specific process of calculating the instantaneous frequency of each eigenmode function (IMF) after the decomposition in step S4 is as follows: and performing Hilbert transform on each IMF to obtain instantaneous phase information of each IMF, performing time derivation to obtain instantaneous frequency information of each IMF, finding IMFs near the meshing frequency, and further performing subsequent analysis.
Further, the specific process of obtaining the rotation speed of the sun gear shaft in step S5 is as follows: the normalized signal based on the IMF selected in step S4 is as follows:
Figure BDA0003201450590000043
Figure BDA0003201450590000044
in the formula (I), the compound is shown in the specification,
Figure BDA0003201450590000045
is ximf(t) Hilbert transform. The instantaneous phase and sun shaft speed information is calculated as follows:
Figure BDA0003201450590000046
Figure BDA0003201450590000047
in the formula, ZrI is the number of teeth of the ring gear, and i is the transmission ratio of the planetary gearbox.
Further, step S6 establishes x of the shaft vibration signal at the time-varying rotation speedshThe specific process of (t) is as follows:
Figure BDA0003201450590000048
the initial phase theta is found by a small increment from 0 to 2 pi
Figure BDA0003201450590000049
And xsh(t) phase angle estimation of the minimum mean square error.
Further, the specific process of establishing the shaft vibration signal y (t) at the constant rotation speed in step S7 is as follows:
y(t)＝cos(2πfst+β)
speed f at constant rotational speedsFrom x in step S6sh(t) is determined by the number of peaks p, fsHas a value range of
Figure BDA00032014505900000410
To
Figure BDA00032014505900000411
Wherein t iszongFor the total sampling time, the initial phase information beta is increased from 0 to 2 pi by a small increment, and y (t) and x are foundsh(t) phase angle estimation of the minimum mean square error.
Further, step S8 is implemented by using FDTW algorithm to obtain shaft vibration signal x at time-varying rotation speedshThe specific process of the optimal alignment path of the shaft reference signal y (t) under the rotation speed and the (t) is as follows: x is the number ofshThe lengths of (t) and y (t) are M and N, respectively:
xsh(t)＝x1,x2,...,xi,...,xM
y(t)＝y1,y2,...,yi,...,yN
twisted path, W:
W＝w1,w2,...,wk,...,wK
k is the length of the bend, and the kth element of the curved path is Wk(i, j). Wherein i is from xsh(t), j is the time series index from y (t), and if the curved path is (i, j), then it represents the time series xsh(t) aligns with the jth element of the time series y (t).
Searching a shortest path by using Euclidean distance to obtain xsh(t) and y (t) the shortest distance between the two time sequences is calculated as follows:
Figure BDA0003201450590000051
d(xsh(i),y(j))＝(xsh(i)-y(j))2
where d represents the distance between two points, an M × N distance matrix is obtained by calculating the distance between different points of the two time series:
Figure BDA0003201450590000052
selecting a minimum twist cost:
Figure BDA0003201450590000053
the twist cost distance in the iterative process is as follows:
Figure BDA0003201450590000054
the initial value is D (2,2) ═ D (x)sh(1),y(1))+d(xsh(2),y(2))
Then, a matrix D of M × N size is obtained, and a warp path with the shortest distance is found from the matrix D as W.
The warp path must satisfy the following condition:
(1)w1at the beginning of each time series, (1,1) the warp path must start;
(2)wK(N, M), the warped path must end at the end of both time series;
and the fast dynamic time warping is to reduce the iteration range and complexity in the iteration process.
Further, the specific process of step S9 distorting the planetary gearbox fault signal x (n) using the optimal distortion path W and resampling to restore the distorted signal to the original length is: obtaining a path W (i) from the optimal distortion path W obtained in step S8, and distorting the planetary gearbox fault signal x (n) by using the path to obtain a distortion signal x1(k) The resampled signal z (n). The principle of resampling to the original signal length is to start from W (j) in the warped path W to the length K of the warped path W, if Wk-1J "and j are not equal, then z (n) x1(k) Otherwise, the largest of the two is taken.
Further, the specific process of step S10 is: and for the resampling signal z (n), carrying out fast Fourier transform to obtain an order spectrum, thereby judging fault information.
Compared with the prior art, the invention has the advantages that the invention has the following characteristics:
(1) the invention not only gives play to the advantage that the POVMD algorithm self-adaptive parameter optimizes and decomposes signals to obtain single-component signals and further estimates the shaft-to-frequency, but also eliminates the phenomenon of frequency spectrum ambiguity caused by rotation speed fluctuation from a time domain by utilizing the FDTW algorithm, and solves the problems of interpolation error of equal-angle resampling, extra need of a rotation speed measuring device and need of accurate rotation speed information in the traditional method.
(2) The parameter optimization in the POVMD algorithm of the invention overcomes the problems that the number of modal number K is too large, so that the central frequency appears repeatedly, the decomposition is incomplete if the number of the modal number K is too small, and the identification of modal components is influenced because the value of a penalty factor alpha is too large, the bandwidth of the frequency band of the components is narrowed, and the bandwidth of the frequency band is widened if the number of the modal number K is too small.
(3) The FDTW algorithm of the present invention differs from previous thought of order tracking technology without a rotational speed measuring device in that the original time scale is projected onto a new time scale by distorting the vibration signal at the measured time-varying rotational speed, wherein the time-varying speed of the shaft is squeezed towards a constant reference speed, so that on this transformed time scale the blurred spectrum whose frequency varies with speed is also squeezed into individual peaks.
(4) The invention not only measures the vibration signal of the gearbox at the time-varying rotation speed, but also needs approximate rotation speed information of the shaft, and the required rotation speed information is not as accurate as that required by the traditional order analysis.
Drawings
FIG. 1 is a flow chart of a diagnostic method of the present invention;
FIG. 2 is a schematic diagram of a planetary gear box fault test stand at a time-varying rotation speed in an embodiment of the invention;
FIG. 3 is a time domain waveform of a sun gear tooth breakage fault vibration signal and its Fourier spectrum in an embodiment of the present invention;
FIG. 4 is a POVMD decomposition of a sun fault signal in an embodiment of the present invention;
FIG. 5 is a POVMD decomposition of the sun fault signal for the instantaneous frequency of each IMF in an embodiment of the present invention;
FIG. 6 is the result of an estimation of the sun gear shaft speed in an embodiment of the present invention;
FIG. 7 is the result of a vibration signal at a time varying rotational speed of the sun gear shaft in an embodiment of the present invention;
FIG. 8 is a schematic illustration of reference signals in an embodiment of the present invention;
FIG. 9 is a graph of the warped path after implementation of the FDTW algorithm in an embodiment of the present invention;
FIG. 10 is a graph of a signal after the FDTW algorithm has been warped and resampled in an embodiment of the present invention;
FIG. 11 is the order spectrum of the last diagnosis in the example of the present invention.
In the figure: 1 driving a motor; 2, a frequency converter; 3, a torque sensor; 4, bearing seats; 5, collecting cards; 6 a planetary gearbox; 7 a tension controller; 8 magnetic powder loader.
Detailed Description
The following describes the specific implementation process of the present invention with reference to the drawings and examples.
The invention provides a fault diagnosis method of a planetary gear box under a time-varying rotating speed based on POVMD-FDTW, a flow chart is shown in figure 1, and the main steps are as follows:
step S1: an acceleration sensor is arranged right above a planetary gearbox shell to collect fault vibration signals x (t) (t is a time sequence);
step S2: optimizing the modal number K and the penalty factor alpha for VMD decomposition of the fault signal x (t) of the planetary gearbox;
step S3: performing POVMD decomposition on the planetary gearbox fault signal x (t) based on the parameters in step S2;
step S4: calculating instantaneous frequency of each eigenmode function (IMF) after decomposition, and selecting mesh frequency fmNearby IMFs;
step S5: performing Hilbert transform based on the IMF selected in step S4, and normalizing the selected IMF to
Figure BDA0003201450590000071
Obtaining an instantaneous phase phiinst(t) information, deriving time t, dividing by number of teeth Z of ring gearrObtaining the absolute rotation frequency of the sun wheel according to the transmission ratio i
Figure BDA0003201450590000072
Information;
step S6: based on the sun gear shaft rotation speed in step S5
Figure BDA0003201450590000073
Information, establishing a shaft vibration signal x at a time-varying rotational speedsh(t)；
Step S7: a shaft reference signal y (t) at a fixed rotational speed is established.
Step S8: for those obtained in step S6Shaft vibration signal x under time-varying rotating speedsh(t) and the shaft reference signal y (t) at the constant rotation speed obtained in step S7, implementing FDTW algorithm to achieve the optimal alignment therebetween, obtaining the optimal twisted path W, and being capable of eliminating the phase difference between the two signals caused by the rotation speed variation;
step S9: distorting the fault signal x (n) of the planetary gearbox based on the optimal distortion path W obtained in the step S8 to obtain a distortion signal x1(k) And resampled so that x1(k) Restoring to the length of the original signal to obtain a re-picked signal z (n);
step S10: performing fast Fourier transform on the resampled signal z (n) obtained in the step S9 to obtain an order spectrum, and judging fault information;
example 1:
and (3) building a fault diagnosis laboratory table of the planetary gearbox, as shown in figure 2. The experiment table mainly comprises a driving motor 1, a planetary gear box 6, a magnetic powder loader 8 and the like. The structural parameters of the planetary gearbox are shown in the table 1. A portion is cut out of a certain tooth of the sun gear as a failure by a wire cutting technique. During the experiment, the fault sun gear is arranged in the planetary gear box for experimental data acquisition. The acceleration sensor is arranged on the vertical, horizontal and axial measuring points of the shell of the planetary gear box. The input rotating speed of the motor is linearly increased by 5-8Hz, the data sampling frequency is set to be 12800Hz, and the sampling time duration is 2 s. Under the experimental condition, the corresponding orders of the local faults of the sun gear of the planetary gearbox can be calculated and are shown in the table 2, wherein m and k are positive integers.
TABLE 1 planetary Gear parameters (Unit/Unit)
Figure BDA0003201450590000081
TABLE 2 sun wheel local failure corresponding characteristic orders
Figure BDA0003201450590000082
The flow of the planetary gear box fault diagnosis method based on the POVMD-FDTW time-varying rotation speed is shown in figure 1, and the method specifically comprises the following steps:
the method comprises the following steps: and collecting fault vibration signals of the sun gear of the planetary gear box in the fault diagnosis experiment table under the time-varying rotating speed by using an acceleration sensor. Fig. 3 is a time domain waveform of a sun gear tooth breakage fault signal acquired by the acceleration sensor in the vertical direction and a fourier spectrum thereof. As can be seen from fig. 3, regular periodic impact characteristics are difficult to observe in the time domain waveform of the vibration signal, the frequency spectrum of the sun gear fault characteristic frequency in the frequency spectrum is blurred due to the change of the rotating speed, and the amplitude of noise and other irrelevant interference frequencies is more prominent. Therefore, the characteristic information representing the health state of the sun gear cannot be extracted from the original fault signal by the traditional time domain and frequency domain analysis method.
Step two: firstly, determining the optimal modal number, and decreasing the modal number progressively so as to decompose a modal signal without frequency aliasing under the condition of the narrowest bandwidth; in order to determine the optimal bandwidth of each modal signal, firstly, the number of the modes is decreased once, then, the penalty factor value is decreased progressively, so that the condition of bandwidth aliasing does not occur among the modes, and when the condition of bandwidth aliasing is calculated, the value of one-step penalty factor is returned, and the optimal mode number K and the penalty factor alpha are obtained. The number of modes K resolved by the present invention in this example is 10 and the penalty factor α is 3300.
Step three: based on the parameter information in the second step, performing POVMD decomposition on the fault vibration signal acquired in the first step, as shown in FIG. 4, specifically including the following steps: an equation is given for the optimal solution containing λ and the quadratic penalty term α:
Figure BDA0003201450590000091
solving the above equation using an alternating direction multiplier optimization algorithm:
Figure BDA0003201450590000092
Figure BDA0003201450590000093
step four: and performing Hilbert transform on each IMF in the third step to obtain instantaneous phase information of each IMF, performing time derivation to obtain instantaneous frequency information of each IMF, and finding the IMF near the meshing frequency as shown in FIG. 5 for subsequent analysis.
Step five: obtaining the rotation speed of the sun gear shaft based on the IMF selected in the step four, as shown in fig. 6, the specific process is as follows: the IMF normalized signal is as follows:
Figure BDA0003201450590000094
Figure BDA0003201450590000095
in the formula (I), the compound is shown in the specification,
Figure BDA0003201450590000096
is ximf(t) Hilbert transform. The instantaneous phase and sun shaft speed information is calculated as follows:
Figure BDA0003201450590000097
Figure BDA0003201450590000098
in the formula, ZrI is the number of teeth of the ring gear, and i is the transmission ratio of the planetary gearbox.
Step six: establishing x of shaft vibration signal under time-varying rotating speedsh(t), as shown in fig. 7, the specific process is:
Figure BDA0003201450590000099
the initial phase theta is found by a small increment from 0 to 2 pi
Figure BDA00032014505900000910
And xsh(t) phase angle estimation of the minimum mean square error.
Step seven, establishing a shaft vibration signal y (t) at a fixed rotating speed based on the step six, wherein the specific process is as shown in fig. 8: y (t) cos (2 pi f)st+β)
Speed f at constant rotational speedsFrom step six xsh(t) is determined by the number of peaks p, fsHas a value range of
Figure BDA00032014505900000911
To
Figure BDA00032014505900000912
Finding y (t) and x by a small increment of 0 to 2 pi of the initial phase information betash(t) phase angle estimation of the minimum mean square error.
Step eight, adopting a time-varying rotating speed lower shaft vibration signal x realized by an FDTW algorithmsh(t) and the optimal twist path W of the shaft reference signal y (t) at the rotation speed, as shown in fig. 9, the specific process is as follows: x is the number ofshThe lengths of (t) and y (t) are M and N, respectively:
xsh(t)＝x1,x2,...,xi,...,xM
y(t)＝y1,y2,...,yi,...,yN
twisted path, W:
W＝w1,w2,...,wk,...,wK
k is the length of the bend, and the kth element of the curved path is Wk(i, j). Wherein i is from xsh(t), j is the time series index from y (t), if the curved path is (i, j), thenRepresenting a time series xsh(t) aligns the ith element with the jth element of the time series.
Searching a shortest path by using Euclidean distance to obtain xsh(t) and y (t) the shortest distance between the two time sequences is calculated as follows:
Figure BDA0003201450590000101
d(xsh(i),y(j))＝(xsh(i)-y(j))2
where d represents the distance between two points, an M × N distance matrix is obtained by calculating the distance between different points of the two time series:
Figure BDA0003201450590000102
selecting a minimum twist cost:
Figure BDA0003201450590000103
the twist cost distance in the iterative process is as follows:
Figure BDA0003201450590000104
the initial value is D (2,2) ═ D (x)sh(1),y(1))+d(xsh(2),y(2))
Then, a matrix D of M × N size is obtained, and a warp path with the shortest distance is found from the matrix D as W.
Step nine: and (3) distorting the fault signal x (n) of the planetary gearbox based on the optimal distortion path W obtained in the step eight, and resampling to restore the distorted signal to the original length, as shown in fig. 10, wherein the specific process is as follows: obtaining a path W (i) from the optimal distorted path W obtained in the step eight, and distorting the fault signal x (n) of the planetary gearbox by using the path to obtain a distortion signalNumber x1(k) The resampled signal z (n). The principle of resampling to the original signal length is to start from W (j) in the warped path W to the length K of the warped path W, if Wk-1J "and j are not equal, then z (n) x1(k) Otherwise, the largest of the two is taken.
Performing fast fourier transform analysis on the resampled signal z (n) obtained in the ninth step to obtain an order spectrum, as shown in fig. 11. As can be seen from FIG. 11, the absolute transfer order O of the sun gear and the mesh order O are clearly evident in the order spectrummOrder of side frequency Om±mOsAnd k are positive integers (m ═ 1,2, 3; k ═ 1, 2). The above analysis therefore shows that the sun wheel is malfunctioning, which is consistent with the experimental setup. Thereby verifying the effectiveness and feasibility of the method of the invention.
In conclusion, the fault diagnosis method for the planetary gear box based on the POVMD-FDTW time-varying rotation speed can successfully extract weak fault characteristic information and realize accurate identification and diagnosis of early faults of the planetary gear box. FDTW projects the original time scale onto a new time scale by distorting the vibration signal at the measured time-varying rotational speed, wherein the time-varying speed of the shaft is squeezed towards a constant reference speed, so that on this transformed time scale the blurred spectrum whose frequency varies with speed is also squeezed into individual peaks. In addition to measuring the vibration signal of the gearbox at time varying rotational speeds, the required shaft rotational speed information is approximate and does not need to be as accurate as the rotational speed information required by conventional order analysis.
The above description is only a preferred embodiment of the present invention, and all equivalent modifications made within the scope of the claims of the present invention should be covered by the present invention.

Claims (1)

1. A fault diagnosis method of a planetary gear box under time-varying rotating speed based on POVMD and FDTW is characterized in that: the method comprises the following steps:
step S1: an acceleration sensor is arranged right above a planetary gearbox shell to collect fault vibration signals x (t), wherein t is a time sequence;
step S2: optimizing the modal number K and the penalty factor alpha for VMD decomposition of the fault signal x (t) of the planetary gearbox;
firstly, determining the optimal modal number, and decreasing the modal number to decompose the modal signal without frequency aliasing under the condition of the narrowest bandwidth; in order to determine the optimal bandwidth of each modal signal, firstly, the number of the modes is decreased once, then, the penalty factor value is decreased progressively, so that the condition of bandwidth aliasing does not occur among the modes, and when the condition of bandwidth aliasing is calculated, the value of one-step penalty factor is returned, so that the optimal mode number K and the penalty factor alpha are obtained;
step S3: performing POVMD decomposition on the planetary gearbox fault signal x (t) based on the parameters in step S2;
for the kth mode, calculating a related analysis signal through Hilbert transform to obtain a single-sided spectrum; the frequency spectrum of each mode is moved to a base band through an index, and meanwhile, the estimated center frequency is adjusted to the position of the center of gravity of the power spectrum corresponding to the mode; the bandwidth is measured by gaussian smoothing of the input signal, i.e. the square of the two-norm of the gradient; therefore, the constraint variation problem is given by:
Figure FDA0003201450580000011
wherein f (t) is a given signal, ukIs the kth mode, ωkIs the center frequency, k is the number of modes,
Figure FDA0003201450580000012
is the partial derivative of the function over time t, δ (t) is the dirac distribution, a convolution operator;
introducing a secondary penalty term alpha and a Lagrangian multiplication operator lambda so as to convert the problem into an unconstrained problem; because of the good convergence of the quadratic penalty term under finite weights and the strict enforcement of the λ constraint; giving an equation containing the optimal solution for λ and the quadratic penalty term α:
Figure FDA0003201450580000013
solving the above equation using an alternating direction multiplier optimization algorithm:
Figure FDA0003201450580000014
Figure FDA0003201450580000021
step S4: calculating instantaneous frequency of each eigenmode function (IMF) after decomposition, and selecting mesh frequency fmNearby IMFs;
performing Hilbert transformation on each IMF to obtain instantaneous phase information of each IMF, performing time derivation to obtain instantaneous frequency information of each IMF, finding IMFs near the meshing frequency, and further performing subsequent analysis;
step S5: performing Hilbert transform based on the IMF selected in step S4, and normalizing the selected IMF to
Figure FDA0003201450580000022
Obtaining an instantaneous phase phiinst(t) information, deriving time t, dividing by number of teeth Z of ring gearrObtaining the absolute rotation frequency of the sun wheel according to the transmission ratio i
Figure FDA0003201450580000023
Information;
the normalized signal based on the IMF selected in step S4 is as follows:
Figure FDA0003201450580000024
Figure FDA0003201450580000025
in the formula (I), the compound is shown in the specification,
Figure FDA0003201450580000026
is ximf(t) a hilbert transform; the instantaneous phase and sun shaft speed information is calculated as follows:
Figure FDA0003201450580000027
Figure FDA0003201450580000028
in the formula, ZrThe number of teeth of the gear ring is shown, i is the transmission ratio of the planetary gear box;
step S6: based on the sun gear shaft rotation speed in step S5
Figure FDA0003201450580000029
Information, establishing a shaft vibration signal x at a time-varying rotational speedsh(t) is represented by the following formula:
Figure FDA00032014505800000210
the initial phase theta is found by a small increment from 0 to 2 pi
Figure FDA00032014505800000211
And xsh(t) a phase angle estimate of the minimum mean square error;
step S7: establishing a shaft reference signal y (t) at a fixed rotating speed according to the following formula:
y(t)＝cos(2πfst+β)
speed f at constant rotational speedsFrom x in step S6sh(t) is determined by the number of peaks p, fsHas a value range of
Figure FDA00032014505800000212
To
Figure FDA00032014505800000213
Wherein t iszongFor the total sampling time, the initial phase information beta is increased from 0 to 2 pi by a small increment, and y (t) and x are foundsh(t) a phase angle estimate of the minimum mean square error;
step S8: for the shaft vibration signal x at the time-varying rotation speed acquired in step S6sh(t) and the shaft reference signal y (t) at the constant rotation speed obtained in step S7, implementing FDTW algorithm to achieve the optimal alignment therebetween, obtaining the optimal twisted path W, and eliminating the phase difference between the two signals caused by the rotation speed variation;
shaft vibration signal x realized by adopting FDTW algorithm at time-varying rotating speedshThe specific process of the optimal alignment path of the shaft reference signal y (t) under the rotation speed and the (t) is as follows: x is the number ofshThe lengths of (t) and y (t) are M and N, respectively:
xsh(t)＝x1,x2,...,xi,...,xM
y(t)＝y1,y2,...,yi,...,yN
twisted path, W:
W＝w1,w2,...,wk,...,wK
k is the length of the bend, and the kth element of the curved path is Wk(ii) as (i, j); wherein i is from xsh(t), j is the time series index from y (t), and if the curved path is (i, j), then it represents the time series xsh(t) the ith element is aligned with the jth element of the time series y (t);
searching a shortest path by using Euclidean distance to obtain xsh(t) and y (t) the shortest distance between the two time sequences is calculated as follows:
Figure FDA0003201450580000031
d(xsh(i),y(j))＝(xsh(i)-y(j))2
where d represents the distance between two points, an M × N distance matrix is obtained by calculating the distance between different points of the two time series:
Figure FDA0003201450580000032
selecting a minimum twist cost:
Figure FDA0003201450580000033
the twist cost distance in the iterative process is as follows:
Figure FDA0003201450580000034
the initial value is D (2,2) ═ D (x)sh(1),y(1))+d(xsh(2),y(2))
Obtaining a matrix D with the size of M multiplied by N, and finding out a distortion path with the shortest distance from the matrix D as W;
the warp path must satisfy the following condition:
(1)w1at the beginning of each time series, (1,1) the warp path must start;
(2)wK(N, M), the warped path must end at the end of both time series;
the fast dynamic time warping is to continuously reduce the iteration range in the iteration process and reduce the complexity of iteration;
step S9: distorting the planetary gearbox fault signal x (t) based on the optimal distortion path W obtained in the step S8 to obtain a distortion signal x1(k) And resampled so that x1(k) Restoring to the length of the original signal to obtain a re-picked signal z (n);
the best twisted path W obtained from step S8 is obtainedThe path w (i) is used for twisting the fault signal x (n) of the planetary gearbox to obtain a twisting signal x1(k) The resampled signal z (n); the principle of resampling to the original signal length is to start from W (j) in the warped path W to the length K of the warped path W, if Wk-1J "and j are not equal, then z (n) x1(k) Otherwise, taking the largest of the two;
step S10: and (4) performing fast Fourier transform on the resampled signal z (n) obtained in the step (S9) to obtain an order spectrum, and judging fault information.
CN202110911167.3A 2021-08-09 2021-08-09 Planetary gear box fault diagnosis method based on POVMD and FDTW under time-varying rotation speed Pending CN113702043A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110911167.3A CN113702043A (en) 2021-08-09 2021-08-09 Planetary gear box fault diagnosis method based on POVMD and FDTW under time-varying rotation speed

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110911167.3A CN113702043A (en) 2021-08-09 2021-08-09 Planetary gear box fault diagnosis method based on POVMD and FDTW under time-varying rotation speed

Publications (1)

Publication Number Publication Date
CN113702043A true CN113702043A (en) 2021-11-26

Family

ID=78651999

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110911167.3A Pending CN113702043A (en) 2021-08-09 2021-08-09 Planetary gear box fault diagnosis method based on POVMD and FDTW under time-varying rotation speed

Country Status (1)

Country Link
CN (1) CN113702043A (en)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114216676A (en) * 2021-11-30 2022-03-22 上海海事大学 Compound fault diagnosis method for planetary gearbox without tachometer under time-varying working condition
CN114593908A (en) * 2022-03-01 2022-06-07 西人马(深圳)科技有限责任公司 Gear fault analysis method and device
CN116010806A (en) * 2023-03-28 2023-04-25 湖北工业大学 Time-frequency analysis method for rotary mechanical time-varying multi-component signal
CN117269701A (en) * 2023-11-21 2023-12-22 川力电气有限公司 High-voltage switch cabinet partial discharge positioning method based on artificial intelligence

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107702908A (en) * 2017-10-12 2018-02-16 国网山东省电力公司莱芜供电公司 GIS mechanical oscillation signal Time-Frequency Analysis Methods based on VMD self adapting morphologies
CN108426715A (en) * 2018-06-13 2018-08-21 福州大学 Rolling bearing Weak fault diagnostic method based on PSO-VMD-MCKD
CN109029977A (en) * 2018-07-12 2018-12-18 福州大学 A kind of epicyclic gearbox Incipient Fault Diagnosis method based on VMD-AMCKD
CN110146281A (en) * 2019-06-06 2019-08-20 南京航空航天大学 A kind of epicyclic gearbox method for diagnosing faults based on VMD-SDAE

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107702908A (en) * 2017-10-12 2018-02-16 国网山东省电力公司莱芜供电公司 GIS mechanical oscillation signal Time-Frequency Analysis Methods based on VMD self adapting morphologies
CN108426715A (en) * 2018-06-13 2018-08-21 福州大学 Rolling bearing Weak fault diagnostic method based on PSO-VMD-MCKD
CN109029977A (en) * 2018-07-12 2018-12-18 福州大学 A kind of epicyclic gearbox Incipient Fault Diagnosis method based on VMD-AMCKD
CN110146281A (en) * 2019-06-06 2019-08-20 南京航空航天大学 A kind of epicyclic gearbox method for diagnosing faults based on VMD-SDAE

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
MADHURJYA D. CHOUDHURY 等: "A Methodology to Handle Spectral Smearing in Gearboxes Using Adaptive Mode Decomposition and Dynamic Time Warping", 《IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT》 *
李宏坤 等: "基于POVMD和CAF的低转速齿轮箱故障诊断", 《振动、测试与诊断》 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114216676A (en) * 2021-11-30 2022-03-22 上海海事大学 Compound fault diagnosis method for planetary gearbox without tachometer under time-varying working condition
CN114593908A (en) * 2022-03-01 2022-06-07 西人马(深圳)科技有限责任公司 Gear fault analysis method and device
CN116010806A (en) * 2023-03-28 2023-04-25 湖北工业大学 Time-frequency analysis method for rotary mechanical time-varying multi-component signal
CN116010806B (en) * 2023-03-28 2023-09-01 湖北工业大学 Time-frequency analysis method for rotary mechanical time-varying multi-component signal
CN117269701A (en) * 2023-11-21 2023-12-22 川力电气有限公司 High-voltage switch cabinet partial discharge positioning method based on artificial intelligence
CN117269701B (en) * 2023-11-21 2024-02-02 川力电气有限公司 High-voltage switch cabinet partial discharge positioning method based on artificial intelligence

Similar Documents

Publication Publication Date Title
CN113702043A (en) Planetary gear box fault diagnosis method based on POVMD and FDTW under time-varying rotation speed
CN109520738B (en) Rotating machinery rolling bearing fault diagnosis method based on order spectrum and envelope spectrum
CN105806613A (en) Planetary gear case fault diagnosis method based on order complexity
CN102288286A (en) Method for analyzing and evaluating measure point precision of gearbox in vibration acceleration sensor
CN108801630B (en) Gear fault diagnosis method for single-channel blind source separation
Wei et al. Variational nonlinear component decomposition for fault diagnosis of planetary gearboxes under variable speed conditions
CN112665851B (en) Key-free phase change rotating speed gearbox fault diagnosis method
Shang et al. An intelligent fault diagnosis system for newly assembled transmission
CN110987438A (en) Method for detecting periodical vibration impact signals of hydraulic generator in variable rotating speed process
CN112668518A (en) VMSST time-frequency analysis method for vibration fault signal
CN109596349A (en) A kind of decelerator trouble diagnostic method based on VMD and PCT
Lin et al. A review and strategy for the diagnosis of speed-varying machinery
CN111881736A (en) Rolling bearing early fault diagnosis method based on bandwidth Fourier decomposition
CN114216676A (en) Compound fault diagnosis method for planetary gearbox without tachometer under time-varying working condition
Shi et al. Instantaneous frequency synchronized generalized stepwise demodulation transform for bearing fault diagnosis
CN111610023B (en) Speed reducer noise evaluation method and device and handheld speed reducer noise evaluation instrument
Liu et al. Fault diagnosis of wind turbines under nonstationary conditions based on a novel tacho-less generalized demodulation
CN108398260B (en) Method for quickly evaluating instantaneous angular speed of gearbox based on mixed probability method
CN112362343A (en) Distributed fault feature extraction method for gearbox under variable rotating speed based on frequency modulation dictionary
Zhao et al. Vold-Kalman generalized demodulation for multi-faults detection of gear and bearing under variable speeds
CN111413095A (en) Instantaneous angular velocity-based distributed fault diagnosis and analysis method for planetary bearing
CN113340631A (en) Torsional vibration testing device and signal analysis method
CN104459186A (en) Tachometer-free order-ratio analyzing method based on sparse segmentation fitting and integral approximation
CN107490477A (en) The Fault Diagnosis of Gear Case method compared based on frequency spectrum kernel density function correlation
Wang et al. Application of Improved Particle Swarm Optimization in Gear Fault Diagnosis of Automobile Transmission.

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
RJ01 Rejection of invention patent application after publication

Application publication date: 20211126

RJ01 Rejection of invention patent application after publication