# A Fault Diagnosis Method for One-Dimensional Vibration Signal Based on Multiresolution tlsDMD and Approximate Entropy

Authors: Zhang Dang; Yong Lv; Yourong Li; Guoqian Wei
Journal: Shock and Vibration
 (2019)

---

## Abstract

Dynamic mode decomposition (DMD) has certain advantages compared with the traditional fault signal diagnosis method. By exploiting the strength of DMD algorithm in signal processing, this paper proposes a joint fault diagnosis scheme to extract the spatial and temporal patterns and evaluate them for the complexity to diagnose the fault for one-dimensional mechanical signal. The multiscale method is adopted to decompose the reconstructed matrix of standard DMD modes into multiple scales with a given level parameter. Total least squares DMD algorithm is performed on each level to solve the noise sensitivity problem. Approximate entropy (ApEn) is performed on the grouped multiscale spatiotemporal modes that represent the dynamic characteristic information of the original signal. ApEn values are used as a fault recognizer to identify fault types. By applying the algorithm on three experimental mechanical vibration data, we verify the effectiveness of the proposed method. The result demonstrates that the proposed scheme can effectively recognize different fault forms as a fault diagnosis method.


---

## Body

## 1. Introduction

Fault diagnosis, an important technology to promote the development of national advanced manufacturing level, plays a vital role in structural health monitoring, service life evaluation, product quality assurance, and production efficiency improvement [1]. Mechanical equipment fault has incalculable potential danger for safe and reliable operation of large continuous rolling mill, aero-engine, gas turbine, and other large mechanical equipment. Due to the characteristics of complex structure, variable operating conditions, and long-term online service, early fault signals usually have the characteristics of nonlinear, nonstationary, non-Gaussian, potentiality, weak dynamic response, and multifactor coupling [2, 3]. Generally, the collected amplitude-modulated and/or frequency-modulated (AM-FM) multicomponent signal is inevitably interfered with background noise. Signal analysis methods of one-dimensional vibration noisy signal play a decisive role in fault diagnosis.In the last few decades, researchers have done a lot of work on nonlinear and nonstationary fault signal diagnosis. The short-time Fourier transform (STFT), developed by Gabor [4], compensates the deficiency of Fourier transform (FT) and is capable of analyzing nonstationary signals. However, STFT cannot adapt to the change of time and frequency due to its constant sliding window [5]. Wavelet transform (WT), which was proposed by Mallat [6], overcomes the shortcoming of STFT. Regrettably, the effect of WT is closely related to the choice of wavelet basis function [7]. Wigner–ville distribution (WVD), introduced by Ville [8], shows good time-frequency aggregation in signal processing; however, there is a problem of cross-term interference for multicomponent signals [9]. Subsequently, in order to completely eliminate the limitations of the time-frequency analysis and to obtain an effective result in allusion to the local time and frequency characteristics of nonlinear and nonstationary signal, Huang et al. [10] proposed empirical mode decomposition (EMD), which has the advantages of orthogonality and completeness. While EMD has been widely applied in biomedical, mechanical fault diagnosis, seismic signal analysis, and other fields [11–13], there are still some imperfect problems which need to be solved and improved, such as over-envelope, under-envelope, and modal confusion and endpoint effect [14]. In response to the above problems, Smith [15] proposed local mean decomposition (LMD), which performs better than EMD in reducing the number of iterations and inhibiting the endpoint effect. Though LMD has been well applied in the field of mechanical fault diagnosis [16, 17], it has the imperfections of low efficiency and pattern confusion [16]. Variation mode decomposition (VMD), proposed by Dragomiretskiy and Zosso [18], decomposes the nonstationary signal into intrinsic mode functions (IMFs) adaptively and nonrecursively [19]. As a decomposition method based on spectrum segmentation [20], VMD overcomes the problem of mode aliasing and achieves more effective separation between the IMFs. The selection of initial parameters (e.g. moderate bandwidth constraint parameter, number of modes, and tolerance of convergence criterion) is the prior problem in the practical application of VMD [19]. Proper orthogonal decomposition (POD) [21], a method of the model order reduction technique which is closely related to principle component analysis (PCA), can extract modal contents from nonlinear dynamical systems by generating an orthonormal set of spatial basis functions. Although POD has been widely used in the field of fluid and PCA in the field of mechanical signal diagnosis, they are confronted with the defects of orthogonal basis [22].As a new nonlinear and nonstationary signal processing method, dynamic mode decomposition (DMD) has already gained fruitful achievements in the field of nonlinear dynamics, especially in the fluid field. In 2009, Schmid [23] firstly proposed the standard DMD (sDMD) method, which can extract the spatiotemporal coherent characteristics by decomposing the complex flow field signal into a series of simple expressions. DMD can decompose the time series into a series of single-frequency modal components [24], avoiding the modal aliasing problem in EMD and LMD. At the same time, because of its nature of being equation-free and data-driven [25], DMD has strict mathematical and theoretical foundation, avoiding the drawbacks of EMD. DMD does not depend on any prior assumptions [23] and avoids the advance selection of WT wavelet basis functions. In the algorithm implementation process of sDMD and VMD excepting the parameter of decomposition modes’ number, VMD needs to preset other parameters while sDMD does not. It is more meaningful in physics compared with traditional signal processing methods because DMD decomposes the nonlinear and nonstationary time series into nonorthogonal modal components corresponding to the eigenvalues, and the modal components are related to its inherent spatial and temporal patterns [26]. For a nonlinear and nonstationary system, there is a set of Koopman modes (kernel of DMD algorithm) that can fully characterize its dynamic characteristics [27]. Therefore, DMD avoids the FT’s drawbacks in characterizing the system by harmonic signals. With the advantage of describing the dynamic characteristics in a series of single nonorthogonal frequency modes [24], DMD maximizes the complexity of the original system, while POD decomposes data into a series of multiple orthogonal modal components, lacking the variability in structure and the concept of freedom in signal expression. As a whole, DMD has certain advantages compared with the traditional fault signal diagnosis method. By exploiting the strength of DMD algorithm in signal processing, this paper proposes an improved algorithm to extract the spatial and temporal patterns, and evaluate them the complexity to diagnose the fault type for one-dimensional mechanical signal.Since the sDMD algorithm was proposed, it has been widely applied in the fields of fluid domain [25–30], video processing [31, 32], biometrics [33, 34] and neuroscience [35, 36]. A major issue with sDMD is that the algorithm is sensitive to noisy data. Hemati et al. [37] found there is a systematic bias in the eigenvalue spectrum with the addition of sensor noise in the testing data, and the bias will not be eliminated when increasing the sample size. The influence of small noise on DMD has been carefully explored and characterized in [38]. There are two dominant solutions to denoise or debias the noise’ effect. In the first method, by Dawson et al. [39], the systematic bias was removed by averaging the data in forward and backward time, and it is known as fbDMD. In the second method, by Hemati et al. [37], the DMD algorithm is solved by applying total least squares algorithm instead of the standard least-squares algorithm, and it is known as tlsDMD. Dawson et al. [39] carried out a comprehensive discussion and accessible comparison with the above two methods. Both fbDMD and tlsDMD do not require knowledge of noise characteristics and come to very similar results; however, fbDMD might be an ill-conditioned operation for some data. In this paper, tlsDMD algorithm for mechanical fault signal diagnosis is proposed, solving the noise sensitivity problem of sDMD.More broadly, mechanical vibration signals are often comprised of multiscale temporal and/or spatial features of interest. A significant challenge is making effective and efficient connections between microscale and macroscale effects, such as extracting the magnitude of the complex system’s dynamic modes by orders in space and/or time. STFT and WT are structured to perform such multiresolution analyses. STFT defines a function acting as a time filter for localizing the processing signal over a specific window in time [4]. Frequencies are obtained at each instance of time with the method of integration over the parameter τ, which slides the time filtering window down the entire time domain of the signal. STFT proposes two main principles for time-frequency analysis: translation and scaling of a short-time window. WT firstly splits a signal up into a collection of a series of smaller signals by translating the mother wavelet with translation parameter b over the entire signal, and secondly, WT processes the signals at different frequency resolutions by scaling the wavelet window with the dilation parameter a [6]. As previously mentioned, DMD can decompose a dynamic model into a series of single-frequency nonorthogonal modes, which are regarded as the slow and fast modes and also regarded them as foregrounds and backgrounds. We are inspired to classify the DMD modes in a recursive modality, as the STFT and WT perform the multiresolution analyses. Multiresolution DMD (mrDMD) is developed in this paper to separate the mechanical vibration signal into grouped spatiotemporal DMD modes, we say, multiscale spatiotemporal modes. It is worth noting that in the process of mrDMD algorithm, tlsDMD algorithm is appropriate for eliminating the noise sensitivity.Theoretically, the fault characteristics of a given mechanical vibration fault signal are concealed in one or more mrDMD modes. Careful observations are needed to identify the failure type according to the frequency domain of each mrDMD structure, which brings about great inconvenience. For different vibration signals, the entropy values are also different due to different time complexity [40]. Therefore, the entropy values of mrDMD modes can be used to distinguish different fault types. Here, we take approximate entropy (ApEn) [41] as a vibration signal fault pattern recognizer and setting the entropy values as the input vector of the back propagation neural network (BPNN) [42] to perform fault classification. In summary, a fault diagnosis scheme based on mr-tlsDMD and ApEn is presented for one-dimensional mechanical vibration signal’ fault diagnosis.The layout of this paper is as follows: Section2 introduces the basic theory, including the algorithm of tlsDMD, mrDMD, and ApEn. Section 3 summarizes the current fault diagnosis methods and results based on the public vibration signal datasets and puts forward the technical route of the proposed method. Analysis results of three different experimental mechanical vibration fault signals are described in Section 4. Conclusions are summarized in Section 5.


## 2. Methodology

In this part, three algorithms adopted in this paper are outlined in detail for mechanical vibration signal fault diagnosis. In terms of the algorithmic solutions, tlsDMD is used to eliminate the effect of sensor noise, mrDMD is outlined to define a hierarchical application of mechanical fault signal, and ApEn is employed to calculate the entropy values of the mr-tlsDMD modes which are the inputs of BPNN.


### 2.1. tlsDMD Algorithm for Mechanical Vibration Signal Processing

DMD stems from the fact that it is an equation-free, data-driven method capable of providing potential dynamic characteristics of a complex high-dimensional system [23]. Essentially, it is a kind of order reduction methods based on the theory of singular value decomposition (SVD) and mode decomposition. The detailed algorithm process of sDMD is described in [23, 24]. Currently, DMD algorithm and multiple variations of DMD algorithm have been successfully applied in fluid dynamics, enhancing the capabilities for predication, state estimation, and control of complex system [26–29, 31, 37–39, 43]. Among that, tlsDMD showed its better ability to extract real and accurate dynamic features from noise-corrupted data. As we introduced in the introduction, Dawson et al. [39] proposed the tlsDMD algorithm and proved that it is insensitive to noise in the fluid domain. Here, we generalize the tlsDMD algorithm framework for one-dimensional mechanical vibration signal processing to eliminate the effect of sensor noise.Suppose that a one-dimensional time series ofN sample points is collected from a mechanical system at equal interval intervals. S=x1,x2,…,xi,…,xN,xi∈R, while the parameter of the given time interval Δt is defined as Δt=ti+1−ti. Therefore, the sampling frequency of the mechanical system is typically calculated as Fs=1/Δt. We transformed the 1×N observed signal S into m×n shift-stack Hankel matrix S′ by(1)S′=x1x2…xnx2x3…xn+1⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅xixi+1…xn+i−1⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅xmxm+1…xm+n−1=S1S2…Sn,S′∈Rm×n.It is possible to divide the matrixS′ into two m×n−1 data matrices X and Y.(2)X=S1S2…Sn−1,Y=S2S3…Sn,X,Y∈Rm×n−1.DMD is algorithmically a regression by assuming the so-called Koopman operatorA, mapping the current data X to the subsequent data Y, defined as Y=AX. Though the Koopman operator A is a linear, infinite-dimensional operator, it can represent the action of a nonlinear dynamical system in the Hilbert space [44]. Apparently, the evolution of the sequences X and Y is determined by the eigenvalues of A.Vibration signals of the mechanical system collected by sensors are inevitably mixed with noises [11]. In other words, the above matrices X and Y both contain noise components. The algorithm of sDMD takes the influence of noises in the matrix X into consideration, assuming that the first r-order eigenvalues of X represent the actual vibration signals and the subsequent m−r order eigenvalues as noises. Here, we consider that m≤n−1, and X is a full rank matrix, which is consistent with the actual mechanical vibration noisy signal. That is, sDMD identifies dynamics with(3)A:Y=AX¯+EX,argminAEXF,where X¯ and EX, respectively, represent actual vibration signals and noise components of the matrix X. ⋅F represents the Frobenius norm.While tlsDMD performs a single algorithm that finds total least squares solutions for the noise error on bothX and Y, we seek(4)A:Y¯+EY=AX¯+EX,argminAEXEYF.Similarly,Y¯ and EY represent actual vibration signals and noise components of the matrix Y, respectively.In what follows, we outline the core steps of tlsDMD algorithm for signal processing:(1)
Firstly, construct an augmented matrixZ=XY and compute the SVD on Z=U∑V∗, where Z∈R2m×n−1.U and V are orthonormal and are called the left and right singular vectors, respectively. U∗U=I,V∗V=I. The symbol ∗ denotes the complex conjugate transpose. ∑∈Rp×p contains a number of nonzero singular values σ1,σ2,…,σp by descending sequence in its diagonal, p=min2m,n−1.It is well known that when the noisy matrixZ is close to the square matrix, the maximum characteristic frequency components can be obtained with SVD. For the sake of obtaining the maximum spatial and temporal complexity of the original noisy signal, the product of 2m and n−1 should be as large as possible. According to the principle of inequality, the product achieves the maximum when 2m and n−1 are equal or close to each other. Then, the dimension m of Equation (1) is defined as follows:(5)m=N4,ifN=2z,N−14,ifN=2z−1,where z∈R is the positive integer sequence that z=1,2,3,….We choose a given number of truncated ranksr and retain the first r columns of V. Similarity matrices X˜ and Y˜ can be obtained by projecting Vr onto proper orthogonal decomposition modes. Then, the solution of the total least squares problem comes to two “corrected” matrices to account for noise.(6)X˜=XVrVr∗,Y˜=YVrVr∗,where X˜,Y˜∈Rr×r. Vr∗ is the complex conjugate transpose of the matrix Vr.Thus, tlsDMD identifies the mechanical signal system’s dynamics with(7)A:Y=AX,argminAEXEYF,where A˜ is the similarity matrix of A, defined as A˜=UrAUr∗, A˜∈Rr×r.The eigenvalues and eigenvectors ofA are then represented by those of A˜ as they process the same dynamical features.(2)
The matrixA˜ of Equation (4) may be obtained by the pseudoinverse of X after applying SVD on it.(8)A=YX−1=YVr∑r−1Ur∗.(3)
Eigenvalue decomposition is applied on the similarity matrixA˜:(9)A˜=WΛW−1,where W=ω1,ω2,…,ωr∈Rr×r is the eigenvector, and Λ=diagλ1,λ2,…,λr∈Rr×r is the corresponding diagonal eigenvalue matrix.(4)
Compute the approximate solution of the reconstruction matrixSDMD, representing the dynamic information of the original mechanical signal.(10)SDMD=∑i=1rϕiexpϖiΔtbi=ΦexpΩtb,SDMD∈Rm×n−1,where Φ is a matrix consisted of DMD modes ϕi, ϕi=YVi∑i−1Wi. Ω is a diagonal matrix whose entries are the dynamic system’s eigenvalues ϖi, ϖi=lnλi/Δt. t is a time series t=Δt,2Δt,…,N∗Δt. b is a vector containing the initial amplitude of each mode bi, b=ΦΓX, and Γ denotes the Moore–Penrose pseudoinverse.Similar to sDMD algorithm, tlsDMD is based on SVD. Determining how many singular values to select (truncated rankr) is of the most choice. In fact, the selection depends on a number of factors including the quality of the signals collected by the sensor, the signal-to-noise ratio (SNR), and the energy level excited by the fault characteristics. In the early stage of mechanical signal failure, the excited energy is low, leading to the fault characteristics that are not obvious. Moreover, it is usually unavoidable, and the experimental mechanical signal may be corrupted with noise, making it unclear to determine where to truncate. In the past, many researchers have done a lot of work in the field of truncation criteria. The proposed methods include hard threshold [45], soft threshold [46], and the nonconvex optimization algorithm [47]. The above methods either need to know the SNR in advance or add additional parameters. In this paper, we take the mechanical early failure signal into consideration, and choose the truncated rank r by retaining 99% of the singular value energy. In other words, the mechanical signal is slightly denoised to ensure that the fault signals are in the reconstructed matrix SDMD. The implementation process of tlsDMD for one-dimensional simulation signal is shown in Figure 1.Figure 1
tlsDMD implementation process for one-dimensional simulation signal.


### 2.2. mrDMD Framework

The reconstructed matrixSDMD contains all the dynamic characteristic information of the original system, including the increased/decreased information. In fluid and image processing, the dynamic information is usually divided into two separations: foreground and background [29]. The background refers to the invariant inherent characteristic, and we say, it represents slow modes in a mechanical signal. The foreground refers to the additional dynamic characteristic, which represents the fast modes in a mechanical signal. It is unclear that the fault signal accounts for the energy level for a given system. That is to say, in the series of single-frequency modes Φ, the fault characteristic frequency may exist in the slow modes component or/and in the fast modes component, making the fault extraction more difficult. In this paper, multiscale method is adopted to decompose the reconstructed signal into multiple scales with different levels and different time-frequency resolutions.In Equation (10), the reconstruction matrix of SDMD is obtained by a superposition of ϕiexpϖiΔtbi. Every single ϖi corresponds to a DMD mode ϕi. It becomes apparent that any period of the signal, which does not change or changes very slowly with time, must have an associated Fourier mode with the frequency ϖi⟶0. We regard the slowest-rank components as background modes by setting a threshold ξ1 and producing a representation of the reconstructed matrix SDMD into the background modes reconstruction matrix S11DMDϖi≤ξ1 and the remaining modes reconstruction matrix S12DMDϖi>ξ1. Equation (10) can be expressed as(11)SDMD=S1DMD=S11DMDϖi≤ξ1+S12DMDϖi>ξ1=∑ϖi≤ξ1ϕiexpϖiΔtbi+∑ϖi>ξ1ϕiexpϖiΔtbi.Mathematically, we divide the solution ofSDMD into slow modes S11DMDϖi≤ξ1 and fast modes S12DMDϖi>ξ1, where the symbol 1 at the bottom of the capitalized letter S represents the first level of reconstruction, the symbol 11 represents first reconstruction level’s slow modes, and the symbol 12 represents first reconstruction level’s fast modes. The tlsDMD algorithm outlined in the previous section can now be performed once again on the data matrix sm/2′. Thus the matrix sm/2′ can be separated again into slow modes and fast modes by setting another threshold ξ2.(12)S2DMD=S21DMDϖi≤ξ2+S22DMDϖi>ξ2=∑ϖi≤ξ2ϕiexpϖiΔtbi+∑ϖi>ξ2ϕiexpϖiΔtbi.This recursive process continues repeatedly until a desired termination. The dynamic characteristics of the original signal are represented by multiple levels of spatiotemporal modes. Accordingly, a hierarchical DMD has been introduced to split the original reconstructed signal matrixSDMD with several levels. After applying arbitrary L-th level with mrDMD, the matrix SDMD can then be reconstructed as follows:(13)SDMD=XL1DMD+XL2DMD+⋯+XL2L−1DMD.Figure2 illustrates the process in terms of the mrDMD solution. Multiple scale signals jointly represent the dynamic characteristic information of the original signal.Figure 2
Illustration of mrDMD solution. Here, the slow modes and fast modes are divided by splitting the number ofϖi at its median.


### 2.3. Basic Theory of ApEn

So far,2L−1 modes with different spatiotemporal resolutions are obtained by performing tlsDMD in each decomposition level. These multiresolution modes are considered similar to the IMFs of EMD. All the modes together represent the dynamic characteristics of the original one-dimensional signal, and each scale mode represents part of the inherent characteristics. Fault features may be concealed in one and/or more modes, making it inconvenient to identify fault types efficiently. For different mechanical vibration signals, due to different time complexity, their entropy values are also different. So, different signal types can be distinguished according to entropy values. In this paper, ApEn values of the 2L−1 modes are used as a fault recognizer to identify the fault types.ApEn, proposed by Pincus [41] in 1991, is an algorithm to measure the complexity and regularity of time series without coarsening. Owing to its advantages of strong anti-interference ability and short required data (least 1000 data), ApEn has been widely used in either mechanical signal or medical signal processing [48, 49]. One disadvantage of ApEn is that it measures the complexity of time series on a single scale. Considering that we have conducted multiscale decomposition of the original one-dimensional mechanical signal, multiscale entropy is not performed in this work for simplifying the overall algorithm. Here, we make full use of ApEn’s short requirement for time series and briefly describe its algorithm process [41, 50].For a given time seriesX,kk=1,2,…,N, firstly determine the pattern dimension τ and extract the elements in Xk sequentially, forming a set of vectors.(14)X=ixi,xi+1,…,xi+τ−1,i=1,2,…,N−τ+1.And determine the distance between any two vectorsXi and Xj by calculating the absolute value of the maximum difference between their corresponding elements.(15)dX,iXj=maxxi+k−xj+k,k=0,1,…,τ−1;i,j=1,2,…,N−τ+1;i≠j.Secondly, count the number of statistical distance which is smaller than a certain thresholdγ, and calculate the ratio with the number of vector Xi.(16)Ciτγ=ηN−τ+1,where η is the number of dX,iXj which is smaller than the threshold γ. And γ, the similarity tolerance, needs to be defined in advance, the same as the parameter τ.Subsequently, define the self-correlation degreeΦγτ of the sequence Xi.(17)Φγτ=1N−τ+1∑i=1N−τ+1lnCiτγ.Then, repeat the above steps by changing the pattern dimensionτ into τ+1, getting another self-correlation degree Φγτ+1.Finally, ApEn is defined for the original given time series.(18)ApEnτ,γ=Φγτ−Φγτ+1.The value of ApEn is related to the data lengthN, pattern dimension τ, and similarity tolerance γ. When τ=1or2 and γ=0.1∼0.25SD, (SD is the standard deviation of the original data sequence), the ApEn calculated has a reasonable statistical characteristic [41, 50]. In this paper, we choose τ=2 to obtain more information during the joint reconstruction. For the similarity tolerance, too much statistical information will be lost and too little will result in unsatisfactory statistical effect [41]; therefore, we set γ=0.2SD.


## 2.1. tlsDMD Algorithm for Mechanical Vibration Signal Processing

DMD stems from the fact that it is an equation-free, data-driven method capable of providing potential dynamic characteristics of a complex high-dimensional system [23]. Essentially, it is a kind of order reduction methods based on the theory of singular value decomposition (SVD) and mode decomposition. The detailed algorithm process of sDMD is described in [23, 24]. Currently, DMD algorithm and multiple variations of DMD algorithm have been successfully applied in fluid dynamics, enhancing the capabilities for predication, state estimation, and control of complex system [26–29, 31, 37–39, 43]. Among that, tlsDMD showed its better ability to extract real and accurate dynamic features from noise-corrupted data. As we introduced in the introduction, Dawson et al. [39] proposed the tlsDMD algorithm and proved that it is insensitive to noise in the fluid domain. Here, we generalize the tlsDMD algorithm framework for one-dimensional mechanical vibration signal processing to eliminate the effect of sensor noise.Suppose that a one-dimensional time series ofN sample points is collected from a mechanical system at equal interval intervals. S=x1,x2,…,xi,…,xN,xi∈R, while the parameter of the given time interval Δt is defined as Δt=ti+1−ti. Therefore, the sampling frequency of the mechanical system is typically calculated as Fs=1/Δt. We transformed the 1×N observed signal S into m×n shift-stack Hankel matrix S′ by(1)S′=x1x2…xnx2x3…xn+1⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅xixi+1…xn+i−1⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅xmxm+1…xm+n−1=S1S2…Sn,S′∈Rm×n.It is possible to divide the matrixS′ into two m×n−1 data matrices X and Y.(2)X=S1S2…Sn−1,Y=S2S3…Sn,X,Y∈Rm×n−1.DMD is algorithmically a regression by assuming the so-called Koopman operatorA, mapping the current data X to the subsequent data Y, defined as Y=AX. Though the Koopman operator A is a linear, infinite-dimensional operator, it can represent the action of a nonlinear dynamical system in the Hilbert space [44]. Apparently, the evolution of the sequences X and Y is determined by the eigenvalues of A.Vibration signals of the mechanical system collected by sensors are inevitably mixed with noises [11]. In other words, the above matrices X and Y both contain noise components. The algorithm of sDMD takes the influence of noises in the matrix X into consideration, assuming that the first r-order eigenvalues of X represent the actual vibration signals and the subsequent m−r order eigenvalues as noises. Here, we consider that m≤n−1, and X is a full rank matrix, which is consistent with the actual mechanical vibration noisy signal. That is, sDMD identifies dynamics with(3)A:Y=AX¯+EX,argminAEXF,where X¯ and EX, respectively, represent actual vibration signals and noise components of the matrix X. ⋅F represents the Frobenius norm.While tlsDMD performs a single algorithm that finds total least squares solutions for the noise error on bothX and Y, we seek(4)A:Y¯+EY=AX¯+EX,argminAEXEYF.Similarly,Y¯ and EY represent actual vibration signals and noise components of the matrix Y, respectively.In what follows, we outline the core steps of tlsDMD algorithm for signal processing:(1)
Firstly, construct an augmented matrixZ=XY and compute the SVD on Z=U∑V∗, where Z∈R2m×n−1.U and V are orthonormal and are called the left and right singular vectors, respectively. U∗U=I,V∗V=I. The symbol ∗ denotes the complex conjugate transpose. ∑∈Rp×p contains a number of nonzero singular values σ1,σ2,…,σp by descending sequence in its diagonal, p=min2m,n−1.It is well known that when the noisy matrixZ is close to the square matrix, the maximum characteristic frequency components can be obtained with SVD. For the sake of obtaining the maximum spatial and temporal complexity of the original noisy signal, the product of 2m and n−1 should be as large as possible. According to the principle of inequality, the product achieves the maximum when 2m and n−1 are equal or close to each other. Then, the dimension m of Equation (1) is defined as follows:(5)m=N4,ifN=2z,N−14,ifN=2z−1,where z∈R is the positive integer sequence that z=1,2,3,….We choose a given number of truncated ranksr and retain the first r columns of V. Similarity matrices X˜ and Y˜ can be obtained by projecting Vr onto proper orthogonal decomposition modes. Then, the solution of the total least squares problem comes to two “corrected” matrices to account for noise.(6)X˜=XVrVr∗,Y˜=YVrVr∗,where X˜,Y˜∈Rr×r. Vr∗ is the complex conjugate transpose of the matrix Vr.Thus, tlsDMD identifies the mechanical signal system’s dynamics with(7)A:Y=AX,argminAEXEYF,where A˜ is the similarity matrix of A, defined as A˜=UrAUr∗, A˜∈Rr×r.The eigenvalues and eigenvectors ofA are then represented by those of A˜ as they process the same dynamical features.(2)
The matrixA˜ of Equation (4) may be obtained by the pseudoinverse of X after applying SVD on it.(8)A=YX−1=YVr∑r−1Ur∗.(3)
Eigenvalue decomposition is applied on the similarity matrixA˜:(9)A˜=WΛW−1,where W=ω1,ω2,…,ωr∈Rr×r is the eigenvector, and Λ=diagλ1,λ2,…,λr∈Rr×r is the corresponding diagonal eigenvalue matrix.(4)
Compute the approximate solution of the reconstruction matrixSDMD, representing the dynamic information of the original mechanical signal.(10)SDMD=∑i=1rϕiexpϖiΔtbi=ΦexpΩtb,SDMD∈Rm×n−1,where Φ is a matrix consisted of DMD modes ϕi, ϕi=YVi∑i−1Wi. Ω is a diagonal matrix whose entries are the dynamic system’s eigenvalues ϖi, ϖi=lnλi/Δt. t is a time series t=Δt,2Δt,…,N∗Δt. b is a vector containing the initial amplitude of each mode bi, b=ΦΓX, and Γ denotes the Moore–Penrose pseudoinverse.Similar to sDMD algorithm, tlsDMD is based on SVD. Determining how many singular values to select (truncated rankr) is of the most choice. In fact, the selection depends on a number of factors including the quality of the signals collected by the sensor, the signal-to-noise ratio (SNR), and the energy level excited by the fault characteristics. In the early stage of mechanical signal failure, the excited energy is low, leading to the fault characteristics that are not obvious. Moreover, it is usually unavoidable, and the experimental mechanical signal may be corrupted with noise, making it unclear to determine where to truncate. In the past, many researchers have done a lot of work in the field of truncation criteria. The proposed methods include hard threshold [45], soft threshold [46], and the nonconvex optimization algorithm [47]. The above methods either need to know the SNR in advance or add additional parameters. In this paper, we take the mechanical early failure signal into consideration, and choose the truncated rank r by retaining 99% of the singular value energy. In other words, the mechanical signal is slightly denoised to ensure that the fault signals are in the reconstructed matrix SDMD. The implementation process of tlsDMD for one-dimensional simulation signal is shown in Figure 1.Figure 1
tlsDMD implementation process for one-dimensional simulation signal.


## 2.2. mrDMD Framework

The reconstructed matrixSDMD contains all the dynamic characteristic information of the original system, including the increased/decreased information. In fluid and image processing, the dynamic information is usually divided into two separations: foreground and background [29]. The background refers to the invariant inherent characteristic, and we say, it represents slow modes in a mechanical signal. The foreground refers to the additional dynamic characteristic, which represents the fast modes in a mechanical signal. It is unclear that the fault signal accounts for the energy level for a given system. That is to say, in the series of single-frequency modes Φ, the fault characteristic frequency may exist in the slow modes component or/and in the fast modes component, making the fault extraction more difficult. In this paper, multiscale method is adopted to decompose the reconstructed signal into multiple scales with different levels and different time-frequency resolutions.In Equation (10), the reconstruction matrix of SDMD is obtained by a superposition of ϕiexpϖiΔtbi. Every single ϖi corresponds to a DMD mode ϕi. It becomes apparent that any period of the signal, which does not change or changes very slowly with time, must have an associated Fourier mode with the frequency ϖi⟶0. We regard the slowest-rank components as background modes by setting a threshold ξ1 and producing a representation of the reconstructed matrix SDMD into the background modes reconstruction matrix S11DMDϖi≤ξ1 and the remaining modes reconstruction matrix S12DMDϖi>ξ1. Equation (10) can be expressed as(11)SDMD=S1DMD=S11DMDϖi≤ξ1+S12DMDϖi>ξ1=∑ϖi≤ξ1ϕiexpϖiΔtbi+∑ϖi>ξ1ϕiexpϖiΔtbi.Mathematically, we divide the solution ofSDMD into slow modes S11DMDϖi≤ξ1 and fast modes S12DMDϖi>ξ1, where the symbol 1 at the bottom of the capitalized letter S represents the first level of reconstruction, the symbol 11 represents first reconstruction level’s slow modes, and the symbol 12 represents first reconstruction level’s fast modes. The tlsDMD algorithm outlined in the previous section can now be performed once again on the data matrix sm/2′. Thus the matrix sm/2′ can be separated again into slow modes and fast modes by setting another threshold ξ2.(12)S2DMD=S21DMDϖi≤ξ2+S22DMDϖi>ξ2=∑ϖi≤ξ2ϕiexpϖiΔtbi+∑ϖi>ξ2ϕiexpϖiΔtbi.This recursive process continues repeatedly until a desired termination. The dynamic characteristics of the original signal are represented by multiple levels of spatiotemporal modes. Accordingly, a hierarchical DMD has been introduced to split the original reconstructed signal matrixSDMD with several levels. After applying arbitrary L-th level with mrDMD, the matrix SDMD can then be reconstructed as follows:(13)SDMD=XL1DMD+XL2DMD+⋯+XL2L−1DMD.Figure2 illustrates the process in terms of the mrDMD solution. Multiple scale signals jointly represent the dynamic characteristic information of the original signal.Figure 2
Illustration of mrDMD solution. Here, the slow modes and fast modes are divided by splitting the number ofϖi at its median.


## 2.3. Basic Theory of ApEn

So far,2L−1 modes with different spatiotemporal resolutions are obtained by performing tlsDMD in each decomposition level. These multiresolution modes are considered similar to the IMFs of EMD. All the modes together represent the dynamic characteristics of the original one-dimensional signal, and each scale mode represents part of the inherent characteristics. Fault features may be concealed in one and/or more modes, making it inconvenient to identify fault types efficiently. For different mechanical vibration signals, due to different time complexity, their entropy values are also different. So, different signal types can be distinguished according to entropy values. In this paper, ApEn values of the 2L−1 modes are used as a fault recognizer to identify the fault types.ApEn, proposed by Pincus [41] in 1991, is an algorithm to measure the complexity and regularity of time series without coarsening. Owing to its advantages of strong anti-interference ability and short required data (least 1000 data), ApEn has been widely used in either mechanical signal or medical signal processing [48, 49]. One disadvantage of ApEn is that it measures the complexity of time series on a single scale. Considering that we have conducted multiscale decomposition of the original one-dimensional mechanical signal, multiscale entropy is not performed in this work for simplifying the overall algorithm. Here, we make full use of ApEn’s short requirement for time series and briefly describe its algorithm process [41, 50].For a given time seriesX,kk=1,2,…,N, firstly determine the pattern dimension τ and extract the elements in Xk sequentially, forming a set of vectors.(14)X=ixi,xi+1,…,xi+τ−1,i=1,2,…,N−τ+1.And determine the distance between any two vectorsXi and Xj by calculating the absolute value of the maximum difference between their corresponding elements.(15)dX,iXj=maxxi+k−xj+k,k=0,1,…,τ−1;i,j=1,2,…,N−τ+1;i≠j.Secondly, count the number of statistical distance which is smaller than a certain thresholdγ, and calculate the ratio with the number of vector Xi.(16)Ciτγ=ηN−τ+1,where η is the number of dX,iXj which is smaller than the threshold γ. And γ, the similarity tolerance, needs to be defined in advance, the same as the parameter τ.Subsequently, define the self-correlation degreeΦγτ of the sequence Xi.(17)Φγτ=1N−τ+1∑i=1N−τ+1lnCiτγ.Then, repeat the above steps by changing the pattern dimensionτ into τ+1, getting another self-correlation degree Φγτ+1.Finally, ApEn is defined for the original given time series.(18)ApEnτ,γ=Φγτ−Φγτ+1.The value of ApEn is related to the data lengthN, pattern dimension τ, and similarity tolerance γ. When τ=1or2 and γ=0.1∼0.25SD, (SD is the standard deviation of the original data sequence), the ApEn calculated has a reasonable statistical characteristic [41, 50]. In this paper, we choose τ=2 to obtain more information during the joint reconstruction. For the similarity tolerance, too much statistical information will be lost and too little will result in unsatisfactory statistical effect [41]; therefore, we set γ=0.2SD.


## 3. Fault Diagnosis Methods for Mechanical Vibration Signals

### 3.1. A Brief Overview of Previous Approaches to Fault Diagnosis

There are many types of vibration signals in the field of mechanical industry, especially in the vital equipment with an online monitoring system. Usually, transmission system, which is prone to failure, is the main monitoring object. Therefore, the vibration signals usually come from the rolling bearings or/and the gears. Fault diagnosis is an important research field with the development of intelligent information technology, and a large number of research reports based on bearing and gear vibration signals fault diagnosis are published every year. Here, we present a part of the review literatures [51–55] on fault identification of mechanical vibration signals reported in the last three years. Support vector machines (SVM) [56] and artificial neural networks (ANN) [57] have been reported as the most popular techniques used as fault diagnosis methods in the view of the review literatures.However, due to the lack of general and public mechanical databases, it is difficult to compare the fault diagnosis techniques proposed in the literatures [58]. Currently available open source datasets include the NASA Prognostics Data [59] and Case Western Reserve University (CWRU) bearing signals [60]. Accordingly, the fault diagnosis methods and the algorithm performance based on SVM and ANN for the two datasets are shown in Table 1.Table 1
Fault diagnosis methods and performance for NASA and CWRU signals based on SVM and ANN.
Ref.
Main features
Database, bearing states
Classifier
Accuracy (%), state samples
[61]
EMD energy entropy of the first eight IMFs
NASA, seven
ANN
93. Total 5394 records divided into five folds, training: Four folds, test: one fold
[62]
Time and frequency domain features
NASA, seven
Linear SVM and quadratic SVM
99.4 of linear SVM, 99.3 of quadratic SVM. Training: 80%, test: 20%
[63]
Time domain features
NASA, two
Classical SVM (CSVM), incremental SVM (ISVM)
Outer: 91.1 CSVM and 98.7 ISVM, inner: 92.0 CSVM and 94.5 ISVM. Training: 70, test: 30
[64]
FFT
NASA, two
1D convolutional neural networks (1D CNN)
97.1 of 1D CNN, 94.5 of FFT-SVM
[65]
LMD, sample entropy, and energy ratio
CWRU, four
SVM
100. Training: 60, test: 20
[66]
Time and frequency, EMD energy entropy
CWRU, even
Adaptive neuro-fuzzy inference system (ANFIS)
94.7 (average). Training: 140, test: 70
[67]
DWT (cluster-based feature extraction)
CWRU, ten
Probabilistic neural network
98.2 (maximum). Training: 168, test: 60
[68]
EMD sample entropy of the first ten IMFs
CWRU, six
Improved shuffled frog leaping algorithm (ISFLA)
100 for H, others 95.4 (maximum). Training: 140, test: 70
[69]
LMD-SVD
CWRU, ten
BPNN SVM, extreme learning machine (ELM)
97.7 (average) for BP, 98.8 (average) for SVM, 99.3 (average) for ELM. Test: 228
[70]
Continuous wavelet transform (CWT)
CWRU, ten
Convolutional neural network (CNN)
99.7 of CNN, 99.7 of CNN, 85.1 of BPNN training: 200, Test: 200
[71]
Transfer learning
CWRU, six
Neural networks
91.8 (total). Total: 4832, training: 1208
[72]
Time domain features
CWRU, ten
Hierarchical adaptive deep CNN (ADCNN)
99.7 (average). Training: 500, test: 500
[73]
SVD, the singular values
CWRU, four
ANN
95.1, training: 336, test: 144Intelligent fault diagnosis methods usually include three key steps: (1) data preprocessing; (2) feature extraction and selection; and (3) fault classification. Data preprocessing is typically processed to eliminate outliers [70], such as sensor noise. Traditional feature extraction methods are performed on the basis of time domain characteristics [63, 72], frequency domain characteristics [64, 67, 70], time and frequency domain characteristics [62, 66], and mode characteristics [10, 61, 65, 66, 73]. Feature selection technology, significantly reducing feature dimension and improving diagnosis efficiency while retaining important and representative features, plays a crucial role in the following step of fault classification [70]. In terms of the fault classification method, Kankar et al. [74] compared the performances of ANN and SVM for fault diagnosis of rotor bearing systems: the experimental results indicate that ANN has a higher classification accuracy than SVM in the cases considered in the study.


### 3.2. Fault Diagnosis Method in This Paper

As stated in Section2, a novel feature extraction and selection approach is defined with mr-tlsDMD. A hierarchical application of the one-dimensional mechanical fault signal is adopted to achieve multiple scale DMD modes in different spatial and temporal resolutions. In each layer of decomposition, tlsDMD is performed to eliminate the effect of sensor noise. Then, ApEn is implemented on each multiresolution mode, and a vector of 2L−1 number of entropy values can be obtained. This fault extraction and selection method is similar to the main features of the literatures [61, 65, 66, 68], in which the dynamic characteristics of the original mechanical system were obtained by EMD/LMD/VMD and entropy algorithm. Supposing that there is a bunch of one-dimensional fault signals with different types, we take a part of the ApEn vectors for which we knew the fault type (actually they can be obtained through industrial field accumulation) as the input vectors of the BPNN for training and take the remaining for testing to perform fault classification.In the implementation of our proposed algorithm, several parameters need to be defined in advance, including the truncated rankr in tlsDMD, the decomposition layers L and the thresholds of slow modes and fast modes in each layer of mrDMD, pattern dimension τ, and similarity tolerance γ of ApEn. For a given number of n periods of mechanical signals, the flowchart of the proposed algorithm is shown in Figure 3. The procedure of the fault diagnosis method proposed in this paper is as follows:(1)
Each period of the mechanical signal is decomposed by mrDMD in the decomposition layerl=1,2,…,L, and in each decomposition level, tlsDMD is firstly applied to eliminate noise interference with the method r rank truncating in SVD (as in Equation (8)). A total number of 2L−1 multiresolution modes characterizing the dynamic information of the original signal are obtained (as in Equation (13)).(2)
ApEn with the parametersτ and γ is performed on each multiresolution mode, resulting in an entropy value. The multiresolution modes of each signal correspond to 2L−1 entropy values, forming a fault feature vector. Thus, there are totally a number of n fault feature vectors for the n periods of mechanical signals.(3)
Based on the fault feature vectors, the condition of the rolling bearings is identified by BPNN. All the vectors are randomly divided into two parts, one part is taken as the training data of BPNN and the remaining as testing data to perform fault classification. A comparison with other previous fault diagnosis methods is made in terms of classification accuracy.Figure 3
A flowchart of the proposed algorithm.


## 3.1. A Brief Overview of Previous Approaches to Fault Diagnosis

There are many types of vibration signals in the field of mechanical industry, especially in the vital equipment with an online monitoring system. Usually, transmission system, which is prone to failure, is the main monitoring object. Therefore, the vibration signals usually come from the rolling bearings or/and the gears. Fault diagnosis is an important research field with the development of intelligent information technology, and a large number of research reports based on bearing and gear vibration signals fault diagnosis are published every year. Here, we present a part of the review literatures [51–55] on fault identification of mechanical vibration signals reported in the last three years. Support vector machines (SVM) [56] and artificial neural networks (ANN) [57] have been reported as the most popular techniques used as fault diagnosis methods in the view of the review literatures.However, due to the lack of general and public mechanical databases, it is difficult to compare the fault diagnosis techniques proposed in the literatures [58]. Currently available open source datasets include the NASA Prognostics Data [59] and Case Western Reserve University (CWRU) bearing signals [60]. Accordingly, the fault diagnosis methods and the algorithm performance based on SVM and ANN for the two datasets are shown in Table 1.Table 1
Fault diagnosis methods and performance for NASA and CWRU signals based on SVM and ANN.
Ref.
Main features
Database, bearing states
Classifier
Accuracy (%), state samples
[61]
EMD energy entropy of the first eight IMFs
NASA, seven
ANN
93. Total 5394 records divided into five folds, training: Four folds, test: one fold
[62]
Time and frequency domain features
NASA, seven
Linear SVM and quadratic SVM
99.4 of linear SVM, 99.3 of quadratic SVM. Training: 80%, test: 20%
[63]
Time domain features
NASA, two
Classical SVM (CSVM), incremental SVM (ISVM)
Outer: 91.1 CSVM and 98.7 ISVM, inner: 92.0 CSVM and 94.5 ISVM. Training: 70, test: 30
[64]
FFT
NASA, two
1D convolutional neural networks (1D CNN)
97.1 of 1D CNN, 94.5 of FFT-SVM
[65]
LMD, sample entropy, and energy ratio
CWRU, four
SVM
100. Training: 60, test: 20
[66]
Time and frequency, EMD energy entropy
CWRU, even
Adaptive neuro-fuzzy inference system (ANFIS)
94.7 (average). Training: 140, test: 70
[67]
DWT (cluster-based feature extraction)
CWRU, ten
Probabilistic neural network
98.2 (maximum). Training: 168, test: 60
[68]
EMD sample entropy of the first ten IMFs
CWRU, six
Improved shuffled frog leaping algorithm (ISFLA)
100 for H, others 95.4 (maximum). Training: 140, test: 70
[69]
LMD-SVD
CWRU, ten
BPNN SVM, extreme learning machine (ELM)
97.7 (average) for BP, 98.8 (average) for SVM, 99.3 (average) for ELM. Test: 228
[70]
Continuous wavelet transform (CWT)
CWRU, ten
Convolutional neural network (CNN)
99.7 of CNN, 99.7 of CNN, 85.1 of BPNN training: 200, Test: 200
[71]
Transfer learning
CWRU, six
Neural networks
91.8 (total). Total: 4832, training: 1208
[72]
Time domain features
CWRU, ten
Hierarchical adaptive deep CNN (ADCNN)
99.7 (average). Training: 500, test: 500
[73]
SVD, the singular values
CWRU, four
ANN
95.1, training: 336, test: 144Intelligent fault diagnosis methods usually include three key steps: (1) data preprocessing; (2) feature extraction and selection; and (3) fault classification. Data preprocessing is typically processed to eliminate outliers [70], such as sensor noise. Traditional feature extraction methods are performed on the basis of time domain characteristics [63, 72], frequency domain characteristics [64, 67, 70], time and frequency domain characteristics [62, 66], and mode characteristics [10, 61, 65, 66, 73]. Feature selection technology, significantly reducing feature dimension and improving diagnosis efficiency while retaining important and representative features, plays a crucial role in the following step of fault classification [70]. In terms of the fault classification method, Kankar et al. [74] compared the performances of ANN and SVM for fault diagnosis of rotor bearing systems: the experimental results indicate that ANN has a higher classification accuracy than SVM in the cases considered in the study.


## 3.2. Fault Diagnosis Method in This Paper

As stated in Section2, a novel feature extraction and selection approach is defined with mr-tlsDMD. A hierarchical application of the one-dimensional mechanical fault signal is adopted to achieve multiple scale DMD modes in different spatial and temporal resolutions. In each layer of decomposition, tlsDMD is performed to eliminate the effect of sensor noise. Then, ApEn is implemented on each multiresolution mode, and a vector of 2L−1 number of entropy values can be obtained. This fault extraction and selection method is similar to the main features of the literatures [61, 65, 66, 68], in which the dynamic characteristics of the original mechanical system were obtained by EMD/LMD/VMD and entropy algorithm. Supposing that there is a bunch of one-dimensional fault signals with different types, we take a part of the ApEn vectors for which we knew the fault type (actually they can be obtained through industrial field accumulation) as the input vectors of the BPNN for training and take the remaining for testing to perform fault classification.In the implementation of our proposed algorithm, several parameters need to be defined in advance, including the truncated rankr in tlsDMD, the decomposition layers L and the thresholds of slow modes and fast modes in each layer of mrDMD, pattern dimension τ, and similarity tolerance γ of ApEn. For a given number of n periods of mechanical signals, the flowchart of the proposed algorithm is shown in Figure 3. The procedure of the fault diagnosis method proposed in this paper is as follows:(1)
Each period of the mechanical signal is decomposed by mrDMD in the decomposition layerl=1,2,…,L, and in each decomposition level, tlsDMD is firstly applied to eliminate noise interference with the method r rank truncating in SVD (as in Equation (8)). A total number of 2L−1 multiresolution modes characterizing the dynamic information of the original signal are obtained (as in Equation (13)).(2)
ApEn with the parametersτ and γ is performed on each multiresolution mode, resulting in an entropy value. The multiresolution modes of each signal correspond to 2L−1 entropy values, forming a fault feature vector. Thus, there are totally a number of n fault feature vectors for the n periods of mechanical signals.(3)
Based on the fault feature vectors, the condition of the rolling bearings is identified by BPNN. All the vectors are randomly divided into two parts, one part is taken as the training data of BPNN and the remaining as testing data to perform fault classification. A comparison with other previous fault diagnosis methods is made in terms of classification accuracy.Figure 3
A flowchart of the proposed algorithm.


## 4. Fault Diagnosis of Experimental Mechanical Vibration Signals

Next, we implement three experiments to verify the effectiveness of the proposed method. The first two sets of rolling bearing’s signals come from CWRU and NASA, and the third set of gear’s signal comes from the data tested by authors’ research team.


### 4.1. Fault Diagnosis of Rolling Bearing Signals of CWRE

The rolling bearing fault simulation systems of CWRU mainly consists of a motor, a torque transducer, couplings, electronic control equipment, and other related devices. Please refer to [60] for details of the bearing experimental device and the test vibration signal. From the point that the signals adopted in this paper of view, four different forms (normal, inner race fault, ball fault, and outer race fault) of SKF bearings 6205-2RS were installed at the motor drive end. The fault bearings are machined with a single point of failure using the electro discharge technique with different sizes. A sensor, which is marked in blue triangle, was placed on the pedestal with a magnetic base. Digital data was collected by sampling with the frequency of 48 kHz in normal vibration and 12 kHz in fault vibration. Both the motor rotation speed was 1750 r/min. We used the test signals with motor power of 2 HP and bearing failure size of diameter 0.021″ to verify the effectiveness of the proposed algorithm. The normal bearing signal, inner race fault signal, ball fault signal, and outer race fault signal come from the database of 99, 211, 224, and 236, respectively. The sampling time of each group lasts for 10 seconds. A total of 120,000 data points were adopted in our experiment. Each time series is divided into 50 sections at intervals of 0.2 second, with the duration of 0.5 second. That is to say, the first section is taken from 0 to 0.5 seconds, the second from 0.2 to 0.7 seconds, and so on. Therefore, the total number of bearing state is 200 with four bearing states, and the number of data points per bearing state is 24,000.All bearing states were firstly decomposed by mr-tlsDMD. We define the parameter ofr by setting that the energy of the truncated rank orders reaches to 99% of the total energy in the algorithm of tlsDMD. The decomposition layers L is set to 4, and the threshold ξl,l=1,2,…,L of slow modes and fast modes are defined by splitting the number of ϖi at its median (the number of slow modes equals the number of fast modes). In other words, there are eight ApEn values for each bearing state when applying the algorithm of ApEn on the multiresolution modes later. Figure 4 shows the original one-dimensional signals with only one representation per bearing state listed. Here, we plot the signals with the first line of Hankel matrices, by using Equations (1) and (5). Figure 5 shows the FT spectrum of the signals in Figure 4. The specific parameters of the experimental bearing 6205-2RS are shown in Table 2, and the calculated frequency characteristics of the fault (refer to the formulas in [35]) are shown in Table 3. Multiscale decomposition modes of each signal are shown in Figures 6–9, and their FT spectrums are shown in Figures 10–13.Figure 4
The time domain diagram of CWRU signals, where (a)–(d) represent normal bearing signal, inner race fault signal, ball fault signal, and outer race fault signal, coming from 99, 211, 224, and 236 of CWRU, respectively. Each signal is listed with the sample number 6000.
(a)
(b)
(c)
(d)Figure 5
The Frequency domain diagram of CWRU signals, where (a)–(d) are corresponding to signals in Figure5, respectively. The characteristic frequency distinction of each fault form is quite obvious, especially in (a) and (b) for the normal signal and the inner race fault signal. It is easy to find out the AM frequency in (b). However, there are a lot of interference frequency components in the spectrum of the roller fault signal and outer race fault signal, as in (c) and (d). Moreover, the FM components (by the cage frequency) of the ball frequency cannot be detected.
(a)
(b)
(c)
(d)Table 2
Rolling element bearing parameters of 6205-2RS.
Inner diameter
Outer diameter
Ball number
Ball diameter
Contact angle
Pitch diameter
25 mm
52 mm
9
7.9 mm
0 rad.
46.4 mmTable 3
Calculated characteristic frequencies of 6205-2RS (Hz).
Rotation frequencyfr
Inner racefi
Outer racefo
Ball frequencyfb
Cage frequencyfc
29.53
159.92
105.87
139.20
11.69Figure 6
Multiresolution DMD modes of a normal roller bearing signal. Here, we select the first column ofXijDMD that represents the fourth layer slow or fast modes. (a)–(h) represent the modes of X41DMD, X42DMD, X43DMD, X44DMD, X45DMD, X46DMD, X47DMD, and X48DMD, respectively. The symbols in the Figures 8–10 share the same meaning.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 7
Multiresolution DMD modes of roller bearing’s inner race fault signal (CWRU).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 8
Multiresolution DMD modes of roller bearing’s ball fault signal (CWRU).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 9
Multiresolution DMD modes of roller bearing’s outer race fault signal (CWRU).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 10
The frequency domain diagram of multiresolution DMD modes of the normal roller bearing signal, where (a)–(h) are corresponding to signals in Figure6, respectively. The characteristic frequencies of the normal bearing signal are distributed in multiple modes, and these multiple modes jointly represent the features of the normal signals.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 11
The frequency domain diagram of multiresolution DMD modes of the inner race fault signal, where (a)–(h) are corresponding to signals in Figure7, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 12
The frequency domain diagram of multiresolution DMD modes of the roller fault signal, where (a)–(h) are corresponding to signals in Figure8, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 13
The frequency domain diagram of multiresolution DMD modes of the outer race fault signal, where (a)–(h) are corresponding to signals in Figure9, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Subsequently, ApEn algorithm (τ=2, γ=0.2SD, and SD is the standard deviation of data sequence) is applied on the multiresolution DMD modes, and a total of 200 vectors are obtained. Each vector contains eight elements corresponding to each bearing state. 30 vectors of each bearing state were randomly selected as the training data of BPNN and the remaining 20 vectors as the test data. Numbers of records per bearing state and ANN outputs are summarized in Table 4. The parameters of the classifier are set as follows: the input layer has 20 nodes, the hidden layer has 30 nodes, and the output layer has 4 nodes. The classification results of the training and testing signals are shown in Figure 14. The category labels 1, 2, 3, and 4 represent normal, inner race fault, ball fault, and outer race fault, respectively. Note that the classification accuracy is the ratio between the total numbers of correctly classified test samples to the total number of test samples. It is clear that the accuracy of the proposed algorithm is 100%, indicating that its performance is significant.(19)classification accuracy=Number of correctly classified samplesTotal number of classified in testing samples×100.Table 4
Number of records per bearing state and BPNN outputs.
Bearing state
Normal
Inner race fault
Ball fault
Outer race fault
Total records
50
50
50
50
Training records
30
30
30
30
Testing records
20
20
20
20
BPNN outputs
1000
0100
0010
0001Figure 14
The classification results of roller bearing signals.


### 4.2. Fault Diagnosis of Rolling Bearing Signals of NASA

In Section4.1, although we have successfully classified and identified four kinds of bearing faults, which proved the effectiveness of the proposed algorithm, the fault types can also be obtained from the spectrum diagram of the original signals. Next, we apply the algorithm on the early signals from NASA’s run-to failure bearing experiment to show the proposed algorithm’s effectiveness on weak fault signal.The used signals were generated by Intelligent Maintenance Systems of Cincinnnati, and the data set is available on the NASA website [59]. Four Rexnord ZA-2115 double row bearings were installed on a shaft, which was driven by an AC motor at the speed of 2,000 RPM via sever rub belts. The specific parameters of the experimental bearing are shown in Table 5, and the calculated frequency characteristics of the fault are shown in Table 6.Table 5
Rolling element bearing parameters of 6205-2RS.
Roller number
Roller diameter
Contact angle
Pitch diameter
16
8.407 mm
0.265 rad.
71.5 mmTable 6
Calculated characteristic frequencies of 6205-2RS (Hz).
Rotation frequencyfr
Inner racefi
Outer racefo
Roller frequencyfb
Cage frequencyfc
33.73
300.6
236.4
141.6
14.8High sensitivity quartz ICP accelerometers (PCB 353B33) were installed on the bearing housing, receiving the vibration signals, and NI DAQ Card 6062E was used to collect the signals. Each data set consists of individual files that are one second vibration signal snapshots. Each snapshot, recorded at every 10 minutes’ interval, consists of 20,480 points with the sampling frequency 20 kHz. Three tests were made. The first test was carried out for 35 days until inner race defect occurred in bearing 3 and roller element defect in bearing 4. The whole experiment was discontinuous, and the continuous running time of the bearings from normal to severe failure was about 15 days. Senanayaka et al. [62] calculated the root mean square (RMS) values of the whole snapshots. They concluded that the first 12 days of the bearings’ running duration were the normal service period, since then the bearing 3 (inner race fault) was subjected to a 2.5 days’ predegradation period, and severe failure of the bearing occurred in the last 6 hours. In our experiment, we use the first channel data of bearing 1 as normal bearing signals and the fifth channel data of bearing 3 as inner race fault signals. The normal bearing signals come from the snapshot 2003.10.31.08.59.46 (the snapshot after the bearings running for 2 days) and the following 100 snapshots of the “1st_test.” The inner fault signals come from the snapshot 2003.11.22.11.06.56 (the snapshot after the bearings running for 12.5 days) and the following 100 snapshots of the “1st_test.” The second test was carried out for 8 days until outer race failure occurred in bearing 1. Senanayaka et al. [62] concluded that the first 3.5 days of the bearings’ running duration were the normal service period, since then bearing 1 (outer race fault) was subjected to a 3.5 days’ pre-degradation period, and the severe failure of the bearing occurred in the last 3.5 hours. We use the first channel data of the snapshots from the snapshot of 2004.02.15.10.32.39 (the snapshot after the bearings running for 3 days) and the following 100 snapshots of the “2st_test” as outer fault signals. Hence, there are totally 300 samples representing three types of incipient failures (normal bearing, inner race fault, and outer race fault). Figure 15 shows the three kinds of signals with the first line of the Hankel matrix, derived from the snapshots 2003.10.31.08.59.46, 2003.11.22.11.06.56, and 2004.02.15.10.32.39, respectively. Figure 16 shows the FT spectrum of the signals in Figure 15. Multiscale decomposition modes of each signal are shown in Figures 17–19, and their FT spectrums are shown in Figures 20–22.Figure 15
The time domain diagram of the training signals, where (a), (b), and (c) represent normal bearing signal, inner race fault signal, and outer race fault signal come from the snapshots 2003.10.31.08.59.46, 2003.11.22.11.06.56, and 2004.02.15.10.32.39, respectively. Each signal is listed with the sample number 5120.
(a)
(b)
(c)Figure 16
The frequency domain diagram of the signals in Figure15.
(a)
(b)
(c)Figure 17
Multiresolution DMD modes of the normal roller bearing signal (NASA, snapshot 2003.10.31.08.59.46). Here, we select the first column ofXijDMD that represents the fourth layer slow or fast modes. (a)–(h) represent the modes of X41DMD, X42DMD, X43DMD, X44DMD, X45DMD, X46DMD, X47DMD, and X48DMD, respectively. The symbols in the Figures 18 and 19 share the same meaning.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 18
Multiresolution DMD modes of roller bearing’s inner race fault signal (NASA, snapshots 2003.11.22.11.06.56).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 19
Multiresolution DMD modes of roller bearing’s outer race fault signal (NASA, snapshot 2004.02.15.10.32.39).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 20
The Frequency domain diagram of multiresolution DMD modes of normal roller bearing signal, where (a)–(h) are corresponding to signals in Figure17, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 21
The Frequency domain diagram of multiresolution DMD modes of the inner race fault signal, where (a)–(h) are corresponding to signals in Figure18, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 22
The Frequency domain diagram of multiresolution DMD modes of the outer race fault signal, where (a)–(h) are corresponding to signals in Figure19, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)ApEn algorithm (τ=2, γ=0.2SD) is applied on the multi-resolution DMD modes, and a total of 300 vectors are obtained for the three bearing fault states. 60 vectors of each state are randomly selected as the inputs of BPNN, the remaining 40 vectors are selected for testing. Numbers of records per bearing state and ANN outputs are summarized in Table 7. The parameters of the classifier are set as the same as 4.1. The classification results of the training and testing signals are shown as in Figure 23. The category labels 1, 2, and 3 represent normal bearing, inner race fault, outer race fault, respectively. The average test accuracy of our algorithm is 98.7%.Table 7
Number of records per bearing state and BPNN outputs.
Bearing state
Normal
Inner race fault
Outer race fault
Total records
100
100
100
Training records
60
60
60
Testing records
40
40
40
BPNN outputs
100
010
001Figure 23
The classification results of roller bearing signals.


### 4.3. Fault Diagnosis of Gear Vibration Signals

Figure24 demonstrates the gear signal testing device consisting of an AC motor, couplings, dynamometers, a magnetic powder brake, and a gearbox. A pair of meshing gears is internally fitted in the gearbox. The numbers of teeth for the driving and driven wheel are 20 and 37, and the gear modulus is equal to 3. The motor speed is maintained at 1154 r/min. A sensor is displayed on the upside of the gearbox cover. Vibration signals of four different faults (normal gear, broken teeth fault, pitch errors fault, and wear fault) are collected with the sampling frequency 2000 Hz. 19 group data for each meshing pair were gathered under the brake load of 10 N/m, and each data sample has 12,000 data points.Figure 24
Gear signal testing device. (a) Physical object of the device and (b) structure diagram of test device: (1) AC motor; (2) coupling; (3) gearbox; (4) dynamometers; and (5) magnetic powder brake.
(a)
(b)Figure25 shows the signals with the first line of Hankel matrices, with only one representation per fault form listed. Figure 26 shows the FT spectrum of the signals in Figure 25. Multiscale decomposition modes of them are presented from Figures 27–30, and their FT spectrums are shown in Figures 31–34. All parameters in the computational process are selected the same as the roller bearing’s signal experiment, except for the sample numbers.Figure 25
The time domain diagram of the collected signals, where (a)–(d) represent normal gear signal, broken teeth fault signal, pitch errors signal, and wear gear signal, respectively. Each signal is listed with 3000 points.
(a)
(b)
(c)
(d)Figure 26
The frequency domain diagram of the signals in Figure25, where FGM is the gear mesh frequency.
(a)
(b)
(c)
(d)Figure 27
Multiresolution DMD modes of the normal gear signal. We select the first column ofXijDMD that represents the fourth layer slow or fast modes. (a)–(h) represent the modes of X41DMD, X42DMD, X43DMD, X44DMD, X45DMD, X46DMD, X47DMD, and X48DMD, respectively. The symbols in Figures 28–30 share the same meanings.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 28
Multiresolution DMD modes of the broken teeth fault signal.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 29
Multiresolution DMD modes of the pitch error gear fault signal.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 30
Multiresolution DMD modes of the wear gear fault signal.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 31
The frequency domain diagram of multiresolution DMD modes of the normal gear signal, where (a)–(h) are corresponding to signals in Figure27, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 32
The frequency domain diagram of multiresolution DMD modes of the broken teeth fault signal, where (a)–(h) are corresponding to signals of Figure28, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 33
The Frequency domain diagram of multiresolution DMD modes of the pitch errors signal, where (a)–(h) are corresponding to signals in Figure29, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 34
The frequency domain diagram of multiresolution DMD modes of and wear gear signal, where (a)–(h) are corresponding to signals in Figure30, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)As in Figure26, although the spectrums of the four gears’ fault signals are visually different, it is difficult to identify the characteristic frequencies of different faults with the perspective of amplitude and AM/FM components. In Figures 31–34, the gear meshing frequencies are distributed in multiple multiresolution modes, and similarly, the series of AM/FM components cannot be distinguished. It is hard to judge the gear fault type either by the FT spectrums or the multiresolution DMD modes. Therefore, we apply ApEn on the multiresolution modes to calculate their complexity, and take the ApEn vector as the characteristic information of each signal.ApEn algorithm (τ=2, γ=0.2SD) is applied on the multiresolution DMD modes, and a total of 76 vectors are obtained for the four gear states. 9 vectors of each gear state type are randomly selected as the inputs of BPNN, the remaining are selected for testing. Numbers of records per gear state and ANN outputs are summarized in Table 8. The parameters of the classifier are set the same as the parameters of the roller bearing’s experiment. The classification results of the training and testing signals are shown as in Figure 35. The category labels 1, 2, 3, and 4 represent normal gear, broken teeth fault, pitch errors fault, and wear fault, respectively. The average test accuracy of our algorithm is 97.5%.Table 8
Number of records per gear state and BPNN outputs.
Bearing state
Normal
Broken teeth
Pitch errors
Wear gear
Total records
19
19
19
19
Training records
9
9
9
9
Testing records
10
10
10
10
BPNN outputs
1000
0100
0010
0001Figure 35
The classification results of gear signals.


### 4.4. Result Discussion

Table1 summarizes the fault diagnosis methods and the algorithm performance based on SVM and ANN using data from the CWRU and NASA bearing databases. Since mr-tlsDMD is essentially an algorithm based on SVD and mode decomposition, we only compare the proposed method with the literatures that adopted mode decomposition algorithm in the step of feature extraction and selection. The comparisons of the CWRU database are shown in Table 9. In Ref. [65], LMD and sample entropy were used to extract features and SVM was applied to classify the four bearing states, and finally 100% accuracy was achieved. Tian et al. [69] integrated LMD, SVD, and BPNN/SVM/ELM to distinguish ten bearing states and got the testing accuracy with 97.7% (BPNN), 98.8% (SVM), and 99.3 (ELM), respectively. Zhao et al. [68] used EMD sample entropy of the first eight IMFs to get features and ISFLA to classify six bearing states and achieved 94.7% accuracy. Lei et al. [66] proposed a method combining time-frequency domain features and EMD energy entropy to extract features and ANFIS to diagnose faults, obtaining a classification accuracy of 95.4% for faulty bearing and 100% for healthy bearing. The comparisons of the NASA database are shown in Table 10. Jaouher et al. [53] introduced a feature extraction method based on EMD energy entropy. The first eight IMFs’ energy entropies were used to extract features, ANN was applied to classify seven bearing states, and finally 93% accuracy was achieved.Table 9
Classification comparisons of CWRU database-adopted mode decomposition algorithm.
Ref.
Main features
Bearing states
Training samples
Testing samples
Test accuracy (%)
[65]
LMD, sample entropy, energy ratio + SVM
Four
60
20
100
[69]
LMD-SVD + BPNN/SVM/ELM
Ten
172
228
97.7, 98.8, 99.3
[68]
EMD sample entropy of first ten IMFs + ISFLA
Six
140
70
100 for health, others 95.4
[66]
Time and frequency, EMD energy entropy + ANFIS
Even
140
70
94.7
Proposed
mr-tlsDMD, ApEn + BPNN
Four
120
80
100Table 10
Classification comparisons of NASA database-adopted mode decomposition algorithm.
Ref.
Main features
Bearing states
Training samples
Testing samples
Test accuracy (%)
[53]
Eight IMFs of EMD, energy entropy + SVM
Seven
2728
2666
93
Proposed
mr-tlsDMD, ApEn + BPNN
Three
180
120
98.7Ge et al. [75] present a joint fault diagnosis scheme via tensor nuclear norm canonical polyadic decomposition (TNNCPD) and multiscale permutation entropy (MSPE) to classify four gear states. BPNN isused to perform fault classification. The results illustrate that their proposed scheme can accurately recognize different gear working conditions, while the performance of WT is worse. It is worth noting that Ge et al. [75] adopt the same set of signals in the experiment as we used in this paper. The comparisons between the performance in [75] and our proposed method are shown in Table 11. The test accuracy of our result is not as perfect as that in [75], nevertheless it is significantly better than WT.Table 11
Classification comparisons of the gear database (brake load: 10 N/m).
Ref.
Main features
Bearing states
Training samples
Testing samples
Accuracy (%)
Train
Test
[75]
TNNCPD, MSPE + BPNN
Four
5
14
100
100
[75]
WT, MSPE + BPNN
Four
5
14
90
76.8
Proposed
mr-tlsDMD, ApEn + BPNN
Four
10
9
100
97.5The computational efficiency is as important as the result’s accuracy for an emerging algorithm, especially when the algorithm is adopted on the online monitoring system. Therefore, the computational efficiency of the proposed algorithm is compared with that of other mode-based methods. Only the computation time of mode extraction procedure is concerned, without considering the computation time of ApEn and BPNN. The experiments were completed on a desktop computer with Matlab version 2014a, and the basic configuration of the computer was Intel® core™ i7-4770, CPU, 3.4 GHz, and memory capacity 32G. In the process of experiments, all comment language and result processing statements were removed. Table12 shows the computational efficiency of several mode-based algorithms for NASA snapshot 2003.10.31.08.59.46, which was used as a normal bearing signal in Section 4.2. It can be seen from the table that the calculation time of EMD and LMD is significantly shorter, that of VMD is 51.67 seconds, while that of mr-tlsDMD is 445.23 seconds. It is obvious that the proposed method is not superior in computational efficiency. In the conclusion part, we will analyze the reasons for the long calculation time.Table 12
Computational efficiency of several mode-based algorithms for the NASA snapshot (s.).
Method
EMD
LMD
VMD
mr-tlsDMD (L=4)
Time
0.87
1.43
51.67
445.23Three experiments are implemented to recognize the fault type of mechanical vibration signals with the proposed algorithm. The first two sets of experiment’s signals come from CWRU and NASA, which are public bearing database, and the third come from a gear test bench. Experiments, particularly on the early bearing fault signals of Section4.2 and complex gear fault signals of Section 4.3, show that the proposed fault diagnosis method has excellent classification performance. The proposed algorithm is robust whether dealing with bearing signals or gear signals; however, it is not as efficient as other modal decomposition algorithms.


## 4.1. Fault Diagnosis of Rolling Bearing Signals of CWRE

The rolling bearing fault simulation systems of CWRU mainly consists of a motor, a torque transducer, couplings, electronic control equipment, and other related devices. Please refer to [60] for details of the bearing experimental device and the test vibration signal. From the point that the signals adopted in this paper of view, four different forms (normal, inner race fault, ball fault, and outer race fault) of SKF bearings 6205-2RS were installed at the motor drive end. The fault bearings are machined with a single point of failure using the electro discharge technique with different sizes. A sensor, which is marked in blue triangle, was placed on the pedestal with a magnetic base. Digital data was collected by sampling with the frequency of 48 kHz in normal vibration and 12 kHz in fault vibration. Both the motor rotation speed was 1750 r/min. We used the test signals with motor power of 2 HP and bearing failure size of diameter 0.021″ to verify the effectiveness of the proposed algorithm. The normal bearing signal, inner race fault signal, ball fault signal, and outer race fault signal come from the database of 99, 211, 224, and 236, respectively. The sampling time of each group lasts for 10 seconds. A total of 120,000 data points were adopted in our experiment. Each time series is divided into 50 sections at intervals of 0.2 second, with the duration of 0.5 second. That is to say, the first section is taken from 0 to 0.5 seconds, the second from 0.2 to 0.7 seconds, and so on. Therefore, the total number of bearing state is 200 with four bearing states, and the number of data points per bearing state is 24,000.All bearing states were firstly decomposed by mr-tlsDMD. We define the parameter ofr by setting that the energy of the truncated rank orders reaches to 99% of the total energy in the algorithm of tlsDMD. The decomposition layers L is set to 4, and the threshold ξl,l=1,2,…,L of slow modes and fast modes are defined by splitting the number of ϖi at its median (the number of slow modes equals the number of fast modes). In other words, there are eight ApEn values for each bearing state when applying the algorithm of ApEn on the multiresolution modes later. Figure 4 shows the original one-dimensional signals with only one representation per bearing state listed. Here, we plot the signals with the first line of Hankel matrices, by using Equations (1) and (5). Figure 5 shows the FT spectrum of the signals in Figure 4. The specific parameters of the experimental bearing 6205-2RS are shown in Table 2, and the calculated frequency characteristics of the fault (refer to the formulas in [35]) are shown in Table 3. Multiscale decomposition modes of each signal are shown in Figures 6–9, and their FT spectrums are shown in Figures 10–13.Figure 4
The time domain diagram of CWRU signals, where (a)–(d) represent normal bearing signal, inner race fault signal, ball fault signal, and outer race fault signal, coming from 99, 211, 224, and 236 of CWRU, respectively. Each signal is listed with the sample number 6000.
(a)
(b)
(c)
(d)Figure 5
The Frequency domain diagram of CWRU signals, where (a)–(d) are corresponding to signals in Figure5, respectively. The characteristic frequency distinction of each fault form is quite obvious, especially in (a) and (b) for the normal signal and the inner race fault signal. It is easy to find out the AM frequency in (b). However, there are a lot of interference frequency components in the spectrum of the roller fault signal and outer race fault signal, as in (c) and (d). Moreover, the FM components (by the cage frequency) of the ball frequency cannot be detected.
(a)
(b)
(c)
(d)Table 2
Rolling element bearing parameters of 6205-2RS.
Inner diameter
Outer diameter
Ball number
Ball diameter
Contact angle
Pitch diameter
25 mm
52 mm
9
7.9 mm
0 rad.
46.4 mmTable 3
Calculated characteristic frequencies of 6205-2RS (Hz).
Rotation frequencyfr
Inner racefi
Outer racefo
Ball frequencyfb
Cage frequencyfc
29.53
159.92
105.87
139.20
11.69Figure 6
Multiresolution DMD modes of a normal roller bearing signal. Here, we select the first column ofXijDMD that represents the fourth layer slow or fast modes. (a)–(h) represent the modes of X41DMD, X42DMD, X43DMD, X44DMD, X45DMD, X46DMD, X47DMD, and X48DMD, respectively. The symbols in the Figures 8–10 share the same meaning.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 7
Multiresolution DMD modes of roller bearing’s inner race fault signal (CWRU).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 8
Multiresolution DMD modes of roller bearing’s ball fault signal (CWRU).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 9
Multiresolution DMD modes of roller bearing’s outer race fault signal (CWRU).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 10
The frequency domain diagram of multiresolution DMD modes of the normal roller bearing signal, where (a)–(h) are corresponding to signals in Figure6, respectively. The characteristic frequencies of the normal bearing signal are distributed in multiple modes, and these multiple modes jointly represent the features of the normal signals.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 11
The frequency domain diagram of multiresolution DMD modes of the inner race fault signal, where (a)–(h) are corresponding to signals in Figure7, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 12
The frequency domain diagram of multiresolution DMD modes of the roller fault signal, where (a)–(h) are corresponding to signals in Figure8, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 13
The frequency domain diagram of multiresolution DMD modes of the outer race fault signal, where (a)–(h) are corresponding to signals in Figure9, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Subsequently, ApEn algorithm (τ=2, γ=0.2SD, and SD is the standard deviation of data sequence) is applied on the multiresolution DMD modes, and a total of 200 vectors are obtained. Each vector contains eight elements corresponding to each bearing state. 30 vectors of each bearing state were randomly selected as the training data of BPNN and the remaining 20 vectors as the test data. Numbers of records per bearing state and ANN outputs are summarized in Table 4. The parameters of the classifier are set as follows: the input layer has 20 nodes, the hidden layer has 30 nodes, and the output layer has 4 nodes. The classification results of the training and testing signals are shown in Figure 14. The category labels 1, 2, 3, and 4 represent normal, inner race fault, ball fault, and outer race fault, respectively. Note that the classification accuracy is the ratio between the total numbers of correctly classified test samples to the total number of test samples. It is clear that the accuracy of the proposed algorithm is 100%, indicating that its performance is significant.(19)classification accuracy=Number of correctly classified samplesTotal number of classified in testing samples×100.Table 4
Number of records per bearing state and BPNN outputs.
Bearing state
Normal
Inner race fault
Ball fault
Outer race fault
Total records
50
50
50
50
Training records
30
30
30
30
Testing records
20
20
20
20
BPNN outputs
1000
0100
0010
0001Figure 14
The classification results of roller bearing signals.


## 4.2. Fault Diagnosis of Rolling Bearing Signals of NASA

In Section4.1, although we have successfully classified and identified four kinds of bearing faults, which proved the effectiveness of the proposed algorithm, the fault types can also be obtained from the spectrum diagram of the original signals. Next, we apply the algorithm on the early signals from NASA’s run-to failure bearing experiment to show the proposed algorithm’s effectiveness on weak fault signal.The used signals were generated by Intelligent Maintenance Systems of Cincinnnati, and the data set is available on the NASA website [59]. Four Rexnord ZA-2115 double row bearings were installed on a shaft, which was driven by an AC motor at the speed of 2,000 RPM via sever rub belts. The specific parameters of the experimental bearing are shown in Table 5, and the calculated frequency characteristics of the fault are shown in Table 6.Table 5
Rolling element bearing parameters of 6205-2RS.
Roller number
Roller diameter
Contact angle
Pitch diameter
16
8.407 mm
0.265 rad.
71.5 mmTable 6
Calculated characteristic frequencies of 6205-2RS (Hz).
Rotation frequencyfr
Inner racefi
Outer racefo
Roller frequencyfb
Cage frequencyfc
33.73
300.6
236.4
141.6
14.8High sensitivity quartz ICP accelerometers (PCB 353B33) were installed on the bearing housing, receiving the vibration signals, and NI DAQ Card 6062E was used to collect the signals. Each data set consists of individual files that are one second vibration signal snapshots. Each snapshot, recorded at every 10 minutes’ interval, consists of 20,480 points with the sampling frequency 20 kHz. Three tests were made. The first test was carried out for 35 days until inner race defect occurred in bearing 3 and roller element defect in bearing 4. The whole experiment was discontinuous, and the continuous running time of the bearings from normal to severe failure was about 15 days. Senanayaka et al. [62] calculated the root mean square (RMS) values of the whole snapshots. They concluded that the first 12 days of the bearings’ running duration were the normal service period, since then the bearing 3 (inner race fault) was subjected to a 2.5 days’ predegradation period, and severe failure of the bearing occurred in the last 6 hours. In our experiment, we use the first channel data of bearing 1 as normal bearing signals and the fifth channel data of bearing 3 as inner race fault signals. The normal bearing signals come from the snapshot 2003.10.31.08.59.46 (the snapshot after the bearings running for 2 days) and the following 100 snapshots of the “1st_test.” The inner fault signals come from the snapshot 2003.11.22.11.06.56 (the snapshot after the bearings running for 12.5 days) and the following 100 snapshots of the “1st_test.” The second test was carried out for 8 days until outer race failure occurred in bearing 1. Senanayaka et al. [62] concluded that the first 3.5 days of the bearings’ running duration were the normal service period, since then bearing 1 (outer race fault) was subjected to a 3.5 days’ pre-degradation period, and the severe failure of the bearing occurred in the last 3.5 hours. We use the first channel data of the snapshots from the snapshot of 2004.02.15.10.32.39 (the snapshot after the bearings running for 3 days) and the following 100 snapshots of the “2st_test” as outer fault signals. Hence, there are totally 300 samples representing three types of incipient failures (normal bearing, inner race fault, and outer race fault). Figure 15 shows the three kinds of signals with the first line of the Hankel matrix, derived from the snapshots 2003.10.31.08.59.46, 2003.11.22.11.06.56, and 2004.02.15.10.32.39, respectively. Figure 16 shows the FT spectrum of the signals in Figure 15. Multiscale decomposition modes of each signal are shown in Figures 17–19, and their FT spectrums are shown in Figures 20–22.Figure 15
The time domain diagram of the training signals, where (a), (b), and (c) represent normal bearing signal, inner race fault signal, and outer race fault signal come from the snapshots 2003.10.31.08.59.46, 2003.11.22.11.06.56, and 2004.02.15.10.32.39, respectively. Each signal is listed with the sample number 5120.
(a)
(b)
(c)Figure 16
The frequency domain diagram of the signals in Figure15.
(a)
(b)
(c)Figure 17
Multiresolution DMD modes of the normal roller bearing signal (NASA, snapshot 2003.10.31.08.59.46). Here, we select the first column ofXijDMD that represents the fourth layer slow or fast modes. (a)–(h) represent the modes of X41DMD, X42DMD, X43DMD, X44DMD, X45DMD, X46DMD, X47DMD, and X48DMD, respectively. The symbols in the Figures 18 and 19 share the same meaning.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 18
Multiresolution DMD modes of roller bearing’s inner race fault signal (NASA, snapshots 2003.11.22.11.06.56).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 19
Multiresolution DMD modes of roller bearing’s outer race fault signal (NASA, snapshot 2004.02.15.10.32.39).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 20
The Frequency domain diagram of multiresolution DMD modes of normal roller bearing signal, where (a)–(h) are corresponding to signals in Figure17, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 21
The Frequency domain diagram of multiresolution DMD modes of the inner race fault signal, where (a)–(h) are corresponding to signals in Figure18, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 22
The Frequency domain diagram of multiresolution DMD modes of the outer race fault signal, where (a)–(h) are corresponding to signals in Figure19, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)ApEn algorithm (τ=2, γ=0.2SD) is applied on the multi-resolution DMD modes, and a total of 300 vectors are obtained for the three bearing fault states. 60 vectors of each state are randomly selected as the inputs of BPNN, the remaining 40 vectors are selected for testing. Numbers of records per bearing state and ANN outputs are summarized in Table 7. The parameters of the classifier are set as the same as 4.1. The classification results of the training and testing signals are shown as in Figure 23. The category labels 1, 2, and 3 represent normal bearing, inner race fault, outer race fault, respectively. The average test accuracy of our algorithm is 98.7%.Table 7
Number of records per bearing state and BPNN outputs.
Bearing state
Normal
Inner race fault
Outer race fault
Total records
100
100
100
Training records
60
60
60
Testing records
40
40
40
BPNN outputs
100
010
001Figure 23
The classification results of roller bearing signals.


## 4.3. Fault Diagnosis of Gear Vibration Signals

Figure24 demonstrates the gear signal testing device consisting of an AC motor, couplings, dynamometers, a magnetic powder brake, and a gearbox. A pair of meshing gears is internally fitted in the gearbox. The numbers of teeth for the driving and driven wheel are 20 and 37, and the gear modulus is equal to 3. The motor speed is maintained at 1154 r/min. A sensor is displayed on the upside of the gearbox cover. Vibration signals of four different faults (normal gear, broken teeth fault, pitch errors fault, and wear fault) are collected with the sampling frequency 2000 Hz. 19 group data for each meshing pair were gathered under the brake load of 10 N/m, and each data sample has 12,000 data points.Figure 24
Gear signal testing device. (a) Physical object of the device and (b) structure diagram of test device: (1) AC motor; (2) coupling; (3) gearbox; (4) dynamometers; and (5) magnetic powder brake.
(a)
(b)Figure25 shows the signals with the first line of Hankel matrices, with only one representation per fault form listed. Figure 26 shows the FT spectrum of the signals in Figure 25. Multiscale decomposition modes of them are presented from Figures 27–30, and their FT spectrums are shown in Figures 31–34. All parameters in the computational process are selected the same as the roller bearing’s signal experiment, except for the sample numbers.Figure 25
The time domain diagram of the collected signals, where (a)–(d) represent normal gear signal, broken teeth fault signal, pitch errors signal, and wear gear signal, respectively. Each signal is listed with 3000 points.
(a)
(b)
(c)
(d)Figure 26
The frequency domain diagram of the signals in Figure25, where FGM is the gear mesh frequency.
(a)
(b)
(c)
(d)Figure 27
Multiresolution DMD modes of the normal gear signal. We select the first column ofXijDMD that represents the fourth layer slow or fast modes. (a)–(h) represent the modes of X41DMD, X42DMD, X43DMD, X44DMD, X45DMD, X46DMD, X47DMD, and X48DMD, respectively. The symbols in Figures 28–30 share the same meanings.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 28
Multiresolution DMD modes of the broken teeth fault signal.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 29
Multiresolution DMD modes of the pitch error gear fault signal.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 30
Multiresolution DMD modes of the wear gear fault signal.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 31
The frequency domain diagram of multiresolution DMD modes of the normal gear signal, where (a)–(h) are corresponding to signals in Figure27, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 32
The frequency domain diagram of multiresolution DMD modes of the broken teeth fault signal, where (a)–(h) are corresponding to signals of Figure28, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 33
The Frequency domain diagram of multiresolution DMD modes of the pitch errors signal, where (a)–(h) are corresponding to signals in Figure29, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 34
The frequency domain diagram of multiresolution DMD modes of and wear gear signal, where (a)–(h) are corresponding to signals in Figure30, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)As in Figure26, although the spectrums of the four gears’ fault signals are visually different, it is difficult to identify the characteristic frequencies of different faults with the perspective of amplitude and AM/FM components. In Figures 31–34, the gear meshing frequencies are distributed in multiple multiresolution modes, and similarly, the series of AM/FM components cannot be distinguished. It is hard to judge the gear fault type either by the FT spectrums or the multiresolution DMD modes. Therefore, we apply ApEn on the multiresolution modes to calculate their complexity, and take the ApEn vector as the characteristic information of each signal.ApEn algorithm (τ=2, γ=0.2SD) is applied on the multiresolution DMD modes, and a total of 76 vectors are obtained for the four gear states. 9 vectors of each gear state type are randomly selected as the inputs of BPNN, the remaining are selected for testing. Numbers of records per gear state and ANN outputs are summarized in Table 8. The parameters of the classifier are set the same as the parameters of the roller bearing’s experiment. The classification results of the training and testing signals are shown as in Figure 35. The category labels 1, 2, 3, and 4 represent normal gear, broken teeth fault, pitch errors fault, and wear fault, respectively. The average test accuracy of our algorithm is 97.5%.Table 8
Number of records per gear state and BPNN outputs.
Bearing state
Normal
Broken teeth
Pitch errors
Wear gear
Total records
19
19
19
19
Training records
9
9
9
9
Testing records
10
10
10
10
BPNN outputs
1000
0100
0010
0001Figure 35
The classification results of gear signals.


## 4.4. Result Discussion

Table1 summarizes the fault diagnosis methods and the algorithm performance based on SVM and ANN using data from the CWRU and NASA bearing databases. Since mr-tlsDMD is essentially an algorithm based on SVD and mode decomposition, we only compare the proposed method with the literatures that adopted mode decomposition algorithm in the step of feature extraction and selection. The comparisons of the CWRU database are shown in Table 9. In Ref. [65], LMD and sample entropy were used to extract features and SVM was applied to classify the four bearing states, and finally 100% accuracy was achieved. Tian et al. [69] integrated LMD, SVD, and BPNN/SVM/ELM to distinguish ten bearing states and got the testing accuracy with 97.7% (BPNN), 98.8% (SVM), and 99.3 (ELM), respectively. Zhao et al. [68] used EMD sample entropy of the first eight IMFs to get features and ISFLA to classify six bearing states and achieved 94.7% accuracy. Lei et al. [66] proposed a method combining time-frequency domain features and EMD energy entropy to extract features and ANFIS to diagnose faults, obtaining a classification accuracy of 95.4% for faulty bearing and 100% for healthy bearing. The comparisons of the NASA database are shown in Table 10. Jaouher et al. [53] introduced a feature extraction method based on EMD energy entropy. The first eight IMFs’ energy entropies were used to extract features, ANN was applied to classify seven bearing states, and finally 93% accuracy was achieved.Table 9
Classification comparisons of CWRU database-adopted mode decomposition algorithm.
Ref.
Main features
Bearing states
Training samples
Testing samples
Test accuracy (%)
[65]
LMD, sample entropy, energy ratio + SVM
Four
60
20
100
[69]
LMD-SVD + BPNN/SVM/ELM
Ten
172
228
97.7, 98.8, 99.3
[68]
EMD sample entropy of first ten IMFs + ISFLA
Six
140
70
100 for health, others 95.4
[66]
Time and frequency, EMD energy entropy + ANFIS
Even
140
70
94.7
Proposed
mr-tlsDMD, ApEn + BPNN
Four
120
80
100Table 10
Classification comparisons of NASA database-adopted mode decomposition algorithm.
Ref.
Main features
Bearing states
Training samples
Testing samples
Test accuracy (%)
[53]
Eight IMFs of EMD, energy entropy + SVM
Seven
2728
2666
93
Proposed
mr-tlsDMD, ApEn + BPNN
Three
180
120
98.7Ge et al. [75] present a joint fault diagnosis scheme via tensor nuclear norm canonical polyadic decomposition (TNNCPD) and multiscale permutation entropy (MSPE) to classify four gear states. BPNN isused to perform fault classification. The results illustrate that their proposed scheme can accurately recognize different gear working conditions, while the performance of WT is worse. It is worth noting that Ge et al. [75] adopt the same set of signals in the experiment as we used in this paper. The comparisons between the performance in [75] and our proposed method are shown in Table 11. The test accuracy of our result is not as perfect as that in [75], nevertheless it is significantly better than WT.Table 11
Classification comparisons of the gear database (brake load: 10 N/m).
Ref.
Main features
Bearing states
Training samples
Testing samples
Accuracy (%)
Train
Test
[75]
TNNCPD, MSPE + BPNN
Four
5
14
100
100
[75]
WT, MSPE + BPNN
Four
5
14
90
76.8
Proposed
mr-tlsDMD, ApEn + BPNN
Four
10
9
100
97.5The computational efficiency is as important as the result’s accuracy for an emerging algorithm, especially when the algorithm is adopted on the online monitoring system. Therefore, the computational efficiency of the proposed algorithm is compared with that of other mode-based methods. Only the computation time of mode extraction procedure is concerned, without considering the computation time of ApEn and BPNN. The experiments were completed on a desktop computer with Matlab version 2014a, and the basic configuration of the computer was Intel® core™ i7-4770, CPU, 3.4 GHz, and memory capacity 32G. In the process of experiments, all comment language and result processing statements were removed. Table12 shows the computational efficiency of several mode-based algorithms for NASA snapshot 2003.10.31.08.59.46, which was used as a normal bearing signal in Section 4.2. It can be seen from the table that the calculation time of EMD and LMD is significantly shorter, that of VMD is 51.67 seconds, while that of mr-tlsDMD is 445.23 seconds. It is obvious that the proposed method is not superior in computational efficiency. In the conclusion part, we will analyze the reasons for the long calculation time.Table 12
Computational efficiency of several mode-based algorithms for the NASA snapshot (s.).
Method
EMD
LMD
VMD
mr-tlsDMD (L=4)
Time
0.87
1.43
51.67
445.23Three experiments are implemented to recognize the fault type of mechanical vibration signals with the proposed algorithm. The first two sets of experiment’s signals come from CWRU and NASA, which are public bearing database, and the third come from a gear test bench. Experiments, particularly on the early bearing fault signals of Section4.2 and complex gear fault signals of Section 4.3, show that the proposed fault diagnosis method has excellent classification performance. The proposed algorithm is robust whether dealing with bearing signals or gear signals; however, it is not as efficient as other modal decomposition algorithms.


## 5. Conclusions

This paper proposes a joint fault diagnosis scheme for fault classification, which exploits the strengths of DMD’s nature of being equation-free and data-driven. Multiresolution DMD algorithm decomposes the original signal into multiscale spatiotemporal modes. These slow and fast modes jointly represent the dynamic characteristics of the original system. ApEn is performed to calculate the complexity of the multiscale spatiotemporal modes. In processing mechanical signals of different fault forms, part of entropy vectors is randomly selected as training data for the categorizer of BPNN and the rest as test data. The proposed algorithm was applied to the roller bearing signals of CWRU, NASA, and the gear signals of authors’ apparatus, and the total test accuracy of the proposed algorithm was 100%, 98.7%, and 97.5%, respectively. The results show that the proposed algorithm has excellent performance in fault diagnosis.This paper proposes an algorithm for fault diagnosis based on mr-tlsDMD and ApEn. In the implementation of the algorithm, multiple parameters need to be defined in advance, such as the number of multiscale decomposition layers, the selection of truncated rank, and the frequency threshold in the process of each layer’s decomposition. In order to minimize manual intervention, we adopt fixed parameters. The proposed algorithm with the chosen parameters was demonstrated to be effective. Different parameter selection methods according to different signal characteristics are the primary research direction in the future study. We have started to study the application of the algorithm to the piezoelectric ceramic signal and acoustic emission signal in the health monitoring of bolted structure, and the work of signal acquiring is in progress.Since DMD algorithm is essentially based on SVD and mode decomposition, the process of matrix decomposition with SVD is time consuming. Furthermore, mrDMD decomposition is proposed to obtain the dynamic characteristic information of the original mechanical signal with L-th SVD algorithm. Therefore, the method proposed in this paper does not have an advantageous computational time. Another future research may focus on the acceleration algorithm-based mr-tlsDMD, such as adopting the method of compressed sensing and sub-Nyquist sampling.Last but not the least, as it was stated in the introduction, DMD algorithm is a mode decomposition technique, similar to the algorithm of EMD/LMD/VMD. The remarkable advantage of DMD is that it can obtain a series of single-frequency modes corresponding to time scale. There is determinately no mode aliasing in DMD. However, in the practical application of the algorithm, it needs some strategies to identify the fault frequency in the series of modes, especially when the geometric parameters of bearings or/and gears are inaccurate. This paper adopts the multiresolution representation strategy to decompose the original one-dimensional signals into multiple groups of signals according to certain rules, then calculate the complexity (ApEn) of the multiple groups, respectively, and finally take the ApEn vector as the statistical characteristics of the original signals to distinguish the fault type by BPNN. Research on the more efficient strategy to make full use of the series of single-frequency modes is the core problem to adopt DMD in fault diagnosis, especially in early fault diagnosis.



---

*Source: 3262818-2019-10-30.xml*