A method for determining parameters of variational mode decomposition

Abstract

The present invention discloses a method for determining parameters of variational mode decomposition, which comprises the following steps: first, determining the search range of K and α, determining the total number of iterations N T , the current number of iterations N C Set to 0; Step 2, select P different variational mode decomposition parameters (K, α) within the search range; Step 3, for all P variational mode decomposition parameters, calculate the objective function of each variational mode decomposition parameter; Step 4, retain the variational mode decomposition parameters that obtain the smallest objective function result, and use the optimization algorithm for each variational mode decomposition parameter that obtains the largest objective function result to obtain the updated parameters; Step 5, N C Add 1; if N T Greater than N C , go to step 3; if N T Equal to N C , the automatic optimization of parameters is completed, and the parameters that obtain the minimum objective function result are the optimal parameters. This method automatically optimizes K and α of variational mode decomposition, solving the problem of incorrect decomposition results caused by manual setting of two parameters.

Description

Technical Field

[0001] The invention belongs to the technical field of signal processing and relates to a method for determining parameters of variational mode decomposition. Background Art

[0002] Modal decomposition is a method for decomposing and detecting complex signals. Currently, decomposition algorithms represented by variational modal decomposition have become an important method for signal extraction.

[0003] Variational mode decomposition (Dragomiretskiy K, Zosso D. "Variational mode decomposition." IEEE Transactions on Signal Processing, 2014; Vol. 62, No. 3, pp. 531-544) can adaptively decompose complex time series signals into simple modal components, and achieve effective separation of intrinsic mode functions and frequency domain division of signals by determining the optimal center frequency and finite bandwidth of each mode. Variational mode decomposition has good noise resistance and can overcome the frequency aliasing problem when using empirical mode decomposition. Variational mode decomposition has received widespread attention since it was proposed, and has been widely used in financial data analysis, seismic data analysis, medical disease detection, mechanical fault diagnosis and other fields.

[0004] Variational mode decomposition is an adaptive signal processing method that has good processing effects on non-stationary and nonlinear signals. It can use iteration to find the optimal value in the variational mode. It can divide the sequence data in the time dimension into a set of intrinsic mode functions (IMFs) with limited width. Each intrinsic mode function has a limited bandwidth and a relatively concentrated frequency. Therefore, it can adaptively change the optimal center frequency and bandwidth between different intrinsic mode functions.

[0005] In variational mode decomposition, it is assumed that the original signal s is composed of many eigenmode functions, and each eigenmode function is regarded as an amplitude-frequency modulated signal:

[0006] u k (t) = A k (t)cos(φ k (t)) (1);

[0007] Where t represents time; let the total number of intrinsic mode function signals be K, and the kth intrinsic mode function signal be u k (t), k takes the value of 1, 2, ..., K; A k (t) is the instantaneous amplitude of the kth intrinsic mode function signal changing with time t; φ k (t) is the instantaneous phase of the kth intrinsic mode function signal changing with time t. For the intrinsic mode function signal uk (t), satisfying: φ k (t) is a non-decreasing function whose derivative is non-negative; A k (t) is non-negative.

[0008] Variational mode decomposition is a process of solving variational problems. Variational problems can be viewed as obtaining K intrinsic mode function signals u k (t), in order to make each eigenmode function signal have a modal component of finite bandwidth at the center frequency, and at the same time minimize the sum of the estimated bandwidths of each eigenmode function signal, the constraint condition is that the sum of all eigenmode function signals is equal to the original signal, then the corresponding process includes three steps:

[0009] (1) Using Hilbert transform, find the analytical signal of each intrinsic mode function signal;

[0010] (2) Mix each eigenmode function signal with the complex exponential of its center frequency estimate, using the exponential Modify the signal to shift its spectrum to the appropriate baseband;

[0011] (3) Estimate the bandwidth of each intrinsic mode function signal: that is, the square of the L2 norm of the gradient.

[0012] The above process is expressed in mathematical formula as follows:

[0013]

[0014] Where s is the original signal, u k is the kth intrinsic mode function signal, the intrinsic mode function signal u k The center frequency is ω k , δ(t) is the impulse function, j is the imaginary part of the complex number, t is time, π is the circumference of the circle, e is the natural exponential, * indicates convolution, ∑ indicates summation, min indicates minimum value, and ||||2 indicates L2 norm. It means to find the square of L2 norm. It means partial derivative with respect to t.

[0015] In order to find the optimal solution to the variational constrained problem, the Lagrange multiplier λ(t) and the penalty factor α are introduced to transform the variational constrained problem into a variational unconstrained problem. The Lagrange multiplier λ(t) can ensure the strictness of the constraints, and the penalty factor α can make the signal reconstruction more accurate under Gaussian noise conditions. The variational unconstrained problem is expressed as follows:

[0016]

[0017] Among them, λ(t) is the Lagrange multiplier, α is the penalty factor, s(t) is the original signal, It means to find λ(t) and The inner product of .

[0018] To solve this constraint problem, we can use the alternating direction multiplier method by updating and λ n+1 Solve the saddle point of the Lagrangian function in the above problem, use the Parseval / Plancherel Fourier isometry transform to transform the problem into the frequency domain, and get the update method as follows:

[0019]

[0020]

[0021] in, and are s(t) and u respectively. k The Fourier transform of (t) and λ(t), n represents the number of iterations, Represents the nth iteration Represents the n+1th iteration represents the ω of the n+1th iteration k For a given error ε, until The iteration ends when K intrinsic mode function signals are obtained, and ε takes the empirical value of 1×10 -9 .

[0022] In actual application, variational mode decomposition requires setting the total number K of intrinsic mode function signals and the penalty factor α. However, both have a great influence on the actual results of variational mode decomposition. When the value of K is too large, false modal components are generated; when the value of K is too small, it will cause modal aliasing. Improper α value will also produce incorrect decomposition results. Moreover, the characteristics of the signal to be decomposed affect the selection of parameters, and the use of empirical values ​​of the parameters cannot obtain correct results. In order to better apply variational mode decomposition to the signal, it is necessary to determine the appropriate parameters K and α.

[0023] Gravitational Search Algorithm (GSA) is a group search algorithm that is derived from Newton's law of universal gravitation and the law of motion. Its main ideas include the following aspects: First, in space, any two particles will be affected by universal gravitation, the direction of the force is the line connecting the centers of mass of the two particles, and the magnitude of the gravitational force conforms to the universal gravitation formula. Second, the gravitational force produces acceleration on the particles.

[0024] In the gravity search algorithm, a particle is abstracted into a point in the solution space, called an individual. There is a gravitational force between individuals, and the magnitude of this gravitational force depends on the mass of the two individuals and the distance between the two individuals. The mass of the individual is determined by the value of the objective function. In the solution space, the acceleration of each individual is calculated by dividing the sum of the attractive forces of other individuals by its mass. The smaller the value of the objective function, the greater the mass, and the greater the acceleration provided by the individual, so that the individual moves in the direction of a better solution. As the movement continues, individuals with large mass occupy an advantageous position, and eventually the entire group will gather towards the individual with the largest mass in the population, thereby finding the optimal solution to the problem.

[0025] The gravity search algorithm has the following advantages: First, the gravity search algorithm is a group intelligence optimization algorithm, which is characterized by emphasizing the interaction and influence between individuals. Second, in the gravity search algorithm, gravity is an important tool for information transmission between individuals, thereby realizing the sharing of optimization information between individuals. Third, the gravity search algorithm enables all individuals in the group to exchange information, and individuals will change their search behavior based on the information obtained, which makes the entire group have some abilities that independent individuals do not have, so it can make the group more effective than a single individual in finding the global optimal solution. Fourth, the gravity search algorithm requires fewer parameters to be set, the process is simple, and it is easy to implement and easy to apply in engineering. Fifth, it does not require the objective function to have mathematical properties such as differentiability or continuity, which makes it easy to apply to parameter optimization of variational mode decomposition. Summary of the invention

[0026] The purpose of the present invention is to provide a method for determining parameters of variational mode decomposition, which automatically optimizes the input parameters K and α of variational mode decomposition and solves the problem that manual setting of two parameters leads to incorrect decomposition results.

[0027] The objective of the present invention is achieved through the following technical solutions:

[0028] A method for determining parameters of variational mode decomposition comprises the following steps:

[0029] The first step is to determine the total number of intrinsic mode function signals K and the search range of the penalty factor α, where K is a positive integer; determine the total number of iterations N T ; Current number of iterations N C Set to 0;

[0030] Step 2: Select P different variational mode decomposition parameters (K, α) within the search range;

[0031] The third step is to calculate the objective function of each variational mode decomposition parameter for all P variational mode decomposition parameters. The objective function L of the nth variational mode decomposition parameter is n The calculation formula is as follows:

[0032] L n =α+10 4 ×n max (6);

[0033] Where n=1,2,…,P; n max The maximum number of maximum points in the amplitude spectrum of each intrinsic mode function signal after the variational mode decomposition is completed using the parameters;

[0034] Step 4: retain the parameters of the variational mode decomposition that obtains the small objective function result, and use the optimization algorithm for each variational mode decomposition parameter that obtains the large objective function result to obtain the updated parameters;

[0035] Step 5: Current number of iterations N C Add 1; if the total number of iterations is N T Greater than the current number of iterations N C , go to step 3; if the total number of iterations N T Equal to the current number of iterations N C , then the automatic optimization of parameters ends, and the parameters that obtain the minimum objective function result are the optimal parameters.

[0036] Compared with the prior art, the present invention has the following advantages:

[0037] The present invention provides a method for automatically optimizing the parameters of variational mode decomposition, and continuously updates the positions of the parameters through an optimization algorithm, gradually reduces the loss function, and finally obtains the global optimal parameters; at the same time, in the optimization process, the present invention adopts two types of loss functions, namely, under-decomposition loss and energy loss, to limit the decomposition parameters from two aspects: whether the decomposition is insufficient and whether the lost energy is too large, so that the decomposition parameters can give full play to the performance of variational mode decomposition, and the result is optimal. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] Figure 1 It is the main process of the parameter determination method of variational mode decomposition of the present invention;

[0039] Figure 2 is the time domain waveform and amplitude spectrum of the original signal s(t);

[0040] Figure 3 For the case of under-decomposition;

[0041] Figure 4 It is the statistical maximum point of the amplitude spectrum after two morphological filtering;

[0042] Figure 5 for Figure 4 Local magnification of

[0043] Figure 6 is the result of using the variational mode decomposition of the present invention;

[0044] Figure 7 It is the result of variational mode decomposition in reference 1;

[0045] Figure 8 It is the result of variational mode decomposition in reference 2;

[0046] Fig. 9 This is the result of variational mode decomposition in reference 3. DETAILED DESCRIPTION

[0047] The technical solution of the present invention is further described below in conjunction with the accompanying drawings, but is not limited thereto. Any modification or equivalent replacement of the technical solution of the present invention without departing from the spirit and scope of the technical solution of the present invention should be included in the protection scope of the present invention.

[0048] The present invention provides a method for determining parameters of variational mode decomposition. The method designs a loss function to limit the decomposition parameters from two aspects: whether the decomposition is insufficient and whether the lost energy is too large, so that the decomposition parameters can give full play to the performance of variational mode decomposition. The position of the parameters is continuously updated through the optimization algorithm, and the loss function is gradually reduced, and finally the global optimal parameters are obtained. At the same time, the amplitude spectrum of each intrinsic mode function is subjected to two morphological filtering and threshold setting methods to count the extreme points, thereby eliminating the influence of noise, so as to ensure that accurate optimal parameters are finally obtained. Figure 1 As shown, the following steps are included:

[0049] The first step is to determine the total number of intrinsic mode function signals K and the search range of the penalty factor α, where K is a positive integer; determine the total number of iterations N T ; Current number of iterations N C Set to 0.

[0050] The second step is to select P different variational mode decomposition parameters (K, α) within the search range.

[0051] The third step is to calculate the objective function of each variational mode decomposition parameter for all P variational mode decomposition parameters. The objective function L of the nth variational mode decomposition parameter is n The calculation formula is as follows:

[0052] L n =α+10 4 ×n max (6);

[0053] Objective function Ln The smaller the result, the better the variational mode decomposition effect is and the better the parameter is; where n = 1, 2, …, P; n max It is the maximum value of the number of maximum points contained in the amplitude spectrum of each intrinsic mode function signal after completing the variational mode decomposition using parameters.

[0054] In this step, n max The calculation method of is as follows: for each intrinsic mode function signal amplitude spectrum after the variational mode decomposition is completed using the parameter, find the number of maximum points it contains, and the number of the largest maximum points is n max .

[0055] In this step, the signal amplitude spectrum of each intrinsic mode function is used to find the number of maximum points contained in it using a mathematical morphology method. The specific method is as follows:

[0056] Step (1) performing a morphological closing operation on the amplitude spectrum to obtain a main peak envelope, wherein: a structure element of the morphological closing operation is a circle, and the width of the structure element is one percent of the width of the amplitude spectrum;

[0057] Step (2) performing a morphological closing operation on the main peak envelope obtained in step (1) to obtain a final envelope of the amplitude spectrum, wherein: the structural element of the morphological closing operation is a straight line, and the width of the structural element is 2% of the width of the amplitude spectrum;

[0058] Step (3) counts the number of maximum points on the final envelope of the amplitude spectrum obtained in step (2), and the specific steps are as follows:

[0059] Step (31) sets a maximum amplitude threshold: the maximum amplitude threshold is one thirtieth of the maximum value of the amplitude spectrum;

[0060] Step (32) sets the maximum value interval threshold: the maximum value interval threshold is one percent of the amplitude spectrum width;

[0061] Step (33) For each frequency position in the amplitude spectrum, the frequency position that satisfies the following conditions is a maximum point, and is counted in the number of maximum points:

[0062] If there is a value on the left side of the position, the value of the position is greater than the value on the left; if there is a value on the right side of the position, the value of the position is greater than the value on the right; the value of the position is greater than the maximum amplitude threshold; the distance between the position and other maximum points on the left and right sides is greater than the maximum interval threshold.

[0063] The fourth step is to retain the variational mode decomposition parameters that obtain a small objective function result, and use an optimization algorithm for each variational mode decomposition parameter that obtains a large objective function result to obtain updated parameters.

[0064] In this step, the method of using the optimization algorithm to obtain the updated parameters is as follows:

[0065] Suppose the parameter (K p ,α p ), p = 1, 2, ..., P; get the updated variational mode decomposition parameters (K ′ p ,α′ p ), using the gravity search optimization algorithm, the following calculation process is used:

[0066] Nth C At the iteration, the gravitational constant G is:

[0067]

[0068] Among them, e is a natural constant;

[0069] Among all P variational mode decomposition parameters, the quality of the nth variational mode decomposition parameter is denoted by M n , the calculation formula is as follows:

[0070]

[0071] Among them, the intermediate variable m corresponding to the parameter of the nth variational mode decomposition n The calculation formula is as follows:

[0072]

[0073] Among them, L b is the minimum value of the objective function obtained by all P variational mode decomposition parameters; L w is the maximum value of the objective function obtained by all P variational mode decomposition parameters, L n is the value of the objective function obtained by the parameters of the nth variational mode decomposition.

[0074] Similarly, the intermediate variable m corresponding to the parameter of the jth variational mode decomposition is j The calculation formula is as follows:

[0075]

[0076] Among them, L j is the value of the objective function obtained by the parameters of the j-th variational mode decomposition, j = 1, 2, …, P.

[0077] Assume that the parameters of the ith variational mode decomposition are (K i ,α i ), the parameter of the jth variational mode decomposition is (K j ,α j), is the parameter K of the ith variational mode decomposition i and the parameter K of the jth variational mode decomposition j The force between them is calculated as follows:

[0078]

[0079] Among them, R ij is the parameter of the ith variational mode decomposition (K i ,α i ) and the parameters of the jth variational mode decomposition (K j ,α j ), the calculation formula is as follows:

[0080]

[0081] is the parameter α of the ith variational mode decomposition i and the parameter α of the jth variational mode decomposition j The force between them is calculated as follows:

[0082]

[0083] For the pth variational mode decomposition parameter (K p ,α p ), parameter K p The total force F K for:

[0084]

[0085] Where R is a random number between 0 and 1;

[0086] Similarly, the parameter α p The total force F α for:

[0087]

[0088] Then the parameter K p The acceleration a K for:

[0089]

[0090] Parameter α p The acceleration a α for:

[0091]

[0092] Since the search range of K is small and the search range of α is large, the parameter K p and parameter α p The corresponding speed is calculated using different strategies, and the parameter K p The speed is:

[0093] v K =R+a K (17);

[0094] Parameter α p The speed is:

[0095] v α =R×α+a α (18);

[0096] The new variational mode decomposition parameters are as follows:

[0097]

[0098] α′ p =α p +v α (20);

[0099] in, Indicates K p +v K The calculation result is rounded to the nearest integer, so the final K′ p is an integer.

[0100] When the variational mode decomposition parameter (K′ p ,α′ p ) and also satisfy:

[0101] Condition 1: K′ p and α′ p Within the search range of parameters K and α respectively;

[0102] Condition 2: K′ p ≠K p , and α′ p ≠α p ; When the pth variational mode decomposition parameter (K p ,α p ) is updated to (K′ p ,α′ p );

[0103] Otherwise, repeat the above process of using the gravity search optimization algorithm to obtain the variational mode decomposition parameters until (K′ p ,α′ p ) satisfies the above two conditions.

[0104] Step 5: Current number of iterations NC Add 1; if the total number of iterations is N T Greater than the current number of iterations N C , go to step 3; if the total number of iterations N T Equal to the current number of iterations N C , then the automatic optimization of parameters ends, and the parameters that obtain the minimum objective function result are the optimal parameters.

[0105] Example:

[0106] In this embodiment, Figure 2 As shown in FIG. 1 , the original vibration signal s(t) of the hydropower equipment collected for condition monitoring contains 4650 sampling points, lasts for 4.65 seconds, and has a sampling rate of 1000 Hz. 101 is a time domain diagram, and 102 is a frequency domain diagram.

[0107] After using variational mode decomposition to obtain each eigenmode function, the number of peaks in the amplitude spectrum of each function can reflect the degree of under-decomposition. If the parameters are set improperly, an eigenmode function contains part of other eigenmode functions, and the amplitude spectrum of this eigenmode function will contain multiple peaks. For example, Figure 3 In the eigenmode function shown in the figure, in addition to the main peak 103, there are two peaks 104 and 105. Obviously, the two peaks 104 and 105 can be decomposed into other eigenmode functions again. Therefore Figure 3 The intrinsic mode functions of are under-resolved.

[0108] The signal that actually needs to be decomposed contains noise, and the extreme points in the amplitude spectrum of each intrinsic mode function are dense. Only counting the number of extreme points cannot produce meaningful results. It is necessary to extract the envelope of the amplitude spectrum and count the number of maximum points in the envelope to determine the number of peaks. Morphological filtering is a global nonlinear filtering method that changes the local morphology of the signal by sliding the structure element globally.

[0109] This embodiment first performs a morphological closing operation on the amplitude spectrum 201, and a circular structure element is used as the structural element, and the width of the structural element is 1% of the width of the amplitude spectrum, the width of the amplitude spectrum is 500, and the width of the structural element is 5, to fill the part of the amplitude spectrum with small amplitude and violent fluctuations, and obtain the main peak envelope 202; then the main peak envelope 202 obtained is further subjected to a morphological closing operation, and a straight line structure element is used as the structural element, and the width of the structural element is 2% of the width of the amplitude spectrum, the width of the amplitude spectrum is 500, and the width of the structural element is 10, and the secondary peak generated by the first closing operation is smoothed to obtain the final envelope 203 of the amplitude spectrum; finally, when counting the number of maximum points, the maximum amplitude threshold is set to one thirtieth of the maximum value of the amplitude spectrum, the maximum value of the amplitude spectrum is 0.26, and the maximum amplitude threshold is 0.00867; the maximum interval threshold is set to one% of the width of the amplitude spectrum, the width of the amplitude spectrum is 500, and the maximum interval threshold is 5; by setting the threshold, the extreme points with small amplitudes are ignored, thereby achieving the purpose of removing noise. Figure 4 and Figure 5 As shown, Figure 4 This is a schematic diagram of the statistical maximum points of the amplitude spectrum after two morphological filtering. Figure 5 yes Figure 4 In the figure, 201 is the original amplitude spectrum, 202 is the result of the first closed operation, and 203 is the result of the second closed operation. Then the maximum point is found, and the result is the solid points 204, 205 and 206 in the figure. The number of maximum points is 3, indicating that 3 modal components are decomposed into the intrinsic mode function, and there is under-decomposition. It can be seen that the method of this embodiment can effectively remove the influence of noise and accurately count the number of maximum points.

[0110] When the original signal contains noise, variational mode decomposition does not need to include these noises in the reconstructed signal. At this time, the Lagrange multiplier remains unchanged during the iteration. Let the Lagrange multiplier be initially 0, then equation (3) only contains two parts. The first part constrains the severity of the change of each intrinsic mode function through regularization, which plays a role in bandwidth limitation; the second part limits the reconstructed signal from being too distorted compared with the original signal. The weight distribution of the two is determined by the input parameter α. The larger α is, the simpler each intrinsic mode function is, but the larger the residual between the reconstructed signal and the original signal is; the smaller α is, the more complex each intrinsic mode function is, and the smaller the energy loss of the reconstructed signal is compared with the original signal. Therefore, the energy loss can be reduced as much as possible under the premise of ensuring that each intrinsic mode function does not have under-decomposition. The objective function (6) of the present invention incorporates the input parameter α into the objective function. The objective function is the weighted sum of the number of extreme points obtained after filtering each intrinsic mode function and α. The weight makes the magnitude of the under-decomposition loss function much larger than the magnitude of the energy loss. At the beginning of the optimization, each eigenmode function contains a large number of maximum points. After several iterations, the number of maximum points is reduced to a small enough level. At this time, the optimization algorithm searches so that the number of extreme points does not increase while reducing α as much as possible.

[0111] This embodiment adopts the gravity search optimization algorithm. The first step is to determine the search range of parameters K and α: the search range of K is an integer between 2 and 20, and the search range of α is 10 to 2000; determine the total number of iterations N T is 50, the current number of iterations N C is 0. Step 2: Set P to 40, and select 40 different variational mode decomposition parameters (K, α) within the search range. Step 3: For all 40 variational mode decomposition parameters, calculate the objective function of each variational mode decomposition parameter. Step 4: Keep 10 variational mode decomposition parameters that obtain small objective function results; for each variational mode decomposition parameter of the remaining 30 that obtain large objective function results, use the optimization algorithm to obtain the updated parameters. Step 5: Current number of iterations N C Add 1; if the total number of iterations is N T Greater than the current number of iterations N C , go to step 3; if the total number of iterations N T Equal to the current number of iterations N C , then the automatic optimization of parameters ends, and the parameters that obtain the minimum objective function result are the optimal parameters. Finally, the optimal parameters K=11, α=24 are obtained, and the decomposition results are as follows Figure 6 shown.

[0112] In order to further illustrate the effect of this embodiment, the results of this embodiment are compared with those of literature 1 (Tang Guiji, Wang Xiaolong. Application of parameter optimization variational mode decomposition method in early fault diagnosis of rolling bearings. Journal of Xi'an Jiaotong University, 2015, Vol. 49, No. 5, pp. 73-8) and literature 2 (Zhang Shuqing, Li Jun, Jiang Anqi, Huang Jiao, Liu Haitao, Ai Hongke. A new two-stage short-term power load forecasting based on FPA-VMD and BiLSTM neural network. Power Grid Technology The method described in reference 1 is used for comparison. The optimization is performed according to the method described in reference 1, and the decomposition results are as follows: Figure 7 . While obtaining the eigenmode function with the smallest envelope entropy, no restrictions are imposed on other eigenmode functions, resulting in over-decomposition and under-decomposition. The optimization method described in Reference 2 is used to find the decomposition result. Figure 8 This method lacks the penalty for over-decomposition, and an inappropriately large K value is obtained. The optimization is performed according to the method described in Reference 3, and the decomposition result is as follows Fig. 9. Since the modal components with smaller permutation entropy tend to have simpler spectra and require larger α, the method minimizes the permutation entropy of each intrinsic mode function to maximize the decomposition parameter α, resulting in serious distortion of the reconstructed signal, and the high-frequency components in the spectrum are discarded, resulting in serious under-decomposition results. It can be seen that the method of this embodiment is better than the methods in the three existing literatures.

Claims

1. A method for determining parameters of variational mode decomposition, characterized in that The method is used for status monitoring of the original vibration signal of hydropower equipment, and comprises the following steps: The first step is to determine the total number of intrinsic mode function signals K and the search range of the penalty factor α, where K is a positive integer; determine the total number of iterations N T ; Current number of iterations N C Set to 0; Step 2: Select P different variational mode decomposition parameters (K, α) within the search range; The third step is to calculate the objective function of each variational mode decomposition parameter for all P variational mode decomposition parameters, where: The objective function L of the nth variational mode decomposition parameter is n The calculation formula is as follows: L n =a+10 4 ×n max ; Where n=1,2,…,P; n max The maximum number of maximum points contained in the amplitude spectrum of each intrinsic mode function signal after the variational mode decomposition is completed using the parameters; The max The calculation method is as follows: For each intrinsic mode function signal amplitude spectrum after the variational mode decomposition is completed using the parameters, find the number of maximum points it contains, where the largest number of maximum points is n max ; The signal amplitude spectrum of each intrinsic mode function is calculated by using a mathematical morphology method to find the number of maximum points contained therein. The specific method is as follows: Step (1) performing a morphological closing operation on the amplitude spectrum to obtain a main peak envelope, wherein: a structure element of the morphological closing operation is a circle, and the width of the structure element is one percent of the width of the amplitude spectrum; Step (2) performs a morphological closing operation on the main peak envelope obtained in step (1) to obtain a final envelope of the amplitude spectrum, wherein: the structural element of the morphological closing operation in step (2) is a straight line, and the width of the structural element is 2% of the width of the amplitude spectrum; Step (3) counting the number of maximum points on the final envelope of the amplitude spectrum obtained in step (2); Step 4: retain the parameters of the variational mode decomposition that obtains the small objective function result, and use the optimization algorithm for each variational mode decomposition parameter that obtains the large objective function result to obtain the updated parameters; Step 5: Current number of iterations N C Add 1; if the total number of iterations is N T Greater than the current number of iterations N C , go to step 3; if the total number of iterations N T Equal to the current number of iterations N C , then the automatic optimization of parameters ends, and the parameters that obtain the minimum objective function result are the optimal parameters.

2. The method for determining parameters of variational mode decomposition according to claim 1, characterized in that In the fourth step, the method of using the optimization algorithm to obtain the updated parameters is as follows: Suppose the parameter (K p ,α p ), p = 1, 2, ..., P; get the updated variational mode decomposition parameters (K ′ p ,α′ p ), using the gravity search optimization algorithm, the following calculation process is used: Nth C At the iteration, the gravitational constant G is: Among them, e is a natural constant; Among all P variational mode decomposition parameters, the quality of the nth variational mode decomposition parameter is denoted by M n ; Assume that the parameters of the ith variational mode decomposition are (K i ,α i ), the parameter of the jth variational mode decomposition is (K j ,α j ), is the parameter K of the ith variational mode decomposition i and the parameter K of the jth variational mode decomposition j The force between is the parameter α of the ith variational mode decomposition i and the parameter α of the jth variational mode decomposition j The force between For the pth variational mode decomposition parameter (K p ,α p ), parameter K p The total force F K and parameter α p The total force F α They are: Where R is a random number between 0 and 1; Then the parameter K p The acceleration a K and parameter α p The acceleration a α They are: Parameter K p The speed and parameter α p The speeds are: v K =R+a K ; in α =R×α+a α ; The new variational mode decomposition parameters are as follows: K p ′=[K p +v K ]; α′ p ＝α p +v α ； Among them, [K p +v K ] indicates K p +v K The calculation result is rounded off; When the variational mode decomposition parameter (K′ p ,α′ p ) and also satisfy: Condition 1: K′ p and α′ p Within the search range of parameters K and α respectively; Condition 2: K′ p ≠K p , and α′ p ≠α p When , the parameter of the pth variational mode decomposition (K p ,α p ) is updated to (K′ p ,α′ p ); Otherwise, repeat the above process of using the gravity search optimization algorithm to obtain the variational mode decomposition parameters until (K′ p ,α′ p ) satisfies the above two conditions.

3. The method for determining parameters of variational mode decomposition according to claim 2, characterized in that The M n The calculation formula is as follows: Among them, the intermediate variable m corresponding to the parameter of the nth variational mode decomposition n The calculation formula is as follows: Among them, L b is the minimum value of the objective function obtained by all P variational mode decomposition parameters; L w is the maximum value of the objective function obtained by all P variational mode decomposition parameters, L n is the value of the objective function obtained by the parameters of the nth variational mode decomposition; The intermediate variable m corresponding to the parameter of the jth variational mode decomposition j The calculation formula is as follows: Among them, L j is the value of the objective function obtained by the parameters of the j-th variational mode decomposition, j = 1, 2, …, P.

4. The method for determining parameters of variational mode decomposition according to claim 2, characterized in that Said and The calculation formula is as follows: Among them, R ij is the parameter of the ith variational mode decomposition (K i ,α i ) and the parameters of the jth variational mode decomposition (K j ,α j ) between the two.

5. The method for determining parameters of variational mode decomposition according to claim 4, characterized in that The R ij The calculation formula is as follows:

Images