Rolling bearing fault diagnosis method, medium and equipment

Abstract

The invention provides a rolling bearing fault diagnosis method, a medium and equipment, and relates to the technical field of large-scale unit bearing fault diagnos.The method comprises the steps that collected original life cycle vibration data of a rolling bearing are obtained; performing early warning processing on the vibration data by using a Beta early warning algorithm, and determining and obtaining a required diagnosis signal according to an early warning interval; determining an input parameter of variational mode decomposition by using the diagnosis signal, performing optimal parameter selection on the variational mode decomposition by combining an improved pelican optimization algorithm, and decomposing the diagnosis signal by combining the optimal parameter with the variational mode decomposition to obtain a plurality of intrinsic mode functions; and evaluating the plurality of intrinsic mode functions by using the correlation entropy in combination with the kurtosis so as to select a proper intrinsic mode function for reconstruction, and obtaining a final diagnosis result according to envelope spectrum analysis of a reconstructed signal. According to the method, the rolling bearing fault can be early warned and accurately diagnosed, and the accuracy, the real-time performance and the intelligent level of fault diagnosis are improved.

Description

Technical Field

[0001] The present invention relates to the technical field of bearing fault diagnosis for large units, and particularly to a rolling bearing fault diagnosis method, medium and device. Background Art

[0002] As a key component in large units, the health state of a rolling bearing directly affects the working performance and operation safety of the equipment. In industrial production, rolling bearings often bear large loads and operate in complex environments, making them prone to wear, cracks or other damages. Since rolling bearing faults usually have a long latency period, early fault warning and diagnosis can not only effectively reduce equipment downtime, but also avoid the serious consequences of sudden failures, reduce maintenance costs, and improve the overall operation efficiency of the equipment. Therefore, research on early fault warning and diagnosis of rolling bearings has important engineering practical significance. Traditional rolling bearing fault diagnosis methods rely on manual experience, with limited diagnostic accuracy and real-time performance, and are difficult to handle complex working conditions and multi-dimensional data analysis. They lack intelligence and adaptability, have a long processing cycle, and cannot effectively conduct early fault warning, making it difficult to meet the requirements of modern equipment for efficient, accurate and intelligent fault diagnosis. Summary of the Invention

[0003] The purpose of the present invention is to provide a rolling bearing fault diagnosis method, medium and device, aiming to efficiently extract the characteristic signals of rolling bearings, overcome the influence of noise interference and working condition changes, realize early warning and accurate diagnosis of rolling bearing faults, and improve the accuracy, real-time performance and intelligent level of fault diagnosis. The specific technical solutions are as follows:

[0004] A rolling bearing fault diagnosis method, the method comprising the following steps:

[0005] S100. Obtain the collected original life cycle vibration data of the rolling bearing;

[0006] S200. Use the Beta warning algorithm to perform warning processing on the vibration data, and determine the required diagnostic signal according to the warning interval;

[0007] S300. Use the diagnostic signal to determine the input parameters of variational mode decomposition, and combine the improved pelican optimization algorithm to select the optimal parameters for the variational mode decomposition. Decompose the diagnostic signal through the optimal parameters combined with variational mode decomposition to obtain multiple intrinsic mode functions;

[0008] S400. Use correlation entropy combined with kurtosis to evaluate multiple intrinsic mode functions to select appropriate intrinsic mode functions for reconstruction, and obtain the final diagnostic result according to the envelope spectrum analysis of the reconstructed signal.

[0009] Furthermore, the probability density function of the random variable is as follows:

[0010]

[0011] where x is the vibration data, a and b are the upper and lower boundaries of the x value, and the shape parameters γ and η are μ x the mean of x, and σ x is the standard deviation of x. If γ > 0, η > 0, and z = (x - a) / (b - a), then x is said to follow a Beta distribution, denoted as x ~ β(γ, η);

[0012] The specific steps of S200 are as follows:

[0013] S210. Normalize the vibration data, then solve the probability density function f(x, γ, η, a, b) of the vibration data, and use the least squares method to determine the two shape parameters γ and η of the Beta distribution;

[0014] S220. Draw the Beta distribution according to the obtained shape parameters γ and η, thereby calculating the two-sided α1 quantiles λ1 and λ2 of the Beta distribution, and further calculating the threshold interval expression including:

[0015]

[0016] l1 = λ1(max(x) - min(x)) + min(x),

[0017] l2 = λ2(max(x) - min(x)) + min(x)

[0018] where the two-sided α1 quantile λ1 is the 0.5α1 quantile of, and λ2 is expressed as the 1 - 0.5α1 quantile of, and α1 is defined as the error caused by noise interference during data acquisition, represents the random fault signal of the rolling bearing, represents the probability that the fault signal is less than a certain value λ1, and [l1, l2] is the threshold interval;

[0019] S230. By comparing and judging the input vibration time series data with the threshold interval in turn, output the vibration data point x t that first exceeds the threshold interval [l1, l2], which is the warning moment, and then intercept the vibration time series data within one cycle according to the warning moment as the diagnostic signal for the next fault diagnosis and identification.

[0020] Furthermore, in the step S300, the variational mode decomposition includes the following steps:

[0021] S310. According to different situations, variational mode decomposition decomposes the diagnostic signal, i.e., the time-domain signal f(t), into K intrinsic mode functions u k (t) with different center frequencies and finite bandwidths. The expression of variational mode decomposition is as follows:

[0022] u k (t) = A k (t) cos(φ k (t))

[0023] where t represents the time series, u k (t) represents the k-th mode component, A k (t) represents the instantaneous amplitude of the k-th mode component, φ k (t) represents the instantaneous phase of the k-th mode component, and cos(φ k (t)) represents the amplitude-modulated sine signal;

[0024] S320. The constrained variational model expression is:

[0025]

[0026] where f(t) represents the input signal, represents the time derivative operator, δ(t) represents the Dirac function, represents the core part of the Hilbert transform, w k represents the center frequency of the k-th mode component, {u k} represents the set of all mode components, {w k} represents the center frequencies of each mode component. In the formula, {u k}, {w k} represent the variable sets to be solved in the optimization process, and the two variables are adjusted to minimize the error;

[0027] S330. By introducing the quadratic penalty factor α and the Lagrange multiplier operator λ(t), the constrained variational model is transformed into an unconstrained variational model. The unconstrained variational model expression is:

[0028]

[0029] where L({u k}, {w k}, λ) represents the Lagrangian function, μ k represents the k-th mode component factor, λ represents the Lagrangian function factor, α represents the penalty factor, λ(t) represents the Lagrange multiplier, represents the inner product between the Lagrangian and the signal difference;

[0030] S340. Continuously update \(u\) in the above formula using the alternating direction multiplier algorithm k , \(w\) k such that the original signal \(f(t)\) is effectively decomposed into \(K\) intrinsic mode functions.

[0031] Furthermore, taking the minimum value of the modal weighted entropy as the fitness function, an improved pelican optimization algorithm is used to autonomously optimize the variational mode decomposition parameters to obtain the optimal mode functions and penalty factors. The diagnostic signal is subjected to variational mode decomposition using the optimal mode functions and penalty factors to obtain \(K\) intrinsic mode functions.

[0032] Furthermore, the expression of the improved pelican optimization algorithm includes:

[0033]

[0034] where \(SampEn(m, r, N)\) represents the modal sample entropy, \(EnpEn\) represents the modal envelope entropy, \(N\) represents the length of \(u(t)\) obtained by decomposing the diagnostic sequence, \(m\) is the embedding dimension, \(r\) is the similarity tolerance, \(B\) represents the number of template pairs with a template length of \((m + 1)\) and that match each other, \(C\) represents the number of template pairs with a template length of \(m\) and that match each other; \(p\) k represents the normalized probability distribution of the signal \(u(j)\), \(\chi(j)\) k,j represents the envelope signal obtained by Hilbert transforming the signal \(u(j)\); k k k k k The expression of the fitness function is:

[0035] F(u

[0036] ) = \(\varepsilon_1 SampEn(m, r, N)+\varepsilon_2 EnpEn\) k

[0037] where \(F(u\) k ) represents the fitness function value, \(\varepsilon_1,\varepsilon_2\) represent the weight coefficients of the two, and satisfy \(\varepsilon_1+\varepsilon_2 = 1\);

[0038] The expression for initializing the input random variables of the pelican optimization algorithm using the good point set theory is:

[0039] r d = 2cos(2πd / g), 1 ≤ d ≤ s

[0040] where \(d\) represents a real number within the interval range, \(g\) represents the smallest prime number point such that \((g - 3) / 2\geq s\), \(s\) represents the dimension in Euclidean space, \(r\) d represents the point after passing through the good point set;

[0041] The calculation expression in the exploration stage of the improved pelican optimization algorithm is:

[0042]

[0043] Among them, represents the position of the \(p_i\)-th pelican in the \(p_j\)-th dimension after the exploration phase update. \(I\) represents a random integer 1 or 2, and \(o\) pj represents the position of the prey in the \(p_j\)-th dimension, and \(F\) pi represents the function value of the \(p_i\)-th pelican, and \(F\) p represents the objective function value of the \(p_i\)-th pelican, represents the non-linear increasing inertia weight factor. \(N_t\) represents the current iteration number, and \(N_T\) represents the maximum iteration number, and \(h\) pi,pj represents the position of the \(p_i\)-th pelican in the \(p_j\)-th dimension before the exploration phase update, and \(x\) pi,pj represents the value of the variable of the \(p_i\)-th pelican in the \(p_j\)-th dimension before the exploration phase update. \(rand(0, 1)\) represents a random number in the interval \([0, 1]\);

[0044] The position update expression in the exploration phase is as follows:

[0045]

[0046] Among them, represents the updated position of the \(p_i\)-th pelican in the exploration phase, represents the objective function value of the \(p_i\)-th pelican at this position in the exploration phase, and \(o\) pi represents the position when the objective function value of the \(p_i\)-th pelican after the update in the exploration phase is not greater than that before the update, and \(O\) pi represents the position of the \(p_i\)-th pelican after meeting the condition;

[0047] The calculation expression in the water surface flight phase in the improved pelican optimization algorithm is as follows:

[0048]

[0049] Among them, represents the position of the \(i\)-th pelican in the \(p_j\)-th dimension after the update in the water surface flight phase, and \(o\) pi, x j represents the position of the \(p_i\)-th pelican in the \(p_j\)-th dimension, and \(R\) represents a random integer 1 or 2;

[0050] The position update expression in the water surface flight phase is as follows:

[0051]

[0052] Among them, represents the updated position of the \(p_i\)-th pelican in the water surface flight phase, The objective function value of the pi-th pelican at this position representing the water surface flight phase.

[0053] Furthermore, using the correlation entropy combined with kurtosis to evaluate multiple intrinsic mode functions to select appropriate intrinsic mode functions for reconstruction includes the following steps:

[0054] S410. For the intrinsic mode function u k (t) obtained by decomposition, calculate its absolute value component y(t) = |u k (t)|, and use the kernel density estimation method to solve the probability density function F(y) of y(t). The calculation formula is:

[0055]

[0056] where d is the bandwidth parameter in kernel density estimation;

[0057] S420. Using the range of y(t) for which the probability density function is greater than the critical probability φ, 0 ≤ φ ≤ 1, the expression for the minimum value of y(t) within this range is:

[0058]

[0059] where y low (t) represents the minimum value within the range, and f(s) represents the absolute value component of u k (t);

[0060] S430. Using the linear interpolation algorithm to replace the points greater than the said y low to obtain a new modal component

[0061] S440. Calculate the kurtosis value according to the new modal component . The calculation formula is:

[0062]

[0063] where K adj represents the kurtosis value calculated according to the newly obtained modal component;

[0064] S450. Using the kurtosis value combined with the correlation entropy to calculate an improved kurtosis index. The expression includes:

[0065]

[0066] where MCC(e i ) represents the maximum correlation entropy, CEIK(i) is the improved kurtosis index, and e i = x i - y i represents the error signal, x i , yi represent two different modal components, κ σ is a Gaussian kernel function, and MCC'(i) represents the transposed matrix of the maximum correlation entropy;

[0067] S460. Taking the preset percentage of the maximum value component of the improved kurtosis index as the threshold, select the intrinsic mode functions greater than the threshold among multiple intrinsic mode functions for reconstruction to obtain the reconstructed signal.

[0068] Furthermore, the preset percentage is 70%.

[0069] Furthermore, perform Hilbert transform analysis on the reconstructed signal to obtain the envelope spectrum, and diagnose the final fault category according to whether the peak spectral values on the envelope spectrum are consistent with the theoretical frequencies of different fault characteristics of the rolling bearing.

[0070] The present invention also provides a computer-readable storage medium, on which a computer program is stored, and when the computer program is executed by a processor, the steps of the rolling bearing fault diagnosis method described above are implemented.

[0071] The present invention also provides a computer device, including a memory, a processor, and a computer program stored on the memory and executable on the processor. When the processor executes the program, the steps of the rolling bearing fault diagnosis method described above are implemented.

[0072] A rolling bearing fault diagnosis method, medium and device provided by the present invention have the following beneficial effects:

[0073] By acquiring the original life cycle vibration data of the collected rolling bearing; using the Beta warning algorithm to perform warning processing on the vibration data, and determining the required diagnostic signal according to the warning interval; using the diagnostic signal to determine the input parameters of the variational mode decomposition, and combining the improved pelican optimization algorithm to perform optimal parameter selection for the variational mode decomposition, decomposing the diagnostic signal through the optimal parameters combined with the variational mode decomposition to obtain multiple intrinsic mode functions; using correlation entropy combined with kurtosis to evaluate multiple intrinsic mode functions to select appropriate intrinsic mode functions for reconstruction, and obtaining the final diagnostic result according to the envelope spectrum analysis of the reconstructed signal; it can efficiently extract the characteristic signals of the rolling bearing, overcome the influence of noise interference and working condition changes, realize the early warning and accurate diagnosis of rolling bearing faults, improve the accuracy, real-time and intelligent level of fault diagnosis, and has strong adaptability and broad application prospects. Description of the Drawings

[0074] Figure 1 is a schematic flow chart of a rolling bearing fault diagnosis method provided by the present invention;

[0075] Figure 2It is the overall flowchart of the embodiment of the present invention;

[0076] Figure 3 It is the vibration data time series diagram provided by the embodiment of the present invention;

[0077] Figure 4 It is the warning point diagram obtained by evaluating and diagnosing the Beta warning algorithm provided by the embodiment of the present invention;

[0078] Figure 5 It is the vibration data time series diagram within one cycle at the warning point moment intercepted after passing through the Beta warning algorithm provided by the embodiment of the present invention;

[0079] Figure 6 It is the envelope spectrum diagram before the Beta-IPOA-VMD fault diagnosis provided by the embodiment of the present invention;

[0080] Figure 7 It is the convergence comparison curve diagram of the fitness function of the IPOA algorithm provided by the embodiment of the present invention;

[0081] Figure 8 It is the comparison diagram of the improved kurtosis index quantity of the correlation entropy and the kurtosis index quantity of the correlation entropy of each modal component obtained by VMD decomposition provided by the embodiment of the present invention;

[0082] Figure 9 It is the envelope spectrum diagram after the Beta-IPOA-VMD fault diagnosis provided by the embodiment of the present invention;

[0083] Figure 10 It is the structural block diagram of the computer device provided by the embodiment of the present invention. Detailed implementation manners

[0084] Next, in combination with the accompanying drawings provided by the present invention, the technical solutions in the embodiments of the present invention will be clearly and completely described. According to the following description, the advantages and features of the present invention will be clearer. It should be noted that the accompanying drawings are all in a very simplified form and use non-precise scales, only for the purpose of facilitating and clearly assisting in explaining the purpose of the embodiments of the present invention.

[0085] Embodiment 1

[0086] This embodiment provides a rolling bearing fault diagnosis method. Referring to Figure 1 、 2 as shown, the method includes the following steps:

[0087] S100. Obtain the original life cycle vibration data of the collected rolling bearing.

[0088] Specifically, according to the failure mode during the operation of the rolling bearing, a high-precision sensor is installed at a suitable position to collect the original life cycle vibration data of the rolling bearing.

[0089] The accelerated life test data of a rolling bearing from a domestic university are used to verify the effectiveness of the early warning and diagnosis methods proposed in this paper. The tested bearing is an LDK UER204 rolling bearing. The test sampling frequency is 25.6 kHz, the sampling duration for each time is 1.28 s, and the sampling interval is 60 s. The full-life cycle test data of the outer ring fault of the drive-end bearing are selected for verification. The bearing speed is 2100 r / min, and a radial force of 12 kN is applied to the bearing by the hydraulic loading system. The tested bearing runs until complete failure, and a total of 123 samples are obtained. The time series diagram of the vibration signal is as Figure 3 shown, Figure 3 which represents the amplitude diagram of the vibration signal changing with time. The abscissa represents the number of points, and the ordinate represents the amplitude. It can be clearly seen from the figure that the data starts to fluctuate greatly at [a certain point]. Therefore, it can be preliminarily guessed that the rolling bearing is in a fault state at this point. The specific technical parameter information of the tested bearing is shown in Table 1 below:

[0090] Table 1 Specific technical parameter table of the bearing

[0091]

[0092] S200. Use the Beta early warning algorithm to perform early warning processing on the vibration data, and determine the required diagnostic signal according to the early warning interval.

[0093] Traditional early warning processing algorithms have certain requirements for system hardware in application, and there may be subjectivity in judging the early warning points of data faults. Especially when setting the early warning threshold interval, it is easily affected by human interference, resulting in inaccurate and unreliable early warning results. To eliminate subjective factors and overcome the deficiencies of traditional early warning algorithms, the present invention uses the Beta early warning algorithm to perform early warning processing on rolling bearing data. This algorithm automatically evaluates the early warning interval, and intelligently determines the fault early warning point according to the historical data and real-time operating status of the rolling bearing, so as to accurately identify the time series signal that needs to be diagnosed. In this way, it can effectively reduce the error of manually setting the threshold, avoid the uncertainty brought by subjective judgment, and improve the accuracy, real-time performance and robustness of fault early warning. In addition, the application of the Beta early warning algorithm can adaptively adjust the early warning model according to different working conditions and equipment operating conditions, further improving the reliability and practicality of the early warning system.

[0094] The steps of the Beta warning algorithm are as follows: First, fit the Beta distribution model based on historical data to determine the probability distributions of the system in normal and faulty states; then, automatically set the warning threshold interval according to the distribution characteristics and evaluate the matching degree between the device operation data and historical data in real time; when the current data exceeds the warning interval, the algorithm determines it as a potential fault point and issues a warning signal, indicating that the device may malfunction. This process can automatically adjust parameters, reduce subjective intervention, and improve the accuracy and real-time performance of warnings.

[0095] In one embodiment, the probability density function of the random variable is:

[0096]

[0097] where x is the vibration data, a and b are the upper and lower boundaries of the x value, and the shape parameters γ and η are μ x the mean of x, and σ x is the standard deviation of x, γ > 0, η > 0, z = (x - a) / (b - a), then x is said to follow a Beta distribution, denoted as x ~ β(γ, η);

[0098] The specific steps of S200 are as follows:

[0099] S210. Normalize the vibration data, then solve the probability density function f(x, γ, η, a, b) of the vibration data, and use the least squares method to determine the two shape parameters γ and η of the Beta distribution;

[0100] S220. Draw the Beta distribution according to the obtained shape parameters γ and η, thereby calculate the two-sided α1 quantiles λ1 and λ2 of the Beta distribution, and further calculate the threshold interval expression including:

[0101]

[0102] l1 = λ1(max(x) - min(x)) + min(x),

[0103] l2 = λ2(max(x) - min(x)) + min(x)

[0104] where the two-sided α1 quantile λ1 is the 0.5α1 quantile of, and λ2 is expressed as the 1 - 0.5α1 quantile of, and α1 is defined as the error caused by noise interference during data acquisition, represents the random fault signal of the rolling bearing, represents the probability that the fault signal is less than a certain value λ1; [l1, l2] is the threshold interval.

[0105] S230. By sequentially comparing the input vibration time series data with the threshold interval, output the vibration data point x that first exceeds the threshold interval [l1, l2]. t , which is the warning moment. Then, intercept the vibration time series data within one cycle according to the warning moment as the diagnostic signal for the next fault diagnosis and identification.

[0106] In a preferred embodiment, α1 = 5%.

[0107] As Figure 4 shown, where Figure 4 represents the vibration data time series diagram after the Beta warning algorithm. The two red dotted lines in the figure are the upper and lower limits of the warning interval. It can be clearly seen from the figure that the vibration data can exceed the warning interval at 7.209×10 5 . Therefore, it can be determined that the vibration data starts to fail here.

[0108] As Figure 5 shown, where Figure 5 represents the time series diagnostic signal of one cycle after the fault point. The signal data in the figure fluctuates randomly and is difficult to observe. It is not easy to diagnose from the envelope spectrum Figure 6 of the diagnostic signal either. Figure 6 is the envelope spectrum diagram of the diagnostic signal, which represents the envelope spectrum diagram drawn after the diagnostic signal undergoes Hilbert-Huang transform. In the figure, the 1× frequency outer race fault signal frequency and the 2× frequency outer race fault signal frequency with relatively prominent peaks can be faintly identified, and the fluctuation range of other frequencies is relatively large and difficult to distinguish. Therefore, further diagnosis is required.

[0109] Through the above processing steps, the parameters can be automatically adjusted, subjective intervention can be reduced, and the accuracy and real-time performance of the warning can be improved.

[0110] S300. Use the diagnostic signal to determine the input parameters of variational mode decomposition, combine the improved pelican optimization algorithm to select the optimal parameters for the variational mode decomposition, and decompose the diagnostic signal through the optimal parameters combined with variational mode decomposition to obtain multiple intrinsic mode functions (IMFs).

[0111] Variational Mode Decomposition (VMD) is a signal decomposition method based on variational optimization theory, aiming to decompose complex signals into multiple sub-signals with different frequency characteristics. Its brief steps include: First, through an iterative optimization process, the original signal is decomposed into several mode functions, and each mode function corresponds to a specific frequency band. Then, VMD automatically adjusts the frequency range of each mode according to the set central frequency and bandwidth to accurately extract different frequency band information in the signal, avoiding the problems of frequency overlap or mode aliasing in traditional signal processing methods.

[0112] The main advantage of VMD is its ability to process non-stationary signals and effectively reduce noise, thereby retaining the key features in the signal. Compared with traditional signal decomposition methods, VMD has strong self-adaptability and can automatically select suitable decomposition parameters, avoiding errors caused by human intervention. In addition, when dealing with complex signals, VMD can accurately separate each frequency band, making it suitable for application fields such as fault diagnosis and condition monitoring, especially for signals with non-linearity and non-stationarity.

[0113] In one embodiment, the variational mode decomposition includes the following steps:

[0114] S310. According to different situations, variational mode decomposition decomposes the diagnostic signal, that is, the time-domain signal f(t), into K intrinsic mode functions u k (t) with different central frequencies and finite bandwidths. The expression of variational mode decomposition is:

[0115] u k (t) = A k (t)cos(φ k (t))

[0116] Among them, t represents the time series, u k (t) represents the k-th mode component, A k (t) represents the instantaneous amplitude of the k-th mode component, φ k (t) represents the instantaneous phase of the k-th mode component, and cos(φ k (t)) represents an amplitude-modulated sine signal;

[0117] S320. The constrained variational model expression is:

[0118]

[0119] Among them, f(t) represents the input signal, represents the time derivative operator, δ(t) represents the Dirac function, represents the core part of the Hilbert transform, w kdenotes the central frequency of the k-th modal component, {u k} denotes the set of all modal components, {w k} denotes the central frequency of each modal component. In the formula, {u k}, {w k} denote the set of variables to be solved in the optimization process. By adjusting these two variables, the error is minimized;

[0120] S330. Introduce the quadratic penalty factor α and the Lagrange multiplier operator λ(t) to transform the constrained variational model into an unconstrained variational model. The unconstrained variational model is expressed as:

[0121]

[0122] where L({u k}, {w k}, λ) represents the Lagrangian function, μ k denotes the k-th modal component factor, λ represents the Lagrangian function factor, α represents the penalty factor, and λ(t) represents the Lagrange multiplier. represents the inner product between the Lagrangian and the signal difference; By adjusting the Lagrange multiplier, the constraint condition is incorporated into the optimization process and it is ensured that this constraint condition is satisfied when minimizing the objective function;

[0123] S340. Continuously update u k , w k of the above formula using the alternating direction multiplier algorithm, so that the original signal f(t) is effectively decomposed into K intrinsic mode functions.

[0124] In one embodiment, with the minimum value of modal weighted entropy as the fitness function, the improved pelican optimization algorithm (IPOA) is used to perform autonomous optimization of the variational mode decomposition parameters to obtain the optimal modal function and penalty factor. The diagnostic signal is subjected to variational mode decomposition using the optimal modal function and penalty factor to obtain K intrinsic mode functions.

[0125] The Pelican Optimization Algorithm (POA) is a swarm intelligence optimization algorithm based on the foraging behavior of pelicans, which simulates the process of pelicans finding food through cooperation and competition during predation. The algorithm designs a fitness function to simulate the search behavior of individual pelicans, uses swarm collaboration and information sharing for global search, and gradually approaches the optimal solution of the problem. The advantages of the Pelican Optimization Algorithm lie in its strong global search ability and fast convergence speed, which can effectively avoid the problem of local optimal solutions. In addition, the algorithm is simple and easy to implement, with relatively few parameter adjustments, strong adaptability, and can achieve good performance in a variety of complex optimization problems, especially suitable for dealing with high-dimensional, non-linear, and multi-objective optimization problems.

[0126] In one embodiment, the expression of the improved Pelican Optimization Algorithm includes:

[0127]

[0128] Among them, SampEn(m, r, N) represents the modal sample entropy, EnpEn represents the modal envelope entropy, N represents the length of u k (t) obtained by decomposing the diagnostic sequence, m is the embedding dimension, r is the similarity tolerance, B represents the number of template pairs with a template length of (m + 1) and mutual matching, C represents the number of template pairs with a template length of m and mutual matching; p k,j represents the normalized probability distribution of the signal u k (j), χ k (j) represents the envelope signal obtained by Hilbert transform of the signal u k (j);

[0129] The expression of the fitness function is:

[0130] F(u k ) = ε1SampEn(m, r, N)) + ε2EnpEn

[0131] Among them, F(u k ) represents the fitness function value, ε1, ε2 represent the weight coefficients of the two, and satisfy ε1 + ε2 = 1;

[0132] The expression for initializing the input random variables of the Pelican Optimization Algorithm using the good point set theory is:

[0133] r d = 2cos(2πd / g), 1 ≤ d ≤ s

[0134] Among them, d represents a real number within the interval range, g represents the smallest prime number point where (g - 3) / 2 ≥ s, s represents the dimension in Euclidean space, r d represents the point after passing through the good point set;

[0135] In the improved pelican optimization algorithm, the calculation expression in the exploration stage is as follows:

[0136]

[0137] Wherein, represents the position of the pi-th pelican in the pj-th dimension after update in the exploration stage, I represents a random integer 1 or 2, o pj represents the position of the prey in the pj-th dimension, F pi represents the function value of the pi-th pelican, F p represents the objective function value of the pi-th pelican, represents the non-linear increasing inertia weight factor, Nt represents the current iteration number, NT represents the maximum iteration number, h pi,pj represents the position of the pi-th pelican in the pj-th dimension before update in the exploration stage, x pi,pj represents the value of the variable of the pi-th pelican in the pj-th dimension before update in the exploration stage, rand(0, 1) represents a random number in the interval [0, 1];

[0138] The position update expression in the exploration stage is as follows:

[0139]

[0140] Wherein, represents the updated position of the pi-th pelican in the exploration stage, represents the objective function value of the pi-th pelican at this position in the exploration stage, o pi represents the position when the objective function value of the pi-th pelican after update in the exploration stage is not greater than that before update, O pi represents the position of the pi-th pelican after meeting the condition;

[0141] In the improved pelican optimization algorithm, the calculation expression in the surface flight stage is as follows:

[0142]

[0143] Wherein, represents the position of the i-th pelican in the pj-th dimension after update in the surface flight stage, o pi,pj represents the position of the pi-th pelican in the pj-th dimension, R represents a random integer 1 or 2;

[0144] The position update expression in the surface flight stage is as follows:

[0145]

[0146] Wherein, The updated position of the pi-th pelican representing the water surface flight stage, The objective function value of the pi-th pelican at this position representing the water surface flight stage.

[0147] Through the improved pelican optimization algorithm, a better fitness value is obtained, as Figure 7 shown. The figure shows the curve of the fitness value change of the improved pelican optimization algorithm with the change of the number of iterations. At the same time, the pelican optimization algorithm and the particle swarm optimization algorithm with good performance are compared. Among them, it can be seen from the figure that the improved pelican optimization algorithm can reach the minimum fitness value faster and has already been the minimum in the second iteration; while the other compared algorithms need 4 iterations and the fitness value is not the minimum. Therefore, it can be obtained that the performance of the improved pelican optimization algorithm is better.

[0148] S400. Use correlation entropy combined with kurtosis to evaluate multiple intrinsic mode functions to select appropriate intrinsic mode functions for reconstruction, and obtain the final diagnosis result according to the envelope spectrum analysis of the reconstructed signal.

[0149] In one embodiment, using correlation entropy combined with kurtosis to evaluate multiple intrinsic mode functions to select appropriate intrinsic mode functions for reconstruction includes the following steps:

[0150] S410. For the intrinsic mode function u k (t) obtained by decomposition, calculate its absolute value component y(t) = |u k (t)|, and use the kernel density estimation method to solve the probability density function F(y) of y(t). The calculation formula is:

[0151]

[0152] where d is the bandwidth parameter in the kernel density estimation;

[0153] S420. Use the range of y(t) for which the probability density function is greater than the critical probability φ, 0 ≤ φ ≤ 1. The minimum value expression of y(t) within this range is:

[0154]

[0155] where y low (t) represents the minimum value within the range, and f(s) represents the absolute value component of u k (t);

[0156] S430. Use the linear interpolation algorithm to replace the points greater than the said y low to obtain a new modal component

[0157] S440. According to the new modal component Calculate the kurtosis value, and the calculation formula is:

[0158]

[0159] Among them, K adj represents the kurtosis value calculated according to the newly obtained modal component;

[0160] The kurtosis is used to measure the kurtosis and tail thickness of the data distribution, and is used to describe the distribution form, detect outliers, and evaluate risks or non-Gaussian characteristics.

[0161] S450. Calculate the improved kurtosis index (CEIK) by combining the kurtosis value with the correlation entropy. The expression includes:

[0162]

[0163] Among them, MCC(e i ) represents the maximum correlation entropy, CEIK(i) is the improved kurtosis index, and e i =x i -y i represents the error signal, x i , y i represent two different modal components, κ σ is the Gaussian kernel function, and MCC'(i) represents the transpose matrix of the maximum correlation entropy;

[0164] The CEIK index is used to evaluate the modal components. The larger the index, the more effective information the corresponding intrinsic mode function component contains. Use the CEIK index to evaluate the modal components. As Figure 8 shown in the figure, the CEIK indexes of different modal components are shown in the figure, and the CEK indexes are compared. Among them, the CEIK indexes of IMF3, IMF6, and IMF7 are significantly larger than those of other IMFs and are relatively easy to distinguish; while the CEK indexes of IMF1, IMF2, IMF3, IMF4, IMF5, and IMF6 are almost close and difficult to distinguish, and only the CEK index of IMF7 is larger, resulting in difficult distinction or omission of effective modal components. Therefore, it can be shown that the CEIK index is more effective in screening modal components.

[0165] S460. Take the preset percentage of the component with the maximum improved kurtosis index as the threshold, and select the intrinsic mode functions greater than the threshold among multiple intrinsic mode functions for reconstruction to obtain the reconstructed signal.

[0166] In a preferred embodiment, the preset percentage is 70%.

[0167] In one embodiment, the final diagnostic result obtained by envelope spectrum analysis of the reconstructed signal includes the following steps:

[0168] S470. Perform Hilbert transform analysis on the reconstructed signal to obtain the envelope spectrum, and diagnose the final fault category based on whether the peak spectral values on the envelope spectrum are consistent with the theoretical frequencies of different fault characteristics of the rolling bearing.

[0169] As Figure 9 shown, the figure shows the envelope spectrogram of the diagnostic signal after the method of the present invention. It can be clearly seen from the figure that the 1-fold and 2-fold rotation frequencies f r of the rolling bearing are relatively prominent, the 1-fold, 2-fold, 3-fold, and 4-fold outer ring fault diagnosis frequencies f o and the modulation frequency components f o - f r and f o + f r centered on the theoretical fault frequency and with the rotation frequency as the sideband, and the noise frequency fluctuation is significantly weakened. Therefore, the fault occurring in the diagnostic signal can be judged as an outer ring fault.

[0170] A rolling bearing fault diagnosis method based on Beta - IPOA - VMD provided by an embodiment of the present invention. To verify the effectiveness of this method, we took the accelerated life test data of a rolling bearing in a domestic university as the background and obtained sample data under different working conditions. Compare the iterative effects of the traditional optimization algorithm and the IPOA - VMD optimization algorithm in the present invention. The results show that the method of the present invention is superior to the traditional optimization algorithm in terms of the accuracy and convergence speed of fault diagnosis. During the iteration process, the IPOA - VMD method effectively optimizes the parameter selection of VMD through the improved pelican optimization algorithm, significantly improving the accuracy of modal decomposition and the reliability of the diagnosis result. The specific comparison results are as Figure 7 shown. The iterative curve in the figure shows that this method can reach a lower fitness value within fewer iteration times, indicating a better optimization effect. In addition, this method also shows good fault warning ability and signal denoising effect in the case of large noise. The result comparison is as Figure 6 and Figure 9 shown.

[0171] The experimental results show that the rolling bearing fault diagnosis method based on Beta-IPOA-VMD of the present invention has remarkable effectiveness in the fault diagnosis of bearings in large units. By introducing the Beta early warning algorithm, the early warning of bearings is successfully realized, and the improved pelican optimization algorithm (IPOA) combined with the correlated entropy improved kurtosis (CEIK) index optimizes the selection of VMD parameters. Compared with the traditional POA and PSO algorithms, this method has obvious advantages in optimization performance, with the number of iterations reduced by about 20% and the fitness value decreased by about 15%. The CEIK index effectively improves the signal-to-noise ratio (SNR) of the signal by about 25% and reduces the root mean square error (RMSE) by 5%. The experimental results also show that in the noise background, this method can give early warning in time and effectively reduce noise, extract obvious fault features, significantly superior to the VMD method, with the SNR increased by 18% and the RMSE reduced by 3%. These results indicate that this method has good prospects in the application of rolling bearing fault early warning and diagnosis.

[0172] In summary, the fault diagnosis method proposed by the present invention effectively improves the accuracy and reliability of diagnosis by introducing the Beta early warning algorithm and the IPOA-VMD noise reduction technology. The experimental results verify the superior performance of this method in complex environments, showing its advantages in fault early warning and signal processing, and having broad application potential and good development prospects.

[0173] Embodiment 2

[0174] This embodiment provides a computer-readable storage medium, on which a computer program is stored, and when the computer program is executed by a processor, the steps of the rolling bearing fault diagnosis method described above are implemented.

[0175] Among them, the storage medium can be a magnetic disk, an optical disk, a read-only memory (ROM), a random access memory (RAM), a flash memory, a hard disk drive (HDD) or a solid-state drive (SSD), etc.; the storage medium can also include a combination of the above types of memories.

[0176] Embodiment 3

[0177] This embodiment provides a computer device, which includes: a memory, a processor, and a computer program stored on the memory and executable on the processor, and when the processor executes the program, the steps of the rolling bearing fault diagnosis method described above are implemented.

[0178] Such as Figure 10As shown, the computer device may include: at least one processor 71, such as a CPU (Central Processing Unit), at least one communication interface 73, a memory 74, and at least one communication bus 72. Among them, the communication bus 72 is used to implement the connection and communication between these components. Among them, the communication interface 73 may include a display screen and a keyboard. Optionally, the communication interface 73 may also include a standard wired interface and a wireless interface. The memory 74 may be a high-speed RAM memory (Random Access Memory, volatile random access memory), or a non-volatile memory, such as at least one disk memory. Optionally, the memory 74 may also be at least one storage device located far from the aforementioned processor 71. Among them, an application program is stored in the memory 74, and the processor 71 calls the program code stored in the memory 74 to execute any of the above method steps.

[0179] Among them, the communication bus 72 may be a peripheral component interconnect (PCI) bus or an extended industry standard architecture (EISA) bus, etc. The communication bus 72 can be divided into an address bus, a data bus, a control bus, etc. For the sake of convenience of representation, Figure 10 only a thick line is used to represent it in the figure, but it does not mean that there is only one bus or one type of bus.

[0180] Among them, the memory 74 may include a volatile memory, such as a random-access memory (abbreviation: RAM); the memory may also include a non-volatile memory, such as a flash memory, a hard disk drive (abbreviation: HDD), or a solid-state drive (abbreviation: SSD); the memory 74 may also include a combination of the above types of memories.

[0181] Among them, the processor 71 may be a central processing unit (abbreviation: CPU), a network processor (abbreviation: NP), or a combination of a CPU and an NP.

[0182] Among them, the processor 71 may further include a hardware chip. The above-mentioned hardware chip may be an application-specific integrated circuit (ASIC), a programmable logic device (PLD), or a combination thereof. The above-mentioned PLD may be a complex programmable logic device (CPLD), a field-programmable gate array (FPGA), a generic array logic (GAL), or any combination thereof.

[0183] Optionally, the memory 74 is further configured to store program instructions. The processor 71 may call the program instructions to implement the rolling bearing fault diagnosis method of the present invention.

[0184] Those skilled in the art of the present technology should understand that the present invention may be implemented in many other specific forms without departing from the spirit and scope of the present invention. Based on the embodiments of the present invention, any changes and modifications made by those of ordinary skill in the art of the present invention according to the above disclosure are within the protection scope of the claims.

Claims

1. A rolling bearing fault diagnosis method, characterized in that: The method comprises the following steps: S100, acquiring collected original life cycle vibration data of the rolling bearing; S200, using the Beta warning algorithm to perform warning processing on the vibration data, and obtaining a required diagnostic signal according to the warning interval; S300, using the diagnostic signal to determine input parameters of variational mode decomposition, combining the improved Pelican optimization algorithm to select optimal parameters for the variational mode decomposition, and decomposing the diagnostic signal by combining the optimal parameters with variational mode decomposition to obtain multiple intrinsic mode functions; S400, using correlation entropy combined with kurtosis to evaluate multiple intrinsic mode functions to select a suitable intrinsic mode function for reconstruction, and obtaining a final diagnosis result according to envelope spectrum analysis of the reconstructed signal.

2. The rolling bearing fault diagnosis method according to claim 1, characterized in that: The density function of a random variable is: Among them, x is the vibration data, a and b are the upper and lower boundaries of x, and the shape parameters γ and η are μ x is the mean of x, σ x is the standard deviation of x, Then x is said to obey Beta distribution, denoted as x~β(γ,η); The S200 specifically includes the following steps: S210, normalizing the vibration data, then solving the probability density function f(x, γ, η, a, b) of the vibration data, and determining two shape parameters γ and η of the Beta distribution using the least squares method; S220. Draw the Beta distribution according to the obtained shape parameters γ and η, so as to calculate the bilateral α1 quantiles λ1 and λ2 of the Beta distribution, and then calculate the threshold interval expression including: l1=λ1(max(x)-min(x))+min(x), l2=λ2(max(x)-min(x))+min(x) Among them, the bilateral α1 quantile λ1 is The 0.5α1 quantile, λ2 is expressed as 1-0.5α1 quantile, α1 is defined as the error caused by noise interference during data acquisition. Indicates random fault signals of rolling bearings, represents the probability that the fault signal is less than a certain value λ1, [l1,l2] is the threshold interval; S230, by comparing the input vibration time series data with the threshold interval in turn, outputting the vibration data point x that exceeds the threshold interval [l1, l2] for the first time t , which is the warning moment, then the vibration time series data within a cycle is intercepted according to the warning moment as the diagnostic signal for the next step of fault diagnosis and identification.

3. The rolling bearing fault diagnosis method according to claim 1, characterized in that: In step S300, the variational mode decomposition includes the following steps: S310, variational mode decomposition: according to different situations, the diagnostic signal, that is, the time domain signal f(t), is decomposed into K eigenmode functions u with different center frequencies and limited bandwidths. k (t), the expression of variational mode decomposition is: u k (t)=A k (t)cos(φ k (t)) Among them, t represents the time series, u k (t) represents the kth modal component, A k (t) represents the instantaneous amplitude of the kth modal component, φ k (t) represents the instantaneous phase of the kth modal component, cos(φ k (t)) represents an amplitude modulated sinusoidal signal; S320, the constrained variational model expression is: Where f(t) represents the input signal, represents the time derivative operator, δ(t) represents the Dirac function, represents the core part of Hilbert transform, w k represents the center frequency of the kth modal component, {u k } represents the set of all modal components, {w k } represents the center frequency of each modal component, where {u k },{w k } represents the set of variables that need to be solved in the optimization process, and the error is minimized by adjusting these two variables; S330, introducing a quadratic penalty factor α and a Lagrange multiplication operator λ(t), and transforming the constrained variational model into an unconstrained variational model. The unconstrained variational model is expressed as: Among them, L({u k },{w k },λ) represents the Lagrangian function, μ k represents the kth modal component factor, λ represents the Lagrangian function factor, α represents the penalty factor, λ(t) represents the Lagrangian multiplier, represents the inner product between the Lagrangian and the signal difference; S340, using the alternating direction multiplier algorithm to continuously update u in the above formula k ,w k , so that the original signal f(t) is effectively decomposed into K intrinsic mode functions.

4. The rolling bearing fault diagnosis method according to claim 3, characterized in that: The minimum value of modal weighted entropy is taken as the fitness function, and the improved Pelican optimization algorithm is used to autonomously optimize the variational modal decomposition parameters to obtain the optimal modal function and penalty factor. The optimal modal function and penalty factor are used to perform variational modal decomposition on the diagnostic signal to obtain K intrinsic mode functions.

5. The rolling bearing fault diagnosis method according to claim 4, characterized in that: The expression of the improved Pelican optimization algorithm includes: Among them, SampEn(m,r,N) represents the modal sample entropy, EnpEn represents the modal envelope entropy, and N represents the u obtained by decomposing the diagnostic sequence. k (t), m is the embedding dimension, r is the similarity tolerance, B is the number of template pairs with a template length of (m+1) and matching each other, and C is the number of template pairs with a template length of m and matching each other; p k,j Represents the signal u k The normalized probability distribution of (j), χ k (j) represents the signal u k (j) Envelope signal obtained by Hilbert transform; The fitness function expression is: F(u k )=ε1SampEn(m,r,N)+ε2EnpEn Among them, F(u k ) represents the fitness function value, ε1, ε2 represent the weight coefficients of the two, and satisfy ε1+ε2=1; The expression for initializing the input random variables of the Pelican optimization algorithm using the good point set theory is: r d =2cos(2πd / g),1≤d≤s Among them, d represents a real number in the interval, g represents the smallest prime number point (g-3) / 2≥s, s represents the dimension in Euclidean space, and r d represents the point after passing through the good point set; The calculation expression of the exploration phase in the improved Pelican optimization algorithm is: in, represents the position of the pi-th pelican in the pj-th dimension after the exploration phase update, I represents a random integer 1 or 2, o pj represents the position of the prey in the pjth dimension, F pi represents the function value of the pi-th pelican, F p represents the objective function value of the pi-th pelican, represents the nonlinear increasing inertia weight factor, Nt represents the current number of iterations, NT represents the maximum number of iterations, and h pi,pj represents the position of the pi-th pelican in the pj-th dimension before the exploration phase update, x pi,pj represents the value of the variable of the pjth dimension before the update of the pith pelican in the exploration phase, and rand(0,1) represents a random number in the interval [0,1]; The position update expression in the exploration phase is: in, represents the updated position of the pi-th pelican during the exploration phase, represents the objective function value of the pi-th pelican at this location during the exploration phase, o pi represents the position of the pi-th pelican in the exploration phase when the objective function value after the updated position is not greater than the objective function value before the update, O pi Indicates the position of the pi-th pelican after the conditions are met; The calculation expression of the surface flight phase in the improved Pelican optimization algorithm is: in, represents the pj-th dimension position of the i-th pelican after updating during the surface flight phase, o pi,pj represents the position of the pi-th pelican in the pj-th dimension, and R represents a random integer 1 or 2; The position update expression during the surface flight phase is: in, represents the updated position of the pi-th pelican during the surface flight phase, Represents the objective function value of the pi-th pelican at this position during the surface flight phase.

6. The rolling bearing fault diagnosis method according to claim 5, characterized in that: In step S400, using correlation entropy combined with kurtosis to evaluate multiple intrinsic mode functions to select a suitable intrinsic mode function for reconstruction includes the following steps: S410, for the decomposed intrinsic mode function u k (t), calculate its absolute value component y(t) = |u k (t)|, use the kernel density estimation method to solve the probability density function F(y) of y(t), the calculation formula is: Where d is the bandwidth parameter in the kernel density estimation; S420, using the range of y(t) where the probability density function is greater than the critical probability φ, 0≤φ≤1, the minimum value expression of y(t) within the range is: Among them, y low (t) represents the minimum value in the range, and f(s) represents u k The absolute value component of (t); S430, using a linear interpolation algorithm to low The points are replaced to obtain new modal components S440, according to the new modal components Calculate the kurtosis value using the following formula: Among them, K adj represents the kurtosis value calculated based on the newly obtained modal components; S450, using the kurtosis value in combination with the correlation entropy to calculate an improved kurtosis index, the expression includes: Among them, MCC(e i ) represents the maximum correlation entropy, CEIK(i) is the improved kurtosis index, e i =x i -y i represents the error signal, x i ,y i represents two different modal components, κ σ is the Gaussian kernel function, MCC'(i) represents the transposed matrix of the maximum correlation entropy; S460, taking a preset percentage of the maximum value component of the improved kurtosis index as a threshold, selecting an eigenmode function greater than the threshold from among a plurality of eigenmode functions for reconstruction to obtain a reconstructed signal.

7. The rolling bearing fault diagnosis method according to claim 6, characterized in that: The default percentage is 70%.

8. The rolling bearing fault diagnosis method according to claim 5, characterized in that: The reconstructed signal is subjected to Hilbert transform analysis to obtain the envelope spectrum. The final fault category is diagnosed based on whether the peak spectrum value on the envelope spectrum is consistent with the theoretical frequency of different fault characteristics of rolling bearings.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the steps of the rolling bearing fault diagnosis method as described in any one of claims 1 to 8 are implemented.

10. A computer device comprising a memory, a processor and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the program, the steps of the rolling bearing fault diagnosis method according to any one of claims 1 to 8 are implemented.

Images