Method for predicting residual service life of rotary machine

Abstract

The invention belongs to the technical field of residual life prediction, and discloses a method for predicting the residual service life of a rotating machine. According to the method, linear multi-fraction Levy stable motion is adopted to describe long-range dependence, non-Gaussian and local and global irregularities in the rotating mechanical property degradation process. The BGTLLM model uses a BGRU drift network as a drift component in the LMSM model. And the prior selection problem of the degradation path in the traditional random process theory is eliminated. By adopting the similarity transfer learning strategy based on the DTW-Wasserstein distance, historical data can be fully applied, and the RUL prediction performance under the current part of information is improved. According to the BGTLLM adaptive hybrid RUL prediction model provided by the method, the neural network and the random process are combined and have complementary advantages.

Description

Technical Field

[0001] The present invention belongs to the technical field of remaining service life prediction, and specifically relates to a method for predicting the remaining service life of a rotating machinery based on the construction of a multi-domain dynamic entropy weighted minimum quantization error health index, a combination of multi-fractional-order Lévy stable motion and a BGRU neural network. Background Art

[0002] Due to the influence of internal and external environments, the performance of rotating machinery will gradually degrade with the increase of service time and eventually fail. Its degradation or failure will directly affect the performance and reliability of the equipment. By accurately capturing the performance changes of the equipment and estimating the remaining life, a reasonable maintenance strategy can be formulated in a targeted manner, which can greatly reduce the losses caused by failure and improve the operational reliability of rotating machinery.

[0003] Patent specification CN 111159846 A discloses a method for predicting the remaining useful life of a rolling bearing. The method is a method for predicting the remaining effective life of a rolling bearing based on fractional Lévy stable motion. The method uses a variational mode decomposition method to obtain the intrinsic mode function, and selects the optimal mode through monotonicity, robustness and trend indicators. By setting a failure threshold, the degradation process of the bearing is modeled based on the fractional Lévy stable motion, and the remaining effective life of the bearing is predicted using the Monte Carlo method. The shortcomings of this method are that, first, the fractional Lévy stable motion (FLSM) uses a fixed Hurst index, which reflects the global irregularity of the degradation process. Due to the influence of changes in the external environment and measurement noise, local irregularities may occur in the degradation process of the equipment. The degradation model constructed by FLSM usually requires the degradation path of the equipment to be determined in advance based on experience or knowledge. However, due to the individuality of the degradation process, the degradation process is often unknowable. In addition, this method and similar random process methods can often only use the degradation information of the current data for prediction, and lack the effective use of the same type of historical data. Therefore, it is necessary to propose a new prediction method for the remaining useful life of rotating machinery. Summary of the invention

[0004] The purpose of the present invention is to address the technical drawbacks of the existing methods for predicting the remaining service life of rotating machinery, and to propose a method for predicting the remaining service life of rotating machinery based on random processes and neural networks.

[0005] The technical solution of the present invention:

[0006] A method for predicting the remaining service life of a rotating machinery comprises the following steps:

[0007] Step 1: Use an acceleration vibration sensor to collect high-frequency vibration signals of rotating machinery, obtain historical data of the entire life cycle of similar rotating machinery and historical data of currently predicted rotating machinery, and select and extract characteristic indicators reflecting the performance degradation of rotating machinery based on time domain, frequency domain, time-frequency domain, entropy domain and chaos domain indicators of vibration signals;

[0008] Step 2: Using monotonicity, trend, and robustness indicators, a comprehensive evaluation method that reflects the quality of characteristic indicators is established. Multiple characteristic indicators extracted in step 1 are scored and evaluated, and the top 15 characteristic indicators with the highest scores are selected as preferred indicators. Based on the dynamic entropy weight method, the preferred indicators are weighted to characterize their contribution to degradation and obtain contribution weights. Based on the SOM neural network, the preferred indicators and the contribution weights, the preferred indicators are fused by dimensionality reduction to obtain the entropy weighted minimum quantization error (DEWMQE) as a preliminary health indicator, which is calculated as shown in the following formula:

[0009]

[0010] in, represents the kth optimal index of the Best Matching Unit (BMU), represents the value of the kth optimal index at time m, w k is the contribution weight of the kth preferred indicator, and H represents the number of preferred indicators;

[0011] Step 3: A windowed inertial smoothing mechanism is used to eliminate local trends and outliers in the entropy weighted minimum quantization error (DEWMQE), and a unified benchmark failure threshold is defined based on a scale correction mechanism;

[0012] First, the entropy-weighted minimum quantization error is defined as the initial health indicator (IHI); next, the average value of IHI during the stable operation phase is calculated to obtain the relative value of IHI RHI; secondly, a sliding window of length n is defined, and the least squares-based linear regression method is applied to fit within the sliding window; let WRHI(j) = RHI(j)-a·x j -b to eliminate local trends, define {WRHI(j)} max =μ+3σ and {WRHI(j)} min =μ-3σ to eliminate outliers; the construction of HI based on window inertial smoothing is described by the following formula:

[0013]

[0014] in, s represents the number of IHI in the healthy stage, assuming s = 100; y = a·x i +b by solving Obtained, a, b are the slope and intercept in the linear regression fitting model, x i represents the time series in the sliding window, y refers to the collective name of the values ​​obtained by linear regression fitting on the sliding window, y i refers to the specific i-th fitting value; WRHI(j) represents the IHI after eliminating the local trend, i represents the i-th point in HI, and j represents the j-th point in the sliding window; and are the mean and standard deviation of WRHI respectively; {WRHI(j)} min and {WRHI(j)} max Respectively represent the maximum and minimum limits for filtering outliers;

[0015] A scale correction mechanism is used to normalize different HIs in the same domain; assuming that the normalized fault threshold is expressed as 1, the correction parameter CP is determined by the following formula:

[0016] CP=1 / ((T max +T min ) / twenty three)

[0017] Among them, T max represents the maximum failure threshold in the health indicator, T min Represents the lowest fault threshold in the health indicator; after obtaining the correction parameter CP, all health indicators are corrected to obtain the final health indicator

[0018] Health indicators and achieve a unified fault threshold; as shown in the following formula:

[0019] HI=CP*SHI(4)

[0020] HI represents the final health index obtained after the scale correction mechanism;

[0021] Step 4: Use the multi-fractional Lévy stable motion (LMSM) Y(u,H(u)) and the BGRU drift neural network Ψ(u; θ Ψ ) is used to model the degradation process of rotating machinery; wherein, the degradation model based on LMSM is modeled as shown in the following formula:

[0022] Z(u)=Z(0)+Ψ(u;θ Ψ )+λY(u,H(u))(5)

[0023] In the formula, Ψ(u; θ Ψ ) represents the BGRU drift neural network, θ Ψrepresents the parameters in the drift neural network Ψ(.), Y(u,H(u)) represents the LMSM, λ is the diffusion coefficient, Z(0) is the initial value in the degradation process, and Z(0) = 0;

[0024] The incremental differential iteration of the degradation process is shown in the following formula:

[0025]

[0026] Among them, S α represents the Lévy stable distribution; H(u) represents the Holder index; α is the stability index in the characteristic function of the Lévy stable distribution, which represents the degree of influence of the local pulse on the distribution; Δu is the time increment; the characteristic function of the Lévy stable distribution is defined as:

[0027]

[0028] Among them, β is the skewness parameter, which is used to describe the skewness of the distribution, -1≤β≤1,β≠0; the position of the distribution is controlled by the parameter γ, The parameter δ ≥ 0 indicates the degree of deviation of the random sequence from the normal level; κ is a real number; the Lévy stable distribution uses α and β to describe its shape, and the linear transformation is controlled by μ and δ;

[0029] Step 5: Based on multiple health indicators of similar rotating machinery, the BGRU drift neural network is trained and used as the drift term of the LMSM degradation model;

[0030] Step 6: Use a DTW-Wasserstein comprehensive distance measurement method to find the historical degradation data closest to the current degradation state of the rotating machinery, and fine-tune the parameters of the current BGRU drift neural network using a transfer learning method based on distance similarity;

[0031] When fine-tuning the BGRU drift neural network, the shared GRU layer parameters remain unchanged, while the FCN layer parameters are retrained to adapt to the characteristics of the transferable samples;

[0032] Among them, the calculation of the DTW-Wasserstein comprehensive distance measurement method is shown in the following formula:

[0033]

[0034] in, represents the DTW distance between the target health indicator and the jth historical health indicator, represents the Wasserstein distance between the target health indicator and the jth historical health indicator, ρ1 and ρ2 correspond to the distance weights between the two; D jThe minimum value D corresponds to the historical health indicator that is most similar to the target health indicator;

[0035] Step 7: Substitute the updated GRU drift neural network into the LMSM-based degradation model in step (4), and use the parameter estimation method to estimate the diffusion coefficient and other parameters in the LMSM-based degradation model, and use the local stationary method to estimate the Holder index;

[0036] The process of parameter estimation is as follows:

[0037] Assume that the sampling time of the degradation data is {u0,u1,u2,…,u k}, the sampling interval is T; the degenerate sampling data is Then the degenerate difference vector is: where the pth element is defined as z p According to formula 7,

[0038] In the parameter estimation process, α, β, μ, and δ are the stability index, skewness parameter, location parameter, and scale parameter, respectively. are the estimated stability index, location parameter, scale parameter and drift coefficient respectively.

[0039] Based on the characteristic function of Levy motion, let κ = 1; we can get the estimated value of δ

[0040]

[0041] Define κ = κ0, and estimate α as follows:

[0042]

[0043] in,

[0044] Using the difference sequence z 1:k The complex field estimates the location parameter γ of the Levy motion:

[0045]

[0046] Since the mathematical relationship between the scaling function δ and the diffusion coefficient λ is The estimated value of the diffusion coefficient λ is therefore:

[0047]

[0048] The Holder index is estimated as follows:

[0049]

[0050] in, k is the largest integer within S / n. S is the number of data points, Y(j,H u ) represents the local sequence of the degradation process.

[0051] Step 8: Based on the first arrival time, failure threshold and BGTLLM model, the remaining life of the equipment is predicted using Monte Carlo simulation (MC);

[0052] Specifically, the first arrival time when the equipment performance degrades to the failure threshold ω is defined as:

[0053]

[0054] The steps to predict the remaining life of equipment using Monte Carlo simulation (MC) are as follows:

[0055] 1) Define the number of MC simulation experiments as M and set the starting position of the current RUL prediction as u f , initialize RUL value l=1;

[0056] 2) Sampling random variables

[0057] 3) Determine whether the random variable R is greater than [Ψ(u+Δu; θ Ψ )-Ψ(u;θ Ψ )], if yes, the current l is considered to be the RUL value of the current starting point; if it is less than [Ψ(u+Δu; θ Ψ )-Ψ(u;θ Ψ )], then let l = l + 1, and return to step 2); after M simulations, count the frequency corresponding to each l to form a vector l 1:M , whose PDF is defined as P(l 1:M ), so at the starting position u f The RUL prediction value of the BGTLLM model can be calculated as E(l 1:M );

[0058] 4) Plot the probability density of RUL.

[0059] Beneficial effects of the present invention: The multi-domain dynamic entropy weighted minimum quantization error health index proposed in the present invention overcomes the shortcomings of single-domain feature extraction, extracts features of reaction equipment degradation in multi-domain and multi-dimensional, and weighs the sensitivity and advantages of different features, and realizes effective dimension reduction and fusion of features based on the SOM algorithm; the window inertia smoothing mechanism is used to remove the local trend and outliers of the initial HI, and the scale correction mechanism is used to construct a unified failure threshold. The proposed BGTLLM adaptive hybrid RUL prediction model combines neural networks and random processes to complement each other. Specifically, the linear multi-fractional Lévy stable motion (LMSM) is used to describe the long-range dependence, non-Gaussian, local and global irregularities of the degradation process of rotating machinery performance. The BGTLLM model uses the GRU drift network as the drift component in the LMSM model. The problem of prior selection of degradation paths in traditional random process theory is eliminated. The similarity transfer learning strategy based on DTW-Wasserstein distance can fully apply historical data and improve the RUL prediction performance under the current partial information. BRIEF DESCRIPTION OF THE DRAWINGS

[0060] Figure 1 It is the overall framework diagram for RUL prediction of rotating machinery based on new health indicators and BGTLLM method.

[0061] Figure 2 This is the basic principle diagram of GRU-TL for the trend term of random process.

[0062] Figure 3 It is a diagram of health index and failure threshold, where (a) is before scale correction and (b) is after scale correction. Figure 4 It is the RUL prediction result diagram of the bearings, where (a) is the RUL prediction of bearing 1-4, and (b) is the RUL prediction of bearing 1-7. DETAILED DESCRIPTION

[0063] The specific implementation of the present invention is further described below in conjunction with the accompanying drawings and technical solutions.

[0064] A method for predicting the remaining service life of a rotating machinery comprises the following steps:

[0065] (1) The vibration signal of the rotating machinery is collected in real time by using sensors, and multi-domain indicators are extracted from the vibration signal to characterize the degradation state of the system. The selected time domain and frequency domain indicators are shown in Table 1. Among them, for the bearing, the four fault characteristic frequencies of the bearing are selected as features. For the gear, first, TSA is used to eliminate the random noise generated during the rotation of the gear. The GMF component and its harmonic components unrelated to the fault are removed from the denoised signal to obtain the residual signal, from which the quasi-normalized kurtosis (NA4) is extracted. On this basis, the sideband frequency components around the GMF at a certain distance from the rotation frequency of the shaft are further removed to obtain the differential signal, and then the kurtosis (FM4), sixth-order moment (M6A) and eighth-order moment (M8A) of the signal are extracted.

[0066] In the time-frequency domain, wavelet packet decomposition is selected to extract degradation features. Morlet is selected as the mother wavelet to extract the impact characteristics inside the signal, with 4 decomposition layers and the first 8 normalized energies as the time-frequency domain feature parameters. In the entropy domain, sample entropy and distribution entropy of the extracted signal are selected. In the chaos domain, Lyapunov exponent spectrum is selected to extract the chaotic features of the vibration signal.

[0067] Table 1 Time domain and frequency domain characteristics

[0068]

[0069]

[0070] (2)) Using monotonicity, trend, and robustness indicators, a comprehensive evaluation method that reflects the quality of characteristic indicators is established, and multiple characteristic indicators extracted in step 1 are scored and evaluated. The top 15 characteristic indicators with the highest scores are taken as preferred indicators. Next, the dynamic entropy weight method is used to evaluate the weights of the preferred indicators to characterize their contribution to degradation and obtain the contribution weights. Next, based on the SOM neural network, the preferred index and the contribution weights are used to reduce the dimension of the above 15 preferred indicators by calculating the distance between the characteristic matrix and the BMU, and the dynamic entropy weighted minimum quantization error (DEWMQE) is obtained as the initial health indicator, which is calculated as shown in the following formula:

[0071]

[0072] in, represents the kth weight of BMU, represents the value of the kth preferred feature at time m, w k is the contribution weight of the kth feature. H represents the number of preferred features.

[0073] (3) A windowed inertial smoothing mechanism is used to eliminate local trends and outliers in DEWMQE. Furthermore, a scale correction mechanism is used to unify the fault threshold.

[0074] Specifically, first, DEWMQE is defined as the initial HI (IHI). Next, the average value of IHI during the stable operation stage is calculated to obtain the relative value of HI (RHI). This method helps to initially alleviate the individual differences between multiple HIs. Secondly, on the obtained RHI, a sliding window of length k is defined, and a least squares-based linear regression method is applied within the window for fitting. Let WRHI(j) = RHI(j)-a·x j -b to eliminate local trends. Define {WRHI(j)} max =μ+3σ and {WRHI(j)} min =μ-3σ to eliminate outliers. The construction of HI based on window inertial smoothing is described by the following formula:

[0075]

[0076] in, y=a·x i +b,s represents the number of IHI in the healthy stage. Here, s = 100. y = a·x i +b, by solving get. and are the mean and standard deviation of WRHI, respectively.

[0077] Even between devices of the same type and under the same conditions, HI can vary significantly due to variables such as individual differences, environmental factors, and measurement errors. Therefore, a scale correction mechanism is used to normalize different HIs within the same domain. Assuming that the normalized fault threshold is expressed as 1, the correction parameter CP is determined by the following formula:

[0078] CP=1 / ((T max +T min ) / twenty four)

[0079] Among them, T max represents the maximum fault threshold in the health indicator sample, T min Represents the lowest fault threshold in the health indicator sample. After obtaining the correction parameter CP, all health indicators are corrected to obtain the final health indicator and achieve a unified fault threshold. As shown in the following formula:

[0080] HI=CP*HI′(5)

[0081] Among them, HI is the final health index, and HI′ represents the health index when the scale correction mechanism is not performed.

[0082] (4) The degradation process of rotating machinery is modeled using multi-fractional Lévy stable motion (LMSM) and GRU drift neural network.

[0083] The degradation sequence of rotating machinery {Z(u),u≥0} represents a one-dimensional random process, which describes the degradation of rotating machinery as a time-dependent function. If the degradation process of rotating machinery has long-range dependence, incremental non-Gaussian and multi-fractal characteristics, the degradation process based on LMSM and BGRU drift neural network can be defined as follows:

[0084] Z(u)=Z(0)+Ψ(u;θ Ψ )+λY(u,H(u))(6)

[0085] Among them, Ψ(u; θ Ψ ) represents the BGRU drift neural network, θ Ψ represents the parameters in the network, Y(u,H(u)) represents the LMSM. λ is the drift coefficient. Z(0) is the initial value in the degradation process. In order not to lose generality, Z(0)=0.

[0086] The incremental differential iteration of the degradation process is shown in the following formula:

[0087]

[0088] Among them, S α represents the Lévy stable distribution. H(u) represents the Holder index, and α is the stability index in the characteristic function of the Levy stable distribution, which represents the degree of influence of the local pulse on the distribution. The characteristic function of the Lévy stable distribution is defined as:

[0089]

[0090] Among them, H (-1≤β≤1,β≠0) is the skewness parameter, which is used to describe the skewness of the distribution. The location of the distribution is determined by the parameter To control, and the parameter δ ≥ 0 indicates the degree of deviation of the random sequence from the normal level. The shape of the Lévy stable distribution is usually described by α and β, while the linear transformation is controlled by γ and δ.

[0091] (5) Based on the historical health indicators of the same type of rotation, the GRU drift neural network is trained and used as the drift term of the LMSM degradation model;

[0092] (6) Next, a DTW-Wasserstein comprehensive distance measurement method is used to find the historical health indicator that is closest to the health indicator of the current device from multiple historical health indicators of similar devices. A distance similarity-based transfer learning method is used to fine-tune the parameters of the current BGRU drift network. Specifically, during fine-tuning, the shared BGRU layer parameters remain unchanged, while the fully connected (FCN) layer parameters are retrained to adapt to the characteristics of the transferable samples. Figure 2 shown.

[0093] Among them, the calculation of the DTW-Wasserstein comprehensive distance measurement method is shown in the following formula:

[0094]

[0095] in, represents the DTW distance between the target health indicator and the jth historical health indicator; represents the Wasserstein distance between the target health indicator and the jth historical health indicator, ρ1 and ρ2 are distance weights. j The minimum value D corresponds to the historical health indicator that is most similar to the target health indicator.

[0096] (7) Substitute the updated BGRU drift network into the degradation model based on LMSM in step (4), and use the parameter estimation method to estimate the diffusion coefficient and other parameters in LMSM, and use the local stability method to estimate the Holder index.

[0097] The process of parameter estimation is as follows:

[0098] Assume that the sampling time of the degradation data is {u0,u1,u2,…,u k}. The sampling interval is T. The degraded sampling data is Then the degenerate difference vector is: where the qth element is defined as z q According to formula 7,

[0099] In the parameter estimation process, α, β, γ, and δ are the stability index, skewness parameter, location parameter, and scale parameter, respectively. are the estimated stability index, location parameter, scale parameter and drift coefficient respectively.

[0100] Based on the characteristic function of Levy motion, let κ = 1; we can get the estimated value of δ

[0101]

[0102] Define κ = κ0, and estimate α as follows:

[0103]

[0104] in,

[0105] Using the difference sequence z 1:k The complex field estimates the location parameter γ of the Levy motion:

[0106]

[0107] Since the mathematical relationship between the scaling function δ and the diffusion coefficient λ is The estimated value of the diffusion coefficient λ is therefore:

[0108]

[0109] The Holder index is estimated as follows:

[0110]

[0111] in, k is the largest integer within S / n. S is the number of data points, Y(j,H u ) represents the local sequence of the degradation process.

[0112] (8) Based on the first arrival time, failure threshold and BGTLLM model, the remaining life of the equipment is predicted using Monte Carlo simulation (MC).

[0113] Specifically, the first arrival time when the equipment performance degrades to the failure threshold ω is defined as:

[0114]

[0115] The steps to predict the remaining life of equipment using Monte Carlo simulation (MC) are as follows:

[0116] 1) Define the number of MC simulation experiments as M and set the starting position of the current RUL prediction as u f , initialize RUL value l=1;

[0117] 2) Sample random variable R,

[0118] 3) Determine whether the random variable R is greater than [Ψ(u+Δu; θ Ψ )-Ψ(u;θ Ψ )], if yes, then the current l is considered to be the RUL value of the current starting point. If it is less than [Ψ(u+Δu;θ Ψ )-Ψ(u;θ Ψ)], then let l = l + 1, and return to step 2). After the M simulations are completed, count the frequency corresponding to each l to form a vector l 1:M , its PDF is defined as P(l 1:M ), therefore, at the starting position u f The RUL prediction value of the BGTLLM model can be calculated as E(l 1:M );

[0119] 4) Plot the probability density of RUL.

[0120] Implementation Cases:

[0121] This embodiment verifies the method proposed in the present invention with the rolling bearing accelerated life test data on the rolling bearing test platform. The test bearing is a rolling bearing, and the experimental data of working condition 1 is used. The speed is set to 1800r / min and the radial force is 4000N. The vibration signal of the rolling bearing is collected using an acceleration sensor, with a sampling frequency of 25.6kHz, a sampling interval of 10s, and a sampling time of 0.1s. The implementation framework of the BGTLLM hybrid model for RUL prediction is as follows Figure 1 shown.

[0122] The HI and failure thresholds constructed by multi-domain feature extraction, dynamic entropy weight method, SOM dimensionality reduction fusion, window inertial smoothing and scale correction are as follows: Figure 3 shown.

[0123] Bearings 1-4 and 1-7 are used as test bearings. Their predicted remaining life is as follows Figure 4 shown.

[0124] Practical application effect: The present invention can be applied to the fault prediction and health management of rotating machinery such as gears and bearings, and realize the online and offline prediction of the remaining service life of the equipment. It provides theoretical basis and technical support for the intelligent operation and maintenance decision-making of mechanical systems. The technical application of this invention can improve the safety, reliability and availability of enterprise mechanical equipment and reduce the maintenance cost of equipment.

Claims

1. A method for predicting the remaining useful life of a rotating machine, characterized in that: The following steps are involved: Step 1: Use an acceleration vibration sensor to collect high-frequency vibration signals of rotating machinery, obtain historical data of the entire life cycle of similar rotating machinery and historical data of currently predicted rotating machinery, and select and extract characteristic indicators reflecting the performance degradation of rotating machinery based on time domain, frequency domain, time-frequency domain, entropy domain and chaos domain indicators of vibration signals; Step 2: Using monotonicity, trend, and robustness indicators, a comprehensive evaluation method that reflects the quality of characteristic indicators is established. The multiple characteristic indicators extracted in step 1 are scored and evaluated, and the top 15 characteristic indicators with the highest scores are selected as preferred indicators. Based on the dynamic entropy weight method, the preferred indicators are weighted to characterize their contribution to degradation, and the contribution weights are obtained. Based on the SOM neural network, the optimal indicators and contribution weights, the optimal indicators are fused with dimensionality reduction to obtain the minimum quantization error of entropy weight as the preliminary health indicator; Step 3: Use the window inertial smoothing mechanism to eliminate local trends and outliers in the entropy weighted minimum quantization error, and define a unified benchmark failure threshold based on the scale correction mechanism; Step 4: Use the multi-fractional Lévy stable motion (LMSM) Y(u,H(u)) and the BGRU drift neural network Ψ(u; θ Ψ ) Modeling the degradation process of rotating machinery; Step 5: Based on multiple health indicators of similar rotating machinery, the BGRU drift neural network is trained and used as the drift term of the LMSM degradation model; Step 6: Use a DTW-Wasserstein comprehensive distance measurement method to find the historical degradation data closest to the current degradation state of the rotating machinery, and fine-tune the parameters of the current BGRU drift neural network using a transfer learning method based on distance similarity; When fine-tuning the BGRU drift neural network, the shared BGRU layer parameters remain unchanged, while the fully connected FCN layer parameters are retrained to adapt to the characteristics of the transferable samples; Step 7: Substitute the updated GRU drift neural network into the LMSM-based degradation model in step (4), and use the parameter estimation method to estimate the diffusion coefficient and other parameters in the LMSM-based degradation model, and use the local stationary method to estimate the Holder index; Step 8: Based on the first arrival time, failure threshold and BGTLLM model, the remaining life of the equipment is predicted using Monte Carlo simulation (MC).

2. The method for predicting the remaining useful life of a rotating machine according to claim 1, characterized in that: In step 2, based on the SOM neural network, the preferred indicators and the contribution weights, the preferred indicators are fused with dimensionality reduction to obtain the minimum quantization error of the entropy weight as the preliminary health indicator, which is calculated as shown in the following formula: in, represents the kth optimal index of the Best Matching Unit (BMU), represents the value of the kth optimal index at time m, w k is the contribution weight of the kth preferred indicator, and H represents the number of preferred indicators.

3. The method for predicting the remaining useful life of a rotating machine according to claim 1, characterized in that: The specific steps of step 3 are as follows: First, the entropy-weighted minimum quantization error is defined as the initial health index, i.e., IHI. Next, the average value of IHI during the stable operation phase is calculated to obtain the relative value RHI of IHI. Secondly, a sliding window of length n is defined, and a linear regression method based on least squares is applied to fit the sliding window; let WRHI(j) = RHI(j)-a·x j -b to eliminate local trends, define {WRHI(j)} max =μ+3σ and {WRHI(j)} min =μ-3σ to eliminate outliers; the construction of HI based on window inertial smoothing is described by the following formula: in, s represents the number of IHI in the healthy stage, assuming s = 100; y = a·x i +b by solving Obtained, a, b are the slope and intercept in the linear regression fitting model, x i represents the time series in the sliding window, y refers to the collective name of the values ​​obtained by linear regression fitting on the sliding window, y i refers to the specific i-th fitting value; WRHI(j) represents the IHI after eliminating the local trend, i represents the i-th point in HI, and j represents the j-th point in the sliding window; and are the mean and standard deviation of WRHI respectively; {WRHI(j)} min and {WRHI(j)} max Respectively represent the maximum and minimum limits for filtering outliers; A scale correction mechanism is used to normalize different HIs in the same domain; assuming that the normalized fault threshold is expressed as 1, the correction parameter CP is determined by the following formula: CP=1 / ((T max +T min ) / 2) (3) Among them, T max represents the maximum failure threshold in the health indicator, T min Represents the lowest fault threshold in the health index; after obtaining the correction parameter CP, all health indicators are corrected to obtain the final health index and achieve a unified fault threshold; as shown in the following formula: HI=CP*SHI (4) HI represents the final health index obtained after the scale correction mechanism.

4. The method for predicting the remaining useful life of a rotating machine according to claim 1, characterized in that: In step 4, the degradation model based on LMSM is modeled as shown in the following formula: Z(u)=Z(0)+Ψ(u;θ Ψ )+λY(u,H(u)) (5) In the formula, Ψ(u; θ Ψ ) represents the BGRU drift neural network, θ Ψ represents the parameters in the drift neural network Ψ(.), Y(u,H(u)) represents the LMSM, λ is the diffusion coefficient, Z(0) is the initial value in the degradation process, and Z(0) = 0; The incremental differential iteration of the degradation process is shown in the following formula: Among them, S α represents the Lévy stable distribution; H(u) represents the Holder index; α is the stability index in the characteristic function of the Lévy stable distribution, which represents the degree of influence of the local pulse on the distribution; Δu is the time increment; the characteristic function of the Lévy stable distribution is defined as: Among them, β is the skewness parameter, which is used to describe the skewness of the distribution, -1≤β≤1,β≠0; the position of the distribution is controlled by the parameter γ, The parameter δ ≥ 0 indicates the degree of deviation of the random sequence from the normal level; κ is a real number; the Lévy stable distribution uses α and β to describe its shape, and the linear transformation is controlled by μ and δ.

5. The method for predicting the remaining useful life of a rotating machine according to claim 1, characterized in that: In step 6, the calculation of the DTW-Wasserstein comprehensive distance metric method is shown in the following formula: D=min(D j ) (9) in, represents the DTW distance between the target health indicator and the jth historical health indicator, represents the Wasserstein distance between the target health indicator and the jth historical health indicator, ρ1 and ρ2 correspond to the distance weights between the two; D j The minimum value D corresponds to the historical health indicator that is most similar to the target health indicator.

6. The method for predicting the remaining useful life of a rotating machine according to claim 1, characterized in that: In step 7, the parameter estimation process is as follows: Assume that the sampling time of the degradation data is {u0,u1,u2,…,u k }, the sampling interval is T; the degenerate sampling data is Then the degenerate difference vector is: where the pth element is defined as z p According to formula 7, In the parameter estimation process, α, β, μ, and δ are stability index, skewness parameter, location parameter, and scale parameter, respectively; They are the estimated stability index, location parameter, scale parameter and drift coefficient respectively; based on the characteristic function of Levy motion, let κ = 1; we can get the estimated value of δ Define κ = κ0, and estimate α as follows: in, Using the difference sequence z 1:k The complex field estimates the location parameter γ of the Levy motion: Since the mathematical relationship between the scaling function δ and the diffusion coefficient λ is The estimated value of the diffusion coefficient λ is therefore: The Holder index is estimated as follows: in, k is the largest integer within S / n; S is the number of data points, Y(j,H u ) represents the local sequence of the degradation process.

7. The method for predicting the remaining useful life of a rotating machine according to claim 1, characterized in that: The specific process of step 7 is as follows: The first arrival time when the equipment performance degrades to the failure threshold ω is defined as: The steps to predict the remaining life of equipment using Monte Carlo simulation (MC) are as follows: 1) Define the number of MC simulation experiments as M and set the starting position of the current RUL prediction as u f , initialize RUL value l=1; 2) Sample random variable R, 3) Determine whether the random variable R is greater than [Ψ(u+Δu; θ Ψ )-Ψ(u;θ Ψ )], if yes, the current l is considered to be the RUL value of the current starting point; if it is less than [Ψ(u+Δu; θ Ψ )-Ψ(u;θ Ψ )], then let l = l + 1, and return to step 2); after M simulations, count the frequency corresponding to each l to form a vector l 1:M , whose PDF is defined as P(l 1:M ), so at the starting position u f The RUL prediction value of the BGTLLM model can be calculated as E(l 1:M ); 4) Plot the probability density of RUL.

Images