A method for extracting characteristic frequency of rotating machinery under strong interference

Abstract

The application belongs to the field of big data learning models and provides a rotating machinery characteristic frequency extraction method under strong interference, which is based on the optimal balance parameter alpha in VMD decomposition roughly positioned according to the set reference modal correlation 最优 The balance parameter alpha of each decomposition is locally optimized by combining the sparrow search optimization algorithm. The modulus K of decomposition is adaptively determined by the iteration decomposition times, and the improved recursive VMD method avoids the influence of the inaccurate preset decomposition number and the balance parameter on the decomposition effect. The application is applied to the constructed simulation signal to achieve the optimal decomposition effect that can be achieved by the VMD method, and is applied to the pump cavitation flow-induced vibration signal processing to successfully achieve the effective extraction of the fluid machinery flow-induced vibration characteristic frequency under the condition of low signal-to-noise ratio.

Description

Technical Field

[0001] This invention belongs to the field of big data learning models, specifically relating to a method for extracting characteristic frequencies of rotating machinery under strong interference. Background Technology

[0002] Pumps, turbines, and propellers are typical examples of hydraulic rotating machinery, and flow-induced faults are inevitable during their operation. Abnormal flows such as tip vortex cavitation in propellers, tip clearance cavitation in axial-flow pumps, and vortex bands in turbine tailraces often induce accompanying phenomena like vibration. Vibration signals carry a wealth of information about flow-induced faults. The shaft frequency and blade frequency signals excited during the operation of hydraulic rotating machinery exhibit typical low-frequency characteristics. Under complex operating conditions, these characteristics are often contaminated by strong background noise and consist of superimposed multi-frequency features. Therefore, it is necessary to extract fault features through signal decomposition. Accurate demodulation of vibration signals under low signal-to-noise ratio conditions is a crucial step in fault diagnosis and target identification of hydraulic rotating machinery.

[0003] Various methods have been designed to address this problem, such as wavelet decomposition, wavelet packet decomposition, empirical mode decomposition (EMD), and local mean decomposition (LMD). However, wavelet decomposition and wavelet packet decomposition are non-adaptive signal analysis methods because the wavelet basis functions are pre-selected. Although EMD and LMD are adaptive signal processing methods, their application is limited due to mode mixing. Noise-aided techniques such as ensemble EMD and ensemble LMD alleviate the mode mixing problem to some extent, but their computational complexity increases dramatically, and they cannot effectively eliminate added white noise.

[0004] Variational Mode Decomposition (VMD) is a novel adaptive signal decomposition method proposed in recent years. With its rigorous mathematical guidance, fast convergence speed, strong noise robustness, and effective avoidance of endpoint effects, over-envelope, and under-envelope in decomposition, it has been widely applied in fault diagnosis. However, the effectiveness of VMD decomposition is highly dependent on the decomposition modulus K and the balancing parameter alpha. Currently, there is no unified method for determining these two parameters.

[0005] Although some studies have proposed improved VMD methods that combine optimization algorithms and fitness functions, most methods optimize alpha and K simultaneously, resulting in all submodes sharing the same alpha value. Since different submodes have different bandwidth characteristics, sharing alpha leads to unreasonable under-decomposition or over-decomposition. Summary of the Invention

[0006] To overcome the difficulties in setting the number of decomposition modes and the problem of all modes sharing the same equilibrium parameter in VMD applications, this invention draws on the idea of ​​recursive decomposition in EMD to adaptively determine the decomposition mode K, and combines the sparrow search optimization algorithm to locally optimize the equilibrium parameter a for each decomposition. Combining the two, this invention proposes an optimized fluid mechanical flow-induced vibration demodulation method for recursive VMD.

[0007] To achieve the above objectives, the present invention provides the following technical solution:

[0008] This invention provides a method for extracting characteristic frequencies of rotating machinery under strong interference, comprising:

[0009] S1: Set the acquired vibration signal as a residual signal and perform VMD operation on the residual signal;

[0010] S2: Set reference mode u ref And calculate the reference mode u ref cliff

[0011] S3: Set calibration test points for the balance parameter α, and after performing VMD decomposition on each calibration test point, calculate the signal sequence u of its decomposed mode u′ and the reference mode u′. ref The relevant parameters between them;

[0012] S4: By judging the correlation coefficient Reducing the number of calibration test points for the equilibrium parameter α further narrows down the optimal equilibrium parameter α. 最优 Within the specified range, reduce unnecessary VMD operations;

[0013] S5: Refine the optimization interval of the balance parameter α based on the characteristics of the decomposed signal;

[0014] S6: Determine the fitness function for optimization;

[0015] S7: Within the optimization range of the refined equilibrium parameter α, the Sparrow Search Optimization Algorithm (SSA) is used to select the optimal equilibrium parameter α for the target mode. 最优 ;

[0016] S8: Utilize the optimal equilibrium parameter α found in step S7 最优 The value of u is set and the decomposition modulus K=1 is used for VMD operation to extract the unique decomposition mode u. e ′;

[0017] S9: The unique decomposition mode u extracted in the residual signal removal step S8 e A new residual signal is then formed;

[0018] S10: Utilize the unique decomposition mode u extracted in step S8 e Perform signal reconstruction;

[0019] S11: Determine whether to continue the next mode decomposition iteration based on the power spectrum similarity between the reconstructed signal and the original signal.

[0020] Furthermore, in step (2), a reference mode u is set. ref The steps are as follows:

[0021] Initialize the equilibrium parameter α 初 The value of is set to c, the decomposition modulus K is set to 1, and the extracted unique mode is used as the reference mode u. ref ;

[0022] Reference mode u ref cliff The calculation formula is:

[0023]

[0024] Among them, u ref N, These are the signal sequence of the reference mode, the number of samples of the mode, and the average time-domain signal value of the reference mode, respectively.

[0025] Furthermore, the specific process of step S3 is as follows:

[0026] Initialize the equilibrium parameter α 初 When the value of α is set to c, the influence of the balancing parameter α on the signal decomposition result is greater when the value is within c. Therefore, a larger balancing parameter α interval is selected in the range of c to 10*c; a smaller balancing parameter α interval is selected within c, forming a calibration set α′. Each time a different value of the balancing parameter α is selected for VMD decomposition, the decomposition modulus K is always set to 1, and the decomposition mode is denoted as u′. After each decomposition, the relevant parameters of the decomposition mode u′ are calculated, including:

[0027] The correlation coefficient between the decomposed mode u′ and the reference mode

[0028]

[0029] Kur of decomposition mode u′ u′ :

[0030]

[0031] kurtosis of the envelope spectrum of decomposed mode u′ u′ :

[0032]

[0033] The ratio of kurtosis to envelope entropy of the decomposed mode u′ (KSE)u′ :

[0034]

[0035] in,

[0036] In is the logarithm to the base e;

[0037] It is the reference mode u ref The average value of the signal sequence;

[0038] It is the average value of the decomposed mode u′;

[0039] ES u′ (ω) is the envelope spectrum of u′ at frequency ω.

[0040] E u′ (j) is the envelope signal sequence E obtained after Hilbert demodulation of the decomposed mode u′. u′ discrete points;

[0041] P j It is E u′ The normalized form of (j);

[0042] J represents the envelope signal sequence E. u′ Sampling points;

[0043] j is the envelope signal sequence E u′ The sampling point located at the value of j.

[0044] Furthermore, in step S4, the optimal equilibrium parameter α is reduced. 最优 The process within the interval is as follows:

[0045] After performing VMD operation on the calibration test point of a certain equilibrium parameter α, if the correlation coefficient... If the value is greater than 0.95, the calibration test point of the balance parameter α is retained; otherwise, it is considered that the decomposed signal has a completely different center frequency from the reference mode. The calibration test point of the balance parameter α serves as one of the boundaries of the optimization interval of the refined balance parameter α. Subsequently, VMD operations will no longer be performed on other calibration test points that cross the boundary. The remaining calibration test points within the boundary form a new optimization interval.

[0046] The calibration test points for eliminating redundant equilibrium parameters α start from the calibration test points close to the equilibrium parameter α = c, and continue until the conditions in S3 are met. Stop. This ensures that the optimization interval of the selected equilibrium parameter α includes both the globally optimal equilibrium parameter α and... 最优 This can reduce unnecessary VMD operations and improve the computational performance of the algorithm.

[0047] Furthermore, the rule for refining the optimization interval of the balance parameter α in step S5 is as follows:

[0048] If the reference mode u ref cliff The extracted decomposition mode u′ is considered to have significant impact characteristics. In this case, the envelope kurtosis KSE of the decomposition mode u′ corresponding to the calibration test points of the remaining equilibrium parameter α within the boundary is... u′ Find the maximum value among them;

[0049] If the reference mode u ref cliff Then, the envelope spectrum kurtosis ESK of the decomposition mode u′ corresponding to the remaining calibration test points of the equilibrium parameter α within the boundary is... u′ Find the maximum value; with the value of the balance parameter α corresponding to the maximum value as the center, move one position to the left and one position to the right in the calibration test point set α′. The values ​​of the two corresponding balance parameters α form the optimization interval of the final refined balance parameter α.

[0050] Furthermore, the specific process of step S6 is as follows:

[0051] if The extracted decomposition mode u′ is considered to have significant impact characteristics, and the maximum KSE is selected in this case. u′ As a fitness function for optimization;

[0052] if Select the largest ESK u′ As a fitness function for optimization;

[0053] The optimal fitness function ESK for:

[0054]

[0055] Among them, ESK(u k ) is the decomposed mode u obtained by optimization within the optimization space. k The corresponding envelope kurtosis; u k It is a decomposition mode obtained by optimization within the refined optimization space;

[0056] KSE(u k ) is the decomposed mode u obtained by optimization within the optimization space. k The ratio of the corresponding kurtosis to the envelope entropy.

[0057] Furthermore, the specific process of step S7 is as follows:

[0058] S71: Randomly initialize the sparrow population and define relevant parameters, and define the maximum number of iterations.

[0059]

[0060] Among them, d represents the dimension of the optimization problem variables, n is the number of sparrows; X represents the position of the sparrow population;

[0061] S72: Calculate the fitness of the initial population, sort it, and then select the current optimal value and the worst value.

[0062]

[0063] Among them, f is determined by the fitness function determined in step S6;

[0064] S73: Update the position of the discoverer, and the formula is as follows:

[0065]

[0066] Among them, m represents the m-th sparrow;

[0067] t represents the current iteration number, n = 1, 2, 3,..., d; iter max is the set maximum number of iterations; X m，n represents the position information of the m-th sparrow in the n-th dimension; σ ∈ (0, 1] is a random number; R2 and ST respectively represent the warning value and the safety value, where R2 ∈ (0, 1], ST ∈ (0.5, 1]; Q represents a random number obeying the normal distribution; L represents a 1×d matrix, and all elements in this matrix are 1;

[0068] When R2 < ST, it means that there is no predator around the foraging environment at this time, and the discoverer can search in a larger range; <00003​​​​​​​​​​​​​​​​​​​​-1 ;

[0073] When m>n / 2, it indicates that the mth entrant with a lower fitness value needs to fly to other places to forage in order to obtain more energy.

[0074] S75: Update the location of sparrows that are aware of danger, using the following formula:

[0075]

[0076] Among them, X best It represents the current global optimal position; β, as the step size control parameter, is a random number between 0 and 1 following a normal distribution; Q is a random number between -1 and 1, fun m This is the fitness value of the current individual sparrow. g and fun w These are the current best and worst fitness values ​​globally, respectively; ε is a minimum constant to avoid the denominator being zero.

[0077] when fun i >fun g This indicates that the sparrows are currently on the edge of their population and extremely vulnerable to predators; X best Indicates the safest location in the population; fun i =fun g When Q indicates that the sparrow in the middle of the population is aware of the danger and needs to move closer to other sparrows to minimize their risk of being preyed upon; Q represents the direction of the sparrow's movement and is also a step size control parameter.

[0078] S76: Obtain the current optimal value. If the current optimal value is better than the optimal value of the previous iteration, perform an update operation. Otherwise, do not perform an update operation and continue the iteration operation until the condition is met. Finally, obtain the global optimal value and the best fitness value.

[0079] Furthermore, the signal reconstruction function f′ is:

[0080]

[0081] Furthermore, the judgment criteria for step S11 are as follows:

[0082] Set the power spectrum of the original signal to PS. f The formula is:

[0083] PS f =|f| 2 ;

[0084] Set the reconstructed signal power spectrum to PS. f′ The formula is:

[0085] PS f′ =|f′| 2 ;

[0086] Where f is the time series of the original signal;

[0087] Power spectral similarity coefficient between the reconstructed signal and the original signal for:

[0088]

[0089] in, is the mean of the original signal power spectrum; I is the total number of iterations;

[0090] q is the number of iterations;

[0091] if If the value exceeds the threshold of 0.9, the iteration stops.

[0092] if If the value is less than or equal to the threshold of 0.9, the iteration continues, and we return to step S1 to continue.

[0093] The threshold of 0.9 is derived from experience, and the specific threshold is not limited. The threshold of 0.9 is the optimal state.

[0094] The present invention has the following beneficial effects:

[0095] (1) In order to overcome the difficulties in setting the number of decomposition modes and the problem of all modes sharing the same equilibrium parameters in VMD applications, this invention proposes a method for extracting characteristic frequencies of rotating machinery under strong interference by drawing on the idea of ​​EMD recursive decomposition.

[0096] (2) The present invention roughly locates the optimal equilibrium parameter α in VMD decomposition based on the correlation of the set reference modes. 最优 Within the given interval, the balance parameter α for each decomposition is locally optimized using the sparrow search optimization algorithm.

[0097] (3) In selecting the optimal balance parameter α, this invention... 最优 The optimization interval of the selected balance parameter α is determined by starting from the calibration test point and then performing one test on each side of the selected calibration test point. This ensures that the optimization interval of the selected balance parameter α can include the globally optimal balance parameter, while reducing unnecessary VMD operations and improving the computational performance of the algorithm.

[0098] (4) The present invention adaptively determines the decomposition modulus K by iterative decomposition number, and the improved recursive VMD method adopted avoids the influence of inaccurate pre-setting of decomposition number and balance parameter on decomposition effect.

[0099] (5) The present invention was applied to the constructed simulation signal to achieve the best decomposition effect that the VMD method can achieve. At the same time, it was applied to the pump cavitation flow-induced vibration signal processing and successfully achieved the effective extraction of the characteristic frequency of fluid mechanical flow-induced vibration under low signal-to-noise ratio conditions. Attached Figure Description

[0100] Figure 1 This is a flowchart illustrating the extraction method of Example 1.

[0101] Figure 2 The spectrum diagram of the simulated harmonic signal x1(t) constructed in Example 1.

[0102] Figure 3 The spectrum diagram of the simulated variable frequency multi-component modulation signal x2(t) constructed in Example 1.

[0103] Figure 4 The spectrum diagram of the simulated random pulse c1(t) signal constructed in Example 1.

[0104] Figure 5 The spectrum diagram of the simulated periodic transient pulse c2(t) caused by the fault constructed in Example 1.

[0105] Figure 6 The spectrum of the simulated random pulse c1(t) signal constructed after adding Gaussian white noise in Example 2 is shown.

[0106] Figure 7 The spectrum of the simulated periodic transient pulse c2(t) caused by the fault constructed after adding Gaussian white noise in Example 2 is shown.

[0107] Figure 8 The simulated spectrum of the harmonic signal x1(t) constructed after adding Gaussian white noise in Example 2 is shown.

[0108] Figure 9 The above is a simulated spectrum of the variable-frequency multi-component modulated signal x2(t) constructed by adding Gaussian white noise in Example 2.

[0109] Figure 10 The spectrum diagram of the simulated signal constructed after adding Gaussian white noise in Example 2 is shown.

[0110] Figure 11 The spectrum diagram of the simulated signal constructed after adding Gaussian white noise in Example 2 is shown.

[0111] Figure 12 To demodulate the noisy new signal in Example 2 using the spectral kurtosis method.

[0112] Figure 13To apply the traditional VMD method with K=4 and initial equilibrium parameter α 初 The first mode in the decomposition result of =2000.

[0113] Figure 14 To apply the traditional VMD method with K=4 and initial equilibrium parameter α 初 The second mode in the decomposition result of =2000.

[0114] Figure 15 To apply the traditional VMD method with K=4 and initial equilibrium parameter α 初 The third mode in the decomposition result of =2000.

[0115] Figure 16 To apply the traditional VMD method with K=4 and initial equilibrium parameter α 初 The fourth mode in the decomposition result of =2000.

[0116] Figure 17 This is a diagram of vibration signals collected during pump cavitation.

[0117] Figure 18 This is a graph of vibration signals collected when the pump is not cavitating.

[0118] Figure 19 To apply the method proposed in this invention, the frequency domain and corresponding envelope spectrum of the vibration signal corresponding to the sub-mode were obtained under the pump air-conditioning condition in Example 3.

[0119] Figure 20 To apply the method proposed in this invention, the frequency domain and corresponding envelope spectrum of the corresponding sub-modes were obtained by performing VMD iteration 1 decomposition on the vibration signal collected under the pump air-conditioning condition in Example 3.

[0120] Figure 21 To apply the method proposed in this invention, the frequency domain and corresponding envelope spectrum of the corresponding sub-modes were obtained by performing VMD iteration 2 times on the vibration signal collected under the pump air-conditioning condition in Example 3.

[0121] Figure 22 To apply the method proposed in this invention, the frequency domain and corresponding envelope spectrum of the corresponding sub-mode were obtained by performing VMD iteration 3 times on the vibration signal collected under the pump air-conditioning condition in Example 3.

[0122] Figure 23 To apply the method proposed in this invention, the frequency domain and corresponding envelope spectrum of the corresponding sub-mode were obtained by performing VMD iteration 4 times on the vibration signal collected under the pump air-conditioning condition in Example 3.

[0123] Figure 24To apply the method proposed in this invention, the frequency domain and corresponding envelope spectrum of the corresponding sub-mode were obtained by performing VMD iteration 5 times on the vibration signal collected under the pump air-conditioning condition in Example 3.

[0124] Figure 25 To apply the method proposed in this invention, the frequency domain and corresponding envelope spectrum of the corresponding sub-modes were obtained by performing VMD iteration 6 times on the vibration signal collected under the pump air-conditioning condition in Example 3.

[0125] Figure 26 To apply the traditional VMD method to analyze the frequency domain and corresponding envelope spectrum of the vibration signal corresponding to the sub-modes collected under the pump-airing condition in Example 3.

[0126] Figure 27 To apply the traditional VMD method to decompose the vibration signal collected under the pump-airing condition in Example 3 by VMD iteration 1, the frequency domain of the corresponding sub-mode and the corresponding envelope spectrum are obtained.

[0127] Figure 28 To apply the traditional VMD method to perform two iterations of VMD decomposition on the vibration signal collected under the pump-airing condition in Example 3, the frequency domain of the corresponding sub-mode and the corresponding envelope spectrum are obtained.

[0128] Figure 29 To apply the traditional VMD method to decompose the vibration signal collected under the pump-airing condition in Example 3 by VMD iteration 3 times, the frequency domain of the corresponding sub-mode and the corresponding envelope spectrum are obtained.

[0129] Figure 30 To apply the traditional VMD method to decompose the vibration signal collected under the pump-airing condition in Example 3 by VMD iteration 4 times, the frequency domain of the corresponding sub-mode and the corresponding envelope spectrum are obtained.

[0130] Figure 31 To apply the traditional VMD method to decompose the vibration signal collected under the pump-airing condition in Example 3 by VMD iteration 5 times, the frequency domain of the corresponding sub-mode and the corresponding envelope spectrum are obtained.

[0131] Figure 32 To apply the traditional VMD method to decompose the vibration signal collected under the pump-airing condition in Example 3 by VMD iteration 6 times, the frequency domain of the corresponding sub-mode and the corresponding envelope spectrum are obtained. Detailed Implementation

[0132] The specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. It should be noted that the embodiments are only specific illustrations of the invention and should not be regarded as limitations on the invention. The purpose of the embodiments is to enable those skilled in the art to better understand and reproduce the technical solution of the present invention. The scope of protection of the present invention should still be determined by the scope defined in the claims.

[0133] To overcome the difficulties in setting the number of decomposition modes and the problem of all modes sharing the same equilibrium parameter in VMD applications, this invention draws on the idea of ​​EMD recursive decomposition to adaptively determine the decomposition mode K, and combines the sparrow search optimization algorithm to locally optimize the equilibrium parameter alpha for each decomposition. Combining the two, this invention proposes an optimized fluid mechanical flow-induced vibration demodulation method for recursive VMD.

[0134] Example 1

[0135] like Figure 1 As shown, the present invention provides a method for extracting characteristic frequencies of rotating machinery under strong interference, comprising:

[0136] S1: Set the acquired vibration signal as a residual signal and perform VMD operation on the residual signal;

[0137] like Figure 2-5 As shown, a multi-component amplitude-frequency modulation signal is constructed by inputting a vibration signal f(t), illustrating the advantages of this invention in anti-mode aliasing characteristics and noise robustness in signal decomposition. Meanwhile, in order to demonstrate the advantages of the improved VMD method in decomposing frequency-modulated multi-component modulation signals, other algorithms and the method of this invention are used to decompose the same signal respectively.

[0138] The formula for constructing a signal is as follows:

[0139]

[0140] Where x1(t) represents the harmonic signal, f1 represents the 25Hz shaft frequency; x2(t) represents the multi-component modulated signal with variable frequency; f2 is the signal frequency of 200Hz; c1(t) is used to simulate the random pulses caused by the bursting of bubbles near the pump body; f c1 The resonant frequency is 3000Hz; A is a random number between 0 and 2; c2(t) represents a periodic transient pulse that may be caused by a failure of the shaft and bearing system; the sampling frequency and the number of signals are set to 20kHz and 1s, respectively; the characteristic frequencies of the designed periodic harmonics and impulse signals are both 25Hz, while the characteristic frequency of the random pulse signal is 200Hz.

[0141] S2: Set reference mode u ref And calculate the reference mode u ref cliff

[0142] Initialize the equilibrium parameter α 初 The value of is set to 2000, the decomposition modulus K is set to 1, and the extracted unique mode is used as the reference mode u. ref This facilitates testing of each subsequent decomposed mode u′ and the reference mode u. refThe correlation;

[0143] Reference mode u ref cliff The calculation formula is:

[0144]

[0145] Among them, u ref N, These are the signal sequence of the reference mode, the number of samples of the mode, and the average time-domain signal value of the reference mode, respectively.

[0146] S3: Set the calibration test points for the equilibrium parameter α, and calculate the relevant parameters of the decomposed mode u′ after performing VMD decomposition on each calibration test point;

[0147] Initialize the equilibrium parameter α 初 When the value of α is set to 2000, the influence of the balance parameter α on the signal decomposition result is greater when the value is within 2000. Therefore, a larger interval of the balance parameter α is selected in the range of 2000-20000, specifically α = 6000, α = 10000, α = 15000, and α = 20000. Within 2000, a smaller interval of the balance parameter α is selected, specifically α = 1000, α = 600, and α = 200, and the decomposition results are tested respectively. The calibration set α′ = [200, 600, 1000, 2000, 6000, 10000, 15000, 20000]. Each time a different value of the balance parameter α is selected for VMD decomposition, the decomposition modulus K is always set to 1, and the decomposition mode is denoted as u′. After each decomposition, the relevant parameters of the decomposition mode u′ are calculated, including:

[0148] The correlation coefficient between the decomposed mode u′ and the reference mode

[0149]

[0150] Kur of decomposition mode u′ u′ :

[0151]

[0152] kurtosis of the envelope spectrum of decomposed mode u′ u′ :

[0153]

[0154] The ratio of kurtosis to envelope entropy of the decomposed mode u′ (KSE) u′ :

[0155]

[0156] in,

[0157] In is the logarithm to the base e;

[0158] It is the reference mode u ref The average value of the signal sequence;

[0159] It is the average value of the decomposed mode u′;

[0160] ES u′ (ω) is the envelope spectrum of u′ at frequency ω;

[0161] E u′ (j) is the envelope signal sequence E obtained after Hilbert demodulation of the decomposed mode u′. u′ discrete points;

[0162] P j It is E u′ The normalized form of (j);

[0163] J represents the envelope signal sequence E. u′ Sampling points;

[0164] j is the envelope signal sequence E u′ The sampling point located at the value of j.

[0165] S4: By judging the correlation coefficient Reducing the number of calibration test points for α further narrows down the optimal equilibrium parameter α. 最优 Within the specified range, reduce unnecessary VMD operations;

[0166] After performing VMD operation on the calibration test point of a certain equilibrium parameter α, if the correlation coefficient... If the value is greater than 0.95, the calibration test point of the balance parameter α is retained; otherwise, it is considered that the decomposed signal has a completely different center frequency from the reference mode. The calibration test point of the balance parameter α serves as one of the boundaries of the optimization interval of the refined balance parameter α. Subsequently, VMD operations will no longer be performed on other calibration test points that cross the boundary. The remaining calibration test points within the boundary form a new optimization interval.

[0167] The calibration test points for eliminating redundant balance parameters α start from the calibration test point close to the balance parameter α = 2000, that is, by performing one test on each side, gradually increasing the optimal balance parameter α from the balance parameter α = 2000. 最优 The search range continues until the conditions in S3 are met. Stop; the specific order is α = [2000, 1000, 6000, 600, 10000, 200, 15000, 20000]. This ensures that the optimization interval of the selected balance parameter α includes the globally optimal balance parameter α. 最优 This can reduce unnecessary VMD operations and improve the computational performance of the algorithm.

[0168] Through calculation, the remaining calibration test points for the equilibrium parameter α are α = [600, 1000, 2000, 6000, 10000, 15000, 20000], and α = 600 is eliminated;

[0169] S5: Refine the optimization interval of the balance parameter α based on the characteristics of the decomposed signal;

[0170] If the reference mode u ref cliff The extracted decomposition modes are considered to have significant impact characteristics. In this case, the envelope kurtosis KSE of the decomposition modes corresponding to the calibration test points of the remaining equilibrium parameter α within the boundary is... u′ Find the maximum value among them;

[0171] If the reference mode u ref cliff Then, the envelope spectrum kurtosis (ESK) of the decomposition modes corresponding to the remaining calibration test points of the equilibrium parameter α within the boundary is... u′ Find the maximum value; with the value of the balance parameter α corresponding to the maximum value as the center, move one position to the left and one position to the right in the calibration test point set α′. The values ​​of the two corresponding balance parameters α form the optimization interval of the final refined balance parameter α.

[0172] Calculated: The largest ESK u′ =1456. Since the balance parameter α = 6000, the optimization space of the refined balance parameter α is α = [2000, 10000].

[0173] S6: Determine the fitness function for optimization;

[0174] if The extracted decomposition modes are considered to have significant impact characteristics, and the maximum KSE is selected in this case. u′ As a fitness function for optimization;

[0175] if Select the largest ESK u′ As the fitness function for optimization; fitness function for optimization ESK for:

[0176]

[0177] Among them, ESK(u k ) is the decomposed mode u obtained by optimization within the optimization space. k The corresponding envelope kurtosis; u k It is a decomposition mode obtained by optimization within the refined optimization space;

[0178] KSE(u k ) is the decomposed mode u obtained by optimization within the optimization space. k The ratio of the corresponding kurtosis to the envelope entropy.

[0179] because Select the largest ESK u′ As a fitness function for optimization.

[0180] S7: Within the optimization range of the refined equilibrium parameter α, the Sparrow Search Optimization Algorithm (SSA) is used to select the optimal equilibrium parameter α for the target mode. 最优 ;

[0181] S71: Randomly initialize a sparrow population and define relevant parameters, including the maximum number of iterations;

[0182]

[0183] Where d represents the dimension of the optimization problem variables, n is the number of sparrows, and X represents the location of the sparrow population;

[0184] S72: Calculate the fitness of the initial population and sort it to select the current best and worst values.

[0185] Where f is determined by the fitness function determined in step S6;

[0186] S73: Update the discoverer's location, using the following formula:

[0187]

[0188] Where m represents the m-th sparrow;

[0189] t represents the current iteration number, n = 1, 2, 3, ..., d; iter max The maximum number of iterations is set; X m，n Let represent the position information of the m-th sparrow in the n-th dimension; σ∈(0,1] is a random number; R2 and ST represent the warning value and the safety value, respectively, where R2∈(0,1] and ST∈(0.5,1]; Q represents a random number that follows a normal distribution; L represents a 1×d matrix, where each element in the matrix is ​​1;

[0190] When R2 < ST, it indicates that there are no predators around the foraging environment at this time, and the discoverers can search within a larger range;

[0191] When R2 ≥ ST, it means that some sparrows in the population have discovered predators and sent out alarms to other sparrows in the population. At this time, all sparrows need to quickly fly to other safe places to forage;

[0192] S74: Update the position of the joiner, and the formula is as follows:

[0193]

[0194] where, X p is the optimal position currently occupied by the discoverer, X worst represents the current globally worst position; A represents a 1×d matrix, where each element is randomly assigned 1 or -1, and it is stipulated that A + = A T (AA T ) -1 ;

[0195] When m > n / 2, it indicates that the mth joiner with a lower fitness value needs to fly to other places to forage in time to obtain more energy.

[0196] S75: Update the position of the sparrows aware of danger, and the formula is as follows:

[0197]

[0198] where, X best is the current globally optimal position; β is the step size control parameter, which is a random number between 0 and 1 that follows a normal distribution; Q is a random number between -1 and 1, fun m is the fitness value of the current sparrow individual, fun g and fun w are the current globally best and worst fitness values respectively; ε is the smallest constant to avoid the denominator being zero;

[0199] When fun i > fun g , it indicates that the sparrow is at the edge of the population at this time and is extremely vulnerable to predators; X best represents the safest position in the population; fun i = fun g when, this indicates that the sparrows in the middle of the population are aware of danger and need to get closer to other sparrows to minimize their risk of being preyed upon; Q represents the direction of the sparrow's movement and is also the step size control parameter;

[0200] S76: Obtain the current optimal value. If the current optimal value is better than the optimal value of the previous iteration, perform an update operation; otherwise, do not perform an update operation and continue the iteration operation until the condition is met; finally, obtain the global optimal value and the best fitness value.

[0201] Within the optimization range of the refined equilibrium parameter α, an optimization procedure is used to select the optimal equilibrium parameter α for the target mode. 最优 The optimal equilibrium parameter α is obtained. 最优 =4192.

[0202] S8: Utilize the optimal equilibrium parameter α found in step S7 最优 =4192, and set the decomposition modulus K=1 to perform VMD operation, extracting the unique decomposition mode as u e ′;

[0203] S9: The unique decomposition mode u extracted in the residual signal removal step S8 e A new residual signal is then formed;

[0204] S10: Utilize the unique decomposition mode u extracted in step S8 e Signal reconstruction is performed, and the signal reconstruction function f′ is:

[0205]

[0206] S11: Determine whether to continue the next mode decomposition iteration based on the power spectrum similarity between the reconstructed signal and the original signal.

[0207] Set the power spectrum of the original signal to PS. f The formula is:

[0208] PS f =|f| 2 ;

[0209] Set the reconstructed signal power spectrum to PS. f′ The formula is:

[0210] PS f′ =|f′| 2 ;

[0211] Where f is the time series of the original signal;

[0212] Power spectral similarity coefficient between the reconstructed signal and the original signal for:

[0213]

[0214] in, It is the mean of the original signal power spectrum;

[0215] if If the value exceeds the threshold of 0.9, the iteration stops.

[0216] if If the value is less than the threshold of 0.9, the iteration continues, and the process returns to step S1. This is achieved through calculation... If the result is less than the threshold of 0.9, return to S1 for the second iteration decomposition.

[0217] The optimal equilibrium parameter α is decomposed in the second iteration. 最优 =200; after iterative decomposition If the value exceeds the threshold of 0.9, the iterative decomposition stops.

[0218] Example 2

[0219] Based on the signal x1(t)+x2(t)+c1(t)+c2(t) constructed in Example 1, different levels of Gaussian white noise are added, such as... Figure 6-12 As shown, the formula for constructing a new noisy signal is as follows:

[0220]

[0221] In the formula, n(t) represents the added noise signal, and η represents Gaussian white noise.

[0222] The above-mentioned noisy new signal x(t) is iteratively decomposed according to the method provided in this invention.

[0223] To quantify the decomposition accuracy of five methods—the improved recursive VMD method, the traditional VMD method, the CEEMD method, the VMD method that optimizes by finding the maximum correlation coefficient between the decomposed mode u′ and the input mode, and the SK method—under different signal-to-noise ratios, the correlation coefficient between the original signal and the reconstructed signal was used for characterization. Table 1 summarizes the comparison results.

[0224] Table 1. Comparison of results between the method of the present invention and the conventional method.

[0225]

[0226]

[0227] like Figure 13-16As shown, except for the SK method which failed, the other methods are very close to 1 under noise-free conditions; as the signal-to-noise ratio decreases, the decomposition ability of all methods weakens when processing analog signals. The improved recursive VMD method achieves the best decomposition results compared to the traditional VMD method, and in many cases, the improved recursive VMD method yields the best results that can be obtained by applying VMD. Although in reality, the sub-decomposition modes of the signal are unknown before decomposition, it is impossible to determine the optimal decomposition parameters by finding the maximum correlation coefficient between the decomposed sub-modes and the sub-signals of the original signal. The improved recursive VMD method proposed in this invention avoids this problem and can intelligently determine the parameters. Obviously, compared with the traditional VMD, CEEMD, and SK methods, the improved recursive VMD method still has satisfactory applicability in processing composite signals with different components, even in noisy environments. In summary, the above simulations fully highlight the anti-aliasing characteristics of the improved recursive VMD method and the advantages of the proposed method in noise robustness in processing multi-component non-stationary signals.

[0228] Where u1 is the first mode obtained by applying the decomposition method;

[0229] u2 is the second mode obtained by applying the decomposition method;

[0230] x1, x2, c1, and c2 are the simulated sub-signals of the harmonic signal x1(t), the multi-component modulated signal x2(t), the random pulse c1(t), and the periodic transient pulse c2(t), respectively.

[0231] Example 3

[0232] S1: Acquire pump cavitation signal, such as Figure 17-18 As shown; the above vibration signal is decomposed using the method of this invention; the signal is set as a residual signal, and a VMD operation is performed on the residual signal; the equilibrium parameter α is initialized. 初 Set to 2000, and set the decomposition modulus to 1;

[0233] S2: The extracted unique mode is used as the reference mode to facilitate testing the correlation between each subsequent decomposition and it; at the same time, the reference mode u is calculated. ref cliff

[0234] S3: Set up calibration test points for the balance parameter α and perform VMD operations according to its parameters. The balance parameter α has a greater impact on the signal decomposition results when its value is below 2000; therefore, select larger intervals of the balance parameter α within the range of 2000-20000, specifically α = 6000, α = 10000, α = 15000, and α = 20000. Within 2000, select α = 1000, α = 600, and α = 200 to test the decomposition results respectively. Form a calibration set α′ = [200, 600, 1000, 2000, 6000, 10000, 15000, 20000]. Each time a different value of the balance parameter α is selected for VMD decomposition, the decomposition modulus K is always set to 1; the decomposition mode is denoted as u′; the test starts from the calibration test point closest to 2000; after each decomposition, calculate the similarity between the decomposition mode u′ and the reference mode. Kur of decomposition mode u′ u ′、Envelope spectrum kurtosis ESK of decomposed mode u′ u′ The ratio of kurtosis to envelope entropy of the decomposed mode u′, KSE u′ Four parameters.

[0235] S4: Reduce the number of calibration test points for the equilibrium parameter α to minimize unnecessary VMD operations. When performing VMD operations on calibration test points for a certain equilibrium parameter α, it is found that if the correlation coefficient... If the value is greater than 0.95, the calibration test point of the balance parameter α will be retained; otherwise, it will be considered that the decomposed signal has a completely different center frequency from the decomposed signal under the reference standard. The calibration test point of the balance parameter α will be used as one of the boundaries of the optimization interval of the refined balance parameter α. VMD operation will no longer be performed on other calibration test points that cross the boundary.

[0236] S5: Refine the optimization interval of the balance parameter α based on the characteristics of the decomposed signal;

[0237] if The extracted decomposed modes are considered to have significant impact characteristics. In this case, the KSE corresponding to the remaining few α calibration test points within the boundary... u′ Find the maximum value among them;

[0238] if Then, the ESK corresponding to the remaining calibration test points of the equilibrium parameter α within the boundary. u′ Find the maximum value. Centering on the value of the balance parameter α corresponding to the maximum value, shift one position to the left and one position to the right in the calibration set α′. The values ​​of the two corresponding balance parameters α form the optimization interval of the final refined balance parameter α.

[0239] S6: Determine the fitness function for optimization;

[0240] if The extracted decomposition modes are considered to have significant impact characteristics, and the maximum KSE is selected in this case. u′ As a fitness function for optimization.

[0241] if Select the largest ESK u′ As a fitness function for optimization.

[0242] S7: Utilize the optimal equilibrium parameter α found in S6 最优的 The value and K=1 are used for VMD operation; the unique decomposition mode is extracted and denoted as u. i .

[0243] S8: The residual signal is formed by removing the decomposition mode extracted in S7;

[0244] S9: Use the decomposed modes extracted in S7 to reconstruct the signal until the correlation coefficient between the reconstructed signal and the original signal is greater than 0.95.

[0245] Figure 26-32 To apply the traditional VMD method to the vibration signal collected under the pump-airing condition in Example 3, the frequency domain and corresponding envelope spectrum of the sub-modes were obtained after VMD iterative decomposition. The analysis of the sub-mode envelope spectrum shows that the characteristic frequencies shaft frequency (BPF), blade frequency (RF), as well as BPF-RF and BPF+RF, were effectively extracted.

[0246] Figure 19-25 To apply the method proposed in this invention, the frequency domain and corresponding envelope spectra of the sub-modes corresponding to the vibration signals collected under pump-air conditions in Example 3 were obtained after VMD iterative decomposition. Analysis of the sub-mode envelope spectra showed that the characteristic frequencies shaft frequency (BPF), blade frequency (RF), as well as BPF-RF and BPF+RF, were effectively extracted. Through comparison... Figure 19-25 and Figure 26-32 Among modes with the same frequency, the modes decomposed by the method proposed in this invention have higher kurtosis, which shows that the improved recursive VMD method can effectively extract the obvious shaft frequency and blade frequency caused by fluid mechanical flow-induced vibration, proving the superiority of the method proposed in this invention.

[0247] Although preferred embodiments of this application have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments as well as all changes and modifications falling within the scope of this application.

Claims

1. A method for extracting characteristic frequencies of rotating machinery under strong interference, characterized in that, Including: S1: Set the collected vibration signal as the residual signal, and perform VMD operation on the residual signal; S2: Set reference mode u ref And calculate the reference mode u ref cliff S3: Set calibration test points for the balance parameter α, and after performing VMD decomposition on each calibration test point, calculate the signal sequence u of its decomposed mode u′ and the reference mode u′. ref The relevant parameters between them; S4: By judging the correlation coefficient Reducing the number of calibration test points for the equilibrium parameter α further narrows down the optimal equilibrium parameter α. 最优 Within the specified range, reduce unnecessary VMD operations; S5: Refine the optimization interval of the balance parameter α according to the characteristics of the decomposed signal; S6: Determine the fitness function for optimization; S7: Within the optimization range of the refined equilibrium parameter α, a sparrow search optimization algorithm is used to select the optimal equilibrium parameter α for the target mode. 最优 ; S8: Utilize the optimal equilibrium parameter α found in step S7 最优 The value of u is set and the decomposition modulus K=1 is used for VMD operation to extract the unique decomposition mode u. e ′; S9: The unique decomposition mode u extracted in the residual signal removal step S8 e A new residual signal is then formed; S10: Utilize the unique decomposition mode u extracted in step S8 e Perform signal reconstruction; S11: Judge whether to continue the next iteration of modal decomposition based on the power spectrum similarity between the reconstructed signal and the original signal.

2. The method for extracting characteristic frequencies of rotating machinery under strong interference according to claim 1, characterized in that, In step (2), the reference mode u is set. ref The steps are as follows: Initialize the equilibrium parameter α 初 The value of is set to c, the decomposition modulus K is set to 1, and the extracted unique mode is used as the reference mode u. ref ; Reference mode u ref cliff The calculation formula is: Among them, u ref N, These are the signal sequence of the reference mode, the number of samples of the mode, and the average time-domain signal value of the reference mode, respectively.

3. The method for extracting characteristic frequencies of rotating machinery under strong interference according to claim 1, characterized in that, The specific process of step S3 is as follows: Initialize the equilibrium parameter α 初 When the value of α is set to c, the influence of the balance parameter α on the signal decomposition result is greater when the value of α is within c. Therefore, a larger balance parameter α interval is selected in the range of c to 10*c; and a smaller balance parameter α interval is selected when the value is within c. When different values of the balance parameter α are selected each time to form the calibration set α′ for VMD decomposition, the decomposition modulus K is always set to 1, the decomposition mode is denoted as u′, and the relevant parameters of its decomposition mode u′ are calculated after each decomposition, including: The correlation coefficient between the decomposed mode u′ and the reference mode Kur of decomposition mode u′ u′ : kurtosis of the envelope spectrum of decomposed mode u′ u′ : The ratio of kurtosis to envelope entropy of the decomposed mode u′ (KSE) u′ : in, In is the natural logarithm with base e; It is the reference mode u ref The average value of the signal sequence; It is the average value of the decomposed mode u′; ES u′ (ω) is the envelope spectrum of u′ at frequency ω. E u′ (j) is the envelope signal sequence E obtained after Hilbert demodulation of the decomposed mode u′. u′ discrete points; P j It is E u′ The normalized form of (j); J represents the envelope signal sequence E. u′ Sampling points; j is the envelope signal sequence E u′ The sampling point located at the value of j.

4. The method for extracting characteristic frequencies of rotating machinery under strong interference according to claim 3, characterized in that, In step S4, the optimal equilibrium parameter α is reduced. 最优 The process within the interval is as follows: After performing VMD operation on the calibration test point of a certain equilibrium parameter α, if the correlation coefficient... If the value is greater than 0.95, the calibration test point of the balance parameter α is retained; otherwise, it is considered that the decomposed signal has a completely different center frequency from the reference mode. The calibration test point of the balance parameter α serves as one of the boundaries of the optimization interval of the refined balance parameter α. Subsequently, VMD operations will no longer be performed on other calibration test points that cross the boundary. The remaining calibration test points within the boundary form a new optimization interval. The calibration test points for eliminating redundant equilibrium parameter α start from the calibration test points with a value close to c for equilibrium parameter α, and continue until the conditions in S3 are met. stop.

5. The method for extracting characteristic frequencies of rotating machinery under strong interference according to claim 1, characterized in that, The rule for refining the optimization interval of the balance parameter α in step S5 is: If the reference mode u ref cliff The extracted decomposition mode u′ is considered to have significant impact characteristics. In this case, the envelope kurtosis KSE of the decomposition mode u′ corresponding to the calibration test points of the remaining equilibrium parameter α within the boundary is... u′ Find the maximum value among them; If the reference mode u ref cliff Then, the envelope spectrum kurtosis ESK of the decomposition mode u′ corresponding to the remaining calibration test points of the equilibrium parameter α within the boundary is... u′ Find the maximum value; with the value of the balance parameter α corresponding to the maximum value as the center, move one position to the left and one position to the right in the calibration test point set α′. The values ​​of the two corresponding balance parameters α form the optimization interval of the final refined balance parameter α.

6. The method for extracting characteristic frequencies of rotating machinery under strong interference according to claim 1, characterized in that, The specific process of step S6 is as follows: if The extracted decomposition mode u′ is considered to have significant impact characteristics, and the maximum KSE is selected in this case. u′ As a fitness function for optimization; if Select the largest ESK u′ As the fitness function for optimization; fitness function for optimization ESK for: Among them, ESK(u k ) is the decomposed mode u obtained by optimization within the optimization space. k The corresponding envelope kurtosis; u k It is a decomposition mode obtained by optimization within the refined optimization space; ESK(u k ) is the decomposed mode u obtained by optimization within the optimization space. k The ratio of the corresponding kurtosis to the envelope entropy.

7. The method for extracting characteristic frequencies of rotating machinery under strong interference according to claim 1, characterized in that, The specific process of step S7 is as follows: S71: Randomly initialize the sparrow population and define relevant parameters, and define the maximum number of iterations; Among them, d represents the dimension of the optimization problem variable, s is the number of sparrows; X represents the position of the sparrow population; S72: Calculate the fitness of the initial population and sort it to select the current optimal value and the worst value; Among them, fun is determined by the fitness function in step S6; S73: Update the position of the discoverer, and the formula is as follows: Among them, m represents the m-th sparrow; t represents the current iteration number, n = 1, 2, 3, ..., d; iter max The maximum number of iterations is set; X m，n Let represent the position information of the m-th sparrow in the n-th dimension; σ∈(0,1] is a random number; R2 and ST represent the warning value and the safety value, respectively, where R2∈(0,1] and ST∈(0.5,1]; Q represents a random number that follows a normal distribution; L represents a 1×d matrix, where each element in the matrix is ​​1; When R2 < ST, it means that there are no predators around the foraging environment at this time, and the discoverer can search in a larger range; When R2 ≥ ST, it means that some sparrows in the population have discovered predators and issued an alarm to other sparrows in the population. At this time, all sparrows need to quickly fly to other safe places to forage; S74: Update the position of the joiner, and the formula is as follows: Among them, X p This is currently the optimal position occupied by the discoverer, X. worst Let A represent the current worst position globally; let A be a 1×d matrix where each element is randomly assigned a value of 1 or -1, and let A be a 1×d matrix. + =A T (AA T ) -1 ; When m > n / 2, it indicates that the m-th joiner with a lower fitness value needs to fly to other places to forage in time to obtain more energy; S75: Update the position of the sparrows aware of danger, and the formula is as follows: Among them, X best It represents the current global optimal position; β, as the step size control parameter, is a random number between 0 and 1 following a normal distribution; Q is a random number between -1 and 1, fun m This is the fitness value of the current individual sparrow. g and fun w These are the current best and worst fitness values ​​globally, respectively; ε is a minimum constant to avoid the denominator being zero. when fun i >fun g This indicates that the sparrows are currently on the edge of their population and extremely vulnerable to predators; X best Indicates the safest location in the population; fun i =fun g This indicates that the sparrows in the middle of the population are aware of the danger and need to move closer to other sparrows to minimize their risk of being preyed upon; Q represents the direction of the sparrow's movement and is also a step size control parameter. S76: Obtain the current optimal value. If the current optimal value is better than the optimal value of the previous iteration, update it, otherwise do not update it, and continue the iterative operation until the condition is met; finally, obtain the global optimal value and the best fitness value.

8. The method for extracting characteristic frequencies of rotating machinery under strong interference according to claim 1, characterized in that, The signal reconstruction function f′ is:

9. The method for extracting characteristic frequencies of rotating machinery under strong interference according to claim 1, characterized in that, The judgment basis of step S11 is: Set the power spectrum of the original signal to PS. f The formula is: PS f =|f| 2 ; Set the reconstructed signal power spectrum to PS. f′ The formula is: PS f′ =|f′| 2 ; Among them, f is the original signal sequence; Power spectral similarity coefficient between the reconstructed signal and the original signal for: in, It is the mean of the original signal power spectrum; I is the total number of iterations; q is the number of iterations; if If the value exceeds the threshold of 0.9, the iteration stops. if If the value is less than or equal to the threshold of 0.9, the iteration continues, and the process returns to step S1.

Images