1. A rotary machine fault diagnosis method based on iterative extended eigenmode decomposition is characterized by comprising the following steps:

step 1: acquiring a vibration acceleration signal of a rotating machine in the operation of a variable rotating speed working condition;

step 2: carrying out discrete short-time Fourier transform on the vibration acceleration signal to obtain a time-frequency expression;

and step 3: acquiring an initial estimation value of a key time-frequency ridge line with the most significant energy in the time-frequency expression by a bidirectional rapid single-line extraction method;

and 4, step 4: separating and reconstructing a key mode, an instantaneous frequency and an instantaneous amplitude thereof according to the vibration acceleration signal and a key time-frequency ridge line initial estimation value through iterative expansion eigenmode decomposition;

and 5: performing equal-angle resampling on the key modes to obtain steady-state angular domain signals;

step 6: and converting the angle domain into an order domain through discrete Fourier transform, and analyzing each characteristic order to obtain a fault diagnosis result of the rotary machine.

2. The method for fault diagnosis of rotating machinery based on iterative extended eigenmode decomposition according to claim 1, wherein the step 2 comprises the following steps:

step 201: setting the Window Length N_windowAnd the number of Fourier points L_fourierConstructing a corresponding window function wf (t);

step 202: and performing discrete short-time Fourier transform on the vibration acceleration signal x (t) through the window function wf (t) to obtain a time-frequency expression TFR (t, f).

3. The method for diagnosing faults of rotating machinery based on iterative extended eigenmode decomposition according to claim 1, wherein in the step 3, the bidirectional fast single-line extraction method specifically comprises: and positioning an energy maximum value point in the time-frequency expression, searching an energy maximum value point of a frequency dimension in a delta f range in a two-way mode along the time direction, and constructing an initial estimation value of the key time-frequency ridge line, wherein the delta f is preset.

4. The method according to claim 1, wherein the step 3 specifically comprises the following steps:

step 301: initializing a maximum frequency difference value delta f;

step 302: locating the energy maximum point (t) in the time-frequency expression_em,f_em) The positioning expression is as follows:

in the formula, TFR (t, f) is a time-frequency expression,is the starting point of the initial estimated value of the critical time-frequency ridge line;

the energy maximum value point of the frequency dimension in the range of delta f is searched in the time direction in a two-way mode, and the searching expression is as follows:

f_R＝f_em,f_L＝f_em

in the formula (I), the compound is shown in the specification,the initial estimation value of the key time-frequency ridge line is obtained, and N is the number of sampling points;

step 303: obtaining the initial estimation value of the critical time-frequency ridge line

5. The method according to claim 1, wherein the step 4 specifically comprises the following steps:

step 401: initializing parameters: regularization parameter lambda, punishment parameters beta and mu and Fourier expansion order L, and obtaining vibration acceleration signal x (t) and initial estimation value of key time-frequency ridge lineThe calculation expression of the Fourier expansion order L is as follows:

wherein BW is the bandwidth, N is the number of sampling points, f_sIs the sampling frequency;

step 402: performing Hilbert transform on the vibration acceleration signal x (t) to construct a key mode m_key(t) analytical formula

Where k is 0.. end, end is the number of iterations,is the k iterations of the key instantaneous frequency, η (t) is the interference noise, t is the discrete sampling time,to resolve the amplitude:

wherein the content of the first and second substances,in order to be the initial phase position,is the true value of the critical instantaneous frequency, a_key(t) is the instantaneous amplitude of the key mode;

then, the analytic amplitude is establishedRedundant fourier expansion of (a):

wherein the content of the first and second substances, l in (1) is the order of the Fourier series;

step 403: will be described inUnfolding model bring-inConstructing a matrix expression:

K^k＝E^kt is a kernel matrix of N × (2L +1),

in the formula, d is a matrix index position;

the optimization problem is then solved:

in the formula (I), the compound is shown in the specification,is a signal analytic formula;

obtaining a Fourier parameter matrixAnd is reconstituted to obtain

In the formula, I is an identity matrix, and H is a matrix transpose;

step 404: unfolding the reconstruction based on the Euler formulaGetReal part and imaginary part, and estimating error by arc tangent demodulation methodTo pairSmoothing the filter and updating the instantaneous frequency

Step 405: updating the updated values in step 404Substitution step 403 updates the kernel matrix K^k+1Parameter matrixAnd resolving the magnitude matrixRepeatedly calculating the preset end times, stopping, and solving the reconstructed key mode m_key(t) and its amplitude a_key(t)。

6. The method according to claim 5, wherein the method for fault diagnosis of rotating machinery based on iterative extended eigenmode decomposition is characterized in thatThe euler equation of (a) is:

the error isThe estimated expression of (c) is:

in the formula, R^kIs the real part of the complex envelope, I^kIs the imaginary part of the complex envelope.

7. The method according to claim 6, wherein the instantaneous frequency is the instantaneous frequencyThe update expression of (1) is:

wherein the content of the first and second substances,to improve the second order difference operator.

8. The method according to claim 5, wherein the critical mode m is a critical mode m_keyThe computational expression of (t) is:

m_key＝u·A+v·B

the amplitude a_keyThe computational expression of (t) is:

9. the method according to claim 1, wherein the step 5 specifically comprises the following steps:

step 501: for the key mode m_key(t) conducting HillPerforming Bert transform, performing inverse tangent on the ratio of imaginary part to real part, and performing deconvolution processing to obtain instantaneous phase curve

Step 502: establishing the instantaneous phase curve by least squares fittingTime-radian mapping relation psi_key(t)；

Step 503: determining a sampling order O_mAnd equal radian interval delta theta are summed, and the mapping relation is reversely solvedObtaining the equal radian time sequence t_Δθ(k) And radian sequence rad_Δθ(k)；

Step 504: for the time sequence t_Δθ(k) And the vibration acceleration signal x (t) is subjected to cubic spline interpolation to obtain an equal radian amplitude sequence a_Δθ(k)。

10. The method according to claim 1, wherein the step 6 specifically includes the following steps:

step 601: adding a Hanning window w to the steady-state angular domain signal_hanning(k) According to said sampling order O_mCalculating order sequence O (k), calculating amplitude sequence A by discrete Fourier transform_Δθ(k) Transition to the order domain;

step 602: in the order domain, if the first order is the significant energy and the lowest order, directly analyzing each feature order proportional to the first order, including: first order and its multiple order, fault signature order and its multiple order;

in the order domain, if the first order is not energy significant and the lowest order, the energy significant and the lowest order smaller than the first order is located and is inverted to obtain the magnification factor m, the order sequence o (k) is magnified in equal proportion by m, and then step 602 is repeatedly executed.

Background

The dynamic signal of mechanical equipment under variable rotating speed is very complex, and the signal can have non-stationary characteristics such as frequency modulation, amplitude modulation, phase modulation and the like, so that serious frequency spectrum aliasing is caused, and fault diagnosis becomes very difficult. However, in actual engineering, the variable-speed working condition is almost ubiquitous, and especially under the large-speed fluctuation working conditions such as starting and stopping of the vehicle, the variable-speed working condition is difficult to analyze and identify.

Order tracking (Order tracking) is the most direct and effective method for solving the problem of rotation speed fluctuation, and can be divided into three types, namely Hardware Order Tracking (HOT), Calculation Order Tracking (COT) and keyless phase Order tracking (TLOT). The hardware order tracking mainly has the problems of hardware equipment support, difficulty in arrangement and implementation, high cost and the like; the problems that a key phase trigger needs to be installed, the calculation accuracy depends on the quality of a rotating speed signal and the like mainly exist in calculation order tracking; the keyless phase order tracking mainly has the problems of large algorithm calculation consumption, low calculation precision, applicability to weak rotation speed fluctuation, much human intervention and the like.

The invention patent application related to the order tracking technology in China comprises the following steps:

the invention discloses an invention patent named as a gear fault keyless phase angle domain average calculation order analysis method in Chinese patent with the authorization number of CN103499443B, wherein the authorization number is 2016, 1, 20 and 10. The method estimates the instantaneous frequency conversion by utilizing smooth pseudo Wigner-Ville time frequency analysis and a Viterbi optimal path search algorithm, and then performs equal-angle resampling based on an estimated key phase signal, so that keyless phase order tracking is effectively performed on the gear box under the working condition of variable rotating speed, dependence on hardware is eliminated, but the method has the defects that the method is only suitable for weak rotating speed fluctuation, is easy to be interfered by noise and cannot obtain order time domain information.

The patent of the invention is named as 'Vold-Kalman filtering bandwidth optimization method based on order spectrum' in Chinese patent with the authorization date of 2018, 5, month and 1 and the authorization number of CN 105512369B. The order spectrums of confidence intervals are respectively obtained by utilizing calculation order tracking and Vold-Kalman order tracking, the optimal filtering bandwidth in a bandwidth library is selected by calculating the standard deviation of a residual signal, the problem of Vold-Kalman order tracking bandwidth selection is effectively solved, the accuracy of order analysis is improved, but the method has the defects that the method is only suitable for weak rotation speed fluctuation, the calculation consumption is large, and rotation speed information needs to be predicted.

The invention discloses an invention patent named as a calculation order tracking method for self-adaptive noise reduction and order aliasing avoidance in a Chinese patent with the publication number of CN110084208A, wherein the publication number is 8/2/2019. The method utilizes the rotation speed information to define the margin frequency, decomposes and reconstructs signals based on the arrangement entropy PE and the VMD optimized by the differential evolution algorithm, and then performs calculation order tracking, thereby effectively and adaptively reducing noise interference and overcoming the aliasing problem of order tracking under variable working conditions, but the method has the defects of large calculation consumption and dependence of calculation precision on the quality of the rotation speed signals;

the invention discloses a rolling bearing fault diagnosis method based on calculation order tracking and spectral kurtosis in Chinese patent with publication number CN111307460A, namely 6/19/2020. The method utilizes spectral kurtosis and envelope demodulation to preprocess vibration signals, then sets a rotating speed pulse threshold value to carry out calculation order tracking, and effectively extracts bearing fault characteristic orders under variable working conditions, but has the defects of needing a key phase trigger, needing to set the pulse threshold value and having more manual intervention.

Disclosure of Invention

The invention aims to overcome the defects of the prior art that the rotating speed needs to be predicted, the calculation precision is limited, the algorithm calculation consumption is high, the method is only suitable for weak rotating speed fluctuation, multiple human interventions and the like, and provides a rotating machine fault diagnosis method based on iterative expansion eigenmode decomposition, which realizes the fault diagnosis of mechanical equipment under the working condition of variable rotating speed by efficiently analyzing and extracting key eigencomponents in vibration signals.

The purpose of the invention can be realized by the following technical scheme:

a rotary machine fault diagnosis method based on iterative extended eigenmode decomposition comprises the following steps:

step 1: acquiring a vibration acceleration signal of a rotating machine in the operation of a variable rotating speed working condition;

step 2: carrying out discrete short-time Fourier transform on the vibration acceleration signal to obtain a time-frequency expression;

and step 3: acquiring an initial estimation value of a key time-frequency ridge line with the most significant energy in the time-frequency expression by a bidirectional rapid single-line extraction method;

and 4, step 4: separating and reconstructing a key mode, an instantaneous frequency and an instantaneous amplitude thereof according to the vibration acceleration signal and a key time-frequency ridge line initial estimation value through iterative expansion eigenmode decomposition;

and 5: performing equal-angle resampling on the key modes to obtain steady-state angular domain signals;

step 6: and converting the angle domain into an order domain through discrete Fourier transform, and analyzing each characteristic order to obtain a fault diagnosis result of the rotary machine.

Further, the step 2 comprises the following steps:

step 201: setting the Window Length N_windowAnd the number of Fourier points L_fourierConstructing a corresponding window function wf (t);

step 202: and performing discrete short-time Fourier transform on the vibration acceleration signal x (t) through the window function wf (t) to obtain a time-frequency expression TFR (t, f).

Further, in step 3, the bidirectional fast single-line extraction method specifically comprises: and positioning an energy maximum value point in the time-frequency expression, searching an energy maximum value point of a frequency dimension in a delta f range in a two-way mode along the time direction, and constructing an initial estimation value of the key time-frequency ridge line, wherein the delta f is preset.

Further, the step 3 specifically includes the following steps:

step 301: initializing a maximum frequency difference value delta f;

step 302: locating the energy maximum point (t) in the time-frequency expression_em,f_em) The positioning expression is as follows:

in the formula, TFR (t, f) is a time-frequency expression,is the starting point of the initial estimated value of the critical time-frequency ridge line;

and (3) searching energy maximum value points of frequency dimensions in a delta f range in a two-way mode along the time direction: the search expression is:

f_R＝f_em,f_L＝f_em

in the formula (I), the compound is shown in the specification,the initial estimation value of the key time-frequency ridge line is obtained, and N is the number of sampling points;

step 303: obtaining the initial estimation value of the critical time-frequency ridge line

Further, the step 4 specifically includes the following steps:

step 401: initializing parameters: regularization parameter lambda, punishment parameters beta and mu and Fourier expansion order L, and obtaining vibration acceleration signal x (t) and initial estimation value of key time-frequency ridge lineThe calculation expression of the Fourier expansion order L is as follows:

wherein BW is the bandwidth, N is the number of sampling points, f_sIs the sampling frequency;

step 402: subjecting the vibration acceleration signal x (t) to Hilbert transformConstruction of the Key mode m_key(t) analytical formula

Where k is 0.. end, end is the number of iterations,is the k iterations of the key instantaneous frequency, η (t) is the interference noise, t is the discrete sampling time,to resolve the amplitude:

wherein the content of the first and second substances,in order to be the initial phase position,is the true value of the critical instantaneous frequency, a_key(t) is the instantaneous amplitude of the key mode;

then, the analytic amplitude is establishedRedundant fourier expansion of (a):

wherein the content of the first and second substances, l in (1) is the order of the Fourier series;

step 403: will be described inUnfolding model bring-inConstructing a matrix expression:

K^k＝E^kt is a kernel matrix of N × (2L +1),

in the formula, d is a matrix index position;

the optimization problem is then solved:

in the formula (I), the compound is shown in the specification,is a signal analytic formula;

obtaining a Fourier parameter matrixAnd is reconstituted to obtain

In the formula, I is an identity matrix, and H is a matrix transpose;

step 404: unfolding the reconstruction based on the Euler formulaGetReal part and imaginary part, and estimating error by arc tangent demodulation methodTo pairSmoothing the filter and updating the instantaneous frequency

Step 405: updating the updated values in step 404Substitution step 403 updates the kernel matrix K^k+1Parameter matrixAnd resolving the magnitude matrixRepeatedly calculating the preset end times, stopping, and solving the reconstructed key mode m_key(t) and its amplitude a_key(t)。

Further, theThe euler equation of (a) is:

the error isThe estimated expression of (c) is:

in the formula, R^kIs the real part of the complex envelope, I^kIs the imaginary part of the complex envelope.

Further, the instantaneous frequencyThe update expression of (1) is:

wherein the content of the first and second substances,to improve the second order difference operator.

Further, the key modality m_keyThe computational expression of (t) is:

m_key＝u·A+v·B

the amplitude a_keyThe computational expression of (t) is:

further, the step 5 specifically includes the following steps:

step 501: for the key mode m_key(t) performing Hilbert transform, arctangent of the ratio of the imaginary part to the real part, and deconvoluting to obtain an instantaneous phase curve

Step 502: establishing the instantaneous phase curve by least squares fittingTime-radian mapping relation psi_key(t)；

Step 503: determining a sampling order O_mAnd equal radian interval delta theta are summed, and the mapping relation is reversely solvedObtaining the equal radian time sequence t_Δθ(k) Andradian sequence rad_Δθ(k)；

Step 504: for the time sequence t_Δθ(k) And the vibration acceleration signal x (t) is subjected to cubic spline interpolation to obtain an equal radian amplitude sequence a_Δθ(k)。

Further, the step 6 specifically includes the following steps:

step 601: adding a Hanning window w to the steady-state angular domain signal_hanning(k) According to said sampling order O_mCalculating order sequence O (k), calculating amplitude sequence A by discrete Fourier transform_Δθ(k) Transition to the order domain;

step 602: in the order domain, if the first order is the significant energy and the lowest order, directly analyzing each feature order proportional to the first order, including: first order and its multiple order, fault signature order and its multiple order;

in the order domain, if the first order is not energy significant and the lowest order, the energy significant and the lowest order smaller than the first order is located and is inverted to obtain the magnification factor m, the order sequence o (k) is magnified in equal proportion by m, and then step 602 is repeatedly executed.

Compared with the prior art, the invention has the following advantages:

the invention provides a bidirectional fast single-line extraction method (BFSE) and an iterative extended eigenmode decomposition method (IEIMD), which are effectively combined with a keyless phase order tracking technology (TLOT) to realize fault diagnosis of a rotary machine, and carry out signal processing and fault diagnosis on mechanical equipment with variable rotating speed, particularly large rotating speed fluctuation working conditions by accurately extracting and analyzing key intrinsic components in vibration signals, and have the following remarkable advantages different from the common technology:

(1) hardware equipment such as a control circuit and a key phase trigger is not needed, operation and maintenance cost and application difficulty are obviously reduced, meanwhile, the method does not depend on rotating speed information, and calculation precision is not influenced by rotating speed signal quality;

(2) compared with the similar method, the calculation speed is improved by more than 2 times, the instantaneous frequency extraction precision is improved by more than 20%, the modal reconstruction precision is improved by more than 15%, the characteristic order error is less than 1%, the calculation precision is improved by more than 20%, more human intervention is not needed, and the embedded type development and the actual engineering application are facilitated;

(3) the invention is suitable for the condition of large rotation speed fluctuation (the instantaneous rotation frequency change rate is more than 30 percent), and can be used for complex working conditions such as starting and stopping, rapid speed increasing, speed reducing and the like;

(4) the method accurately reconstructs the time domain waveform, the instantaneous frequency and the instantaneous amplitude of the key order, effectively inhibits the noise interference in the original signal, has rich and complete information and is favorable for deeply analyzing the fault.

Drawings

FIG. 1 is a schematic flow chart of a fault diagnosis method for a rotary machine based on iterative extended eigenmode decomposition according to the present invention;

FIG. 2(a) is a time domain diagram of an emulated signal in an embodiment of the present invention;

FIG. 2(b) is a frequency domain diagram of a simulated signal in an embodiment of the present invention;

FIG. 3(a) is a frequency spectrum and a partially enlarged view of an exemplary simulated signal according to the present invention;

FIG. 3(b) is a diagram illustrating an estimation of a critical time-frequency ridge of a BFSE in the simulation embodiment of the present invention;

FIG. 3(c) is a diagram illustrating an extraction of critical time-frequency ridge lines of an IEIMD according to an embodiment of the present invention;

FIG. 3(d) is a graph illustrating a calculated error curve of a critical time-frequency ridge line of an IEIMD according to an embodiment of the present invention;

FIG. 4(a) is a partial enlarged view of a waveform of a critical mode of an emulation signal IEIMD according to an embodiment of the present invention;

FIG. 4(b) is a diagram of a steady-state angular domain of a simulation signal windowing in an embodiment of the present invention;

FIG. 5(a) is a chart of a simulation signal before updating in an embodiment of the present invention;

FIG. 5(b) is a diagram of an updated order spectrogram and a local enlarged view of a simulation signal according to an embodiment of the present invention;

FIG. 6 is a schematic structural diagram of a test bed according to an embodiment of the present invention;

FIG. 7(a) is a frequency spectrum and a partially enlarged view of a test bed signal according to an embodiment of the present invention;

FIG. 7(b) is a diagram illustrating a critical time-frequency ridge extraction of a test-bed signal IEIMD in an embodiment of the present invention;

FIG. 8 is a partial enlarged view of a waveform of a key mode of a test bed signal IEIMD in accordance with an embodiment of the present invention;

FIG. 9(a) is a chart of the test bed signal before updating in an embodiment of the present invention;

FIG. 9(b) is a diagram of an updated order spectrum and a partially enlarged view of a test bed signal according to an embodiment of the present invention.

Detailed Description

In order to make the objects, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are some, but not all, embodiments of the present invention. The components of embodiments of the present invention generally described and illustrated in the figures herein may be arranged and designed in a wide variety of different configurations.

Thus, the following detailed description of the embodiments of the present invention, presented in the figures, is not intended to limit the scope of the invention, as claimed, but is merely representative of selected embodiments of the invention. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.

It should be noted that: like reference numbers and letters refer to like items in the following figures, and thus, once an item is defined in one figure, it need not be further defined and explained in subsequent figures.

As shown in fig. 1, the present invention provides a rotating machine fault diagnosis method based on iterative expansion eigenmode decomposition, assuming that a variable rotation speed vibration signal is x (t), the method specifically includes the following steps:

step 1: installing a vibration acceleration sensor near the rotating machinery, and acquiring a vibration acceleration signal x (t) when the rotating speed is changed to work;

step 2: performing discrete short-time Fourier transform on the collected vibration acceleration signals x (t) to obtain a time-frequency expression, and specifically comprising the following steps:

step 201: setting the Window Length N_windowAnd the number of Fourier points L_fourierConstructing a corresponding window function wf (t);

step 202: and performing discrete short-time Fourier transform (STFT) on the vibration acceleration signal x (t) through the window function wf (t) to obtain a time-frequency expression TFR (t, f).

And step 3: preliminarily estimating the key time-frequency ridge line with the most significant energy in the time-frequency expression by a bidirectional rapid single-line extraction method, and specifically comprising the following steps of:

step 301: initializing a parameter maximum frequency difference value delta f;

step 302: locating an energy maximum point (t) in the TFR (t, f)_em,f_em)：

In the formula, TFR (t, f) is a time-frequency expression,is the starting point of the initial estimated value of the critical time-frequency ridge line;

and searching an energy maximum value point of a frequency dimension in a delta f range in a two-way mode along the time direction:

f_R＝f_em,f_L＝f_em

in the formula (I), the compound is shown in the specification,the initial estimation value N of the key time-frequency ridge line is the number of sampling points.

Step 303: obtaining the initial estimation value of the critical time-frequency ridge line

And 4, step 4: separating and reconstructing a key mode, an instantaneous frequency and an instantaneous amplitude thereof according to the vibration acceleration signal and a key time-frequency ridge line initial value by iterative expansion eigenmode decomposition, and specifically comprising the following steps:

step 401: initializing a regularization parameter lambda, punishment parameters beta and mu and a Fourier expansion order L, and inputting a vibration acceleration signal x (t) and an initial estimation value of a critical time-frequency ridge line

Wherein BW is the bandwidth, N is the number of sampling points, f_sIs the sampling frequency;

step 402: performing Hilbert transform on the x (t) to construct a key mode m_key(t) analytical formula

Where k is 0.. end, end is the number of iterations,is the k iterations of the key instantaneous frequency, η (t) is the interference noise, t is the discrete sampling time,to resolve the amplitude:

wherein the content of the first and second substances,in order to be the initial phase position,is the true value of the critical instantaneous frequency, a_key(t) is the instantaneous amplitude of the key mode; then, the analytic amplitude is establishedRedundant fourier expansion of (a):

wherein the content of the first and second substances, l in (1) is the order of the Fourier series;

step 403: will be described inUnfolding model bring-inConstructing a matrix expression:

K^k＝E^kt isN × (2L +1) kernel matrix:

in the formula, d is a matrix index position;

the optimization problem is then solved:

in the formula (I), the compound is shown in the specification,is a signal analytic formula;

obtaining a Fourier parameter matrixAnd is reconstituted to obtain

Wherein, I is an identity matrix, and H is a matrix transpose;

step 404: unfolding the reconstruction based on the Euler formulaComprises the following steps:

getReal part and imaginary part, and estimating error by arc tangent demodulation method

In the formula, R^kIs the real part of the complex envelope, I^kIs the imaginary part of the complex envelope.

To pairSmoothing the filter and updating the instantaneous frequency

WhereinTo improve the second order difference operator;

step 405: will be described inJump to step 403 to update kernel matrix K^k+1Parameter matrixAnd resolving the magnitude matrixRepeatedly calculating end times and stopping, and then solving and reconstructing the key mode m_key(t) and its amplitude a_key(t)：

m_key＝u·A+v·B

Wherein the content of the first and second substances,

and 5: performing equal-angle resampling on the key mode to obtain a steady-state angular domain signal, and specifically comprising the following steps:

step 501: for the key mode m_key(t) performing Hilbert transform, arctangent of the ratio of the imaginary part to the real part, and deconvoluting to obtain an instantaneous phase curve

Step 502: establishing the instantaneous phase curve by least squares fittingTime-radian mapping relation psi_key(t)；

Step 503: determining a sampling order O_mAnd equal radian interval delta theta are summed, and the mapping relation is reversely solvedObtaining the equal radian time sequence t_Δθ(k) And radian sequence rad_Δθ(k)；

Step 504: for the time sequence t_Δθ(k) And the vibration acceleration signal x (t) is subjected to cubic spline interpolation to obtain an equal radian amplitude sequence a_Δθ(k)；

Step 6: converting the angular domain into an order domain through discrete Fourier transform, and analyzing each characteristic order, wherein the method specifically comprises the following steps:

step 601: for the equal radian amplitude sequence a_Δθ(k) Jiahanning window w_hanning(k) According to said sampling order O_mCalculating order sequence O (k), calculating amplitude sequence A by discrete Fourier transform_Δθ(k) Transition to the order domain;

step 602: in the order domain, if the first order is the significant energy and the lowest order, directly analyzing each feature order proportional to the first order, including: first order and its multiple order, fault signature order and its multiple order;

step 603: in the order domain, if the first order is not significant in energy and is the lowest order, locating the significant and lowest order of energy smaller than the first order, calculating the reciprocal of the order to obtain an amplification factor m, amplifying the order sequence O (k) by m times in an equal proportion, and analyzing each characteristic order according to the step 602 to obtain a fault diagnosis result;

example 1

As shown in FIG. 2(a), the signal-to-noise ratio of the gearbox vibration simulation signal x (t) is-5 dB, the sampling frequency is 4000Hz, and the sampling time is 5 s; as shown in fig. 2(b), serious "spectrum aliasing" occurs in the spectrum x (t), and the failure information cannot be analyzed. As shown in fig. 3(a), the frequency conversion is increased from 400Hz to about 1100Hz within 2.5 seconds, the instantaneous frequency change rate is greater than 30%, the characteristic of obvious large rotation speed fluctuation is achieved, and due to cross term limitation and noise interference, the side frequency band and the time-frequency ridge lines at two ends in the time-frequency spectrogram are fuzzy; as shown in FIG. 3(c), the key time-frequency ridge extracted by IEIMDThe matching degree with the real ridge line is extremely high;as shown in figure 3(d) of the drawings,the calculation error curve is stable and within 0.1 percent, and the deviation of two ends of the Vold-Kalman filtering error curve is large. The Mean Square Error (MSE) and the Relative Error (RE) are shown in table 1, and compared with the similar method, the instantaneous frequency extraction precision is improved by 21%, and especially compared with Vold-kalman filtering, the extraction precision is improved by nearly 3 times.

TABLE 1 MSE and RE indices

As shown in FIG. 4(a), the key modality m obtained by IEIMD reconstruction_key(t) and its instantaneous amplitude a_keyAnd (t) the waveform has good quality and obvious physical characteristics. As shown in Table 2, the reconstruction signal-to-noise ratio (SNR) is improved from original-5 dB to 6.35dB, the interference of random noise is effectively inhibited, compared with the similar method, the modal reconstruction precision is improved by 14%, and especially, compared with the Vold-Kalman filtering, the modal reconstruction precision is improved by 43%.

TABLE 2 SNR index

As shown in fig. 5(a), the energy significant lowest order less than one order is located, and the reciprocal is found as:

key instantaneous frequencyThe second meshing instantaneous frequency is obtained by amplifying the order sequence O (k) by 35.91 timesA standard order spectrogram is obtained, and as shown in fig. 5(b), a significant side band is present around the engagement feature order, and the effective information is abundant. The calculation error of the characteristic order is shown in Table 3, the calculation error of the method is within 0.5%, compared with the similar method, the calculation precision is improved by 19%, and especially compared with the Vold-Kalman filtering, the calculation precision is improved by 40%.

TABLE 3 characteristic order calculation error

Software and hardware environment in the test process, hardware environment: the CPU is Intel (R) core (TM) i5-6300HQ, the main frequency is 2.30GHz, the memory is 8.00GB, and the 64-bit operating system is adopted. Software environment: windows 10 enterprise edition, MATLAB R2018 a. The calculation consumption is shown in table 4, and the calculation speed is improved by about 2 times compared with the similar method.

TABLE 4 calculated consumptions

In conclusion, compared with the similar method, various performance indexes of the method are improved remarkably, and the method is more suitable for the working condition with large rotating speed fluctuation compared with Vold-Kalman order tracking.

Example 2

As shown in fig. 6, the model of the MSF-PK5M mechanical fault simulation test bed is 623C01, the models of the rolling bearings are ER16K, the specific parameters are shown in table 5, and the inner ring fault bearing is selected in the test.

TABLE 5 bearing parameters

As shown in fig. 7(a), the resolution of the time-frequency spectrogram is low due to cross term limitation and noise interference. As shown in FIG. 8, the key modality m obtained by IEIMD reconstruction_key(t) and its amplitude a_keyAnd (t) the waveform is clear, the physical property is good, and the fault deep analysis is facilitated. As shown in fig. 9(a), the energy significant lowest order less than one order is located, and the reciprocal is found as:

key instantaneous frequencyThat is, the inner ring fault instantaneous frequency is obtained by amplifying the order sequence o (k) by 5.42 times to obtain a standard order spectrogram, as shown in fig. 9(b), the fault characteristic order is in a cluster shape, the 2 times and 3 times of orders are significant in energy, and the fault information is rich.

The calculation error of the characteristic order is shown in table 6, the calculation error of the method is within 0.2%, and compared with the similar method, the calculation precision is improved by 38%, and especially compared with COT, the calculation precision is improved by nearly 82%. This is because the COT calculation accuracy is affected by the quality of the rotation speed signal, and therefore the calculation error is the largest.

TABLE 6 characteristic order calculation error

Software and hardware environment in the test process, hardware environment: the CPU is Intel (R) core (TM) i5-6300HQ, the main frequency is 2.30GHz, the memory is 8.00GB, and the 64-bit operating system is adopted. Software environment: windows 10 enterprise edition, MATLAB R2018 a. The calculation consumption of the method is shown in table 7, and compared with the similar method, the calculation speed is improved by about 3 times.

TABLE 7 calculated consumptions

In conclusion, the rotary machine fault diagnosis method based on the keyless phase order tracking technology of the iterative expansion eigenmode decomposition has remarkable advantages, can realize efficient signal processing and fault diagnosis for mechanical equipment with variable rotating speed, particularly large rotating speed fluctuation working conditions, and is beneficial to embedded development and practical engineering application.

The foregoing detailed description of the preferred embodiments of the invention has been presented. It should be understood that numerous modifications and variations could be devised by those skilled in the art in light of the present teachings without departing from the inventive concepts. Therefore, the technical solutions available to those skilled in the art through logic analysis, reasoning and limited experiments based on the prior art according to the concept of the present invention should be within the scope of protection defined by the claims.