# An MCM-Enhanced Compressive Sensing for Weak Fault Feature Extraction of Rolling Element Bearings under Variable Speeds

Authors: Ya He; Kun Feng; Minghui Hu; Jinmiao Cui
Journal: Shock and Vibration
 (2020)

---

## Abstract

The compressive sensing (CS) theory provides a new slight to the big-data problem led by the Shannon sampling theorem in rolling element bearings condition monitoring, where the measurement matrix of CS tends to be designed by the random matrix (RM) to preserve the integrity of signal roughly. However, when the signal to be analyzed is infected with strong noise, not only does the signal become insufficiently sparse, but the randomness of the measurement matrix will bring down the sensing efficiency, resulting in the loss of fault feature. Thus, a sensing-enhanced CS scheme based on a series of modes after VMD decomposition is proposed under this paper. The core of this scheme is as follows: (1) the principal mode of VMD with better sparsity replaces the raw signal for compressive sensing; (2) all these modes contain the time-frequency characteristics of the raw signal; (3) a new measurement matrix called mode-circulant matrix (MCM) is defined by circulating the mode matrix, and when the amount of samples is shrunk, the sensing efficiency can be enhanced greatly. Besides, considering the fault signal of rolling bearings under variable speed, there is a need to use order tracking to overcome the nonstationarity of the signal before applying CS theory. The analysis results of simulation and experiment prove that the VMD- and MCM-based CS can successfully extract the weak fault feature of rolling bearings with operating speed changing.


---

## Body

## 1. Introduction

The condition monitoring of rolling element bearings is a crucial method used in a wide range of helping to prevent unpredictable failure. So far, vibration-based diagnosis is a valid and frequently used pattern of bearing health assessment [1, 2], and the increasing interest in vibration-based diagnosis has heightened the need for vibration analysis approaches suitable to extract the faulty feature of rolling bearings [3]. Of particular concern and significant are the time-domain analysis (time-synchronous averaging [4]; symptom parameters [5]), frequency domain analysis (resonance demodulation [6]; cepstrum analysis [7]), and time-frequency analysis (wavelet transform [8]; empirical mode decomposition [9, 10]).The cited literature demonstrated that the fault diagnosis of rolling bearings is explored increasingly. However, since almost all of the above analysis methods are based on signal sampling without loss, high sampling frequencies must be obliged to ensure the accurate acquisition of bearing fault feature under the Shannon sampling theorem [11]. It subsequently will generate enormous amounts of data exerting more pressure on the data storage, processing, and remote transmission in continuous condition monitoring of operating equipment. Therefore, it would be preferable for its working state to be identified quickly from less transferred data. This means when the equipment fails, it takes less time to diagnose and provide suggestions averting evolving into a fatal accident.Compressive sensing (CS) has provided a fine insight for tackling this challenge. This method was primarily proposed by Candès et al. and Donoho [12, 13], whose main goal was to reconstruct the original signals from a small number of samples using sparse representation and a well-designed measurement matrix. The CS theory has been employed to image processing, nuclear magnetic resonance, remote sensing, radar, and so on. Built on CS, Chan et al. [14] implemented compression and reconstruction of an image under a small number of measurements, and the collection efficiency improved. Lustig et al. [15] studied that the time of data acquisition can be sharply reduced by MRI through CS technology. Xie [16] performed range compression by using the CS technique and a two-dimensional imaging scheme was established.Owing to the excellent effect of CS theory in the application, it has attracted extensive attention from scholars in the field of fault diagnosis. Li et al. [17] attempted to compress and reconstruct monitoring signals of train rolling bearing; the results showed that the fault can be identified from the reconstructed signal. Zhang et al. [18] implemented a bearing fault detection method based on compressed sampling and matching track. As introduced by Chenet al. [19], sparse extraction of impulse by an adaptive dictionary algorithm based on the sparse model of CS can extract impulse components effectively. Tang et al. [20] developed a sparse representation classification strategy for condition evaluation of rotating machinery. Wang et al. [21] utilized a TQWT-based CS method to overcome the problem of inadequate sparsity and extract the faulty feature of rolling bearings. The aforementioned studies, based on CS theory, aimed at promoting the sparsity or reconstruction performance, but few kinds of literature are linked to measurement matrix optimization.The measurement matrix is at the heart of signal compression and reconstruction of CS. Current literature on measurement matrix construction mainly focuses on the deterministic measurement matrix. Depending on the energy of each part of the signal, Guo et al. [22] designed a matching deterministic measurement matrix based on the principle of observation adaptive allocation. Based on the characteristics of the signal, Abtahi and Modarreshashemi [23] used a large number of comparable training signals to guide the construction of deterministic measurement matrixes. To optimize the measurement matrix, S.H.S. Masoumian and Tazehkand [24] proposed that the measurement matrix selected in the next iteration can be given in judging differential entropy reduction speed. In the literature about voice signal and image analysis, the results show that, compared with the random measurement matrix, the sensing efficiency of the deterministic measurement matrix has been greatly improved, which demonstrates the possibility of constructing a deterministic measurement matrix in the application of CS in fault diagnosis. When the CS is applied in fault diagnosis, random measurement matrices are mostly employed as dimension reduction to preserve the integrity of the signal with high probability [25]. However, owing to the randomness, useful and useless information will be acquired with equal probability resulting in low sensing efficiency. Besides, the fault diagnosis of rolling bearing in practical application still has the following challenges [2]: (1) the amplitude of vibration induced by initial faults in rolling bearings is very small; (2) the high natural frequency excited by faults of rolling bearings attenuates rapidly from the faulty bearing to vibration sensors; (3) due to complex transfer paths and complicated working environments, lots of noises are mixed in vibration signals collected in practice. All these factors combined will lead to the vibration caused by bearing failure which is extremely weak, in which the fault feature of bearing may be lost but the real redundant information may be retained with a random measurement matrix.Therefore, a sensing enhanced method named the mode-circulant matrix (MCM) is utilized in this study to increase sensing efficiency on vibration signals, which are advantageous to the weak fault feature of rolling bearings based on CS theory. Assisted by VMD [26, 27], the raw signal can be decomposed into a series of modes that contain the local characteristics of the original signal at different time scales. The MCM is constructed by the modes, and it can achieve higher sensing efficiency than the RM when CS is performed. Furthermore, it is well established that the quality of the CS application can be affected by the sparse representation of the signal. However, the inclusion of strong noise in the signal will seriously decrease the sparsity of the signal in the sparse dictionary. To promote the sparsity, the main mode is determined from a series of modes through VMD based on the principle of energy ratio [28], which can eliminate noise and highlight the faulty feature of bearing. Besides, with the help of order tracking, this method can be extended to extract fault feature of nonstationary vibration signals of rolling bearings running at variable speed. So a weak fault feature extraction scheme based on VMD and CS for rolling element bearings at variable speed is proposed.The structure of this paper is arranged as follows. Section2 presents the theoretical bases of VMD and CS. Section 3 describes the proposed scheme which adopts VMD and CS to construct the MCM and extract the weak fault feature of rolling bearings. Moreover, the detailed procedures of this method are discussed in this section. Section 4 establishes the model of defect signal collected at different locations of rolling bearing and simulates the fault signal of the inner race defect. Then, the effectiveness of the proposed scheme is demonstrated using simulated and experimental signals of rolling bearing faults in Sections 4 and 5, respectively. Section 6 presents the conclusions of the study.


## 2. Basic Concept of the Proposed Feature Extraction Scheme

### 2.1. Variational Mode Decomposition

VMD is an adaptive signal decomposition method. It can nonrecursively decompose the vibration signal into several quasiorthogonal modes. The optimal solution of the variational model can be used to obtain the frequency center and bandwidth of each mode and adaptively separate the mode components [26, 27].To estimate the bandwidth of each mode, the constrained variational problem is constructed as shown in(1)minuk,ωk∑k∂tδt+jπt⋅ukte−jωkt22s.t.∑kuk=f,where uk:=u1,...,uk and ωk:=ω1,...,ωk represent a series of modes and corresponding center frequencies, respectively.The constrained problem can be simplified to the unconstrained variational problem by introducing quadratic penalty factorα and Lagrange multiplication operator λ. The augmented Lagrange multiplier is given as follows:(2)Luk,ωk,λ≔α∑k∂tδt+jπt⋅ukte−jωkt22+ft−∑kukt22+λt,ft−∑kukt.Equation (2) can be solved by the alternate direction method of multipliers. In other words, the optimal solution of the constrained variational problem can be found through updating ukn+1,ωkn+1,λn+1 alternately.The procedure of VMD can be described as follows:(1)
Initializeu^k1,ωk1,λ1,n←0.(2)
Repeatn=n+1.(3)
k=k+1, until k=K. Update uk:(3)u^kn+1ω⟵fω−∑i<ku^inω+λnω/21+2αω−ωkn2.Updateωk:(4)ωkn+1←∫0∞ωu^kn+1ω2dω∫0∞u^kn+1ω2dω.(4)
Updateλ:(5)λn+1ω←λnω+τfω−∑ku^kn+1ω.(5)
Repeat step (2)∼(4) until the convergence∑ku^kn+1−u^kn22/u^nk22<ε. At the end of the cycle, K components are output.


### 2.2. Compressive Sensing

If a sparse signal or the signal is sparse in a transform domain, the original signal can be recovered with high precision through the signal of a small number of samples [13].If an original signal can be expanded under an orthogonal base, it can be shown as(6)x=∑i=1NαiΨi,wherex∈RN×1 denotes the original signal, Ψ=Ψ1,Ψ2,…,ΨN∈RN×Nrepresents the transform domain, and α is a sparse vector representing the transform coefficient corresponding to Ψ. The matrix form of (6) is as follows:(7)x=Ψα.It is assumed that the sparse vector includesS nonzero (or far from zero) coefficients and S is much less than N, and x is S-sparse. S is the sparsity of x in the orthogonal base.x can be projected from high dimension onto low dimension, and a compression signal y can be obtained:(8)y=Φx,where Φ is M×N-dimensional measurement matrix and M is much less than N. N×1-dimensional signal can be down-sampled to M×1 dimension according to (8).Equation (7) is substituted into (8), and then (9) can be obtained:(9)y=Φx=ΦΨα=Aα,where A=ΦΨ∈RM×N is a sensing matrix.As shown in (9), it is a question of solving the linear equation set to recover x from compressed signal y. α can be obtained by solving (9), and x will be known from (7). The solution of (9) can be expressed as(10)argminα0s.t.y=Aα,where α0 is 0-norm of α.As dimensionM of y is much smaller than the total number of samples N, the problem described by (10) is NP-hard. Several feasible theoretical solutions for α may exist, so there is not fast solution of the equation set. In this case, an approximate solution is required. The 0-norm in (10) can be transformed into 1-norm to obtain an accurate sparse solution of α, as shown in(11)argminα1s.t.y=Aα.The sparse vectorα can be solved with little error, but the prerequisite is that A meets the RIP criterion as follows:(12)1−δα22≤Aα22≤1+δα22,where δ∈0,1. That is, the measurement matrix Φ is not related to the sparse dictionary Ψ.


## 2.1. Variational Mode Decomposition

VMD is an adaptive signal decomposition method. It can nonrecursively decompose the vibration signal into several quasiorthogonal modes. The optimal solution of the variational model can be used to obtain the frequency center and bandwidth of each mode and adaptively separate the mode components [26, 27].To estimate the bandwidth of each mode, the constrained variational problem is constructed as shown in(1)minuk,ωk∑k∂tδt+jπt⋅ukte−jωkt22s.t.∑kuk=f,where uk:=u1,...,uk and ωk:=ω1,...,ωk represent a series of modes and corresponding center frequencies, respectively.The constrained problem can be simplified to the unconstrained variational problem by introducing quadratic penalty factorα and Lagrange multiplication operator λ. The augmented Lagrange multiplier is given as follows:(2)Luk,ωk,λ≔α∑k∂tδt+jπt⋅ukte−jωkt22+ft−∑kukt22+λt,ft−∑kukt.Equation (2) can be solved by the alternate direction method of multipliers. In other words, the optimal solution of the constrained variational problem can be found through updating ukn+1,ωkn+1,λn+1 alternately.The procedure of VMD can be described as follows:(1)
Initializeu^k1,ωk1,λ1,n←0.(2)
Repeatn=n+1.(3)
k=k+1, until k=K. Update uk:(3)u^kn+1ω⟵fω−∑i<ku^inω+λnω/21+2αω−ωkn2.Updateωk:(4)ωkn+1←∫0∞ωu^kn+1ω2dω∫0∞u^kn+1ω2dω.(4)
Updateλ:(5)λn+1ω←λnω+τfω−∑ku^kn+1ω.(5)
Repeat step (2)∼(4) until the convergence∑ku^kn+1−u^kn22/u^nk22<ε. At the end of the cycle, K components are output.


## 2.2. Compressive Sensing

If a sparse signal or the signal is sparse in a transform domain, the original signal can be recovered with high precision through the signal of a small number of samples [13].If an original signal can be expanded under an orthogonal base, it can be shown as(6)x=∑i=1NαiΨi,wherex∈RN×1 denotes the original signal, Ψ=Ψ1,Ψ2,…,ΨN∈RN×Nrepresents the transform domain, and α is a sparse vector representing the transform coefficient corresponding to Ψ. The matrix form of (6) is as follows:(7)x=Ψα.It is assumed that the sparse vector includesS nonzero (or far from zero) coefficients and S is much less than N, and x is S-sparse. S is the sparsity of x in the orthogonal base.x can be projected from high dimension onto low dimension, and a compression signal y can be obtained:(8)y=Φx,where Φ is M×N-dimensional measurement matrix and M is much less than N. N×1-dimensional signal can be down-sampled to M×1 dimension according to (8).Equation (7) is substituted into (8), and then (9) can be obtained:(9)y=Φx=ΦΨα=Aα,where A=ΦΨ∈RM×N is a sensing matrix.As shown in (9), it is a question of solving the linear equation set to recover x from compressed signal y. α can be obtained by solving (9), and x will be known from (7). The solution of (9) can be expressed as(10)argminα0s.t.y=Aα,where α0 is 0-norm of α.As dimensionM of y is much smaller than the total number of samples N, the problem described by (10) is NP-hard. Several feasible theoretical solutions for α may exist, so there is not fast solution of the equation set. In this case, an approximate solution is required. The 0-norm in (10) can be transformed into 1-norm to obtain an accurate sparse solution of α, as shown in(11)argminα1s.t.y=Aα.The sparse vectorα can be solved with little error, but the prerequisite is that A meets the RIP criterion as follows:(12)1−δα22≤Aα22≤1+δα22,where δ∈0,1. That is, the measurement matrix Φ is not related to the sparse dictionary Ψ.


## 3. The Proposed Feature Extraction Scheme

Figure1 explains the method proposed in this paper and its principles. The steps of this method are listed as follows.Figure 1
Schematic of the proposed bearing weak fault feature extraction method.


### 3.1. Signal Preprocessing

#### 3.1.1. Order Tracking

Referring to order tracking steps [11], the vibration signal of rolling bearings x(t) is pretreated, by which the time-domain nonstationary signal x(n) can be converted to the angle-domain stationary signal x(θ) under variable working conditions.Since the Fourier dictionary is chosen in the following as an orthogonal dictionary, the input signal length of compressed sensing needs to match the processing capacity of the computer. When it does not fit well, order tracking may have also additional use of compressing signals without loss of bearing fault feature.


### 3.2. Signal Decomposition

#### 3.2.1. Optimization of VMD Parameters

Before the vibration signal is decomposed, we should predefine the balancing parameter of the data-fidelity constraintα and the number of the modes K. As a type of intelligence optimization algorithm, Particle Swarm Optimization (PSO) can fully consider the interaction between two parameters and find the global optimal solution without traversing all the parameter groups in the interval. Therefore, this paper uses PSO to optimize the parameters of VMD in parallel, which can automatically select the best combination of parameters, avoiding the interference of human subjective factors. Therefore, x(θ) is the input of the PSO algorithm, and the energy ratio (ESN) is used as a basis for quantitative comparison in this paper. ESN represents the ratio of energy generated in the case of bearing component failure to the total energy of vibration, and its definition is as follows:(13)ESN=Xkit2+Xkot2+Xkbt2+Xkct2∑i=1NXi2,where N is the length of the signal to be analyzed and X(i) is the amplitude of the envelope order spectrum. X(kit), X(kot), X(kbt), and X(kct) represent the energy of the characteristic order of bearing outer race, inner race, rolling element, and cage fault, respectively, and kit, kot, kbt, and kct, respectively, denote the serial number corresponding to the fault characteristic order of outer race, inner race, rolling element, and cage.


#### 3.2.2. Signal Decomposition with VMD

After the parameters are optimized, according to (1)–(5), x(θ) is decomposed by VMD to obtain a series of modes, mode 1∼mode K. These modes are involved for two uses: (1) The main modes (if mode i) will be obtained by comparing ESN of different modes for noise reduction and sparsity enhancement. If mode i (i = 1, 2, …, K) is the mode with the maximum ESN value, then it is the main mode. Later mode i, also called xd(θ), will take the place of x(θ) to participate in further compressive sensing. (2) Mode 1∼mode K will be used to construct the measurement matrix needed in CS, as detailed below.


### 3.3. MCM Construction

In this paper, the measurement matrix, called the mode-circulant matrix (MCM), is constructed by the modes decomposition of VMD, as described in the following.As is shown in (14), the mode matrix I can be expressed by mode 1∼mode K. n is the length of x(θ). The mode i is a1i,a2i,...,ani, and the superscript i represents the ith mode function:(14)I=mode1mode2⋮modei⋮modeK=a11a21…an−11an1a12a22…an−12an2⋮⋮⋮⋮⋮a1ia2i⋯an−1iani⋮⋮⋮⋮⋮a1Ka2K⋯an−1KanK=h1,h2,h3,…,hn.The measurement matrix Φ can be acquired by circulating the mode matrixI nc times [29], as shown in (15). Φ is the MCM. nc is the number of cycles: (15)The number of rows of Φ can be calculated according to (16), that is, the number of measurements m: (16)m=nc+1×K.


### 3.4. Compressive Sensing

#### 3.4.1. Compressive Sampling

According to (8), the MCM is utilized to reduce the dimensions of x(θ), which compresses the N-dimensional signal x(θ) into M-dimensional y(θ).


#### 3.4.2. Transformation Selection

In this study, the Fourier dictionary is selected for a sparse representation of signals. The reason is that not only is the VMD mode of bearing fault signal sparse in the Fourier domain, but also the sparsity of harmonic signal in the Fourier domain is fixed [12, 13, 21, 30], which is very beneficial to fault feature detection. To investigate the properties of the Fourier dictionary, the inner race fault signal of a rolling bearing was simulated in Section 4.


#### 3.4.3. Fault Feature Detection

The Fourier dictionary Ψ and the MCM Φ make up the sensing matrix A. According to (11), the sparse vector α can be detected from the compressed signal y(θ) by orthogonal match pursuit with the sensing matrix A and sparsity. The sparsity setting should be taken into account in the study.It is known that the typical fault features of rolling bearings are that the fault impact excites resonance of natural frequency, and this frequency is modulated by faulty feature frequency as the carrier frequency. Thus, when the natural frequency and feature frequency (sideband frequency) are detected simultaneously, that is, sparsity is set to 4, the faulty feature frequency can be extracted.


### 3.5. Fault Diagnosis of Rolling Bearings

#### 3.5.1. Feature Visualization

After the crucial frequencies detected, the signal recovery can be achieved with the Fourier dictionary Ψ and the sparse vectorα(θ) using (7). Due to the amplitude modulation of the recovered signal, the envelope demodulation should be applied with Hilbert transform. Then the envelope spectrum is obtained by the Fourier transform of the envelope signal to realize feature extraction.


#### 3.5.2. Fault Identification

Fault diagnosis for rolling bearings can be achieved by comparing the dominant components in the envelope spectrum with the calculation value of the fault characteristic frequency of the bearing.


## 3.1. Signal Preprocessing

### 3.1.1. Order Tracking

Referring to order tracking steps [11], the vibration signal of rolling bearings x(t) is pretreated, by which the time-domain nonstationary signal x(n) can be converted to the angle-domain stationary signal x(θ) under variable working conditions.Since the Fourier dictionary is chosen in the following as an orthogonal dictionary, the input signal length of compressed sensing needs to match the processing capacity of the computer. When it does not fit well, order tracking may have also additional use of compressing signals without loss of bearing fault feature.


## 3.1.1. Order Tracking

Referring to order tracking steps [11], the vibration signal of rolling bearings x(t) is pretreated, by which the time-domain nonstationary signal x(n) can be converted to the angle-domain stationary signal x(θ) under variable working conditions.Since the Fourier dictionary is chosen in the following as an orthogonal dictionary, the input signal length of compressed sensing needs to match the processing capacity of the computer. When it does not fit well, order tracking may have also additional use of compressing signals without loss of bearing fault feature.


## 3.2. Signal Decomposition

### 3.2.1. Optimization of VMD Parameters

Before the vibration signal is decomposed, we should predefine the balancing parameter of the data-fidelity constraintα and the number of the modes K. As a type of intelligence optimization algorithm, Particle Swarm Optimization (PSO) can fully consider the interaction between two parameters and find the global optimal solution without traversing all the parameter groups in the interval. Therefore, this paper uses PSO to optimize the parameters of VMD in parallel, which can automatically select the best combination of parameters, avoiding the interference of human subjective factors. Therefore, x(θ) is the input of the PSO algorithm, and the energy ratio (ESN) is used as a basis for quantitative comparison in this paper. ESN represents the ratio of energy generated in the case of bearing component failure to the total energy of vibration, and its definition is as follows:(13)ESN=Xkit2+Xkot2+Xkbt2+Xkct2∑i=1NXi2,where N is the length of the signal to be analyzed and X(i) is the amplitude of the envelope order spectrum. X(kit), X(kot), X(kbt), and X(kct) represent the energy of the characteristic order of bearing outer race, inner race, rolling element, and cage fault, respectively, and kit, kot, kbt, and kct, respectively, denote the serial number corresponding to the fault characteristic order of outer race, inner race, rolling element, and cage.


### 3.2.2. Signal Decomposition with VMD

After the parameters are optimized, according to (1)–(5), x(θ) is decomposed by VMD to obtain a series of modes, mode 1∼mode K. These modes are involved for two uses: (1) The main modes (if mode i) will be obtained by comparing ESN of different modes for noise reduction and sparsity enhancement. If mode i (i = 1, 2, …, K) is the mode with the maximum ESN value, then it is the main mode. Later mode i, also called xd(θ), will take the place of x(θ) to participate in further compressive sensing. (2) Mode 1∼mode K will be used to construct the measurement matrix needed in CS, as detailed below.


## 3.2.1. Optimization of VMD Parameters

Before the vibration signal is decomposed, we should predefine the balancing parameter of the data-fidelity constraintα and the number of the modes K. As a type of intelligence optimization algorithm, Particle Swarm Optimization (PSO) can fully consider the interaction between two parameters and find the global optimal solution without traversing all the parameter groups in the interval. Therefore, this paper uses PSO to optimize the parameters of VMD in parallel, which can automatically select the best combination of parameters, avoiding the interference of human subjective factors. Therefore, x(θ) is the input of the PSO algorithm, and the energy ratio (ESN) is used as a basis for quantitative comparison in this paper. ESN represents the ratio of energy generated in the case of bearing component failure to the total energy of vibration, and its definition is as follows:(13)ESN=Xkit2+Xkot2+Xkbt2+Xkct2∑i=1NXi2,where N is the length of the signal to be analyzed and X(i) is the amplitude of the envelope order spectrum. X(kit), X(kot), X(kbt), and X(kct) represent the energy of the characteristic order of bearing outer race, inner race, rolling element, and cage fault, respectively, and kit, kot, kbt, and kct, respectively, denote the serial number corresponding to the fault characteristic order of outer race, inner race, rolling element, and cage.


## 3.2.2. Signal Decomposition with VMD

After the parameters are optimized, according to (1)–(5), x(θ) is decomposed by VMD to obtain a series of modes, mode 1∼mode K. These modes are involved for two uses: (1) The main modes (if mode i) will be obtained by comparing ESN of different modes for noise reduction and sparsity enhancement. If mode i (i = 1, 2, …, K) is the mode with the maximum ESN value, then it is the main mode. Later mode i, also called xd(θ), will take the place of x(θ) to participate in further compressive sensing. (2) Mode 1∼mode K will be used to construct the measurement matrix needed in CS, as detailed below.


## 3.3. MCM Construction

In this paper, the measurement matrix, called the mode-circulant matrix (MCM), is constructed by the modes decomposition of VMD, as described in the following.As is shown in (14), the mode matrix I can be expressed by mode 1∼mode K. n is the length of x(θ). The mode i is a1i,a2i,...,ani, and the superscript i represents the ith mode function:(14)I=mode1mode2⋮modei⋮modeK=a11a21…an−11an1a12a22…an−12an2⋮⋮⋮⋮⋮a1ia2i⋯an−1iani⋮⋮⋮⋮⋮a1Ka2K⋯an−1KanK=h1,h2,h3,…,hn.The measurement matrix Φ can be acquired by circulating the mode matrixI nc times [29], as shown in (15). Φ is the MCM. nc is the number of cycles: (15)The number of rows of Φ can be calculated according to (16), that is, the number of measurements m: (16)m=nc+1×K.


## 3.4. Compressive Sensing

### 3.4.1. Compressive Sampling

According to (8), the MCM is utilized to reduce the dimensions of x(θ), which compresses the N-dimensional signal x(θ) into M-dimensional y(θ).


### 3.4.2. Transformation Selection

In this study, the Fourier dictionary is selected for a sparse representation of signals. The reason is that not only is the VMD mode of bearing fault signal sparse in the Fourier domain, but also the sparsity of harmonic signal in the Fourier domain is fixed [12, 13, 21, 30], which is very beneficial to fault feature detection. To investigate the properties of the Fourier dictionary, the inner race fault signal of a rolling bearing was simulated in Section 4.


### 3.4.3. Fault Feature Detection

The Fourier dictionary Ψ and the MCM Φ make up the sensing matrix A. According to (11), the sparse vector α can be detected from the compressed signal y(θ) by orthogonal match pursuit with the sensing matrix A and sparsity. The sparsity setting should be taken into account in the study.It is known that the typical fault features of rolling bearings are that the fault impact excites resonance of natural frequency, and this frequency is modulated by faulty feature frequency as the carrier frequency. Thus, when the natural frequency and feature frequency (sideband frequency) are detected simultaneously, that is, sparsity is set to 4, the faulty feature frequency can be extracted.


## 3.4.1. Compressive Sampling

According to (8), the MCM is utilized to reduce the dimensions of x(θ), which compresses the N-dimensional signal x(θ) into M-dimensional y(θ).


## 3.4.2. Transformation Selection

In this study, the Fourier dictionary is selected for a sparse representation of signals. The reason is that not only is the VMD mode of bearing fault signal sparse in the Fourier domain, but also the sparsity of harmonic signal in the Fourier domain is fixed [12, 13, 21, 30], which is very beneficial to fault feature detection. To investigate the properties of the Fourier dictionary, the inner race fault signal of a rolling bearing was simulated in Section 4.


## 3.4.3. Fault Feature Detection

The Fourier dictionary Ψ and the MCM Φ make up the sensing matrix A. According to (11), the sparse vector α can be detected from the compressed signal y(θ) by orthogonal match pursuit with the sensing matrix A and sparsity. The sparsity setting should be taken into account in the study.It is known that the typical fault features of rolling bearings are that the fault impact excites resonance of natural frequency, and this frequency is modulated by faulty feature frequency as the carrier frequency. Thus, when the natural frequency and feature frequency (sideband frequency) are detected simultaneously, that is, sparsity is set to 4, the faulty feature frequency can be extracted.


## 3.5. Fault Diagnosis of Rolling Bearings

### 3.5.1. Feature Visualization

After the crucial frequencies detected, the signal recovery can be achieved with the Fourier dictionary Ψ and the sparse vectorα(θ) using (7). Due to the amplitude modulation of the recovered signal, the envelope demodulation should be applied with Hilbert transform. Then the envelope spectrum is obtained by the Fourier transform of the envelope signal to realize feature extraction.


### 3.5.2. Fault Identification

Fault diagnosis for rolling bearings can be achieved by comparing the dominant components in the envelope spectrum with the calculation value of the fault characteristic frequency of the bearing.


## 3.5.1. Feature Visualization

After the crucial frequencies detected, the signal recovery can be achieved with the Fourier dictionary Ψ and the sparse vectorα(θ) using (7). Due to the amplitude modulation of the recovered signal, the envelope demodulation should be applied with Hilbert transform. Then the envelope spectrum is obtained by the Fourier transform of the envelope signal to realize feature extraction.


## 3.5.2. Fault Identification

Fault diagnosis for rolling bearings can be achieved by comparing the dominant components in the envelope spectrum with the calculation value of the fault characteristic frequency of the bearing.


## 4. Simulation

To verify the properties of the Fourier dictionary and MCM in this method, the fault signal of the inner race of the rolling bearing is constructed under speed-up conditions. The simulations are as follows [31, 32].


### 4.1. Rolling Bearings Fault Signal Modeling

The simulated signal of rolling bearings’ fault model can be expressed strictly as follows:(17)xt=∑i=1NAist−Ti−τi,where Ai is the amplitude of the ith impulse and s(t-Ti-τi) is the ith impulse oscillation generated by the fault with the natural frequency fn. Generally, s(t) is an exponential decay process. Under variable speed conditions, s(t) and Ai can be refined into(18)st=e−Bt−Ti−τi⋅sin2πfnt−Ti−τi,(19)Ai=Acos2πfrt+φA+CA.The bearing fault simulation signal can be concluded as follows:(20)xt=∑i=1NAcos2πfmt+φA+CA⋅e−Bt−Ti−τicos2πfnt−Ti−τi,where B is an attenuation coefficient, Ti is the corresponding time or the ith fault impulse, τi=λi⋅Ti is the minor random fluctuation around Ti, and n(t) represents random noise.The calculation process ofTi under variable speed condition is as follows.In this simulation, the fault characteristic order of bearing inner race isOi = 14.9; that is, the fault impact caused by the rotor every rotation cycle is 14.9. Also, it is assumed that the variation equation of rotor speed is fr(t) = at + b. The relationship between the rotating angle of the inner race and the time coordinate corresponding to the fault impulse can be expressed by(21)∫0tfrtdt=a2t2+bt=θi,where θi = i/Oi (i = 1, 2, …, K). The time coordinate Ti corresponding to the impact can be obtained by solving the above equation.This study assumesA = 2, CA = 0.5, ΦA = 30°, B = 800, fn = 3000, λi = 0.01, a = 4, b = 15, and Oi = 14.9. Besides, the fault characteristic orders of bearing inner race, rolling elements, and cage are 9.92, 2.4, and 0.4, respectively. In this paper, the sampling frequency of the simulation signal is 12.8 kHz, and the sampling length is 8192. The simulation signal of inner race defect without noise is shown, respectively, in Figure 2.Figure 2
The inner race fault simulation signal.The fault vibration signal of rolling bearing measuredxb(t) is composed of the fault signal x(t) and some random noise n(t). It can be expressed as follows:(22)xbt=xt+nt.Gaussian noise is added to the simulated signal in Section4.1 to hide the fault features in the noise. To quantify the added noise influence, two means of calculating signal-to-noise ratio (SNR) are used, one for gauging the local peaks of the random noise (referred to as SNR1) [33] and another for further gauging the noise level across fault characteristic order (referred to as SNR2) [34]. SNR1 is defined as(23)SNR1=20log10PsAn,where Ps denotes the RMS value of the signal x(t) and An represents the peak value of the noise. According to (23), the SNR1 of the fault simulation signal xb(t) is −22.7 dB.SNR2 is defined as follows:(24)SNR=10∗log10SOtNOt,SOt=SO1t+SO2t+⋯+SOit+⋯+SOnt,n=OmaxOtSOit=∑j=kit−lkit+lXj2,i=1,n,l≠0,NOt=∑Xj2−SOt,where Omax is the maximum order to be analyzed and X(j) is the amplitude of the envelope order spectrum. S(Ot) represents the total energy of bearing fault characteristic order and its multiple in the range of 0 to Omax. S(Oit) is the energy which is in the vicinity of iOt − lΔo and iOt + lΔo and it centers on iOt. Oit is the i time of theoretical fault characteristic order, and kit is the serial number of Oit. N(Ot) denotes the energy of noise in the same range and it excludes S(Ot). According to (24), the SNR2 of the fault simulation signal is xb(t) − 16.25 dB.The signal after being affected by noise is shown in Figure3(a). Additionally, Figure 3(b) provides the corresponding speed trend.Figure 3
Simulated signal of bearing fault with noise. (a) The waveform of vibration and (b) speed trend.
(a)(b)


### 4.2. Optimization of VMD Parameters

According to the steps of the proposed method in Section3, the signals Figure 3(a) are resampled in the angle domain, respectively, and the obtained angular waveform is as shown in Figure 4.Figure 4
The angular signal with noise.Particle Swarm Optimization has been applied to analyze the noisy signal in angle domain. A fitness curve can be obtained, as being presented in Figure5. The setting of the main parameters and optimization results is given in Table 1. The maximum value of the local energy ratio appears at 0.1475 in the tenth generation, and the best combination of decomposition parameters is [α, K] = [664.7, 6].Figure 5
The trend of local characteristic order amplitude with evolution algebra.Table 1
Parameters setting and optimization results of PSO.
ParametersαKRange of values[0, 2000][1,10]Maximum evolution50Maximum population10Local search speed1Global search speed1.2Optimization results664.76Subsequently, the modes are generated by the angular signal’s VMD decomposition in whichK = 6 and α = 664.7 are used. Table 2 shows the fitness values for each mode, and mode 4 has the maximum fitness value, almost 3∼5 times that of other modes. So the mode 4 is the main mode, that is, denoised signal for the simulation signal with noise.Table 2
Fitness values obtained by analyzing simulation signal with VMD.
Modes obtained by VMD and their order spectraFitness (%)3.74.264.214.754.184.33


### 4.3. Properties Verification

Based on Section4.2, this part will further validate three properties: (1) after noise reduction of rolling bearing fault signal, the sparsity is promoted in Fourier dictionary; (2) the sensing efficiency of MCM is higher than that of Gaussian random matrix (GRM); (3) MCM-based CS (MCM-CS) is more effective than Narrowband Envelope Analysis (NEA) in the feature extraction of bearing fault.


#### 4.3.1. Sparsity Analysis of Bearing Fault Signal

In the sparse transform domain, onlyK nonzero elements are sparse (where K is far less than the signal length). Generally, the signal does not necessarily satisfy the above requirement of absolute sparsity. The signal can be regarded as an approximate sparse signal when the sparse transform coefficients of the signal decrease exponentially to zero [12, 13, 30]. Therefore, we can judge whether the rolling bearing fault signal in 4.1 is sparse under the Fourier dictionary by analyzing the decreasing trend of Fourier coefficients.In this paper, firstly, the exponential function is utilized to fit the attenuation curve. Secondly, we calculate the sum of squares error (SSE) and coefficient of determination (R-squared) between the original curve and the fitted curve to assess whether the attenuation of the sparse coefficient is close to the exponential attenuation. Here, SSE is the square sum of the errors of the corresponding points of the fitting data and the original data. R-square is used to represent the quality of the fitting by the change of data. The closer it is to 1, the better the data fitting is. The closer SSE is to 0, the more successful the fitting is.The angular waveform without noise is shown in Figure4 (a), and the sparse coefficients obtained by the Fourier transform are shown in Figure 6(a). The absolute values of Fourier transform coefficients are sorted from large to small, and the sorted result is shown in Figure 6(b).Figure 6
In Fourier dictionary, (a) sparse coefficients of pure simulation signal; (b) attenuation curve of (a).
(a)(b)The attenuation curve of Figure6(b) is fitted by using the least-square method, through the curve fitting toolbox of MATLAB, and the formula obtained is fx=0.2e−0.059x+0.03e−0.0018x. The actual attenuation curve and the fitting curve are as presented in Figure 7, where SSE is 0.011 and R-square is 0.984, which are both close to the ideal value. Furthermore, from the detailed area, it is obvious that the sparse coefficients decrease to zero. So it can be considered that the Fourier coefficients decay to zero exponentially, indicating that the vibration signal of rolling bearing is almost sparse in the Fourier domain.Figure 7
The fitting result of coefficients attenuation curve in Figure6(b): SSE = 0.011; R-square = 0.984.The existence of noise in the vibration signal will increase the difficulty of sparse representation in the Fourier domain. Sparse coefficients of the noisy signal and the sequencing result are shown in Figures8(a) and 8(b), respectively. The fitting result of the attenuation curve is shown in Figure 9, and the formula of the fitted curve is fx=0.065e−0.00047x−5.5×10−12e−0.0052x. Due to the influence of noise, SSE is as high as 0.18 and R-square is only 0.852, both of which deviate slightly from the ideal value. What is more, the attenuation curvature of sparse coefficients of the noisy signal is increasing, not decreasing as exponential attenuation. Thus, sparse coefficients in Figure 8(b) cannot be concluded to attenuate to zero exponentially, which means the polluted signal has inadequate sparsity affected by full-band noise.Figure 8
In Fourier dictionary, (a) sparse coefficients of the noisy signal; (b) attenuation curve of (a).
(a)(b)Figure 9
The fitting result of coefficients attenuation curve in Figure8(b): SSE = 0.18; R-square = 0.852.Finally, the sparsity of mode 4 in Table2 is analyzed in the Fourier domain. Sparse coefficients of the denoised signal and the sorting result are shown in Figures 10(a) and 10(b), respectively. The attenuation curve of Figure 10(b) is further fitted, and the formula obtained is fx=0.16e−0.091x+0.056e−0.0031x. In Figure 11, SSE is 0.007, R-square is 0.994, and the actual attenuation curve almost overlaps the fitted curve. It can be inferred that the sparse coefficients of the denoised signal can be decreased exponentially to zero under the Fourier transform. Compared with the noisy signal, the sparsity of the denoised signal is greatly enhanced.Figure 10
In Fourier dictionary, (a) sparse coefficients of the mode 3 and (b) attenuation curve of (a).
(a)(b)Figure 11
The fitting result of coefficient attenuation curve in Figure10(b): SSE = 0.007; R-square = 0.994.


#### 4.3.2. High Sensing Efficiency of MCM

Gaussian random matrix (GRM) is frequently used as a measurement matrix of CS for a typical one-dimensional vibration signal [18–21, 25]. The reason is just simple and convenient to construct, and most importantly, the measurement matrix constructed by the basis vector obeying the Gauss distribution can be roughly guaranteed to be irrelevant to any sparse basis.In this study, a new compressive sensing measurement matrix is constructed in the form of the circulant matrix by using the VMD modes. Its sensing efficiency is much higher than GRM, which means the extraction accuracy of the MCM can be achieved by fewer measurements than GRM. The following parts will demonstrate it.This verification firstly used the modes in Table2 to construct the MCM and then applied it and GRM, respectively, to compress the denoised vibration signal. Fault feature extraction is considered successful when the fault feature order appears in the final envelope spectrum.Figure12(a) and 12(b) are the final envelope order spectra obtained by using MCM and GRM, respectively, when the number of measurements is 600. The dominant order in Figure 12(a) is consistent with the inner simulated fault order O = 14.9, representing a 100% success rate of fault feature extraction, while Figure 12(b) fails; that is, the success rate is 0. According to (24), the SNR2 of Figure 12(a) and 12(b) are 33.4 dB and −48.7 dB, respectively, and bearing fault can be directly identified with the former.Figure 12
The envelope order spectra obtained by using (a) MCM: the spectrum is clear, and SNR2 is 33.4 dB; (b) GRM: the spectrum is clear, but SNR2 is only −44.7 dB.
(a)(b)Moreover, the change in success rate is shown in Figure13 as the number of measurements increases from 20 to 1000. From Figure 13, when the number of measurements is only 60, the MCM can successfully extract bearing fault features, and the results are stable and consistent within a varied number of measurements. By contrast, the GRM has obvious randomness and poor robustness although it can occasionally extract fault features under certain measurements. Therefore, sensing efficiency of the MCM is higher than that of the GRM, which means that the beneficial information of the signal can be measured more effectively.Figure 13
Comparison of detection effect of two measuring matrices.


#### 4.3.3. Validation of NEA Method

Section4.3.2 has verified that MCM has higher sensing efficiency than GRM. Combined with CS, MCM-CS can effectively extract bearing fault features under a fewer number of measurements. To further highlight the significance of MCM-CS for bearing fault compressive detection, the Narrowband Envelope Analysis (NEA) is selected for comparison.Before analyzing the bearing fault signal to NEA, the most fundamental issue to be solved first is the setting of the center frequency and bandwidth of the band-pass filter. However, it fails to extract effective features based on experience frequently. In recent years, SK has been extensively used in the selection of demodulation band parameters and gained good reputation [35]. Thus, SK is used to design the band-pass filter in the NEA method in this section. Besides, considering variable speed conditions, NEA is carried out as follows:(1)
Even angle resampling in the speed domain, order tracking is applied to the original signal, achieving a stationary signal in the angle domain instead of a nonstationary signal in the time domain.(2)
The optimal envelope center frequency and bandwidth are determined by SK, based on which a band-pass filter is designed to filter the resampled signal.(3)
The envelope demodulation method based on the Hilbert transform is applied to achieve the envelope order spectrum.Figure14 shows the analysis results for the fault simulation signal shown in Figure 3(a) using the NEA based on SK and order tracking. Figure 14(a) is the kurtogram. The SK shows that the center frequency of the optimum demodulation band is 2933 Hz, which is extremely near to the resonance frequency of the simulated signal. This proves that SK is highly significant. Figure 14(b) shows the obscure envelope order spectrum, where the amplitude of 2.099 order is the highest. An order of 14.79, corresponding to the inner fault, can be seen but it is not the dominant component in this spectrum.Figure 14
The analysis results using NEA: (a) kurtogram; (b) envelope order spectrum: the spectrum is obscure, and SNR2 is only −9.07 dB.
(a)(b)Besides, according to (24), the SNR2 of the filtered signal is −9.07 dB, which is far from the recovered signal through MCM-CS. As observed in the above, although NEA is indeed simpler than MCM-CS, it is not only unable to achieve fault compressive detection but also inferior in weak fault feature extraction.The facts of Section4.3 prove that the VMD mode of fault vibration signals of rolling bearings can be sparsely represented in the Fourier dictionary. What is more, the MCM-enhanced compressive sensing is endowed with higher sensing efficiency. These are the necessary conditions for CS to extract fault feature of rolling bearings in Section 5.


## 4.1. Rolling Bearings Fault Signal Modeling

The simulated signal of rolling bearings’ fault model can be expressed strictly as follows:(17)xt=∑i=1NAist−Ti−τi,where Ai is the amplitude of the ith impulse and s(t-Ti-τi) is the ith impulse oscillation generated by the fault with the natural frequency fn. Generally, s(t) is an exponential decay process. Under variable speed conditions, s(t) and Ai can be refined into(18)st=e−Bt−Ti−τi⋅sin2πfnt−Ti−τi,(19)Ai=Acos2πfrt+φA+CA.The bearing fault simulation signal can be concluded as follows:(20)xt=∑i=1NAcos2πfmt+φA+CA⋅e−Bt−Ti−τicos2πfnt−Ti−τi,where B is an attenuation coefficient, Ti is the corresponding time or the ith fault impulse, τi=λi⋅Ti is the minor random fluctuation around Ti, and n(t) represents random noise.The calculation process ofTi under variable speed condition is as follows.In this simulation, the fault characteristic order of bearing inner race isOi = 14.9; that is, the fault impact caused by the rotor every rotation cycle is 14.9. Also, it is assumed that the variation equation of rotor speed is fr(t) = at + b. The relationship between the rotating angle of the inner race and the time coordinate corresponding to the fault impulse can be expressed by(21)∫0tfrtdt=a2t2+bt=θi,where θi = i/Oi (i = 1, 2, …, K). The time coordinate Ti corresponding to the impact can be obtained by solving the above equation.This study assumesA = 2, CA = 0.5, ΦA = 30°, B = 800, fn = 3000, λi = 0.01, a = 4, b = 15, and Oi = 14.9. Besides, the fault characteristic orders of bearing inner race, rolling elements, and cage are 9.92, 2.4, and 0.4, respectively. In this paper, the sampling frequency of the simulation signal is 12.8 kHz, and the sampling length is 8192. The simulation signal of inner race defect without noise is shown, respectively, in Figure 2.Figure 2
The inner race fault simulation signal.The fault vibration signal of rolling bearing measuredxb(t) is composed of the fault signal x(t) and some random noise n(t). It can be expressed as follows:(22)xbt=xt+nt.Gaussian noise is added to the simulated signal in Section4.1 to hide the fault features in the noise. To quantify the added noise influence, two means of calculating signal-to-noise ratio (SNR) are used, one for gauging the local peaks of the random noise (referred to as SNR1) [33] and another for further gauging the noise level across fault characteristic order (referred to as SNR2) [34]. SNR1 is defined as(23)SNR1=20log10PsAn,where Ps denotes the RMS value of the signal x(t) and An represents the peak value of the noise. According to (23), the SNR1 of the fault simulation signal xb(t) is −22.7 dB.SNR2 is defined as follows:(24)SNR=10∗log10SOtNOt,SOt=SO1t+SO2t+⋯+SOit+⋯+SOnt,n=OmaxOtSOit=∑j=kit−lkit+lXj2,i=1,n,l≠0,NOt=∑Xj2−SOt,where Omax is the maximum order to be analyzed and X(j) is the amplitude of the envelope order spectrum. S(Ot) represents the total energy of bearing fault characteristic order and its multiple in the range of 0 to Omax. S(Oit) is the energy which is in the vicinity of iOt − lΔo and iOt + lΔo and it centers on iOt. Oit is the i time of theoretical fault characteristic order, and kit is the serial number of Oit. N(Ot) denotes the energy of noise in the same range and it excludes S(Ot). According to (24), the SNR2 of the fault simulation signal is xb(t) − 16.25 dB.The signal after being affected by noise is shown in Figure3(a). Additionally, Figure 3(b) provides the corresponding speed trend.Figure 3
Simulated signal of bearing fault with noise. (a) The waveform of vibration and (b) speed trend.
(a)(b)


## 4.2. Optimization of VMD Parameters

According to the steps of the proposed method in Section3, the signals Figure 3(a) are resampled in the angle domain, respectively, and the obtained angular waveform is as shown in Figure 4.Figure 4
The angular signal with noise.Particle Swarm Optimization has been applied to analyze the noisy signal in angle domain. A fitness curve can be obtained, as being presented in Figure5. The setting of the main parameters and optimization results is given in Table 1. The maximum value of the local energy ratio appears at 0.1475 in the tenth generation, and the best combination of decomposition parameters is [α, K] = [664.7, 6].Figure 5
The trend of local characteristic order amplitude with evolution algebra.Table 1
Parameters setting and optimization results of PSO.
ParametersαKRange of values[0, 2000][1,10]Maximum evolution50Maximum population10Local search speed1Global search speed1.2Optimization results664.76Subsequently, the modes are generated by the angular signal’s VMD decomposition in whichK = 6 and α = 664.7 are used. Table 2 shows the fitness values for each mode, and mode 4 has the maximum fitness value, almost 3∼5 times that of other modes. So the mode 4 is the main mode, that is, denoised signal for the simulation signal with noise.Table 2
Fitness values obtained by analyzing simulation signal with VMD.
Modes obtained by VMD and their order spectraFitness (%)3.74.264.214.754.184.33


## 4.3. Properties Verification

Based on Section4.2, this part will further validate three properties: (1) after noise reduction of rolling bearing fault signal, the sparsity is promoted in Fourier dictionary; (2) the sensing efficiency of MCM is higher than that of Gaussian random matrix (GRM); (3) MCM-based CS (MCM-CS) is more effective than Narrowband Envelope Analysis (NEA) in the feature extraction of bearing fault.


### 4.3.1. Sparsity Analysis of Bearing Fault Signal

In the sparse transform domain, onlyK nonzero elements are sparse (where K is far less than the signal length). Generally, the signal does not necessarily satisfy the above requirement of absolute sparsity. The signal can be regarded as an approximate sparse signal when the sparse transform coefficients of the signal decrease exponentially to zero [12, 13, 30]. Therefore, we can judge whether the rolling bearing fault signal in 4.1 is sparse under the Fourier dictionary by analyzing the decreasing trend of Fourier coefficients.In this paper, firstly, the exponential function is utilized to fit the attenuation curve. Secondly, we calculate the sum of squares error (SSE) and coefficient of determination (R-squared) between the original curve and the fitted curve to assess whether the attenuation of the sparse coefficient is close to the exponential attenuation. Here, SSE is the square sum of the errors of the corresponding points of the fitting data and the original data. R-square is used to represent the quality of the fitting by the change of data. The closer it is to 1, the better the data fitting is. The closer SSE is to 0, the more successful the fitting is.The angular waveform without noise is shown in Figure4 (a), and the sparse coefficients obtained by the Fourier transform are shown in Figure 6(a). The absolute values of Fourier transform coefficients are sorted from large to small, and the sorted result is shown in Figure 6(b).Figure 6
In Fourier dictionary, (a) sparse coefficients of pure simulation signal; (b) attenuation curve of (a).
(a)(b)The attenuation curve of Figure6(b) is fitted by using the least-square method, through the curve fitting toolbox of MATLAB, and the formula obtained is fx=0.2e−0.059x+0.03e−0.0018x. The actual attenuation curve and the fitting curve are as presented in Figure 7, where SSE is 0.011 and R-square is 0.984, which are both close to the ideal value. Furthermore, from the detailed area, it is obvious that the sparse coefficients decrease to zero. So it can be considered that the Fourier coefficients decay to zero exponentially, indicating that the vibration signal of rolling bearing is almost sparse in the Fourier domain.Figure 7
The fitting result of coefficients attenuation curve in Figure6(b): SSE = 0.011; R-square = 0.984.The existence of noise in the vibration signal will increase the difficulty of sparse representation in the Fourier domain. Sparse coefficients of the noisy signal and the sequencing result are shown in Figures8(a) and 8(b), respectively. The fitting result of the attenuation curve is shown in Figure 9, and the formula of the fitted curve is fx=0.065e−0.00047x−5.5×10−12e−0.0052x. Due to the influence of noise, SSE is as high as 0.18 and R-square is only 0.852, both of which deviate slightly from the ideal value. What is more, the attenuation curvature of sparse coefficients of the noisy signal is increasing, not decreasing as exponential attenuation. Thus, sparse coefficients in Figure 8(b) cannot be concluded to attenuate to zero exponentially, which means the polluted signal has inadequate sparsity affected by full-band noise.Figure 8
In Fourier dictionary, (a) sparse coefficients of the noisy signal; (b) attenuation curve of (a).
(a)(b)Figure 9
The fitting result of coefficients attenuation curve in Figure8(b): SSE = 0.18; R-square = 0.852.Finally, the sparsity of mode 4 in Table2 is analyzed in the Fourier domain. Sparse coefficients of the denoised signal and the sorting result are shown in Figures 10(a) and 10(b), respectively. The attenuation curve of Figure 10(b) is further fitted, and the formula obtained is fx=0.16e−0.091x+0.056e−0.0031x. In Figure 11, SSE is 0.007, R-square is 0.994, and the actual attenuation curve almost overlaps the fitted curve. It can be inferred that the sparse coefficients of the denoised signal can be decreased exponentially to zero under the Fourier transform. Compared with the noisy signal, the sparsity of the denoised signal is greatly enhanced.Figure 10
In Fourier dictionary, (a) sparse coefficients of the mode 3 and (b) attenuation curve of (a).
(a)(b)Figure 11
The fitting result of coefficient attenuation curve in Figure10(b): SSE = 0.007; R-square = 0.994.


### 4.3.2. High Sensing Efficiency of MCM

Gaussian random matrix (GRM) is frequently used as a measurement matrix of CS for a typical one-dimensional vibration signal [18–21, 25]. The reason is just simple and convenient to construct, and most importantly, the measurement matrix constructed by the basis vector obeying the Gauss distribution can be roughly guaranteed to be irrelevant to any sparse basis.In this study, a new compressive sensing measurement matrix is constructed in the form of the circulant matrix by using the VMD modes. Its sensing efficiency is much higher than GRM, which means the extraction accuracy of the MCM can be achieved by fewer measurements than GRM. The following parts will demonstrate it.This verification firstly used the modes in Table2 to construct the MCM and then applied it and GRM, respectively, to compress the denoised vibration signal. Fault feature extraction is considered successful when the fault feature order appears in the final envelope spectrum.Figure12(a) and 12(b) are the final envelope order spectra obtained by using MCM and GRM, respectively, when the number of measurements is 600. The dominant order in Figure 12(a) is consistent with the inner simulated fault order O = 14.9, representing a 100% success rate of fault feature extraction, while Figure 12(b) fails; that is, the success rate is 0. According to (24), the SNR2 of Figure 12(a) and 12(b) are 33.4 dB and −48.7 dB, respectively, and bearing fault can be directly identified with the former.Figure 12
The envelope order spectra obtained by using (a) MCM: the spectrum is clear, and SNR2 is 33.4 dB; (b) GRM: the spectrum is clear, but SNR2 is only −44.7 dB.
(a)(b)Moreover, the change in success rate is shown in Figure13 as the number of measurements increases from 20 to 1000. From Figure 13, when the number of measurements is only 60, the MCM can successfully extract bearing fault features, and the results are stable and consistent within a varied number of measurements. By contrast, the GRM has obvious randomness and poor robustness although it can occasionally extract fault features under certain measurements. Therefore, sensing efficiency of the MCM is higher than that of the GRM, which means that the beneficial information of the signal can be measured more effectively.Figure 13
Comparison of detection effect of two measuring matrices.


### 4.3.3. Validation of NEA Method

Section4.3.2 has verified that MCM has higher sensing efficiency than GRM. Combined with CS, MCM-CS can effectively extract bearing fault features under a fewer number of measurements. To further highlight the significance of MCM-CS for bearing fault compressive detection, the Narrowband Envelope Analysis (NEA) is selected for comparison.Before analyzing the bearing fault signal to NEA, the most fundamental issue to be solved first is the setting of the center frequency and bandwidth of the band-pass filter. However, it fails to extract effective features based on experience frequently. In recent years, SK has been extensively used in the selection of demodulation band parameters and gained good reputation [35]. Thus, SK is used to design the band-pass filter in the NEA method in this section. Besides, considering variable speed conditions, NEA is carried out as follows:(1)
Even angle resampling in the speed domain, order tracking is applied to the original signal, achieving a stationary signal in the angle domain instead of a nonstationary signal in the time domain.(2)
The optimal envelope center frequency and bandwidth are determined by SK, based on which a band-pass filter is designed to filter the resampled signal.(3)
The envelope demodulation method based on the Hilbert transform is applied to achieve the envelope order spectrum.Figure14 shows the analysis results for the fault simulation signal shown in Figure 3(a) using the NEA based on SK and order tracking. Figure 14(a) is the kurtogram. The SK shows that the center frequency of the optimum demodulation band is 2933 Hz, which is extremely near to the resonance frequency of the simulated signal. This proves that SK is highly significant. Figure 14(b) shows the obscure envelope order spectrum, where the amplitude of 2.099 order is the highest. An order of 14.79, corresponding to the inner fault, can be seen but it is not the dominant component in this spectrum.Figure 14
The analysis results using NEA: (a) kurtogram; (b) envelope order spectrum: the spectrum is obscure, and SNR2 is only −9.07 dB.
(a)(b)Besides, according to (24), the SNR2 of the filtered signal is −9.07 dB, which is far from the recovered signal through MCM-CS. As observed in the above, although NEA is indeed simpler than MCM-CS, it is not only unable to achieve fault compressive detection but also inferior in weak fault feature extraction.The facts of Section4.3 prove that the VMD mode of fault vibration signals of rolling bearings can be sparsely represented in the Fourier dictionary. What is more, the MCM-enhanced compressive sensing is endowed with higher sensing efficiency. These are the necessary conditions for CS to extract fault feature of rolling bearings in Section 5.


## 4.3.1. Sparsity Analysis of Bearing Fault Signal

In the sparse transform domain, onlyK nonzero elements are sparse (where K is far less than the signal length). Generally, the signal does not necessarily satisfy the above requirement of absolute sparsity. The signal can be regarded as an approximate sparse signal when the sparse transform coefficients of the signal decrease exponentially to zero [12, 13, 30]. Therefore, we can judge whether the rolling bearing fault signal in 4.1 is sparse under the Fourier dictionary by analyzing the decreasing trend of Fourier coefficients.In this paper, firstly, the exponential function is utilized to fit the attenuation curve. Secondly, we calculate the sum of squares error (SSE) and coefficient of determination (R-squared) between the original curve and the fitted curve to assess whether the attenuation of the sparse coefficient is close to the exponential attenuation. Here, SSE is the square sum of the errors of the corresponding points of the fitting data and the original data. R-square is used to represent the quality of the fitting by the change of data. The closer it is to 1, the better the data fitting is. The closer SSE is to 0, the more successful the fitting is.The angular waveform without noise is shown in Figure4 (a), and the sparse coefficients obtained by the Fourier transform are shown in Figure 6(a). The absolute values of Fourier transform coefficients are sorted from large to small, and the sorted result is shown in Figure 6(b).Figure 6
In Fourier dictionary, (a) sparse coefficients of pure simulation signal; (b) attenuation curve of (a).
(a)(b)The attenuation curve of Figure6(b) is fitted by using the least-square method, through the curve fitting toolbox of MATLAB, and the formula obtained is fx=0.2e−0.059x+0.03e−0.0018x. The actual attenuation curve and the fitting curve are as presented in Figure 7, where SSE is 0.011 and R-square is 0.984, which are both close to the ideal value. Furthermore, from the detailed area, it is obvious that the sparse coefficients decrease to zero. So it can be considered that the Fourier coefficients decay to zero exponentially, indicating that the vibration signal of rolling bearing is almost sparse in the Fourier domain.Figure 7
The fitting result of coefficients attenuation curve in Figure6(b): SSE = 0.011; R-square = 0.984.The existence of noise in the vibration signal will increase the difficulty of sparse representation in the Fourier domain. Sparse coefficients of the noisy signal and the sequencing result are shown in Figures8(a) and 8(b), respectively. The fitting result of the attenuation curve is shown in Figure 9, and the formula of the fitted curve is fx=0.065e−0.00047x−5.5×10−12e−0.0052x. Due to the influence of noise, SSE is as high as 0.18 and R-square is only 0.852, both of which deviate slightly from the ideal value. What is more, the attenuation curvature of sparse coefficients of the noisy signal is increasing, not decreasing as exponential attenuation. Thus, sparse coefficients in Figure 8(b) cannot be concluded to attenuate to zero exponentially, which means the polluted signal has inadequate sparsity affected by full-band noise.Figure 8
In Fourier dictionary, (a) sparse coefficients of the noisy signal; (b) attenuation curve of (a).
(a)(b)Figure 9
The fitting result of coefficients attenuation curve in Figure8(b): SSE = 0.18; R-square = 0.852.Finally, the sparsity of mode 4 in Table2 is analyzed in the Fourier domain. Sparse coefficients of the denoised signal and the sorting result are shown in Figures 10(a) and 10(b), respectively. The attenuation curve of Figure 10(b) is further fitted, and the formula obtained is fx=0.16e−0.091x+0.056e−0.0031x. In Figure 11, SSE is 0.007, R-square is 0.994, and the actual attenuation curve almost overlaps the fitted curve. It can be inferred that the sparse coefficients of the denoised signal can be decreased exponentially to zero under the Fourier transform. Compared with the noisy signal, the sparsity of the denoised signal is greatly enhanced.Figure 10
In Fourier dictionary, (a) sparse coefficients of the mode 3 and (b) attenuation curve of (a).
(a)(b)Figure 11
The fitting result of coefficient attenuation curve in Figure10(b): SSE = 0.007; R-square = 0.994.


## 4.3.2. High Sensing Efficiency of MCM

Gaussian random matrix (GRM) is frequently used as a measurement matrix of CS for a typical one-dimensional vibration signal [18–21, 25]. The reason is just simple and convenient to construct, and most importantly, the measurement matrix constructed by the basis vector obeying the Gauss distribution can be roughly guaranteed to be irrelevant to any sparse basis.In this study, a new compressive sensing measurement matrix is constructed in the form of the circulant matrix by using the VMD modes. Its sensing efficiency is much higher than GRM, which means the extraction accuracy of the MCM can be achieved by fewer measurements than GRM. The following parts will demonstrate it.This verification firstly used the modes in Table2 to construct the MCM and then applied it and GRM, respectively, to compress the denoised vibration signal. Fault feature extraction is considered successful when the fault feature order appears in the final envelope spectrum.Figure12(a) and 12(b) are the final envelope order spectra obtained by using MCM and GRM, respectively, when the number of measurements is 600. The dominant order in Figure 12(a) is consistent with the inner simulated fault order O = 14.9, representing a 100% success rate of fault feature extraction, while Figure 12(b) fails; that is, the success rate is 0. According to (24), the SNR2 of Figure 12(a) and 12(b) are 33.4 dB and −48.7 dB, respectively, and bearing fault can be directly identified with the former.Figure 12
The envelope order spectra obtained by using (a) MCM: the spectrum is clear, and SNR2 is 33.4 dB; (b) GRM: the spectrum is clear, but SNR2 is only −44.7 dB.
(a)(b)Moreover, the change in success rate is shown in Figure13 as the number of measurements increases from 20 to 1000. From Figure 13, when the number of measurements is only 60, the MCM can successfully extract bearing fault features, and the results are stable and consistent within a varied number of measurements. By contrast, the GRM has obvious randomness and poor robustness although it can occasionally extract fault features under certain measurements. Therefore, sensing efficiency of the MCM is higher than that of the GRM, which means that the beneficial information of the signal can be measured more effectively.Figure 13
Comparison of detection effect of two measuring matrices.


## 4.3.3. Validation of NEA Method

Section4.3.2 has verified that MCM has higher sensing efficiency than GRM. Combined with CS, MCM-CS can effectively extract bearing fault features under a fewer number of measurements. To further highlight the significance of MCM-CS for bearing fault compressive detection, the Narrowband Envelope Analysis (NEA) is selected for comparison.Before analyzing the bearing fault signal to NEA, the most fundamental issue to be solved first is the setting of the center frequency and bandwidth of the band-pass filter. However, it fails to extract effective features based on experience frequently. In recent years, SK has been extensively used in the selection of demodulation band parameters and gained good reputation [35]. Thus, SK is used to design the band-pass filter in the NEA method in this section. Besides, considering variable speed conditions, NEA is carried out as follows:(1)
Even angle resampling in the speed domain, order tracking is applied to the original signal, achieving a stationary signal in the angle domain instead of a nonstationary signal in the time domain.(2)
The optimal envelope center frequency and bandwidth are determined by SK, based on which a band-pass filter is designed to filter the resampled signal.(3)
The envelope demodulation method based on the Hilbert transform is applied to achieve the envelope order spectrum.Figure14 shows the analysis results for the fault simulation signal shown in Figure 3(a) using the NEA based on SK and order tracking. Figure 14(a) is the kurtogram. The SK shows that the center frequency of the optimum demodulation band is 2933 Hz, which is extremely near to the resonance frequency of the simulated signal. This proves that SK is highly significant. Figure 14(b) shows the obscure envelope order spectrum, where the amplitude of 2.099 order is the highest. An order of 14.79, corresponding to the inner fault, can be seen but it is not the dominant component in this spectrum.Figure 14
The analysis results using NEA: (a) kurtogram; (b) envelope order spectrum: the spectrum is obscure, and SNR2 is only −9.07 dB.
(a)(b)Besides, according to (24), the SNR2 of the filtered signal is −9.07 dB, which is far from the recovered signal through MCM-CS. As observed in the above, although NEA is indeed simpler than MCM-CS, it is not only unable to achieve fault compressive detection but also inferior in weak fault feature extraction.The facts of Section4.3 prove that the VMD mode of fault vibration signals of rolling bearings can be sparsely represented in the Fourier dictionary. What is more, the MCM-enhanced compressive sensing is endowed with higher sensing efficiency. These are the necessary conditions for CS to extract fault feature of rolling bearings in Section 5.


## 5. Experiments

To further verify the effectiveness of the weak fault feature extraction method based on the sensing-enhanced CS, the weak fault signals of the rolling bearings under variable speed conditions are collected in the test rig described below.


### 5.1. Experimental Setup

To simulate the weak fault of rolling elements to a large extent, the dual-rotor test rig is selected, which is designed and built to simulate the working state of the intershaft bearing on the actual double rotor equipment. As is illustrated in Figure15, it consists of two driving motors, #1 intershaft bearing, #2 bearing, #2 bearing housing, and dual rotors. In particular, #2 bearings are fitted with thin-walled supports.Figure 15
Dual-rotor vibration test rig.In actual dual-rotor equipment, since the inner and outer rings of intershaft bearings are rotating, there is no fixed bearing housing for the installation of sensors, which must require indirect measurement of vibration. The indirectly measured vibration response needs to be transmitted to the adjacent bearing via the shaft and then to sensor position through the elastic support and complicated thin-wall paths so that the intershaft bearing signal is buried in the strong vibration and interference.Similar to the actual dual-rotor equipment, #1 intershaft bearing, in the dual-rotor test rig, has not directly linked bearing housing. Thus, although the inner and outer races of #1 intershaft bearing are machined defectively, the accelerometer is installed at #2 bearing housing. Figure16 is the transfer path of vibration of #1 intershaft bearing in the dual-rotor test rig, it can be seen that the test rig reproduces the running status of an actual intershaft bearing: the vibration response of #1 intershaft bearing needs to be transferred to adjacent bearing via the inner rotor; then, via elastic supports and thin-walled structure, the response is transmitted to the #2 bearing housing. The complicated transfer path indicates that the detected bearing signal is severely attenuated during the transfer. Therefore, compared with traditional rolling bearings used with a single rotor and fixed housing, the fault signal of #1 intershaft bearing collected in this experiment is relatively weak.Figure 16
Transfer path of vibration signal of #1 intershaft bearing.The faults of inner and outer races of the intershaft bearing were simulated as follows. As for the outer race defect, a 1 mm deep and 1 mm wide groove was machined across the axis direction on the inner surface of the outer race. The outer surface of the inner race was also machined with the same size defects to simulate inner race fault.The vibration signal is acquired by the BK4519 acceleration sensor, a key phase sensor, charge amplifier, and LMS SCADAS [36] data acquisition system. During the experiment, only the drive motor of the inner race rotor of 1 # intermediate bearing is started, while the outer ring rotor is fixed.The 1# intershaft bearing fault vibration data used in this study are divided into the following two groups: (1) outer race fault signal where the inner race speed is increased from 625 r/min to 766 r/min; the sampling frequency is 25.6 kHz, and the number of vibration signals sampled for analysis is 102400; (2) inner race fault signal where rotational speed of inner race increased from 738 r/min to 825 r/min; sampling frequency is 102.4 kHz, and the amount of sample is 204800. The vibration waveforms and speed trends are shown in Figures17 and 18.Figure 17
Fault in the outer race: (a) waveform of vibration; (b) speed trend.
(a)(b)Figure 18
Fault in the inner race: (a) waveform of vibration; (b) speed trend.
(a)(b)


### 5.2. Fault Characteristic Frequencies of Rolling Bearing

In the experiment, geometrical parameters of the faulty bearing are listed in Table3, and the fault characteristic orders are shown in Table 4.Table 3
Geometrical parameters of 1# intershaft bearing.
Pitch diameter D (mm)Ball diameter d (mm)Contact angleα (°)Number of rolling elementsz1258034Table 4
Fault characteristic orders of 1# intershaft bearing.
ComponentOuter raceInner raceRolling elementsCageOrder15.918.17.80.5


### 5.3. Compressive Extraction of Fault Feature by MCM-CS

#### 5.3.1. Outer Race Fault Detection

The proposed method is first utilized to extract the faulty feature of 1# intershaft bearings with an outer race fault. As presented in Figure17, the original vibration signal has intense time-domain nonstationarity and the number of vibration signals for fault extraction is a bit larger, which increases the difficulty of fault feature extraction based on CS. Firstly, the original vibration signal is order tracked to obtain the stationary signal in the angle domain. Compared with the original signal, the stationarity of the angular signal in Figure 19 is greatly enhanced. Although there are still more than 23000, the sample is compressed by 77.6% while retaining the faulty feature of the bearing, which eases the pressure of further compression by CS.Figure 19
The angular signal of outer race fault.PSO is utilized to analyze the outer race fault signal in Figure19. The maximum value of local energy ratio changes with the evolution algebra as shown in Figure 20. The main parameters’ setting and optimization results are presented in Table 5. The optimal combination of [α, K] obtained by optimization is [65, 3].Figure 20
The trend of local energy ratio with evolution algebra.Table 5
Parameters setting and optimization results of PSO.
ParametersαKRange of values[0, 2000][1, 10]Maximum evolution50Maximum population10Local search speed1Global search speed1.2Optimization results653The angular signal is decomposed by VMD, where the penalty factor and decomposition numbers are artificially set to be 65 and 3, respectively. Table6 shows three modes obtained by VMD decomposition, the fitness value of each mode. It indicates that mode 3 is the best mode corresponding to the maximum energy ratio.Table 6
Fitness values obtained by analyzing the outer race fault signal with VMD.
Modes obtained by VMD and their order spectraFitness (%)5.263.26.28Therefore, mode 3, as depicted in Figure21, will replace the angular signal shown in Figure 19 for compression detection of fault feature.Figure 21
The denoised signal of outer race fault.Modes obtained by VMD are continued to be used to construct the MCM, in the detection of the outer-race faulty bearing based on CS. The study selects the number of measurements to be 1500. Next, the amount of the analyzed signal is compressed to 1500 sampling points assisted by MCM and GRM, respectively.Based on the constructed MCM, the fault feature is detected by orthogonal matching pursuit with the help of the Fourier dictionary and setting sparsityK = 4. As provided in Table 4, the theoretical fault feature order of the outer race fault is 15.9. The detected result is illustrated in Figure 22(a), where 239.7 and 255.6 are the sparse coefficients under the Fourier dictionary, and the difference between them is 15.9.Figure 22
The extracted result of outer race fault using CS-based on MCM: (a) the detected sparse coefficients; (b) the visualized envelope order spectrum.
(a)(b)Furthermore, the Fourier dictionary is used to transform the detected result in Figure22(a) into the recovered signal, and envelope demodulation is carried out to obtain the order spectrum which can directly recognize the fault feature, as presented in Figure 22(b). Thus, since the extracted result is equal to the theoretical one, a conclusion can be drawn that there is a fault in the outer race of 1# intershaft bearing.


#### 5.3.2. Inner Race Fault Detection

To fully validate the effectiveness of the proposed fault feature extraction strategy, 1# intershaft bearing with a faulty inner race was employed. Firstly, the raw vibration signal of inner race fault in Figure18 is processed by order tracking to obtain the angular vibration signal, the waveform of which is as shown in Figure 23.Figure 23
The angular signal of inner race fault.PSO is utilized to analyze the inner race fault signal in Figure24. A fitness curve can be obtained, as shown in Figure 24. The main parameters’ setting and optimization results are shown in Table 7. The optimal combination of [α, K] obtained by optimization is [109, 5].Figure 24
The trend of local energy ratio with evolution algebra.Table 7
Parameters setting and optimization results of PSO.
ParametersαKRange of values[0, 2000][1, 10]Maximum evolution50Maximum population10Local search speed1Global search speed1.2Optimization results1095Five modes are obtained by VMD decomposition with penalty factor 109 and the decomposition numbers 5, which are prepared for denoising and constructing the MCM. Among the modes, mode 2 is the best component corresponding to the maximum energy ratio, so mode 2 is the main mode (Table8. The waveform of mode 2 is shown in Figure 25.Table 8
Fitness values obtained by analyzing inner race fault signal with VMD.
Modes obtained by VMD and their order spectraFitness (%)0.534.291.351.370.77Figure 25
The denoised signal of inner race fault.Next, the MCM and GRM are used as a measurement matrix to acquire low-dimensional signals with the measurement numbers set to 1500, respectively. When the sparsity is determined as 4, the coefficients corresponding to the faulty feature can be detected from the compressed signal by orthogonal matching pursuit with the sensing matrix constructed by MCM and Fourier dictionary. Next, the detected result, as shown in Figure26(a), can be recovered by the Fourier dictionary. Finally, the recovered signal can be demodulated by envelope analysis to obtain the envelope order spectrum, as presented in Figure 26(b). The previous contains the modulated frequency order 162.8 and sideband 18.1 which are almost equal to the theoretical one of the inner fault feature. Moreover, the faulty feature order 18.09 and its twice value are extracted directly in the envelope order spectrum. Thus, the judgment that there is an inner-race fault in 1# intershaft bearing can be concluded.Figure 26
The extracted result of inner race fault using CS-based on MCM: (a) the detected sparse coefficients and (b) the visualized envelope order spectrum.
(a)(b)


### 5.4. Effect Comparison of Different Methods

#### 5.4.1. Compressive Extraction of Fault Feature by GRM-CS

To further highlight the high sensing efficiency of MCM and the advantages of MCM-CS in signal compression and feature extraction, the GRM-based CS (GRM-CS) is selected as a comparison. The following steps can be used:(1)
Signal preprocessing and sparsity enhancement: the same is as in Sections3.1 and 3.2.(2)
Compressive sensing: the GRM is utilized to reduce the dimensions of the denoised signal, and the Fourier dictionary and GRM make up the sensing matrix; finally, the sparse vector is searched by orthogonal matching pursuit.(3)
Fault feature extraction: the signal recovery can be achieved using the Fourier dictionary and the sparse vector. Because of amplitude modulation of the recovered signal, the envelope demodulation is applied to achieve the envelope order spectrum.Thus, the envelope order spectrum is obtained, as shown in Figures27 and 28.Figure 27
The envelope order spectrum of outer race fault using GRM-CS.Figure 28
The analysis results for fault data of outer race by using NEA based on SK: (a) kurtogram; (b) envelope order spectrum.Figure27 shows the envelope order spectrum of the fault signal of the intershaft bearing’s outer race, in which the dominant order is 58.18. An order of 15.9, corresponding to a fault in the outer race, can be seen but it is not the dominant component.Figure29 shows the envelope order spectrum of the fault signal of the intershaft bearing’s inner race. However, the order component 18.09 corresponding to the inner-race fault of 1# intershaft bearing cannot appear in Figure 29, which means the fault feature cannot be extracted through the GRM-based CS.Figure 29
The envelope order spectrum of inner race fault using GRM-CS.
(a)(b)The analysis of experimental data proves that the proposed MCM-based CS method has higher sensing efficiency than GRM-based CS again. Under the same measurement number, the former can successfully extract the weak fault feature of rolling bearings.


#### 5.4.2. Fault Feature Extraction by NEA Based on SK

To assess the value of the proposed method properly and show that its novelty can be compensated by performance, following the steps in Section4.3.3, we use NEA based on SK to analyze experimental signals for comparison. However, the following analysis shows that it is hard to achieve significant effects with the SK-based NEA.Figure28 shows the analysis results for the fault signal of the 1# intershaft bearing’s outer race using NEA based on SK. Figure 28(a) is the kurtogram, where the center frequency of the optimum band-pass filter is 2000 Hz and the bandwidth is 800 Hz. That is, the optimum demodulation band is [1600, 2400] Hz. Figure 28(b) shows the envelope order spectrum. In this spectrum, an order of 2.7 has the highest amplitude. The outer race fault feature order 15.9 is fully buried in noise.Similarly, Figure30 shows the analysis results for the fault signal of the inner race. Figure 30(a) shows that the optimal demodulation frequency band is [0, 25600] Hz. Figure 30(b) is the envelope order spectrum; an order of 1 and its harmonic are the dominant components, in which three times harmonic is particularly prominent. An order of 18.1, corresponding to the fault of the bearing inner race, is not effectively extracted.Figure 30
The analysis results for fault data of inner race by using NEA based on SK: (a) kurtogram; (b) envelope order spectrum.
(a)(b)It can be seen that the order of prominence in Figures28 and 30 is the result of other impacts in the rotor system. Owing to its existence, NEA based on SK fails to extract the features generated by bearing fault impact. Combined with the analysis of Section 4.3.3, Section 5.3, and Section 5.4.2, it is indicated that SK-based NEA fails to extract the weak fault features of bearings in the presence of random noise and impact noise, which exactly highlights the unique advantages of MCM-CS in face of such difficulties.


## 5.1. Experimental Setup

To simulate the weak fault of rolling elements to a large extent, the dual-rotor test rig is selected, which is designed and built to simulate the working state of the intershaft bearing on the actual double rotor equipment. As is illustrated in Figure15, it consists of two driving motors, #1 intershaft bearing, #2 bearing, #2 bearing housing, and dual rotors. In particular, #2 bearings are fitted with thin-walled supports.Figure 15
Dual-rotor vibration test rig.In actual dual-rotor equipment, since the inner and outer rings of intershaft bearings are rotating, there is no fixed bearing housing for the installation of sensors, which must require indirect measurement of vibration. The indirectly measured vibration response needs to be transmitted to the adjacent bearing via the shaft and then to sensor position through the elastic support and complicated thin-wall paths so that the intershaft bearing signal is buried in the strong vibration and interference.Similar to the actual dual-rotor equipment, #1 intershaft bearing, in the dual-rotor test rig, has not directly linked bearing housing. Thus, although the inner and outer races of #1 intershaft bearing are machined defectively, the accelerometer is installed at #2 bearing housing. Figure16 is the transfer path of vibration of #1 intershaft bearing in the dual-rotor test rig, it can be seen that the test rig reproduces the running status of an actual intershaft bearing: the vibration response of #1 intershaft bearing needs to be transferred to adjacent bearing via the inner rotor; then, via elastic supports and thin-walled structure, the response is transmitted to the #2 bearing housing. The complicated transfer path indicates that the detected bearing signal is severely attenuated during the transfer. Therefore, compared with traditional rolling bearings used with a single rotor and fixed housing, the fault signal of #1 intershaft bearing collected in this experiment is relatively weak.Figure 16
Transfer path of vibration signal of #1 intershaft bearing.The faults of inner and outer races of the intershaft bearing were simulated as follows. As for the outer race defect, a 1 mm deep and 1 mm wide groove was machined across the axis direction on the inner surface of the outer race. The outer surface of the inner race was also machined with the same size defects to simulate inner race fault.The vibration signal is acquired by the BK4519 acceleration sensor, a key phase sensor, charge amplifier, and LMS SCADAS [36] data acquisition system. During the experiment, only the drive motor of the inner race rotor of 1 # intermediate bearing is started, while the outer ring rotor is fixed.The 1# intershaft bearing fault vibration data used in this study are divided into the following two groups: (1) outer race fault signal where the inner race speed is increased from 625 r/min to 766 r/min; the sampling frequency is 25.6 kHz, and the number of vibration signals sampled for analysis is 102400; (2) inner race fault signal where rotational speed of inner race increased from 738 r/min to 825 r/min; sampling frequency is 102.4 kHz, and the amount of sample is 204800. The vibration waveforms and speed trends are shown in Figures17 and 18.Figure 17
Fault in the outer race: (a) waveform of vibration; (b) speed trend.
(a)(b)Figure 18
Fault in the inner race: (a) waveform of vibration; (b) speed trend.
(a)(b)


## 5.2. Fault Characteristic Frequencies of Rolling Bearing

In the experiment, geometrical parameters of the faulty bearing are listed in Table3, and the fault characteristic orders are shown in Table 4.Table 3
Geometrical parameters of 1# intershaft bearing.
Pitch diameter D (mm)Ball diameter d (mm)Contact angleα (°)Number of rolling elementsz1258034Table 4
Fault characteristic orders of 1# intershaft bearing.
ComponentOuter raceInner raceRolling elementsCageOrder15.918.17.80.5


## 5.3. Compressive Extraction of Fault Feature by MCM-CS

### 5.3.1. Outer Race Fault Detection

The proposed method is first utilized to extract the faulty feature of 1# intershaft bearings with an outer race fault. As presented in Figure17, the original vibration signal has intense time-domain nonstationarity and the number of vibration signals for fault extraction is a bit larger, which increases the difficulty of fault feature extraction based on CS. Firstly, the original vibration signal is order tracked to obtain the stationary signal in the angle domain. Compared with the original signal, the stationarity of the angular signal in Figure 19 is greatly enhanced. Although there are still more than 23000, the sample is compressed by 77.6% while retaining the faulty feature of the bearing, which eases the pressure of further compression by CS.Figure 19
The angular signal of outer race fault.PSO is utilized to analyze the outer race fault signal in Figure19. The maximum value of local energy ratio changes with the evolution algebra as shown in Figure 20. The main parameters’ setting and optimization results are presented in Table 5. The optimal combination of [α, K] obtained by optimization is [65, 3].Figure 20
The trend of local energy ratio with evolution algebra.Table 5
Parameters setting and optimization results of PSO.
ParametersαKRange of values[0, 2000][1, 10]Maximum evolution50Maximum population10Local search speed1Global search speed1.2Optimization results653The angular signal is decomposed by VMD, where the penalty factor and decomposition numbers are artificially set to be 65 and 3, respectively. Table6 shows three modes obtained by VMD decomposition, the fitness value of each mode. It indicates that mode 3 is the best mode corresponding to the maximum energy ratio.Table 6
Fitness values obtained by analyzing the outer race fault signal with VMD.
Modes obtained by VMD and their order spectraFitness (%)5.263.26.28Therefore, mode 3, as depicted in Figure21, will replace the angular signal shown in Figure 19 for compression detection of fault feature.Figure 21
The denoised signal of outer race fault.Modes obtained by VMD are continued to be used to construct the MCM, in the detection of the outer-race faulty bearing based on CS. The study selects the number of measurements to be 1500. Next, the amount of the analyzed signal is compressed to 1500 sampling points assisted by MCM and GRM, respectively.Based on the constructed MCM, the fault feature is detected by orthogonal matching pursuit with the help of the Fourier dictionary and setting sparsityK = 4. As provided in Table 4, the theoretical fault feature order of the outer race fault is 15.9. The detected result is illustrated in Figure 22(a), where 239.7 and 255.6 are the sparse coefficients under the Fourier dictionary, and the difference between them is 15.9.Figure 22
The extracted result of outer race fault using CS-based on MCM: (a) the detected sparse coefficients; (b) the visualized envelope order spectrum.
(a)(b)Furthermore, the Fourier dictionary is used to transform the detected result in Figure22(a) into the recovered signal, and envelope demodulation is carried out to obtain the order spectrum which can directly recognize the fault feature, as presented in Figure 22(b). Thus, since the extracted result is equal to the theoretical one, a conclusion can be drawn that there is a fault in the outer race of 1# intershaft bearing.


### 5.3.2. Inner Race Fault Detection

To fully validate the effectiveness of the proposed fault feature extraction strategy, 1# intershaft bearing with a faulty inner race was employed. Firstly, the raw vibration signal of inner race fault in Figure18 is processed by order tracking to obtain the angular vibration signal, the waveform of which is as shown in Figure 23.Figure 23
The angular signal of inner race fault.PSO is utilized to analyze the inner race fault signal in Figure24. A fitness curve can be obtained, as shown in Figure 24. The main parameters’ setting and optimization results are shown in Table 7. The optimal combination of [α, K] obtained by optimization is [109, 5].Figure 24
The trend of local energy ratio with evolution algebra.Table 7
Parameters setting and optimization results of PSO.
ParametersαKRange of values[0, 2000][1, 10]Maximum evolution50Maximum population10Local search speed1Global search speed1.2Optimization results1095Five modes are obtained by VMD decomposition with penalty factor 109 and the decomposition numbers 5, which are prepared for denoising and constructing the MCM. Among the modes, mode 2 is the best component corresponding to the maximum energy ratio, so mode 2 is the main mode (Table8. The waveform of mode 2 is shown in Figure 25.Table 8
Fitness values obtained by analyzing inner race fault signal with VMD.
Modes obtained by VMD and their order spectraFitness (%)0.534.291.351.370.77Figure 25
The denoised signal of inner race fault.Next, the MCM and GRM are used as a measurement matrix to acquire low-dimensional signals with the measurement numbers set to 1500, respectively. When the sparsity is determined as 4, the coefficients corresponding to the faulty feature can be detected from the compressed signal by orthogonal matching pursuit with the sensing matrix constructed by MCM and Fourier dictionary. Next, the detected result, as shown in Figure26(a), can be recovered by the Fourier dictionary. Finally, the recovered signal can be demodulated by envelope analysis to obtain the envelope order spectrum, as presented in Figure 26(b). The previous contains the modulated frequency order 162.8 and sideband 18.1 which are almost equal to the theoretical one of the inner fault feature. Moreover, the faulty feature order 18.09 and its twice value are extracted directly in the envelope order spectrum. Thus, the judgment that there is an inner-race fault in 1# intershaft bearing can be concluded.Figure 26
The extracted result of inner race fault using CS-based on MCM: (a) the detected sparse coefficients and (b) the visualized envelope order spectrum.
(a)(b)


## 5.3.1. Outer Race Fault Detection

The proposed method is first utilized to extract the faulty feature of 1# intershaft bearings with an outer race fault. As presented in Figure17, the original vibration signal has intense time-domain nonstationarity and the number of vibration signals for fault extraction is a bit larger, which increases the difficulty of fault feature extraction based on CS. Firstly, the original vibration signal is order tracked to obtain the stationary signal in the angle domain. Compared with the original signal, the stationarity of the angular signal in Figure 19 is greatly enhanced. Although there are still more than 23000, the sample is compressed by 77.6% while retaining the faulty feature of the bearing, which eases the pressure of further compression by CS.Figure 19
The angular signal of outer race fault.PSO is utilized to analyze the outer race fault signal in Figure19. The maximum value of local energy ratio changes with the evolution algebra as shown in Figure 20. The main parameters’ setting and optimization results are presented in Table 5. The optimal combination of [α, K] obtained by optimization is [65, 3].Figure 20
The trend of local energy ratio with evolution algebra.Table 5
Parameters setting and optimization results of PSO.
ParametersαKRange of values[0, 2000][1, 10]Maximum evolution50Maximum population10Local search speed1Global search speed1.2Optimization results653The angular signal is decomposed by VMD, where the penalty factor and decomposition numbers are artificially set to be 65 and 3, respectively. Table6 shows three modes obtained by VMD decomposition, the fitness value of each mode. It indicates that mode 3 is the best mode corresponding to the maximum energy ratio.Table 6
Fitness values obtained by analyzing the outer race fault signal with VMD.
Modes obtained by VMD and their order spectraFitness (%)5.263.26.28Therefore, mode 3, as depicted in Figure21, will replace the angular signal shown in Figure 19 for compression detection of fault feature.Figure 21
The denoised signal of outer race fault.Modes obtained by VMD are continued to be used to construct the MCM, in the detection of the outer-race faulty bearing based on CS. The study selects the number of measurements to be 1500. Next, the amount of the analyzed signal is compressed to 1500 sampling points assisted by MCM and GRM, respectively.Based on the constructed MCM, the fault feature is detected by orthogonal matching pursuit with the help of the Fourier dictionary and setting sparsityK = 4. As provided in Table 4, the theoretical fault feature order of the outer race fault is 15.9. The detected result is illustrated in Figure 22(a), where 239.7 and 255.6 are the sparse coefficients under the Fourier dictionary, and the difference between them is 15.9.Figure 22
The extracted result of outer race fault using CS-based on MCM: (a) the detected sparse coefficients; (b) the visualized envelope order spectrum.
(a)(b)Furthermore, the Fourier dictionary is used to transform the detected result in Figure22(a) into the recovered signal, and envelope demodulation is carried out to obtain the order spectrum which can directly recognize the fault feature, as presented in Figure 22(b). Thus, since the extracted result is equal to the theoretical one, a conclusion can be drawn that there is a fault in the outer race of 1# intershaft bearing.


## 5.3.2. Inner Race Fault Detection

To fully validate the effectiveness of the proposed fault feature extraction strategy, 1# intershaft bearing with a faulty inner race was employed. Firstly, the raw vibration signal of inner race fault in Figure18 is processed by order tracking to obtain the angular vibration signal, the waveform of which is as shown in Figure 23.Figure 23
The angular signal of inner race fault.PSO is utilized to analyze the inner race fault signal in Figure24. A fitness curve can be obtained, as shown in Figure 24. The main parameters’ setting and optimization results are shown in Table 7. The optimal combination of [α, K] obtained by optimization is [109, 5].Figure 24
The trend of local energy ratio with evolution algebra.Table 7
Parameters setting and optimization results of PSO.
ParametersαKRange of values[0, 2000][1, 10]Maximum evolution50Maximum population10Local search speed1Global search speed1.2Optimization results1095Five modes are obtained by VMD decomposition with penalty factor 109 and the decomposition numbers 5, which are prepared for denoising and constructing the MCM. Among the modes, mode 2 is the best component corresponding to the maximum energy ratio, so mode 2 is the main mode (Table8. The waveform of mode 2 is shown in Figure 25.Table 8
Fitness values obtained by analyzing inner race fault signal with VMD.
Modes obtained by VMD and their order spectraFitness (%)0.534.291.351.370.77Figure 25
The denoised signal of inner race fault.Next, the MCM and GRM are used as a measurement matrix to acquire low-dimensional signals with the measurement numbers set to 1500, respectively. When the sparsity is determined as 4, the coefficients corresponding to the faulty feature can be detected from the compressed signal by orthogonal matching pursuit with the sensing matrix constructed by MCM and Fourier dictionary. Next, the detected result, as shown in Figure26(a), can be recovered by the Fourier dictionary. Finally, the recovered signal can be demodulated by envelope analysis to obtain the envelope order spectrum, as presented in Figure 26(b). The previous contains the modulated frequency order 162.8 and sideband 18.1 which are almost equal to the theoretical one of the inner fault feature. Moreover, the faulty feature order 18.09 and its twice value are extracted directly in the envelope order spectrum. Thus, the judgment that there is an inner-race fault in 1# intershaft bearing can be concluded.Figure 26
The extracted result of inner race fault using CS-based on MCM: (a) the detected sparse coefficients and (b) the visualized envelope order spectrum.
(a)(b)


## 5.4. Effect Comparison of Different Methods

### 5.4.1. Compressive Extraction of Fault Feature by GRM-CS

To further highlight the high sensing efficiency of MCM and the advantages of MCM-CS in signal compression and feature extraction, the GRM-based CS (GRM-CS) is selected as a comparison. The following steps can be used:(1)
Signal preprocessing and sparsity enhancement: the same is as in Sections3.1 and 3.2.(2)
Compressive sensing: the GRM is utilized to reduce the dimensions of the denoised signal, and the Fourier dictionary and GRM make up the sensing matrix; finally, the sparse vector is searched by orthogonal matching pursuit.(3)
Fault feature extraction: the signal recovery can be achieved using the Fourier dictionary and the sparse vector. Because of amplitude modulation of the recovered signal, the envelope demodulation is applied to achieve the envelope order spectrum.Thus, the envelope order spectrum is obtained, as shown in Figures27 and 28.Figure 27
The envelope order spectrum of outer race fault using GRM-CS.Figure 28
The analysis results for fault data of outer race by using NEA based on SK: (a) kurtogram; (b) envelope order spectrum.Figure27 shows the envelope order spectrum of the fault signal of the intershaft bearing’s outer race, in which the dominant order is 58.18. An order of 15.9, corresponding to a fault in the outer race, can be seen but it is not the dominant component.Figure29 shows the envelope order spectrum of the fault signal of the intershaft bearing’s inner race. However, the order component 18.09 corresponding to the inner-race fault of 1# intershaft bearing cannot appear in Figure 29, which means the fault feature cannot be extracted through the GRM-based CS.Figure 29
The envelope order spectrum of inner race fault using GRM-CS.
(a)(b)The analysis of experimental data proves that the proposed MCM-based CS method has higher sensing efficiency than GRM-based CS again. Under the same measurement number, the former can successfully extract the weak fault feature of rolling bearings.


### 5.4.2. Fault Feature Extraction by NEA Based on SK

To assess the value of the proposed method properly and show that its novelty can be compensated by performance, following the steps in Section4.3.3, we use NEA based on SK to analyze experimental signals for comparison. However, the following analysis shows that it is hard to achieve significant effects with the SK-based NEA.Figure28 shows the analysis results for the fault signal of the 1# intershaft bearing’s outer race using NEA based on SK. Figure 28(a) is the kurtogram, where the center frequency of the optimum band-pass filter is 2000 Hz and the bandwidth is 800 Hz. That is, the optimum demodulation band is [1600, 2400] Hz. Figure 28(b) shows the envelope order spectrum. In this spectrum, an order of 2.7 has the highest amplitude. The outer race fault feature order 15.9 is fully buried in noise.Similarly, Figure30 shows the analysis results for the fault signal of the inner race. Figure 30(a) shows that the optimal demodulation frequency band is [0, 25600] Hz. Figure 30(b) is the envelope order spectrum; an order of 1 and its harmonic are the dominant components, in which three times harmonic is particularly prominent. An order of 18.1, corresponding to the fault of the bearing inner race, is not effectively extracted.Figure 30
The analysis results for fault data of inner race by using NEA based on SK: (a) kurtogram; (b) envelope order spectrum.
(a)(b)It can be seen that the order of prominence in Figures28 and 30 is the result of other impacts in the rotor system. Owing to its existence, NEA based on SK fails to extract the features generated by bearing fault impact. Combined with the analysis of Section 4.3.3, Section 5.3, and Section 5.4.2, it is indicated that SK-based NEA fails to extract the weak fault features of bearings in the presence of random noise and impact noise, which exactly highlights the unique advantages of MCM-CS in face of such difficulties.


## 5.4.1. Compressive Extraction of Fault Feature by GRM-CS

To further highlight the high sensing efficiency of MCM and the advantages of MCM-CS in signal compression and feature extraction, the GRM-based CS (GRM-CS) is selected as a comparison. The following steps can be used:(1)
Signal preprocessing and sparsity enhancement: the same is as in Sections3.1 and 3.2.(2)
Compressive sensing: the GRM is utilized to reduce the dimensions of the denoised signal, and the Fourier dictionary and GRM make up the sensing matrix; finally, the sparse vector is searched by orthogonal matching pursuit.(3)
Fault feature extraction: the signal recovery can be achieved using the Fourier dictionary and the sparse vector. Because of amplitude modulation of the recovered signal, the envelope demodulation is applied to achieve the envelope order spectrum.Thus, the envelope order spectrum is obtained, as shown in Figures27 and 28.Figure 27
The envelope order spectrum of outer race fault using GRM-CS.Figure 28
The analysis results for fault data of outer race by using NEA based on SK: (a) kurtogram; (b) envelope order spectrum.Figure27 shows the envelope order spectrum of the fault signal of the intershaft bearing’s outer race, in which the dominant order is 58.18. An order of 15.9, corresponding to a fault in the outer race, can be seen but it is not the dominant component.Figure29 shows the envelope order spectrum of the fault signal of the intershaft bearing’s inner race. However, the order component 18.09 corresponding to the inner-race fault of 1# intershaft bearing cannot appear in Figure 29, which means the fault feature cannot be extracted through the GRM-based CS.Figure 29
The envelope order spectrum of inner race fault using GRM-CS.
(a)(b)The analysis of experimental data proves that the proposed MCM-based CS method has higher sensing efficiency than GRM-based CS again. Under the same measurement number, the former can successfully extract the weak fault feature of rolling bearings.


## 5.4.2. Fault Feature Extraction by NEA Based on SK

To assess the value of the proposed method properly and show that its novelty can be compensated by performance, following the steps in Section4.3.3, we use NEA based on SK to analyze experimental signals for comparison. However, the following analysis shows that it is hard to achieve significant effects with the SK-based NEA.Figure28 shows the analysis results for the fault signal of the 1# intershaft bearing’s outer race using NEA based on SK. Figure 28(a) is the kurtogram, where the center frequency of the optimum band-pass filter is 2000 Hz and the bandwidth is 800 Hz. That is, the optimum demodulation band is [1600, 2400] Hz. Figure 28(b) shows the envelope order spectrum. In this spectrum, an order of 2.7 has the highest amplitude. The outer race fault feature order 15.9 is fully buried in noise.Similarly, Figure30 shows the analysis results for the fault signal of the inner race. Figure 30(a) shows that the optimal demodulation frequency band is [0, 25600] Hz. Figure 30(b) is the envelope order spectrum; an order of 1 and its harmonic are the dominant components, in which three times harmonic is particularly prominent. An order of 18.1, corresponding to the fault of the bearing inner race, is not effectively extracted.Figure 30
The analysis results for fault data of inner race by using NEA based on SK: (a) kurtogram; (b) envelope order spectrum.
(a)(b)It can be seen that the order of prominence in Figures28 and 30 is the result of other impacts in the rotor system. Owing to its existence, NEA based on SK fails to extract the features generated by bearing fault impact. Combined with the analysis of Section 4.3.3, Section 5.3, and Section 5.4.2, it is indicated that SK-based NEA fails to extract the weak fault features of bearings in the presence of random noise and impact noise, which exactly highlights the unique advantages of MCM-CS in face of such difficulties.


## 6. Conclusions

In this paper, a VMD-based CS method is proposed to extract the weak fault feature of rolling element bearings under variable speeds to overcome the inadequate sensing efficiency. First of all, order tracking is employed to eliminate the nonstationarity of the signal caused by the speed change. Secondly, through VMD, the angular vibration signal from the previous step can be decomposed into a series of modes, before which Particle Swarm Optimization algorithm is involved in the optimization of VMD parameters to retain bearing fault signals as much as possible. The maximum energy ratio mode is a great substitute for the next compression detection, because of the enhanced sparsity. In the application of CS, the mode matrix is used to construct MCM as the measurement matrix to compress the denoised signal. Based on the compressed signal, the orthogonal matching pursuit is employed to detect the sparse coefficients corresponding to fault feature where the Fourier dictionary is selected to construct sensing matrix and the sparsity is determined as 4 by signal feature. Finally, for feature visualization, the recovered signal is obtained from the detected result using the Fourier dictionary, and the demodulation is employed to extract the fault feature which characterizes the performance of rolling element bearings. The simulated signals and experimental data are utilized to validate that the sensing efficiency of the proposed method is significantly promoted. Despite the strong background noise, the proposed scheme can successfully extract the weak fault feature of rolling element bearings operating at variable speed.



---

*Source: 1745184-2020-08-13.xml*