© AFM, EDP Sciences 2019

1 Introduction

Roller sleeve as an important component is widely used in crusher, the status of its running performance directly affects the health of the whole equipment [1]. However, due to the complex working load, dust and other harsh working conditions, the service life of crusher roller sleeve is not long, and the accurate life of high speed spindle in all kinds of crushers is only thousands of hours. Once the working hours exceed the service life limit, the operation precision of roller will drop sharply and further cause that the machine cannot work properly. So it is very important to improve the reliability, safety and work efficiency of the crusher roller sleeve by means of prognostics and health management (PHM). It is an important part of PHM to predict the RUL of equipment and evaluate the performance of devices [2].

It is critical to establish an appropriate model in the life prediction process. In brief, condition-based monitoring is becoming more and more significant, especially in the RUL prediction. Vibration signal on-line monitoring is one of the most effective methods to monitor the state of crusher roller sleeve health condition (SOH) [3]. The RUL prediction of crusher roller sleeve based on vibration monitoring data is divided into two steps: firstly, construct an indicator to accurately assess the performance degradation of crusher roll sleeve; subsequently, establish an effective model to predict the RUL of crusher roller sleeve.

How to establish an appropriate model under partially observable state is the key to predict the RUL accurately, and it is also the urgent demand for the industrial production. SVR-PF is such a machine learning algorithm to make classification and prediction under small samples [4]. This method based on the statistics theory has been successfully applied to the prediction in the financial, electric and other systems [5,6].

In this paper, a method using acceleration and temperature data is proposed firstly to solve the challenge of low RUL prediction precision based on single data source. However, the noise and vibration interferences caused by other mechanical and systems may severely obscure the roller sleeve signal collected from sensors and make it very challenging to reliably detect the effective components. For the above reasons, a variety of signal analysis methods have been proposed by researchers, such as time-domain, frequency-domain and time-frequency technique. Wavelet analysis is such a widely accepted approach. Then it selects the sensitive features of two signals as input and constructs the SVR-PF model to solve the problem that is difficult to predict with finite state data. The proposed method is evaluated using experimental data respectively. Finally, after assessing the prediction result errors, the conclusion is given in Table 5.

2 Theory introduction of D-S data fusion and SVR-PF

2.1 Theory introduction of D-S data fusion

D-S theory fuses data from different sources through the basic probability assignment (BPA) function, and analyses the belief of all the possible propositions in the identification framework, so as to achieve the goal of data fusion [7–10].

To set the identification framework consists of evidence B and C, m₁ and m₂ are two BPA functions in the identification framework, m_1,2 is the fused BPA function, then the D-S data fusion can be expressed as:$$m_{1,2}(\mathrm{Ø})=0$$(1) $$m_{1,2}(A)=\frac{1}{1-K}{\displaystyle {\displaystyle \sum }_{B\cap C=A}}{m}_1(B)m_2(C)$$(2)where K is the degree of conflict between the evidences B and C $$K={\displaystyle {\displaystyle \sum }_{B\cap C=\mathrm{\varnothing }}}{m}_1(B)m_2(C).$$(3)

D-S data fusion theory combines with the same view for different sources of the same problem and eliminates the all conflicting views at once, so that a more reliable fused posterior BPA function can be obtained.

The RUL prediction of crusher roll sleeve based on D-S data fusion has two data sources: (1) the RUL prediction based on temperature data; (2) the RUL prediction based on acceleration data. In this paper, the proposed prediction method based on D-S data fusion and SVR-PF fuses the results of two prediction methods to gain the fused RUL prediction result.

The whole identified framework is set as Ω$$\mathrm{\Omega }=\{T,a\}.$$(4)

Because there's no intersection between T and a, prediction by acceleration data and prediction by temperature data are independent events, then the power set can be expressed as 2^Ω.$$2^{\Omega }=\{\mathrm{\varnothing },\{T\},\{a\},\{T\hspace{0.17em}\cup \hspace{0.17em}a\left\}\right\}.$$(5)

The meaning of all the propositions in the power set 2^Ω is explained as follows.

• {T} represents the RUL prediction credibility obtained by temperature data;

• {a} represents the RUL prediction credibility obtained by acceleration data;

•  {T ∪ a } represents the RUL prediction credibility obtained by acceleration or temperature data.

Meanwhile, the BPA functions m₁ and m₂ defined in the power set 2^Ω mean that:

• m₁ represents the prediction credibility distribution obtained by temperature data in the power set 2^Ω;

• m₂ represents the prediction credibility distribution obtained by acceleration data in the power set 2^Ω.

The combination of the BPA function based on data fusion is shown in Table 1.

From Table 1, the posterior BPA function based on data fusion can be expressed as:$$m\left(T\right)=\frac{1}{1-K}{\displaystyle {\displaystyle \sum }_{B\cap C=T}}{m}_1(B)m_2(C)=\frac{b}{1-K}$$(6) $$m(a)=\frac{1}{1-K}{\displaystyle {\displaystyle \sum }_{B\cap C=a}}{m}_1(B)m_2(C)=\frac{c}{1-K}$$(7)where$$b=m_1(T)m_2(T)+m_1(T\hspace{0.17em}\cup \hspace{0.17em}a)m_2(T)+m_1(T)m_2(T\hspace{0.17em}\cup \hspace{0.17em}a)$$(8) $$c=m_1(a)m_2(a)+m_1(T\hspace{0.17em}\cup \hspace{0.17em}a)m_2(a)+m_1(a)m_2(T\hspace{0.17em}\cup \hspace{0.17em}a)$$(9) $$K={\displaystyle {\displaystyle \sum }_{B\cap C=\varnothing }}{m}_1(B)m_2(C)=m_1(a)m_2(T)+m_1(T)m_2(a).$$(10)

That's the proposed RUL prediction method for crusher roller sleeve, it uses the posterior fusion BPA function and combines the prediction results of two prediction methods to gain a more accurate prediction result.

2.2 Basic theory of SVR-PF

2.2.1 Basic theory of particle filter

On the basis of the recursive Bayesian estimation [11–13], particle filter becomes a universal algorithm drawing samples from posterior distributions and assigns weights to all the particles by using the Monte Carlo method [13–16].

Particle filter has more excellent performance on nonlinear and non-Gaussian system than Kalman filter which only has good performance on liner and Gaussian system [17].

The particle filter system state space model can be described as:$$\{\begin{array}{l}{x}_k=f(x_{k-1},v_{k-1})\\ z_k=h(x_k,n_k)\end{array}$$(11)where x_k is the system state, z_k is either the system output or the measurement, v_k−1 is the system noise, and n_k is the measurement noise.

We assume that the prior distribution $p(x_{0:k-1}^i|z_{1:k-1})$ of system is known and N samples from the posterior distribution of system (11) have drawn. The posterior distribution can be approximately described as:$$p(x_{0:k}|z_{1:k})\approx {\displaystyle \sum _{i=1}^N}{w}_k^i\delta (x_{0:k}-x_{0:k}^i)$$(12)where $\left\{x_k^i\right\}$ is the sample, $\left\{w_k^i\right\}$ is the sample weight, which have ${\sum }_i^N{w}_k^i=1$. The higher is the weight, the higher is the sample probability. δ(⋅) represents the Dirac-Delta function.

In order to solve the problem that is very difficult to sample directly from a posterior distribution, a good deal of the problem is the importance sampling technique. It can draw samples directly from the importance distribution. The importance distribution can be described as:$$q(x_{0:k}|z_{1:k})\approx {\displaystyle \sum _{i=1}^N}\delta (x_{0:k}-x_{0:k}^i).$$(13)Plugging the importance distribution (13) into (12), then the weight can be updated:$$w_k^i=\frac{p(z_k|x_k^i)p(x_k^i|x_{k-1}^i)p(x_{0:k-1}^i|z_{1:k-1})}{q(x_k^i|x_{0:k-1}^i,z_{1:k})q(x_{0:k-1}^i|z_{1:k-1})}=w_{k-1}^i\frac{p(z_k|x_k^i)p(x_k^i|x_{k-1}^i)}{q(x_k^i|x_{0:k-1}^i,z_{1:k})}$$(14)where $p(z_k|x_k^i)$ is the likelihood function,$p(x_k^i|x_{k-1}^i)$ is the state transfer distribution. If system (11) subjects to the Markov process, the weight update equation (14) can be reduced to:$$w_k^i=w_{k-1}^i\frac{p(z_k|x_k^i)p(x_k^i|x_{k-1}^i)}{q(x_k^i|x_{k-1}^i,z_k)}.$$(15)

We set state transfer distribution as the importance distribution:$$q(x_k^i|x_{k-1}^i,z_k)=p(x_k^i|x_{k-1}^i).$$(16)

If the likelihood function $p(z_k|x_k^i)$ and the prior weights are used to update the new weights [15], the weight renew equation can be reduced to equation (17):$$w_k^i=w_{k-1}^{i}p(z_k|x_k^i).$$(17)

There is a wider problem of PF, which is known as degeneracy phenomenon. In order to avoid the problem, resampling is a suitable method. If the system iterates without resampling, the weight of some particles will tend to zero, and all efforts for the weights calculation become meaningless.

The standard method to avoid the degeneracy phenomenon is to renormalize the distribution by removing the small weight particles and duplicating the large weight particles. The weights of all the particles are set to 1/N (N is the number of particles). The resampling algorithm of the standard PF is shown above.$$N_{eff}=\frac{N}{1+\mathrm{\backslash rm\; var}(w_k^i)}\approx \frac{1}{{\displaystyle \sum _{i=1}^N}{(w_k^i)}^2}$$(18)where N_eff is the threshold of resampling.

2.2.2 Support vector regression-particle filter

The standard PF algorithm eliminates the small weight particles and duplicates the large weight particles to avoid the degeneracy phenomenon that would cause the loss of particle diversity. Which would make most particles aggregate around the larger weighted ones, so the degeneracy phenomenon still exists. In view of this problem, a new resampling algorithm known as SVR is introduced to rebuild a posterior distribution [18], which has an extremely fast learning speed and advantageous generalization capability. Moreover, SVR has a commendable performance in both classification and regression with a simple structure. Compared with other methods, SVR can avoid the degeneracy phenomenon and keep the diversity of particles in the case of limited samples. What's more, the training speed of SVR is much faster while obtaining better generalization. In view of these advantages, the SVR is selected in this paper to establish RUL prediction model. The application of the SVR is detailed in some studies [19,20].

The fundamental principle of SVR is known as an optimization problem expressed by a regularized functional with constraints [21], the form can be described as:$$\{\begin{array}{l}\Omega ={(f,f)}_H\\ s.t.{\displaystyle \underset{x}{\sup }}\left|F(x)-F_l(x)\right|={\displaystyle \underset{x}{\sup }}\left|F_l(x)-{\displaystyle \underset{-\infty }{\overset{x}{\int }}}f(t)dt\right|={\sigma }_l<ϵ\end{array}$$(19)where the regularized functional defined in Hilbert space and generated by σ_l is represented by Ω = (f, f)_H. The error between the distribution functions F(x) and their estimation F_l(x) is represented by σ_l. The constraint is represented by ϵ.

The estimated probability density function (PDF) of distribution F_l(x) is F(x). Only the points x_i(i = 1, 2,…, m) in the particle set should be considered, so equation (19) can be reduced to:$${\displaystyle \underset{i}{\max }}{\left|F_l(x)-{\displaystyle \underset{-\infty }{\overset{x}{\int }}}f(t)dt\right|}_{x=x_i}={\sigma }_l<ϵ.$$(20)If the PDF f(x) is described by kernel functions:$$f(x)={\displaystyle \sum _{i=1}^m}{\beta }_{i}K(x_i,x)$$(21)

Kernel function K(x_i , x) = ϕ _{^T  (xi)ϕ(x)} satisfies Mercer's condition. Then the regularized functional can be described as:$$\Omega (f)={(f,f)}_H={\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^m}{\beta }_i{\beta }_{j}K(x_i,x_j).$$(22)

The posterior distribution prediction can be described as an optimization problem with constraints:$$\{\begin{array}{l}\min w_p(\beta )={\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^m}{\beta }_i{\beta }_{j}K(x_i,x_j)\\ s.t.{\displaystyle \underset{i}{\max {\left|F_i(x)-{\displaystyle \sum _{j=1}^m}{\beta }_j{\displaystyle \underset{-\infty }{\overset{x}{\int }}}K(x_j,t)dt\right|}_{x=x_j}}}={\sigma }_l& .\end{array}$$(23)

Set y_i =F_i(x_i)_{, w =}[_{β1, β2,…, βm]  ^T} , $z_j(x)={\int }_{-\infty }^{x_i}K(x,t)dt$, z_i  = (z_i (x₁),z_i (x₂),…, z_i (x_m ))_{, ξi} and ${\xi }_i^{*}$ are non-negative slack variables, then equation (23) can be reduced to a quadratic programming problem:$$\{\begin{array}{l}\min J(w,{\xi }_i,{\xi }_i^{*})=\frac{1}{2}{w}^{T}w+C({\displaystyle \sum _{i=1}^m}{\xi }_i+{\displaystyle \sum _{i=1}^m}{\xi }_i^{*})\\ s.t.w^T{z}_i-y_i\le {\sigma }_l+{\xi }_i\\ y_i-w^T{z}_i\le {\sigma }_l+{\xi }_i^{*}\\ {\xi }_i,{\xi }_i^{*}\ge 0,i=1,2,\dots ,m\end{array}$$(24)where C is the penalty coefficient. By introducing Lagrange coefficients a_{i ,} $a_i^{*}$ to equation (24), we get:$$\{\begin{array}{l}\max w(a_i,a_i^{*})=-\frac{1}{2}{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^m}(a_i^{*}-a_i)(a_j^{*}-a_j)(z_i^T{z}_j)-{\sigma }_l{\displaystyle \sum _{i=1}^m}(a_i^{*}+a_i)\\ +{\displaystyle \sum _{i=1}^m}{y}_i(a_i^{*}-a_i)\\ s.t.{\displaystyle \sum _{i=1}^m}(a_i^{*}-a_i)=0,0\le a_i,a_i\le C,i=1,2,\dots ,m& .\end{array}$$(25)

Now the solution of equation (25) can be described as:$${\beta }_j={\displaystyle \sum _{i=1}^m}(a_i^{*}-ai)z_i(x_j).$$(26)

In equation (26), x_i is the support vector and the corresponding parameter of non-zero coefficients $a_i^{*}$, ai. Substituting equation (26) into (21), the solution can be transformed into a posterior distribution estimation of an optimization problem.

As discussed above, the PF algorithm can be modified into a new PF algorithm by integrating SVR, which can be described as follows.

Resampling of the posterior distribution starts once the effective sample N_eff below the threshold. The two training groups are particle $x_k^i$ and corresponding weight $w_k^i=F_l(x_k^i)$. The resampling posterior distribution is rebuild by these groups. The flow chart of the SVR-PF algorithm is shown in Figure 1.

In Figure 1, the rebuilt particles and weights are represented by ${\tilde{{x}}}_k^1,\dots ,{\tilde{\mathrm{x}}}_k^{\mathit{m}}$ _{and ${\tilde{\mathrm{w}}}_k^1,\dots ,{\tilde{\mathrm{w}}}_k^{\mathit{m}}$.}

Fig. 1 Fundamental illustration of SVR-PF.

3 Proposed prediction method

The proposed method mainly consists of three parts, feature construction, feature signal processing and RUL prediction, see Table 2 for details.

3.1 Feature construction

Feature signals are extracted from respective original vibration and temperature signal of crusher roller sleeves. Since the definition of the original signal in different stage is relatively vague, it is crucial to select a significant sensitive feature that can fully reflect the degradation of roller sleeve. The proposed method evaluates the degradation of roller sleeves by calculating the tendency degree between each feature and running time, which is defined as the Karl Pearson coefficient of the feature.

The Karl Pearson coefficient uses the rank to evaluate the tendency degree of a feature. It cannot only evaluate the nonlinear relationship but the monotonicity of the features$$R=\frac{{\displaystyle \sum _{i=1}^{N_s}}\left(x_i-{\displaystyle \bar{x}}\right)\left(y_i-{\displaystyle \bar{y}}\right)}{\sqrt{{{\displaystyle \sum _{i=1}^{N_s}}\left(x_i-{\displaystyle \bar{x}}\right)}^2{\left(y_i-{\displaystyle \bar{y}}\right)}^2}}$$(27)where x_i and y_i are the ranks of the time t_i, and the ith feature, respectively. N_i is the length of the time sequence.  $\overline{\mathit{x}}$ and $\overline{\mathit{y}}$ are the means of x_i and y_i, respectively.

The sensitive feature is chosen from the features which has the highest tendency value, i.e., the original feature has the most obvious monotonic trends.

3.2 Feature signal processing

In the signal processing part, original vibration and temperature signals are usually formed by the superposition of the characteristic signal and the noise signal, and the random disturbance signal in the noise affects the precision of the prediction result deeply. So, it is important to process the signal for a more accurate prediction result. Signal processing contains the removal of outliers, eliminates the trend item and denoising.

In general, the detection of outliers in signal is based on the previous normal monitoring data. The least squares polynomial is established to estimate the value of the observation data at the next moment, the absolute value of the estimated value subtracts the actual data at current time and the difference is further determined whether the difference is more than a given threshold. If the difference is more than the threshold, it is considered that the observation data are outlier, otherwise they are considered normal data.

In measurement process $\widehat{\mathit{x}}(n)=\mathrm{x}(\mathrm{n}-1)+\frac{1}{2}$ or signal processing, set x (n − 4)_{, x (n − 3), x (n − 2)}, x (n − 1) are four consecutive data of signal x(n) before time point n. The estimated value of the current time $\widehat{\mathit{x}}(n)$ can be obtained by the linear extrapolation of the least square estimation [22–25].$$\widehat{\mathit{x}}(n)=\mathit{x}(\mathit{n}-1)+\frac{1}{2}\mathit{x}(\mathit{n}-2)-\frac{1}{2}\mathrm{x}(\mathit{n}-4).$$(28)

Calculating the absolute value of $\widehat{\mathit{x}}(n)$ subtracts the measured value, and compares it with the threshold value δ, i.e.$$|\widehat{\mathit{x}}(\mathit{n})-\mathit{x}(\mathit{n})|\le \mathit{\delta }$$(29)where x(n) are current data, $\widehat{\mathit{x}}(n)$ are the estimated value of the current data obtained by the linear extrapolation through the least square estimation; σ is the standard deviation of measured data residuals. If the equation (29) is established, x(n) is the normal value, otherwise, it is the outlier.

The trend term in the measurement signal is the frequency component of the signal, which is larger than the sampling length of the signal. It is generally the result of a slow change of the time sequence in the measurement system. Except the working frequency of the original signal collected by the sensor, there are some random interference signals. The existence of these trends, will cause great error in the correlation analysis or power spectrum analysis in space domain, even distort the low frequency completely. If the measurement signal without removing the trend term is directly used to predict the RUL of roller sleeve, it will directly affect the forecast results, make inappropriate judgments and conclusions. So the extraction and elimination of the measurement signal trend term is an important part of tested data processing.

The original signal is x(t), which can get a discrete time series x(n) by uniformly-spaced sampling. The least square method is used to construct a pth-order polynomial [26–28].$$y(t)=a_0+a_{1}t+a_2{t}^2+\cdots +a_p{t}^p={\displaystyle \sum _{k=0}^p}{a}_k{t}^k$$(30)where p is a positive integer, means the order of polynomial, and the selection of p value is based on the estimation of the signal trend. If the trend of the signal is linear, choose p = 1. With y(t), we subtract the original signal x(t) by the polynomial trend term y(t), i.e.$$\widehat{\mathrm{y}}(\mathit{t})-\mathit{x}(\mathit{t})-\mathit{y}(\mathit{t})$$(31)where $\widehat{\mathrm{y}}(\mathit{t})$ is the signal that removes the trend item.

Because of the existence of so lot of noise, the field monitoring signal may submerge in other vibration signals and random noise, which can further causes great impact on the online monitoring. In this paper, wavelet is used to denoise the signals. The basic idea of wavelet denoising is to decompose and reconstruct the signal. Because signal and noise at different wavelet spectrum scales have different expressions, we remove the spectral components especially the dominant portions generated by noise at different scales. The wavelet spectrum preserved in this way is the wavelet spectrum of the original signal, basically. We reconstruct the original signal using the reconstruction algorithm of wavelet transform at last [29–32].

Signal x(t) can be expressed as:$$x(t)={\displaystyle \sum _{k=-\infty }^{+\infty }}{c}_{j,k}\varphi (t-k)+{\displaystyle \sum _{k=-\infty }^{+\infty }}{\displaystyle \sum _{j=0}^{+\infty }}{d}_{j,k}{c}_k\varphi (2^{j}t-k)$$(32)where c_j,k =⟨ x(t), ϕ_j,k(t) ⟩ is the scale coefficient, d_j,k =⟨ x(t), ϕ_j,k(t) ⟩ is the wavelet coefficient.

In the multi-scale decomposition process, x(t) is always progressively decomposed to two subspaces V_j and W_j from the space V_{j−1 .} According to the two-scale equation, we can get the fast recursive algorithm about projection coefficient from c_j−1,k of x(t) in V_j−1 to c_j,k and d_j,k of x(t) in V_j and W_{j .} $$c_{j,k}={\displaystyle \sum _{m\in z}}h\left(m-2k\right)c_{j-1,m}$$(33) $$d_{j,k}={\displaystyle \sum _{m\in z}}g\left(m-2k\right)c_{j-1,m.}$$(34)

On the contrary, c_j−1,k also can be reconstructed by c_j,k and d_j,k, and the reconstruction formula is as follows.$$c_{j-1,k}={\displaystyle \sum _m}{c}_{j,m}h\left(k-2m\right)+{\displaystyle \sum _m}{d}_{j,m}g\left(k-2m\right).$$(35)

3.3 RUL prediction

3.3.1 Initial state of fusion prediction

Set N as the initial time of prediction, a_N is the acceleration in the Nth time point, $a_{T,N}^{*}$ is the acceleration prediction obtained by temperature data in the Nth time point,$a_{a,N}^{*}$ is the acceleration prediction obtained by acceleration data in the Nth time point.

The first step is calculating the initial value of m_1,i(T) BPA function. From the central limit theorem, a large number of temperature data measurement errors obeys normal distribution, i.e.$$T_{T,N}\sim N\left({\mu }_{T_{T,N}},{\sigma }_T^2\right)$$(36)where ${\mu }_{T_{T,N}}$ is the mean value of temperature, ${\sigma }_T^2$ is the variance.

From the experimental result we can see that there is a strong linear relationship between acceleration signal and temperature signal. So, the a_T,N estimation also obeys normal distribution, i.e.$$a_{T,N}\sim N({\mu }_{a_{T,N}},{\sigma }_T^2)$$(37)where$${\mu }_{a_t,N}={\alpha }_N*{\mu }_{T,N}+{\beta }_N=a_N$$(38) $${\sigma }_T^2={\alpha }_N^2*{\sigma }_T^2.$$(39)

Thus the initial value of the BPA function m_1,i(T) can be described as m_1,N+1(T)$$m_{1,N+1}(T)=\frac{1}{\sqrt{2\pi }{\sigma }_T}\exp \left[-\frac{{\left(a_{T,N}^{*}-a_N\right)}^2}{2{\sigma }_T^2}\right].$$(40)

The second step is calculating the initial value of the BPA function m_2,i(a). From the central limit theorem, a large number of acceleration data measurement errors obeys the normal distribution, i.e.$$a_{a,N}\sim N({\mu }_{a_a,N},{\sigma }_a^2)$$(41)where μ_a,N is the mean value of acceleration; ${\sigma }_a^2$ is the variance. If μ_a,N meets:$${\mu }_{a_a,N}=a_N.$$(42)

The initial value of the BPA function m_2,i(a) can be described as m_2,N+1(a)$$m_{2,N+1}(a)=\frac{1}{\sqrt{2\pi }{\mathit{\sigma }}_{\mathit{a}}}\exp [-\frac{({\mathit{a}}_{\mathit{a},\mathit{N}}^{*}-{\mathit{a}}_{\mathit{N}}{)}^2}{{2\mathrm{\sigma }}_{\mathrm{a}}2}]\mathrm{.}$$(43)

After getting the value of BPA function m_1,N+1(a) and m_2,N+1(a), then we calculate the value of the other BPA functions. Because there is no correlation between the two kinds of prediction methods, it's easy to know that:$$m_{1,N+1}(a)=m_{2,N+1}(T)=0.$$(44)

According to the properties of BPA function:$$m_{1,N+1}(T\hspace{0.17em}\cup \hspace{0.17em}a)=1-m_{1,N+1}(T)-\hspace{0.17em}{m}_{1,N+1}(a)=1-m_{1,N+1}(T)$$(45) $$m_{2,N+1}(T\hspace{0.17em}\cup \hspace{0.17em}a)=1-m_{2,N+1}(T)-m_{2,N+1}(a)=1-m_{2,N+1}(a).$$(46)

After determining all the values of the BPA function, the third step is calculating the posterior fusion BPA function. The fusion BPA function can be described as m, from formulas (6) and (7), it's easy to know that:$$\begin{array}{ll}{m}_{N+1}(T)\hfill & =\frac{1}{1-K_{N+1}}{\displaystyle {\displaystyle \sum }_{B\cap C=T}}{m}_{1,N+1}(B)m_{2,N+1}(C)\hfill \\ \hfill & =\frac{m_{1,N+1}(T)m_{2,N+1}(T\hspace{0.17em}\cup \hspace{0.17em}a)}{1-K_{N+1}}\hfill \end{array}$$(47) $$\begin{array}{ll}{m}_{N+1}(a)\hfill & =\frac{1}{1-K_{N+1}}{\displaystyle {\displaystyle \sum }_{B\cap C=a}}{m}_{1,N+1}(B)m_{2,N+1}(C)\hfill \\ \hfill & =\frac{m_{1,N+1}(T\hspace{0.17em}\cup \hspace{0.17em}a)m_{2,N+1}(a)}{1-K_{N+1}}\hfill \end{array}$$(48)where$$K_{N+1}={\displaystyle {\displaystyle \sum }_{B\cap C=\varnothing }}{m}_{1,N+1}(B)m_{2,N+1}(C)=m_{1,N+1}(a)m_{2,N+1}(T)+m_{1,N+1}(T)m_{2,N+1}\left(a\right)=m_{1,N+1}(T)m_{2,N+1}(a).$$(49)

With the fusion posterior BPA function, the acceleration can be predicted in the N+1th time point. Set $a_{T,N+1}^{*}$ is the acceleration prediction obtained by temperature data in the N+1th time point,  $a_{a,N+1}^{*}$ is the acceleration prediction obtained by acceleration data in the N+1 ^th time point, $a_{N+1}^{*}$ is the acceleration prediction obtained by data fusion method in the N+1th time point. Based on formulas (47)–(49), $a_{N+1}^{*}$ meets:$$a_{N+1}^{*}=m_{N+1}(T)a_{T,N+1}^{*}+m_{N+1}(a)a_{a,N+1}^{*}.$$(50)

Thus we get the initial state of RUL prediction based on D-S data fusion and SVR-PF:$$\text{Initial}={[m_{1,N+1}(T),m_{1,N+1}(a),m_{1,N+1}(T\hspace{0.17em}\cup \hspace{0.17em}a),m_{2,N+1}(T),m_{2,N+1}(a)\cdots m_{2,N+1}(T\hspace{0.17em}\cup \hspace{0.17em}a),m_{N+1}(T),m_{N+1}(a),a_{N+1}^{*}]}^T.$$(51)

3.3.2 Fusion prediction process

Step 1: calculating the BPA functions. Set N_EOL is the end of roller sleeve working life, $a_k^{*}$ _{is the} acceleration prediction obtained by data fusion in kth time point, $a_{T,k}^{*}$ is the acceleration prediction obtained by the temperature data, $a_{a,k}^{*}$ is the acceleration prediction obtained by the acceleration data. Where N + 1 ≤ k ≤ N_EOL, the BPA function can be described as m_1,k+1(T) and m_2,k+1(a).$$m_{1,k+1}(T)=\frac{1}{\sqrt{2\pi }{\sigma }_T}\exp \left[-\frac{({\mathit{a}}_{\mathit{T},\mathit{k}}^{*}-{\mathit{a}}_{\mathit{k}}^{*}{)}^2}{\mathit{2}\mathit{\sigma }{\mathit{T}}^2}\right]$$(52) $$m_{2,k+1}(a)=\frac{1}{\sqrt{2\pi }{\sigma }_a}\exp [-\frac{({\mathit{a}}_{\mathit{a},\mathit{k}}^{*}-{\mathit{a}}_{\mathit{k}}^{*})}{{2\mathrm{\sigma }}_{\mathit{a}}^{*}}].$$(53)

Because there is no correlation between the two prediction models, it's easy to know that:$$m_{1,k+1}(a)=m_{2,k+1}(T)=0.$$(54)

Based on the properties of BPA functions, there are:$$m_{1,k+1}(T\cup a)=1-m_{1,k+1}(T)$$(55) $$m_{2,k+1}(T\cup a)=1-m_{2,k+1}(a).$$(56)

Step 2: calculating the acceleration prediction in the k+1th time point. From the formulas (6) and (7), formulas (52)–(56), the posterior fusion BPA function m_k+1(T) and m_k+1(a) can be described as:$$\begin{array}{ll}{m}_{k+1}(T)\hfill & =\frac{1}{1-K_{k+1}}{\displaystyle {\displaystyle \sum }_{B\cap C=T}}{m}_{1,k+1}(B)m_{2,k+1}(C)\hfill \\ \hfill & =\frac{m_{1,k+1}(T)m_{2,k+1}(T\cup a)}{1-k_{k+1}}\hfill \end{array}$$(57) $$\begin{array}{l}{m}_{k+1}(a)=\frac{1}{1-K_{k+1}}{\displaystyle {\displaystyle \sum }_{B\cap C=a}}{m}_{1,k+1}(B)m_{2,k+1}(C)\\ =\frac{m_{1,k+1}(T\cup a)m_{2,k+1}(a)}{1-k_{k+1}}.\end{array}$$(58)

Thus we obtain the acceleration prediction $a_{k+1}^{*}$ in the K+1th time point through D-S data fusion.$$a_{k+1}^{*}=m_{k+1}(a)a_{T,k+1}^{*}+m_{k+1}(a)a_{a,k+1}^{*}.$$(59)

Step 3: determining whether the acceleration reaches the threshold. If not, then go back to step 1 and continue the prediction; otherwise, calculating the RUL prediction $\overline{\mathit{L}}*$:$$\overline{\mathit{L}}*={\mathit{N}}_{\mathit{E}\mathit{O}\mathit{L}}-\mathit{N}=(\mathit{k}+1)-\mathit{N}.$$(60)

In summary, the proposed RUL prediction method based on D-S data fusion and SVR-PF can be described as flow chart in Figure 2.

A prediction model based on D-S data fusion and SVR-PF is established as:$$\{\begin{array}{c}\hfill X_{k+1}=f\left(X_k,V_k\right)\hfill \\ \hfill Y_k=g\left(X_k,N_k\right)\hfill \end{array}$$(61)where X_k is the prediction state, X_k, V_k are noise.$$X_k={[X_{,k}^T,X_{a,k}^T,X_{DS,k}^T]}^T$$(62)where $X_{T,k}^T$ is the state obtained by the analysis of temperature data; $X_{a,k}^T$ is the state obtained by the analysis of acceleration data; $X_{DS,k}^T$ is the state obtained by the analysis of D-S data fusion, and there are:$$X_{T,k}={[{\lambda }_{T,k}^{*},T_{T,k}^{*},a_{T,k}^{*}]}^T$$(63) $$X_{a,k}={[{\lambda }_{T_{T,k}}^{*},T_{T,k}^{*},a_{a,k}^{*}]}^T$$(64) $$X_{DS,k}={[m_{1,k}(T),m_{1,k}(T\cup a),m_{2,k}(a),m_{2,k}(T\cup a)\dots m_k(T),m_k(a),a_k^{*}]}^T.$$(65)

In formulas (63)–(65): λ_T,k is the degradation parameter of temperature in the kth time point [33]; T_T,k, ${\lambda }_{T_{T,k}}$ are the temperature degradation parameters in the kth time point [34]; * is the prediction of each corresponding variable.

Thus, in combination with the RUL prediction model of literature [33,34] we can obtain the state equation of prediction model (61), the partial state equation of prediction by temperature data can be expressed as follows.$$\{\begin{array}{l}{\lambda }_{T_{T,k+1}}^{*}={\lambda }_{T_{T,k}}^{*}+v_{a,k}\\ T_{T,k+1}^{*}=T_{T,k}^{*}\exp \left({\lambda }_{T_{T,k^{\Delta k}}}^{*}\right)+v_{b,k}\\ a_{T,k+1}^{*}={\alpha }_N*T_{T,k}+{\beta }_N+v_{c,k}\end{array}$$(66)where α_N, β_N is the degradation parameter prediction of acceleration in the Nth time point through the analysis of the temperature data. The partial state equation of prediction by acceleration data can be described as follows.$$\{\begin{array}{l}{\lambda }_{T_{T,k}}^{*}={\lambda }_{T,k}^{*}+v_{1,k}^{*}\\ T_{T,k+1}^{*}=T_{T,k}^{*}\exp \left({\lambda }_{T_T,k}^{*}\right)+v_{2,k}^{*}\\ a_{a,k+1}^{*}=T_{T,k+1}^{*}\exp \left({\lambda }_{T_{T,k+1}}^{*}\right)+v_{3,k}^{*}\end{array}$$(67)

The partial state equation of prediction by data fusion can be described as follows.$$\{\begin{array}{l}{m}_{1,k+1}\left(T\right)=\frac{1}{\sqrt{2\pi }{\sigma }_T}\exp \left[-\frac{({\mathit{a}}_{\mathit{T},\mathit{k}}^{*}-{\mathit{a}}_{\mathit{k}}^{*}\mathrm{^2}}{{2\mathrm{\sigma }}_{\mathit{T}}^2}\right]\\ m_{2,k+1}(a)=\frac{1}{\sqrt{2\pi }{\sigma }_a}\exp \left[-\frac{({\mathit{a}}_{\mathit{a},\mathit{k}}^{*}-{\mathit{a}}_{\mathit{k}}^{*}\mathrm{^2}}{{2\mathrm{\sigma }}_{\mathit{a}}^2}\right]\\ m_{1,k+1}(T\hspace{0.17em}\cup \hspace{0.17em}a)=1-m_{1,k+1}(T)\\ m_{2,k+1}(T\hspace{0.17em}\cup \hspace{0.17em}a)=1-m_{2,k+1}(a)\\ K_{k+1}=m_{1,k+1}(T)m_{2,k+1}(a)\\ m_{k+1}(T)=\frac{m_{1,k+1}(T)m_{2,k+1}(T\hspace{0.17em}\cup \hspace{0.17em}a)}{1-K_{k+1}}\\ m_{k+1}(a)=\frac{m_{1,k+1}(T\hspace{0.17em}\cup \hspace{0.17em}a)m_{2,k+1}(a)}{1-K_{k+1}}\\ {\mathit{a}}_{\mathit{k}+1}^{*}=m_{k+1}(T){\mathit{a}}_{\mathit{T},\mathit{k}+1}^{*}+m_{k+1}(a){\mathit{a}}_{\mathit{a},\mathit{k}+1}^{*}\end{array}$$(68)

Combining formulas (66)–(68) can obtain the state equation of prediction model (61), the measurement equation (61) can be described as follows.$$\{\begin{array}{l}{\widehat{\mathrm{T}}}_{\mathit{T},\mathit{k}}^{*}={\mathit{T}}_{T,k}^{*}+{\mathit{n}}_{\mathit{a},\mathit{k}}\\ {\mathit{T}}_{\mathit{a},\mathit{k}}=T_{T,k}^{*}+{\mathit{n}}_{1,k}^{*}\\ a_k^{*}={\mathit{m}}_{\mathit{k}}(\mathrm{T}){\mathit{a}}_{\mathit{T},\mathit{k}}^{*}+{\mathit{m}}_{\mathit{k}}(\mathit{a}){\mathit{a}}_{\mathit{a},k}^{*}\end{array}$$(69)where  ${\widehat{T}}_{T,k}$ is the temperature prediction obtained by temperature data, and T_a,k is the degradation parameter prediction of acceleration obtained by temperature data.

Fig. 2 Flowchart of the proposed method.

4 Experimental demonstration

4.1 Introduction to the data acquisition platform

The test platform layout is shown in Figure 3, named PRONOSTIA [35]. The testing platform is designed by the AS2M Department of the FEMTO-ST Association, the full life test of the roller sleeve is carried out on the data acquisition platform of the rolling bearing, the vibration signal is collected by the 3035B Dytran acceleration sensor (the maximum acquisition range is 50 g), temperature signal acquisition using JCJ100TLB temperature sensor (maximum acquisition range is 200 °C). Because the acceleration signal is more severe than the temperature signal, so the full-life test process stops if the acceleration signal amplitude is found to exceed 20 g. Even if the roller sleeve does not out of work, in order to avoid the test platform damage caused by the roller sleeve, we determine the failure of the roller sleeve, stop testing. The acceleration test sampling frequency is 25.6 kHz, each 10 s stores a set of data, each group of data 2560 points, the temperature test sampling frequency is 10 Hz, and each 10 s stores a set of data, each group of data 100 points.

In the test, the testing roller sleeve is 22,324 tapered roller bearing, the roller life-test is carried out 4 times, each time one test roll is damaged. The 1st test and the 2nd test are worked under the condition of radial load 4000 kN, speed 1800 rpm/min; the 3rd test and the 4th test are worked under the condition of radial load 4200 kN, speed 1650 rpm/min. The test results are shown in Table 3.

Due to the different working state of the 4 roller sleeves and the different structure of the roller sleeves, the experimental results are different, which is conform to the actual engineering facts. Then we use the measured experimental data to verify the performance of the proposed RUL prediction method. Because the 4 groups tests are carried out on the same platform, and the analysis methods of each roller sleeve is same, below takes the 1st roller sleeve as the research object to make explanation. The measured temperature and vibration data are shown in Figures 4 and 5.

Fig. 3 Overview of the data acquisition platform.

Fig. 4 Original temperature data.

Fig. 5 Original acceleration data.

4.2 Feature construction

In order to compare the prediction performance of proposed D-S data fusion and SVR-PF prediction method with the prediction method using the single acceleration data and the prediction method with the temperature data with finite data available, the first key step is to select a good feature signal.

In this paper, the features are selected by calculating the Karl Pearson correlation coefficient between the time-domain characteristics of temperature and acceleration and RUL. As an indicator the features whose correlation coefficient is highest are selected. As a result, the root mean square (RMS) feature of vibration and the absolute mean value feature of temperature are selected. The Karl Pearson correlation coefficient result is shown in Table 4. The feature signal of acceleration and temperature are shown in Figures 6 and 7.

Fig. 6 Temperature feature signal.

Fig. 7 Acceleration feature signal.

4.3 Feature signal processing

4.3.1 Removal of outliers

The first step of signal processing is the removal of outliers. For the nonlinear and non-stationary signal, the existence of outliers can produce spurious harmonic components, further can influence the prediction accuracy. According to the statistical properties of the original data, 3σ criterion is used to remove the outliers here. If the residuals in equation (29) exceed 3σ the outliers can be eliminated. The temperature and vibration signals after removing the outliers are shown in Figures 8 and 9.

Fig. 8 Temperature signal after removing outliers.

Fig. 9 Acceleration signal after removing outliers.

4.3.2 Remove the trend of a smooth

Because of the zero drift of the amplifier caused by temperature variation, the performance of low frequency which exceeds the frequency range of the sensor is not stable with ambient interference around the sensor etc., which caused the collected data of vibration signal and temperature signal in life-test will often deviate from the baseline, and even the degree of the deviation from the baseline will vary over time. The whole process of the deviation from the baseline directly affects the correctness of the signal and should be removed as the trend term. This paper from the perspective of engineering application adopts a simple and practical method to remove the trend items − the modified function method. The temperature and vibration signals after removing the trend term are shown in Figures 10 and 11.

Fig. 10 Smoothed temperature signal.

Fig. 11 Smoothed acceleration signal.

4.3.3 Denoising

Wavelet analysis is known as the microscope of signal processing. The key of wavelets analysis is the selection of wavelet basis and the decomposition level. Decomposition layer has great influence on the effect of denoising. The more is the decomposition layer, the lower is the noise–signal ratio. Meanwhile, when the layer increases, the processing becomes slow. Although few decomposition layer has high noise-signal ratio, the signal is decomposed to very small frequency bandwidths. Only the high frequency coefficients can be processed to remove the corresponding noise, while the corresponding low frequency noise is all reserved. Therefore, the choice of wavelet decomposition layer should be neither too large as considering the improvement of the noise-signal ratio nor too small as considering the suppression of low frequency noise. The purpose of denoising is to get the useful feature signal, so wavelet coefficients can reflect the minimum frequency components in the useful signal. The wavelet decomposition is to decompose the signal into various independent bands, high detail coefficient reflects the low-frequency part of the signal. So this paper is based on the minimum frequency signal to determine the maximum level of wavelet decomposition. In this paper, sym8 wavelet is chosen as the wavelet base and soft thresholds are used to denoise. The temperature and vibration signal denoising processes are shown in Figures 12–17.

Fig. 12 Low frequency temperature signal.

Fig. 13 High frequency temperature signal.

Fig. 14 Low frequency acceleration signal.

Fig. 15 High frequency acceleration signal.

Fig. 16 Denoising temperature signal.

Fig. 17 Denoising acceleration signal.

4.4 RUL prediction

We compare the prediction performance of proposed prediction method which based on D-S data fusion and SVR-PF with prediction method which uses single data source and other prediction methods. Roller sleeve 1 was used to predict in three cases, predicted by acceleration data, predicted by temperature data and predicted by fused data.

From Figures 18–20, it's easy to know that the results predicted by the proposed prediction method are more accurate than other prediction methods.

Fig. 18 Predicted acceleration by temperature.

Fig. 19 Predicted acceleration by acceleration.

Fig. 20 Fusion prediction.

5 Conclusion

In view of the engineering problem that is difficult to accurately predict the remaining useful life under partially observed state, a new method based on D-S data fusion and SVR-PF is proposed.

From Table 5, we can see that compared to other prediction methods such as obtained by the temperature or acceleration-based data-driven, the prediction accuracy of proposed method is significantly improved. Meanwhile, it provides a basis to make maintenance decision for the equipment working under severe conditions, further reduces the maintenance cost, improves the utilization rate and the reliability of the equipment, which has good practicability and popularization value.

Nomenclature

SOH: State of health

PHM: Prognostics and health management

BPA: Basic probability assignment

m₁(·) : BPA function 1 under the identification framework

m₂(·): BPA function 2 under the identification framework

m_1,2(·): BPA function after fusion m₁ and m₂ under the identification framework

K : Degree of conflict between the two evidences

T : Temperature magnitude

a: Acceleration magnitude

a_N : Acceleration magnitude at time N

$a_{T,N}^{*}$ : Prediction acceleration magnitude at time N obtained by temperature data

a_T,N : Acceleration magnitude at time N obtained by temperature data

$a_{a,N}^{*}$ : Prediction acceleration magnitude at time N obtained by acceleration data

R : Value of Karl Pearson coefficient

^α_N,β_N : Acceleration degradation parameters

m_1,i(T): BPA function of temperature at time i

$a_{T,N+1}^{*}$ : Prediction acceleration magnitude at time N+1 obtained by temperature data

K_N+1 : Degree of conflict between the two evidences at time N+1

N_EOL : Threshold of life cycle N

$\widehat{\mathrm{L}}*$ : Prediction of remaining useful life

X_k : State of prediction at time k

V_k,N_k : Measurement noise at time k

f(⋅),g(⋅): Transition function and measurement function

$X_{T,k}^T$ : State of temperature at time k

$X_{a,k}^T$ : State of acceleration at time k

$X_{DS,k}^T$ : State of fusion at time k

λ_T,k : Temperature degradation parameters obtained by temperature at time k

${\mathrm{\lambda }}_{\mathit{T},\mathit{k}}^{*}$ : Prediction temperature magnitude obtained by temperature at time k

$T_{T,k}^{*}$ : Prediction temperature magnitude obtained by temperature at time k

$a_{T,k}^{*}$ : Prediction acceleration magnitude obtained by temperature at time k

X_a,k : State of acceleration at time k

${\lambda }_{T_1,k}^{*}$,$T_{1,k}^{*}$ : Acceleration degradation parameters obtained by acceleration at time k

${\lambda }_{T,k+1}^{*}$ : Prediction temperature degradation parameters obtained by acceleration at time k+1

v_1,k : λ degradation parameters obtained by temperature at time k

${\lambda }_{T,k^{\Delta k}}^{*}$ : Distribution state of λ at time k

v_2,k : State noise of prediction temperature obtained by temperature at time k

v_3,k : State noise of prediction acceleration obtained by temperature at time k

$v_{1,k}^{*}$ : λ degradation parameters obtained by acceleration at time k

$v_{2,k}^{*}$ : State noise of prediction acceleration obtained by temperature at time k

$v_{3,k}^{*}$ : Measurement noise of prediction acceleration obtained by acceleration at time k

n_1,k : Measurement noise of prediction temperature at time k

$n_{1,k}^{*}$ : Measurement noise of prediction acceleration degradation parameters at time k

References

All Tables

All Figures

	Fig. 1 Fundamental illustration of SVR-PF.
In the text

	Fig. 2 Flowchart of the proposed method.
In the text

	Fig. 3 Overview of the data acquisition platform.
In the text

	Fig. 4 Original temperature data.
In the text

	Fig. 5 Original acceleration data.
In the text

	Fig. 6 Temperature feature signal.
In the text

	Fig. 7 Acceleration feature signal.
In the text

	Fig. 8 Temperature signal after removing outliers.
In the text

	Fig. 9 Acceleration signal after removing outliers.
In the text

	Fig. 10 Smoothed temperature signal.
In the text

	Fig. 11 Smoothed acceleration signal.
In the text

	Fig. 12 Low frequency temperature signal.
In the text

	Fig. 13 High frequency temperature signal.
In the text

	Fig. 14 Low frequency acceleration signal.
In the text

	Fig. 15 High frequency acceleration signal.
In the text

	Fig. 16 Denoising temperature signal.
In the text

	Fig. 17 Denoising acceleration signal.
In the text

	Fig. 18 Predicted acceleration by temperature.
In the text

	Fig. 19 Predicted acceleration by acceleration.
In the text

	Fig. 20 Fusion prediction.
In the text