METHOD FOR DIAGNOSING AND PREDICTING OPERATION CONDITIONS OF LARGE-SCALE EQUIPMENT BASED ON FEATURE FUSION AND CONVERSION

Abstract

A method for diagnosing and predicting operation conditions of large-scale equipment based on feature fusion and conversion, including: collecting a vibration signal of each operating condition of the equipment, and establishing an original vibration acceleration data set of the vibration signal; performing noise reduction on the original vibration acceleration data set, and calculating a time domain parameter; performing EMD on a de-noised vibration acceleration and calculating a frequency domain parameter; constructing a training sample data set through the time domain parameter and the frequency domain parameter; establishing a GBDT model, and inputting the training sample data set into the GBDT model; extracting a leaf node number set from a trained GBDT model; performing one-hot encoding on the leaf node number set to obtain a sparse matrix; and inputting the sparse matrix into a factorization machine to obtain a prediction result.

Claims

1. A method for diagnosing and predicting operation conditions of large-scale equipment based on feature fusion and conversion, comprising: (S1) collecting a vibration signal of individual operating conditions of the large-scale equipment, and establishing an original vibration acceleration data set of the vibration signal; (S2) performing noise reduction on the original vibration acceleration data set to obtain a de-noised vibration acceleration; (S3) calculating a time domain parameter based on the de-noised vibration acceleration; (S4) performing empirical mode decomposition (EMD) on the de-noised vibration acceleration to obtain a plurality of components; (S5) calculating a frequency domain parameter based on the plurality of components; (S6) constructing a training sample data set based on the time domain parameter and the frequency domain parameter; (S7) establishing a gradient boosting decision tree (GBDT) model, and inputting the training sample data set into the GBDT model for training to obtain a trained GBDT model; (S8) extracting a leaf node number set from the trained GBDT model; (S9) performing one-hot encoding on the leaf node number set to obtain a sparse matrix; and (S10) inputting the sparse matrix into a factorization machine to obtain a prediction result; wherein the step (S2) is performed through steps of: (S2-1) removing a direct-current (DC) component of a h.sup.th original vibration acceleration g.sub.h(i) in the original vibration acceleration data set according to formula (1) to obtain a h.sup.th first vibration acceleration g′.sub.h(i) with the DC component removed: $\begin{array}{cc}{g}_{h\left(i\right)}^{\prime }=g_{h\left(i\right)}+\frac{{.Math.}_{i=1}^n{g}_{h\left(i\right)}}{n};& \left(1\right)\end{array}$ wherein n is a data length of the h.sup.th original vibration acceleration; and (S2-2) removing a trend of the h.sup.th first vibration acceleration according to formulas (2)-(4) to obtain a h.sup.th de-noised vibration acceleration: $\begin{array}{cc}{g}_{h\left(i\right)}^{″}=g_{h\left(i\right)}^{\prime }-z;& \left(2\right)\end{array}$ $\begin{array}{cc}z={\alpha }^{\prime }+{\beta }^{\prime }t;& \left(3\right)\end{array}$ $\begin{array}{cc}\left({\alpha }^{\prime },{\beta }^{\prime }\right)=\arg \underset{\left(\alpha ,\beta \right)}{\min }\left[\stackrel{n}{\underset{i=1}{.Math.}}{\left(\alpha +\beta t-g_{h\left(i\right)}^{\prime }\right)}^2\right];& \left(4\right)\end{array}$ wherein z is the fitting trend linear function when a α=α′ and β=β′; α and β are parameters of a fitting trend linear function; (α′, β′) represents values of α and β corresponding to a minimum value of $\stackrel{n}{\underset{i=1}{.Math.}}{\left(\alpha +\beta t-g_{h\left(i\right)}^{\prime }\right)}^2;$ arg min[⋅] indicates minimization operation; and t represents time within a range of one sample; a sampling time t.sub.1 corresponding to the h.sup.th de-noised vibration acceleration is expressed as:
t.sub.1=e+f×i i≤n   (5); wherein e is an origin time of vibration acceleration; and f is a sampling frequency; in step (S3), the time domain parameter comprises an effective value, kurtosis, variance and crest factor of the de-noised vibration acceleration; the step (S5) is performed through steps of: (S5-1) obtaining an energy value P.sub.k of a k.sup.th component according to formula (6): $\begin{array}{cc}{P}_k=\underset{\xi =1}{\overset{N_1}{.Math.}}{h}_k^2\left(t^{*}\right)\times \mathrm{dT};& \left(6\right)\end{array}$ wherein N.sub.1 is the total number of discrete data points in the k.sup.th component h.sub.k(t*); dT represents a sampling frequency; and i.sub.1 represents a serial number of the discrete data points in the k.sup.th component h.sub.k(t*); (S5-2) sorting energy values of the plurality of components; selecting X components with the largest energy value; and subjecting the X components to fast Fourier transform to obtain X initial spectra; (S5-3) extracting a maximum amplitude P.sub.0.sup.x of each of the X initial spectra; (S5-4) setting a threshold; and selecting frequencies respectively corresponding to an amplitude higher than the threshold to establish X frequency-amplitude data sets; (S5-5) calculating an average amplitude P.sub.i.sup.x of the X frequency-amplitude data sets according to formula (7): $\begin{array}{cc}{P}_1^X=\frac{{.Math.}_{i_2=1}^{N_2}.Math.f_{i_2}.Math.}{N_2}{N}_2\le N_1;& \left(7\right)\end{array}$ wherein N.sub.2 is the number of selected frequencies; and f.sub.1.sub.2 is an amplitude corresponding to a i.sub.2.sup.th frequency; (S5-6) calculating a root mean square value P.sub.2.sup.X of the X frequency-amplitude data sets according to formula (8): $\begin{array}{cc}{P}_2^X=\sqrt{\frac{{.Math.}_{i_2=1}^{N_2}.Math.f_{i_2}^2.Math.}{4{\pi }^2{f}_{i_2}^2}};& \left(8\right)\end{array}$ wherein π represents a 180-degree angle; and (S5-7) selecting a frequency P.sub.3.sup.X corresponding to a maximum amplitude in the X frequency-amplitude data sets; in step (S7), the GBDT model is established through steps of: (S7-1) initializing a first weak learner G.sub.0(D) according to formula (9): $\begin{array}{cc}{G}_0\left(D\right)=\log \frac{P\left(y_k=0|D\right)}{1-P\left(y_k=0|D\right)}h\in \left\{1,.Math.,q\right\};& \left(9\right)\end{array}$ wherein P(y.sub.h=0|D) represents a probability that a h.sup.th tagged value y.sub.h in the training sample data set D is marked as 0 in an integrated feature real number set; and log is a base-10 logarithmic function; (S7-2) establishing S trees, and calculating a pseudo-residual rss.sub.s,h of a s.sup.th tree G′.sub.s(x.sub.h) by using a log-likelihood function, expressed as: $\begin{array}{cc}{\mathrm{YSS}}_{s,h}=y_h-\frac{1}{1+e^{-G_s^{\prime }\left(x_h\right)}}s\in \left\{0,1,.Math.S-1\right\};& \left(10\right)\end{array}$ wherein x.sub.h is a h.sup.th integrated feature real number in the training sample data set D; and e represents natural logarithm; (S7-3) fitting data (x.sub.h,rss.sub.s,h) using a classification and regression tree (Cart) to obtain a s.sup.th regression tree G.sub.s(x.sub.h) as a s.sup.th weak learner; (S7-4) calculating an optimal negative gradient fitted value c.sub.s,ks of each leaf node in the s.sup.th regression tree according to formula (11): $\begin{array}{cc}{c}_{s,\mathrm{ks}}=.Math.\frac{{.Math.}_{x_k\in R_{s,\mathrm{ks}}}{\mathrm{rss}}_{s,h}}{{.Math.}_{x_j\in R_{s,\mathrm{ks}}}\left(y_k-{\mathrm{rss}}_{s,h}\right)\left(1-y_h+{\mathrm{rss}}_{s,h}\right)}.Math.;& \left(11\right)\end{array}$ wherein R.sub.s,ks is a leaf node region of the s.sup.th regression tree; and ks represents a k.sup.th leaf node of the s.sup.th regression tree; (S7-5) updating a strong learner G .sub.s+1(x.sub.h) according to formula (12): $\begin{array}{cc}{G}_{s+1}\left(x_h\right)=G_s\left(x_h\right)+\underset{\mathrm{kz}=1}{\overset{K_s}{.Math.}}{c}_{z,\mathrm{ks}}I\left(x_k\in R_{s,\mathrm{ks}}\right);& \left(12\right)\end{array}$ wherein I(⋅) is a weighted count of x.sub.h within the leaf node region; and K.sub.s is the total number of leaf nodes of the s.sup.th regression tree; and (S7-6) obtaining a final strong learner G.sub.s+1′(x.sub.h) as the trained GBDT model based on an output c.sub.s+1,ks of an updated strong learner according to formula (13): $\begin{array}{cc}{G}_{s+1}^{\prime }\left(x_k\right)=G_0\left(x_h\right)+\underset{s=2}{\overset{S}{.Math.}}\underset{k=1}{\overset{K_s}{.Math.}}{c}_{s-1,k}I\left(x\in R_{s,\mathrm{ks}}\right).& \left(13\right)\end{array}$

2. The method of claim 1, wherein the step (S4) is performed through steps of: (S4-1) finding all extreme points of the de-noised vibration acceleration through a find_peaks function; wherein the find_peaks function is configured to find peaks; (S4-2) connecting local maximum points m.sub.up in all extreme points into an upper envelope n.sub.up , and connecting local minimum points m.sub.up in all the extreme points into a lower envelope n.sub.down through a cubic spline curve, expressed as follows:
n.sub.up=a.sub.up+b.sub.upm.sub.up+c.sub.upm.sub.up.sup.2+d.sub.upm.sub.up.sup.3 up≤n n.sub.down=a.sub.down+b.sub.downm.sub.down+.sub.downm.sub.down.sup.2+d.sub.downm.sub.down.sup.3 down≤n   (14); wherein a.sub.up, b.sub.up, c.sub.up and d.sub.up are coefficients of the upper envelope; and a.sub.down, b.sub.down, c.sub.down and d.sub.down are coefficients of the lower envelope; (S4-3) determining the coefficients of the upper envelope and the coefficients of the lower envelope based on all extreme points of the de-noised vibration acceleration through an interpolation algorithm; (S4-4) obtaining a component to be judged h.sub.1(t*)at moment t* according to formula (15): $\begin{array}{cc}{m}_1\left(t^{*}\right)=\frac{n_{\mathrm{sp}}+n_{\mathrm{down}}}{2}& \left(15\right)\end{array}$ $h_1\left(t^{*}\right)=g_h^{″}\left(t^{*}\right)-m_1\left(t^{*}\right);$ wherein m.sub.1(t*) is an average value of the upper envelope and the lower envelope at the moment t* ; and g′.sub.h(t*) is a time domain expression of a h.sup.th de-noised vibration acceleration; (S4-5) determining whether h.sup.1(t*) satisfies an intrinsic mode function; if yes, taking h.sub.1(t*) as a first component of the h.sup.th de-noised vibration acceleration, and proceeding to step (S4-6); otherwise, taking h.sub.1(t*) as an input signal, and repeating steps (S4-1)-(S4-5) until h.sub.1(t*) satisfies the intrinsic mode function, and proceeding to step (S4-6); (S4-6) obtaining a total set R.sub.1(t*) of remaining signals after a first decomposition according to formula (16):
R.sub.1(t*)=g′.sub.h(t*)−h(t*)   (16); and (S4-7) repeating steps (S4-1)-(S4-6) based on the total set of remaining signals to obtain the k.sup.th component h.sub.k(t*) .

3. The method of claim 1, wherein the step (S6) is performed through steps of: (S6-1) combining the frequency domain parameter into a frequency-domain feature real number set x.sub.2.sup.h={P.sub.0.sup.x, P.sub.1.sup.x, P.sub.2.sup.x, P.sub.3.sup.x}.sub.h, and forming the integrated feature real number set x.sub.h={x.sub.1.sup.h, x.sub.2.sup.h} in combination with the time domain parameter; wherein h is a h.sup.th data in the original vibration acceleration data set; and x.sub.1.sup.h is a time domain feature real number set; and (S6-2) tagging the integrated feature real number set to obtain the training sample data set D expressed as: D={(x.sub.1, y.sub.1), (x.sub.2, Y.sub.2), . . .(x.sub.q, Y.sub.q)}; wherein q represents the total amount of data in the original vibration acceleration data set; y.sub.q is a tagged value, and Y.sub.g ∈ {0,1}; y.sub.q=0 indicates that data in the integrated feature real number set is in a normal operating state; and y.sub.q=1 indicates that the data in the integrated feature real number set is in an abnormal operating state.

4. The method of claim 1, wherein the step (S8) is performed through steps of: extracting the leaf node number set T output by the trained GBDT model according to formula (16): $\begin{array}{cc}T={\left[\begin{array}{ccc}{l}_{\left(1,1\right)}& .Math.& l_{\left(1,t\right)}\\ .Math.& \ddots & .Math.\\ l_{\left(q,1\right)}& .Math.& l_{\left(q,t\right)}\end{array}\right]}_{q,S};& \left(16\right)\end{array}$ wherein the leaf node number set T has a size of q×S; l.sub.(q,s) represents a serial number of a leaf node on the s.sup.th regression tree where an output feature of a q.sup.th training sample after being trained; q is the total number of training samples; and S is the total number of regression trees.

5. The method of claim 4, wherein the step (S9) is performed through steps of: (S9-1) performing one-hot encoding conversion on the leaf node number set T; marking a leaf node where an output of a regression tree corresponding to a tree formed by each training sample in the training sample data set is located as 1, and marking remaining leaf nodes as 0 to obtain a sub-matrix; and (S9-2) subjecting each training sample to sub-matrix transformation followed by combining into the sparse matrix with a size of (S×K.sub.s, q) herein S is the total number of regression trees; K.sub.s is the total number of leaf nodes on the s.sup.th regression tree corresponding to the tree formed by each training sample; and q is the total number of training samples.

6. The method of claim 5, wherein the step (S10) is performed through steps of: (S10-1) based on decomposability of a real symmetric positive definite matrix, factorizing the sparse matrix to calculate a feature latent vector of the sparse matrix; (S10-2) obtaining a second-order polynomial model ŷ(x) formed by combination of a linear model and a cross term according to formula (17): $\begin{array}{cc}\hat{y}\left(x\right)={\omega }_0+\underset{u=1}{\overset{S\times K_s}{.Math.}}{w}_u{x}_u+\underset{v=1}{\overset{S\times K_z}{.Math.}}\underset{u=v=1}{\overset{S\times K_z}{.Math.}}.Math.V_u,V_v.Math.x_u{x}_v& \left(17\right)\end{array}$ $.Math.V_u,V_v.Math.=\underset{l}{\overset{d}{.Math.}}{v}_{u,l}\times v_{v,l}$ wherein ω.sub.0 is a bias term of the linear model; w.sub.u is a weight of a u.sup.th feature in the sparse matrix; x.sub.u is the u.sup.th feature in the sparse matrix; x.sub.v is a v.sup.th feature of the sparse matrix; V.sub.u, V.sub.v indicates conversion of a coefficient of the cross term into a latent vector inner product; V.sub.u is a feature latent vector of the u.sup.th feature of the sparse matrix; V.sub.v is a feature latent vector of the v.sup.th feature of the sparse matrix; V.sub.u,l is a l.sup.th element obtained from factorization of the u.sup.th feature; and V.sub.v,l is a l.sup.th element obtained from factorization of the v.sup.th feature; and (S10-3) optimizing the second-order polynomial model according to back propagation of a Hinge loss function in the factorization machine to obtain a prediction result.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0083] FIG. 1 is a flow chart of a method for diagnosing and predicting operation conditions of large-scale equipment based on feature fusion and conversion according to an embodiment of the present disclosure;

[0084] FIG. 2 is a receiver operating characteristic (ROC) curve of a gradient boosting decision tree (GBDT) and a factorization machine (FM) integrated training set A according to an embodiment of the present disclosure;

[0085] FIG. 3 shows a first component obtained from empirical mode decomposition (EMD) according to an embodiment of the present disclosure;

[0086] FIG. 4 shows a second component obtained from the EMD according to an embodiment of the present disclosure;

[0087] FIG. 5 shows a third component obtained from the EMD according to an embodiment of the present disclosure;

[0088] FIG. 6 shows a fourth component obtained from the EMD according to an embodiment of the present disclosure;

[0089] FIG. 7 shows a fifth component obtained from the EMD according to an embodiment of the present disclosure;

[0090] FIG. 8 shows a sixth component obtained from the EMD according to an embodiment of the present disclosure;

[0091] FIG. 9 shows a seventh component obtained from the EMD according to an embodiment of the present disclosure;

[0092] FIG. 10 shows an eighth component obtained from the EMD according to an embodiment of the present disclosure; and

[0093] FIG. 11 shows the rest of EMD according to an embodiment of the present disclosure.

DETAILED DESCRIPTION OF EMBODIMENTS

[0094] The disclosure will be described completely and clearly below with reference to the accompanying drawings and embodiments to make the object, technical solutions, and beneficial effects of the present disclosure clearer. Obviously, provided below are merely some embodiments of the disclosure, which are not intended to limit the disclosure. It should be understood that all other embodiments obtained by those skilled in the art based on the content disclosed herein without paying any creative effort should fall within the scope of the present disclosure.

[0095] Illustrated in FIG. 1 is a method for diagnosing and predicting operation conditions of large-scale equipment based on feature fusion and conversion, which is performed through the following steps.

[0096] (S1) A vibration signal of individual operating conditions of the large-scale equipment is collected, and an original vibration acceleration data set of the vibration signal is established.

[0097] (S2) Noise reduction is performed on the original vibration acceleration data set to obtain a de-noised vibration acceleration.

[0098] (S3) A time domain parameter based on the de-noised vibration acceleration is calculated.

[0099] (S4) Empirical mode decomposition (EMD) is performed on the de-noised vibration acceleration to obtain a plurality of components.

[0100] (S5) A frequency domain parameter based on the plurality of components is calculated.

[0101] (S6) A training sample data set based on the time domain parameter and the frequency domain parameter is constructed.

[0102] (S7) A gradient boosting decision tree (GBDT) model is established, and the training sample data set is input into the GBDT model for training to obtain a trained GBDT model.

[0103] (S8) A leaf node number set is extracted from the trained GBDT model.

[0104] (S9) One-hot encoding is performed on the leaf node number set to obtain a sparse matrix.

[0105] (S10) The sparse matrix is input into a factorization machine to obtain a prediction result.

[0106] The step (S2) is performed through the following steps.

[0107] (S2-1) A direct-current (DC) component of a h.sup.th original vibration acceleration g.sub.h(i) in the original vibration acceleration data set is removed according to formula (1) to obtain a h.sup.th first vibration acceleration g′.sub.h(i) with the DC component removed:

[00015]$\begin{array}{cc}{g}_{h\left(i\right)}^{\prime }=g_{h\left(i\right)}+\frac{{.Math.}_{i=1}^n{g}_{h\left(i\right)}}{n};& \left(1\right)\end{array}$

[0108] where n is a data length of the h.sup.th original vibration acceleration.

[0109] (S2-2) A trend of the first vibration acceleration is removed according to formulas (2)-(4) to obtain a h.sup.th de-noised vibration acceleration:

[00016]$\begin{array}{cc}{g}_{h\left(i\right)}^{″}=g_{h\left(i\right)}^{\prime }-z;& \left(2\right)\end{array}$$\begin{array}{cc}z={\alpha }^{\prime }+{\beta }^{\prime }t;\mathrm{and}& \left(3\right)\end{array}$$\begin{array}{cc}\left({\alpha }^{\prime },{\beta }^{\prime }\right)=\arg \underset{\left(\alpha ,\beta \right)}{\min }\left[\underset{i=1}{\overset{n}{.Math.}}{\left(\alpha +\beta t-g_{h\left(i\right)}^{\prime }\right)}^2\right];& \left(4\right)\end{array}$

[0110] where z is the fitting trend linear function when α=α′ and β=β′; α and β are parameters of a fitting trend linear function; (α′, β′) represents values of α and β corresponding to a minimum value of

[00017]$\underset{i=1}{\overset{n}{.Math.}}{\left(\alpha +\beta t-g_{h\left(i\right)}^{\prime }\right)}^2;$

arg min[⋅] indicates minimization operation; and t represents time within a range of one sample.

[0111] A sampling time t.sub.i corresponding to the h.sup.th de-noised vibration acceleration is expressed as:

t.sub.i=e+f×i i≤n   (5)

[0112] where e is an origin time of vibration acceleration; and f is a sampling frequency.

[0113] In step (S3), the time domain parameter includes an effective value, kurtosis, variance and crest factor of the de-noised vibration acceleration.

[0114] The step (S4) is performed through the following steps.

[0115] (S4-1) All extreme points of the de-noised vibration acceleration are found through a find_peaks function. The find_peaks function is configured to find peaks.

[0116] (S4-2) Local maximum points m.sub.up in all extreme points are connected into an upper envelope n.sub.up , and local minimum points m.sub.down in all the extreme points are connected into a lower envelope n.sub.down through a cubic spline curve, expressed as follows:

n.sub.up=a.sub.up+b.sub.upm.sub.up+c.sub.upm.sub.up.sup.2+d.sub.upm.sub.up.sup.3d up≤n n.sub.down=a.sub.down+b.sub.downm.sub.down+c.sub.downm.sub.down.sup.2+d.sub.downm.sub.down.sup.3 down≤n   (14)

[0117] where a.sub.up, b.sub.up, c.sub.up and d.sub.up, are coefficients of the upper envelope; and a .sub.down, b.sub.down, C.sub.down and d.sub.down are coefficients of the lower envelope.

[0118] (S4-3) The coefficients of the upper envelope and the coefficients of the lower envelope based on all extreme points of the de-noised vibration acceleration are determined through an interpolation algorithm.

[0119] (S4-4) A component to be judged h.sub.1(t*) at a moment t* is obtained according to formula (15):

[00018]$\begin{array}{cc}{m}_1\left(t^{*}\right)=\frac{n_{\mathrm{up}}+n_{\mathrm{down}}}{2}& \left(15\right)\end{array}$$h_1\left(t^{*}\right)=g_h^{″}\left(t^{*}\right)-m_1\left(t^{*}\right);$

[0120] where m.sub.1(t*) is an average value of the upper envelope and the lower envelope at the moment t*; and g″.sub.h(t*) is a time domain expression of a h.sup.th de-noised vibration acceleration.

[0121] (S4-5) Whether h.sub.1(t*) satisfies an intrinsic mode function is determined.

[0122] If yes, h.sub.1(t*) is taken as a first component of the h.sup.th de-noised vibration acceleration, and step (S4-6) is performed.

[0123] Otherwise, h.sup.1t(t*) is taken as an input signal, and steps (S4-1)-(S4-5) are repeated until h.sub.1(t*) satisfies the intrinsic mode function, and step (S4-6) is performed.

[0124] (S4-6) A total set R.sub.1(t*) of remaining signals after a first decomposition is obtained according to formula (16):

R.sub.1(t*)=g″.sub.h(t*)   (16)

[0125] (S4-7) Steps (54-1)-(54-6) are repeated based on the total set of remaining signals to obtain the k.sup.th component h.sub.k(t*).

[0126] The step (S5) is performed through the following steps.

[0127] (S5-1) An energy value P.sub.k of a k.sup.th component is obtained according to formula (6):

[00019]$\begin{array}{cc}{P}_k=\underset{i_1=1}{\overset{N_1}{.Math.}}{h}_k^2\left(t^{*}\right)\times \mathrm{dT};& \left(6\right)\end{array}$

[0128] where N.sub.1 is the total number of discrete data points in the k.sup.th component h.sub.k(t*); dT represents a sampling frequency; and i.sub.1 represents a serial number of the discrete data points in the k.sup.th component h.sub.k(t*).

[0129] (S5-2) Energy values of the plurality of components are sorted. X components with the largest energy value are selected. The X components are subjected to fast Fourier transform to obtain X initial spectra.

[0130] (S5-3) A maximum amplitude P.sub.0.sup.x of each of the X initial spectra is extracted.

[0131] (S5-4) A threshold is set. Frequencies respectively corresponding to an amplitude higher than the threshold is selected to establish X frequency-amplitude data sets.

[0132] (S5-5) An average amplitude P.sub.1.sup.X of the X frequency-amplitude data sets is calculated according to formula (7):

[00020]$\begin{array}{cc}{P}_1^X=\frac{{.Math.}_{t_2=1}^{N_2}.Math.f_{t_2}.Math.}{N_2}{N}_2\le N_1;& \left(7\right)\end{array}$

[0133] where N.sub.2 is the number of selected frequencies; and f.sub.i.sub.2 is an amplitude corresponding to a i.sub.2.sup.th frequency.

[0134] (S5-6) A root mean square value P.sub.2.sup.X of the X frequency-amplitude data sets is calculated according to formula (8):

[00021]$\begin{array}{cc}{P}_2^X=\sqrt{\frac{{.Math.}_{t_2=1}^{N_2}.Math.f_{i_2}^2.Math.}{4{\pi }^2{f}_{i_2}^2}};& \left(8\right)\end{array}$

[0135] where π represents a 180-degree angle. (S5-7) A frequency P.sub.3.sup.x corresponding to a maximum amplitude in the X frequency-amplitude data sets is selected.

[0136] The step (S6) is performed through the following steps.

[0137] (S6-1) Frequency domain parameters are combined into a frequency-domain feature real number set x.sub.2.sup.h={P.sub.0.sup.x, P.sub.1.sup.x, P.sub.2.sup.x, P.sub.3.sup.x}.sub.h, and an integrated feature real number set x.sub.h={x.sub.1.sup.h, x.sub.2.sup.h} is formed in combination with the time domain parameters; where h is a h.sup.th data in the original vibration acceleration data set; and x.sub.1.sup.h is a time domain feature real number set.

[0138] (S6-2) The integrated feature real number set is tagged to obtain the training sample data set D expressed as: D={(x.sub.1,y.sub.1), (x.sub.2, y.sub.2), . . . (x.sub.q, y.sub.q)}.

[0139] where q represents the total amount of data in the original vibration acceleration data set; y.sub.q is a tagged value, and y.sub.q ∈ {0,1}; y.sub.q=0 indicates that data in the integrated feature real number set is in a normal operating state; and y.sub.q=1 indicates that the data in the integrated feature real number set is in an abnormal operating state.

[0140] In step (S7), the GBDT model is established through the following steps. (S7-1) A first weak learner G.sub.0(D) is initialized according to formula (9):

[00022]$\begin{array}{cc}{G}_0\left(D\right)=\log \frac{P\left(y_h=0|D\right)}{1-P\left(y_k=0|D\right)}h\in \left\{1,.Math.,q\right\};& \left(9\right)\end{array}$

[0141] where P(y.sub.h=0|D) represents a probability that a h.sup.th tagged value y.sub.h in the training sample data set D is marked as 0 in the integrated feature real number set; and log is a base-10 logarithmic function.

[0142] (S7-2) S trees are established. A pseudo-residual rss.sub.s,h of a s.sup.th tree G′.sub.s(x.sub.h) is calculated by using a log-likelihood function, expressed as:

[00023]$\begin{array}{cc}{\mathrm{rss}}_{z,h}=y_h-\frac{1}{1+e^{-G_z^{\prime }\left(x_k\right)}}s\in \left\{0,1,.Math.S-1\right\}.& \left(10\right)\end{array}$

[0143] x.sub.h is a h.sup.th integrated feature real number in the training sample data set D; and e represents natural logarithm.

[0144] (S7-3) Data (x.sub.h,rss.sub.s,h) is fit by using a classification and regression tree (Cart) to obtain a silt regression tree G.sub.s(x.sub.h) as a s.sub.th weak learner.

[0145] (S7-4) An optimal negative gradient fitted value c.sub.s,ks of each leaf node in the s.sup.th regression tree is calculated according to formula (11):

[00024]$\begin{array}{cc}{c}_{s,\mathrm{kz}}=.Math.\frac{{.Math.}_{x_h\in R_z,\mathrm{ks}}{\mathrm{rss}}_{s,h}}{{.Math.}_{x_s\in R_{z,\mathrm{ks}}}\left(y_h-{\mathrm{rss}}_{s,h}\right)\left(1-y_h+{\mathrm{rss}}_{s,h}\right)}.Math.;& \left(11\right)\end{array}$

[0146] where R.sub.s,ks is a leaf node region of the s.sup.th regression tree; and ks represents a k.sup.th leaf node of the s.sup.th regression tree.

[0147] (S7-5) A strong learner G.sub.s+1(x.sub.h) is updated according to formula (12):

[00025]$\begin{array}{cc}{G}_{z-1}\left(x_h\right)=G_s\left(x_h\right)+\underset{\mathrm{kz}=1}{\overset{K_s}{.Math.}}{c}_{z,\mathrm{ks}}I\left(x_h\in R_{s,\mathrm{ks}}\right);& \left(12\right)\end{array}$

[0148] where I(⋅) is a weighted count of x.sub.h within the leaf node region; and K.sub.s is the total number of leaf nodes of the s.sup.th regression tree.

[0149] (S7-6) A final strong learner G.sub.s+1′(x.sub.h) as the trained GBDT model based on an output c.sub.s+1, ks of an updated strong learner is obtained according to formula (13):

[00026]$\begin{array}{cc}{G}_{s+1}^{\prime }\left(x_k\right)=G_0\left(x_h\right)+\stackrel{S}{\underset{s=2}{.Math.}}\underset{k=1}{\overset{K_s}{.Math.}}{c}_{s-1,\mathrm{ks}}I\left(x\in R_{s,\mathrm{ks}}\right).& \left(13\right)\end{array}$

[0150] The step (S8) is performed through the following steps.

[0151] The leaf node number set T output by the trained GBDT model is extracted according to formula (16):

[00027]$\begin{array}{cc}T={\left[\begin{array}{ccc}{l}_{\left(1,1\right)}& .Math.& l_{\left(1,s\right)}\\ .Math.& \ddots & .Math.\\ l_{\left(q,1\right)}& .Math.& l_{\left(q,s\right)}\end{array}\right]}_{q,S};& \left(16\right)\end{array}$

[0152] where the leaf node number set T has a size of q×S ; l.sub.(q,s) represents a serial number of a leaf node on the s.sup.th regression tree where an output feature of a q.sup.th training sample after being trained; q is the total number of training samples; and S is the total number of regression trees.

[0153] The step (S9) is performed through the following steps. 6p (S9-1) One-hot encoding conversion is performed on the leaf node number set T A leaf node where an output of a regression tree corresponding to a tree formed by each training sample in the training sample data set is located is marked as 1, and remaining leaf nodes are marked as 0 to obtain a sub-matrix.

[0154] (S9-2) Each training sample is subjected to sub-matrix transformation, and combined into the sparse matrix with a size of (S×K.sub.s, q) where S is the total number of regression trees; K.sub.s is the total number of leaf nodes on the s.sup.th regression tree corresponding to the tree formed by each training sample; and q is the total number of training samples.

[0155] The step (S10) is performed through steps of:

[0156] (S10-1) Based on decomposability of a real symmetric positive definite matrix, the sparse matrix is factorized to calculate a feature latent vector of the sparse matrix.

[0157] (S10-2) A second-order polynomial model ŷ(x) formed by combination of a linear model and a cross term is obtained according to formula (17):

[00028]$\begin{array}{cc}\hat{y}\left(x\right)={\omega }_0+\underset{u=1}{\overset{S\times K_s}{.Math.}}{w}_u{x}_u+\underset{v=1}{\overset{S\times K_s}{.Math.}}\underset{u=v=1}{\overset{S\times K_s}{.Math.}}.Math.V_u,V_v.Math.x_u{x}_v& \left(17\right)\end{array}$$.Math.V_u,V_v.Math.=\stackrel{d}{\underset{l}{.Math.}}{v}_{u,i}\times v_{v,i};$

[0158] where ω.sub.0 is a bias term of the linear model; w.sub.u is a weight of a u.sup.th feature in the sparse matrix; x.sub.u is the u.sup.th feature in the sparse matrix; x.sub.v is a with feature of the sparse matrix; V.sub.u, V.sub.v indicates a conversion of a coefficient of the cross term into a latent vector inner product; V.sub.u is a feature latent vector of the u.sup.th feature of the sparse matrix; V.sub.v is a feature latent vector of the with feature of the sparse matrix; V.sub.u,l is a l.sup.th element obtained from factorization of the u.sup.th feature; and V.sub.,l is a l.sup.th element obtained from factorization of the v.sup.th feature.

[0159] (S10-3) The second-order polynomial model is optimized according to back propagation of a Hinge loss function in the factorization machine to obtain a prediction result.

[0160] A total set R.sub.k(t*) of remaining signals after a k.sup.th decomposition in step (S4-7) is expressed as:

R.sub.k(t*)=R.sub.k-1(t*)−h.sub.k(t*)   (18);

[0161] where h.sub.k(t*) is the k.sup.th component; R.sub.k−1(t*) is a total set of remaining signals after the (k−1).sup.th decomposition.

[0162] A value of the weighted count I(x.sub.h ∈ R.sub.s,ks) in step (S7-5) is expressed as:

[00029]$I\left(x_h\in R_{s,\mathrm{ks}}\right)=\{\begin{array}{cc}1,& \mathrm{if}{x}_h\mathrm{is}\mathrm{in}a{\mathrm{ks}}^{\mathrm{th}}\mathrm{leaf}\mathrm{node}\mathrm{on}\mathrm{the}{s}^{\mathrm{th}}\mathrm{regression}\mathrm{tree}\\ 0,& \mathrm{if}{x}_h\mathrm{is}\mathrm{not}\mathrm{in}a{\mathrm{ks}}^{\mathrm{th}}\mathrm{leaf}\mathrm{node}\mathrm{on}\mathrm{the}{s}^{\mathrm{th}}\mathrm{regression}\mathrm{tree}\end{array}.$

[0163] In an embodiment, a high-frequency cut-off frequency of a low-pass filter is set to 40 Hz. The original vibration acceleration data set is subjected to low-pass filtering to remove the DC component. The threshold in step (S5-4) is set to 0.2 times of the maximum amplitude. The integrated feature real number set is divided into a training set and a test set in a ratio of 7:3 for model training, and the training set is divided into a GBDT model training set B and a GBDT+factorization machine (FM) integrated training set A in a ratio of 1:1. A learning rate is 0.01. A sub-tree depth is 3. The maximum number of weak learners is 40. The GBDT model training set B is used in step (S7) for training, and the GBDT+FM integrated training set A is used in step (S8) and later to build the data model. In step (S10-3), the factorization machine performs 300 iterations, and a bivariate parameter dimension is 50.

[0164] As shown in FIG. 2, provided herein is a receiver operating characteristic (ROC) curve of the GBDT+FM integrated training set A used in step (S10); where TPR is a proportion of correctly identified positive data in the total positive data; and FPR is a proportion of predicting negative data as positive data. The test set is input into the model to be tested and predicted, and the results are analyzed. AUC (area under the ROC curve) is 0.93, which shows that the method in the present disclosure has a high prediction accuracy.

[0165] As shown in FIGS. 3-11, provided herein are the component diagrams obtained by decomposition in step (S4), where time series is the time series; IMF1-8 are the components obtained from the decomposition in an order of frequency from high to low; and Res is the remaining signals obtained from the decomposition.

[0166] Regarding the method provided herein, time domain features and frequency domain features are combined, so that the physical characteristics of original signals can be directly reflected, and the initial feature fusion can also be fully utilized to obtain reliable signal data.

[0167] Empirical mode decomposition (EMD) is adopted to de-noise and decompose the signal after the frequency domain conversion, and the selected frequency domain feature parameters can effectively represent the characteristics of the original signals and further relieve the noise interference while retaining the characteristics of original signals, so as to prepare for the next algorithm processing.

[0168] The diagnostic model involves the combination of a gradient boosting decision tree (GBDT) and a factorization machine, where by means of the GBDT, effective features can be found to obtain a local optimal solution; based on the ultra-high dimensional sparse characteristic of the factorization machine, the GBDT model is trained to obtain a global optimal solution. The combination of the two allows a higher fault resolution without manual feature selection.