1. The bearing early fault detection method based on health state index construction and fault early warning limit self-learning comprises the following steps of:

s1: time-frequency domain fault feature extraction

S11: acquiring time-frequency domain information in a non-stationary original vibration signal of a rolling bearing, namely the change condition of signal frequency along with a time sequence, and instantaneous frequency and amplitude thereof at each moment;

s12: carrying out gray processing on the obtained time-frequency image so as to facilitate subsequent model input, processing the three-channel image into a single-channel gray image, calling an RGB2gray (RGB) function in MATLAB, eliminating hue and saturation information, simultaneously keeping brightness, and converting the RGB image into a gray image; by calculating R, G a weighted sum of the B components: 0.2989 × R +0.5870 × G +0.1140 × B, converting RGB values into gray values, and processing image pixels into 32 × 32;

s2: MVGCDBN-based health state characterization index construction

Adopting a matrix variable Gaussian convolution depth confidence network MVGCDBN for constructing the health state index:

the model adopts a visual layerV and a hidden layer H; the visible layer is an N_V×N_VA binary cell matrix including K N_W×N_W(N_W＝N_V-N_H+1) convolution kernel W_k(ii) a The hidden layer is composed of K N_H×N_HThe feature mapping surface composition of the binary unit matrix marks the k-th feature mapping surface H^kHas a convolution kernel of W^kOffset is b_kAll visible layers share bias a;

assuming that each cell of the visible layer and the hidden layer obeys a variance of σ₁ ²And gamma₁ ²The energy function of the model is as follows (1):

wherein v and h represent the state vectors of all visible layers and hidden layers respectively; v. of_i,jIs an element of the ith row and jth column of the matrix v,is the element of the ith row and jth column of the kth hidden unit,is an element of the kth filter row r column s;is thatThe matrix of the flip-flop of (c),denotes convolution operation, denotes the product of elements and then sums, i.e. a · B ═ trA^TB；

The hidden layer and visible layer unit conditional activation probabilities are respectively:

the purpose of model training is to determine an optimal parameter set, and parameters of the model are updated by adopting a contrast divergence CD algorithm based on Gibbs sampling; sampling hidden layer nodes according to the calculated hidden layer conditional probability, and taking the average value of the hidden layer nodes asVariance is σ₁ ²The Gaussian distribution is sampled, the computed probability value is not compared with a random number from 0 to 1 to obtain a value of 0 or 1, and the sampled data obtained by the method is not 0 or 1 but a real number sampled near the mean value; the sampling method of the visible layer node is the same as that of the hidden layer node;

inputting: training samples, namely a visual layer vector v, the number K of filters in a hidden layer, a learning rate, a regularization constant zeta, a sparsity coefficient rho and a maximum iteration number Q;

and (3) outputting: connection weight parameter W^kA visible layer bias a and a hidden layer bias b_k；

Random initialization: connection weight parameter W^kA visible layer bias a and a hidden layer bias b_kIs a random smaller number;

for q＝1 to Q do

for k＝1 to K do

step 1 updating weight parameters

Step 2 updating weights with gradients of regularization term

Step 3, repeating Step 1 and Step 2 until convergence or Q ═ Q

end for

end for

S3: self-learning fault early warning control limit

S31: fitting of distributions

(1) Normalizing the HI data in the 'normal' stage;

(2) calculating the shape parameters of normal distribution and gamma distribution of HI data in a 'normal' stage by adopting maximum likelihood estimation;

if the random variable x (hi) obeys a probability distribution with a position parameter μ and a scale parameter σ, its probability density function is:

the random variable x is then said to obey the parameter μ, σ²Normal distribution of (d) is recorded as: x to N (mu, sigma)²)；

If the random variable x obeys a probability distribution with a shape parameter of α and an inverse scale parameter of β, its probability density function is:

let x-gamma (alpha, beta), and letNamely, it isWherein alpha is more than 0, beta is more than 0; then the random variable X is said to obey the gamma distribution with parameters α, β, written as: x to G (. alpha.,. beta.);

s32: and (3) testing the goodness of fit: checking the fitted distribution to determine a distribution to which the HI obeys;

s33: and calculating a threshold corresponding to the quantile on the HI to determine a normalized control limit, and then carrying out reverse normalization on the normalized control limit to obtain a fault early warning control limit, so as to realize early prediction of a fault point.

2. The method for early fault detection of a bearing based on health status indicator construction and fault early warning limit self-learning of claim 1, wherein: step S11, acquiring time-frequency domain information in the non-stationary original vibration signal of the rolling bearing by using short-time fourier transform STFT, that is, the change of the signal frequency with the time series, and the instantaneous frequency and amplitude thereof at each moment; let g_τν(t)Is a localized standard single window function centered around (tau, v), i.e. low-pass filtering, centered in the interval tau-deltag, tau + deltag in the time domain]Wherein e is^iνtIs a frequency shift factor;

g_τν(t)＝g(t-τ)e^iνtτ,ν∈R(a)

in the formula T_x(τ, ν) represents the fourier transformation of the signal within the window, x (t) is the original signal that needs to be transformed,as a part of the signal x (t) is intercepted by a window of width 2 Δ g centered on τ, while as τ moves, different parts of the signal will be taken into the window in sequence; STFT is interpreted as: the signal is first frequency shifted by x (t) e^-iνtThen go through g_τν(t)Low-pass filtering to obtain a time-domain representation at a frequency ν; wherein, a Hamming window function with the width of 25ms is adopted, and the displacement of the frame is 10 ms.

3. The method for early fault detection of a bearing based on health status indicator construction and fault early warning limit self-learning of claim 1, wherein: step S11 performs wavelet transform on the original signal x (t) to obtain time-frequency domain information in the non-stationary original vibration signal of the rolling bearing, that is, the change of the signal frequency along with the time series, and the instantaneous frequency and amplitude thereof at each time, specifically as follows:

in the formula, a₀Is a scale factor, a₀＞0；b₀Is a time-shift factor;is a wavelet basis function; the essence of wavelet transform is to use wavelet basis functionsThe signal x (t) is decomposed into sub-signals of different frequency bands by inner product operation to obtain the frequency spectrum.

4. The method for early fault detection of a bearing based on state of health indicator construction and fault early warning limit self-learning according to any one of claims 1-3, wherein: in the step S32, the distribution of the fit is checked by adopting K-S check to determine the distribution obeyed by HI; here, the significance level a₁0.05, the confidence is 1- α₁＝95％；

(1) Let f (x) be the frequency distribution, g (x) the theoretical distribution, assuming H₀: the data conforms to theoretical distribution;

(2) calculating D ═ max | f (x) -g (x) |, when the actual observed value D is>D(n，α₁) Then refuse H₀Otherwise, H is accepted₀It is assumed that, i.e. when h is 0,the original hypothesis is accepted.

5. The method for early fault detection of a bearing based on state of health indicator construction and fault early warning limit self-learning according to any one of claims 1-3, wherein: the goodness-of-fit test in the step S32 adopts Chi²The fit distribution is tested by testing, and the distribution obeyed by HI is determined; here, the significance level a₁0.05, the confidence levelIs 1-alpha₁＝95％，

In the formula (f)₀Actual observation frequency representing sample observations of the ith group; f. of_eFrequency that is theoretically expected; nbins is the number of groups.

Background

Predictive maintenance is a potential explosion point in the field of smart manufacturing. The early failure prediction of the equipment is taken as a key ring, becomes a hotspot of research in the engineering field in recent years, and develops a series of failure prediction methods with great application prospects. However, the vibration signals collected by the mechanical system are complex, non-stationary and high-background-noise signals, and the health state evolution trend is usually non-monotonous or even abrupt, so that it is a great challenge to effectively learn the health characterization indexes capable of characterizing the equipment from the vibration signals so as to detect and alarm the impending slight or abnormal fault symptoms of the equipment in advance and guide the predictive maintenance of the equipment.

Unlike many available fault diagnosis methods, early fault prediction regarding mechanical parts is in the beginning stage, and related research is also rapidly developing. Currently, early failure prediction models can be divided into physical models and data-driven models. The model-based method needs to build an accurate physical model to predict the health condition of the equipment according to the characteristics of the equipment. The limitations are the very high requirement on the prior knowledge of the study subject and the poor generalization ability. And it is very difficult to establish a complete physical failure prediction model for complex mechanical equipment which is installed in complex structure, has bad operation environment and has large mutual influence between components, so that the use and popularization of the method have great limitation.

The data-driven prediction method obtains a potential degradation rule of the equipment by analyzing sensor monitoring data of the equipment, so that the health state evolution trend of the equipment is predicted. It need only collect enough performance degradation data and convert the collected data into relevant state of health features or performance degradation models to predict future data. The data-driven prediction method has lower requirements on the relevant prior knowledge of the equipment, can well solve the prediction problem when the failure mechanism of the equipment is complex and an accurate model is difficult to establish, has high efficiency and good generalization capability of the model, and obtains more and more attention and research. Data-driven prediction methods can be further divided into three categories: the first type is a statistical Model prediction method, such as a Gray Model (GM), a Particle Filter (PF), a Bayesian Prediction Model (BPM), and the like; the second type is a shallow Machine learning prediction method, such as an Artificial Neural Network (ANN) and a Support Vector Machine (SVM); the third category is Deep learning prediction methods, such as Deep Belief Networks (DBN) and Recurrent Neural Networks (RNN).

Although the first two types of data-driven models have found widespread use in health state prediction, they still suffer from the following inherent disadvantages: SVMs and ANN are considered to have shallow structures, and it is difficult to effectively learn the highly nonlinear relationship in the health state evolution process. Typically, due to changes in operating conditions, bearing state-of-health evolution is a non-linear, non-stationary time sequence, and worse, state-of-health evolution information tends to be overwhelmed by noise during early stages of bearing failure. In such a case, it is a great challenge to use a shallow architecture for accurate fault prediction. Therefore, it is necessary to deeply study the complex nonlinear relationship in the bearing failure prediction problem.

The deep learning model is a deep neural network model with a plurality of nonlinear mapping levels, can more effectively and flexibly express any highly-changed nonlinear function than a shallow learning model, and can abstract and extract characteristics of input signals layer by layer. In other words, the deep learning method can get rid of the dependence on manual signal processing techniques and domain experts. Deep Belief Networks (DBNs), Recurrent Neural Networks (RNNs), and variants thereof are increasingly being used in machine health state prediction as a popular deep learning prediction method. The DBN is a generative neural network and has strong unsupervised feature learning capability.

Even though the past failure prediction achieves certain results, most of the attention is paid to the improvement of a prediction algorithm, and the most important premise of the failure prediction based on deep learning is to extract a characteristic index (HI) capable of comprehensively representing the health state from a vibration signal, so that the conventional method has the following defects: the existing health state assessment indexes are mostly constructed by manually extracting partial time-frequency domain statistical characteristics or entropy characteristics from original vibration signals and performing data fusion, only the characteristic quantity of a certain domain or certain domains is used as a bearing health state assessment index, the comprehensive description of the degradation information of the bearing health state assessment index cannot be met, and each original characteristic is only sensitive to a certain fault at a certain degradation stage and is difficult to accurately capture the fault information of the bearing degradation process. (2) In the vibration characteristic signal, some key characteristics are independent of each other, are sensitive to certain specific fault information of the bearing, can well describe the degradation state of the bearing, and are fully utilized. However, there are also some redundant or irrelevant features that are not sensitive enough to the degradation state of the bearing, which can reduce the computational efficiency of the fault state trend prediction and even interfere with the evaluation and prediction results. Therefore, before the feature signal set is input into the health state prediction model, sensitive features need to be screened from the feature set by using a feature evaluation method, irrelevant redundant information in the feature set is removed, and the existing feature selection and evaluation depend on the engineering experience of experts to a great extent. (3) The health status indicator is generally constructed for a specific degradation process and therefore does not generalize well to other degradation processes.

The Convolutional Neural Network (CNN) has the advantages that the CNN has the characteristic of learning capability and displacement invariance, and the trained model has stronger generalization capability through a weight sharing strategy of Convolutional connection. CNN is generally a supervised training discriminant model. The DBN is a generating neural network, has strong unsupervised feature learning capability, but requires that input data have translational invariance, and the replacement of matrix multiplication by discrete convolution of a small kernel is a standard way for solving the input problem with spatial translational invariance or a time structure. Although a CDBN (Convolutional Deep belief network) method in an unsupervised learning mode has good effects on the aspects of mechanical state monitoring and fault diagnosis, when the CDBN is applied to the construction of the bearing health state index in the prior art, a main problem still exists that hidden layer nodes take 0 or 1 value, the information of original data is lost in the hidden layer output obtained by the model to a certain extent, and the health state cannot be comprehensively represented. Therefore, the improvement of the CDBN model and the construction of the health status index have become an urgent task.

Disclosure of Invention

The technical problem to be solved by the invention is as follows: the rolling bearing early failure prediction method based on the automatic construction of the health state indexes and the early warning control limit self-learning overcomes the defects of the existing problems, and can effectively solve the problems that the failure characteristics of the bearing degradation process are difficult to completely capture, the robustness of a prediction model is not strong, the alarm leakage rate and the false alarm rate are high and the like in the existing rolling bearing failure prediction research.

In order to solve the technical problems, the technical scheme adopted by the invention is as follows: the bearing early fault detection method based on health state index construction and fault early warning limit self-learning comprises the following steps of:

s1: time-frequency domain fault feature extraction

S11: acquiring time-frequency domain information in a non-stationary original vibration signal of a rolling bearing, namely the change condition of signal frequency along with a time sequence, and instantaneous frequency and amplitude thereof at each moment;

s12: carrying out gray processing on the obtained time-frequency image so as to facilitate subsequent model input, processing the three-channel image into a single-channel gray image, calling an RGB2gray (RGB) function in MATLAB, eliminating hue and saturation information, simultaneously keeping brightness, and converting the RGB image into a gray image; by calculating R, G a weighted sum of the B components: 0.2989 × R +0.5870 × G +0.1140 × B, converting RGB values into gray values, and processing image pixels into 32 × 32;

s2: MVGCDBN-based health state characterization index construction

Adopting a matrix variable Gaussian convolution depth confidence network MVGCDBN for constructing the health state index:

the model adopts a two-layer structure consisting of a visible layer V and a hidden layer H; the visible layer is an N_V×N_VA binary cell matrix including K N_W×N_W(N_W＝N_V-N_H+1) convolution kernel Wk; the hidden layer is composed of K N_H×N_HThe feature mapping surface composition of the binary unit matrix marks the k-th feature mapping surface H^kHas a convolution kernel of W^kOffset is b_kAll visible layers share bias a;

assuming that each cell of the visible layer and the hidden layer obeys a variance of σ₁ ²And gamma₁ ²The energy function of the model is as follows (1):

wherein v and h represent the state vectors of all visible layers and hidden layers respectively; v. of_i,jIs an element of the ith row and jth column of the matrix v,is the element of the ith row and jth column of the kth hidden unit,is an element of the kth filter row r column s;is thatThe matrix of the flip-flop of (c),denotes convolution operation, denotes the product of elements and then sums, i.e. a · B ═ trA^TB；

The hidden layer and visible layer unit conditional activation probabilities are respectively:

the purpose of model training is to determine an optimal parameter set, and parameters of the model are updated by adopting a contrast divergence CD algorithm based on Gibbs sampling; sampling hidden layer nodes according to the calculated hidden layer conditional probability, and taking the average value of the hidden layer nodes asVariance is σ₁ ²The Gaussian distribution is sampled, the computed probability value is not compared with a random number from 0 to 1 to obtain a value of 0 or 1, and the sampled data obtained by the method is not 0 or 1 but a real number sampled near the mean value; the sampling method of the visible layer node is the same as that of the hidden layer node;

inputting: training samples, namely a visual layer vector v, the number K of filters in a hidden layer, a learning rate, a regularization constant zeta, a sparsity coefficient rho and a maximum iteration number Q;

and (3) outputting: connection weight parameter W^kA visible layer bias a and a hidden layer bias b_k；

Random initialization: connection weight parameter W^kA visible layer bias a and a hidden layer bias b_kIs a random smaller number;

s3: self-learning fault early warning control limit

S31: fitting of distributions

(1) Normalizing the HI data in the 'normal' stage;

(2) calculating the shape parameters of normal distribution and gamma distribution of HI data in a 'normal' stage by adopting maximum likelihood estimation;

if the random variable x (hi) obeys a probability distribution with a position parameter μ and a scale parameter σ, its probability density function is:

the random variable x is then said to obey the parameter μ, σ²Normal distribution of (d) is recorded as: x to N (mu, sigma)²)；

If the random variable x obeys a probability distribution with a shape parameter of α and an inverse scale parameter of β, its probability density function is:

let x-gamma (alpha, beta), and letNamely, it isWherein alpha is more than 0, beta is more than 0; then the random variable X is said to obey the gamma distribution with parameters α, β, written as: x to G (. alpha.,. beta.);

s32: and (3) testing the goodness of fit: checking the fitted distribution to determine a distribution to which the HI obeys;

s33: and calculating a threshold corresponding to the quantile on the HI to determine a normalized control limit, and then carrying out reverse normalization on the normalized control limit to obtain a fault early warning control limit, so as to realize early prediction of a fault point.

As a preferable scheme, the step S11 adopts short-time fourier transform STFT to obtain time-frequency domain information in the non-stationary original vibration signal of the rolling bearing, that is, the situation that the signal frequency changes along with the time series, and the instantaneous frequency and amplitude thereof at each moment; let g_τν(t)Is a localized standard single window function centered around (tau, v), i.e. low-pass filtering, centered in the interval tau-deltag, tau + deltag in the time domain]Wherein e is^iνtIs a frequency shift factor;

g_τν(t)＝g(t-τ)e^iνtτ,ν∈R (a)

in the formula T_x(τ, ν) represents the fourier transformation of the signal within the window, x (t) is the original signal that needs to be transformed,as a part of the signal x (t) is intercepted by a window of width 2 Δ g centered on τ, while as τ moves, different parts of the signal will be taken into the window in sequence; STFT is interpreted as: the signal is first frequency shifted by x (t) e^-iνtThen go through g_τν(t)Low-pass filtering to obtain a time-domain representation at a frequency ν; wherein, a Hamming window function with the width of 25ms is adopted, and the displacement of the frame is 10 ms.

As another preferable scheme, in step S11, the wavelet transform is performed on the original signal x (t) to obtain time-frequency domain information in the non-stationary original vibration signal of the rolling bearing, that is, the case that the signal frequency changes with the time series and the instantaneous frequency and amplitude thereof at each time, specifically as follows:

in the formula, a₀Is a scale factor, a₀＞0；b₀Is a time-shift factor;is a wavelet basis function; the essence of wavelet transform is to use wavelet basis functionsThe signal x (t) is decomposed into sub-signals of different frequency bands by inner product operation to obtain the frequency spectrum.

As a preferred scheme, the goodness-of-fit test of step S32 uses K-S test to test the fitted distribution, and determines the distribution to which HI obeys; here, the significance level a₁0.05, the confidence is 1- α₁＝95％；

(1) Let f (x) be the frequency distribution, g (x) the theoretical distribution, assuming H₀: the data conforms to theoretical distribution;

(2) calculating D ═ max | f (x) -g (x) |, when the actual observed value D > D (n, alpha)₁) Then refuse H₀Otherwise, H is accepted₀Suppose thatThe original hypothesis is accepted.

As another preferable scheme, the goodness-of-fit test in step S32 uses χ²The fit distribution is tested by testing, and the distribution obeyed by HI is determined; here, the significance level a₁0.05, the confidence is 1- α₁＝95％，

In the formula (f)₀Actual observation frequency representing sample observations of the ith group; f. of_eFrequency that is theoretically expected; nbins is the number of groups; as can be seen from equation (d), χ is the observed result if it is consistent with the expected result²When the fitting degree is 0, the fitting degree is in a perfect condition; when x²If the accumulated contribution exceeds the set confidence range, rejecting the original hypothesis that the data conforms to the theoretical distribution, χ²The larger the size the less acceptable the original assumption.

The invention has the beneficial effects that:

the invention provides a bearing early fault prediction method based on health state index construction and fault early warning control limit self-learning, which is characterized in that an early fault prediction model capable of running at the end of edge computing equipment is constructed, only one section of historical data of a normal running state is needed, early fault information can be accurately obtained without external expert intervention and manual analysis, and technical support is provided for realizing predictive maintenance in an information physical system environment; the problem that the established health state evaluation index cannot meet the comprehensive description of degradation information of the health state evaluation index in the existing research and can hardly and accurately capture fault information of a bearing degradation process is solved.

Compared with the traditional fixed threshold limit alarm mode, the invention can greatly advance the fault prediction time, reduce the false alarm rate and the missing alarm rate, has better generalization, and has important value for realizing predictive maintenance, planning the supply and reserve of spare parts and avoiding accident loss possibly caused by unplanned shutdown.

Drawings

FIG. 1 is a diagram of the variation trend of the horizontal vibration acceleration data, kurtosis and root mean square of 5 bearing full life cycles

FIG. 2 is a diagram of horizontal vibration acceleration data, kurtosis and root mean square variation trends of the other 5 bearings in the whole life cycle

FIG. 3 is a diagram of the horizontal vibration acceleration data, kurtosis and root mean square variation trends of the other 5 bearings in the whole life cycle

FIG. 4 is a graph of horizontal vibration acceleration data, kurtosis and root mean square variation trends of the other 2 bearings in the full life cycle

FIG. 5 is a detailed flow chart of an embodiment

FIG. 6 is HI values for constructing bearings based on MVGCDBN model

FIG. 7 shows the failure prediction result (a) gamma distribution (B) normal distribution of the bearing B1-3

FIG. 8 shows the failure prediction result (a) gamma distribution (B) normal distribution of the bearing B2-3

FIG. 9 shows the failure prediction result (a) gamma distribution (B) normal distribution of the bearing B3-3

FIG. 10 is a comparison of different failure prediction methods

Detailed Description

The present invention will be described in detail below with reference to the accompanying drawings 1 to 10, tables 1 to 4, and specific examples. It is to be understood that the embodiments described are only a few embodiments of the present invention, and not all embodiments. All other embodiments, which can be obtained by a person skilled in the art without any inventive step based on the embodiments of the present invention, shall fall within the scope of protection of the present invention.

The invention discloses an early fault prediction method verified by an experimental standard data set of an accelerated life of a rolling bearing actually measured by a PRONOSTIA test bed. The standard data set collects data under three working conditions, the sampling frequency is 25.6kHz, the data are collected every 10 seconds, the sampling duration lasts 0.1s every time, and 2560 acceleration data points can be obtained every time of sampling. The data set comprises vibration monitoring signals of 17 rolling bearings which are obtained in accelerated degradation experiments and run to failure, working condition information is shown in table 1, and bearing samples are described in table 2:

TABLE 1 bearing condition table

TABLE 2 bearing sample description table

In the rolling bearing accelerated life test, the early failure time of each natural deterioration of the bearing performance is not clearly described. However, sensitive monitoring indicators can still detect significant changes in dynamic behavior to a large extent. Kurtosis and Root Mean Square (RMS) are two of the most commonly used statistical parameters that quantify a bearing from healthy to failed, and are often used as early failure predictors for the bearing. Calculating the kurtosis and the root mean square of each bearing, drawing a time domain parameter curve of each bearing, and drawing figures 1-4 are horizontal direction vibration acceleration data, kurtosis and root mean square change trend graphs of the whole life cycles of 17 bearings. As can be seen from fig. 1-4:

the kurtosis of B1-1 fluctuates greatly in the life cycle, and no obvious growth trend or change occurs; RMS presents an increasing trend over the full life cycle, but no sudden increase points. Neither kurtosis nor RMS is therefore able to determine the early point of failure of bearing B1-1. The kurtosis of B1-2 fluctuates greatly in the life cycle, and no obvious growth trend or change occurs; the RMS was in a plateau phase before point 827 and showed a marked increase and change after point 827. Therefore, the kurtosis could not be used to determine early failure of bearing B1-2, and RMS suggested that bearing B1-2 failed early at 827. The kurtosis of B1-3 fluctuates greatly in the life cycle, and no obvious growth trend or change occurs; RMS presented an increasing trend over the full life cycle, but without abrupt increase points, B1-3 was not judged to have failed early. The kurtosis of B1-4 was in a more steady state before point 1083, and showed a significant trend and change after point 1083; the RMS was in a relatively steady state before point 1090 and showed a significant trend and change after point 1090. Thus kurtosis believes that bearing B1-4 failed early at point 1083 and RMS believes that bearing B1-4 failed early at point 1090. The kurtosis of B1-5 was in a more stable stage before point 2412, and showed a clear increasing trend and change after point 2412; the RMS was in a more steady state before point 2450 and showed a significant trend and change after point 2450. Kurtosis and RMS therefore predict early failure of bearing B1-4 at points 2412 and 2450, respectively. The kurtosis and RMS of B1-6 were in a more plateau phase before point 1631, and showed a significant trend and change in growth after point 1631. Both kurtosis and RMS therefore suggest that bearing B1-6 failed early at point 1631. The kurtosis and RMS of B1-7 were at a more plateau before 2210, and showed a significant trend and change in growth after 2210. Both kurtosis and RMS therefore suggest that bearing B1-7 failed early at 2210.

The kurtosis of B2-1 fluctuates greatly in the life cycle, and no obvious growth trend or change occurs; the RMS tended to fluctuate in the early, middle and late stages, showing a significant trend and change after 875 points. The kurtosis therefore failed early in bearing B2-1, and RMS assumed that bearing B2-1 failed early at point 875. The kurtosis of B2-2 is in a relatively steady state before point 752 except for unstable data of early operation, and shows obvious growth trend and change after point 752; the RMS showed a significant trend and change after 752 points. Both kurtosis and RMS therefore suggest that bearing B2-2 failed early at point 752. The kurtosis of B2-3 fluctuates greatly in the whole life cycle, and no obvious growth trend or change occurs; the RMS was in a relatively smooth stage before point 1947 except for rough running data at early times, and showed a significant trend and change after point 1947, so the kurtosis failed to determine the early failure point of bearing B2-3, whereas the RMS considered that bearing B2-3 failed early at point 1947. The kurtosis of B2-4 is in a relatively smooth stage before point 262, greatly fluctuates between points 262 and 410, is in a relatively smooth stage after point 410, and fluctuates a plurality of times from point 603 to the failure stage without obvious growth trend and change; the RMS was in a more steady state before 744 and showed a significant trend and change after 744, except for early running unstable data. The kurtosis therefore fails to determine the early failure point of bearing B2-4, and RMS assumes that an early failure occurred at 744. The kurtosis of B2-5, in addition to the earlier running jerky data, was at a more jerky stage before point 2269, showing a significant trend and change of growth after point 2269; the RMS fluctuated very well before point 530 and then was relatively smooth but showed no significant growth trend or change. The kurtosis therefore suggests that bearing B2-5 failed early at point 2269, and the RMS was unable to determine the early failure point of bearing B12. The kurtosis of B2-6 is in a relatively steady stage before 688 except for unstable data in early operation, and shows obvious growth trend and change after 688; the RMS was in a more steady state before point 688, with a marked increase and variation after point 688. Thus, both kurtosis and RMS suggest that bearing B2-6 failed early at 688. The kurtosis of B2-7 has many abrupt increase points in the life cycle, but has no obvious growth trend and change; the RMS was in a relatively steady state before 226 and showed a significant trend and change after 226. Therefore, the kurtosis does not determine the early failure point of bearing B2-7, and RMS believes that bearing B2-7 failed early at 226.

The kurtosis and RMS of B3-1 were in a more plateau phase before point 491, and showed a marked trend and change in growth after point 491. Both kurtosis and RMS therefore suggest that bearing B3-1 failed early at point 491. The kurtosis of B3-2 fluctuates at many places, and there is a clear increasing trend at point 1433; the RMS was in a more steady state before point 1614 and showed a significant increase and change after point 1614. Thus, kurtosis considers bearing B3-2 to fail early at 1433, and RMS considers bearing B3-2 to fail early at 1614. The kurtosis of B3-3 has a plurality of sudden increase points in the life cycle, but no obvious growth trend and change appear; the RMS was in a relatively steady state before point 310 and showed a significant trend and change after point 310. Thus, the kurtosis does not determine the early failure point of bearing B3-3, and the RMS assumes that bearing B3-3 failed early at point 310.

Due to the randomness of the degradation process, there is some difference in the full-life vibration signals of the multiple bearings. In order to avoid the contingency and the particularity of the results, the first 2 bearings (namely the bearings numbered 1-1, 1-2, 2-1, 2-2, 3-1 and 3-2) under three different working conditions are selected as a training set, and the rest bearings are selected as a testing set.

The specific flowchart of an embodiment of the early failure prediction method provided by the invention is shown in fig. 5, and mainly comprises the following steps:

s1: time-frequency domain fault feature extraction

S11: the STFT was used to obtain the time-frequency domain information (i.e., the case where the signal frequency varies with the time series, the instantaneous frequency at each time and its amplitude) in the non-stationary raw vibration signals of 17 rolling bearings. Let g_τν(t)Is a localized standard single window function centered around (tau, v), i.e. low-pass filtering, centered in the interval tau-deltag, tau + deltag in the time domain]Wherein e is^iνtIs the frequency shift factor.

g_τν(t)＝g(t-τ)e^iνtτ,ν∈R (1)

In the formula T_x(τ, ν) represents the fourier transformation of the signal within the window, x (t) is the original signal that needs to be transformed,it can be seen that a portion of the signal x (t) is intercepted by a window of width 2 Δ g centered on τ, and as τ moves, different portions of the signal will be sequentially captured in the window. STFT can be interpreted as: the signal is first frequency shifted by x (t) e^-iνtThen go through g_τν(t)Low-pass filtering is performed resulting in a time-domain representation at the frequency v. Wherein the content of the first and second substances,the hamming window function with a width of 25ms is used and the frame displacement is 10 ms.

S12: and carrying out gray processing on the obtained time-frequency image so as to facilitate subsequent model input, processing the three-channel image into a single-channel gray image, calling an RGB2gray (RGB) function in MATLAB, eliminating hue and saturation information, simultaneously keeping brightness, and converting the RGB image into a gray image. By calculating R, G a weighted sum of the B components: 0.2989 × R +0.5870 × G +0.1140 × B, the RGB values are converted into grayscale values, and the image pixels are processed to 32 × 32. And inputting the extracted time-frequency domain fault characteristics of the bearing into a GMVCDBN model to construct HI of the time-frequency domain fault characteristics.

S2: MVGCDBN-based health state characterization index construction

The standard CDBN is a new hierarchical generative model constructed from a number of convolution limited Boltzmann machines (CRBMs). The CRBM is an expansion of a standard Restricted Boltzmann Machines (RBM), adopts a local connection and weight sharing mode, considers the two-dimensional (2-D) structure and the periodic characteristic of input data, is suitable for spatial expansion of the data, overcomes the defect that all visual units are associated with all hidden units through different weight parameters, can directly acquire effective representation from original vibration data through alternate use of a convolution layer and a pooling layer, automatically extracts local features hidden in the data, combines and abstracts the local features into high-level features, and constructs HI through a full connection layer.

The standard CDBN is a two-layer structure consisting of a visible layer V and a hidden layer H. The visible layer is an N_V×N_VA binary cell matrix including K N_W×N_W(N_W＝N_V-N_H+1) convolution kernel W_k(i.e., a filter). The hidden layer is composed of K N_H×N_HThe feature mapping surface composition of the binary unit matrix marks the k-th feature mapping surface H^kHas a convolution kernel of W^kOffset is b_kAll visible layers share bias a.

If the visual unit is binary, the energy function can be defined as equation (3):

if the visual unit is real, the energy function can be defined as equation (4):

where v and h represent the state vectors of all visible and hidden layers, respectively. v. of_i,jIs an element of the ith row and jth column of the matrix v,is the element of the ith row and jth column of the kth hidden unit,is the element of the kth filter row r column s.Is thatThe matrix of the flip-flop of (c),denotes convolution operation, denotes the product of elements and then sums, i.e. a · B ═ trA^TB。

Mapping a surface H on a feature^kIn the above, the conditional activation probability of the hidden layer unit in the mth row and nth column can be expressed as formula (5):

similarly, the conditional activation probability of a visual layer unit can be expressed as formula (6):

in the method, a new one-dimensional CDBN model is constructed by adopting the Gaussian visual unit and the binary hidden layer, and the hidden layer output obtained by the model loses the information of the original data to a certain extent. The distribution of hidden layer nodes and visible layer nodes in a CRBM may be any exponential family distribution as indicated in a paper published by Welling et al in 2004 in the journal of the Neural Information Processing Systems. According to the method, the visual layer and the hidden layer both meet Gaussian distribution, the visual layer input is assumed to be Gaussian distribution by the model, input data can be fitted more accurately, and meanwhile richer information of an original sample can be extracted. In addition, when high-dimensional data such as images and videos is input, the input data needs to be vectorized first, which destroys the spatial internal structure of the sample. Therefore, by using the idea of Matrix variable limited Boltzmann Machine (MVRBM), a Matrix variable Gaussian Convolutional conditional Deep Belief Network (MVGCDBN) is proposed for the construction of the health status indicator. Assuming that each cell of the visible layer and the hidden layer obeys a variance of σ₁ ²And gamma₁ ²The energy function of the improved model is as follows, equation (7):

the hidden layer and visible layer unit conditional activation probabilities are respectively:

the purpose of model training is to determine the bestThe optimal parameter set, similar to the standard CDBN, updates the parameters of the model using a Contrast Divergence (CD) algorithm based on Gibbs sampling. Sampling hidden layer nodes according to the calculated hidden layer conditional probability, and taking the average value of the hidden layer nodes asVariance is σ₁ ²The Gaussian distribution is sampled, the computed probability value is not compared with a random number from 0 to 1 to obtain a value of 0 or 1, and the sampled data obtained by the method is not 0 or 1 but a real number sampled near the mean value. The sampling of the nodes of the visual layer is the same as the hidden layer sampling method.

Inputting: training samples (visual layer vectors) v, the number K of filters in a hidden layer, a learning rate, a regularization constant zeta, a sparsity coefficient rho and a maximum iteration number Q.

And (3) outputting: connection weight parameter W^kA visible layer bias a and a hidden layer bias b_k。

Random initialization: connection weight parameter W^kA visible layer bias a and a hidden layer bias b_kAre random smaller numbers.

The architecture of the model is determined through experiments, and the main parameters of the MVGCDBN model are shown in Table 3:

TABLE 3 parameters for the construction of a Rolling bearing HI

Note: the symbols of the first and second CRBMs are not distinguished

The results of the test bearing data set health curves for the MVGCDBN configuration are shown in fig. 6.

S3: self-learning fault early warning control limit

S31: fitting of distributions

(1) And normalizing the HI data in the normal stage.

(2) The shape parameters of the normal distribution and the gamma distribution of the HI data in the 'normal' stage are calculated by adopting maximum likelihood estimation.

If the random variable x (hi) obeys a probability distribution with a position parameter μ and a scale parameter σ, its probability density function is:

the random variable x is then said to obey the parameter μ, σ²Normal distribution of (a), generally noted as: x to N (mu, sigma)²)。

If the random variable x (HI data) obeys a probability distribution with a shape parameter α and an inverse scale parameter β, its probability density function is:

let x-gamma (alpha, beta), and letNamely, it isWherein alpha is more than 0 and beta is more than 0. The random variable X is then said to obey a gamma distribution with parameters α, β, commonly denoted as: x to G (. alpha.,. beta.).

S32: goodness of fit test

The fit distribution was examined using a K-S test to determine the distribution to which the HI obeys. Here, the significance level a₁＝0.05。

(1) Let f (x) be the frequency distribution, g (x) the theoretical distribution, assuming H₀: the data fit the theoretical distribution.

(2) Calculating D ═ max | f (x) -g (x) |, when the actual observed value D > D (n, alpha)₁) Then refuse H₀Otherwise, H is accepted₀It is assumed that, i.e. when h is 0,the original hypothesis is accepted.

S33: and calculating a threshold corresponding to the quantile on the HI to determine a normalized control limit, and then carrying out reverse normalization on the normalized control limit to obtain a fault early warning control limit, so as to realize early prediction of a fault point.

For reasons of space, only the calculation results of the early warning control limit of the first test data set (i.e. bearings B1-3, B2-3, B3-3) under each working condition are shown. To verify the superiority and generalization of the invented method, the invented method was compared with the usual statistical parameters (kurtosis and RMS) quantifying the bearing from health to failure, and other results were statistically and contrastively analyzed in table 4.

In the self-learning control limit construction process, the peak error generated in the acquisition process by external influence is generally 5%, so that the quantile on both sides is 0.05. As shown in fig. 7, is the result of the self-learning bearing B1-3 fault early warning control limit. And performing goodness-of-fit test on the gamma distribution and the normal distribution respectively by adopting a K-S test. The gamma distribution was first examined, and since h is 0 and p is 0.5486, the original assumption was accepted at a significance level of 0.05, i.e. the HI of B1-3 was considered to obey the gamma distribution with parameters α 5.6695 and β 0.034572. Normal distributions were re-examined and the original hypothesis was accepted at significance level 0.05, since h is 0 and p is 0.6321, considering the HI obeying parameters μ 0.19601 for B1-3, σ²Normal distribution of 0.082291. In summary, the HI compliance parameter for B1-3 was considered to be μ -0.19601, σ²Normal distribution of 0.082291. As shown in fig. 7(B), the normal data and the fault data are separated, the data of the previous 2110 groups are below the self-learning control limit (HI ═ 0.33) except the unstable data just started to operate, and the 2210 point is judged to be the early fault point, i.e. the bearing is operated for about 368 minutes (the total life of the bearing B1-3 is about 396 minutes). Based on the data of the example B1-3, as shown in FIGS. 1-4, the failure warning point cannot be determined by calculating the kurtosis and the root mean square.

As shown in fig. 8, is a self-learning shaftAnd carrying out fault early warning control limit result B2-3. The HI compliance parameter of B2-3 is determined to be mu-0.63004 and sigma through K-S test²Normal distribution of 0.012808. As shown in fig. 8(B), the normal data and the failure data have been separated, the previous 411 sets of data are below the self-learning control limit (0.84) except for the unstable data just started to operate, and the point 411 is judged as the early failure point, that is, the bearing is operated for 69 minutes (about 326 minutes for the total life of the bearing B2-3). Based on the data of the example B2-3, as shown in FIGS. 1-4, the RMS determination 1947 point is a failure warning point, and the kurtosis cannot determine the failure warning point.

As shown in fig. 9, is the result of the self-learning bearing B3-3 fault early warning control limit. The HI of B3-3 was considered by K-S test to comply with the gamma distribution with the parameters α 3.0943 and β 0.10952. As shown in fig. 9(a), the normal data and the failure data have been separated, the former 310 data is below the self-learning control limit (0.70) except for the unstable data just started to operate, and the early failure point is judged at 310, that is, the bearing is operated for about 52 minutes (about 72 minutes for the total life of the bearing B3-3). Based on the data of the example B3-3, as shown in FIGS. 1-4, the RMS certainty 310 is a failure warning point, and the kurtosis is unable to determine the failure warning point.

TABLE 4 comparison of results of different failure prediction methods

Note: "/" indicates that early failure cannot be predicted

As can be seen from table 4 and fig. 10, the overall performance of the proposed method in early failure prediction is better than kurtosis and RMS, which is particularly shown in:

(1) early failure of bearings B1-1, B1-3, both kurtosis and RMS could not be predicted; the invented method judges that the early failure points of B1-1 and B1-3 are respectively: 2433 points (i.e., bearings run to about 406 minutes), 2210 points (i.e., bearings run to about 368 minutes).

(2) For bearings B1-2, B2-1, B2-3, B2-4, B2-7 and B3-3, the kurtosis cannot judge the early failure points, and the RMS judges the early failure points to be respectively: 827 point (bearing run to about 138 minutes), 875 point (bearing run to about 146 minutes), 1947 point (bearing run to about 325 minutes), 744 point (bearing run to about 124 minutes), 226 point (bearing run to about 38 minutes), 310 point (bearing run to about 52 minutes); the invented method judges that their early failure points are respectively: 720 points (bearing run to about 120 minutes), 790 points (bearing run to about 132 minutes), 411 points (bearing run to about 69 minutes), 720 points (bearing run to about 120 minutes), 180 points (bearing run to about 30 minutes), 310 points (bearing run to about 52 minutes). Compared with the RMS method, the invented method is consistent with the RMS detection results of the B3-3 bearing. For the remaining bearings, the inventive method predicted early failure of bearings B1-2, B2-1, B2-4 and B2-7 18 minutes, 14 minutes, 4 minutes, 8 minutes earlier than RMS, respectively.

(3) For bearing B2-5, the RMS was unable to determine its early failure point, and both the inventive method and kurtosis were able to determine its early failure point to be 2269 (the bearing ran for about 378 minutes).

(4) For bearings B2-2 and B2-6, the kurtosis, RMS and method of the present invention simultaneously judged their early failure points to be 752 (bearing run to about 125 minutes) and 688 (bearing run to about 115 minutes), respectively.

(5) For bearings B1-4, B1-5, B1-6, B1-7, B3-1 and B3-2, the kurtosis judges that the early failure points are respectively as follows: 1083 points (bearing run to about 181 minutes), 2412 points (bearing run to about 402 minutes), 1631 points (bearing run to about 272 minutes), 2210 points (bearing run to about 368 minutes), 491 points (bearing run to about 82 minutes), 1433 points (bearing run to about 239 minutes); RMS judges that their early failure points are: 1090 point (bearing run to about 182 minutes), 2450 point (bearing run to about 408 minutes), 1631 point (bearing run to about 272 minutes), 2210 point (bearing run to about 368 minutes), 491 point (bearing run to about 82 minutes), 1614 point (bearing run to about 269 minutes); the method judges that the early fault points are respectively as follows: point 1030 (bearing run to about 172 minutes), point 2410 (bearing run to about 402 minutes), point 1600 (bearing run to about 267 minutes), point 2136 (bearing run to about 356 minutes), point 480 (bearing run to about 80 minutes), point 1330 (bearing run to about 223 minutes), all preceded the first two methods. Compared with kurtosis, for bearings (B1-4, B1-5, B1-6, B1-7, B3-1 and B3-2), the invented method can predict early failure thereof 9 minutes, 0 minutes, 5 minutes, 12 minutes, 2 minutes and 16 minutes ahead. Compared to RMS, the inventive method can predict early failure 10, 6, 5, 12, 2, and 16 minutes ahead.

The above-mentioned embodiments are merely illustrative of the principles and effects of the present invention, and some embodiments may be used, not restrictive; it should be noted that, for those skilled in the art, various changes and modifications can be made without departing from the inventive concept of the present invention, and these changes and modifications belong to the protection scope of the present invention.