Trivariate Empirical Mode Decomposition via Convex Optimization for Rolling Bearing Condition Identification

As a multichannel signal processing method based on data-driven, multivariate empirical mode decomposition (MEMD) has attracted much attention due to its potential ability in self-adaption and multi-scale decomposition for multivariate data. Commonly, the uniform projection scheme on a hypersphere is used to estimate the local mean. However, the unbalanced data distribution in high-dimensional space often conflicts with the uniform samples and its performance is sensitive to the noise components. Considering the common fact that the vibration signal is generated by three sensors located in different measuring positions in the domain of the structural health monitoring for the key equipment, thus a novel trivariate empirical mode decomposition via convex optimization was proposed for rolling bearing condition identification in this paper. For the trivariate data matrix, the low-rank matrix approximation via convex optimization was firstly conducted to achieve the denoising. It is worthy to note that the non-convex penalty function as a regularization term is introduced to enhance the performance. Moreover, the non-uniform sample scheme was determined by applying singular value decomposition (SVD) to the obtained low-rank trivariate data and then the approach used in conventional MEMD algorithm was employed to estimate the local mean. Numerical examples of synthetic defined by the fault model and real data generated by the fault rolling bearing on the experimental bench are provided to demonstrate the fruitful applications of the proposed method.


Introduction
With the development of industrial production, the key equipment and structure always have the characteristics of complex structure, variable operation status and continuous online service. Especially in the metallurgical industry, most parts such as rolling bearing are working under an environment of high speed, overloading and high temperatures. As an important approach of structural health monitoring, the damage feature extraction algorithm for rolling bearings based on signal processing technology has important significance [1][2][3]. Naturally, the vibration signal is expressed as nonlinear, non-stationary and time-varying [4][5][6][7][8]. Moreover, the noise components directly affect the performance of the signal analysis [9][10][11]. Thus, improved signal processing methods are required.
In order to improve the signals to noise ratio (SNR) of actual measured vibration signal, many scholars have done a lot of research work [12][13][14][15][16]. Short-time Fourier transform (STFT) to verify the validity of the method, the simulation signal model and measured experimental data of fault rolling bearing are conducted to analysis. This paper was organized as follows: the basic principles and characteristics of the trivariate empirical mode decomposition via convex optimization were introduced in the second chapter. In the third section, the numerical simulation analysis was conducted to confirm the effectiveness of the algorithm. In the fourth part, the validity of the proposed method was verified by the bearing data of the experimental bench. The conclusions of the study were given in section fifth.

Theoretical Descriptions
Given a trivariate signal Y(t) = {y 1 (t), y 2 (t), y 3 (t)} generated by the mechanical equipment using triaxial acceleration sensor or three sensors installed in different positions, it is a common fact that the actual measured trivariate signal is composed by useful signal and noise components. The problem of estimating a low-rank matrix X from its noisy observation Y is considered as follows: where W represents zero-mean additive white Gaussian noise. It was reported that the non-convex regularization methods outperformed the convex regularization methods on matrix approximation and could enhance sparsity [32]. Thus, in this paper, a low-rank matrix approximation [33] for the trivariate signal Y can be expressed as the following model: where k = min(m, n), φ is the parameterized non-convex penalty function, σ i (X) represents the singular values of the matrix X. a 0 and a 1 represent the parameters of non-convex penalty functions, which are associated with the convexity of the objective function. The weights of two parameterized non-convex penalty functions are adjusted by the regularization parameters λ 0 and λ 1 . Obviously, if λ 1 = 0 and φ(x, a) = |x|, Equation (2) is transformed into a typical nuclear norm minimization problem. Typical non-convex penalty functions should be 2nd order differentiable, continuous and symmetric. An example of a non-convex penalty function is the rational penalty function, which is defined as: Equation (3) is parameterized by the parameter a ≥ 10. The proximal operator of φ(x, a) is expressed as: It is noted that the proximity operator is associated with the non-convex function φ(x, a). In order to ensure that the objective function ψ(X) is strictly convex, the parameter a 0 and a 1 should meet the following requirement: Then, the alternating direction method of multipliers (ADMM) algorithm in conjunction with variable splitting is employed to solve the trivariate data matrix of the low-rank approximation problem. The ADMM algorithm is applied to Equation (2), and yields the following iterative procedure. Firstly, initialize the model with Z = 0 and D = 0. Then, the iterative optimization is conducted as following: where µ > 1 is the scalar augmented Lagrangian parameter to ensure the strictly convex of the algorithm; SVD indicates the operator of singular value decomposition. The trivariate data preprocessing method based on low-rank matrix approximation via convex optimization requires the specification of two regularization parameters λ 0 and λ 1 , two penalty parameters a 0 and a 1 , and the scalar augmented Lagrangian parameter µ. Generally, the regularization parameters λ 0 and λ 1 are determined by the follow formulation: where β i is chosen so as to maximize the signal-to-noise ratio. σ is corresponding to the spectral noise variance, which can be estimated by the wavelet coefficients by using the local variance analysis [34]. The values of two penalty parameters a 0 and a 1 can be simplified chosen by Equation (5). It's also noted that the value of µ does not affect the solution due to the iterative algorithm converges, and µ = 1.5 is always suitable for algorithm convergence. Through the processing of the proposed method, a low-rank approximation matrix X is obtained from the observed spectral data matrix X. Then, the singular value decomposition (SVD) of the covariance matrix to X is performed, E XX T = VΛV T , where V is the eigenvector matrix and Λ = diag{λ 1 , λ 2 , λ 3 } is the eigenvalues corresponding to the eigenvalue matrix. The corresponding azimuth angle φ H and inclination angle θ H of the Hammerseley projections are obtained by uniformly sampling a sphere using the Hammerseley sequence [29]. It is aimed to identify the Cartesian coordinates of the uniformly sampled sphere. Thus, the non-uniform sample e p on a sphere can be determined by constructing an ellipsoid with the following parameters: The direction of the highest curvature is sampled by rotating the ellipsoid e l . Therefore, the novel non-uniform samples scheme can be defined as e p = Ve l = {e 1 , e 2 , e 3 }. As a critical step of signal decomposition, the local mean estimation according to the conventional MEMD algorithm using the non-uniform samples e p is implemented and the computational procedure of trivariate empirical mode decomposition via convex optimization can be described in Algorithm 1. Perform the low-rank matrix approximation via convex optimization framework to a trivariate signal Y(t) and the new observed signal X(t) can be obtained.

2.
The projection p k (t)(k = 1, 2, . . . , K) can be calculated of the input low-rank trivariate signal X(t) along the direction vector e 1 k , e 2 k , e 3 k . It should also be noted that K is the number of the directional vector sets.

3.
The time instants t k m corresponding to the maxima of the set of projected signals p k (t) is determined.

4.
Interpolate [t k m , X(t k m )] to obtain multivariate envelope curves E k (t). Then, the envelop mean can be calculated M(t) = 1

5.
Calculate the residual by R(t) = X(t) − M(t). If the stopping criterion condition of iteration can be satisfied, then R(t) is set as one IMF and repeat the above steps to X(t) − R(t) until the next IMF is isolated. If it does not satisfy the stopping criterion, then repeat the above steps to R(t) until it meets the criterion. Cauchy-type convergence condition is chosen as the stopping criterion of the sifting iterative process, which indicates the iterative will stop under the circumstance that the discrepancy between adjacent sifting results is less than a threshold with a range 0.2-0.3. For clarify, K = 32 is always suitable for the application of trivariate signal empirical mode decomposition.

Simulation Signal Analysis
Rolling bearing is an important transmission component supporting the rotating part of the mechanical equipment. If a bearing with an outer race that is fixed to the structure, a typical fault rolling bearing model can be simplified as follows [35]: where f i is the inner ring failure frequency and f r is corresponding to rotational frequency, respectively. Then, the following three source signals X = {x 1 , x 2 , x 3 } are employed to simulate the collected vibration signals in this section: where the characteristics frequency about three simulation source signals are chosen as f 1 = 10 Hz, f 2 = 45 Hz, f i = 110 Hz, f r = 10 Hz. The sampling number and sampling frequency are set as N = 1024 and f s = 1024 Hz, respectively. In the analog process of sensor collecting signal, the vibration signal will be collected by any sensor simultaneously, and the trivariate signal is the instantaneous mixed signal composed by three above mentioned simulation signal. Considering the influence of noise, the zero-mean additive white Gaussian noise with variance 0.5 is added to the observation signal, and it is denoted as S N . A random matrix of 3 × 3 is chosen to mix the simulation source signal together and it can be expressed as follows: Then, the instantaneous mixed signal model can be described as Y = AX + S N and Y = {y 1 , y 2 , y 3 }, which indicates trivariate signal with noisy generated by three source signal using three sensors. The analysis results of trivariate signal in time-domain and frequency spectrum provided by FFT are shown in Figures 1 and 2.
Then, the instantaneous mixed signal model can be described as = + N Y AX S and , which indicates trivariate signal with noisy generated by three source signal using three sensors. The analysis results of trivariate signal in time-domain and frequency spectrum provided by FFT are shown in Figures 1 and 2.  As shown in Figures 1 and 2, it is hard to identify the frequency characteristics and modulation phenomenon of the observed trivariate signals with noisy. Only the inner race fault frequency i f in channel#3 can be inspected, and other frequency components are not obvious for all the channels. It is illustrated that the influence of noise should not be neglected and more advanced methods are required. Then, the conventional MEMD methods are employed to deal with the simulation signal and the result is plotted in Figure 3. From the analysis result of the three channels shown in Figure 3, it is proved that MEMD has advantages in guaranteeing the stability of the decomposition results.  , , y y y = Y , which indicates trivariate signal with noisy generated by three source signal using three sensors. The analysis results of trivariate signal in time-domain and frequency spectrum provided by FFT are shown in Figures 1 and 2.  As shown in Figures 1 and 2, it is hard to identify the frequency characteristics and modulation phenomenon of the observed trivariate signals with noisy. Only the inner race fault frequency i f in channel#3 can be inspected, and other frequency components are not obvious for all the channels. It is illustrated that the influence of noise should not be neglected and more advanced methods are required. Then, the conventional MEMD methods are employed to deal with the simulation signal and the result is plotted in Figure 3. From the analysis result of the three channels shown in Figure 3, it is proved that MEMD has advantages in guaranteeing the stability of the decomposition results. As shown in Figures 1 and 2, it is hard to identify the frequency characteristics and modulation phenomenon of the observed trivariate signals with noisy. Only the inner race fault frequency f i in channel#3 can be inspected, and other frequency components are not obvious for all the channels. It is illustrated that the influence of noise should not be neglected and more advanced methods are required. Then, the conventional MEMD methods are employed to deal with the simulation signal and the result is plotted in Figure 3. From the analysis result of the three channels shown in Figure 3, it is proved that MEMD has advantages in guaranteeing the stability of the decomposition results. However, the characteristic frequencies f 1 , f 2 and the frequency modulation f i ± f r still cannot be identified. However, the characteristic frequencies 1 f , 2 f and the frequency modulation i r f f ± still cannot be identified. Eventually, the proposed method is applied to simulate trivariate signal and the result is drawn in Figure 4. For clarify the denosing performance, the 3th, 4th, and 7th IMFs are chosen according to the maximum similarity with the original signal shown in Figure 5. It is obvious that the frequency modulation i r f f ± can be found in the 3th IMF. In addition, the characteristic frequencies 1 f , 2 f also can be viewed in the 4th and 7th IMF, respectively. The result provided by the proposed method is coincides with practical situation and the effective can be demonstrated.  Eventually, the proposed method is applied to simulate trivariate signal and the result is drawn in Figure 4. For clarify the denosing performance, the 3th, 4th, and 7th IMFs are chosen according to the maximum similarity with the original signal shown in Figure 5. It is obvious that the frequency modulation f i ± f r can be found in the 3th IMF. In addition, the characteristic frequencies f 1 , f 2 also can be viewed in the 4th and 7th IMF, respectively. The result provided by the proposed method is coincides with practical situation and the effective can be demonstrated.  Eventually, the proposed method is applied to simulate trivariate signal and the result is drawn in Figure 4. For clarify the denosing performance, the 3th, 4th, and 7th IMFs are chosen according to the maximum similarity with the original signal shown in Figure 5. It is obvious that the frequency modulation i r f f ± can be found in the 3th IMF. In addition, the characteristic frequencies 1 f , 2 f also can be viewed in the 4th and 7th IMF, respectively. The result provided by the proposed method is coincides with practical situation and the effective can be demonstrated.

Analysis of the Roll Bearing on the Experimental Bench
Commonly, the measured vibration signal is more complicated than the simulation signal in the process of actual signal analysis. In order to verify the effectiveness of the proposed algorithm of trivariate empirical mode decomposition via convex optimization, the outer ring failure signal of rolling bearing in the experimental device is analyzed. The experimental apparatus and its structure diagram are drawn in Figure 6, where bench comprises a drive shaft which is driven by 550 W (220 V/50 Hz) AC motor. The electric spark machining method is used to carry out pitting treatment on the outer ring of the replaceable bearing to simulate the outer faults. A PCB-352C33 acceleration sensor is used to

Analysis of the Roll Bearing on the Experimental Bench
Commonly, the measured vibration signal is more complicated than the simulation signal in the process of actual signal analysis. In order to verify the effectiveness of the proposed algorithm of trivariate empirical mode decomposition via convex optimization, the outer ring failure signal of rolling bearing in the experimental device is analyzed. The experimental apparatus and its structure diagram are drawn in Figure 6, where bench comprises a drive shaft which is driven by 550 W (220 V/50 Hz) AC motor.

Analysis of the Roll Bearing on the Experimental Bench
Commonly, the measured vibration signal is more complicated than the simulation signal in the process of actual signal analysis. In order to verify the effectiveness of the proposed algorithm of trivariate empirical mode decomposition via convex optimization, the outer ring failure signal of rolling bearing in the experimental device is analyzed. The experimental apparatus and its structure diagram are drawn in Figure 6, where bench comprises a drive shaft which is driven by 550 W (220 V/50 Hz) AC motor. The electric spark machining method is used to carry out pitting treatment on the outer ring of the replaceable bearing to simulate the outer faults. A PCB-352C33 acceleration sensor is used to  The electric spark machining method is used to carry out pitting treatment on the outer ring of the replaceable bearing to simulate the outer faults. A PCB-352C33 acceleration sensor is used to collect the acceleration signals. The installation position of sensor is located in the horizontal direction, vertical direction and axial direction of the bearing on the right side of the experimental platform. The bearing type is deep groove ball bearing with model number of 6207. The external diameter of removable bearing was D = 72 mm and the inner diameter was d = 35 mm, respectively. The testing flow chart for experimental bench is demonstrated in Figure 7 and the sampling frequency is determined as 16,000 Hz. Vibration detection in horizontal (x), vertical (y) and axial (z) directions should be considered as far as possible, which is aimed to collect the trivariate signal. Suppose there is no relative sliding between the raceway surface and the rolling elements, and the outer ring is fixed. Then, the outer fault frequency 87.01 Hz can be obtained. The detailed experimental parameters are shown in Table 1.  Figure 7 and the sampling frequency is determined as 16,000 Hz. Vibration detection in horizontal (x), vertical (y) and axial (z) directions should be considered as far as possible, which is aimed to collect the trivariate signal. Suppose there is no relative sliding between the raceway surface and the rolling elements, and the outer ring is fixed. Then, the outer fault frequency 87.01 Hz can be obtained. The detailed experimental parameters are shown in Table 1 Figures 8 and 9, respectively.  It can be seen from Table 1 that the rotation frequency is f r = 24.17 HZ. According to the theory of bearing fault diagnosis, the frequency of outer ring fault signal of rolling bearing is calculated as f o = 87.01 HZ. Then, the time domain and frequency domain plots of the original measured vibration signal are shown in Figures 8 and 9, respectively. Sensors 2018, 18, x FOR PEER REVIEW 10 of 14  The spectral analysis by FFT is applied to the measured trivariate signal in Figure 9 and the result shows that the signal components are complex. In addition, the fault characteristic frequency and frequency modulation phenomenon cannot be identified. Then, the conventional MEMD method is used to analysis the fault signal and the result in frequency-domain is plotted in Figure 10. It is also noted that all the analyses employed 32 projection vectors to capture the direction of highest curvature of trivariate signals. By observing the frequency features of each IMFs, the fault characteristic frequency still cannot be inspected.  The spectral analysis by FFT is applied to the measured trivariate signal in Figure 9 and the result shows that the signal components are complex. In addition, the fault characteristic frequency and frequency modulation phenomenon cannot be identified. Then, the conventional MEMD method is used to analysis the fault signal and the result in frequency-domain is plotted in Figure 10. It is also noted that all the analyses employed 32 projection vectors to capture the direction of highest curvature of trivariate signals. By observing the frequency features of each IMFs, the fault characteristic frequency still cannot be inspected. The spectral analysis by FFT is applied to the measured trivariate signal in Figure 9 and the result shows that the signal components are complex. In addition, the fault characteristic frequency and frequency modulation phenomenon cannot be identified. Then, the conventional MEMD method is used to analysis the fault signal and the result in frequency-domain is plotted in Figure 10. It is also noted that all the analyses employed 32 projection vectors to capture the direction of highest curvature of trivariate signals. By observing the frequency features of each IMFs, the fault characteristic frequency still cannot be inspected. Finally, the proposed method is employed to the measured trivariate vibration signal and the result in time-domain is shown in Figure 11. For the sake of clarity, the FFT operators to the IMF3 for all the channels are conducted, and the corresponding results are plotted in Figure 12 can also be identified. Thus, we can firmly believe that the fault is located on the outer ring, which is consistent with the actual situation shown in Figure 13.  Finally, the proposed method is employed to the measured trivariate vibration signal and the result in time-domain is shown in Figure 11. For the sake of clarity, the FFT operators to the IMF3 for all the channels are conducted, and the corresponding results are plotted in Figure 12. Fortunately, the outer ring fault frequency f o and its doubling 2 f o can be found. Additionally, the special frequency modulation phenomenon of fault rolling bearing such as f o ± f r and 2 f o ± f r can also be identified. Thus, we can firmly believe that the fault is located on the outer ring, which is consistent with the actual situation shown in Figure 13. Finally, the proposed method is employed to the measured trivariate vibration signal and the result in time-domain is shown in Figure 11. For the sake of clarity, the FFT operators to the IMF3 for all the channels are conducted, and the corresponding results are plotted in Figure 12 can also be identified. Thus, we can firmly believe that the fault is located on the outer ring, which is consistent with the actual situation shown in Figure 13.   By comparing the results respectively provided by MEMD and the method proposed in this paper, we can make a conclusion that the proposed method has better performance in trivariate signal mode decomposition and fault feature extraction.

Conclusions
Based on the traditional MEMD method, a novel method of trivariate empirical mode decomposition via convex optimization for rolling bearing condition identification is proposed and tested in this paper. The main research works are listed as follows: (1) a low-rank matrix approximation via convex optimization framework is proposed for the trivariate signal denoising. By introducing parameterized non-convex function, its performance is improved; (2) the principal component obtained by SVD operator is employed to indicate the sample directions, thus the non-uniform sample scheme is presented to meet the requirement of inhomogeneity of the data distribution in high-dimensional space; (3) through the analysis of the simulation signal and the measured vibration signal, it is demonstrated that the proposed method is superior to MEMD in multi-scale feature extraction. It can be regarded as powerful tool in trivariate signal processing.   By comparing the results respectively provided by MEMD and the method proposed in this paper, we can make a conclusion that the proposed method has better performance in trivariate signal mode decomposition and fault feature extraction.

Conclusions
Based on the traditional MEMD method, a novel method of trivariate empirical mode decomposition via convex optimization for rolling bearing condition identification is proposed and tested in this paper. The main research works are listed as follows: (1) a low-rank matrix approximation via convex optimization framework is proposed for the trivariate signal denoising. By introducing parameterized non-convex function, its performance is improved; (2) the principal component obtained by SVD operator is employed to indicate the sample directions, thus the non-uniform sample scheme is presented to meet the requirement of inhomogeneity of the data distribution in high-dimensional space; (3) through the analysis of the simulation signal and the measured vibration signal, it is demonstrated that the proposed method is superior to MEMD in multi-scale feature extraction. It can be regarded as powerful tool in trivariate signal processing. By comparing the results respectively provided by MEMD and the method proposed in this paper, we can make a conclusion that the proposed method has better performance in trivariate signal mode decomposition and fault feature extraction.

Conclusions
Based on the traditional MEMD method, a novel method of trivariate empirical mode decomposition via convex optimization for rolling bearing condition identification is proposed and tested in this paper. The main research works are listed as follows: (1) a low-rank matrix approximation via convex optimization framework is proposed for the trivariate signal denoising. By introducing parameterized non-convex function, its performance is improved; (2) the principal component obtained by SVD operator is employed to indicate the sample directions, thus the non-uniform sample scheme is presented to meet the requirement of inhomogeneity of the data distribution in high-dimensional space; (3) through the analysis of the simulation signal and the measured vibration signal, it is demonstrated that the proposed method is superior to MEMD in multi-scale feature extraction. It can be regarded as powerful tool in trivariate signal processing.
Author Contributions: Y.L. and C.Y. conceived and designed the experiments; Y.L. and H.Z. performed the experiments; H.Z. and C.Y. analyzed the data; C.Y. contributed to materials/analysis tools development; C.Y. and Y.L. wrote the paper.