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

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

Info

Publication number
CN103499443B
CN103499443B CN201310416504.7A CN201310416504A CN103499443B CN 103499443 B CN103499443 B CN 103499443B CN 201310416504 A CN201310416504 A CN 201310416504A CN 103499443 B CN103499443 B CN 103499443B
Authority
CN
China
Prior art keywords
signal
frequency
omega
time
domain average
Prior art date
Application number
CN201310416504.7A
Other languages
Chinese (zh)
Other versions
CN103499443A (en
Inventor
訾艳阳
何正嘉
曹宏瑞
万志国
冯超亮
何水龙
Original Assignee
西安交通大学
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 西安交通大学 filed Critical 西安交通大学
Priority to CN201310416504.7A priority Critical patent/CN103499443B/en
Publication of CN103499443A publication Critical patent/CN103499443A/en
Application granted granted Critical
Publication of CN103499443B publication Critical patent/CN103499443B/en

Links

Abstract

The invention provides a kind of gear distress without key phase angular domain average computation order analysis method, first the vibration acceleration signal of collection is carried out calculating its time-frequency distributions by smoothed pseudo wigner ville disstribution after low pass protects phase filtering, then the instantaneous frequency of gear case rotating shaft is estimated by Viterbi optimum route search algorithm, recycling key signal estimation model carries out point-by-point integration to instantaneous frequency and obtains estimating key signal, finally combine angularly resampling and angle domain average technology carries out calculating order analysis to vibration acceleration signal, the order obtained based on instantaneous Frequency Estimation is composed.This order spectrogram fully can reflect the characteristic information of gearbox fault.The present invention combines the advantage of smoothed pseudo wigner ville disstribution, Viterbi optimum route search algorithm, angle domain average technology, calculating order analysis, effectively can carry out fault diagnosis to the gear case under variable speed operating condition.

Description

A kind of gear distress is without key phase angular domain average computation order analysis method
Technical field
The invention belongs to gear distress analytical approach field, be specifically related to a kind of gear distress without key phase angular domain average computation order analysis method.
Background technology
Traditional gear case vibration-testing and the analysis of test signal, test are carried out under constant rotational speed, constant duty, and gear case in real work runs under variable speed, variable working condition load often, the such as gear case of wind power generator group, Helicopter gearboxes, car deceleration device, crane etc.For crane, handling steel ladle is responsible for by crane, the process that during work, gear case experience raising speed promotes, reduction of speed is toppled over.At present, under gear case variable speed operating condition condition, in calculating order analysis method used, the key signal of rotating shaft be carry out gearbox fault calculate order analysis institute requisite.But, for some casing relative closures (rotating machinery is non-disconnectable) or can not, at axle head sensor installation gear case, be difficult to extract key signal.Particularly some gear reduction design time inconsiderate or assembling time compact conformation, cannot fitting key phase pulse transducer, to be difficult under variable rotate speed realization and to wait periodic sampling, to cause in practical work process to gear case vibration-testing and to the analysis of test signal, test is more difficult completes.Therefore; research effective method; for the gear case under variable rotate speed and rotary speed information cannot gather when monitor the running status of gear case; arrange rational predetermined spare part, shut down replacing construction; before maintenance of equipment is arranged in the generation of equipment failure accident, reduces burst shutdown event to greatest extent, give full play to equipment effectiveness; extension device serviceable life, significant to the security of national product, reliability.
Summary of the invention
The object of the present invention is to provide a kind of gear distress without key phase angular domain average computation order analysis method, the method can solve the difficulties carrying out fault diagnosis under variable rotate speed in the not collectable situation of key signal.
To achieve these goals, the technical solution used in the present invention is:
A kind of gear distress, without key phase angular domain average computation order analysis method, comprises the following steps:
1) instantaneous Frequency Estimation, specifically comprises following three steps:
1-1) gather vibration acceleration signal, low pass is carried out to the vibration acceleration signal gathered and protects phase filtering, obtain low frequency signal x (t) reflecting rotor turns;
The time-frequency distributions of low frequency signal x (t) 1-2) calculated by smoothed pseudo wigner ville disstribution, obtains time frequency distribution map;
In the time frequency distribution map obtained, 1-3) utilize Viterbi optimum route search algorithm to estimate the instantaneous frequency (InstantaneousFrequency, IF) of gear case rotating shaft; .
2) order based on instantaneous Frequency Estimation is composed, and specifically comprises following four steps:
2-1) utilizing key signal estimation model to step 1-3) instantaneous frequency that obtains carries out point-by-point integration, obtains estimating key signal;
2-2) utilize the estimation key signal that obtains to carry out angularly resampling to vibration acceleration signal, obtain the signal after angularly resampling;
2-3) angle domain average process is carried out to the signal after angularly resampling, obtain the signal after angle domain average process;
2-4) Fast Fourier Transform (FFT) (FFT) is carried out to the signal after angular domain average treatment, obtain reflecting that the order based on instantaneous Frequency Estimation of gearbox fault characteristic information is composed, gear distress is analyzed.
Described step 1-2) concrete operations be:
Step 1-1) the smoothed pseudo wigner ville disstribution SPW (t, ω) of low frequency signal x (t) that obtains is by shown in formula (1):
SPW ( t , ω ) = 1 2 π ∫ - ∞ + ∞ ∫ - ∞ + ∞ x ( t + 1 2 τ ) x * ( t - 1 2 τ ) h ( τ ) g ( u ) e - iτω dτdu - - - ( 1 )
Namely time frequency distribution map is obtained by formula (1); In formula (1), * represents conjugation, and t represents the time, and ω represents frequency, and h (τ) and g (u) is window function, and τ, u represent the shift factor of time shaft and frequency axis respectively.
Described window function comprises Gauss window and Hamming window.
Described step 1-3) concrete operations be:
Get the time series [n on the smoothed pseudo wigner ville disstribution of x (t) 1, n 2], time point n ∈ [n 1, n 2], the frequency values on each time point n is SPW n, SPW nby the array after descending sort for being designated as SPW ' n;
From 0 s, press increment 1 successively by order from big to small and distinguish assignment again, the time-frequency distributions sequence after assignment be designated as f (n, SPW ' n);
Instantaneous frequency is the path making formula (2) minimum:
ω ^ ( n ) = arg min [ Σ n = n 1 n 2 g ( ω n , ω n + 1 ) + Σ n = n 1 n 2 f ( n , SPW n ′ ) ] - - - ( 2 )
Wherein, g (ω n, ω n+1) be continuity in order to meet frequency change and the penalty factor that larger frequency change of adjusting the distance adds.
Described penalty factor g (ω n, ω n+1) such as formula shown in (3):
g ( ω n , ω n + 1 ) = 0 , | ω n - ω n + 1 | ≤ δ c ( | ω n - ω n + 1 | - δ ) , | ω n - ω n + 1 | > δ - - - ( 3 )
Wherein δ is the maximum expected value of the Instantaneous frequency variations of 2 continuity points.
Described step 2-1) concrete operations be:
Work under a rotating shaft is in variable rotate speed, rotation instantaneous frequency is f (t); Elapsed time t n, rotating shaft turns over n circle, and rotating shaft turns over angle such as formula shown in (4):
∫ 0 t n 2 πf ( t ) dt = 2 πn - - - ( 4 )
To f (t) point-by-point integration, record integrated value is the t of integer n, at t nthe key signal of a step is inserted in position, is namely the estimation key signal obtained by instantaneous frequency.
Described step 2-2) concrete operations be:
Timing node t when obtaining angularly θ resampling according to formula (5), carries out angularly resampling, obtains the signal after angularly resampling;
t = 1 2 b 2 [ 4 b 2 ( θ - b 0 ) + b 1 2 - b 1 ] - - - ( 5 )
B in formula (5) 0, b 1, b 2obtain according to formula (6);
Wherein
Described step 2-3) concrete operations be:
Burst after angularly resampling is designated as the burst y after angle domain average process is obtained by the formula formula (7) of angle domain average process n,
y n = 1 P Σ P = 0 P - 1 y ~ n + PN - - - ( 7 )
Wherein, P is the hop count signal averaging after angularly resampling be divided into, and N is each section of identical sampling number.
Relative to prior art, beneficial effect of the present invention is:
The gear distress that this patent proposes is without key phase angular domain average computation order analysis method method, the rotary speed information of rotating shaft can be extracted from vibration acceleration signal, effectively solve the problem being difficult to diagnosis without gear case fault under variable speed operating condition time key signal (key signal can not gather), namely can without the diagnosis completed during key signal and under variable rotate speed gearbox fault.The method combines smoothed pseudo wigner ville disstribution, Viterbi optimum route search algorithm, key signal estimation, angle domain average process and calculates the advantage of order analysis, has following differences in the significant advantage of classic method:
1) in instantaneous Frequency Estimation, present invention incorporates the time-frequency focusing that smoothed pseudo wigner ville disstribution is good, and Viterbi optimum route search algorithm follows the feature of free-running frequency change continuity principle, can estimate the instantaneous frequency of gear case rotating shaft effectively accurately.
2) compose in calculating at the order based on instantaneous Frequency Estimation, utilize key signal estimation model to carry out point-by-point integration to instantaneous frequency and estimated key signal accurately, then combine angularly resampling technique and angle domain average treatment technology and calculating order analysis is carried out to vibration acceleration signal, effectively resist the interference of noise and unrelated signal, obtained reflecting that the order based on instantaneous Frequency Estimation of gearbox fault characteristic information is composed.
3) on the whole, method provided by the invention combines smoothed pseudo wigner ville disstribution, Viterbi optimum route search algorithm, key signal estimation, angle domain average treatment technology and calculates the advantage of order analysis, utilize the present invention effectively can carry out fault diagnosis to the gear case that cannot gather key signal under variable rotate speed, improve the accuracy rate of fault diagnosis.
Accompanying drawing explanation
Fig. 1 (a) is the schematic diagram of Viterbi optimum route search algorithm sweep forward;
Fig. 1 (b) is the schematic diagram of Viterbi optimum route search algorithm sweep backward;
Fig. 1 (c) is the schematic diagram in the local optimum path of Viterbi optimum route search algorithm;
Fig. 1 (d) is the schematic diagram in the global optimum path of Viterbi optimum route search algorithm;
Fig. 2 is the schematic diagram of the rotational model of axle;
Fig. 3 (a) is emulation vibration signal for faster figure;
Fig. 3 (b) is the interpolation resampling figure to emulation vibration signal for faster;
Fig. 3 (c) is the estimation key signal figure of emulation vibration signal for faster;
Fig. 3 (d) carries out the signal graph after angularly resampling to estimation key signal;
Fig. 4 tests the vibration acceleration signal figure recorded;
Fig. 5 vibrates the signal graph after signal for faster low-pass filtering to gained;
Fig. 6 is the SPWVD time frequency distribution map of signal after low-pass filtering, and wherein a is that 1X turns frequently, and b is that 2X turns frequently;
Fig. 7 is the instantaneous frequency figure utilizing Viterbi optimum route search algorithm search to obtain;
Fig. 8 is the estimation key signal figure utilizing key signal estimation model to obtain;
Fig. 9 is the signal graph after angularly resampling;
Figure 10 is the order spectrogram based on instantaneous Frequency Estimation;
Figure 11 is that the epicyclic gearbox centre wheel flank of tooth evenly puts corrosion figure.
Embodiment
Below in conjunction with accompanying drawing, content of the present invention is described in further detail.
The invention discloses a kind of gear distress without key phase angular domain average computation order analysis method (Angle-domain-averageAndNon-keyphasorCalculatedOrderTrack ing, AN-COT), specifically implement according to the following steps:
1) instantaneous Frequency Estimation.Low pass is carried out to the vibration acceleration signal gathered and protects phase filtering, obtain the low frequency signal reflecting rotor turns.Then the time-frequency distributions of low frequency signal is calculated by smoothed pseudo wigner ville disstribution, the frequency division Butut obtained utilize Viterbi optimum route search algorithm estimate the instantaneous frequency (InstantaneousFrequency, IF) of gear case rotating shaft.Specifically comprise following 3 steps:
1-1) gather vibration acceleration signal, low pass is carried out to the vibration acceleration signal gathered and protects phase filtering, extract reflection turning axle and first the low pass guarantor phase filtering that input shaft turns signal is frequently comprised to signal, remove the interference of other orders in signal, the vibration signal (low frequency signal of reflection rotor turns) that to obtain with input shaft vibration signal and low power order thereof be fundamental component.The upper cut-off frequency of low-pass filter is according to the highest running speed setting in machine operation procedure.Protect phase filtering very important, if otherwise have the changing of the relative positions of phase place that IF can be caused to estimate and vibration signal on a timeline asynchronous after filtering simultaneously.
1-2) time frequency analysis: the time-frequency distributions of the low frequency signal calculated by smoothed pseudo wigner ville disstribution, obtains time frequency distribution map.
The Wigner-Ville distribution of low frequency signal x (t) is defined as,
WVD ( t , ω ) = 1 2 π ∫ - ∞ + ∞ x * ( t - τ / 2 ) x ( t + τ / 2 ) e - iτω dτ
Wigner-Ville has numerous advantages, particularly time-frequency locality, but can introduce cross term.In order to overcome the deficiency of cross term, can in time domain or time domain, frequency domain windowing simultaneously to overcome cross term, Here it is pseudo-Winger distribution (PWVD) and Smoothing Pseudo Winger distribution (SPWVD).Wigner-Ville distribution after time-domain windowed is called pseudo NM-algebra, is defined as follows:
PW ( t , ω ) = 1 2 π ∫ - ∞ + ∞ h ( τ ) x * ( t - 1 2 τ ) x ( t + 1 2 τ ) e - iτω dτ
Wherein, in formula, * represents conjugation.H (t) is window function.Conventional window function comprises Gauss window, Hamming window etc.
Then can at time and frequency both direction suppressing crossterms simultaneously in time domain, frequency domain windowing simultaneously.The Wigner-Ville distribution of this both direction windowing process is called smoothed pseudo wigner ville disstribution (SPWVD), and the smoothed pseudo wigner ville disstribution of x (t) is defined as follows:
SPW ( t , ω ) = 1 2 π ∫ - ∞ + ∞ ∫ - ∞ + ∞ x ( t + 1 2 τ ) x * ( t - 1 2 τ ) h ( τ ) g ( u ) e - iτω dτdu - - - ( 1 )
Wherein, t represents the time, and ω represents frequency, and g (u) is used to do level and smooth window function in frequency axis direction, and τ, u represent the shift factor of time shaft and frequency axis respectively.Namely time frequency distribution map is obtained by formula (1).
In the time frequency distribution map obtained, 1-3) utilize Viterbi optimum route search algorithm to estimate the instantaneous frequency (InstantaneousFrequency, IF) of gear case rotating shaft.
At 1-2) step tries to achieve in time frequency distribution map, adopts Viterbi algorithm search spectrum optimum time frequency path, peak.Spectrum blob detection just searches for the peak value in time-frequency distributions successively, and actual IF is as a continually varying physical quantity, and should consider that the change between two adjacent time points can not be excessive, Viterbi algorithm just in time considers this point.Time-frequency distributions proposes during Searching I F following two hypothesis,
(1) if on certain time point, the maximal value of time-frequency distributions is not the IF of signal, then the IF value of signal is that the probability of SPWVD n-th (n is less) individual maximal value is very large;
The IF change of (2) 2 continuous time points is not very greatly.
Only seeing that Article 1 is supposed, is in fact spectrum blob detection method.So the crucial part Article 2 hypothesis of Viterbi algorithm, what this hypothesis reflected is the continuity that in rotating machinery variable parameter operation engineering, free-running frequency changes, and can not occur very large frequency agility within the extremely short time.
Viterbi algorithm ascribes an extreme-value problem to, namely gets the time series [n in SPWVD distribution 1, n 2] between time-frequency distributions analyze, time point n ∈ [n 1, n 2], the frequency values on each time point n is SPW n, namely
SPW n=[SPW(n,ω 1),SPW(n,ω 2),…,SPW(n,ω M)]
M is that the frequency comprised is counted.
SPW nbe SPW ' by the array define after descending sort n, then f (n, SPW ' n) be defined as the time-frequency distributions sequence after assignment,
f(n,SPW′ n)=j-1(j=0,1,…,M)
The meaning of above formula be exactly the value at all Frequency point places on SPWVD distribution time point n by order from big to small from 0 successively assignment by increment 1 assignment respectively.
Instantaneous frequency is the path making formula (2) minimum,
ω ^ ( n ) = arg min [ Σ n = n 1 n 2 g ( ω n , ω n + 1 ) + Σ n = n 1 n 2 f ( n , SPW n ′ ) ] - - - ( 2 )
In formula, f (n, SPW ' n) be construct according to hypothesis (1), g (ω n, ω n+1) according to hypothesis (2) structure.G (ω n, ω n+1) be continuity in order to meet frequency change and penalty factor that larger frequency change of adjusting the distance adds, can ensure that IF path meets above hypothesis (2).G (ω n, ω n+1) be defined as follows:
g ( ω n , ω n + 1 ) = 0 , | ω n - ω n + 1 | ≤ δ c ( | ω n - ω n + 1 | - δ ) , | ω n - ω n + 1 | > δ - - - ( 3 )
In formula, the optimal selection of δ is the maximum expected value that 2 continuity point IF change, and the cost function namely changing (| x-y|≤δ) for less IF is 0, chooses less δ (as δ=3) in an experiment and can obtain good effect.
With δ=1, c=5 is the implementation procedure that example illustrates Viterbi optimum route search algorithm, and get certain local time frequency plane and analyze, time frequency plane has 3 time points, 8 Frequency points, the concrete steps of Viterbi optimum route search algorithm are:
(1) according to formula f (n, SPW ' n)=j-1 (j=0,1 ..., M) and calculate each point (n, ω i) corresponding f (n, ω i) functional value, determine the local optimum path of all frequencies of each time point n and all frequencies of its previous time point n-1 and record, choose f (n, SPW ' n)+g (ω n, ω n-1) minimum path as the local optimum path (left arrow line segment) between two time points, as shown in Fig. 1 (a);
(2) the local optimum path of all Frequency points in like manner, determining each time n and all Frequency points of a time point n+1 thereafter is also recorded (to the right arrow line segment), as shown in Fig. 1 (b);
(3) according to formula (3), with the local optimum path between two time points for benchmark, first along the functional value in all paths of local optimum path computing determined by step (1) straight line, the minimum path of record functional value is as shown in arrow left in Fig. 1 (c); Then the local optimum path searched for backward is determined after the same method, as shown in arrow to the right in Fig. 1 (c);
(4) compare to get in Fig. 1 (c) forward, the size of searching route functional value backward, the path selecting minimum function value corresponding is optimal path, as shown in Fig. 1 (d).
In addition, due to the edge effect of time-frequency distributions, there is very big error in the rim value that IF estimates, can remove some data of the right and left in time frequency distribution map.Simultaneously for the feature of multi-components composition in actual condition middle gear case vibration signal, need to go out on the basis of an IF in above-mentioned Viterbi algorithm search, estimate other IF.Suppose now to estimate it is the signal that respective signal is the strongest.Will value zero setting in scope, forms new time-frequency representation, and ξ is the constant of a setting.Above Viterbi searching algorithm step is repeated to new time-frequency representation, until estimated all to estimate by the IF of component each in signal.
2) order based on instantaneous Frequency Estimation is composed.Utilize key signal estimation model to carry out point-by-point integration to instantaneous frequency and estimated key signal accurately.Then carry out calculating order analysis in conjunction with angle domain average technology to vibration acceleration signal, the order obtained based on instantaneous Frequency Estimation is composed.Specifically comprise following four steps:
2-1) utilizing key signal estimation model to step 1-3) instantaneous frequency that obtains carries out point-by-point integration, obtains estimating key signal.On the basis obtaining IF, in order to computation order analysis, need to estimate key signal with IF.Suppose that a rotating shaft works under being in variable rotate speed, rotating speed is ω (t), as shown in Figure 2.The pass of ω (t) and rotation instantaneous frequency f (t) is ω (t)=2 π f (t).Getting the angle that certain time period in Fig. 2 turns over is Δ θ, and Δ θ can be expressed as:
Δθ=ω(t)dt=2πf(t)dt
Now suppose to add up from the measurement starting point t of a point as the time, every a bit of Δ θ, then elapsed time t 1, rotor forwards again a point to, and the angle that rotor turns over is 2 π.Obtain following equation thus:
∫ 0 t 1 2 πf ( t ) dt = 2 π
If through t n, rotating shaft turns over n circle, and rotating shaft turns over angle such as formula shown in (4):
∫ 0 t n 2 πf ( t ) dt = 2 πn - - - ( 4 )
In formula (4), t nit is the unknown quantity needing to solve.To f (t) point-by-point integration, record integrated value is the t of integer n, at t nthe key signal of a step is inserted in position, is namely the estimation key signal obtained by IF.
2-2) utilize the estimation key signal that obtains to carry out angularly resampling to vibration acceleration signal, obtain the signal after angularly resampling.
Utilize spindle keyway to mark or a pulse can only be collected in each turn of reflecting piece, angularly resampling needs the work done to be exactly the pulse utilizing a series of mark each turn of 2 π radians gathered, utilize corresponding computing method interpolation to be divided into N section, the radian of each section is for calculating the angularly resampling moment, suppose that rotating shaft rotating shaft within the time period that each is little is made even angular acceleration and run.
The uniform angular velocity rotation equation of shaft rotary corner θ can be described as:
θ(t)=b 0+b 1t+b 2t 2
In formula: b 0, b 1, b 2for unknown parameter, t is time point.
Get the time point t of three continuous impulses that tacho-pulse sensor gathers 1, t 2, t 3, with t 1as the start index of angle step, then t 1and t 2, t 1and t 3angle step between time point is respectively following equation can be listed:
B can be solved by equation 0, b 1, b 2, shown in (6):
Meanwhile, by equation θ (t)=b 0+ b 1t+b 2t 2solve t:
t = 1 2 b 2 [ 4 b 2 ( θ - b 0 ) + b 1 2 - b 1 ] - - - ( 5 )
Formula (6) is substituted into timing node t when formula (5) can obtain angularly θ resampling, then carry out angularly resampling, obtain the signal after angularly resampling.The schematic diagram of an angularly resampling is given shown in Fig. 3.In the drawings, in every two key signals collected every individual insertion resample points (namely often turning sampling 8 points), these time points inserted, from time shaft, are not constant durations, but in angular domain read fortune every with equally spaced
2-3) angle domain average process is carried out to the signal after angularly resampling, obtain the signal after angle domain average process.
Burst after angularly resampling is defined as identical with Synchronous time average thought, the algorithmic formula of angular domain synchronized averaging is formula (7), obtains the burst y after angle domain average process by formula (7) n.
y n = 1 P Σ P = 0 P - 1 y ~ n + PN - - - ( 7 )
Wherein,
P: the hop count that the signal averaging after angularly resampling is divided into;
N: each section identical sampling number.
Periodic mechanical signal sequences y (n) is had for what obtained by angular domain resampling in rotating machinery or reciprocating machine operational process, suppose that y (n) is made up of periodic signal f (n) and random noise s (n), that is:
x(n)=f(n)+s(n)
Go to intercept signal x (n) with the cycle T of f (n), intercept altogether M section, be then added, due to the irrelevance of random noise, can obtain:
x ( n i ) = Mf ( n i ) + M s ( n i )
Again to x (n i) be averaging, obtain average after output signal y (n i):
y ( n i ) = f ( n i ) + 1 M s ( n i )
Now random noise is the random noise in original input signal x (t) thus signal to noise ratio (S/N ratio) improves doubly.
Therefore, the wave filter of angular domain synchronized sampling is similar with the filter characteristic of time domain synchronized sampling, can stress release treatment and unrelated signal effectively.
2-4) Fast Fourier Transform (FFT) (FFT) is carried out to the signal after angular domain average treatment, namely obtain reflecting that the order based on instantaneous Frequency Estimation of gearbox fault characteristic information is composed.
Case study on implementation
The validity of this invention of signal authentication is moved with certain planetary reduction gear testing table institute vibration measuring.This testing table comprises motor, speed reduction unit, generator, control desk four part.Planetary reduction gear is a secondary precision speed reduction device, and reduction gear ratio is 32, and the structural parameters of planetary reduction gear are as shown in table 1 below.
Table 1 star reducer structure parameter
According to the parameter of table 1 and the meshing frequency of planetary gear meshing frequency computing formula calculating first order transmission sun gear and planetary gear.Planetary gear meshing frequency computing formula is as follows:
F engagement=(n 1-n h) × z 1
In formula, n 1: sun gear rotating speed; n h: planet carrier rotating speed; z 1: the sun gear number of teeth.
F is obtained by above formula engagement=10.5n 1.10.5 is important parameters, contrasts previously described order concept, J n=10.5 is planetary gear order parameters when breaking down.
Manual adjustments control box, makes the speed of input motor be a moderating process, in time domain, gathers vibration signal respectively as shown in Figure 4.Sample frequency is 10240Hz, meets sample frequency and arranges requirement.
Setting maximum analysis order according to feature order is 40, and 80 points of namely angularly sampling within the time that axle often circles, the sampling angle of equiangular sampling is 2 π/80.
Carry out to signal the zero phase filtering that cutoff frequency is 50Hz, filtered signal as shown in Figure 5.
SPWVD is utilized to ask for the time-frequency distributions of signal.Getting window function is Hamming window, and length of window is 1/10 of signal length, after obtaining time-frequency distributions, then removes the time-frequency distributions of 1s duration at time-frequency distributions two ends, to reduce end effect.The time-frequency distributions obtained is illustrated in fig. 6 shown below.Can see and time-frequency figure obviously has input shaft to turn 1 times frequently, the time-frequency distributions of 2 times.
Utilize Viterbi algorithm, time-frequency figure searches out 1 times of instantaneous turn of turning frequently respectively to be estimated in as Fig. 7 shown in heavy line frequently, 2 times turn instantaneous turn frequently and frequently estimate that (2 times of optimal paths turning IF frequently of extraction are divided by 2, then obtain 1 times of IF estimation turned frequently) as shown in Fig. 7 dotted line, true instantaneous frequency profile is as shown in Fig. 7 fine line.As can be seen from the figure the result that estimates of instantaneous turn of frequency and true instantaneous frequency profile can extraordinaryly coincide.
First turn IF frequently with 1 times to estimate, to conjunction key phase algorithm for estimating, namely can obtain the key signal estimated, as shown in Figure 8.
Utilize the key signal estimated, interpolation resampling is carried out to original signal and can obtain equal angular resampling signal, as shown in Figure 9.
Angle domain average process is carried out to the signal after angular domain resampling, and carries out FFT conversion., obtain order spectrum as shown in Figure 10.Spectrogram can be seen significantly relevant with the order 10.5 in the transmission of epicyclic gearbox first order spectrum peak, as figure J nshown in mark, both sides is with obvious fundamental frequency characteristic edge frequency band simultaneously, and namely sideband is rotating speed order 1.At 2 times of meshing frequency 2J nplace, the frequency conversion band being spaced apart fundamental frequency is high and narrow, belongs to the even pitting fault feature of typical gear.Open this epicyclic gearbox, each flank of tooth of discovery first order transmission centre wheel has the even point corrosion of the flank of tooth shown in Figure 11 (in oval circle).
The AN-COT algorithm that this patent proposes, under the strong phase signals of nothing, the Accurate Diagnosis flank of tooth of gear case under the condition under variable rotate speed even pitting fault feature, improves the accuracy of fault diagnosis, ensures effectively reliably for the Fault Diagnosis of Gear Case under variable speed condition provides.

Claims (8)

1. gear distress is without a key phase angular domain average computation order analysis method, it is characterized in that, comprises the following steps:
1) instantaneous Frequency Estimation:
1-1) gather vibration acceleration signal, low pass is carried out to the vibration acceleration signal gathered and protects phase filtering, obtain low frequency signal x (t) reflecting rotor turns;
The time-frequency distributions of low frequency signal x (t) 1-2) calculated by smoothed pseudo wigner ville disstribution, obtains time frequency distribution map;
In the time frequency distribution map obtained, 1-3) utilize Viterbi optimum route search algorithm to estimate the instantaneous frequency of gear case rotating shaft;
2) order based on instantaneous Frequency Estimation is composed:
2-1) utilizing key signal estimation model to step 1-3) instantaneous frequency that obtains carries out point-by-point integration, obtains estimating key signal;
2-2) utilize the estimation key signal that obtains to carry out angularly resampling to vibration acceleration signal, obtain the signal after angularly resampling;
2-3) angle domain average process is carried out to the signal after angularly resampling, obtain the signal after angle domain average process;
2-4) Fast Fourier Transform (FFT) is carried out to the signal after angular domain average treatment, obtain reflecting that the order based on instantaneous Frequency Estimation of gearbox fault characteristic information is composed, gear distress is analyzed.
2. gear distress according to claim 1 is without key phase angular domain average computation order analysis method, it is characterized in that: described step 1-2) concrete operations be:
Step 1-1) the smoothed pseudo wigner ville disstribution SPW (t, ω) of low frequency signal x (t) that obtains is by shown in formula (1):
S P W ( t , ω ) = 1 2 π ∫ - ∞ + ∞ ∫ - ∞ + ∞ x ( t + 1 2 τ ) x * ( t - 1 2 τ ) h ( τ ) g ( u ) e - i τ ω d τ d u - - - ( 1 )
Namely time frequency distribution map is obtained by formula (1); In formula (1), * represents conjugation, and t represents the time, and ω represents frequency, and h (τ) and g (u) is window function, and τ, u represent the shift factor of time shaft and frequency axis respectively.
3. gear distress according to claim 2 is without key phase angular domain average computation order analysis method, it is characterized in that: described window function comprises Gauss window and Hamming window.
4. according to the gear distress in claim 1-3 described in any one without key phase angular domain average computation order analysis method, it is characterized in that: described step 1-3) concrete operations be:
Get the time series [n on the smoothed pseudo wigner ville disstribution of x (t) 1, n 2], time point n ∈ [n 1, n 2], the frequency values on each time point n is SPW n, SPW nby the array after descending sort for being designated as SPW ' n;
From 0 s, press increment 1 successively by order from big to small and distinguish assignment again, the time-frequency distributions sequence after assignment be designated as f (n, SPW ' n);
Instantaneous frequency is the path making formula (2) minimum:
ω ^ ( n ) = arg min [ Σ n = n 1 n 2 g ( ω n , ω n + 1 ) + Σ n = n 1 n 2 f ( n , SPW n ′ ) ] - - - ( 2 )
Wherein, g (ω n, ω n+1) be continuity in order to meet frequency change and the penalty factor that larger frequency change of adjusting the distance adds.
5. gear distress according to claim 4 is without key phase angular domain average computation order analysis method, it is characterized in that: described penalty factor g (ω n, ω n+1) such as formula shown in (3):
g ( ω n , ω n + 1 ) = 0 , | ω n - ω n + 1 | ≤ δ c ( | ω n - ω n + 1 | - δ ) , | ω n - ω n + 1 | > δ - - - ( 3 )
Wherein δ is the maximum expected value of the Instantaneous frequency variations of 2 continuity points.
6. gear distress according to claim 4 is without key phase angular domain average computation order analysis method, it is characterized in that: described step 2-1) concrete operations be:
Work under a rotating shaft is in variable rotate speed, rotation instantaneous frequency is f (t); Elapsed time t n, rotating shaft turns over n circle, and rotating shaft turns over angle such as formula shown in (4):
∫ 0 t n 2 π f ( t ) d t = 2 π n - - - ( 4 )
To f (t) point-by-point integration, record integrated value is the t of integer n, at t nthe key signal of a step is inserted in position, is namely the estimation key signal obtained by instantaneous frequency.
7. gear distress according to claim 4 is without key phase angular domain average computation order analysis method, it is characterized in that: described step 2-2) concrete operations be:
Timing node t when obtaining angularly θ resampling according to formula (5), carries out angularly resampling, obtains the signal after angularly resampling;
t = 1 2 b 2 [ 4 b 2 ( θ - b 0 ) + b 1 2 - b 1 ] - - - ( 5 )
B in formula (5) 0, b 1, b 2obtain according to formula (6);
Wherein
8. gear distress according to claim 4 is without key phase angular domain average computation order analysis method, it is characterized in that: described step 2-3) concrete operations be:
Burst after angularly resampling is designated as i=1,2,3 ..., N, obtains the burst y after angle domain average process by the formula (7) of angle domain average process n,
y n = 1 P Σ P = 0 P - 1 y ~ n + P N - - - ( 7 )
Wherein, P is the hop count signal averaging after angularly resampling be divided into, and N is each section of identical sampling number.
CN201310416504.7A 2013-09-12 2013-09-12 A kind of gear distress is without key phase angular domain average computation order analysis method CN103499443B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201310416504.7A CN103499443B (en) 2013-09-12 2013-09-12 A kind of gear distress is without key phase angular domain average computation order analysis method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201310416504.7A CN103499443B (en) 2013-09-12 2013-09-12 A kind of gear distress is without key phase angular domain average computation order analysis method

Publications (2)

Publication Number Publication Date
CN103499443A CN103499443A (en) 2014-01-08
CN103499443B true CN103499443B (en) 2016-01-20

Family

ID=49864674

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201310416504.7A CN103499443B (en) 2013-09-12 2013-09-12 A kind of gear distress is without key phase angular domain average computation order analysis method

Country Status (1)

Country Link
CN (1) CN103499443B (en)

Families Citing this family (17)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104266747B (en) * 2014-06-09 2017-10-24 中能电力科技开发有限公司 A kind of method for diagnosing faults based on vibration signal order analysis
CN104077474A (en) * 2014-06-23 2014-10-01 华南理工大学 Meshing frequency and spectrum correction technology based wind power gear box order tracking method
CN104655380B (en) * 2015-03-16 2017-10-24 北京六合智汇技术有限责任公司 A kind of rotating machinery fault signature extracting method
CN104819841B (en) * 2015-05-05 2017-04-19 西安交通大学 Built-in-coding-information-based single sensing flexible angle-domain averaging method
CN105277362B (en) * 2015-11-23 2017-09-12 西安交通大学 Gear distress detection method based on encoder multidigit angular signal
CN105547698B (en) * 2015-12-31 2017-11-14 新疆金风科技股份有限公司 The method for diagnosing faults and device of rolling bearing
US10094743B2 (en) * 2016-03-14 2018-10-09 Epro Gmbh Order analysis system
CN106644467B (en) * 2016-12-27 2019-05-14 华南理工大学 A kind of gear-box non-stationary signal fault signature extracting method
CN107292067B (en) * 2017-08-17 2021-02-02 湖南纬拓信息科技有限公司 Gear fault diagnosis method based on compressed sensing and bispectrum analysis
CN107907324A (en) * 2017-10-17 2018-04-13 北京信息科技大学 A kind of Fault Diagnosis of Gear Case method composed based on DTCWT and order
CN107941510B (en) * 2017-10-19 2019-07-19 西安交通大学 Extracting method based on the angularly Rolling Bearing Fault Character of dual sampling
CN108225764A (en) * 2017-12-05 2018-06-29 昆明理工大学 It is a kind of based on the high-precision of envelope extraction without key signal Order Tracking and system
CN108362942B (en) * 2018-02-11 2020-11-20 中国铁道科学研究院集团有限公司 Time-frequency spectrum obtaining method and device for multi-component signals
CN108871742A (en) * 2018-05-03 2018-11-23 西安交通大学 A kind of improved no key phase fault feature order extracting method
CN109813546B (en) * 2019-03-20 2020-09-15 苏州微著设备诊断技术有限公司 Off-line detection method for abnormal knocking sound of gear box
CN110686768A (en) * 2019-10-17 2020-01-14 昆明理工大学 Improved rotating machinery nonstationary vibration signal calculation order ratio analysis method
CN110907162A (en) * 2019-12-13 2020-03-24 北京天泽智云科技有限公司 Rotating machinery fault feature extraction method without tachometer under variable rotating speed

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102175393A (en) * 2011-01-04 2011-09-07 广东电网公司电力科学研究院 Imbalance phase estimation method based on procession decomposition technology
US8024120B2 (en) * 2008-05-16 2011-09-20 Turner Larry A Complex phase locked loop
FR2967541A1 (en) * 2010-11-17 2012-05-18 Centre Nat Etd Spatiales Oversampled digital signal i.e. split phase-level/phase modulation signal, demodulating method, involves calculating frequency estimation from channel estimation, extracting data, and compensating carrier phase with channel estimation
CN102998110A (en) * 2012-11-29 2013-03-27 西安交通大学 Rotary machine fault characteristic extraction method based on order-holospectrum principle
CN103196547A (en) * 2013-03-11 2013-07-10 安徽新力电业科技咨询有限责任公司 Method achieving rotary machinery vibration signal synchronization order ratio tracing analysis
CN103234627A (en) * 2013-04-17 2013-08-07 国家电网公司 Complete alternation synchronous sampling and analyzing method for rotating machinery vibration signals
CN103278235A (en) * 2013-06-03 2013-09-04 合肥伟博测控科技有限公司 Novel transient oscillation signal angular domain order tracking sampling and analytical method

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8024120B2 (en) * 2008-05-16 2011-09-20 Turner Larry A Complex phase locked loop
FR2967541A1 (en) * 2010-11-17 2012-05-18 Centre Nat Etd Spatiales Oversampled digital signal i.e. split phase-level/phase modulation signal, demodulating method, involves calculating frequency estimation from channel estimation, extracting data, and compensating carrier phase with channel estimation
CN102175393A (en) * 2011-01-04 2011-09-07 广东电网公司电力科学研究院 Imbalance phase estimation method based on procession decomposition technology
CN102998110A (en) * 2012-11-29 2013-03-27 西安交通大学 Rotary machine fault characteristic extraction method based on order-holospectrum principle
CN103196547A (en) * 2013-03-11 2013-07-10 安徽新力电业科技咨询有限责任公司 Method achieving rotary machinery vibration signal synchronization order ratio tracing analysis
CN103234627A (en) * 2013-04-17 2013-08-07 国家电网公司 Complete alternation synchronous sampling and analyzing method for rotating machinery vibration signals
CN103278235A (en) * 2013-06-03 2013-09-04 合肥伟博测控科技有限公司 Novel transient oscillation signal angular domain order tracking sampling and analytical method

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
无转速计旋转机械升降速振动信号零相位阶比跟踪滤波;郭瑜等;《机械工程学报》;20040331;第40卷(第3期);50-54 *

Also Published As

Publication number Publication date
CN103499443A (en) 2014-01-08

Similar Documents

Publication Publication Date Title
Feng et al. Iterative generalized synchrosqueezing transform for fault diagnosis of wind turbine planetary gearbox under nonstationary conditions
Abboud et al. The spectral analysis of cyclo-non-stationary signals
Borghesani et al. A new procedure for using envelope analysis for rolling element bearing diagnostics in variable operating conditions
Chow et al. Induction machine fault diagnostic analysis with wavelet technique
Chen et al. Generator bearing fault diagnosis for wind turbine via empirical wavelet transform using measured vibration signals
US20160033580A1 (en) Detecting Faults in Turbine Generators
Wilkinson et al. Condition monitoring of generators & other subassemblies in wind turbine drive trains
Li et al. New procedure for gear fault detection and diagnosis using instantaneous angular speed
CN102937668B (en) Electric system low-frequency oscillation detection method
CN102901929B (en) For calculating device and the battery impedance measuring system of cell impedance
CN101221066B (en) Engineering nonlinear vibration detecting method
CN103321854B (en) A kind of vibration control method for wind generator set tower
Li et al. A new noise-controlled second-order enhanced stochastic resonance method with its application in wind turbine drivetrain fault diagnosis
Wang et al. Rolling element bearing fault diagnosis via fault characteristic order (FCO) analysis
Zhao et al. A tacho-less order tracking technique for large speed variations
CN102435844B (en) Sinusoidal signal phasor calculating method being independent of frequency
He et al. Multi-scale stochastic resonance spectrogram for fault diagnosis of rolling element bearings
Tsao et al. An insight concept to select appropriate IMFs for envelope analysis of bearing fault diagnosis
Sawalhi et al. Gear parameter identification in a wind turbine gearbox using vibration signals
Blunt et al. Detection of a fatigue crack in a UH-60A planet gear carrier using vibration analysis
CN102483368B (en) Method for detecting structural defect in mechanical assembly including rotary member
CN102156043B (en) Online state monitoring and fault diagnosis system of wind generator set
Feng et al. Fault diagnosis of wind turbine planetary gearbox under nonstationary conditions via adaptive optimal kernel time–frequency analysis
Barszcz et al. Application of spectral kurtosis for detection of a tooth crack in the planetary gear of a wind turbine
Feng et al. Joint envelope and frequency order spectrum analysis based on iterative generalized demodulation for planetary gearbox fault diagnosis under nonstationary conditions

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant