Bearing life prediction method

Abstract

The invention provides a bearing life prediction method, and belongs to the technical field of equipment life prediction, and the method comprises the steps: obtaining a bearing health state index, and comprehensively calculating a comprehensive health score CHS of a certain time step; based on the comprehensive health score CHS, constructing a prediction model of the bearing health state to the remaining use time; training the constructed prediction model by using a comprehensive health score (CHS); and the trained prediction model is used to identify the dynamic characteristics of the bearing health state to the remaining use time, and the bearing health state is predicted according to the trend of historical health state scoring so as to complete the prediction of the bearing life. According to the method, the problems of high model complexity, insufficient explanation, poor applicability and universality and large over-fitting risk in the prior art are solved, so that high precision, high robustness and high reliability of bearing life prediction are realized.

Description

Technical Field

[0001] The present invention belongs to the technical field of equipment life prediction, and particularly relates to a bearing life prediction method. Background Art

[0002] With the continuous advancement of industrial automation and intelligence, the reliability and maintenance efficiency of equipment have become the core issues of concern for enterprises. As an important component of mechanical equipment, the performance of bearings directly affects the operation stability and safety of the entire system. Therefore, accurately predicting the remaining useful life (RUL) of bearings is crucial for preventing failures, reducing maintenance costs, and improving production efficiency. In recent years, a variety of technologies have been applied in the field of bearing life prediction, including physics-based methods, data-driven methods, and hybrid model methods. However, these technologies still have certain limitations in practical applications and need to be further optimized.

[0003] The physics-based method describes the failure mechanism of bearings by constructing a mathematical model and combines empirical knowledge with the defect growth equation for life prediction. Although this method is effective under specific conditions, the model is usually complex and highly dependent on bearing types and usage environments, so it is difficult to be widely applied. The data-driven method uses historical data and real-time monitoring information and adopts machine learning and statistical techniques to establish a prediction model. However, this method has high requirements for data quality and quantity. Especially in the case of insufficient data or high noise, it is prone to overfitting problems, which will reduce the prediction accuracy. In recent years, with the progress of deep learning technology, recurrent models such as long short-term memory networks (LSTM) and convolutional neural networks (CNN) have been widely used in equipment health state prediction. LSTM performs well in capturing long-term sequence dependence features, but its model architecture is relatively simple, and there are still limitations in applications in complex industrial environments. CNN has advantages in extracting multi-dimensional features, but its overly deep network structure may lead to difficulties in training. To this end, introducing a residual network can effectively alleviate this problem.

[0004] In terms of bearing health state feature extraction, the health state is related to various indicators such as vibration, temperature, current, and acoustics. Although many studies focus on the time-frequency feature extraction of vibration signals, in actual production, different bearings have different dependencies on each feature. Therefore, it is necessary to distinguish the influence degree of each index on the health state through a weighting method. Summary of the Invention

[0005] Aiming at the above deficiencies in the prior art, the present invention provides a bearing life prediction method, which overcomes the problems of high model complexity, insufficient interpretability, poor applicability and generality, and high overfitting risk in the prior art, thereby achieving high-precision, high-robustness, and high-reliability bearing life prediction.

[0006] To achieve the above object, the technical solution adopted by the present invention is as follows: A bearing life prediction method, comprising the following steps:

[0007] S1. Obtain bearing health state indicators and comprehensively calculate the comprehensive health score CHS at a certain time step;

[0008] S2. Based on the comprehensive health score CHS, construct a prediction model for the remaining service life of the bearing health state;

[0009] S3. Use the comprehensive health score CHS to train the constructed prediction model;

[0010] S4. Use the trained prediction model to identify the dynamic characteristics of the bearing health state with respect to the remaining service life, and predict the bearing health state based on the trend of historical health state scores, thereby completing the prediction of the bearing life.

[0011] The beneficial effect of the present invention is that the present invention proposes a bearing life prediction method based on a hybrid data-driven model (CRN-Autoformer) of a convolutional residual network and Autoformer. Specifically, the present invention combines a deep convolutional residual network with the Autoformer model evolved from Transformer, and through multi-dimensional state evaluation and feature extraction, comprehensively analyzes the operating conditions of the bearing. The present invention overcomes the problems of high model complexity, insufficient interpretability, poor applicability and generality, and large overfitting risk in the prior art, thereby realizing high-precision, high-robustness and high-reliability of bearing life prediction.

[0012] Further, in the step S1, specifically:

[0013] Obtain bearing health state indicators including vibration index VI and temperature index TI, and comprehensively calculate the comprehensive health score CHS at a certain time step.

[0014] The beneficial effect of the above further solution is that this solution more comprehensively evaluates the health state of the bearing. The vibration index can reflect the running smoothness of mechanical components, while the temperature index can capture abnormal conditions such as overheating. The comprehensive health score calculated by combining these two indicators can be regarded as a comprehensive and intuitive view of the bearing health state.

[0015] Still further, the expression of the comprehensive health score CHS is as follows:

[0016] CHS = w V ·VI + w T ·TI

[0017] Wherein, w V and w T both represent learnable parameters.

[0018] The beneficial effects of the above further solution are as follows: By introducing learnable parameters and assigning different weights to the indicators, the system can automatically adjust the importance of these two indicators according to the actual application scenario and device characteristics, so as to more accurately reflect the health status of the bearing. These parameters are continuously optimized through training data to better adapt to different types of bearings and working environments. This adaptive ability enables the Comprehensive Health Score (CHS) to more accurately reflect the actual health status.

[0019] Furthermore, the process of obtaining the vibration index VI is as follows:

[0020] A1. Collect vibration signals within the complete life cycle from the bearing sensor, and preprocess the vibration signals to obtain the original vibration signals for extracting time-domain features;

[0021] A2. Use the Hilbert-Huang Transform (HHT) to perform non-linear and non-stationary signal analysis on the original vibration signals, and obtain the vibration index VI based on the Hilbert spectrum generated from the analysis results. Among them, the Hilbert-Huang Transform (HHT) includes Empirical Mode Decomposition (EMD) and Hilbert Spectrum Analysis (HAS). The Empirical Mode Decomposition (EMD) decomposes the original vibration signals into several Intrinsic Mode Functions (IMFs) to capture the inherent features of the vibration signals; the Hilbert Spectrum Analysis (HAS) performs the Hilbert transform on each Intrinsic Mode Function (IMF) and generates the Hilbert spectrum by calculating the instantaneous amplitude and instantaneous frequency to reveal the energy distribution of the signals in time and frequency.

[0022] The beneficial effects of the above further solution are as follows: By introducing the Hilbert-Huang Transform (HHT) to extract the time-frequency domain features of the bearing vibration signals, the present invention can adapt to non-linear and non-stationary signal characteristics. Especially during the fault development process, it can accurately capture the dynamic features of the vibration signals changing with time and frequency, providing high-quality input data for the model.

[0023] Furthermore, the specific process of A2 is as follows:

[0024] B1. Obtain all the maximum and minimum values of the original vibration signal segment, and use the cubic spline interpolation method to generate the upper envelope and the lower envelope respectively;

[0025] B2. Calculate the average value of the upper envelope and the lower envelope;

[0026] B3. Subtract the average value from the original vibration signal to obtain h1(t), where h1(t) represents the vibration signal after normalization processing;

[0027] B4. If h1(t) satisfies the conditions of the Intrinsic Mode Function (IMF), then h1(t) is taken as the first Intrinsic Mode Function (IMF), where the conditions for the Intrinsic Mode Function (IMF) are that the number of extreme points and zero-crossing points are equal or differ by one, and the average of the envelopes of the local maxima and minima is zero. The extreme points include the maxima and minima of the original vibration signal segment;

[0028] B5. Subtract the first Intrinsic Mode Function (IMF) from the original vibration signal to obtain the residual signal r1(t);

[0029] B6. Determine whether the residual signal r1(t) is a monotonic function. If so, use the following formula to represent the original vibration signal as the sum of several Intrinsic Mode Functions (IMFs) and a residual term, and proceed to step B7. Otherwise, return to step B1 until the residual signal r n (t) is a monotonic function:

[0030]

[0031] where x(t) represents the original vibration signal, i represents the index variable, n' represents the upper limit value of the division, and IMF i (t) represents the i-th Intrinsic Mode Function (IMF);

[0032] B7. Perform the Hilbert transform on the i-th Intrinsic Mode Function (IMF) using the following formula to obtain the analytic signal z i (t):

[0033] z i (t) = IMF i (t) + jH[IMF i (t)]

[0034]

[0035] where j represents the imaginary unit, H[] represents the Hilbert transform, P.V. represents the Cauchy principal value function, τ represents the integration variable, t represents the time parameter of the signal, and d represents the differential symbol;

[0036] B8. Based on the analytic signal z i (t), use the following formula to obtain the instantaneous amplitude a i (t) and the instantaneous phase φ i (t):

[0037]

[0038] B9. According to the instantaneous phase φ i (t), use the following formula to obtain the instantaneous frequency ω i (t):

[0039]

[0040] B10. According to the instantaneous amplitude a i (t) and the instantaneous frequency ω i (t), the Hilbert spectrum H(ω, t) is obtained by the following formula. The Hilbert spectrum H(ω, t) represents the distribution of the instantaneous amplitude a i (t) at different frequencies and times:

[0041]

[0042] where δ() represents the Dirac function, ω represents the frequency variable, which is used to represent the distribution of the signal at different frequencies;

[0043] B11. From each Hilbert spectrum H(ω, t), the total energy, frequency center, frequency bandwidth, and maximum energy frequency are obtained by the following formula:

[0044]

[0045] where E i represents the energy of the i-th component, T represents the upper limit of the time variable, d represents the differential symbol, ω represents the frequency variable, CenterFreq i represents the center frequency of the i-th component, Bandwidth i represents the bandwidth of the i-th component, MaxEnergyFreq i represents the maximum energy frequency of the i-th component;

[0046] B12. Using the following formula, the vibration index VI is obtained by fusing the results of B11:

[0047] VI = [E1, CenterFreq1, Bandwidth1, MaxEnergyFreq1,...,

[0048] E n , CenterFreq n , Bandwidth n , MaxEnergyFreq n

[0049] where E n represents the energy of the last component, CenterFreq n represents the center frequency of the i-th component, Bandwidth n represents the bandwidth of the i-th component, MaxEnergyFreq n represents the maximum energy frequency of the i-th component.

[0050] ​The beneficial effects of the above further solution are as follows: The vibration index calculated by integrating the information of multiple frequency bands can more comprehensively reflect the characteristics of the vibration signal, which helps to improve the accuracy and reliability of fault detection.

[0051] Furthermore, the expression of the temperature index TI is as follows:

[0052]

[0053] Wherein, represents the average temperature in the nth time period, represents the standard deviation in the nth time period, represents the maximum temperature in the nth time period, represents the minimum temperature in the nth time period, represents the average temperature in the i'th time period, T i' represents the length of the i'th time period, T(t) represents the temperature at time t, d represents the differential symbol, t represents the time variable, t i' represents the start time of the i'th time period, Δt represents the time interval, represents the standard deviation in the i'th time period, represents the maximum temperature in the i'th time period, represents the minimum temperature in the i'th time period.

[0054] The beneficial effects of the above further solution are as follows: By integrating the information of multiple time periods (including average temperature, standard deviation, maximum temperature, and minimum temperature), this solution can more comprehensively reflect the characteristics of the temperature signal. This helps to improve the accuracy and reliability of fault detection. At the same time, by analyzing each time period separately, abnormal conditions in different time periods can be captured. This refined method enables the system to identify small fluctuations or high-frequency noises that may be masked in the overall signal.

[0055] Furthermore, the prediction model includes:

[0056] An input data preprocessing module, which is used to perform normalization processing on the comprehensive health score CHS and cut it into time series segments of a fixed size;

[0057] A CNN residual network, which is used to extract low-level features based on the time series segments using an initial convolutional layer, and based on the extracted low-level features, use two residual blocks to extract high-level features by increasing the number of channels, and use a global average pooling layer to compress the feature dimension to extract the spatial features of the bearing health state data;

[0058] The Autoformer module is used to separate the trend term and the periodic term based on the spatial features of the bearing health state data, and generate the bearing life prediction result by aggregating the time series features;

[0059] The prediction output module is used to output the bearing life prediction result.

[0060] The beneficial effects of the above further solution are as follows: The present invention uses the convolutional residual network (CRN) as the feature extraction module, and effectively alleviates the problem of gradient disappearance in the training of deep networks through residual connections, significantly improving the extraction ability of the temporal features of the bearing vibration signal and the spatial features of the temperature data, so as to comprehensively characterize the bearing health state. In addition, the introduced Autoformer module combines the time series decomposition mechanism to specifically identify the trend and seasonal features, significantly enhancing the modeling ability of long-distance dependencies in time series data, and enabling the bearing life prediction to maintain high accuracy and stability over a long time span. In summary, the prediction model uses the convolutional network to efficiently extract the multi-dimensional data features of the bearing, and combines the Autoformer to comprehensively integrate the extracted temporal features, so as to achieve accurate prediction of the bearing RUL.

[0061] Furthermore, the comprehensive health score CHS is normalized and cut into time series segments of a fixed size, specifically:

[0062] C1. Normalize the comprehensive health score CHS;

[0063] C2. Use a Gaussian filter to remove the data noise in the normalized comprehensive health score CHS;

[0064] C3. Use the following formula to cut the comprehensive health score CHS processed by C2 into time series segments of a fixed size:

[0065] X i' = x[(i'-1)S:(i'-1)S+L]

[0066] where X i represents the i'-th time series segment, x represents the original time series data, S represents the step size, and L represents the fixed length of each segment.

[0067] The beneficial effects of the above further solution are as follows: Preprocessing the data is beneficial to improving the training quality of the neural network and avoiding gradient disappearance and explosion during the training process. The generalization ability and robustness are improved by maintaining the numerical stability of the neural network.

[0068] Furthermore, based on the features after compressing the feature dimension, the trend term and the periodic term are separated, and the bearing life prediction result is generated by aggregating the time series features, specifically:

[0069] D1. Based on the spatial characteristics of the bearing health state data, perform time series decomposition using the following formula:

[0070] X t = T t + S t + R t

[0071] where X t represents the original time series value at time t, T t represents the trend term, reflecting the long-term change trend of the time series, S t represents the seasonal term, reflecting the periodic change of the time series, and R t represents the random term, reflecting the random fluctuation of the time series;

[0072] D2. Based on the time series decomposition result, perform trend extraction analysis using the following formula:

[0073]

[0074] where k represents the moving average window size, and X t+i represents the original time series value at time t + i';

[0075] D3. Based on the trend extraction analysis result, use the autocorrelation mechanism to capture the periodicity and autocorrelation in the time series;

[0076] D4. Based on the autocorrelation, construct an attention mechanism to focus on the bearing health state time points with strong correlation;

[0077] D5. Based on the attention result, perform aggregation of the time series features using the following formula:

[0078]

[0079] where F t represents the aggregated time series features, S represents the number of time scales considered, w s represents the weights of each time scale, X t-s:t represents the time window data from t - s to t, and θ s represents the parameter of the one-dimensional convolution;

[0080] D6. Based on the aggregated time series features, generate the bearing life prediction result using the following formula:

[0081]

[0082] where The bearing life prediction result at time t is denoted as, φ represents the activation function of the output layer, both W1 and W2 represent weight matrices, both b1 and b2 represent bias terms, and σ represents the activation function of the intermediate layer.

[0083] The beneficial effects of the above further solution are as follows: Autoformer learns the time operating state of the bearing by identifying the periodic patterns of the time series. The autocorrelation coefficient explains the patterns of the time series, and these patterns may be key indicators for early fault warning. At the same time, by utilizing the autocorrelation information, a more accurate time series prediction model can be constructed. High autocorrelation means that future values largely depend on past values, which helps to improve the prediction accuracy. When the actually observed autocorrelation structure does not match the expectation, it may indicate the occurrence of abnormal conditions or impending failures, thus allowing preventive measures to be taken and improving the accuracy of early warning. Brief Description of the Drawings

[0084] Figure 1 It is the flowchart of the method of the present invention.

[0085] Figure 2 It is the block diagram of the bearing life prediction of the present invention.

[0086] Figure 3 It is the schematic diagram of the Autoformer module architecture of the present invention. Detailed Embodiments

[0087] The following describes the detailed embodiments of the present invention to facilitate those skilled in the art of the present technology to understand the present invention. However, it should be clear that the present invention is not limited to the scope of the detailed embodiments. For those of ordinary skill in the art of the present technology, as long as various changes are within the spirit and scope of the present invention defined and determined by the appended claims, these changes are obvious, and all inventions and creations using the concept of the present invention are within the scope of protection.

[0088] Embodiment

[0089] Before describing the present invention, the following terms are explained:

[0090] Autoformer (Adaptive Transformer): An innovative time series prediction model based on Transformer proposed by the research team of Alibaba Cloud.

[0091] The present invention is based on a hybrid data-driven model (CRN-Autoformer) of a Convolutional Residual Network and Autoformer, and evaluates the operating state of a bearing and predicts its remaining service life through sensor data of bearing health state indicators. Aiming at the problem that the life data of bearings in the actual use environment is easy to obtain, while the reference life provided by the manufacturer does not fully consider environmental differences, the present invention introduces the idea of transfer learning and fine-tunes the weights and parameters of the pre-trained network model to learn the time-series characteristics of the health state of bearings at different time nodes.

[0092] As Figure 1 and Figure 2 shown, the present invention provides a bearing life prediction method, and its implementation method is as follows:

[0093] S1. Obtain bearing health state indicators and comprehensively calculate the comprehensive health score CHS at a certain time step, which is specifically:

[0094] Obtain bearing health state indicators including the vibration index VI and the temperature index TI, and comprehensively calculate the comprehensive health score CHS at a certain time step;

[0095] In this embodiment, the process of obtaining the vibration index VI is as follows:

[0096] A1. Collect vibration signals within the complete life cycle from a bearing sensor, and preprocess the vibration signals to obtain the original vibration signals for extracting time-domain features;

[0097] A2. Use the Hilbert-Huang transform HHT to perform non-linear and non-stationary signal analysis on the original vibration signals, and obtain the vibration index VI based on the Hilbert spectrum generated from the analysis results. Among them, the Hilbert-Huang transform HHT includes empirical mode decomposition EMD and Hilbert spectrum analysis HAS. The empirical mode decomposition EMD decomposes the original vibration signals into several intrinsic mode functions IMF to capture the intrinsic characteristics of the vibration signals; the Hilbert spectrum analysis HAS performs the Hilbert transform on each intrinsic mode function IMF and generates the Hilbert spectrum by calculating the instantaneous amplitude and instantaneous frequency to reveal the energy distribution of the signals in time and frequency, which is specifically:

[0098] B1. Obtain all the maximum and minimum values of the original vibration signal segment, and use the cubic spline interpolation method to generate the upper envelope and the lower envelope respectively;

[0099] B2. Calculate the average value of the upper envelope and the lower envelope;

[0100] B3. Subtract the average value from the original vibration signal to obtain h1(t), where h1(t) represents the vibration signal after normalization processing;

[0101] B4. If \(h_1(t)\) satisfies the conditions of the Intrinsic Mode Function (IMF), then \(h_1(t)\) is taken as the first Intrinsic Mode Function (IMF). The conditions for the Intrinsic Mode Function (IMF) are that the number of extreme points and zero-crossing points are equal or differ by one, and the average of the envelopes of the local maxima and minima is zero. The extreme points include the maxima and minima of the original vibration signal segment.

[0102] B5. Subtract the first Intrinsic Mode Function (IMF) from the original vibration signal to obtain the remaining signal \(r_1(t)\).

[0103] B6. Determine whether the remaining signal \(r_1(t)\) is a monotonic function. If so, represent the original vibration signal as the sum of several Intrinsic Mode Functions (IMFs) and a residual term, and proceed to step B7. Otherwise, return to step B1 until the remaining signal \(r(t)\) is a monotonic function. n (t) is a monotonic function.

[0104] B7. Perform the Hilbert transform on the \(i\)-th Intrinsic Mode Function (IMF) to obtain the analytic signal \(z(t)\). i (t);

[0105] B8. Based on the analytic signal \(z(t)\), obtain the instantaneous amplitude \(a(t)\) and the instantaneous phase \(\varphi(t)\). i (t), obtain the instantaneous amplitude \(a(t)\) i (t) and the instantaneous phase \(\varphi\) i (t);

[0106] B9. According to the instantaneous phase \(\varphi(t)\), obtain the instantaneous frequency \(\omega(t)\). i (t), obtain the instantaneous frequency \(\omega\) i (t);

[0107] B10. According to the instantaneous amplitude \(a(t)\) and the instantaneous frequency \(\omega(t)\), obtain the Hilbert spectrum \(H(\omega,t)\). The Hilbert spectrum \(H(\omega,t)\) represents the distribution of the instantaneous amplitude \(a(t)\) at different frequencies and times. i (t) and the instantaneous frequency \(\omega\) i (t), obtain the Hilbert spectrum \(H(\omega,t)\). The Hilbert spectrum \(H(\omega,t)\) represents the distribution of the instantaneous amplitude \(a\) i (t) at different frequencies and times;

[0108] B11. From each Hilbert spectrum \(H(\omega,t)\), obtain the total energy, frequency center, frequency bandwidth, and maximum energy frequency.

[0109] B12. Fuse the results of B11 to obtain the vibration index VI.

[0110] In this embodiment, health state indicators of the bearing in multiple aspects are collected, including the vibration index (VI) and the temperature index (TI), and the comprehensive health score (CHS) at a certain time step is calculated comprehensively. The comprehensive health score CHS is the weighted sum of each health indicator. Among them, the vibration index VI and the temperature index TI have a significant impact on the health state in most scenarios. The formula for evaluating the bearing state at a certain time step is as follows:

[0111] CHS = w V ·VI + w T ·TI

[0112] where w V and w T both represent learnable parameters.

[0113] In this embodiment, the overall extraction process of the vibration index VI is as follows:

[0114] First, vibration signal acquisition and preprocessing are carried out. Vibration signals within the complete life cycle are collected from the bearing sensor, and preprocessing such as denoising and filtering is performed to extract time-domain features reflecting the basic characteristics of the signal. Subsequently, time-frequency domain feature extraction is carried out. The Hilbert-Huang Transform (HHT) is used to analyze the non-linear and non-stationary vibration signals. The Hilbert-Huang Transform HHT consists of the following two parts: Empirical Mode Decomposition (EMD) decomposes the original signal into several Intrinsic Mode Functions (IMFs) to capture the intrinsic characteristics of the vibration signal; Hilbert Spectrum Analysis performs the Hilbert transform on each IMF, calculates the instantaneous amplitude and frequency, generates the Hilbert spectrum, and reveals the energy distribution of the signal in time and frequency.

[0115] In this embodiment, the processing process of the Empirical Mode Decomposition EMD is as follows:

[0116] For the original signal x(t), all the maximum and minimum values of the original signal segment are found, and the upper envelope u(t) and the lower envelope l(t) are respectively generated using the cubic spline interpolation method. The calculation formulas are as follows:

[0117] u(t) = spline(t, maxima)

[0118] l(t) = spline(t, minima)

[0119] Among them, \(t\) represents the time variable, maxima represents all the maximum points of the original signal \(x(t)\), minima represents all the minimum points of the original signal \(x(t)\), and spline represents the cubic spline interpolation function, which is used to generate a smooth envelope curve.

[0120] The average value \(m1(t)\) of the upper envelope and the lower envelope calculated by the following formula:

[0121]

[0122] Subsequently, the first Intrinsic Mode Function (IMF) is extracted. The IMF is several Intrinsic Mode Functions (IMFs) obtained by decomposing the original signal. The IMFs have the following characteristics: within the entire data segment, the number of extreme points (maximum and minimum points) is equal to or at most one different from the number of zero-crossing points; at any time point, the average value of the envelope of the local maximum and the envelope of the local minimum is zero.

[0123] Extract the first Intrinsic Mode Function (IMF), subtract the local mean from the original signal to obtain \(h1(t)\):

[0124] \(h1(t)=x(t)-m1(t)\)

[0125] Check whether \(h1(t)\) is an Intrinsic Mode Function (IMF): If \(h1(t)\) meets the conditions of the Intrinsic Mode Function (IMF) (the number of extreme points and zero-crossing points is equal or different by one, and the average value of the envelopes of the local maximum and local minimum is zero), then it is taken as the first Intrinsic Mode Function (IMF); otherwise, repeat the steps until the conditions are met.

[0126] After obtaining the first Intrinsic Mode Function (IMF), subtract the first Intrinsic Mode Function (IMF) from the original signal to obtain the remaining signal \(r1(t)\):

[0127] \(r1(t)=x(t)-IMF1(t)\)

[0128] Subsequently, iterate the above process until the remaining signal \(r\) n (t) becomes a monotonic function. Finally, the original signal \(x(t)\) can be expressed as the sum of several Intrinsic Mode Functions (IMFs) and a residual term

[0129]

[0130] Among them, \(x(t)\) represents the original vibration signal, \(i\) represents the index variable, \(n'\) represents the upper limit value of the division, and \(IMF\) i (t) represents the \(i\)-th Intrinsic Mode Function (IMF).

[0131] Next, perform Hilbert spectral analysis, that is, Hilbert transform. The calculation process is as follows:

[0132] Perform the Hilbert transform on the i-th intrinsic mode function IMF, IMFi(t), to obtain its analytic signal z i (t):

[0133] z i (t) = IMF i (t) + jH[IMF i (t)]

[0134]

[0135] where j represents the imaginary unit, H[] represents the Hilbert transform, P.V. represents the Cauchy principal value function, τ represents the integration variable, t represents the time parameter of the signal, and d represents the differential symbol.

[0136] Subsequently, the instantaneous amplitude a i (t) and the instantaneous phase φ i (t) can be obtained. First, represent the obtained analytic signal in exponential form:

[0137]

[0138] where the instantaneous amplitude ai(t) and the instantaneous phase φi(t) are respectively:

[0139]

[0140] Furthermore, differentiate the instantaneous phase to obtain the instantaneous frequency:

[0141]

[0142] From the above, the Hilbert spectrum H(ω,t) is obtained. The Hilbert spectrum H(ω,t) is the distribution of the instantaneous amplitude a i (t) at different frequencies and times. The calculation formula is as follows:

[0143]

[0144] where δ() represents the Dirac function, ω represents the frequency variable, which is used to represent the distribution of the signal at different frequencies.

[0145] In summary, n intrinsic mode functions (IMFs) have been obtained through empirical mode decomposition (EMD), and the Hilbert transform has been performed on each intrinsic mode function IMF to obtain the Hilbert spectrum H(ω,t). Construct the eigenvector of the vibration index VI as follows:

[0146] The total energy, frequency center, frequency bandwidth, and maximum energy frequency are obtained from each Hilbert spectrum, and the calculation formulas are as follows:

[0147]

[0148] where E i represents the energy of the i-th component, T represents the upper limit of the time variable, d represents the differential symbol, ω represents the frequency variable, and CenterFreq i represents the center frequency of the i-th component, Bandwidth i represents the bandwidth of the i-th component, and MaxEnergyFreq i represents the maximum energy frequency of the i-th component.

[0149] The vibration index VI is obtained by fusing the above values:

[0150] VI = [E1, CenterFreq1, Bandwidth1, MaxEnergyFreq1,...,

[0151] E n , CenterFreq n , Bandwidth n , MaxEnergyFreq n

[0152] where E n represents the energy of the last component, CenterFreq n represents the center frequency of the i-th component, Bandwidth n represents the bandwidth of the i-th component, and MaxEnergyFreq n represents the maximum energy frequency of the i-th component.

[0153] In this embodiment, for the feature extraction of the temperature index, for n time windows, 4 features are extracted from each time window to keep the dimension of the temperature index TI consistent with that of the vibration index VI. The average temperature, temperature standard deviation, maximum temperature, and minimum temperature of each time window are calculated respectively. The temperature quadruple of each time window is calculated as follows:

[0154]

[0155] where, represents the average temperature in the n-th time period, represents the standard deviation in the n-th time period, represents the maximum temperature in the n-th time period, represents the minimum temperature in the n-th time period, ​represents the average temperature within the \(i^{th}\) time period, \(T\) i' represents the length of the \(i^{th}\) time period, \(T(t)\) represents the temperature at time \(t\), \(d\) represents the differential symbol, \(t\) represents the time variable, \(t\) i' represents the start time of the \(i^{th}\) time period, \(\Delta t\) represents the time interval represents the standard deviation within the \(i^{th}\) time period represents the maximum temperature within the \(i^{th}\) time period represents the minimum temperature within the \(i^{th}\) time period

[0156] Combine the feature vectors of all time windows into a total feature vector \(TI\):

[0157]

[0158] Finally, by weighted fusion of the vibration index \(VI\) and the temperature index \(TI\), the comprehensive bearing condition score can be obtained

[0159] S2. Based on the comprehensive health score \(CHS\), construct a prediction model for the remaining useful life of the bearing health state; the prediction model includes:

[0160] Input data preprocessing module, used to perform normalization processing on the comprehensive health score \(CHS\) and cut it into time series segments of a fixed size, specifically:

[0161] C1. Perform normalization processing on the comprehensive health score \(CHS\);

[0162] C2. Use a Gaussian filter to remove data noise in the normalized comprehensive health score \(CHS\);

[0163] C3. Cut the comprehensive health score \(CHS\) processed by C2 into time series segments of a fixed size;

[0164] CNN residual network, used to extract low-level features based on time series segments using the initial convolutional layer, and based on the extracted low-level features, use two residual blocks to extract high-level features by increasing the number of channels, and use the global average pooling layer to compress the feature dimension to extract the spatial features of the bearing health state data;

[0165] Autoformer module, used to separate the trend term and the periodic term based on the spatial features of the bearing health state data, and generate the bearing life prediction result by aggregating time series features, specifically:

[0166] D1. Perform time series decomposition based on the spatial features of the bearing health state data;

[0167] D2. Perform trend extraction and analysis based on the time series decomposition result;

[0168] D3. Based on the trend extraction analysis results, use the autocorrelation mechanism to capture the periodicity and autocorrelation in the time series;

[0169] D4. Based on the autocorrelation, construct an attention mechanism to focus on the time points of the bearing health state with strong correlation;

[0170] D5. Based on the attention results, perform aggregation of the time series features;

[0171] D6. Based on the aggregated time series features, generate the bearing life prediction results;

[0172] A prediction output module for outputting the bearing life prediction results.

[0173] In this embodiment, the prediction model consists of multiple parts, including a model input data preprocessing module, a CNN residual network module, an Autoformer module, and an output module. The input preprocessing module receives the original bearing health state data, performs weighted fusion on the vibration index and temperature index of the bearing to form a health state vector, which is done using a fully connected layer. Subsequently, the data is normalized, denoised, and cut into time series segments of a fixed size.

[0174] In this embodiment, normalization means scaling the obtained comprehensive health data to improve the convergence speed of the model during training. The formula is as follows:

[0175]

[0176] Among them, x represents the obtained multi-dimensional tensor CHS.

[0177] Subsequently, a Gaussian filter is used to remove the noise in the data. The Gaussian filter is a commonly used smoothing filter and can be implemented through convolution operations. The calculation formula is as follows:

[0178]

[0179] Among them, y(t) represents the signal processed by the Gaussian filter, τ represents the integration variable and the time offset, t represents the time variable and the current time point, d represents the differential symbol for differentiating the variable, and g(t) represents the Gaussian kernel function. The calculation formula is as follows:

[0180]

[0181] Among them, σ represents the standard deviation of the Gaussian kernel.

[0182] The obtained data is cut into time series segments of a fixed size. The cutting formula is as follows:

[0183] X i'= x[(i'-1)S:(i'-1)S + L]

[0184] Wherein, X i represents the i'-th time series segment, x represents the original time series data, S represents the step size, and L represents the fixed length of each segment.

[0185] The convolutional residual network is used to effectively extract the spatial features of bearing health state data. The CNN residual network adopted in the present invention consists of an initial convolutional layer and two residual blocks, a total of three layers. The initial convolutional layer extracts low-level features with 64 3×3 convolutional kernels. Subsequently, it passes through two residual blocks. Each residual block contains two convolutional layers, a BatchNormalization layer, a ReLU activation function, and a residual connection. After passing through the residual blocks, the number of channels is initially increased (from 64 to 128) to extract higher-level features. Finally, global average pooling is used to compress the feature dimension. The detailed calculation process is as follows:

[0186] a. Both the initial convolutional layer and the residual convolution use basic convolutional operations:

[0187]

[0188] Wherein, Y(i, j) represents the value of the output feature map at position (i, j), m and n represent the indices of the convolutional kernel, i and j represent the position indices of the output feature map, X represents the input feature map, K represents the convolutional kernel, Y represents the output feature map, and M and N both represent the size of the convolutional kernel (both are 3 in the present invention).

[0189] b. Perform residual connection. Residual connection can solve the problem of gradient disappearance and gradient explosion in the process of training a neural network, making it easier to train a very deep network architecture where the input of the network is directly added to the output of this layer. Each residual block used in the present invention contains two convolutional layers. Each convolutional layer is followed by a batch normalization layer and a ReLU activation function. The residual connection directly adds the input to the output of this layer. After the first residual block, the number of channels is increased from 64 to 128 through a 1×1 convolution to extract higher-level features. The calculation formula of the ReLU function is as follows:

[0190] f(x) = max(0, x)

[0191] Wherein, x represents the input of the function.

[0192] The formula for the residual process is as follows:

[0193] H(x) = F(x) + x

[0194] H(x) = F(x, {W i}) + W s x

[0195] Among them, H(x) represents the result after residual connection, F(x) represents the residual mapping, x represents the input of the function, and W i represents the set of weights of the residual path, and W s represents the weight matrix of the shortcut mapping.

[0196] c. Subsequently, batch normalization is performed for feature standardization calculation:

[0197]

[0198] Among them, represents the standardized result of the k-th feature channel, x (k) represents the k-th feature channel, E[x (k) represents the expected value of the k-th feature channel, y (k) represents the result after standardization and scaling of the k-th feature channel. Var[] represents the variance, and ò represents a small constant to prevent division by zero, usually randomly initialized. γ (k) and β (k) both represent learnable parameters.

[0199] d. Finally, global average pooling is performed to compress the feature dimension and obtain a fixed-length feature vector. The calculation formula is as follows:

[0200]

[0201] Among them, GAP(X) represents the output of global average pooling, that is, the finally extracted feature. H and W respectively represent the height and width of the feature map, and X ij represents the value at the position (i, j) in the feature map.

[0202] In this embodiment, the Autoformer network is a deep learning architecture for time series prediction. For the model architecture, see Figure 3 . Autoformer embeds a sequence decomposition unit into the deep model to achieve progressive decomposition. During the prediction process, the trend term and the periodic term are gradually separated from the latent variables, and the prediction result optimization and sequence decomposition are alternately performed to promote each other, so as to better handle complex time patterns. Based on the theory of stochastic processes, Autoformer discards the self-attention mechanism of point-to-point connection and designs a self-correlation mechanism of sequence-level connection. This mechanism includes period-based dependence discovery and time-delay information aggregation, which can utilize the inherent periodicity of the sequence, mine the sub-processes between similar phases of different periods, achieve efficient information aggregation, break through the information utilization bottleneck, and have a low complexity. The detailed calculation steps are as follows:

[0203] First, perform time series decomposition on the output obtained from the CNN residual network. The formula is as follows:

[0204] X t = T t + S t + R t

[0205] Among them, X t represents the original time series value at time t, T t represents the trend term, reflecting the long-term change trend of the time series, S t represents the seasonal term, reflecting the periodic change of the time series, R t represents the random term, reflecting the random fluctuation of the time series. This step decomposes the time series into trend term, seasonal term and random term, which is convenient for separate modeling.

[0206] Perform trend extraction analysis on the data:

[0207]

[0208] Among them, k represents the moving average window size, X t+i represents the original time series value at time t + i', and i represents the offset relative to time t.

[0209] Subsequently, perform the calculation steps of the auto-correlation mechanism to capture the periodic patterns and auto-correlations in the time series:

[0210] Auto-correlation calculation:

[0211]

[0212] Among them: μ Q and μ K represent the means of Q and K respectively, σ Q and σ K represent the standard deviations of Q and K respectively, and AutoCorr(Q, K) represents the auto-correlation value between the query matrix Q and the key matrix K.

[0213] Construct an attention mechanism based on auto-correlation to focus on the time points of bearing health status with strong correlations.

[0214]

[0215] Among them, Q represents the query matrix, K represents the key matrix, V represents the value matrix, and d krepresents the scaling factor of the attention mechanism, AutoCorr() represents the autocorrelation operation, and A() represents the output of the attention mechanism.

[0216] Finally, perform the aggregation of temporal features:

[0217]

[0218] Among them, F t represents the aggregated time series features, S represents the number of time scales considered, w s represents the weights of each time scale, X t-s:t represents the time window data from t - s to t, and θ s represents the parameters of the one-dimensional convolution.

[0219] Generate the final prediction based on the extracted features:

[0220]

[0221] Among them, represents the bearing life prediction result at time t, φ represents the activation function of the output layer, both W1 and W2 represent weight matrices, both b1 and b2 represent bias terms, and σ represents the activation function of the intermediate layer.

[0222] S3. Use the comprehensive health score CHS to train the constructed prediction model;

[0223] In this embodiment, using the input data prepared in step S1 and the constructed CRN - Autoformer model, the mean squared error (MSE) is used as the loss function during the training phase, and the Adam optimizer is used to optimize the network parameters. The calculation formula is as follows:

[0224]

[0225] Among them, N represents the number of samples, y i represents the true value of the i - th sample, represents the predicted value of the i - th sample.

[0226]

[0227] Among them, m t represents the first - order moment estimate (i.e., the exponential weighted average of the gradient), β1 and β2 respectively represent the exponential decay rates of the first - order and second - order moment estimates, m t-1 represents the previous first - order moment estimate, g t represents the gradient at the current time step, v t represents the second - order moment estimate (i.e., the exponential weighted average of the gradient square), v t-1Represents the second moment estimate of the previous step, Represents the first moment estimate after bias correction, Represents the t-th power of the exponential decay rate of the first moment estimate, Represents the second moment estimate after bias correction, Represents the t-th power of the exponential decay rate of the second moment estimate, θ t+1 and θ t respectively represent the updated parameter and the parameter before update. α represents the learning rate, set to 0.01, and ò represents a very small value, randomly initialized, used to ensure the stability of the division operation.

[0228] S4. Use the trained prediction model to identify the dynamic characteristics of the bearing health state with respect to the remaining useful life, and predict the bearing health state based on the trend of the historical health state scores, thereby completing the prediction of the bearing life.

[0229] In this embodiment, through the above steps to train the model, the CRN-Autoformer prediction model can identify the dynamic characteristics of the bearing health state with respect to the usage time series, and predict the subsequent bearing health state based on the trend of the historical health state scores, that is, the estimation of the available duration of the bearing.

[0230] In summary, the hybrid data-driven model (CRN-Autoformer) based on convolutional residual network and Autoformer proposed by the present invention has significant advantages in multiple aspects and solves the key problems in the prior art. First, by introducing the Hilbert-Huang transform (HHT) to extract the time-frequency domain features of the bearing vibration signal, it can adapt to the characteristics of non-linear and non-stationary signals. Especially during the fault development process, it can accurately capture the dynamic characteristics of the vibration signal changing with time and frequency, providing high-quality input data for the model. Second, using the convolutional residual network (CRN) as the feature extraction module, the problem of gradient disappearance in the training of deep networks is effectively alleviated through residual connections, significantly improving the extraction ability of the temporal features of the bearing vibration signal and the spatial features of the temperature data, thereby comprehensively characterizing the bearing health state. In addition, the introduced Autoformer module combines the time series decomposition mechanism to specifically identify trend and seasonal features, significantly enhancing the modeling ability for long-range dependencies in time series data, enabling the bearing life prediction to maintain high accuracy and stability over a long time span. At the same time, the present invention utilizes the idea of transfer learning to achieve efficient modeling by fine-tuning the pre-trained model, without uploading a large amount of original data, ensuring data security and reducing computational resource consumption. Combining these technical advantages, the CRN-Autoformer model exhibits high accuracy, strong robustness and excellent generalization ability in complex industrial environments, providing an innovative solution for bearing life prediction.

Claims

1. A bearing life prediction method, characterized in that, It includes the following steps: S1. Obtain the bearing health status indicators and comprehensively calculate the comprehensive health score CHS at a certain time step; S2. Based on the comprehensive health score CHS, construct a prediction model for the remaining useful life of the bearing health status; S3. Use the comprehensive health score CHS to train the constructed prediction model; S4. Use the trained prediction model to identify the dynamic characteristics of the bearing health status with respect to the remaining useful life, and predict the bearing health status based on the trend of the historical health status scores to complete the prediction of the bearing life.

2. The bearing life prediction method according to claim 1, wherein Specifically in step S1: Obtain the bearing health status indicators including the vibration index VI and the temperature index TI, and comprehensively calculate the comprehensive health score CHS at a certain time step.

3. The bearing life prediction method according to claim 2, characterized in that The expression of the comprehensive health score CHS is as follows: CHS = w V ·VI + w T ·TI where, w V and w T both represent learnable parameters.

4. The bearing life prediction method according to claim 2, wherein, The acquisition process of the vibration index VI is as follows: A1. Collect vibration signals during the complete life cycle from the bearing sensor, and preprocess the vibration signals to obtain the original vibration signals for extracting time-domain features; A2. Use the Hilbert-Huang transform HHT to perform non-linear and non-stationary signal analysis on the original vibration signals, and obtain the vibration index VI based on the Hilbert spectrum generated from the analysis results. Among them, the Hilbert-Huang transform HHT includes empirical mode decomposition EMD and Hilbert spectrum analysis HAS. The empirical mode decomposition EMD decomposes the original vibration signals into several intrinsic mode functions IMF to capture the inherent characteristics of the vibration signals; the Hilbert spectrum analysis HAS performs the Hilbert transform on each intrinsic mode function IMF, and generates the Hilbert spectrum by calculating the instantaneous amplitude and instantaneous frequency to reveal the energy distribution of the signals in time and frequency.

5. The bearing life prediction method according to claim 4, characterized in that, The specific process of A2 is as follows: B1. Obtain all the maximum and minimum values of the original vibration signal segment, and use the cubic spline interpolation method to generate the upper envelope and the lower envelope respectively; B2. Calculate the average value of the upper envelope and the lower envelope; B3. Subtract the average value from the original vibration signal to obtain h1(t), where h1(t) represents the vibration signal after normalization processing; B4. When h1(t) meets the conditions of the intrinsic mode function IMF, then take h1(t) as the first intrinsic mode function IMF. Among them, the conditions of the intrinsic mode function IMF are that the number of extreme points and zero-crossing points are equal or differ by one, and the average value of the envelopes of the local maximum and minimum values is zero. The extreme points include the maximum and minimum values of the original vibration signal segment; B5. Subtract the first intrinsic mode function IMF from the original vibration signal to obtain the remaining signal r1(t); B6. Determine whether the remaining signal r1(t) is a monotonic function. If so, use the following formula to represent the original vibration signal as the sum of a number of intrinsic mode functions IMF and a residual term, and proceed to step B7. Otherwise, return to step B1 until the remaining signal r n (t) is a monotonic function: Among them, x(t) represents the original vibration signal, i represents the index variable, n' represents the upper limit value of the division, and IMF i (t) represents the i-th intrinsic mode function IMF; B7. Using the following formula, perform the Hilbert transform on the \(i\)-th intrinsic mode function IMF to obtain the analytic signal \(z(t)\): i (t): z i (t) = IMF i (t) + jH[IMF i (t)] Where j represents the imaginary unit, H[] represents the Hilbert transform, P.V. represents the Cauchy principal value function, τ represents the integration variable, t represents the time parameter of the signal, and d represents the differential symbol; B8. Based on the analytic signal z i (t), the instantaneous amplitude a i (t) and the instantaneous phase φ i (t) are obtained by using the following formula: B9. According to the instantaneous phase φ i (t), the instantaneous frequency ω i (t) is obtained by using the following formula: B10. According to the instantaneous amplitude a i (t) and the instantaneous frequency ω i (t), the Hilbert spectrum H(ω, t) is obtained by using the following formula. The Hilbert spectrum H(ω, t) represents the distribution of the instantaneous amplitude a i (t) at different frequencies and times: Where δ() represents the Dirac function, ω represents the frequency variable, which is used to represent the distribution of the signal at different frequencies; B11. From each Hilbert spectrum H(ω,t), obtain the total energy, frequency center, frequency bandwidth, and maximum energy frequency using the following formula: Among them, E i represents the energy of the i-th component, T represents the upper limit of the time variable, d represents the differential symbol, ω represents the frequency variable, and CenterFreq i represents the center frequency of the i-th component, Bandwidth i represents the bandwidth of the i-th component, and MaxEnergyFreq i represents the maximum energy frequency of the i-th component; B12. Use the following formula to fuse the results of B11 to obtain the vibration index VI: VI = [E1, CenterFreq1, Bandwidth1, MaxEnergyFreq1,... E n , CenterFreq n , Bandwidth n , MaxEnergyFreq n ​ Among them, E n represents the energy of the last component, and CenterFreq n represents the center frequency of the i-th component, Bandwidth n represents the bandwidth of the i-th component, and MaxEnergyFreq n represents the maximum energy frequency of the i-th component.

6. The bearing life prediction method according to claim 2, characterized in that The expression of the temperature index TI is as follows: Among them, represents the average temperature in the nth time period, represents the standard deviation in the nth time period, represents the maximum temperature in the nth time period, represents the minimum temperature in the nth time period, represents the average temperature in the i'th time period, T i' represents the length of the i'th time period, T(t) represents the temperature at time t, d represents the differential symbol, t represents the time variable, t i' represents the start time of the i'th time period, Δt represents the time interval, represents the standard deviation in the i'th time period, represents the maximum temperature in the i'th time period, represents the minimum temperature in the i'th time period.

7. The bearing life prediction method according to claim 1, characterized in that The prediction model includes: An input data preprocessing module, which is used to perform normalization processing on the comprehensive health score CHS and cut it into time series segments of a fixed size; A CNN residual network, which is used to extract low-level features based on the time series segments using an initial convolutional layer, and based on the extracted low-level features, use two residual blocks to extract high-level features by increasing the number of channels, and use a global average pooling layer to compress the feature dimension to extract the spatial features of the bearing health state data; An Autoformer module, which is used to separate the trend term and the periodic term based on the spatial features of the bearing health state data and generate a bearing life prediction result by aggregating the time series features; A prediction output module, which is used to output the bearing life prediction result.

8. The bearing life prediction method according to claim 7, wherein The normalization processing on the comprehensive health score CHS and cutting it into time series segments of a fixed size are specifically as follows: C1. Perform normalization processing on the comprehensive health score CHS; C2. Use a Gaussian filter to remove the data noise in the normalized comprehensive health score CHS; C3. Use the following formula to cut the comprehensive health score CHS processed by C2 into time series segments of a fixed size: X i' = x[(i'-1)S:(i'-1)S+L] Among them, X i represents the i'-th time series segment, x represents the original time series data, S represents the step size, and L represents the fixed length of each segment.

9. The bearing life prediction method according to claim 7, wherein The separation of the trend term and the periodic term based on the features after the compressed feature dimension and the generation of the bearing life prediction result by aggregating the time series features are specifically as follows: D1. Based on the spatial features of the bearing health state data, perform time series decomposition using the following formula: X t = T t + S t + R t where X t represents the original time series value at time t, T t represents the trend term, reflecting the long-term change trend of the time series, S t represents the seasonal term, reflecting the periodic change of the time series, R t represents the random term, reflecting the random fluctuation of the time series; D2. Based on the time series decomposition result, perform trend extraction analysis using the following formula: where k represents the moving average window size, and X t+i represents the original time series value at time t + i'; D3. Based on the trend extraction analysis result, use an autocorrelation mechanism to capture the periodicity and autocorrelation in the time series; D4. Based on the autocorrelation, construct an attention mechanism to focus on the bearing health state time points with strong correlation; D5. Based on the attention result, perform aggregation of the time series features using the following formula: Among them, F t represents the aggregated time series features, S represents the number of time scales considered, w s represents the weights of each time scale, X t-s:t represents the time window data from t - s to t, and θ s represents the parameters of one-dimensional convolution; D6. Based on the aggregated time series features, generate a bearing life prediction result using the following formula: Among them, represents the bearing life prediction result at time t, φ represents the output layer activation function, both W1 and W2 represent weight matrices, both b1 and b2 represent bias terms, and σ represents the intermediate layer activation function.

Images