Current signal iteration Vold-Kalman filtering detection method for rotating machine fault

Abstract

The invention provides a current signal iteration Vold-Kalman filtering detection method for a rotating machinery fault, and belongs to the technical field of rotating machinery, and the method comprises the steps: carrying out the coordinate transformation of a three-phase stator current during the operation of the rotating machinery, and obtaining Id; obtaining a time-frequency spectrum of the original current, constructing a cost function, and extracting an instantaneous frequency ridge line of the fundamental wave of the original current; according to the instantaneous frequency ridge line and the component design parameters, calculating the characteristic frequency when each component breaks down; an iterative Vold-Kalman filter is constructed, the calculated fault frequency and the frequency multiplication of the fault frequency are taken as cut-off frequencies, and fault harmonics in the current signal Id are iteratively checked and extracted through a substitution check method; and carrying out time-frequency analysis on the extracted fault harmonic waves, and identifying the fault type. According to the method, a non-intrusive method is adopted, the non-stationary characteristic of the signal is considered, the calculation complexity is reduced, meanwhile, the fault harmonic extraction precision is improved, and the rotating machine fault can be effectively detected through the current signal under the variable rotating speed working condition.

Description

Technical Field

[0001] The invention belongs to the technical field of rotating machinery, and in particular relates to a current signal iterative Vold-Kalman filtering detection method for rotating machinery faults. Background Art

[0002] Rotating machinery is a key equipment in industrial production and is widely used in high-speed trains, wind power generation, ships, automobiles and other fields. Rotating machinery is often subjected to complex dynamic loads during operation, and is therefore prone to failure. If the failure is not discovered, it is very likely to spread to the matching parts, and even lead to long-term downtime and huge economic losses. Taking the failure data of various systems of Guangzhou Metro trains over a three-year period as an example, running gear failures accounted for 35% of the total number of failures, while rotating machinery failures accounted for more than 78% of running gear failures. Therefore, considering complex working conditions, service environment and other conditions, it is crucial to accurately detect rotating machinery failures and perform maintenance in a timely manner to ensure the safe and reliable operation of equipment.

[0003] Among the current iterative Vold-Kalman filter detection methods for rotating machinery faults, vibration signal analysis is the most common. However, the transmission path from the fault excitation point to the measurement point is complex, and the vibration signal will produce modulation, which increases the difficulty of fault detection. In addition, in actual application scenarios, acceleration sensors are expensive and inconvenient to install. In order to overcome the above-mentioned defects in vibration signal analysis, motor current signal analysis (MCSA) has gradually been used for rotating machinery fault diagnosis. As a non-invasive method, fault diagnosis based on MCSA has obvious advantages: 1. Current signal acquisition will not interfere with system operation; 2. Current signals are not affected by the transmission path; 3. Current sensors are inexpensive and easy to install.

[0004] Although many scholars have proposed a current signal iterative Vold-Kalman filter detection method for rotating machinery faults based on current signals, the existing methods have problems such as low time-frequency resolution, high computational complexity, and low harmonic extraction accuracy. At present, domestic invention patent applications for rotating machinery fault diagnosis using current signals include:

[0005] A Chinese patent with the authorization date of June 20, 2023 and the authorization number CN111680665B discloses an invention patent entitled "A data-driven method for diagnosing motor mechanical faults using current signals". It first reconstructs the stator current signal as a noise signal, then subtracts the reconstructed signal from the original stator current signal to obtain the residual current, and then converts the residual current to the frequency domain and inputs it into a trained machine learning model, effectively identifying the motor state. However, this method requires pre-training of the machine learning model, which takes a long time to train and cannot cover all fault types.

[0006] A Chinese patent with the authorization date of May 3, 2024 and the authorization number CN117112984B discloses an invention patent entitled "A non-invasive lathe belt drive system fault diagnosis method". It uses the square modulation bispectrum to estimate the amplitude and phase of the rotation frequency sideband, effectively realizing the health status monitoring of the lathe spindle motor drive system, but it has the defects of high computational complexity and inability to quickly detect faults.

[0007] The Chinese patent with the publication date of August 2, 2024 and the publication number of CN118425769A discloses an invention patent named "A current signal processing method, a motor operation fault detection method and device". It uses the original current to calculate the instantaneous root mean square value and perform weighted average calculation, which can remove the power frequency carrier signal while retaining the fault characteristic signal. However, there are many types of rotating machinery faults. This method cannot accurately diagnose the type of fault only through time domain analysis, and the diagnostic accuracy needs to be improved. Summary of the invention

[0008] The object of the present invention is to provide a current signal iterative Vold-Kalman filtering detection method for rotating machinery faults, characterized in that it comprises the following steps:

[0009] S1: Acquire the three-phase stator current signal of the driving motor of the rotating machinery when the rotating machinery is running at a variable speed;

[0010] S2: Perform Clark transformation and Park transformation on the three-phase stator current signal to obtain the current signal I in the two-phase rotating coordinate system d ;

[0011] S3: Perform time-frequency analysis on the single-phase current signal of S1 based on short-time Fourier transform to obtain time-frequency energy;

[0012] S4: construct a cost function using time-frequency energy to calculate the instantaneous frequency ridge of the fundamental wave of the current signal;

[0013] S5: Calculate the characteristic frequency f of each component when it fails based on the instantaneous frequency ridge and component design parameters fault ;

[0014] S6: Construct an iterative Vold-Kalman filter, use the fault frequency and its multiples in S5 as the cutoff frequency, and iteratively test and extract the current signal I through the substitution test method d Fault harmonics in

[0015] S7: Perform time-frequency analysis on the extracted fault harmonics to identify the type of rotating machinery fault.

[0016] Furthermore, S2 specifically includes the following steps:

[0017] S21: Clark transformation is performed on the stator current signal collected in the three-phase stationary coordinate system to obtain the current signal in the two-phase stationary coordinate system. The solution formula is expressed as:

[0018]

[0019] Among them, I sa ,I sb and I sc is the stator current signal collected in the three-phase stationary coordinate system; I α and I β is the stator current signal in the two-phase stationary coordinate system;

[0020] S22: stator current signal I in the two-phase stationary coordinate system α ,I β Perform Park transformation to obtain the stator current signal in the two-phase rotating coordinate system. The solution formula is expressed as:

[0021]

[0022] Among them, I d and I q is the stator current signal in the two-phase rotating coordinate system.

[0023] Furthermore, S3 specifically includes the following steps:

[0024] S31: Set the window width W and the number of Fourier points N to construct the short-time Fourier transform function STFT(t);

[0025] S32: Perform discrete short-time Fourier transform on a single-phase stator current signal through a function to obtain time-frequency energy E(t,f).

[0026] Furthermore, S4 specifically includes the following steps:

[0027] S41: Calculate the maximum value of the time-frequency energy in S32, and denote the time and frequency here as t base and f base ;

[0028] S42: With t base is the initial moment, and the cost function is constructed along the time axis to both sides. The expression of the cost function at each moment is expressed as:

[0029]

[0030] Among them, E(t i ,f i ) is the time spectrum t i At the moment iThe energy at the frequency; J is the sampling frequency of the output data; c1 is the coefficient of the time-frequency energy; c2 is the coefficient of the frequency change rate;

[0031] S43: Calculate the minimum value of the cost function at each moment, and use the frequency corresponding to the cost function value as the instantaneous frequency at that moment;

[0032] S44: After calculating the instantaneous frequency at a moment, the instantaneous frequency is recorded as f base , and repeat S43 until the instantaneous frequency of all moments is obtained, and the instantaneous frequencies of all moments are connected in time series, which is the instantaneous frequency ridge line f of several waves of the current signal s .

[0033] Furthermore, S5 specifically includes the following steps:

[0034] S51: Calculate the motor rotation frequency according to the instantaneous frequency ridge line and the number of motor pole pairs. The calculation formula is expressed as:

[0035]

[0036] Among them, f s is the instantaneous frequency ridge of the current fundamental wave; p is the number of motor pole pairs; f r is the motor speed frequency;

[0037] S52: Calculate multiple characteristic frequencies f when each component fails based on the motor frequency and the design parameters of the rotating machinery fault .

[0038] Furthermore, S6 specifically includes the following parameters:

[0039] S61: Construct a Vold-Kalman filter. The state equation of the Vold-Kalman filter is expressed as:

[0040] DA k =ε k

[0041]

[0042] Where D is a second-order n*n subtraction operator; A k is the envelope amplitude matrix of the signal; ε k is a non-homogeneous term;

[0043] The observation equation of the Vold-Kalman filter is expressed as:

[0044] X=A k C k +δ k

[0045] Where X is the n*1 observed current signal matrix; δ k is the n*1 generalized noise matrix; C k =diag[θ k(1) ,θ k(2) ,θ k(3) ,...,θ k(n) ] is the target component X k The carrier matrix of

[0046] According to the state equation, observation equation and the pre-set carrier matrix C k , define the loss function as the non-homogeneous term ε k and generalized noise δ k The weighted energy and loss function are expressed as follows:

[0047]

[0048] Where α is the weighting factor, (·)H represents the conjugate transpose of the matrix;

[0049] When the loss function is minimized, When the target component X is obtained k The amplitude envelope, amplitude envelope A k The expression is expressed as:

[0050]

[0051] Where I is the identity matrix;

[0052] According to the carrier matrix C k and amplitude envelope A k The target component X to be extracted can be reconstructed k , the reconstruction method is as follows:

[0053] X k =A k C k

[0054] S62: Multiple characteristic frequencies f in S52 fault As the cutoff frequency of the Vold-Kalman filter, the current signal I in S2 is d Extract multiple target components (X1, X2, ..., X k );

[0055] S63: testing multiple target components to be fault harmonics or noises by a substitution test method, and determining the null hypothesis of the target component X1 - X1 is noise;

[0056] S64: Generate m groups of substitute data of the target component X1 by iterative amplitude adjustment Fourier transform, and the power spectrum and amplitude distribution of the substitute data are consistent with the target component X1;

[0057] S65: Calculate the envelope energy entropy Q[A] and instantaneous frequency energy entropy Q[F] of the target component and the m groups of surrogate data respectively, adaptively set different weights for Q[A] and Q[F], and construct a statistic of the surrogate test. The calculation formula of the test statistic is expressed as:

[0058]

[0059] Among them, D(ζ a ,ζ f ) is the test statistic; f(x) is the target component or alternative data; Q[f(x)] is the energy entropy of f(x), which is used to measure the order of data. The larger the energy entropy, the greater the uncertainty of the data and the more chaotic the distribution; ζ a represents the weight of Q[A]; ζ f represents the weight of Q[F]; A(X K ) represents the Fourier coefficient of the target component envelope; F(X K ) represents the Fourier coefficient of the instantaneous frequency ridge of the target component; the weights are based on the target component X k Calculation, when the target component changes, the weight of the test statistic is adaptively adjusted;

[0060] S66: When the test statistic of the target component is less than the test statistic of all alternative data, the target component is a fault harmonic, and the null hypothesis is rejected and the fault harmonic is saved; when the test statistic of the target component is greater than or equal to the test statistic of all alternative data, the target component is noise, and the null hypothesis is accepted and the target component is discarded;

[0061] S67: Repeat S63 to S66 to iteratively extract all fault harmonics in the target component in S62.

[0062] Furthermore, S7 specifically includes the following steps:

[0063] S71: Perform time-frequency analysis on the fault harmonics in S67 by short-time Fourier transform;

[0064] S72: Compare the time-frequency analysis result with the fault characteristic frequency of each component of the rotating machinery, so as to identify the fault type of the rotating machinery.

[0065] Compared with the prior art, the beneficial effects of the present invention are mainly reflected in:

[0066] 1. The present invention adopts a non-invasive method, takes into account the non-stationary characteristics of the signal, and can use the current signal to effectively detect rotating machinery faults under variable speed conditions. Compared with the detection method based on vibration signals, the sensor installation scheme of the present invention is easier to implement, has lower hardware costs, and improves the practicality of intelligent health management.

[0067] 2. Each component of the rotating machinery may generate multiple fault harmonics. The present invention does not need to traverse all the fault harmonics of each component in the current signal. By extracting the fault harmonics of other components immediately after extracting the fault harmonics of a certain component, the number of fault harmonics that need to be extracted is reduced, thereby reducing the calculation complexity.

[0068] 3. The present invention proposes a weight-adaptive test statistic, which can accurately identify harmonics even when the signal-to-noise ratio is low, effectively suppresses noise interference, and improves the detection accuracy of fault harmonics. BRIEF DESCRIPTION OF THE DRAWINGS

[0069] Figure 1 It is a schematic diagram of the process of the present invention.

[0070] Figure 2 Schematic diagram of experimental motor speed.

[0071] Figure 3 For missing tooth fault d Time-frequency diagram of fault harmonics.

[0072] Figure 4 For broken tooth fault d Time-frequency diagram of fault harmonics.

[0073] Figure 5 Schematic diagram for comparing fault harmonic extraction methods under missing tooth fault.

[0074] Figure 6 Schematic diagram of the comparison of fault harmonic extraction methods under broken tooth fault.

[0075] Among them, (a) is the working condition 1 of each figure, (b) is the working condition 2 of each figure, (c) is the working condition 3 of each figure, and (d) is the working condition 4 of each figure. DETAILED DESCRIPTION

[0076] The iterative Vold-Kalman filtering detection method for current signals of rotating machinery faults of the present invention will be described in more detail below in conjunction with a schematic diagram, wherein a preferred embodiment of the present invention is shown. It should be understood that those skilled in the art can modify the present invention described herein while still achieving the beneficial effects of the present invention. Therefore, the following description should be understood as being widely known to those skilled in the art and not as a limitation to the present invention.

[0077] like Figure 1As shown, the current signal iterative Vold-Kalman filtering detection method for rotating machinery faults includes the following steps:

[0078] Step 1: Obtain the three-phase stator current signal of the driving motor when the rotating machinery is running at a variable speed.

[0079] Step 2: Perform Clark transformation and Park transformation on the three-phase stator current signal to obtain the current signal I in the two-phase rotating coordinate system d .

[0080] (1) Perform Clark transformation on the stator current signal collected in the three-phase stationary coordinate system to obtain the current signal in the two-phase stationary coordinate system. The solution formula is expressed as:

[0081]

[0082] Among them, I sa ,I sb and I sc is the stator current signal collected in the three-phase stationary coordinate system; I α and I β is the stator current signal in the two-phase stationary coordinate system.

[0083] (2) For the stator current signal I in the two-phase stationary coordinate system α ,I β Perform Park transformation to obtain the stator current signal in the two-phase rotating coordinate system. The solution formula is expressed as:

[0084]

[0085] Among them, I d and I q is the stator current signal in the two-phase rotating coordinate system.

[0086] Step 3: Perform time-frequency analysis on the single-phase current signal of S1 based on short-time Fourier transform to obtain time-frequency energy.

[0087] (1) Set the window width W and the number of Fourier points N to construct the short-time Fourier transform function STFT(t).

[0088] (2) The single-phase stator current signal is subjected to discrete short-time Fourier transform through the function to obtain the time-frequency energy E(t,f).

[0089] Step 4: Use time-frequency energy to construct a cost function and calculate the instantaneous frequency ridge of the fundamental wave of the current signal.

[0090] (1) Calculate the maximum value of the time-frequency energy in step 3(2), and denote the time and frequency here as t base and fbase .

[0091] (2) With t base is the initial moment, and the cost function is constructed along the time axis to both sides. The expression of the cost function at each moment is expressed as:

[0092]

[0093] Among them, E(t i ,f i ) is the time spectrum t i At the moment i The energy at the frequency; J is the sampling frequency of the output data; c1 is the coefficient of time-frequency energy; c2 is the coefficient of the frequency change rate.

[0094] (3) Calculate the minimum value of the cost function at each moment, and take the frequency corresponding to the cost function value as the instantaneous frequency at that moment.

[0095] (4) After calculating the instantaneous frequency at a certain moment, the instantaneous frequency is recorded as f base , and repeat S43 until the instantaneous frequency of all moments is obtained, and the instantaneous frequencies of all moments are connected in time series, which is the instantaneous frequency ridge line f of several waves of the current signal s .

[0096] Step 5: Calculate the characteristic frequency f of each component when it fails based on the instantaneous frequency ridge and component design parameters fault , specifically including the following steps:.

[0097] (1) The motor rotation frequency is calculated based on the instantaneous frequency ridge line and the number of motor pole pairs. The calculation formula is expressed as:

[0098]

[0099] Among them, f s is the instantaneous frequency ridge of the current fundamental wave; p is the number of motor pole pairs; f r Motor speed

[0100] (2) Based on the motor frequency and the design parameters of the rotating machinery, calculate the multiple characteristic frequencies f when each component fails fault .

[0101] Step 6: Construct an iterative Vold-Kalman filter, use the fault frequency in S5 and its multiples as the cutoff frequency, and iteratively test and extract the current signal I through the substitution test method d Fault harmonics in the

[0102] (1) Construct a Vold-Kalman filter. The state equation of the Vold-Kalman filter is expressed as:

[0103]

[0104] Where D is a second-order n*n subtraction operator; A k is the envelope amplitude matrix of the signal; ε k is a non-homogeneous term;

[0105] The observation equation of the Vold-Kalman filter is expressed as:

[0106] X=A k C k +δ k

[0107] Where X is the n*1 observed current signal matrix; δ k is the n*1 generalized noise matrix; C k =diag[θ k(1) ,θ k(2) ,θ k(3) ,...,θ k(n) ] is the target component X k The carrier matrix of

[0108] According to the state equation, observation equation and the pre-set carrier matrix C k , define the loss function as the non-homogeneous term ε k and generalized noise δ k The weighted energy and loss function are expressed as follows:

[0109]

[0110] Where α is the weighting factor, (·)H represents the conjugate transpose of the matrix;

[0111] When the loss function is minimized, When the target component X is obtained k The amplitude envelope, amplitude envelope A k The expression is expressed as:

[0112]

[0113] Where I is the identity matrix;

[0114] According to the carrier matrix C k and amplitude envelope A k The target component X to be extracted can be reconstructed k , the reconstruction method is as follows:

[0115] X k =A k C k

[0116] (2) Multiple characteristic frequencies f in S52 fault As the cutoff frequency of the Vold-Kalman filter, the current signal I in S2 is d Extract multiple target components (X1, X2, ..., X k ).

[0117] (3) Through the substitution test method, multiple target components are tested to see whether they are fault harmonics or noise, and the null hypothesis of the target component X1 is determined - X1 is noise.

[0118] (4) Generate m groups of surrogate data of the target component X1 through iterative amplitude adjustment Fourier transform. The power spectrum and amplitude distribution of the surrogate data are consistent with the target component X1.

[0119] (5) Calculate the envelope energy entropy Q[A] and instantaneous frequency energy entropy Q[F] of the target component and m groups of surrogate data respectively, set different weights for Q[A] and Q[F] adaptively, and construct the statistic of the surrogate test. The calculation formula of the test statistic is expressed as:

[0120]

[0121] Among them, D(ζ a ,ζ f ) is the test statistic; f(x) is the target component or alternative data; Q[f(x)] is the energy entropy of f(x), which is used to measure the order of data. The larger the energy entropy, the greater the uncertainty of the data and the more chaotic the distribution; ζ a represents the weight of Q[A]; ζ f represents the weight of Q[F]; A(X K ) represents the Fourier coefficient of the target component envelope; F(X K ) represents the Fourier coefficient of the instantaneous frequency ridge of the target component; the weights are based on the target component X k Calculate, and when the target component changes, the weight of the test statistic is adaptively adjusted.

[0122] (6) When the test statistic of the target component is less than the test statistic of all alternative data, the target component is a fault harmonic, and the null hypothesis is rejected and the fault harmonic is saved; when the test statistic of the target component is greater than or equal to the test statistic of all alternative data, the target component is noise, and the null hypothesis is accepted and the target component is discarded.

[0123] (7) S63 to S66 are repeatedly executed to iteratively extract all fault harmonics in the target component in S62.

[0124] Step 7: Perform time-frequency analysis on the extracted fault harmonics to identify the type of rotating machinery fault.

[0125] (1) Time-frequency analysis of the fault harmonics in S67 is performed through short-time Fourier transform.

[0126] (2) Compare the time-frequency analysis results with the characteristic frequencies of the faults of various components of the rotating machinery to identify the fault type of the rotating machinery.

[0127] The present invention will be further described below through specific examples.

[0128] A mechanical fault simulation test bench, the current clamp model is iProbe 100, the drive motor is an induction motor, the gearbox is a two-stage parallel shaft gearbox, and the parameters of the motor and the gearbox are shown in Table 1 and Table 2 respectively. In this experiment, the fault is injected into the pinion 2, including two types of faults: missing teeth and broken teeth. Four speed change conditions are set in the experiment, including: (a) acceleration (from 1200rpm to 3000rpm); (b) deceleration (from 3000rpm to 1200rpm); (c) acceleration-uniform speed-deceleration; (d) acceleration-uniform speed-deceleration-uniform speed-acceleration-uniform speed-deceleration. The actual speeds of the four conditions are as follows Figure 2 shown.

[0129]

[0130] Table 1 Induction motor design parameters

[0131]

[0132] Table 2 Gearbox parameters

[0133] After the pinion 2 has a missing tooth fault, the proposed iterative Vold-Kalman filter (IVKF) is used to extract the current signal I d Extract the fault harmonics from the Figure 3 As shown. Figure 3 It can be seen that under condition 1, IVKF can effectively detect I d There is a fault frequency of pinion 2 (f p2 ) of 10 times, 11 times, 12 times, and 13 times; Under working condition 2, IVKF can detect I d There is f p2 9 times, 10 times, and 11 times of the frequency; Under working condition 3, IVKF can detect I d There is f p2 8 times and 11 times of the frequency; Under working condition 4, IVKF can detect I d There is f p2 11 times the frequency.

[0134] Since the missing tooth fault is a more serious fault, the fault harmonics in the current signal are more obvious. In order to explore whether the method proposed in the present invention is still effective under weaker faults, the current data collected under the broken tooth fault is used for verification. Figure 4 The fault harmonic extraction results after the broken tooth fault of pinion 2 are given. Figure 4 It can be seen that under working conditions 1 and 2, IVKF can effectively detect I d There is a fault frequency of pinion 2 (f p2 ) of 10 times and 13 times. Under working conditions 3 and 4, IVKF can detect I d There is f p2 10 times the frequency.

[0135] from Figure 3 and Figure 4 The harmonic extraction results show that: 1. In the case of gear failure, the current I d 1. The proposed IVKF can effectively extract the gear fault harmonics in the current signal, eliminating the influence of speed change, whether it is a simple speed change or a complex speed fluctuation condition; 2. The proposed IVKF can effectively extract the gear fault harmonics in the current signal, eliminating the influence of speed change, whether it is a missing tooth fault with obvious fault characteristics, or a broken tooth fault with weaker fault characteristics. d The time-varying gear fault characteristics can be effectively extracted.

[0136] Table 3 Comparison of the accuracy of three harmonic extraction methods under four working conditions

[0137]

[0138] In order to show the superiority of this method, adaptive chirped mode decomposition (ACMD) and variational nonlinear chirped mode decomposition (VNCMD) are used to compare with this method. d The 11-fold frequency extraction result of the fault frequency of small gear 2 is as follows: Figure 5 As shown in the figure; three methods are used to fix the broken tooth fault d The 10-fold frequency extraction result of the fault frequency of small gear 2 is as follows: Figure 6 As shown in the figure, it can be seen that compared with ACMD and VNCMD, the fault harmonics extracted by the proposed IVKF are closer to the true value. The logarithmic root mean square error (LMSE) is used to quantitatively evaluate the harmonic extraction accuracy of the three methods, and the results are shown in Table 3. It can be seen from the table that whether it is a missing tooth fault or a weaker broken tooth fault, the fault harmonic time-frequency ridge extracted by IVKF is more accurate than the other two methods, and has a stronger ability to suppress noise.

[0139] In summary, the present invention can detect rotating machinery faults by using current signals, and has the advantages of convenient data acquisition, low computational complexity, and high harmonic extraction accuracy. Compared with similar methods, the fault diagnosis accuracy of this method is significantly improved.

[0140] The above is only a preferred embodiment of the present invention and does not limit the present invention in any way. Any technician in the relevant technical field, without departing from the scope of the technical solution of the present invention, makes any form of equivalent replacement or modification to the technical solution and technical content disclosed in the present invention, which does not depart from the content of the technical solution of the present invention and still falls within the protection scope of the present invention.

Claims

1. A current signal iterative Vold-Kalman filtering detection method for rotating machinery faults, characterized in that: The steps include: S1: Acquire the three-phase stator current signal of the driving motor of the rotating machinery when the rotating machinery is running at a variable speed; S2: Perform Clark transformation and Park transformation on the three-phase stator current signal to obtain the current signal I in the two-phase rotating coordinate system d ; S3: Perform time-frequency analysis on the single-phase current signal of S1 based on short-time Fourier transform to obtain time-frequency energy; S4: construct a cost function using time-frequency energy to calculate the instantaneous frequency ridge of the fundamental wave of the current signal; S5: Calculate the characteristic frequency f of each component when it fails based on the instantaneous frequency ridge and component design parameters fault ; S6: Construct an iterative Vold-Kalman filter, use the fault frequency and its multiples in S5 as the cutoff frequency, and iteratively test and extract the current signal I through the substitution test method d Fault harmonics in S7: Perform time-frequency analysis on the extracted fault harmonics to identify the type of rotating machinery fault.

2. The current signal iterative Vold-Kalman filtering detection method for rotating machinery fault according to claim 1, characterized in that: S2 specifically includes the following steps: S21: Perform Clark transformation on the stator current signal collected in the three-phase stationary coordinate system to obtain the stator current signal I in the two-phase stationary coordinate system α and I β ; S22: stator current signal I in the two-phase stationary coordinate system α ,I β Perform Park transformation to obtain the stator current signal I in the two-phase rotating coordinate system d and I q .

3. The current signal iterative Vold-Kalman filtering detection method for rotating machinery fault according to claim 1, characterized in that: S3 specifically includes the following steps: S31: Set the window width W and the number of Fourier points N to construct the short-time Fourier transform function STFT(t); S32: Perform discrete short-time Fourier transform on a single-phase stator current signal through a function to obtain time-frequency energy E(t,f).

4. The current signal iterative Vold-Kalman filtering detection method for rotating machinery fault according to claim 1, characterized in that: S4 specifically includes the following steps: S41: Calculate the maximum value of the time-frequency energy in S32, and denote the time and frequency here as t base and f base ; S42: With t base is the initial moment, and the cost function is constructed along the time axis to both sides. The expression of the cost function at each moment is expressed as: Among them, E(t i ,f i ) is the time spectrum t i At the moment i The energy at the frequency; J is the sampling frequency of the output data; c1 is the coefficient of the time-frequency energy; c2 is the coefficient of the frequency change rate; S43: Calculate the minimum value of the cost function at each moment, and use the frequency corresponding to the cost function value as the instantaneous frequency at that moment; S44: After calculating the instantaneous frequency at a moment, the instantaneous frequency is recorded as f base , and repeat S43 until the instantaneous frequency of all moments is obtained, and the instantaneous frequencies of all moments are connected in time series, which is the instantaneous frequency ridge line f of several waves of the current signal s .

5. The current signal iterative Vold-Kalman filtering detection method for rotating machinery fault according to claim 1, characterized in that: S5 specifically includes the following steps: S51: Calculate the motor rotation frequency f based on the instantaneous frequency ridge line and the number of motor pole pairs r ; S52: Calculate multiple characteristic frequencies f when each component fails based on the motor frequency and the design parameters of the rotating machinery fault .

6. The current signal iterative Vold-Kalman filtering detection method for rotating machinery fault according to claim 1, characterized in that: S6 specifically includes the following parameters: S61: Construct a Vold-Kalman filter. The state equation of the Vold-Kalman filter is expressed as: Where D is a second-order n*n subtraction operator; A k is the envelope amplitude matrix of the signal; ε k is a non-homogeneous term; The observation equation of the Vold-Kalman filter is expressed as: X=A k C k +d k Where X is the n*1 observed current signal matrix; δ k is the n*1 generalized noise matrix; C k =diag[θ k(1) ,θ k(2) ,θ k(3) ,...,θ k(n) ] is the target component X k The carrier matrix of According to the state equation, observation equation and the pre-set carrier matrix C k , define the loss function as the non-homogeneous term ε k and generalized noise δ k The weighted energy and loss function are expressed as follows: Where α is the weighting factor, (·)H represents the conjugate transpose of the matrix; When the loss function is minimized, When the target component X is obtained k The amplitude envelope, amplitude envelope A k The expression is expressed as: Where I is the identity matrix; According to the carrier matrix C k and amplitude envelope A k The target component X to be extracted can be reconstructed k , the reconstruction method is as follows: X k =A k C k S62: Multiple characteristic frequencies f in S52 fault As the cutoff frequency of the Vold-Kalman filter, the current signal I in S2 is d Extract multiple target components (X1, X2, ..., X k ); S63: testing multiple target components to be fault harmonics or noises by a substitution test method, and determining the null hypothesis of the target component X1 - X1 is noise; S64: Generate m groups of substitute data of the target component X1 by iterative amplitude adjustment Fourier transform, and the power spectrum and amplitude distribution of the substitute data are consistent with the target component X1; S65: Calculate the envelope energy entropy Q[A] and instantaneous frequency energy entropy Q[F] of the target component and the m groups of surrogate data respectively, adaptively set different weights for Q[A] and Q[F], and construct a statistic of the surrogate test. The calculation formula of the test statistic is expressed as: D(ζ a ,g f )=ζ a Q[A]+ζ f Q[F] g a =1 / Q[A(X K )] ζ f =1 / Q[F(X K )] Among them, D(ζ a ,ζ f ) is the test statistic; f(x) is the target component or alternative data; Q[f(x)] is the energy entropy of f(x), which is used to measure the order of data. The larger the energy entropy, the greater the uncertainty of the data and the more chaotic the distribution; ζ a represents the weight of Q[A]; ζ f represents the weight of Q[F]; A(X K ) represents the Fourier coefficient of the target component envelope; F(X K ) represents the Fourier coefficient of the instantaneous frequency ridge of the target component; the weights are based on the target component X k Calculation, when the target component changes, the weight of the test statistic is adaptively adjusted; S66: When the test statistic of the target component is less than the test statistic of all alternative data, the target component is a fault harmonic, and the null hypothesis is rejected and the fault harmonic is saved; when the test statistic of the target component is greater than or equal to the test statistic of all alternative data, the target component is noise, and the null hypothesis is accepted and the target component is discarded; S67: Repeat S63 to S66 to iteratively extract all fault harmonics in the target component in S62.

7. The current signal iterative Vold-Kalman filtering detection method for rotating machinery fault according to claim 1, characterized in that: S7 specifically includes the following steps: S71: Perform time-frequency analysis on the fault harmonics in S67 by short-time Fourier transform; S72: Compare the time-frequency analysis result with the fault characteristic frequency of each component of the rotating machinery, so as to identify the fault type of the rotating machinery.

Images