CN113984282B - Rotor dynamic balance fault feature extraction method under keyless phase-change rotating speed - Google Patents

Rotor dynamic balance fault feature extraction method under keyless phase-change rotating speed Download PDF

Info

Publication number
CN113984282B
CN113984282B CN202111188122.4A CN202111188122A CN113984282B CN 113984282 B CN113984282 B CN 113984282B CN 202111188122 A CN202111188122 A CN 202111188122A CN 113984282 B CN113984282 B CN 113984282B
Authority
CN
China
Prior art keywords
frequency
time
phase
spectrum
discrete
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.)
Active
Application number
CN202111188122.4A
Other languages
Chinese (zh)
Other versions
CN113984282A (en
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.)
Xiamen University
Original Assignee
Xiamen University
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 Xiamen University filed Critical Xiamen University
Priority to CN202111188122.4A priority Critical patent/CN113984282B/en
Publication of CN113984282A publication Critical patent/CN113984282A/en
Application granted granted Critical
Publication of CN113984282B publication Critical patent/CN113984282B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01MTESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M1/00Testing static or dynamic balance of machines or structures
    • G01M1/14Determining unbalance
    • G01M1/16Determining unbalance by oscillating or rotating the body to be tested
    • G01M1/22Determining unbalance by oscillating or rotating the body to be tested and converting vibrations due to unbalance into electric variables

Abstract

A method for extracting dynamic balance fault characteristics of a rotor under a keyless phase-change rotating speed comprises the following steps: carrying out band-pass filtering and noise reduction processing on the acquired vibration signals; performing time-frequency analysis on the vibration signals subjected to the noise reduction processing through short-time Fourier transform to obtain time-frequency energy distribution; calculating by using a Sobel operator to obtain an energy gradient, and searching a target ridge line with the minimum energy gradient by using a dynamic programming algorithm to adjust the discontinuity of ridge line change caused by noise and realize the accurate estimation of the instantaneous frequency conversion of the variable-speed vibration signal; and performing order tracking calculation and discrete spectrum correction on the vibration signal subjected to noise reduction processing according to the extracted instantaneous frequency conversion, and extracting the amplitude and phase corresponding to a first order spectrum. The method can overcome the interference of adjacent frequency and strong noise, realize the accurate extraction of the dynamic balance fault characteristics of the variable-speed rotor under the condition of no key phase signal, save the hardware cost and improve the dynamic balance efficiency.

Description

Rotor dynamic balance fault feature extraction method under keyless phase-change rotating speed
Technical Field
The invention relates to the field of rotor dynamic balance fault diagnosis and signal analysis and processing under variable rotating speed, in particular to a rotor dynamic balance fault feature extraction method under keyless phase variable rotating speed.
Background
With the rapid development of production technology and equipment manufacturing, the mechanical equipment structure becomes more and more complex, and in the rotating parts of these complex equipment, the rotor is mostly used as the main working part. Rotor imbalance is becoming a leading cause of rotary machine failure. The vibration signal is the most important information source in the monitoring and diagnosis process of the complex equipment, can directly reflect the running state of the equipment, contains rich fault information, and is usually applied to the diagnosis of the imbalance fault of the rotor. The traditional method for diagnosing the unbalance fault of the rotor is based on a steady state and established on the basis of a steady signal, however, due to the inherent design limitation of a motor, torsional vibration, assembly fault, time-varying torque and the like, the rotor usually runs under the non-steady working condition of variable rotating speed and variable load in the actual running process, at the moment, the rotating speed of the rotor has periodic or quasi-periodic fluctuation around the fundamental frequency rotating speed and is difficult to keep constant, if the signal is artificially assumed to be the steady signal for processing, a serious frequency fuzzy phenomenon is generated as a result, and certain difficulty is brought to the traditional fault diagnosis. How to extract the fundamental frequency signal from the signal measured by the sensor under the condition of variable rotating speed through a signal analysis processing method, and accurately calculate the unbalance amount of the rotor, which becomes the premise and the basis of effective diagnosis of the dynamic balance of the rotor.
In order to solve the problem of difficulty in extracting fault features at variable rotation speed, an angular domain sampling theory and an order tracking method established on the basis of the angular domain sampling theory are developed. Order tracking methods include tachometer-based order tracking methods and keyless phase order tracking methods. The order tracking method based on the tachometer needs additional hardware equipment (an optical encoder, the tachometer, a key phaser and the like), so that the hardware cost is increased, difficulty is brought to installation and debugging, and meanwhile, the method has the defect that the calculation precision is seriously limited by the resolution of the tachometer. For this reason, the keyless phase order tracking method is gradually applied.
The chinese invention patent 201410052317.X proposes a method for extracting a rotor dynamic unbalance signal when the rotor rotates at a non-stationary speed. The method is based on a key phase pulse signal and an interpolation method to obtain an instantaneous frequency estimation value in each revolution, and then rotation speed independence processing is carried out on a vibration signal according to the instantaneous frequency, so that extraction of the amplitude and the phase of an unbalanced signal under a non-stable rotation speed is realized. The method is still applicable to the case of external key phase signal reference. The chinese invention patent 201911285855.2 proposes a method for extracting fault characteristics of a rotating machine without a tachometer under variable rotation speed, which is based on the assumption that the rotation speed of a signal does not change in a short time (such as 1-4 cycles), and realizes the extraction of the instantaneous rotation frequency of the rotating machine with variable rotation speed by means of a matching tracking method, but is not applicable to the case of severe fluctuation of the rotation speed, and even seriously affects the extraction accuracy of the instantaneous rotation frequency. The invention patent 202010401814.1 in China provides a method and a system for diagnosing faults of a transmission shaft under the working condition of rotating speed fluctuation. According to the method, the central main frequency of a signal is extracted from the highest point of time-frequency energy by traversing the time axis in the time-frequency domain image under the condition without a tachometer, so that the time-frequency ridge line is obtained, the instantaneous frequency conversion of the transmission shaft under the fluctuation of the rotating speed is extracted, and the accuracy of the instantaneous frequency conversion is influenced by noise and adjacent frequency interference.
Disclosure of Invention
The invention mainly aims to overcome the defects in the prior art, provides a method for extracting the dynamic balance fault characteristics of the rotor under the condition of the keyless phase change rotating speed, can effectively avoid the problem of 'spectrum ambiguity' caused by adopting a traditional method based on stable signal spectrum analysis, has high extracted instantaneous frequency conversion accuracy of the rotor, has strong robustness under the conditions of strong background noise and adjacent frequency interference, can realize the accurate extraction of the dynamic balance fault characteristics of the variable rotating speed rotor under the condition of the keyless phase signal, saves the hardware cost and improves the dynamic balance efficiency.
The invention adopts the following technical scheme:
a rotor dynamic balance fault feature extraction method under a keyless phase-change rotating speed is characterized by comprising the following steps:
step (1), converting the frequency fcSetting a bandwidth as a center, and then carrying out noise reduction processing on an original vibration signal by adopting a band-pass filter with zero phase offset characteristic;
selecting proper window length and overlapping number, carrying out time-frequency analysis on the vibration signals x (k) subjected to noise reduction treatment by adopting short-time Fourier transform to obtain a time-frequency spectrogram, carrying out logarithm treatment to obtain a time-frequency energy distribution graph, wherein N is the signal length, and k is 1,2, …;
And (3) calculating the obtained time-frequency energy distribution graph through a Sobel operator to obtain an energy gradient, carrying out normalization processing, and normalizingThe energy gradient map after the chemical treatment comprises a plurality of ridge lines, namely corresponding spectral lines of different orders; the extraction problem of the ridge line is changed into an optimization problem, and the frequency conversion f is selectedcThe corresponding points are used as the starting points of the target ridge lines, the target ridge lines with the minimum energy gradient are searched by adopting a dynamic programming algorithm, and the frequency corresponding to each time point of the target ridge lines is the instantaneous frequency conversion of the variable-speed vibration signals in the whole time interval;
and (4) carrying out order tracking calculation on the vibration signal subjected to noise reduction processing in the step (1) according to the extracted instantaneous frequency conversion to obtain an order spectrum, carrying out discrete spectrum correction on the first order spectrum in the order spectrum, and extracting to obtain the amplitude and the phase corresponding to the first order spectrum.
Preferably, the band-pass filter selected in step (1) is required to have a zero phase shift characteristic, for example, a harmonic wavelet filter with a "box" characteristic is adopted, and since the phase difference between the real part and the imaginary part of the harmonic wavelet is 90 °, it means that the phase of the signal remains unchanged before and after the harmonic wavelet transform, and the accuracy of phase extraction in the dynamic balance fault characteristic of the rotor is ensured.
Preferably, the specific steps of obtaining the time-frequency energy distribution map in the step (2) are as follows:
and (2.1) performing discrete short-time Fourier transform on the discrete vibration signal subjected to the noise reduction treatment, wherein the discrete short-time Fourier transform comprises the following steps:
Figure GDA0003654260160000031
where S (i Δ t, f) represents the transformed spectrum, f represents the frequency, Δ t is 1/fsRepresenting the sampling time interval, fsRepresenting a sampling frequency, i represents a time series, k ═ 1,2, …, N represents a sample series, N represents a signal length, x represents a complex conjugate, x (k) represents a discrete form of a noise-reduced vibration signal, and g (k) represents a discrete form of a window function g (t);
and (2.2) converting the time frequency spectrum after the discrete short-time Fourier transform into energy distribution consisting of time discrete points and frequency discrete points, wherein the conversion formula is as follows:
P(i△t,f)=10*log10(S(i△t,f)) (2)。
preferably, the specific steps of obtaining the instantaneous frequency f (i Δ t) of the variable-speed vibration signal in the step (3) are as follows:
step (3.1), regarding the time-frequency energy distribution P (i Δ t, f) as an energy function e (x, y) corresponding to the coordinate point (x, y), calculating the energy gradient of the time-frequency spectrum by using a Sobel operator, and performing normalization processing, wherein the energy gradient G (e (x, y)) is calculated as follows:
G(e(x,y))=abs(filter2(hx,P(i△t,f)))+abs(filter2(hy,P(i△t,f))) (3)
in the formula, hx and hy are Sobel operators in the horizontal direction and the vertical direction respectively,
Figure GDA0003654260160000032
Figure GDA0003654260160000033
filter2(·) represents a two-dimensional digital filter function;
And (3.2) because the local maximum value and the real maximum value in the time frequency spectrum caused by the noise shift, the frequency jump is large and the energy gradient is increased, the ridge line formed by the local maximum values in the time frequency spectrum is discontinuously changed, so that the target ridge line is extracted by searching the minimum value point of the energy gradient to adjust the ridge line discontinuity caused by the noise and the adjacent frequency interference, and the optimal target ridge line is extracted by adopting an optimal decision dynamic programming algorithm. Defining the matrix M (x, y), initializing M (x,1) to G (e (x,1)), the remaining M (x, y) values starting at the second column and ending at the last column are calculated as follows:
M(x,y)=G(e(x,y))+min{M(x-1,y-1),M(x,y-1),M(x+1,y-1)} (4)
step (3.3), selecting the conversion frequency fcCorresponding to the last column of coordinate points (x)1M) as the starting point of the target ridge line, defining a matrix:
Figure GDA0003654260160000041
the previous column target ridge line position is (x)1+P(x1M), m-1), in reverse orderObtaining target ridge lines on the whole time interval from the rest positions of the tracked target ridge lines;
and (3.4) extracting the frequency corresponding to the target ridge line time point, namely the instantaneous frequency f (i delta t) of the variable-speed vibration signal, wherein i represents a time sequence.
Preferably, the step (4) extracts the amplitude and phase corresponding to the first order spectrum, and the specific steps are as follows:
step (4.1), the obtained instantaneous frequency f (i delta t) is integrated to obtain a phase-time mapping relation
Figure GDA0003654260160000042
Where f (Δ t) represents the instantaneous frequency at the first time interval Δ t, as shown in equation (6):
Figure GDA0003654260160000043
step (4.2), a time sequence t (k delta alpha) corresponding to the equal angle is obtained through a polynomial interpolation method, and the formula (7) shows that:
Figure GDA0003654260160000044
wherein, delta alpha is equal angle sampling interval, k is sample sequence, and spline (·) represents interpolation function;
then, according to the obtained equiangular time sequence, a signal x (k Δ α) after equiangular resampling can be obtained through further interpolation, wherein x (i Δ t) represents a vibration signal obtained by equiangular sampling, and the formula (8) is shown as follows:
x(k△α)=spline({i△t,x(i△t)},t(k△α)) (8)
step (4.3), obtaining an order spectrum by carrying out discrete Fourier transform on the signal subjected to equal-angle resampling; when the order spectrum conversion is carried out, certain errors exist in the frequency, amplitude and phase in the order spectrum analysis inevitably due to spectral leakage and barrier effect generated by truncation error and discretization processing, so that the first order spectrum is corrected by adopting a discrete spectrum correction method, and the amplitude and the phase corresponding to the first order spectrum are extracted.
As can be seen from the above description of the present invention, compared with the prior art, the present invention has the following advantages:
the invention provides a method for extracting dynamic balance fault characteristics of a rotor under a keyless phase change rotating speed, which can effectively avoid the problem of 'spectrum ambiguity' caused by the adoption of a traditional stationary signal-based spectrum analysis method.
The present invention is described in further detail below with reference to the accompanying drawings and embodiments, but the method for extracting the dynamic balance fault characteristics of the rotor at the keyless phase-change rotational speed is not limited to the embodiments.
Drawings
FIG. 1 is a flow chart of a method according to an embodiment of the present invention;
FIG. 2 is a time domain waveform diagram of a simulation signal according to an embodiment of the present invention;
FIG. 3 is a STFT time-frequency spectrum of a simulation signal according to an embodiment of the present invention;
FIG. 4 is a comparison diagram of the instantaneous frequency conversion and the actual instantaneous frequency conversion extracted by directly taking the extreme value of the time frequency spectrum of the simulation signal STFT according to the embodiment of the present invention;
FIG. 5 is a diagram of a target ridge line with the minimum energy gradient extracted by the method according to the embodiment of the present invention;
FIG. 6 is a comparison graph of the instantaneous frequency conversion extracted from the simulation signal by the method according to the embodiment of the present invention;
FIG. 7 is a diagram illustrating fundamental frequency amplitude extraction from a simulation signal by a conventional Fourier transform method according to an embodiment of the present invention;
fig. 8 is a diagram of extracting a fundamental phase from a simulation signal by using a conventional fourier transform method according to an embodiment of the present invention;
FIG. 9 is a graph of the first order amplitude extracted from the simulation signal by the present method according to the embodiment of the present invention;
FIG. 10 is a first-order phase diagram extracted from a simulation signal according to an embodiment of the present invention.
Detailed Description
The invention is further described below by means of specific embodiments.
The invention provides a method for extracting dynamic balance fault characteristics of a rotor under a keyless phase-change rotating speed. Based on the method, the dynamic balance fault characteristics of the variable-speed rotor are accurately extracted under the condition of no key phase signal, the hardware cost is saved, and the dynamic balance efficiency is improved.
Firstly, performing band-pass filtering and noise reduction treatment on an acquired vibration signal; performing time-frequency analysis on the vibration signals subjected to the noise reduction treatment through short-time Fourier transform to obtain time-frequency energy distribution; calculating by using a Sobel operator to obtain an energy gradient, and searching for a target ridge line with the minimum energy gradient by using a dynamic programming algorithm to adjust discontinuity of ridge line change caused by noise and adjacent frequency interference and realize accurate estimation of instantaneous frequency conversion of the variable-speed vibration signal; and performing order tracking calculation and discrete spectrum correction on the vibration signal subjected to noise reduction according to the extracted instantaneous frequency conversion, and extracting the amplitude and phase corresponding to a first order spectrum.
Referring to fig. 1, the method for extracting the dynamic balance fault characteristics of the rotor at the keyless phase-change rotating speed provided by the invention comprises the following specific implementation steps:
step (1), converting the frequency fcCentering, after setting a bandwidth, performing noise reduction processing on an original vibration signal by adopting a band-pass filter with a zero phase offset characteristic;
specifically, the harmonic wavelet filter with the box-shaped characteristic is adopted, and the phase difference between the real part and the imaginary part of the harmonic wavelet is 90 degrees, so that the phase of a signal is kept unchanged before and after the harmonic wavelet transform, and the accuracy of phase extraction in the dynamic balance fault characteristic of the rotor is ensured.
Selecting proper window length and overlapping number, performing time-frequency analysis on the vibration signal x (k) subjected to noise reduction treatment by adopting short-time Fourier transform to obtain a time-frequency spectrogram, and performing logarithm treatment to obtain a time-frequency energy distribution map, wherein N is the signal length;
the specific steps of calculating the time-frequency energy distribution map are as follows:
and (2.1) performing discrete short-time Fourier transform on the discrete vibration signal subjected to the noise reduction treatment, wherein the discrete short-time Fourier transform comprises the following steps:
Figure GDA0003654260160000061
where S (i Δ t, f) represents the transformed spectrum, f represents the frequency, Δ t is 1/fsRepresenting the sampling time interval, f sRepresenting the sampling frequency, i represents the time series, k is 1,2, …, N represents the sample series, N represents the signal length, x represents the complex conjugate, x (k) represents the discrete form of the noise-reduced vibration signal, g (k) represents the discrete form of the window function g (t);
and (2.2) converting the time frequency spectrum after the discrete short-time Fourier transform into energy distribution consisting of time discrete points and frequency discrete points, wherein the conversion formula is as follows:
P(i△t,f)=10*log10(S(i△t,f)) (2)
step (3), calculating the obtained time-frequency energy distribution graph through a Sobel operator to obtain an energy gradient, and carrying out normalization processing, wherein the normalized energy gradient graph comprises a plurality of ridges, namely corresponding spectral lines of different orders; converting the extraction problem of the ridge line into an optimization problem, and selecting the frequency conversion fcThe corresponding points are used as the starting points of the target ridge lines, the target ridge lines with the minimum energy gradient are searched by adopting a dynamic programming algorithm, and the frequency corresponding to each time point of the target ridge lines is the instantaneous frequency conversion of the variable-speed vibration signals in the whole time interval;
the method specifically comprises the following steps of accurately estimating the instantaneous frequency f (i delta t) of the variable-speed vibration signal:
step (3.1), regarding the time-frequency energy distribution P (i Δ t, f) as an energy function e (x, y) corresponding to the coordinate point (x, y), calculating the energy gradient of the time-frequency spectrum by using a Sobel operator, and performing normalization processing, wherein the energy gradient G (e (x, y)) is calculated as follows:
G(e(x,y))=abs(filter2(hx,P(i△t,f)))+abs(filter2(hy,P(i△t,f))) (3)
In the formula, hx and hy are Sobel operators in the horizontal direction and the vertical direction respectively,
Figure GDA0003654260160000071
Figure GDA0003654260160000072
filter2(·) represents a two-dimensional digital filter function.
And (3.2) because the local maximum value and the real maximum value in the time frequency spectrum caused by the noise shift, the frequency jump is large and the energy gradient is increased, the ridge line formed by the local maximum values in the time frequency spectrum is discontinuously changed, so that the target ridge line is extracted by searching the minimum value point of the energy gradient to adjust the ridge line discontinuity caused by the noise and the adjacent frequency interference, and the optimal target ridge line is extracted by adopting an optimal decision dynamic programming algorithm. Defining a matrix M (x, y), initializing M (x,1) to G (e (x,1)), the remaining M (x, y) values from the beginning of the second column to the last column are calculated as follows:
M(x,y)=G(e(x,y))+min{M(x-1,y-1),M(x,y-1),M(x+1,y-1)} (4)
step (3.3), selecting the conversion frequency fcCorresponding to the last column of coordinate points (x)1M) as the starting point of the target ridge, defining a matrix:
Figure GDA0003654260160000073
the previous column target ridge position is (x)1+P(x1M), m-1), and sequentially reversely tracking the rest positions of the target ridge line to obtain the target ridge line in the whole time interval.
And (3.4) extracting the frequency corresponding to the target ridge line time point, namely the instantaneous frequency f (i delta t) of the variable-speed vibration signal, wherein i represents a time sequence.
And (4) performing order tracking calculation on the original vibration signal according to the extracted instantaneous frequency conversion, and extracting the amplitude and the phase corresponding to the first-order spectrum.
The specific steps of obtaining the amplitude and the phase corresponding to the first order spectrum through order tracking calculation and extraction are as follows:
step (4.1), the obtained instantaneous frequency f (i delta t) is integrated to obtain a phase-time mapping relation
Figure GDA0003654260160000074
As shown in formula (6)
Figure GDA0003654260160000075
And (4.2) obtaining a time sequence t (k delta alpha) corresponding to the equal angle through a polynomial interpolation method, wherein delta alpha is the equal angle sampling interval, and is shown in a formula (7).
Figure GDA0003654260160000081
And then according to the obtained equiangular time sequence, obtaining a signal x (k delta alpha) after equiangular resampling by further interpolation, wherein the formula is shown in a formula (8).
x(k△α)=spline({i△t,x(i△t)},t(k△α)) (8)
And (4.3) carrying out discrete Fourier transform on the signals subjected to equal-angle resampling to obtain an order spectrum. When the order spectrum conversion is carried out, certain errors exist in the frequency, the amplitude and the phase in the order spectrum analysis due to frequency spectrum leakage and barrier effect which are inevitably generated by truncation errors and discretization processing, and therefore the first-order spectrum is corrected based on the discrete spectrum correction method, and accurate extraction of the amplitude and the phase corresponding to the first-order spectrum is achieved.
The method of the present invention is described below using a set of rotor imbalance simulation data at varying rotational speeds.
The sampling frequency of the simulation data is 2048Hz, the sampling time is 3s, the signals comprise fundamental frequency components of amplitude modulation and phase modulation, adjacent frequency components of frequency multiplication of 2 and 3, and noise components, the fundamental frequency changes along with time, and the simulation signals are shown in a formula (9).
Figure GDA0003654260160000082
Wherein the fundamental frequency is f 015+0.15sin (4 π t), the signal-to-noise ratio η (t) is 10dB of white Gaussian noise. See fig. 2, which is a time domain waveform of the simulated signal.
And filtering and denoising the original simulation signal by adopting harmonic wavelets, so as to improve the signal-to-noise ratio. Then, a hamming window with the length of 513 is selected, the number of overlapping windows is 512, STFT transformation is performed on the denoised signal to obtain a time-frequency spectrogram, see fig. 3, the instantaneous frequency conversion extracted by the traditional extremum extraction method is susceptible to interference of adjacent frequency and noise, and the relative error between the extracted instantaneous frequency conversion and the actual instantaneous frequency conversion is 18.7%, see fig. 4. Calculating the time-frequency spectrogram to obtain a time-frequency energy distribution graph, calculating the obtained time-frequency energy distribution graph through a Sobel operator to obtain an energy gradient, and carrying out normalization processing; a dynamic programming algorithm is used to search for the target ridge with the smallest energy gradient, see fig. 5. The frequency corresponding to each time point of the target ridge line is the instantaneous frequency conversion estimated value in the whole time interval of the variable-speed vibration signal, the relative error between the instantaneous frequency conversion extracted by the method and the actual instantaneous frequency conversion is 0.95%, and the method is shown in figure 6. And integrating the extracted instantaneous frequency to obtain an instantaneous phase, and converting the instantaneous phase into an order spectrum by using an order tracking method. And then, correcting the first order spectrum by adopting a discrete spectrum correction method to obtain the accurate amplitude and phase corresponding to the first order spectrum. Fig. 9 and 10 show the order amplitude spectrum and the order phase spectrum, respectively, and the spectra are clear, and the corresponding amplitude and phase of the first order spectrum are 0.998 and 49.43 ° respectively, which are close to the theoretical values of 1 and 50 °. For comparison, the spectral analysis using the conventional fourier transform method based on stationary signal analysis showed severe "frequency ambiguity", see fig. 7 and 8, where the fundamental amplitude and phase were 0.507 and-174.45 °, respectively, which are very different from the theoretical values of 1 and 50 °. Therefore, the method provided by the invention can realize accurate identification of the rotor unbalance fault characteristics.
The above description is only an embodiment of the present invention, but the design concept of the present invention is not limited thereto, and any insubstantial modifications made by using the design concept should fall within the scope of infringing the present invention.

Claims (4)

1. A rotor dynamic balance fault feature extraction method under a keyless phase-change rotating speed is characterized by comprising the following steps:
step (1), converting the frequency fcCentering, after setting a bandwidth, performing noise reduction processing on an original vibration signal by adopting a band-pass filter with a zero phase offset characteristic;
selecting a proper window length and a proper overlap number, carrying out time-frequency analysis on the vibration signals x (k) after noise reduction treatment by adopting short-time Fourier transform to obtain a time-frequency spectrogram, and carrying out logarithm treatment to obtain a time-frequency energy distribution map, wherein N is the signal length;
step (3), calculating the obtained time-frequency energy distribution graph through a Sobel operator to obtain an energy gradient, and carrying out normalization processing, wherein the energy gradient graph after the normalization processing comprises a plurality of ridges, namely corresponding spectral lines of different orders; converting the extraction problem of the ridge line into an optimization problem, and selecting the frequency conversion fcThe corresponding points are used as the starting points of the target ridge lines, the target ridge lines with the minimum energy gradient are searched by adopting a dynamic programming algorithm, and the frequency corresponding to each time point of the target ridge lines is the instantaneous frequency conversion of the variable-speed vibration signals in the whole time interval;
And (4) performing order tracking calculation on the vibration signal subjected to noise reduction treatment in the step (1) according to the extracted instantaneous frequency conversion to obtain an order spectrum, performing discrete spectrum correction on the first order spectrum in the order spectrum, and extracting to obtain an amplitude and a phase corresponding to the first order spectrum, namely the rotor dynamic balance fault characteristic at variable rotation speed.
2. The method for extracting the dynamic balance fault characteristics of the rotor at the keyless phase-change rotating speed according to claim 1, wherein the step (2) specifically comprises the following steps:
and (2.1) performing discrete short-time Fourier transform on the discrete vibration signal subjected to the noise reduction treatment, wherein the discrete short-time Fourier transform comprises the following steps:
Figure FDA0003654260150000011
where S (i Δ t, f) represents the transformed time-frequency spectrum, f represents the frequency, and Δ t is 1/fsRepresenting the sampling time interval, fsDenotes the sampling frequency, i denotes the time series, k 1, 2., N denotes the sample series, N denotes the signal length, x denotes the complex conjugate, x (k) denotes the discrete form of the noise-reduced vibration signal, g (k) denotes the discrete form of the window function g (t);
and (2.2) converting the time-frequency spectrum after the discrete short-time Fourier transform into energy distribution P (i delta t, f) consisting of time discrete points and frequency discrete points, wherein the conversion formula is as follows:
P(iΔt,f)=10*log10(S(iΔt,f)) (2)。
3. the method for extracting the dynamic balance fault characteristics of the rotor at the keyless phase-change rotating speed according to claim 2, wherein the step (3) specifically comprises:
Step (3.1), regarding the time-frequency energy distribution P (i Δ t, f) as an energy function e (x, y) corresponding to the coordinate point (x, y), calculating the energy gradient of the time-frequency spectrum by using a Sobel operator, and performing normalization processing, wherein the energy gradient G (e (x, y)) is calculated as follows:
G(e(x,y))=abs(filter2(hx,P(iΔt,f)))+abs(filter2(hy,P(iΔt,f))) (3)
in the formula, hx and hy are Sobel operators in the horizontal direction and the vertical direction respectively, and the filter2(·) represents a two-dimensional digital filter function;
step (3.2), defining the matrix M (x, y), initializing M (x, 1) to G (e (x, 1)), then the remaining M (x, y) values from the second column to the last column are calculated as follows:
M(x,y)=G(e(x,y))+min{M(x-1,y-1),M(x,y-1),M(x+1,y-1)} (4)
step (3.3), selecting the conversion frequency fcCorresponding to the last column of coordinate points (x)1M) as the starting point of the target ridge line, defining a matrix:
Figure FDA0003654260150000021
the previous column target ridge line position is (x)1+P(x1M), m-1), sequentially and reversely tracking the rest positions of the target ridge line to obtain the target ridge line in the whole time interval;
and (3.4) extracting the frequency corresponding to the target ridge line time point, namely the instantaneous frequency f (i delta t) of the variable-speed vibration signal, wherein i represents a time sequence.
4. The method for extracting the dynamic balance fault characteristics of the rotor at the keyless phase-change rotating speed according to claim 3, wherein the step (4) specifically comprises the following steps:
step (4.1), the obtained instantaneous frequency f (i delta t) is integrated to obtain a phase-time mapping relation
Figure FDA0003654260150000023
Where f (Δ t) represents the instantaneous frequency at the first time interval Δ t, as shown in equation (6):
Figure FDA0003654260150000022
and (4.2) obtaining a time sequence t (k delta alpha) corresponding to the equal angle by a polynomial interpolation method, wherein the formula (7) is as follows:
Figure FDA0003654260150000031
wherein, delta alpha is an equal angle sampling interval, k is a sample sequence, and spline (·) represents an interpolation function;
then, according to the obtained equiangular time sequence, a signal x (k Δ α) after equiangular resampling can be obtained through further interpolation, wherein x (i Δ t) represents a vibration signal obtained by equiangular sampling, and the formula (8) is shown as follows:
x(kΔα)=spline({iΔt,x(iΔt)},t(kΔα)) (8)
step (4.3), obtaining an order spectrum by performing discrete Fourier transform on the signals subjected to equal-angle resampling; when the order spectrum conversion is carried out, certain errors exist in the frequency, amplitude and phase in the order spectrum analysis inevitably due to spectral leakage and barrier effect generated by truncation error and discretization processing, so that the first order spectrum is corrected by adopting a discrete spectrum correction method, and the amplitude and the phase corresponding to the first order spectrum are extracted.
CN202111188122.4A 2021-10-12 2021-10-12 Rotor dynamic balance fault feature extraction method under keyless phase-change rotating speed Active CN113984282B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111188122.4A CN113984282B (en) 2021-10-12 2021-10-12 Rotor dynamic balance fault feature extraction method under keyless phase-change rotating speed

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111188122.4A CN113984282B (en) 2021-10-12 2021-10-12 Rotor dynamic balance fault feature extraction method under keyless phase-change rotating speed

Publications (2)

Publication Number Publication Date
CN113984282A CN113984282A (en) 2022-01-28
CN113984282B true CN113984282B (en) 2022-07-19

Family

ID=79738268

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111188122.4A Active CN113984282B (en) 2021-10-12 2021-10-12 Rotor dynamic balance fault feature extraction method under keyless phase-change rotating speed

Country Status (1)

Country Link
CN (1) CN113984282B (en)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114674410A (en) * 2022-03-24 2022-06-28 安徽大学 Component number time-varying underwater acoustic signal instantaneous frequency estimation method

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106370403A (en) * 2016-08-22 2017-02-01 南京信息工程大学 Instant frequency estimation method based on edge detection
CN107368456A (en) * 2017-07-24 2017-11-21 潍坊学院 The instantaneous Frequency Estimation method examined based on Sobel operators and t
CN108871742A (en) * 2018-05-03 2018-11-23 西安交通大学 A kind of improved no key phase fault feature order extracting method
CN110887663A (en) * 2019-10-30 2020-03-17 中国石油化工股份有限公司 Bearing fault diagnosis method combining variable working condition calculation order tracking and spectral kurtosis
CN111307452A (en) * 2020-03-05 2020-06-19 江苏天沃重工科技有限公司 Intelligent fault diagnosis method for rotating machinery at time-varying rotating speed
CN112665851A (en) * 2020-10-28 2021-04-16 西安交通大学 Key-free phase change rotating speed gearbox fault diagnosis method
CN113125179A (en) * 2021-03-11 2021-07-16 同济大学 Keyless phase order tracking method for rotating speed fluctuation of rotary machine

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7321809B2 (en) * 2003-12-30 2008-01-22 The Boeing Company Methods and systems for analyzing engine unbalance conditions

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106370403A (en) * 2016-08-22 2017-02-01 南京信息工程大学 Instant frequency estimation method based on edge detection
CN107368456A (en) * 2017-07-24 2017-11-21 潍坊学院 The instantaneous Frequency Estimation method examined based on Sobel operators and t
CN108871742A (en) * 2018-05-03 2018-11-23 西安交通大学 A kind of improved no key phase fault feature order extracting method
CN110887663A (en) * 2019-10-30 2020-03-17 中国石油化工股份有限公司 Bearing fault diagnosis method combining variable working condition calculation order tracking and spectral kurtosis
CN111307452A (en) * 2020-03-05 2020-06-19 江苏天沃重工科技有限公司 Intelligent fault diagnosis method for rotating machinery at time-varying rotating speed
CN112665851A (en) * 2020-10-28 2021-04-16 西安交通大学 Key-free phase change rotating speed gearbox fault diagnosis method
CN113125179A (en) * 2021-03-11 2021-07-16 同济大学 Keyless phase order tracking method for rotating speed fluctuation of rotary machine

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
变转速下机械设备动态信号分析方法的回顾与展望;林京 等;《中国科学》;20151231;第45卷(第7期);669-686 *
齿轮箱关键部件非平稳振动信号分析及诊断方法研究;江星星;《中国博士学位论文全文数据库 工程科技II辑》;20171115(第11期);144-149 *

Also Published As

Publication number Publication date
CN113984282A (en) 2022-01-28

Similar Documents

Publication Publication Date Title
CN111665051A (en) Bearing fault diagnosis method under strong noise variable-speed condition based on energy weight method
CN110617964A (en) Synchronous compression transformation order ratio analysis method for fault diagnosis of rolling bearing
CN108871742B (en) Improved key-phase-free fault feature order extraction method
Wang et al. Rolling element bearing fault diagnosis via fault characteristic order (FCO) analysis
CN110907162B (en) Rotating machinery fault feature extraction method without tachometer under variable rotating speed
Urbanek et al. Comparison of amplitude-based and phase-based method for speed tracking in application to wind turbines
CN110763462B (en) Time-varying vibration signal fault diagnosis method based on synchronous compression operator
Kang et al. Phase difference correction method for phase and frequency in spectral analysis
CN112665851B (en) Key-free phase change rotating speed gearbox fault diagnosis method
CN113125179A (en) Keyless phase order tracking method for rotating speed fluctuation of rotary machine
CN110987167A (en) Fault detection method, device, equipment and storage medium for rotary mechanical equipment
CN113984282B (en) Rotor dynamic balance fault feature extraction method under keyless phase-change rotating speed
Wu et al. A modified tacho-less order tracking method for the surveillance and diagnosis of machine under sharp speed variation
CN112362343A (en) Distributed fault feature extraction method for gearbox under variable rotating speed based on frequency modulation dictionary
Lin et al. A review and strategy for the diagnosis of speed-varying machinery
CN113899444A (en) Resonance frequency measurement method of vibrating wire sensor based on Hanning double windows
CN114509159A (en) Order ratio tracking analysis method, system and computer readable storage medium
CN111323233B (en) Local mean decomposition method for low-speed rotating machine fault diagnosis
CN117419923A (en) Pipelined hardware phase resolving method and resolving system suitable for engine
CN112781723B (en) Harmonic component detection method based on frequency spectrum variance
CN104677486A (en) Aero-engine vibration signal phase measurement method based on revolving speed pulse reconstruction
CN112665712B (en) Wide-area order tracking method and system for monitoring train running gear
CN116256158A (en) Rotary machine instantaneous phase self-adaptive extraction method based on depth signal separation
CN110580471A (en) Mechanical equipment fault diagnosis method based on encoder signal transient characteristics
Zhang et al. Research on the identification of asynchronous vibration parameters of rotating blades based on blade tip timing vibration measurement theory

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
GR01 Patent grant
GR01 Patent grant