Tensor Singular Spectrum Decomposition Algorithm Based on Permutation Entropy for Rolling Bearing Fault Diagnosis

Abstract

Mechanical vibration signal mapped into a high-dimensional space tends to exhibit a special distribution and movement characteristics, which can further reveal the dynamic behavior of the original time series. As the most natural representation of high-dimensional data, tensor can preserve the intrinsic structure of the data to the maximum extent. Thus, the tensor decomposition algorithm has broad application prospects in signal processing. High-dimensional tensor can be obtained from a one-dimensional vibration signal by using phase space reconstruction, which is called the tensorization of data. As a new signal decomposition method, tensor-based singular spectrum algorithm (TSSA) fully combines the advantages of phase space reconstruction and tensor decomposition. However, TSSA has some problems, mainly in estimating the rank of tensor and selecting the optimal reconstruction tensor. In this paper, the improved TSSA algorithm based on convex-optimization and permutation entropy (PE) is proposed. Firstly, aiming to accurately estimate the rank of tensor decomposition, this paper presents a convex optimization algorithm using non-convex penalty functions based on singular value decomposition (SVD). Then, PE is employed to evaluate the desired tensor and improve the denoising performance. In order to verify the effectiveness of proposed algorithm, both numerical simulation and experimental bearing failure data are analyzed.

1. Introduction

Rolling bearing is one of the most widely used and easily damaged rotating machine elements, the operation status of which is directly related to the operation of machinery and production efficiency. Thus, it has a very high practical value to achieve accurate diagnosis and recognition of rolling bearing fault. To maintain smooth operation, the diagnostic object should be changed from the traditional super-fault to early weak fault. The diagnostic methods should also be changed from linear analysis methods to non-traditional analysis methods. Many scholars have proposed different methods to accurately identify and extract the fault characteristics in frequency domain. For instance, Short-time Fourier transform (STFT) [1] is proposed to make up the shortcomings of Fourier transform so that it can be used to analyze non-stationary signal, but the time window width of STFT is fixed and cannot be adapted to the changes in the signal. Other non-stationary signal processing techniques also have a variety of shortcomings. For instance, Empirical Mode Decomposition (EMD) is utilized in the analysis of non-stationary and nonlinear signal [2,3], which easily causes the phenomenon of modal aliasing [4] and edge effect [5]. The selection of decomposition level and wavelet basis function has a great impact on analytical results for Wavelet Transform (WT) [6,7,8], which lacks in adaption to signal processing. Hilbert–Huang transform (HHT) operator [9] owns windowing effects and the demodulated results will be non-transient, so maybe not accurate. Wigner–Ville distribution [10] for multi-component signal will produce crosstalk and spurious frequency components. In addition, there are some time–frequency analysis methods applied to speed signals in rotating systems [11,12,13,14,15,16,17]. In conclusion, these methods have some limitations in non-stationary, non-periodic, non-linear signal processing.

High-dimensional signal contains a lot of information, which is hard to be identified and may not be displayed in one-dimensional space. Therefore, aiming to decompose original signal to extract useful components, high-dimensional signal processing technology may be a breakthrough with great practical engineering significance. Traditional data representational models, such as vectors and matrices, are one-dimensional time series. It is clear that there is some structure in the observed signal allowing one to model it much more compactly. Commonly, the observed signals can be compressible and be expressed as parsimonious models. A higher-order model such as tensor is served as one way to achieve it, which is the most natural representation of high-dimensional data and maintains the intrinsic structural properties of the data. Theoretically, tensors can be seen as higher-order generalizations of vectors and matrices. Tensorization of the observed signal using segmentation and decomposition of the generated tensor are two main aspects of our approach. Tensor decomposition can dig out the potential high-dimensional data structure in the form of low-dimensional subspace, thus be able to reflect the essential characteristics. Generally, the concept of tensor was first proposed by Hitchcock [18]. In 1970, to the psychometrics community, the tensor is developed to the form of CANDECOMP (canonical decomposition) by Carroll [19] and PARAFAC (parallel factors) by Harshman [20]. Then, the research of tensor has gradually been extended to other areas such as signal processing, computer vision, numerical analysis and data mining. CANDECOMP/PARAFAC tensor decomposition, namely CP decomposition model, is a commonly used method. The CP decomposition factorizes a tensor into a sum of R-component rank-one tensors [21]. By the proposed model, three factor matrixes representing the combination of the vectors can be obtained from the rank-one components. Recently, tensor-based singular spectrum algorithm (TSSA) was proposed by Saeid [22] and has been applied to the field of EEG signal processing, which provides an effective way for solving the above problem. The main idea is the one-dimensional times series can be segmented as a matrix using a non-overlapping window. Then each row of matrix can be expressed as a reconstructed attractor matrix through phase space reconstruction [23]. The reconstructed attractor matrix forms the corresponding slice of the tensor, thus we have a 3D tensor to be decomposed. After that, the above-mentioned CP tensor decomposition method is used here. The key step is performed by the alternating least squares method (ALS) [24] to obtain the three factor matrices. TSSA method combines the advantages of phase space reconstruction, singular spectrum analysis [25] as well as tensor decomposition. However, the rank of tensor is difficult to estimate and the optimal selection among the decomposed tensor for reconstruction has a great influence on the final result.

In this paper, an improved TSSA decomposition method is proposed by applying the convex optimization for the rank estimation and permutation entropy (PE) for desired tensor selection. The rank of tensor in TSSA is usually identified by the number of non-zero singular value, while the singular value decomposition (SVD) algorithm is easily affected by the noise interference [26,27]. Consequently, a novel SVD algorithm using convex optimization framework is introduced in this paper. A parameterized non-convex penalty function is put forward to the singular value decomposition and the number of non-zero singular values is used to determine the rank of tensor [28,29,30,31]. Bandt [32] proposed the concept of average entropy parameter in the application of measuring the complexity of one-dimensional time series, namely PE. It is an algorithm study about describing irregular and nonlinear systems [33], which cannot or hard to be quantitatively described, in a relatively simple way. The PE has the advantages of simple calculation, strong anti-noise ability, and high sensitivity to signal change. It can also be applied to the detection of the weak signal [34], and the dynamic mutation of complex system. When an abnormal failure or different kinds of faults occur in gear and bearing parts during operation process, the influence of the nonlinear factors and signal complexity are different. Hence, the vibration signals obtained from the mechanical system will also change, resulting in different values of PE. Thus, PE is employed as novel approach to select the desired tensor for reconstruction in this paper. In order to verify the validity of the proposed method, bearing fault vibration signal and fault simulation test bed are used for analysis. The proposed method is also compared with the traditional singular spectrum analysis (SSA) and EMD.

The structure of this paper is arranged as follows: In Section 2, tensor singular spectrum decomposition algorithm based on PE is briefly described. Then, the performance of abnormal signal detection using PE is focused on. The analysis results of numerical simulation signal and bearing fault signal are, respectively, described in the Section 3 and Section 4. Section 5 introduces the conclusions.

2. Theory Description

2.1. Tensor Singular Spectrum Analysis

The proposed TSSA method in this paper mainly contains two stages. The first stage includes an embedding operation and the second stage is to decompose a 3D tensor. In the stage of embedding progress, a one-dimensional time series $x(i),i=1,2,\dots ,n$ with length n is mapped into a 3D tensor $\underset{\_}{\mathit{X}}$.

There are two steps in the embedding stage. Firstly, $x(i)$ is segmented by a non-overlapping window with the size $L$. A $[n/L]\times L$ matrix $\mathit{T}$ is obtained from $x(i)$ (where $I=n/L$):

$$\mathit{T}=\left(\begin{array}{llll}x(1)& x(2)& \dots & x(L)\\ x(L+1)& x(L+2)& \dots & x(2L)\\ \dots \\ x(I-1)*L& x(I-1)*L+1& \dots & xn\end{array}\right)$$

(1)

Then, the tensor $\underset{\_}{\mathit{X}}$ is obtained by matrix $\mathit{T}$ as demonstrated in Figure 1. Each slice of the tensor can be regarded as a reconstructed attractor matrix by a windowed version of matrix $\mathit{T}$. The segmentation is executed in one direction during this procedure.

The slice ${\underset{\_}{\mathit{X}}}_{i::}$ of tensor $\underset{\_}{\mathit{X}}$ in Figure 1 is formed from the i-th row of matrix $\mathit{T}$ using phase space reconstruction, where $J$ is the length of reconstructed window, $K$ is reconstructed embedding dimension, $\tau$ is the delay time. Consequently, the relationship of $L=(K-1)\times \tau +J$ can be obtained. Our way of converting a matrix to tensor can be described as follows:

$${\underset{\_}{\mathit{X}}}_{ijk}=\mathit{T}(i,(j+(k-1)\times \tau ))i=1,2,\cdots ,I;j=1,2,\cdots ,J;k=1,2,\cdots ,K$$

(2)

Thus, we can obtain a 3D tensor of size $I\times J\times K$. In the second stage, the CP tensor decomposition method for a 3D tensor is performed, which factorizes a tensor into a sum of component rank-one tensors. Naturally, the generalized bilinear principal component analysis is similar with the tensor decomposition [35]. Generally, the detailed CP model can be expressed as [36,37]:

$$\underset{\_}{\mathit{X}}={\displaystyle \sum _{r=1}^R{\mathit{a}}_r\circ {\mathit{b}}_r\circ {\mathit{c}}_r}+\underset{\_}{\mathit{e}}$$

(3)

where $R$ is corresponding to the rank of tensor $\underset{\_}{\mathit{X}}$. ${\mathit{a}}_{\mathit{r}}\in R^{I\times 1}$, ${\mathit{b}}_{\mathit{r}}\in R^{J\times 1}$ and ${\mathit{c}}_{\mathit{r}}\in R^{K\times 1}$ are the vector elements of factor matrix $\mathit{A}\in R^{I\times R}$, $\mathit{B}\in R^{J\times R}$ and $\mathit{C}\in R^{K\times R}$. The tensor of $\underset{\_}{\mathit{e}}\in R^{I\times J\times K}$ is the residual term. Hence, the CP model can be approximately explained by factor matrix $〚\mathit{A},\mathit{B},\mathit{C}〛$:

$$\underset{\_}{\mathit{X}}\approx 〚\mathit{A},\mathit{B},\mathit{C}〛={\displaystyle \sum _{r=1}^R{\mathit{a}}_r\circ {\mathit{b}}_r\circ {\mathit{c}}_r}+\underset{\_}{\mathit{e}}$$

(4)

The above mentioned decomposition model is shown in Figure 2.

Since TSSA algorithm uses the iterative least squares method to seek factor matrix; we can define an error function as follows:

$$f(\mathit{A},\mathit{B},\mathit{C})=\frac{1}{2}{\Vert (\underset{\_}{\mathit{X}}-[[\mathit{A},\mathit{B},\mathit{C}]])\Vert }^2=\frac{1}{2}{\displaystyle \sum _{i=1}^I{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=1}^K{({\underset{\_}{\mathit{X}}}_{ijk}-{\displaystyle \sum _{r=1}^R{a}_{ir}{b}_{jr}{c}_{kr}})}^2}}}$$

(5)

First, the random factor matrix is set as an initial matrix to $\mathit{A}$, $\mathit{B}$, $\mathit{C}$. Then, $\mathit{B}$ and $\mathit{C}$ are fixed to solve for $\mathit{A}$. Similarly, $\mathit{A}$ and $\mathit{C}$ are fixed to solve for $\mathit{B}$, and next $\mathit{A}$ and $\mathit{B}$ are fixed to solve for $\mathit{C}$ in an alternating algorithm until reaching the convergence condition.

2.2. The Rank Estimation of Tensor Based on Convex Optimization

Similar to the matrix rank, the problem of rank determination has an important significance for tensor decomposition. This paper firstly proposed an algorithm based on SVD using convex optimization to estimate the rank of tensor. SVD is a classic subspace decomposition algorithm that constructs trajectory matrix based on the observed time series. The singular value is obtained by the trajectory matrix, which represents the different dynamic characteristics about each component in the original time series. For the one-dimensional time series $x(i)$, based on the theory of phase space reconstruction by selecting the appropriate embedding dimension $m$ and the delay time $\tau$, the neighborhood matrix or trajectory matrix can be expressed as:

$$\mathit{Y}=\left(\begin{array}{cccc}x(1)& x(1+\tau )& \cdots & x(1+(m-1)\tau )\\ x(1)& x(2+\tau )& \cdots & x(2+(m-1)\tau )\\ x(j)& x(j+\tau )& \cdots & x(j+(m-1)\tau )\\ ⋮& ⋮& ⋮\\ x(K)& x(K+\tau )& \cdots & x(K+(m-1)\tau )\end{array}\right)$$

(6)

where $j=1,2,\dots ,K$, $K=n-(m-1)\tau$. Generally, the trajectory matrix $\mathit{Y}$ consists of the useful signal and noise component. Then, SVD is applied to $\mathit{Y}$ by the form of $\mathit{Y}=\mathit{U}\mathit{\Sigma }{\mathit{V}}^{\mathrm{T}}$, where $\Sigma =\left[\begin{array}{ll}S& 0\\ 0& 0\end{array}\right]$ is a diagonal matrix, $\mathit{S}=\mathrm{diag}({\sigma }_1,{\sigma }_2,\dots ,{\sigma }_r)$, $\mathit{U}$ and $\mathit{V}$ represent the left and right feature vector. The singular value of ${\sigma }_1,{\sigma }_2,\dots ,{\sigma }_r$ can be used to identify the useful signals and noise, respectively. If the singular value is arranged in accordance with the descending order, the larger singular values mainly reflect the useful signal and the smaller singular values mainly reflects the noise. The rank is estimated by the eigenvalues of the neighborhood matrix calculated by SVD in the conventional TSSA method. It is common fact that the result provided by SVD is susceptible to noise especially under the strong noise environment. Without loss of generality, it can be considered the problem of estimating a low-rank matrix $\mathit{X}$ from its noisy observation $\mathit{Y}$.

$$\mathit{Y}=\mathit{X}+\mathit{G},\mathit{X},\mathit{Y},\mathit{G}\in R^{m\times n}$$

(7)

where $\mathit{Y}$ represents a noisy neighborhood matrix, $\mathit{X}$ is the expected low rank matrix, and $\mathit{G}$ is the zero-mean additive white Gaussian noise matrix. Note that for unitary matrices $\mathit{U}$ and $\mathit{V}$, $\Phi (\mathit{X})=\Phi (\mathit{U}\mathit{X}\mathit{V})$. Using the unitary invariant property of the Frobenius norm, we can define a typical low-rank matrix approximation (LRMA) issue as follows:

$$\begin{array}{ll}\stackrel{\wedge }{\mathit{X}}& =\arg \underset{\mathit{X}}{\min }\left\{\frac{1}{2}{\Vert \mathit{Y}-\mathit{X}\Vert }_F^2+\lambda \Phi (\mathit{X})\right\}\\ & =\arg \underset{\mathit{X}}{\min }\left\{\frac{1}{2}{\Vert \mathit{\Sigma }-{\mathit{U}}^{\mathrm{T}}\mathit{X}\mathit{V}\Vert }_F^2+\lambda \Phi ({\mathit{U}}^{\mathrm{T}}\mathit{X}\mathit{V})\right\}\\ & =\mathit{U}\cdot \arg \underset{\mathit{X}}{\min }\left\{\frac{1}{2}{\Vert \mathit{\Sigma }-\mathit{X}\Vert }_F^2+\lambda \Phi (\mathit{X})\right\}\cdot {\mathit{V}}^{\mathit{T}}\end{array}$$

(8)

where $\Phi (\mathit{X})={\displaystyle \sum _{i=1}^k\varphi ({\sigma }_i(\mathit{X});a)}$, ${\sigma }_i(\mathit{X})$ is the singular value of matrix $\mathit{X}$, $k=\min (m,n)$, $0\le a<1/\lambda$, $\varphi$ is a sparsity promotion function. When it meets $\varphi (x)=\left|x\right|$, LRMA issue becomes a typical nuclear norm minimization problem (NNM). Note that $\Theta (\mathit{\Sigma };\lambda ,a)=\arg \underset{\mathit{X}}{\min }\left\{\frac{1}{2}{\Vert \mathit{\Sigma }-\mathit{X}\Vert }_F^2+\lambda \Phi (\mathit{X})\right\}$. Therefore, the Equation (8) can be converted to:

$$\stackrel{\wedge }{\mathit{X}}=\mathit{U}\cdot \Theta (\mathit{\Sigma };\lambda ,a)\cdot {\mathit{V}}^{\mathrm{T}}$$

(9)

It is proved that parametric non-convex penalty function can be more accurate in estimating non-zero singular value than the nuclear norm under strong noise background [29]. Therefore, we use the non-convex penalty function defined by Equation (10) to achieve the optimized estimation of singular value.

$$\Theta (\mathit{\Sigma };\lambda ,a)=\min \left\{\left|\mathit{\Sigma }\right|,\max \left\{\left(\left|\mathit{\Sigma }\right|-\lambda \right)/(1-a\lambda ),0\right\}\right\}\mathrm{sign}(\mathit{\Sigma })$$

(10)

In this paper, the number of non-zero singular values obtained by the convex optimization framework is equivalent to the rank $R$.

2.3. The Desired Tensor Selection Based on PE

After the above tensor decomposition and selecting the proper rank of the desired source signal, it is necessary to return it to single channel data, which is likely to be our desired signal components. The desired single channel data from the tensor is obtained by the same hankelisation procedure. Firstly, it is done by hankelizing matrices across the first slab or block hankelizing the tensor. Then, the obtained tensor is hankelized again. Finally, the one-dimensional signal is reconstructed by one-to-one correspondence. It is worth pointing out that the desired components are not a simple superposition of these reconstruction signal. Thus, an effective way to accurately evaluate the useful components should be taken into consideration. Permutation Entropy (PE) is one of the most effective ways to measure the randomness and dynamic mutation of time sequence. Similar to the above operations defined in Section 2.2, the trajectory matrix corresponding to the reconstruction signal is obtained by phase space reconstruction. The matrix in each row is rearranged in ascending order and $j_1,j_2,\dots ,j_m$ indicates the index of each element in the column of reconstruction components. A set of symbols sequence can be obtained from each row in the matrix.

$$S(l)=\left\{j_1,j_2,\dots ,j_m\right\}$$

(11)

where l = 1, 2,…, k, and k ≤ m. The probability of each symbol sequence was calculated and named as P₁, P₂,…, P_k. The value of PE corresponding to different symbol sequences in reconstruction signal is defined as the form of information entropy as follows:

$$H_p(d)={\displaystyle \sum _{j=1}^k{P}_j}\ln (P_j)$$

(12)

It is also noted that the value of $H_p(d)$ represents the degree of randomness of time series, namely, the larger value of $H_p(d)$ declares the more randomness of time series, and vice versa [32]. Commonly, the permutation entropy of harmonic signal and the modulation signal is small, while the permutation entropy for random noise is large. Hence, PE can be used to select the desired tensor for signal reconstruction. Ultimately, the flowchart of the proposed method is shown in Figure 3.

3. Test with Numerical Simulation Signal

3.1. The Performance of Abnormal Signal Detection Using PE

In the process of permutation entropy calculation, there are two parameter values that need to be determined: the embedding dimension $m$ and delay time $\tau$. Under the condition of the above parameters, the following representative signal is considered:

$$x_1=[1+0.7\sin (2\pi f_{1}t)]\sin (2\pi f_{2}t)$$

(13)

$$x_2=0.8\sin (2\pi f_{3}t)$$

(14)

where $f_1=5$, $f_2=10$, $f_3=15$. $x_3$ is the standard Gaussian white noise. Theoretically, the embedding dimension and delay time are chosen as $m=6$ and $\tau =1$ according to the False Nearest Neighbor algorithm (FNN) [38]. According to the definition of PE, the normalized PE values of three simulated signals are calculated as 0.2909, 0.4451 and 0.9644, respectively. The result demonstrated that the PE value corresponding to noisy signal is obviously lager than normal modulating signal and harmonic signal. It is proved that PE is sensitive to noisy signal and can be used to detect the abnormal signal. For the proposed method, the value of PE corresponding to the undesired tensor by hankelisation procedure is greater than the threshold value, which is defined as 0.6.

3.2. The Feature Extraction Result Provided by Proposed Method

Rolling bearings are mainly used to support the rotating part in the mechanical equipment and its vibration signal always contains a lot of information such as fault characteristic and noise component. The key step in fault diagnosis is to extract vibration signal feature in frequency domain. There are many kinds of the simulation signal model about bearing fault and the most typical one is proposed by Randall [39,40]. Without loss of generality, the numerical simulation signal is expressed as follows:

$$x(t)={\displaystyle \sum _{i=0}^M[A_0\cdot \cos (2\pi f_{m}t+{\varphi }_A)+C_A]\cdot e^{-B(t-i\mathrm{T}-{\tau }_i)}\cdot }\cos [2\pi f_n(t-i\mathrm{T}-{\tau }_i)+{\varphi }_w]+n(t)$$

(15)

where $A_0$ is the amplitude of resonance; $f_m$ is the modulation frequency; ${\varphi }_A$, ${\varphi }_w$, and $C_A$ are selected as arbitrary constants; $B$ is the attenuation coefficient; T is defined as the average time between two impacts with $\mathrm{T}=1/f_p$; $f_p$ is fault character frequency; ${\tau }_i$ is regarded as the time lag from its mean period due to the presence of slip; $f_n$ is the resonance frequencies of the bearing system; and $n(t)$ is the additive white Gaussian noise. The three common types of faults on outer race, inner race and rolling elements can be modeled as $f_m=0$, $f_m=f_r$, $f_m=f_{re}$, respectively. It should also be noted that $f_r$ and $f_{re}$ are denoted as rotational frequency and retainer frequency. The fault character frequencies of three common types as f_i, f_o, f_ro are described in

Table 1. The fault character frequency of three common types.

Inner Race fi (Hz)	Outer Race fo (Hz)	Rolling Element fro (Hz)
156	103.9	59.7

. Firstly, the outer race fault feature extraction under noise environment is researched and the parameter selection is listed as

Table 2. The parameter selection of outer race fault simulation signal.

${\mathit{A}}_0$	${\mathit{f}}_{\mathit{r}}$ (Hz)	${\mathit{f}}_{\mathit{m}}$ (Hz)	${\mathit{\tau }}_{\mathit{i}}$	${\mathit{f}}_{\mathit{n}}$ (Hz)	${\mathit{C}}_{\mathit{A}}$	${\mathit{\varphi }}_{\mathit{A}}$	${\mathit{\varphi }}_{\mathit{w}}$	$\mathit{B}$
0.0003	17	0	0.01	2000	1	0	0	800

. Moreover, $n(t)$ is Gaussian white noise with variance 0.5 in the section of simulation signal analysis.

The sampling frequency and the sampling points are set as 10,000 Hz and 2000, respectively. Original signal without noise and noisy synthesis simulation signal in time-domain are shown in Figure 4a,b, respectively. The frequency spectrum about $x(t)$ is conducted in Figure 5.

In all the simulations and experiments, the value of $L$ and $J$ is set to 200 and 50, respectively. Commonly, the delay time $\tau$ is selected as 1 in order to keep the computational complexity low. Then, the parameter of $I$ and $K$ can be obtained by the relation of $I=n/L$, $L=(K-1)\times \tau +J$. In addition, the rank of tensor is determined by the convex optimization framework introduced in Section 2.2, and the desired tensor selection is performed by the index of PE in Section 3.1. According to Figure 4b, it is acceptable that the outer fault characteristics of the original simulation signal in the time domain cannot be clearly indicated in the strong background noise. For the frequency analysis in Figure 5, the outer race fault characteristics frequency is still concealed by noisy components.

In the simulation signal analysis section, in order to accurately evaluate the proposed method in signal reconstruction under noisy conditions, the proposed method, the conventional TSSA method and wavelet transform (WT) are used to the comparative analysis. Figure 6 is the result provided by conventional TSSA model and only the third harmonic ($3f_o$) can be inspected. It is important to note that the “wden” function is chosen as the automatic 1-D denoising function and the wavelet basis function is set as “sym3” for wavelet denoising in the rest of this paper. Figure 7 is a diagram of wavelet denoising in time domain and frequency domain. It is obvious that both the conventional TSSA and WT cannot completely identify the outer race fault characteristic frequency. Figure 8 is the result provided by the proposed method based on PE in frequency domain, which indicates the outer race feature frequency $f_o$ and its frequency multiplication such as $2f_o$, $3f_o$, $4f_o$ can be identified. Then, the proposed method is applied to characteristic frequency extraction for inner race fault and rolling elements fault. The result provided by the proposed method is drawn in Figure 9 and the modulation phenomenon is in keeping with the inner ring and rolling elements fault characteristics. Thus, we can make a conclusion that the proposed method has better performance in denoising and characteristic frequency extraction.

It can be seen from the above analysis result provided by different methods that conventional TSSA and WT cannot completely separate the noisy component from the original signal. However, the proposed tensor singular spectrum algorithm based on PE is more effective in removing noise component, better in reflecting the fault characteristics of the original signal. Through multiple simulations, the improved TSSA method has higher convergence and better reconstruction precision than TSSA.

For a more detailed description of the degree of noise reduction, the paper selects signal-to-noise ratio (SNR) as an index for evaluation, which is critical for assessing the degrees of noise contamination. The SNR of the noisy signal about outer race fault simulation was varied between 2.17 dB and 0.32 dB, while the corresponding reconstruction SNR of the simulation signal was then calculated and the result is plotted in the following Figure 10. The smaller values always indicate the poor reconstruction accuracy. As can be seen in Figure 10, the results also indicate that the proposed method performs well in denoising and signal reconstruction.

4. Applications to Rolling Bearing Fault Feature Extraction

In industrial practice, the measured vibration signal is more complicated than the simulated signal. In order to verify the effectiveness of this approach in the real environment, the experimental device analog failure signal was measured to be analyzed. Simulation apparatus is shown in Figure 11, where bench consists of a 550 W (220 V/50 Hz) AC motor driven drive shaft through the coupling operation. The electric spark machining method was used to carry out pitting treatment on the outer ring of replaceable bearing to simulate the outer faults and the depth of pitting corrosion pit is close to 1.0 mm. The acceleration signal of the experiment was collected in the vertical direction of the bearing on the right side of the experimental platform using a CSI2130 data analyzer (Emerson, Sao Luis, Brazil) and vibration acceleration sensor PCB-352C33 (PCB, New York, NY, USA). The bearing type is deep groove ball bearing with model number of 6207. The medium diameter of removable bearing is D = 53.5 mm, the number of rolling elements is Z = 9, and the roller diameter is d = 11.1 mm. The detailed parameters of the experiment are shown in

Table 3. The experimental parameters and fault frequency.

Rotating Speed r/min	Rotating Frequency/Hz	Sampling Frequency/Hz	Sampling Time/s	Outer Fault Frequency/Hz
1450	24.17	16384	1	87.01

.

It is worth pointing out that the frequency of outer ring fault signal of rolling bearing is $f_o=87.01$ Hz by calculating, and the rotation frequency is $f_r=24.17$ Hz. The original signal in time domain and frequency domain is shown in Figure 12a,b, respectively.

From time-domain waveform in Figure 12a,b, it cannot be determined where the failure occurred. Meanwhile, in order to illustrate the effectiveness the proposed method, the conventional noise reduction method such as wavelet denoising algorithm is applied and the results are shown in Figure 13.

After the wavelet denoising, only the three multiple of outer fault feature frequency $3f_o$ can be inspected, while the other useful components are still difficult to be accurately identified.

Based on the False Nearest Neighbor algorithm (FNN) [38], the embedding dimension and delay time is determined as $m=12$ and $\tau =1$ in the experimental case studies. Finally, the proposed method based on PE is utilized to signal denoising and feature extraction, and the result is shown in Figure 14.

Through the tensor decomposition using PE, the unwanted signal components have been deducted and the feature information can be enhanced. It is obvious that the outer ring fault frequency $f_o$ and its harmonics ($2f_o$, $3f_o$) all can be identified in Figure 14. Hence, we can make a conclusion that the outer ring of the rolling bearing has occurred failure, which is consistent with the actual situation.

5. Conclusions

This paper improves the tensor singular spectrum analysis (TSSA) algorithm to address its problems of inestimable tensor rank and lacking in criterion to achieve the desired tensor selection. The main work of this paper was reflected in the following aspects: (1) the accurate estimation of the rank of tensor was solved by SVD using convex optimization and non-convex penalty function; (2) the desired tensor selection using PE was conducted to improve the denoising performance of the proposed method; and (3) in the processing of measured rolling bearing fault signal, the proposed improved TSSA method based on PE outperformed other methods in extracting weak fault characteristics, which demonstrated that the proposed method was feasible in the future real applications.