# Fault Diagnosis of Variable Load Bearing Based on Quantum Chaotic Fruit Fly VMD and Variational RVM

Authors: Bo Xu; Huipeng Li; Fengxing Zhou; Baokang Yan; Yi Liu; Yajie Ma
Journal: Shock and Vibration
 (2019)

---

## Abstract

Under normal circumstances, bearings generally run under variable loading conditions. Under such conditions, the vibration signals of the bearing malfunctions are often nonstationary signals, which are difficult to process effectively. In order to accurately and effectively diagnose the failure types and damage degree of bearings under variable load conditions, an intelligent diagnostic model based on the variational mode decomposition (VMD) of quantum chaotic fruit fly optimization algorithm (QCFOA) and a multiclassification variational relevance vector machine (VRVM) is proposed. First, the key parameters of the VMD are selected using the QCFOA. Secondly, the known bearing fault signal is decomposed by the optimized VMD, and the center frequency and marginal spectral entropy (MSE) of each natural modal component are extracted to construct two-dimensional MSE. Then, the probit model is used to replace the logistic model, and a simpler and more practical multiclassification model is constructed. The two-dimensional MSE of each intrinsic modal component is used as a learning sample for VRVM. Finally, the bearing fault data under 1 hp load are taken as training samples, and the bearing fault data under two loads of 0 hp and 3 hp are used as test samples to verify the effectiveness of the intelligent diagnosis model. The experimental results show that the intelligent fault diagnosis method proposed in this paper can accurately diagnose the type of fault and the degree of damage and has high robustness.


---

## Body

## 1. Introduction

Rolling bearings, which are the core components of rotating machinery, operate under harsh conditions and are easily damaged. It is necessary to monitor the fault. Once a fault occurs, it may cause huge economic losses and personal safety problems. In recent years, with the continuous update of knowledge in the field of digital signal processing and machine learning, intelligent fault diagnosis technology has become a major development trend [1]. Intelligent diagnosis is essentially a pattern recognition process, including two important steps of fault feature extraction and fault identification. In particular, how to effectively extract weak fault characteristics in the fault signal is key to fault diagnosis.The vibration signal analysis is widely used in the fault diagnosis of bearing. In general, the vibration signals of the rolling bearing fault are mostly nonstationary signals, and the method suitable for processing nonstationary signals should be adopted. Empirical mode decomposition (EMD) [2] as a powerful tool of nonstationary signal processing has received extensive attention from researchers concerned with mechanical fault diagnosis. The method of bearing fault feature extraction based on EMD has been widely applied [3–5]. Inspired by the EMD method, Smith et al in [6] proposed another adaptive signal decomposition method named the “local mean decomposition (LMD)” in 2005, which has aroused the attention of a large number of scholars [7]. The EMD and LMD are excellent self-adaptive processing methods for nonstationary and nonlinear signals. However, they have unavoidable deficiencies, such as the end effect, overshoot, and mode mixing. Recently, many scholars proposed optimization and improvement methods for such deficiencies. However, both of them were regarded as recursive mode decomposition methods, and inherent defects are difficult to solve fundamentally. Thus, a new adaptive signal processing method named the “variation mode decomposition (VMD)” was proposed by Dragomiretskiy and Zosso [8]. The method abandons the recursive decomposition mode and effectively alleviates or avoids a series of deficiencies in the EMD and LMD methods.The excellent features presented by VMD have been used by scholars in the field of fault diagnosis [9]. However, an important feature of VMD is that the number of intrinsic modal components and their penalty parameters must be set in advance; once unsuitable parameter values are selected, the decomposition results will be seriously affected. Therefore, against the selection of key parameter values of VMD, some scholars introduced the particle swarm optimization algorithm to select the key parameters of VMD optimally and achieved certain results [10]. However, the traditional heuristic optimization algorithms, such as particle swarm optimization (PSO) [11], genetic algorithm (GA) [12], ant colony optimization (ACO) [13], and cuckoo algorithm (CA) [14], have some problems such as parameter dependence, computational complexity, convergence speed, and optimization accuracy, which restrict their practical applications. To solve this problem, this paper introduces the fruit fly optimization algorithm (FOA) [15] to optimize the key parameters of VMD. This algorithm is a new swarm intelligence algorithm based on the bionics principle of fruit fly foraging behavior. It is applied in many fields [16–18]. However, its convergence accuracy is very sensitive to the initial value. Once the initial value is not selected properly, the search is likely to fall into a local optimum, and the convergence accuracy is low. This paper combines the quantum logistic chaotic mapping [19] and FOA and propose a quantum chaotic fruit fly algorithm in the three-dimensional space search, which strengthens the ergodicity, avoids the search process falling into the local optimal value, and improves the search efficiency. Then, using the local minimum value of the MSE as the fitness function, the quantum chaotic fruit fly optimization algorithm (QCFOA) is used to search two key parameters of VMD and obtain the optimal combination value of the key parameters. Furthermore, the optimized VMD is used to process the known fault signals, and the effective intrinsic mode function (IMF) component and its MSE are obtained.The traditional diagnosis method usually trains the fault diagnosis model under a certain load and has great limitations in diagnosing the state of equipment under a certain load. In practical applications, many mechanical devices work under variable load conditions. Both the load and the damage degree will affect the amplitude of the fault characteristic frequency of the bearing. Therefore, the single MSE cannot effectively characterize the damage degree of the fault in the variable load condition (the radial load is mainly discussed in this paper). Because the bearing is running under different load conditions, the normal contact load between the rolling body and the raceway will change, which leads to the change in the natural vibration frequency of the bearing. To address these challenges, the center frequency is introduced in this paper, and the one-dimensional MSE is extended into two-dimensional MSE as the learning sample of the variational correlation vector machine. Then, the method is used to identify the various types of faults and the damage degree of the bearing under variable load.The technology of intelligent diagnosis for mechanical faults is changing fast. Artificial neural network (ANN) and support vector machine (SVM), as intelligent recognizers, have received the most extensive attention, and a multitude of recent research efforts have been made to explore the mechanical fault diagnosis. However, the ANN algorithm requires a large number of training samples and has some inherent problems, such as black box operation, low generalization ability, and overlearning [20–24]. Similarly, the SVM also has some unavoidable deficiencies. The method cannot get the probabilistic prediction and the uncertainty in prediction [25–29].RVM is a new machine learning algorithm based on support vector machine (SVM) and Bayesian theory framework [30]. Compared with the SVM method, it can directly give the uncertainty of the result while giving the diagnosis results. The RVM training process needs fewer parameters, and its solution is more sparse [31]. So, the probability output of RVM accords with the actual mechanical fault diagnosis process and has high applied research value. However, when the standard RVM is very limited in the size of the data sample, the computational cost of the training is very high. For this problem, Bishop [32] proposed a method of computing and solving RVM by means of the variational method, named the “VRVM.” In the case of very limited data samples, this method is better than the type II maximum edge likelihood estimation and can give the posterior distribution of parameters and hyperparameters. Compared with the standard RVM, the practicability and performance of the VRVM are better. In addition, VRVM is the same as the standard RVM, and its classification and regression are all mapped by logistic function. Therefore, when the regression problem is converted into a classification problem, the noise variable must be ignored, and the true value of the model cannot be accurately estimated. Therefore, the probit model is used instead of the logistic model in this paper, which makes the classification problem and the regression problem organically combined to avoid the approximate derivation of the logistic model from the continuous output to the discrete output mapping and to reduce the amount of computation.Finally, the proposed method is used to diagnose the variable load fault data collected from the failure platform to verify the effectiveness and robustness of the method.


## 2. The Principle of VMD

In the VMD algorithm, the IMF is redefined as an AM-FM signal, which removes the loop iteration method used by the EMD algorithm. Instead, the signal decomposition process is transferred to the variational structure. By constructing and solving the constrained variational problem and decomposing the original signal into a specified number of IMF components, the construction process of the corresponding variational problem is summarized as follows.For each IMF componentukt, the following analytic signal is obtained through the Hilbert transform:(1)δt+jπt∗ukt.For each analytic signal, a central frequencyωk is estimated, and the frequency of each analytical signal is transformed to the baseband by shifting frequency:(2)δt+jπt∗ukte−jωkt.The Gaussian smoothing index of the frequency shift signal is used to estimate the bandwidth of each IMF component, and then the corresponding constraint variational model is expressed as(3)minuk,ωk∑k∂tδt+jπt∗ukte−jωkt22,s.t.∑kuk=f,where uk=u1,u2,…,uk represents the decomposition of the IMF component, k is the number, ωk=ω1,ω2,…,ωk represents the frequency center of each component, f is the input signal, and e−jωkt is the estimated center frequency. In order to obtain the abovementioned constraint variation problem, the augmented Lagrange function is introduced as follows:(4)Luk,ωk,λ=α∑∂tδt+jπt∗ukte−jωt22+ft−∑kukt22+λt,ft−∑kukt,where α is the quadratic penalty parameter and λ is the Lagrange multiplier. The saddle point of the augmented Lagrange function is obtained by using the alternating direction multiplier algorithm, which is the optimal solution of equation (2) constraint variational model, and the original signal is decomposed into K narrowband IMF components. The solution procedure for the variational model is as follows (Algorithm 1).Algorithm 1:Complete optimization of VMD.
Initialize:u^k1, ωk1, λ^1, n⟵0
repeat:
n
⟵
n
+
1
for k=1:K do
Update u^k for all ω≥0:
u
^
k
n
+
1
ω
⟵
f
^
ω
−
∑
i
<
k
u
^
i
n
+
1
ω
−
∑
i
>
k
u
^
i
n
ω
+
λ
^
n
ω
/
2
/
1
+
2
α
ω
−
ω
k
n
2
Update ωk:
ω
k
n
+
1
⟵
∫
0
∞
ω
u
^
k
n
+
1
ω
2
d
ω
/
∫
0
∞
u
^
k
n
+
1
ω
2
d
ω
end for
Dual ascent for all ω≥0:
λ
^
n
+
1
ω
⟵
λ
^
n
ω
+
τ
f
^
ω
−
∑
k
u
^
k
n
+
1
ω
until convergence: ∑ku^kn+1−u^kn22/u^kn22<εIt is known from the variational model solution process that the performance of VMD is closely related to the decomposition parameters, such as the total number of modalitiesK and the secondary penalty α.The performance of VMD is very sensitive to the value ofK. If the value of K is too small, the data will be undersegmented and some components will be contained in other modalities, and if the value of K is too large, problems such as modal copying will occur. If the value of α is too small, the bandwidth of the modal component will be too large; some components will be included in other modal components, or additional “noise” will be captured; if the value of α is too large, the bandwidth of the modal component will be too small and some components in the original signal will be lost.Therefore, the design of an optimal VMD should be focused on how to obtain the optimal combination value of key parameters of the VMD.


## 3. VMD Optimization Based on Quantum Chaotic FOA

The VMD algorithm needs to set the number of IMF components in advance when processing signals. DifferentK values will result in different decomposition results. In addition, the penalty parameter α also has a great influence on the decomposition result. The smaller the α, the larger the bandwidth of each IMF component; conversely, the bandwidth of the component signal is smaller. Therefore, parameters K and α are two critical parameters that affect VMD performance. In practical application, the components of the bearing fault signal are very complex, and the selection of suitable parameters K and α is key to the effective extraction of bearing fault characteristics by the VMD algorithm. Meanwhile, α and K have interactive effects on the performance of VMD. Therefore, a swarm intelligence algorithm which can optimize parameters α and K is needed to avoid the contingency and blindness of manually setting parameters.


### 3.1. Quantum Chaotic Mapping

Quantum chaotic system has the natural properties of the classical chaotic system [33]. For the same classical chaotic system, different quantitative criteria can produce different quantum chaotic maps. In the work of Goggin [34], the classical logistic system is quantified by the recoil rotor model, and the corresponding quantum logistic mapping is obtained. The definition is as follows:(5)xn+1=μxn−xn2−μyn,yn+1=−yne−2β+e−βμ2−xn−x˜nyn−xnz˜n−x˜nznzn+1=−zne−2β+e−βμ21−x˜nzn−2xnyn−xn,,where μ is a chaos control parameter, β is a dissipative parameter, xn, yn, and zn are the state values of the system, and x˜n and z˜n are complex conjugation of xn and zn, respectively. If the initial value in the chaotic system is real, then the chaotic sequence generated by the system is real, and there are x˜n=xn and z˜n=zn. In [35–37], it is proved that the pseudorandom sequence based on the quantum chaotic mapping not only has all the advantages of the traditional chaotic system but also has a weaker correlation and stronger ergodicity than the traditional chaotic system.


### 3.2. Quantum Chaotic FOA

The advantage of the FOA is that it is easy to understand, simple to search, and easy to implement. Therefore, it is widely used in parameter optimization problems [38, 39]. In this paper, a quantum logistic chaotic map is used to extend the search space of the FOA into a three-dimensional search space, and the location of the fruit fly group is initialized by the better characteristics of nonperiodicity, ergodicity, and class randomness and more sensitivity to system parameters and initial conditions than traditional chaotic systems. The method improves the diversity of population and strengthens the ergodicity of the search, avoiding the premature convergence of the search process to the local optimum and improving search efficiency. The sketch map of its three-dimensional search space foraging is shown in Figure 1.Figure 1
Three-dimensional space foraging map of the fruit fly swarm.The steps of the QCFOA are shown as follows:Step 1.
Initialize the fruit fly swarm locationXaxis,Yaxis,Zaxis=X0,Y0,Z0 = X0,Y0,Z0 with random function rand⋅ randomly and parameters of the first iteration, including the maximum number of iterations T, number of fruit flies N, chaos control parameter μ, and dissipation parameter β.Step 2.
Update the positionXi,Yi,Zi of each fruit fly by using equation (5):(6)Xi=Xaxis+rand1,1,Yi=Yaxis+rand1,1,Zi=Zaxis+rand1,1,where rand⋅∈0,1 represents the random variable of uniform distribution, i=1,2,…,N. Substituting equation (5) into rand⋅ of equation (6), we can get the following:(7)Xi=Xaxis+Xi′,Yi=Yaxis+Yi′,Zi=Zaxis+Zi′,where(8)Xi′=μXaxis−Xaxis2−μYaxis,Yi′=−Yaxise−2β+e−βμ2−2XaxisYaxis−2XaxisZaxis,Zi′=−Zaxise−2β+e−βμ21−XaxisZaxisn−2XaxisYaxis−Xaxis.Step 3.
Each component from Step 2 is mapped to the value of odor concentration judgment via equation (7):(9)Distancei=Xi2+Yi2+Zi2,Si=1Distancei.Step 4.
Calculate odor concentration values based on fitness function:(10)Smelli=fitnessfunctionSi.Step 5.
Select maximum odor concentration:(11)bestSmellIndex=selectionmaxSmell.Step 6.
Update maximum odor concentration:(12)Xi=XbestSmellIndex,Yi=YbestSmellIndex,Zi=ZbestSmellIndex.Step 7.
Check the termination condition: this step compares the current maximum odor concentration to the previous maximum odor concentration. When the current maximum concentration is no longer superior to the previous one, or the current iteration is equal to the maximum number of iterations, the iteration process is terminated. Otherwise, execute Step 2.


### 3.3. VMD Model Based on Quantum Chaotic FOA

#### 3.3.1. The Fitness Function

The decomposition capability of the VMD is heavily determined by the selected parametersα and K. For nonstationary and nonlinear signal processing, it is not feasible to search the optimal parameters of α,K artificially. In the study, the optimization algorithm of quantum chaotic fruit fly is used to search the optimal parameters of VMD. But, before the optimization, a fitness function needs to be determined.Information entropy is a measure of uncertainty in information quality, which represents the average uncertainty of signals. The greater the entropy, the greater the uncertainty of the signal and the more complex the signal. This paper defines the VMD marginal spectral entropyHp of the signal xt based on the information entropy, which can represent the uncertainty of the signal in the frequency domain and measure the complexity of signal frequency.Hilbert transformation for each IMF componentcit of VMD is as follows:(13)Hcit=1π∫−∞+∞ciτt−τdτ.The analytic signalsit is constructed as follows:(14)sit=cit+jHcit=aitejφit,and we can obtain the following:(15)xt=Re∑i=1naitejφit=Re∑i=1naitej2π∫fitdt,where n is the number of IMF components, ait is the instantaneous amplitude, φit is the instantaneous phase, and fit represents the instantaneous frequency.The amplitude and frequency of the Hilbert transform are time-domain functions. The amplitude of the signalxt can be expressed as a function of time and frequency in the three-dimensional space, and it is the Hilbert amplitude spectrum:(16)Hf,t=Re∑i=1naitej2π∫fitdt.The marginal spectrum of the signalxt can be obtained by the time integral of Hf,t:(17)hf=∫−∞+∞Hf,tdt.The variation rule ofxt amplitude with frequency is described by equation (15). The marginal spectrum entropy of the signal xt is defined as follows:(18)Hp=−∑i=1Npilnpi,pi=hi∑i=1Nhi,where pi is the probability of the corresponding amplitude of the ith frequency and hi is the marginal spectrum of the ith IMF component. In order to facilitate analysis, the marginal spectrum entropy value is normalized: HE=Hp/lnL, and the value range is [0,1]. L determines the length of hi sequence.


#### 3.3.2. Proposed Improved Algorithm Framework

When the bearing is in fault, the vibration signal mainly appears as a periodic impact signal. Therefore, when using the VMD algorithm to deal with bearing fault signals, multiple IMF components will be obtained. According to the definition of marginal spectrum entropy in equation (18), if the IMF component contains more noise components, its corresponding MSE value is larger. Conversely, if the IMF component mainly contains the periodic impact component of bearing failure, the marginal spectrum entropy of VMD is very small. Therefore, the minimum value of the marginal spectrum entropy is expressed as a local minimum, and the corresponding IMF component is taken as the optimal component. The local minimum marginal spectral entropy (LMMSE) of the VMD is shown as follows:(19)minLHEIMF=min−∑i=1Lpilnpi.In this study, the LMMSE of the VMD is taken as the fitness value in the optimization process, and the global optimal IMF component is searched. The local minimum value of the marginal spectrum entropy is used as the final optimization target. The optimization process is the same as the FOA, and the process is carried out as follows:Step 8.
Initialize the fruit fly swarm parameters of the first iteration, including the fruit fly group numberN, number of iterations T, and maximum iterations Tmax.Step 9.
The parameters in equation (3) are set as μ = 3.99 and β ≥ 6.Step 10.
Generate two large groups ofUα and UK randomly by using equation (7), and the number of fruit flies in each group is N. The positions UXiα,Yiα,Ziα and UXiK,YiK,ZiK of each fruit fly are then continuously updated using equation (7).Step 11.
CalculateDistancei,α according to equation (5), and it is represented by Di,α and Di,K. Set Si,K=1/Di,K and Si,α=1/Di,α, and Si is represented by Si,α and Si,K. Get the odor concentration Smelli by the fitness value, and obtain the best smell concentration value Smellbest and the corresponding location Si,α,Si,K.Step 12.
Enter iterative optimization to repeat the implementation of Step 3,T=T+1.When the current smell concentration is not superior to the previous iterative smell concentration any more or the iteration number reaches the maximum number of iterations Tmax, the circulation is terminated, and return optimal combination parameters αopt,Kopt. Otherwise, carry out Step 3.


### 3.4. Simulations and Analysis

#### 3.4.1. Simulation Signal

In practice, rolling bearings generate periodic impact signals when pitting or cracking occurs. But in the early stage of failure, the impact signal is very weak and is easily drowned by noise, so it is difficult to find fault characteristic frequency in traditional signal analysis methods. This paper proposes a quantum chaotic fruit fly algorithm to search the optimal parameters of VMD, and the flow chart of algorithm is shown in Figure2. In order to qualitatively and quantitatively analyze the validity and superiority of this method, the simulation signal of the early fault of the bearing under the simulated strong noise background is analyzed [30], and the expression of the simulation signal is as follows:(20)xt=x1t+x2t+xnt,where x1t is an impact simulation signal with periodic pulse attenuation with a frequency of 12 Hz and maximum amplitude of 0.5 V, x2t is a cosine combined signal with a frequency of 35 Hz and 15 Hz, and xnt is a Gauss white noise signal. The time-domain waveform of the simulation signal is shown in Figure 3.Figure 2
Simulation signal of early fault of the bearing.Figure 3
Search process diagram of three optimization algorithms.From the simulation signal of early fault of the bearing in Figure3, it can be seen that the impact attenuation signal in the early stage of simulated bearing failure is almost submerged by low-frequency components and strong noise due to its small amplitude.


#### 3.4.2. Parameter Optimization Analysis

In order to verify the effectiveness and superiority of the improved chaotic FOA parameter optimization proposed in this paper, the performance is compared with that of classical PSO, QPSO, FOA, and CFOA. In this experiment, the number of fruit flies is set asN = 30 and the maximum number of iterations is Tmax = 200. Other parameters are selected according to relevant literature to ensure the best results of each algorithm. The combination parameters αopt,Kopt are optimized by the above method. Figure 4 and Table 1 show that different algorithms have different effects on the solution.Figure 4
Flow chart of the optimal VMD algorithm.Table 1
Computation efficiency comparison among PSO, QPSO, FOA, CFOA, and QCFOA methods.
Algorithms
Global optimal solution
Iterative number
[α, K]
Reconstruction error
PSO
0.0602
110
[1207, 7]
3.01 × 10−4
QPSO
0.0501
94
[1150, 6]
1.47 × 10−4
FOA
0.0454
75
[1040, 6]
1.15 × 10−4
CFOA
0.0443
40
[970, 5]
6.17 × 10−5
QCFOA
0.0236
18
[690, 4]
1.02 × 10−5From Table1, it can be observed that the QCFOA method achieves the minimum number of iterations and the best global optimal solution. The searched VMD parameter combination is [690, 4], and the reconstruction error is 1.02 × 10−5. Figure 4 demonstrates the iterative process using PSO, QPSO, FOA, CFOA, and QCFOA. As shown in Figure 4, the convergence speed of the QCFOA algorithm is the fastest and the reconstruction error of the signal is the smallest.In the search process, the LMMSE changes with the evolution of the population, and the corresponding VMD optimal parameter combination is shown in Table1.From Table1, it can be observed that the QCFOA method achieves the minimum number of iterations and the best global optimal solution. The searched parameter combination of VMD [α,K] is [690, 4], its reconstruction error is 1.02 × 10−5, and it has the least number of iterations. Figure 4 demonstrates the search process by using PSO, QPSO, FOA, CFOA, and QCFOA. The experimental results indicate that the convergence speed of the QCFOA is faster than others. Therefore, the experimental results indicate that the proposed method is more effective and superior than the mentioned optimization algorithms.


#### 3.4.3. Comparison and Analysis of EEMD, LMD, and VMD Methods

In order to validate the effectiveness and superiority of the VMD method based on the QCFOA, EEMD, LMD, and VMD methods are used, respectively, to process the simulation signals in Figure3, and the results are shown in Figure 5.Figure 5
The characteristic components and corresponding marginal spectrum of the simulation signal based on (a) VMD, (b) LMD, and (c) EEMD methods.
(a)
(b)
(c)From Figure5, it can be seen that the EEMD and LMD methods can effectively extract the characteristic frequencies of the low-frequency cosine components in the simulation signals including 25 Hz and 15 Hz but cannot extract the characteristic frequency of the weak impact signal of 12 Hz, whereas there is modal mixing which is an unavoidable deficiency in the decomposition of EEMD and LMD. The low-frequency cosine components of 25 Hz and 15 Hz appear in different components. In addition, the frequency amplitude of the extracted 12 Hz impact signal is also relatively weak. The VMD method based on the QCFOA proposed in this paper not only can effectively extract the characteristic frequency of the low-frequency cosine components of 25 Hz and 15 Hz but also can effectively extract the characteristic frequency of the weak impact signal of 12 Hz and its corresponding doubling frequency. The number of signal components decomposed by the improved VMD method is also significantly less than that obtained by EEMD and LMD. The experimental results show the better effectiveness and superiority of the proposed method.


## 3.1. Quantum Chaotic Mapping

Quantum chaotic system has the natural properties of the classical chaotic system [33]. For the same classical chaotic system, different quantitative criteria can produce different quantum chaotic maps. In the work of Goggin [34], the classical logistic system is quantified by the recoil rotor model, and the corresponding quantum logistic mapping is obtained. The definition is as follows:(5)xn+1=μxn−xn2−μyn,yn+1=−yne−2β+e−βμ2−xn−x˜nyn−xnz˜n−x˜nznzn+1=−zne−2β+e−βμ21−x˜nzn−2xnyn−xn,,where μ is a chaos control parameter, β is a dissipative parameter, xn, yn, and zn are the state values of the system, and x˜n and z˜n are complex conjugation of xn and zn, respectively. If the initial value in the chaotic system is real, then the chaotic sequence generated by the system is real, and there are x˜n=xn and z˜n=zn. In [35–37], it is proved that the pseudorandom sequence based on the quantum chaotic mapping not only has all the advantages of the traditional chaotic system but also has a weaker correlation and stronger ergodicity than the traditional chaotic system.


## 3.2. Quantum Chaotic FOA

The advantage of the FOA is that it is easy to understand, simple to search, and easy to implement. Therefore, it is widely used in parameter optimization problems [38, 39]. In this paper, a quantum logistic chaotic map is used to extend the search space of the FOA into a three-dimensional search space, and the location of the fruit fly group is initialized by the better characteristics of nonperiodicity, ergodicity, and class randomness and more sensitivity to system parameters and initial conditions than traditional chaotic systems. The method improves the diversity of population and strengthens the ergodicity of the search, avoiding the premature convergence of the search process to the local optimum and improving search efficiency. The sketch map of its three-dimensional search space foraging is shown in Figure 1.Figure 1
Three-dimensional space foraging map of the fruit fly swarm.The steps of the QCFOA are shown as follows:Step 1.
Initialize the fruit fly swarm locationXaxis,Yaxis,Zaxis=X0,Y0,Z0 = X0,Y0,Z0 with random function rand⋅ randomly and parameters of the first iteration, including the maximum number of iterations T, number of fruit flies N, chaos control parameter μ, and dissipation parameter β.Step 2.
Update the positionXi,Yi,Zi of each fruit fly by using equation (5):(6)Xi=Xaxis+rand1,1,Yi=Yaxis+rand1,1,Zi=Zaxis+rand1,1,where rand⋅∈0,1 represents the random variable of uniform distribution, i=1,2,…,N. Substituting equation (5) into rand⋅ of equation (6), we can get the following:(7)Xi=Xaxis+Xi′,Yi=Yaxis+Yi′,Zi=Zaxis+Zi′,where(8)Xi′=μXaxis−Xaxis2−μYaxis,Yi′=−Yaxise−2β+e−βμ2−2XaxisYaxis−2XaxisZaxis,Zi′=−Zaxise−2β+e−βμ21−XaxisZaxisn−2XaxisYaxis−Xaxis.Step 3.
Each component from Step 2 is mapped to the value of odor concentration judgment via equation (7):(9)Distancei=Xi2+Yi2+Zi2,Si=1Distancei.Step 4.
Calculate odor concentration values based on fitness function:(10)Smelli=fitnessfunctionSi.Step 5.
Select maximum odor concentration:(11)bestSmellIndex=selectionmaxSmell.Step 6.
Update maximum odor concentration:(12)Xi=XbestSmellIndex,Yi=YbestSmellIndex,Zi=ZbestSmellIndex.Step 7.
Check the termination condition: this step compares the current maximum odor concentration to the previous maximum odor concentration. When the current maximum concentration is no longer superior to the previous one, or the current iteration is equal to the maximum number of iterations, the iteration process is terminated. Otherwise, execute Step 2.


## 3.3. VMD Model Based on Quantum Chaotic FOA

### 3.3.1. The Fitness Function

The decomposition capability of the VMD is heavily determined by the selected parametersα and K. For nonstationary and nonlinear signal processing, it is not feasible to search the optimal parameters of α,K artificially. In the study, the optimization algorithm of quantum chaotic fruit fly is used to search the optimal parameters of VMD. But, before the optimization, a fitness function needs to be determined.Information entropy is a measure of uncertainty in information quality, which represents the average uncertainty of signals. The greater the entropy, the greater the uncertainty of the signal and the more complex the signal. This paper defines the VMD marginal spectral entropyHp of the signal xt based on the information entropy, which can represent the uncertainty of the signal in the frequency domain and measure the complexity of signal frequency.Hilbert transformation for each IMF componentcit of VMD is as follows:(13)Hcit=1π∫−∞+∞ciτt−τdτ.The analytic signalsit is constructed as follows:(14)sit=cit+jHcit=aitejφit,and we can obtain the following:(15)xt=Re∑i=1naitejφit=Re∑i=1naitej2π∫fitdt,where n is the number of IMF components, ait is the instantaneous amplitude, φit is the instantaneous phase, and fit represents the instantaneous frequency.The amplitude and frequency of the Hilbert transform are time-domain functions. The amplitude of the signalxt can be expressed as a function of time and frequency in the three-dimensional space, and it is the Hilbert amplitude spectrum:(16)Hf,t=Re∑i=1naitej2π∫fitdt.The marginal spectrum of the signalxt can be obtained by the time integral of Hf,t:(17)hf=∫−∞+∞Hf,tdt.The variation rule ofxt amplitude with frequency is described by equation (15). The marginal spectrum entropy of the signal xt is defined as follows:(18)Hp=−∑i=1Npilnpi,pi=hi∑i=1Nhi,where pi is the probability of the corresponding amplitude of the ith frequency and hi is the marginal spectrum of the ith IMF component. In order to facilitate analysis, the marginal spectrum entropy value is normalized: HE=Hp/lnL, and the value range is [0,1]. L determines the length of hi sequence.


### 3.3.2. Proposed Improved Algorithm Framework

When the bearing is in fault, the vibration signal mainly appears as a periodic impact signal. Therefore, when using the VMD algorithm to deal with bearing fault signals, multiple IMF components will be obtained. According to the definition of marginal spectrum entropy in equation (18), if the IMF component contains more noise components, its corresponding MSE value is larger. Conversely, if the IMF component mainly contains the periodic impact component of bearing failure, the marginal spectrum entropy of VMD is very small. Therefore, the minimum value of the marginal spectrum entropy is expressed as a local minimum, and the corresponding IMF component is taken as the optimal component. The local minimum marginal spectral entropy (LMMSE) of the VMD is shown as follows:(19)minLHEIMF=min−∑i=1Lpilnpi.In this study, the LMMSE of the VMD is taken as the fitness value in the optimization process, and the global optimal IMF component is searched. The local minimum value of the marginal spectrum entropy is used as the final optimization target. The optimization process is the same as the FOA, and the process is carried out as follows:Step 8.
Initialize the fruit fly swarm parameters of the first iteration, including the fruit fly group numberN, number of iterations T, and maximum iterations Tmax.Step 9.
The parameters in equation (3) are set as μ = 3.99 and β ≥ 6.Step 10.
Generate two large groups ofUα and UK randomly by using equation (7), and the number of fruit flies in each group is N. The positions UXiα,Yiα,Ziα and UXiK,YiK,ZiK of each fruit fly are then continuously updated using equation (7).Step 11.
CalculateDistancei,α according to equation (5), and it is represented by Di,α and Di,K. Set Si,K=1/Di,K and Si,α=1/Di,α, and Si is represented by Si,α and Si,K. Get the odor concentration Smelli by the fitness value, and obtain the best smell concentration value Smellbest and the corresponding location Si,α,Si,K.Step 12.
Enter iterative optimization to repeat the implementation of Step 3,T=T+1.When the current smell concentration is not superior to the previous iterative smell concentration any more or the iteration number reaches the maximum number of iterations Tmax, the circulation is terminated, and return optimal combination parameters αopt,Kopt. Otherwise, carry out Step 3.


## 3.3.1. The Fitness Function

The decomposition capability of the VMD is heavily determined by the selected parametersα and K. For nonstationary and nonlinear signal processing, it is not feasible to search the optimal parameters of α,K artificially. In the study, the optimization algorithm of quantum chaotic fruit fly is used to search the optimal parameters of VMD. But, before the optimization, a fitness function needs to be determined.Information entropy is a measure of uncertainty in information quality, which represents the average uncertainty of signals. The greater the entropy, the greater the uncertainty of the signal and the more complex the signal. This paper defines the VMD marginal spectral entropyHp of the signal xt based on the information entropy, which can represent the uncertainty of the signal in the frequency domain and measure the complexity of signal frequency.Hilbert transformation for each IMF componentcit of VMD is as follows:(13)Hcit=1π∫−∞+∞ciτt−τdτ.The analytic signalsit is constructed as follows:(14)sit=cit+jHcit=aitejφit,and we can obtain the following:(15)xt=Re∑i=1naitejφit=Re∑i=1naitej2π∫fitdt,where n is the number of IMF components, ait is the instantaneous amplitude, φit is the instantaneous phase, and fit represents the instantaneous frequency.The amplitude and frequency of the Hilbert transform are time-domain functions. The amplitude of the signalxt can be expressed as a function of time and frequency in the three-dimensional space, and it is the Hilbert amplitude spectrum:(16)Hf,t=Re∑i=1naitej2π∫fitdt.The marginal spectrum of the signalxt can be obtained by the time integral of Hf,t:(17)hf=∫−∞+∞Hf,tdt.The variation rule ofxt amplitude with frequency is described by equation (15). The marginal spectrum entropy of the signal xt is defined as follows:(18)Hp=−∑i=1Npilnpi,pi=hi∑i=1Nhi,where pi is the probability of the corresponding amplitude of the ith frequency and hi is the marginal spectrum of the ith IMF component. In order to facilitate analysis, the marginal spectrum entropy value is normalized: HE=Hp/lnL, and the value range is [0,1]. L determines the length of hi sequence.


## 3.3.2. Proposed Improved Algorithm Framework

When the bearing is in fault, the vibration signal mainly appears as a periodic impact signal. Therefore, when using the VMD algorithm to deal with bearing fault signals, multiple IMF components will be obtained. According to the definition of marginal spectrum entropy in equation (18), if the IMF component contains more noise components, its corresponding MSE value is larger. Conversely, if the IMF component mainly contains the periodic impact component of bearing failure, the marginal spectrum entropy of VMD is very small. Therefore, the minimum value of the marginal spectrum entropy is expressed as a local minimum, and the corresponding IMF component is taken as the optimal component. The local minimum marginal spectral entropy (LMMSE) of the VMD is shown as follows:(19)minLHEIMF=min−∑i=1Lpilnpi.In this study, the LMMSE of the VMD is taken as the fitness value in the optimization process, and the global optimal IMF component is searched. The local minimum value of the marginal spectrum entropy is used as the final optimization target. The optimization process is the same as the FOA, and the process is carried out as follows:Step 8.
Initialize the fruit fly swarm parameters of the first iteration, including the fruit fly group numberN, number of iterations T, and maximum iterations Tmax.Step 9.
The parameters in equation (3) are set as μ = 3.99 and β ≥ 6.Step 10.
Generate two large groups ofUα and UK randomly by using equation (7), and the number of fruit flies in each group is N. The positions UXiα,Yiα,Ziα and UXiK,YiK,ZiK of each fruit fly are then continuously updated using equation (7).Step 11.
CalculateDistancei,α according to equation (5), and it is represented by Di,α and Di,K. Set Si,K=1/Di,K and Si,α=1/Di,α, and Si is represented by Si,α and Si,K. Get the odor concentration Smelli by the fitness value, and obtain the best smell concentration value Smellbest and the corresponding location Si,α,Si,K.Step 12.
Enter iterative optimization to repeat the implementation of Step 3,T=T+1.When the current smell concentration is not superior to the previous iterative smell concentration any more or the iteration number reaches the maximum number of iterations Tmax, the circulation is terminated, and return optimal combination parameters αopt,Kopt. Otherwise, carry out Step 3.


## 3.4. Simulations and Analysis

### 3.4.1. Simulation Signal

In practice, rolling bearings generate periodic impact signals when pitting or cracking occurs. But in the early stage of failure, the impact signal is very weak and is easily drowned by noise, so it is difficult to find fault characteristic frequency in traditional signal analysis methods. This paper proposes a quantum chaotic fruit fly algorithm to search the optimal parameters of VMD, and the flow chart of algorithm is shown in Figure2. In order to qualitatively and quantitatively analyze the validity and superiority of this method, the simulation signal of the early fault of the bearing under the simulated strong noise background is analyzed [30], and the expression of the simulation signal is as follows:(20)xt=x1t+x2t+xnt,where x1t is an impact simulation signal with periodic pulse attenuation with a frequency of 12 Hz and maximum amplitude of 0.5 V, x2t is a cosine combined signal with a frequency of 35 Hz and 15 Hz, and xnt is a Gauss white noise signal. The time-domain waveform of the simulation signal is shown in Figure 3.Figure 2
Simulation signal of early fault of the bearing.Figure 3
Search process diagram of three optimization algorithms.From the simulation signal of early fault of the bearing in Figure3, it can be seen that the impact attenuation signal in the early stage of simulated bearing failure is almost submerged by low-frequency components and strong noise due to its small amplitude.


### 3.4.2. Parameter Optimization Analysis

In order to verify the effectiveness and superiority of the improved chaotic FOA parameter optimization proposed in this paper, the performance is compared with that of classical PSO, QPSO, FOA, and CFOA. In this experiment, the number of fruit flies is set asN = 30 and the maximum number of iterations is Tmax = 200. Other parameters are selected according to relevant literature to ensure the best results of each algorithm. The combination parameters αopt,Kopt are optimized by the above method. Figure 4 and Table 1 show that different algorithms have different effects on the solution.Figure 4
Flow chart of the optimal VMD algorithm.Table 1
Computation efficiency comparison among PSO, QPSO, FOA, CFOA, and QCFOA methods.
Algorithms
Global optimal solution
Iterative number
[α, K]
Reconstruction error
PSO
0.0602
110
[1207, 7]
3.01 × 10−4
QPSO
0.0501
94
[1150, 6]
1.47 × 10−4
FOA
0.0454
75
[1040, 6]
1.15 × 10−4
CFOA
0.0443
40
[970, 5]
6.17 × 10−5
QCFOA
0.0236
18
[690, 4]
1.02 × 10−5From Table1, it can be observed that the QCFOA method achieves the minimum number of iterations and the best global optimal solution. The searched VMD parameter combination is [690, 4], and the reconstruction error is 1.02 × 10−5. Figure 4 demonstrates the iterative process using PSO, QPSO, FOA, CFOA, and QCFOA. As shown in Figure 4, the convergence speed of the QCFOA algorithm is the fastest and the reconstruction error of the signal is the smallest.In the search process, the LMMSE changes with the evolution of the population, and the corresponding VMD optimal parameter combination is shown in Table1.From Table1, it can be observed that the QCFOA method achieves the minimum number of iterations and the best global optimal solution. The searched parameter combination of VMD [α,K] is [690, 4], its reconstruction error is 1.02 × 10−5, and it has the least number of iterations. Figure 4 demonstrates the search process by using PSO, QPSO, FOA, CFOA, and QCFOA. The experimental results indicate that the convergence speed of the QCFOA is faster than others. Therefore, the experimental results indicate that the proposed method is more effective and superior than the mentioned optimization algorithms.


### 3.4.3. Comparison and Analysis of EEMD, LMD, and VMD Methods

In order to validate the effectiveness and superiority of the VMD method based on the QCFOA, EEMD, LMD, and VMD methods are used, respectively, to process the simulation signals in Figure3, and the results are shown in Figure 5.Figure 5
The characteristic components and corresponding marginal spectrum of the simulation signal based on (a) VMD, (b) LMD, and (c) EEMD methods.
(a)
(b)
(c)From Figure5, it can be seen that the EEMD and LMD methods can effectively extract the characteristic frequencies of the low-frequency cosine components in the simulation signals including 25 Hz and 15 Hz but cannot extract the characteristic frequency of the weak impact signal of 12 Hz, whereas there is modal mixing which is an unavoidable deficiency in the decomposition of EEMD and LMD. The low-frequency cosine components of 25 Hz and 15 Hz appear in different components. In addition, the frequency amplitude of the extracted 12 Hz impact signal is also relatively weak. The VMD method based on the QCFOA proposed in this paper not only can effectively extract the characteristic frequency of the low-frequency cosine components of 25 Hz and 15 Hz but also can effectively extract the characteristic frequency of the weak impact signal of 12 Hz and its corresponding doubling frequency. The number of signal components decomposed by the improved VMD method is also significantly less than that obtained by EEMD and LMD. The experimental results show the better effectiveness and superiority of the proposed method.


## 3.4.1. Simulation Signal

In practice, rolling bearings generate periodic impact signals when pitting or cracking occurs. But in the early stage of failure, the impact signal is very weak and is easily drowned by noise, so it is difficult to find fault characteristic frequency in traditional signal analysis methods. This paper proposes a quantum chaotic fruit fly algorithm to search the optimal parameters of VMD, and the flow chart of algorithm is shown in Figure2. In order to qualitatively and quantitatively analyze the validity and superiority of this method, the simulation signal of the early fault of the bearing under the simulated strong noise background is analyzed [30], and the expression of the simulation signal is as follows:(20)xt=x1t+x2t+xnt,where x1t is an impact simulation signal with periodic pulse attenuation with a frequency of 12 Hz and maximum amplitude of 0.5 V, x2t is a cosine combined signal with a frequency of 35 Hz and 15 Hz, and xnt is a Gauss white noise signal. The time-domain waveform of the simulation signal is shown in Figure 3.Figure 2
Simulation signal of early fault of the bearing.Figure 3
Search process diagram of three optimization algorithms.From the simulation signal of early fault of the bearing in Figure3, it can be seen that the impact attenuation signal in the early stage of simulated bearing failure is almost submerged by low-frequency components and strong noise due to its small amplitude.


## 3.4.2. Parameter Optimization Analysis

In order to verify the effectiveness and superiority of the improved chaotic FOA parameter optimization proposed in this paper, the performance is compared with that of classical PSO, QPSO, FOA, and CFOA. In this experiment, the number of fruit flies is set asN = 30 and the maximum number of iterations is Tmax = 200. Other parameters are selected according to relevant literature to ensure the best results of each algorithm. The combination parameters αopt,Kopt are optimized by the above method. Figure 4 and Table 1 show that different algorithms have different effects on the solution.Figure 4
Flow chart of the optimal VMD algorithm.Table 1
Computation efficiency comparison among PSO, QPSO, FOA, CFOA, and QCFOA methods.
Algorithms
Global optimal solution
Iterative number
[α, K]
Reconstruction error
PSO
0.0602
110
[1207, 7]
3.01 × 10−4
QPSO
0.0501
94
[1150, 6]
1.47 × 10−4
FOA
0.0454
75
[1040, 6]
1.15 × 10−4
CFOA
0.0443
40
[970, 5]
6.17 × 10−5
QCFOA
0.0236
18
[690, 4]
1.02 × 10−5From Table1, it can be observed that the QCFOA method achieves the minimum number of iterations and the best global optimal solution. The searched VMD parameter combination is [690, 4], and the reconstruction error is 1.02 × 10−5. Figure 4 demonstrates the iterative process using PSO, QPSO, FOA, CFOA, and QCFOA. As shown in Figure 4, the convergence speed of the QCFOA algorithm is the fastest and the reconstruction error of the signal is the smallest.In the search process, the LMMSE changes with the evolution of the population, and the corresponding VMD optimal parameter combination is shown in Table1.From Table1, it can be observed that the QCFOA method achieves the minimum number of iterations and the best global optimal solution. The searched parameter combination of VMD [α,K] is [690, 4], its reconstruction error is 1.02 × 10−5, and it has the least number of iterations. Figure 4 demonstrates the search process by using PSO, QPSO, FOA, CFOA, and QCFOA. The experimental results indicate that the convergence speed of the QCFOA is faster than others. Therefore, the experimental results indicate that the proposed method is more effective and superior than the mentioned optimization algorithms.


## 3.4.3. Comparison and Analysis of EEMD, LMD, and VMD Methods

In order to validate the effectiveness and superiority of the VMD method based on the QCFOA, EEMD, LMD, and VMD methods are used, respectively, to process the simulation signals in Figure3, and the results are shown in Figure 5.Figure 5
The characteristic components and corresponding marginal spectrum of the simulation signal based on (a) VMD, (b) LMD, and (c) EEMD methods.
(a)
(b)
(c)From Figure5, it can be seen that the EEMD and LMD methods can effectively extract the characteristic frequencies of the low-frequency cosine components in the simulation signals including 25 Hz and 15 Hz but cannot extract the characteristic frequency of the weak impact signal of 12 Hz, whereas there is modal mixing which is an unavoidable deficiency in the decomposition of EEMD and LMD. The low-frequency cosine components of 25 Hz and 15 Hz appear in different components. In addition, the frequency amplitude of the extracted 12 Hz impact signal is also relatively weak. The VMD method based on the QCFOA proposed in this paper not only can effectively extract the characteristic frequency of the low-frequency cosine components of 25 Hz and 15 Hz but also can effectively extract the characteristic frequency of the weak impact signal of 12 Hz and its corresponding doubling frequency. The number of signal components decomposed by the improved VMD method is also significantly less than that obtained by EEMD and LMD. The experimental results show the better effectiveness and superiority of the proposed method.


## 4. Kernel Parameter Self-Optimization Variational Relevance Vector Machine

### 4.1. Variational Relevance Vector Machine

Many problems in machine learning belong to the category of supervised learning, giving an input vectorX=xnn=1N and the corresponding target output T=tnn=1N. For regression problems, they can be arbitrary values. For classification problems, they are class labels. For regression, tn can be any value, and for classification, tn is the category label.For the standard RVM regression model, the formula is defined as follows:(21)tn=yxn;w=∑n=1NwnKx,xn+w0+εn=wTϕxn+εn,where ϕxn=1,Kxn,x1,…,Kxn,xNT, wn is the weight parameter of the model, Kx,xn is the kernel function, and an RBF kernel function is selected in this section: Kx,xn=expx−xn2/h2, h>0. εn represents the noise, which obeys normal distribution N0,σ2.If the input vectorsX=xnn=1N are independent and identically distributed, then the likelihood function of the corresponding target output T=tnn=1N is defined as(22)PTX,w,σ2=∏n=1NPtnxn,w,σ2,where the parameter w is Gaussian prior distribution and σ2is noise variance:(23)Pwα=∏n=1NNwn0,αn−1,where α=αn is a hyperparameter vector, and each weight value wn is independently assigned a parameter αn. In order to make parameter learning more flexible, ultra-prior distribution is defined for α and noise variance σ2, respectively. The appropriate prior distribution is the gamma distribution:(24)Pα=∏n=0NGammaαnan,bn,Pσ2=Gammaσ2c,d,Gammaαa,b=Γa−1baαa−1e−bα,where Γa=∫0∞ta−1e−tdt is the gamma function, and it is usually defined that the hyperparameter is a very small value such as an = bn = c = d = 10−4; such a hyperparameter before does not provide information for posterior learning, so the posterior depends entirely on the data.When the model is established, the posterior distribution ofw, α, and σ2 can be obtained by the variable Bayesian method [40]. In the process of iterative solution, most of αn tends to infinity and the corresponding wn is zero, realizing the sparse model.RVM classification and regression essentially use the same framework model, except that the conditional distribution of the target value is changed. For the binary classification, the logistic functionσy=1/1+e−y is used in the continuous latent variable yxn;w, and assuming that PTX is the Bernoulli distribution, the likelihood function is as follows [41]:(25)PTX,w=∏n=1Nσyxn;wtn1−σyxn;w1−tn.It is important to note that, in the classification problem without considering the noise variableεn, it is very difficult to directly use the variational method to solve the above model. In [28], a lower bound is introduced by using the inequality:(26)σyxn;wtn1−σyxn;w1−tn=σzn≥σξnexpzn−ξn2−λξnzn2−ξn2,where λξn=1/4ξn×tanhξn/2, zn=2tn−1yxn;w, εn is a variational parameter, and when εn=zn, equation (26) is established. Finally, the variational method is used to solve the lower boundary.The VRVM method is the posterior distribution of RVM model parameters and superparameters by variational Bayesian (VB) function. In the VB method [42] based on the RVM model, the observed variables are X=xnn=1N and T=tnn=1N, the hidden variables are ynn=1N, and the parameter is wpm,αp,σ2p=1,m=1PM; therefore, the logarithmic edge likelihood function can be written as(27)lnPTX=ln∬∬PTX,y,w,a,σ2dydwdadσ2=ln∬∬Qy,w,α,σ2PTX,y,w,a,σ2Qy,w,a,σ2dydwdadσ2≥∬∬Qy,w,a,σ2lnPTX,y,w,a,σ2Qy,w,a,σ2dydwdadσ2,where Qy,w,a,σ2 is the joint probability distribution function between the hidden variable and parameter. The variational Bayesian algorithm assumes that y, w, a, and σ2 are independent of each other; therefore, the joint probability distribution of four variables can be written approximately as(28)Qy,w,a,σ2≈QyQwQaQσ2.Following the assumption of equation (28), equation (21) can be rewritten as follows:(29)lnPTX≥∬∬QyQwQaQσ2,lnPTX,w,a,σ2QyQwQαQσ2dydwdαdσ2=FQy,Qw,Qa,Qσ2.It can be seen from equation (29) that the log-likelihood function has a lower bound and that the real value can be approximated by maximizing the lower bound FQy,Qw,Qa,Qσ2. Then, the lower bound of the log-likelihood is solved using EM (expectation-maximization), and the posterior distribution of the hidden variable and all parameters is obtained. The lower bound expression is as follows:(30)FQy,Qw,Qa,Qσ2=∬∬QyQwQaQσ2lnPty+lnPyXw,σ−2+lnPw0,α−1+lnPaa,b+lnPσ2c,d−lnQy−lnQw−lnQa−lnQσ2dydwdadσ2.Because the variable distributions ofQy, Qw, Qa, and Qσ2 are all conjugate prior distributions, they have the same distribution form as their posterior distribution, and the following is their probability distribution:(31)Qy=∏n=1NNtnyn;m˜n,σ˜n,(32)Qw=Nw;μ˜w,Σ˜w,(33)Qσ2=Gammaσ2;c˜,d˜,(34)Qa∏p=1PGammaαp;a˜p,b˜pM,where Ntn⋅ represents truncated normal distribution, and the direction of truncation N+⋅ or N−⋅ is determined by tn.According to the theory of the EM algorithm, the posterior distribution of the variables is actual expectation of the logarithm of the complete likelihood function with respect to other variables (indicated by⋅), and then the terms related to the variable are extracted. The posterior distributions of y, w, a, and σ2 are, respectively, solved as follows:Step 13.
Extract the term related to variabley in equation (30):(35)F˜Qy=∫QylnPty+∬QwQσ2lnPyXw,σ−2dwdσ2−lnQydy.
Set the derivative ofQy in equation (35) to be 0; hence,(36)lnPty+∬QwQσ2lnPyXw,σ−2dwdσ2−lnQy−1=0.
Solve equation (8) and get the following equation:(37)m˜n=xnw,σ˜n=σ2−1Step 14.
Extract the term related to variablew in equation (30):(38)F˜Qw=∫Qw∬QyQσ2lnPyXw,σ−2dydσ2+∫QalnPw0,α−1da−lnQwdw.
Set the derivative ofQw in equation (38) to be 0. Hence,(39)Σ˜w=diaga+σ2∑n=1NxnTxn−1,μ˜w=σ2Σ˜w∑n=1NxnTyn.Step 15.
Extract the term related to variableσ2 in equation (29):(40)F˜Qσ2=∫Qσ2∬QyQwlnPyXw,σ−2dydw+∫QalnPσ2c,d−lnQσ2dσ2.
Set the derivative ofQσ2 in equation (40) to be 0. Hence,(41)c˜=c+0.5,d=d+1N∑n=1N12yn2−2ynwTxnT+xnwwTxnT.Step 16.
Extract the term related to variablea in equation (29):(42)F˜Qa=∫Qa∫QwlnPw0,a−1dw+lnPaa,b−lnQada.
Set the derivative ofQa in equation (42) to be 0. Hence,(43)α˜p=αp+12,b˜p=bp+1M∑n=1N12wpm2.
Through the process above, the parameter iteration formulas in each variable’s posterior distribution function are obtained, but they are expressed by the expectations of other variables. Therefore, the expectation of these parameters needs to be found [42].
For equation (32), tntn∈−1,+1 determines the truncation direction of the truncation normal distribution; hence,(44)yn=xnμ˜w+tn1σ21/2normpdfxnμ˜w/σ2−1/2normpdftnxnμ˜w/σ2−1/2,yn2=1σ2+xnμ˜w2+tnxnμ˜wσ2−1/2normpdfxnμ˜w/σ21/2normpdftnxnμ˜w/σ2−1/2,where normpdf is a normal probability density function Moreover,(45)w=μ˜w,wwT=Σ˜w+μ˜wμ˜wT,αp=a˜pb˜p,lnαp=Ψa˜p−lnb˜p,σ2=c˜d˜,σ2=Ψc˜−lnd˜.
TheΨ function is defined as follows:(46)Ψa=dlnΓada.
After the model is trained, for a test samplex∗, its probability of prediction can be calculated as follows:(47)Pt∗=1x∗=∫Pt∗=1y∗Py∗x∗,w,σ2dy∗=∫0∞Ny∗;x∗w,σ2−1dy∗=normcdfx∗μ˜wσ2−1/2.
wherenormcdf is a normal cumulative integral distribution function. According to the prediction probability, the test sample can be classified and identified. When Pt∗=1x∗≥0.5, the test sample is judged as a classification; otherwise, it is judged as another class.


### 4.2. Introduction of Probit Model

The binary logistic model of standard RVM is(48)Ptn=1wTϕxn=11+exp−wTϕxn.The logistic function is a mapping from continuous variables to binary outputtn. Logistic mapping function is easy to understand, but it is not a standard probability function; there are many difficulties in the process of reasoning. In addition, the traditional RVM classification model introduces a lower bound to the likelihood function, which is an approximate derivation. Therefore, the real value of the model cannot be estimated accurately. To solve this problem, a mapping method from continuous quantity to discrete quantity through the probit model is defined as follows:(49)Ptn,ynwTϕxn,σ2=PtnynNyn;wTϕxn,σ−2,where the hidden variable ynn=1N is a continuous random variable hidden behind tn. Based on the requirement of the model, the target value tnn=1N is assumed to be −1 or +1. The probability relationship is as follows:(50)Ptnyn=Itn=signyn=1,tn=signyn,0,otherwise,where I⋅ is an indicator function. By integrating yn, it can be found that(51)Ptn,ynwTϕxn,σ2,=∫Ptn=1,ynwTϕxn,σ2dyn,=∫0∞Nyn;wTϕxn,σ−2dyn,=normcdfwTϕxnσ−1,where normcdf⋅ is a normal cumulative distribution function.The advantage of the probit model is to transform the problem of binary-classification output into a regression problem by introducing hidden variables. It makes the models of classification and regression completely equivalent, and the noise variable must be ignored by using the logistic model. It can be seen form Figure6 that the probit model approximates the logistic model very well. Therefore, the reasoning algorithm based on the probit model can be flexibly applied to the classification model directly. In addition, the use of the probit model can also easily extend the binary classification to multiple classifications [43]. Compared with the multivariate logistic model [44], the multivariate probit model can also avoid complex approximate calculations, has more simple and practical characteristics, and can be well approximated to the logistic model.Figure 6
Comparison between the probit model and logistic model (σ2=1).


## 4.1. Variational Relevance Vector Machine

Many problems in machine learning belong to the category of supervised learning, giving an input vectorX=xnn=1N and the corresponding target output T=tnn=1N. For regression problems, they can be arbitrary values. For classification problems, they are class labels. For regression, tn can be any value, and for classification, tn is the category label.For the standard RVM regression model, the formula is defined as follows:(21)tn=yxn;w=∑n=1NwnKx,xn+w0+εn=wTϕxn+εn,where ϕxn=1,Kxn,x1,…,Kxn,xNT, wn is the weight parameter of the model, Kx,xn is the kernel function, and an RBF kernel function is selected in this section: Kx,xn=expx−xn2/h2, h>0. εn represents the noise, which obeys normal distribution N0,σ2.If the input vectorsX=xnn=1N are independent and identically distributed, then the likelihood function of the corresponding target output T=tnn=1N is defined as(22)PTX,w,σ2=∏n=1NPtnxn,w,σ2,where the parameter w is Gaussian prior distribution and σ2is noise variance:(23)Pwα=∏n=1NNwn0,αn−1,where α=αn is a hyperparameter vector, and each weight value wn is independently assigned a parameter αn. In order to make parameter learning more flexible, ultra-prior distribution is defined for α and noise variance σ2, respectively. The appropriate prior distribution is the gamma distribution:(24)Pα=∏n=0NGammaαnan,bn,Pσ2=Gammaσ2c,d,Gammaαa,b=Γa−1baαa−1e−bα,where Γa=∫0∞ta−1e−tdt is the gamma function, and it is usually defined that the hyperparameter is a very small value such as an = bn = c = d = 10−4; such a hyperparameter before does not provide information for posterior learning, so the posterior depends entirely on the data.When the model is established, the posterior distribution ofw, α, and σ2 can be obtained by the variable Bayesian method [40]. In the process of iterative solution, most of αn tends to infinity and the corresponding wn is zero, realizing the sparse model.RVM classification and regression essentially use the same framework model, except that the conditional distribution of the target value is changed. For the binary classification, the logistic functionσy=1/1+e−y is used in the continuous latent variable yxn;w, and assuming that PTX is the Bernoulli distribution, the likelihood function is as follows [41]:(25)PTX,w=∏n=1Nσyxn;wtn1−σyxn;w1−tn.It is important to note that, in the classification problem without considering the noise variableεn, it is very difficult to directly use the variational method to solve the above model. In [28], a lower bound is introduced by using the inequality:(26)σyxn;wtn1−σyxn;w1−tn=σzn≥σξnexpzn−ξn2−λξnzn2−ξn2,where λξn=1/4ξn×tanhξn/2, zn=2tn−1yxn;w, εn is a variational parameter, and when εn=zn, equation (26) is established. Finally, the variational method is used to solve the lower boundary.The VRVM method is the posterior distribution of RVM model parameters and superparameters by variational Bayesian (VB) function. In the VB method [42] based on the RVM model, the observed variables are X=xnn=1N and T=tnn=1N, the hidden variables are ynn=1N, and the parameter is wpm,αp,σ2p=1,m=1PM; therefore, the logarithmic edge likelihood function can be written as(27)lnPTX=ln∬∬PTX,y,w,a,σ2dydwdadσ2=ln∬∬Qy,w,α,σ2PTX,y,w,a,σ2Qy,w,a,σ2dydwdadσ2≥∬∬Qy,w,a,σ2lnPTX,y,w,a,σ2Qy,w,a,σ2dydwdadσ2,where Qy,w,a,σ2 is the joint probability distribution function between the hidden variable and parameter. The variational Bayesian algorithm assumes that y, w, a, and σ2 are independent of each other; therefore, the joint probability distribution of four variables can be written approximately as(28)Qy,w,a,σ2≈QyQwQaQσ2.Following the assumption of equation (28), equation (21) can be rewritten as follows:(29)lnPTX≥∬∬QyQwQaQσ2,lnPTX,w,a,σ2QyQwQαQσ2dydwdαdσ2=FQy,Qw,Qa,Qσ2.It can be seen from equation (29) that the log-likelihood function has a lower bound and that the real value can be approximated by maximizing the lower bound FQy,Qw,Qa,Qσ2. Then, the lower bound of the log-likelihood is solved using EM (expectation-maximization), and the posterior distribution of the hidden variable and all parameters is obtained. The lower bound expression is as follows:(30)FQy,Qw,Qa,Qσ2=∬∬QyQwQaQσ2lnPty+lnPyXw,σ−2+lnPw0,α−1+lnPaa,b+lnPσ2c,d−lnQy−lnQw−lnQa−lnQσ2dydwdadσ2.Because the variable distributions ofQy, Qw, Qa, and Qσ2 are all conjugate prior distributions, they have the same distribution form as their posterior distribution, and the following is their probability distribution:(31)Qy=∏n=1NNtnyn;m˜n,σ˜n,(32)Qw=Nw;μ˜w,Σ˜w,(33)Qσ2=Gammaσ2;c˜,d˜,(34)Qa∏p=1PGammaαp;a˜p,b˜pM,where Ntn⋅ represents truncated normal distribution, and the direction of truncation N+⋅ or N−⋅ is determined by tn.According to the theory of the EM algorithm, the posterior distribution of the variables is actual expectation of the logarithm of the complete likelihood function with respect to other variables (indicated by⋅), and then the terms related to the variable are extracted. The posterior distributions of y, w, a, and σ2 are, respectively, solved as follows:Step 13.
Extract the term related to variabley in equation (30):(35)F˜Qy=∫QylnPty+∬QwQσ2lnPyXw,σ−2dwdσ2−lnQydy.
Set the derivative ofQy in equation (35) to be 0; hence,(36)lnPty+∬QwQσ2lnPyXw,σ−2dwdσ2−lnQy−1=0.
Solve equation (8) and get the following equation:(37)m˜n=xnw,σ˜n=σ2−1Step 14.
Extract the term related to variablew in equation (30):(38)F˜Qw=∫Qw∬QyQσ2lnPyXw,σ−2dydσ2+∫QalnPw0,α−1da−lnQwdw.
Set the derivative ofQw in equation (38) to be 0. Hence,(39)Σ˜w=diaga+σ2∑n=1NxnTxn−1,μ˜w=σ2Σ˜w∑n=1NxnTyn.Step 15.
Extract the term related to variableσ2 in equation (29):(40)F˜Qσ2=∫Qσ2∬QyQwlnPyXw,σ−2dydw+∫QalnPσ2c,d−lnQσ2dσ2.
Set the derivative ofQσ2 in equation (40) to be 0. Hence,(41)c˜=c+0.5,d=d+1N∑n=1N12yn2−2ynwTxnT+xnwwTxnT.Step 16.
Extract the term related to variablea in equation (29):(42)F˜Qa=∫Qa∫QwlnPw0,a−1dw+lnPaa,b−lnQada.
Set the derivative ofQa in equation (42) to be 0. Hence,(43)α˜p=αp+12,b˜p=bp+1M∑n=1N12wpm2.
Through the process above, the parameter iteration formulas in each variable’s posterior distribution function are obtained, but they are expressed by the expectations of other variables. Therefore, the expectation of these parameters needs to be found [42].
For equation (32), tntn∈−1,+1 determines the truncation direction of the truncation normal distribution; hence,(44)yn=xnμ˜w+tn1σ21/2normpdfxnμ˜w/σ2−1/2normpdftnxnμ˜w/σ2−1/2,yn2=1σ2+xnμ˜w2+tnxnμ˜wσ2−1/2normpdfxnμ˜w/σ21/2normpdftnxnμ˜w/σ2−1/2,where normpdf is a normal probability density function Moreover,(45)w=μ˜w,wwT=Σ˜w+μ˜wμ˜wT,αp=a˜pb˜p,lnαp=Ψa˜p−lnb˜p,σ2=c˜d˜,σ2=Ψc˜−lnd˜.
TheΨ function is defined as follows:(46)Ψa=dlnΓada.
After the model is trained, for a test samplex∗, its probability of prediction can be calculated as follows:(47)Pt∗=1x∗=∫Pt∗=1y∗Py∗x∗,w,σ2dy∗=∫0∞Ny∗;x∗w,σ2−1dy∗=normcdfx∗μ˜wσ2−1/2.
wherenormcdf is a normal cumulative integral distribution function. According to the prediction probability, the test sample can be classified and identified. When Pt∗=1x∗≥0.5, the test sample is judged as a classification; otherwise, it is judged as another class.


## 4.2. Introduction of Probit Model

The binary logistic model of standard RVM is(48)Ptn=1wTϕxn=11+exp−wTϕxn.The logistic function is a mapping from continuous variables to binary outputtn. Logistic mapping function is easy to understand, but it is not a standard probability function; there are many difficulties in the process of reasoning. In addition, the traditional RVM classification model introduces a lower bound to the likelihood function, which is an approximate derivation. Therefore, the real value of the model cannot be estimated accurately. To solve this problem, a mapping method from continuous quantity to discrete quantity through the probit model is defined as follows:(49)Ptn,ynwTϕxn,σ2=PtnynNyn;wTϕxn,σ−2,where the hidden variable ynn=1N is a continuous random variable hidden behind tn. Based on the requirement of the model, the target value tnn=1N is assumed to be −1 or +1. The probability relationship is as follows:(50)Ptnyn=Itn=signyn=1,tn=signyn,0,otherwise,where I⋅ is an indicator function. By integrating yn, it can be found that(51)Ptn,ynwTϕxn,σ2,=∫Ptn=1,ynwTϕxn,σ2dyn,=∫0∞Nyn;wTϕxn,σ−2dyn,=normcdfwTϕxnσ−1,where normcdf⋅ is a normal cumulative distribution function.The advantage of the probit model is to transform the problem of binary-classification output into a regression problem by introducing hidden variables. It makes the models of classification and regression completely equivalent, and the noise variable must be ignored by using the logistic model. It can be seen form Figure6 that the probit model approximates the logistic model very well. Therefore, the reasoning algorithm based on the probit model can be flexibly applied to the classification model directly. In addition, the use of the probit model can also easily extend the binary classification to multiple classifications [43]. Compared with the multivariate logistic model [44], the multivariate probit model can also avoid complex approximate calculations, has more simple and practical characteristics, and can be well approximated to the logistic model.Figure 6
Comparison between the probit model and logistic model (σ2=1).


## 5. Analysis of Bearing Fault Diagnosis

### 5.1. Experimental Scheme

In order to further verify the effectiveness of the proposed method for bearing fault diagnosis under unknown variable load conditions, the experiments are conducted to collect vibration signals from bearings operating under variable load conditions. The experimental platform is shown in Figure7; it uses bearing type 6203-2RS JEM SKF deep groove rolling bearing. The bearing inner diameter is 25 mm, outer diameter is 52 mm, and thickness is 15 mm; the rolling diameter is 8.18 mm, and the pitch diameter of the bearing is 44.2 mm; the sampling frequency is 12 kHz, and the sampling data length is 12000. The fault with different etch diameters (simulating varying degrees of damage to the inner ring, the outer ring, and the rolling body) is processed by electric spark. Bearing loads include 3 types: load 0 (0 N, 1797 r/min, simulated no load), load 1 (800 N, 1772 r/min, simulated light load), and load 2 (1600 N, 1750 r/m, simulated heavy load). The pitting diameter of the bearing is 0.1778 mm (it simulates minor fault), 0.3556 mm (it simulates medium fault), and 0.5334 mm (it simulates serious fault), which is used to simulate three different damage degrees of the bearing.Figure 7
Experimental platform.


### 5.2. Fault Feature Extraction

To further verify the effectiveness of the VMD method based on the quantum chaotic FOA in the early bearing fault feature extraction, this paper divides the fault feature extraction into two cases: the fault characteristics of the minor faults under different loads and the fault characteristics of different damage faults under the same load.


#### 5.2.1. Fault Feature Extraction of Weak Damage Degree under Different Loads

Under three different loads, the measured signals of small bearing failure (pitting diameter 0.1778 mm) and small features of the rolling body fault are used as experimental signals. The fault signals are shown in Figure8. The parameter α,Kof the VMD method optimized by the quantum chaotic FOA is used to search the optimal value, and the search results are shown in Table 2.Figure 8
Rolling body minor fault signals under three different loads.Table 2
Optimal search results.
Load
0 hp
1 hp
2 hp
The calculated value of the fault frequency (Hz)
119.39
118.01
116.61
Fault measured value (Hz)
119.4
117.9
116.5
[α, K]
[2000, 6]
[1600, 6]
[1400, 8]
Optimal IMF
IMF4
IMF4
IMF6
Normalized center frequency
0.4947
0.4224
0.3644
LMMSE
0.3228
0.3231
0.3275The optimal IMF component and its corresponding marginal spectrum and center frequency for the small faults of rolling bodies under three different loads are shown in Figures9(a)–9(c), respectively. fr represents the rotor frequency, and fball represents the fault characteristic frequency of the rolling element. It can be seen from Figure 9 that the weak fault characteristic frequencies can be effectively extracted under the three loads, the frequency amplitude decreases with the increase of the load, and the bearing fault feature frequency amplitude under load 2 hp is minimum. In addition, it can be seen from Table 2 that the center frequency of the slight fault of the rolling body is also different under different loads. To further validate the effectiveness of VMD based on the QCFOA, the EEMD and LMD methods are used to deal with the rolling body fault signal with the pitting diameter of 0.1778 mm under load 2 hp, respectively. It can be seen from Figure 9 that the characteristic frequency of the bearing fault obtained by the EEMD method is very weak and almost drowned by other frequency components. The LMD method has improved relative to the EEMD method, but the fault characteristic frequency is still very weak. Compared with the VMD method in Figure 10, the characteristic frequency of the bearing fault is superior than that of EEMD and LMD. Therefore, the above experimental results indicate that the VMD method using the quantum chaotic FOA can accurately extract the characteristic frequency of the weak rolling body fault under heavy loads, and the validity of the proposed method is also verified.Figure 9
The best IMF component and its corresponding marginal spectrum and center frequency. Marginal spectrum of IMF4 components under load: (a) 0 hp; (b) 1 hp; (c) 2 hp.
(a)
(b)
(c)Figure 10
Marginal spectrum of minor rolling body fault under heavy load conditions. (a) EEMD decomposition. (b) LMD decomposition.
(a)
(b)


#### 5.2.2. Fault Feature Extraction of Different Defects under the Same Load

The experiment is conducted on a bearing with different degrees of fault, and the collected vibration signals shown in Figure9 are of a bearing with different defects.It can be seen from Figure11 that the frequency of fault features with different damage degrees under the same load can be effectively extracted, the amplitude of the frequency increases with the increase of the degree of rolling damage, and the frequency amplitude of the minor damage fault features is the least. The center frequencies of rolling body fault with different damage degrees used for the experiment are given in Table 3, and the changes are minor.Figure 11
Vibration signals under different pitting radius.Table 3
Experimental data under different fault degrees.
Damage degree
Minor
Moderate
Severe
Fault frequency calculated value (Hz)
116.61
116.61
116.61
Measured value (Hz)
116.5
116.5
116.5
[α, K]
[1800, 6]
[1600, 6]
[1400, 8]
Optimal IMF
IMF4
IMF4
IMF6
Normalization center frequency
0.4586
0.4624
0.4667
LMMSE
0.4875
0.3994
0.3211


### 5.3. Selection of Eigenvectors

The marginal spectrum can accurately reflect the distribution of the actual frequency components of the signal and the degree of uncertainty of the signal spectrum. The smaller the marginal spectral entropy of the signal is, the greater the concentration of the energy spectrum of the signal is and the more concentrated the marginal spectrum of the signal is. The analysis of the marginal spectrum energy can effectively reflect the working state of the bearing. From Table2, it is seen that there is a small difference in the LMMSE of the same fault type and the fault level at different loads, but the difference in their center frequency values is large. It can be seen from Table 3 that the LMMSE under the same load is significantly different for the same fault type and different fault degrees, but the difference between their central frequencies is small. The experimental data show that the magnitude of load and the degree of fault have a certain effect on the marginal spectrum entropy of failure. Therefore, the marginal spectrum entropy of a single IMF component is difficult to accurately characterize the degree of fault under variable load conditions. For the above reasons, the center frequency is combined with the basis of one-dimensional MSE to extend it into two-dimensional MSE in this paper, and its central frequency can be obtained through equation (2).In order to verify the effectiveness and robustness of the proposed method under variable load conditions, the experimental data under 1 hp load are used as a training sample, with 0 hp and 2 hp representing the unknown load, and the experimental data are used as a test set. The length of each data sample is 4096 points. The detailed data description is shown in Table4.Table 4
Experimental data under variable load conditions.
Training samples (1 hp)
Test samples (0 hp + 2 hp)
Pitting radius (mm)
Bearing state
Fault identifier
40
40
0
Normal
0
80
80
0.1778
Outer race
1
80
80
0.3556
Outer race
2
80
80
0.5334
Outer race
3
80
80
0.1778
Rolling element
4
80
80
0.3556
Rolling element
5
80
80
0.5334
Rolling element
6
80
80
0.1778
Inner race
7
80
80
0.3556
Inner race
8
80
80
0.5334
Inner race
9


### 5.4. Fault Diagnosis

#### 5.4.1. The Process of Intelligent Diagnosis

First, the sample data are decomposed by VMD, and the corresponding IMF components are obtained. The two-dimensional MSE synthesized by marginal spectrum entropy and the central frequency corresponding to each IMF component is taken as the eigenvector. Then, the most commonly used RBF kernel function is selected as the kernel function of VRVM. The parameters of the kernel function are optimized through 10-fold cross-validation. The detailed diagnosis process is shown in Figure12.Figure 12
Variable load fault diagnosis process based on VRVM.In Figure12, three multiclassifiers of VRVM are used. Although the probit model is proposed in this paper, it can easily extend the binary classification to the multiple classification and avoid complex approximate calculation. However, as the number of classes increases, the Hessian matrix in the process of constructing the model will also increase, resulting in an increase in computational complexity. Therefore, in this paper, three VRVM multiclassifiers are constructed and a multiclassification intelligent diagnosis model is constructed based on the combination strategy of “maximum probability win”.


#### 5.4.2. Diagnostic Results

In order to verify the validity and robustness of the proposed method in intelligent fault diagnosis under variable load conditions, the fault data of different damage degrees under the same load and the fault data with slight damage under different radial loads are tested for the diagnosis model in this section. The three different levels of damage are as follows: (1) minor damage, (2) moderate damage, and (3) severe damage. The three different radial loads are 0 hp, 1 hp, and 2 hp, respectively. 0 hp represents the radial load force of 0 N, 1 hp represents the radial load force of 400 N, and 2 hp represents the radial load force of 800 N.According to the proposed bearing fault intelligent diagnosis procedure, the diagnosis results of bearing outer race damage with different degrees under 0 hp load are shown in Figure13(a). It can clearly distinguish the running data of bearing outer race under normal and three different damage degree conditions. Similarly, the diagnosis results of bearing rolling balls under 1 ph load and bearing inner race under 2 hp load are shown in Figures 13(b)–13(c), respectively. Under different loading conditions, the fault diagnosis results of the bearing inner ring, rolling body, and outer ring with slight damage are shown in Figure 14.Figure 13
Intelligent diagnosis of different bearing faults under variable load. Diagnosis results of three different damage degrees of the bearing outer race under 0 hp load (a), bearing rolling balls under 1 hp load (b), and bearing inner race under 2 hp load (c).
(a)
(b)
(c)Figure 14
Intelligent diagnosis of bearing faults with minor damage under three different load conditions. (a) Bearing outer race fault. (b) Bearing rolling ball fault. (c) Bearing inner race fault.
(a)
(b)
(c)According to the proposed procedure shown in Figures13 and 14, fault types and degrees can be effectively identified under variable loads with the intelligent diagnosis method. The fault recognition rate is 100%, and the fault degree recognition rate also reaches 90%; the overall recognition rate is 95%, which proves the validity of the method. Meanwhile, VRVM is compared with other classification methods in terms of classification accuracy and running time. The results are shown in Table 5. The experiment results show that the VRVM method has a good performance in classification accuracy, but the running time is longer. In practice, equipment fault samples are often scarce, and it is even more difficult to get data of different fault types under the same load. It is a practical significance to identify fault types and fault levels under unknown loads by using fault data under known loads.Table 5
Performance comparison of different classification methods.
Algorithms
Classification accuracy (%)
Run times (s)
Testing precision
Training precision
LSSVM
71.0
73.1
27.1
FSVM
80.6
82.1
24.1
KFCM-FSVM
84.3
87.8
36.2
RVM
82.4
84.7
34.8
M-RVM
83.2
88.2
38.7
V-RVM
89.1
95.0
31.6


## 5.1. Experimental Scheme

In order to further verify the effectiveness of the proposed method for bearing fault diagnosis under unknown variable load conditions, the experiments are conducted to collect vibration signals from bearings operating under variable load conditions. The experimental platform is shown in Figure7; it uses bearing type 6203-2RS JEM SKF deep groove rolling bearing. The bearing inner diameter is 25 mm, outer diameter is 52 mm, and thickness is 15 mm; the rolling diameter is 8.18 mm, and the pitch diameter of the bearing is 44.2 mm; the sampling frequency is 12 kHz, and the sampling data length is 12000. The fault with different etch diameters (simulating varying degrees of damage to the inner ring, the outer ring, and the rolling body) is processed by electric spark. Bearing loads include 3 types: load 0 (0 N, 1797 r/min, simulated no load), load 1 (800 N, 1772 r/min, simulated light load), and load 2 (1600 N, 1750 r/m, simulated heavy load). The pitting diameter of the bearing is 0.1778 mm (it simulates minor fault), 0.3556 mm (it simulates medium fault), and 0.5334 mm (it simulates serious fault), which is used to simulate three different damage degrees of the bearing.Figure 7
Experimental platform.


## 5.2. Fault Feature Extraction

To further verify the effectiveness of the VMD method based on the quantum chaotic FOA in the early bearing fault feature extraction, this paper divides the fault feature extraction into two cases: the fault characteristics of the minor faults under different loads and the fault characteristics of different damage faults under the same load.


### 5.2.1. Fault Feature Extraction of Weak Damage Degree under Different Loads

Under three different loads, the measured signals of small bearing failure (pitting diameter 0.1778 mm) and small features of the rolling body fault are used as experimental signals. The fault signals are shown in Figure8. The parameter α,Kof the VMD method optimized by the quantum chaotic FOA is used to search the optimal value, and the search results are shown in Table 2.Figure 8
Rolling body minor fault signals under three different loads.Table 2
Optimal search results.
Load
0 hp
1 hp
2 hp
The calculated value of the fault frequency (Hz)
119.39
118.01
116.61
Fault measured value (Hz)
119.4
117.9
116.5
[α, K]
[2000, 6]
[1600, 6]
[1400, 8]
Optimal IMF
IMF4
IMF4
IMF6
Normalized center frequency
0.4947
0.4224
0.3644
LMMSE
0.3228
0.3231
0.3275The optimal IMF component and its corresponding marginal spectrum and center frequency for the small faults of rolling bodies under three different loads are shown in Figures9(a)–9(c), respectively. fr represents the rotor frequency, and fball represents the fault characteristic frequency of the rolling element. It can be seen from Figure 9 that the weak fault characteristic frequencies can be effectively extracted under the three loads, the frequency amplitude decreases with the increase of the load, and the bearing fault feature frequency amplitude under load 2 hp is minimum. In addition, it can be seen from Table 2 that the center frequency of the slight fault of the rolling body is also different under different loads. To further validate the effectiveness of VMD based on the QCFOA, the EEMD and LMD methods are used to deal with the rolling body fault signal with the pitting diameter of 0.1778 mm under load 2 hp, respectively. It can be seen from Figure 9 that the characteristic frequency of the bearing fault obtained by the EEMD method is very weak and almost drowned by other frequency components. The LMD method has improved relative to the EEMD method, but the fault characteristic frequency is still very weak. Compared with the VMD method in Figure 10, the characteristic frequency of the bearing fault is superior than that of EEMD and LMD. Therefore, the above experimental results indicate that the VMD method using the quantum chaotic FOA can accurately extract the characteristic frequency of the weak rolling body fault under heavy loads, and the validity of the proposed method is also verified.Figure 9
The best IMF component and its corresponding marginal spectrum and center frequency. Marginal spectrum of IMF4 components under load: (a) 0 hp; (b) 1 hp; (c) 2 hp.
(a)
(b)
(c)Figure 10
Marginal spectrum of minor rolling body fault under heavy load conditions. (a) EEMD decomposition. (b) LMD decomposition.
(a)
(b)


### 5.2.2. Fault Feature Extraction of Different Defects under the Same Load

The experiment is conducted on a bearing with different degrees of fault, and the collected vibration signals shown in Figure9 are of a bearing with different defects.It can be seen from Figure11 that the frequency of fault features with different damage degrees under the same load can be effectively extracted, the amplitude of the frequency increases with the increase of the degree of rolling damage, and the frequency amplitude of the minor damage fault features is the least. The center frequencies of rolling body fault with different damage degrees used for the experiment are given in Table 3, and the changes are minor.Figure 11
Vibration signals under different pitting radius.Table 3
Experimental data under different fault degrees.
Damage degree
Minor
Moderate
Severe
Fault frequency calculated value (Hz)
116.61
116.61
116.61
Measured value (Hz)
116.5
116.5
116.5
[α, K]
[1800, 6]
[1600, 6]
[1400, 8]
Optimal IMF
IMF4
IMF4
IMF6
Normalization center frequency
0.4586
0.4624
0.4667
LMMSE
0.4875
0.3994
0.3211


## 5.2.1. Fault Feature Extraction of Weak Damage Degree under Different Loads

Under three different loads, the measured signals of small bearing failure (pitting diameter 0.1778 mm) and small features of the rolling body fault are used as experimental signals. The fault signals are shown in Figure8. The parameter α,Kof the VMD method optimized by the quantum chaotic FOA is used to search the optimal value, and the search results are shown in Table 2.Figure 8
Rolling body minor fault signals under three different loads.Table 2
Optimal search results.
Load
0 hp
1 hp
2 hp
The calculated value of the fault frequency (Hz)
119.39
118.01
116.61
Fault measured value (Hz)
119.4
117.9
116.5
[α, K]
[2000, 6]
[1600, 6]
[1400, 8]
Optimal IMF
IMF4
IMF4
IMF6
Normalized center frequency
0.4947
0.4224
0.3644
LMMSE
0.3228
0.3231
0.3275The optimal IMF component and its corresponding marginal spectrum and center frequency for the small faults of rolling bodies under three different loads are shown in Figures9(a)–9(c), respectively. fr represents the rotor frequency, and fball represents the fault characteristic frequency of the rolling element. It can be seen from Figure 9 that the weak fault characteristic frequencies can be effectively extracted under the three loads, the frequency amplitude decreases with the increase of the load, and the bearing fault feature frequency amplitude under load 2 hp is minimum. In addition, it can be seen from Table 2 that the center frequency of the slight fault of the rolling body is also different under different loads. To further validate the effectiveness of VMD based on the QCFOA, the EEMD and LMD methods are used to deal with the rolling body fault signal with the pitting diameter of 0.1778 mm under load 2 hp, respectively. It can be seen from Figure 9 that the characteristic frequency of the bearing fault obtained by the EEMD method is very weak and almost drowned by other frequency components. The LMD method has improved relative to the EEMD method, but the fault characteristic frequency is still very weak. Compared with the VMD method in Figure 10, the characteristic frequency of the bearing fault is superior than that of EEMD and LMD. Therefore, the above experimental results indicate that the VMD method using the quantum chaotic FOA can accurately extract the characteristic frequency of the weak rolling body fault under heavy loads, and the validity of the proposed method is also verified.Figure 9
The best IMF component and its corresponding marginal spectrum and center frequency. Marginal spectrum of IMF4 components under load: (a) 0 hp; (b) 1 hp; (c) 2 hp.
(a)
(b)
(c)Figure 10
Marginal spectrum of minor rolling body fault under heavy load conditions. (a) EEMD decomposition. (b) LMD decomposition.
(a)
(b)


## 5.2.2. Fault Feature Extraction of Different Defects under the Same Load

The experiment is conducted on a bearing with different degrees of fault, and the collected vibration signals shown in Figure9 are of a bearing with different defects.It can be seen from Figure11 that the frequency of fault features with different damage degrees under the same load can be effectively extracted, the amplitude of the frequency increases with the increase of the degree of rolling damage, and the frequency amplitude of the minor damage fault features is the least. The center frequencies of rolling body fault with different damage degrees used for the experiment are given in Table 3, and the changes are minor.Figure 11
Vibration signals under different pitting radius.Table 3
Experimental data under different fault degrees.
Damage degree
Minor
Moderate
Severe
Fault frequency calculated value (Hz)
116.61
116.61
116.61
Measured value (Hz)
116.5
116.5
116.5
[α, K]
[1800, 6]
[1600, 6]
[1400, 8]
Optimal IMF
IMF4
IMF4
IMF6
Normalization center frequency
0.4586
0.4624
0.4667
LMMSE
0.4875
0.3994
0.3211


## 5.3. Selection of Eigenvectors

The marginal spectrum can accurately reflect the distribution of the actual frequency components of the signal and the degree of uncertainty of the signal spectrum. The smaller the marginal spectral entropy of the signal is, the greater the concentration of the energy spectrum of the signal is and the more concentrated the marginal spectrum of the signal is. The analysis of the marginal spectrum energy can effectively reflect the working state of the bearing. From Table2, it is seen that there is a small difference in the LMMSE of the same fault type and the fault level at different loads, but the difference in their center frequency values is large. It can be seen from Table 3 that the LMMSE under the same load is significantly different for the same fault type and different fault degrees, but the difference between their central frequencies is small. The experimental data show that the magnitude of load and the degree of fault have a certain effect on the marginal spectrum entropy of failure. Therefore, the marginal spectrum entropy of a single IMF component is difficult to accurately characterize the degree of fault under variable load conditions. For the above reasons, the center frequency is combined with the basis of one-dimensional MSE to extend it into two-dimensional MSE in this paper, and its central frequency can be obtained through equation (2).In order to verify the effectiveness and robustness of the proposed method under variable load conditions, the experimental data under 1 hp load are used as a training sample, with 0 hp and 2 hp representing the unknown load, and the experimental data are used as a test set. The length of each data sample is 4096 points. The detailed data description is shown in Table4.Table 4
Experimental data under variable load conditions.
Training samples (1 hp)
Test samples (0 hp + 2 hp)
Pitting radius (mm)
Bearing state
Fault identifier
40
40
0
Normal
0
80
80
0.1778
Outer race
1
80
80
0.3556
Outer race
2
80
80
0.5334
Outer race
3
80
80
0.1778
Rolling element
4
80
80
0.3556
Rolling element
5
80
80
0.5334
Rolling element
6
80
80
0.1778
Inner race
7
80
80
0.3556
Inner race
8
80
80
0.5334
Inner race
9


## 5.4. Fault Diagnosis

### 5.4.1. The Process of Intelligent Diagnosis

First, the sample data are decomposed by VMD, and the corresponding IMF components are obtained. The two-dimensional MSE synthesized by marginal spectrum entropy and the central frequency corresponding to each IMF component is taken as the eigenvector. Then, the most commonly used RBF kernel function is selected as the kernel function of VRVM. The parameters of the kernel function are optimized through 10-fold cross-validation. The detailed diagnosis process is shown in Figure12.Figure 12
Variable load fault diagnosis process based on VRVM.In Figure12, three multiclassifiers of VRVM are used. Although the probit model is proposed in this paper, it can easily extend the binary classification to the multiple classification and avoid complex approximate calculation. However, as the number of classes increases, the Hessian matrix in the process of constructing the model will also increase, resulting in an increase in computational complexity. Therefore, in this paper, three VRVM multiclassifiers are constructed and a multiclassification intelligent diagnosis model is constructed based on the combination strategy of “maximum probability win”.


### 5.4.2. Diagnostic Results

In order to verify the validity and robustness of the proposed method in intelligent fault diagnosis under variable load conditions, the fault data of different damage degrees under the same load and the fault data with slight damage under different radial loads are tested for the diagnosis model in this section. The three different levels of damage are as follows: (1) minor damage, (2) moderate damage, and (3) severe damage. The three different radial loads are 0 hp, 1 hp, and 2 hp, respectively. 0 hp represents the radial load force of 0 N, 1 hp represents the radial load force of 400 N, and 2 hp represents the radial load force of 800 N.According to the proposed bearing fault intelligent diagnosis procedure, the diagnosis results of bearing outer race damage with different degrees under 0 hp load are shown in Figure13(a). It can clearly distinguish the running data of bearing outer race under normal and three different damage degree conditions. Similarly, the diagnosis results of bearing rolling balls under 1 ph load and bearing inner race under 2 hp load are shown in Figures 13(b)–13(c), respectively. Under different loading conditions, the fault diagnosis results of the bearing inner ring, rolling body, and outer ring with slight damage are shown in Figure 14.Figure 13
Intelligent diagnosis of different bearing faults under variable load. Diagnosis results of three different damage degrees of the bearing outer race under 0 hp load (a), bearing rolling balls under 1 hp load (b), and bearing inner race under 2 hp load (c).
(a)
(b)
(c)Figure 14
Intelligent diagnosis of bearing faults with minor damage under three different load conditions. (a) Bearing outer race fault. (b) Bearing rolling ball fault. (c) Bearing inner race fault.
(a)
(b)
(c)According to the proposed procedure shown in Figures13 and 14, fault types and degrees can be effectively identified under variable loads with the intelligent diagnosis method. The fault recognition rate is 100%, and the fault degree recognition rate also reaches 90%; the overall recognition rate is 95%, which proves the validity of the method. Meanwhile, VRVM is compared with other classification methods in terms of classification accuracy and running time. The results are shown in Table 5. The experiment results show that the VRVM method has a good performance in classification accuracy, but the running time is longer. In practice, equipment fault samples are often scarce, and it is even more difficult to get data of different fault types under the same load. It is a practical significance to identify fault types and fault levels under unknown loads by using fault data under known loads.Table 5
Performance comparison of different classification methods.
Algorithms
Classification accuracy (%)
Run times (s)
Testing precision
Training precision
LSSVM
71.0
73.1
27.1
FSVM
80.6
82.1
24.1
KFCM-FSVM
84.3
87.8
36.2
RVM
82.4
84.7
34.8
M-RVM
83.2
88.2
38.7
V-RVM
89.1
95.0
31.6


## 5.4.1. The Process of Intelligent Diagnosis

First, the sample data are decomposed by VMD, and the corresponding IMF components are obtained. The two-dimensional MSE synthesized by marginal spectrum entropy and the central frequency corresponding to each IMF component is taken as the eigenvector. Then, the most commonly used RBF kernel function is selected as the kernel function of VRVM. The parameters of the kernel function are optimized through 10-fold cross-validation. The detailed diagnosis process is shown in Figure12.Figure 12
Variable load fault diagnosis process based on VRVM.In Figure12, three multiclassifiers of VRVM are used. Although the probit model is proposed in this paper, it can easily extend the binary classification to the multiple classification and avoid complex approximate calculation. However, as the number of classes increases, the Hessian matrix in the process of constructing the model will also increase, resulting in an increase in computational complexity. Therefore, in this paper, three VRVM multiclassifiers are constructed and a multiclassification intelligent diagnosis model is constructed based on the combination strategy of “maximum probability win”.


## 5.4.2. Diagnostic Results

In order to verify the validity and robustness of the proposed method in intelligent fault diagnosis under variable load conditions, the fault data of different damage degrees under the same load and the fault data with slight damage under different radial loads are tested for the diagnosis model in this section. The three different levels of damage are as follows: (1) minor damage, (2) moderate damage, and (3) severe damage. The three different radial loads are 0 hp, 1 hp, and 2 hp, respectively. 0 hp represents the radial load force of 0 N, 1 hp represents the radial load force of 400 N, and 2 hp represents the radial load force of 800 N.According to the proposed bearing fault intelligent diagnosis procedure, the diagnosis results of bearing outer race damage with different degrees under 0 hp load are shown in Figure13(a). It can clearly distinguish the running data of bearing outer race under normal and three different damage degree conditions. Similarly, the diagnosis results of bearing rolling balls under 1 ph load and bearing inner race under 2 hp load are shown in Figures 13(b)–13(c), respectively. Under different loading conditions, the fault diagnosis results of the bearing inner ring, rolling body, and outer ring with slight damage are shown in Figure 14.Figure 13
Intelligent diagnosis of different bearing faults under variable load. Diagnosis results of three different damage degrees of the bearing outer race under 0 hp load (a), bearing rolling balls under 1 hp load (b), and bearing inner race under 2 hp load (c).
(a)
(b)
(c)Figure 14
Intelligent diagnosis of bearing faults with minor damage under three different load conditions. (a) Bearing outer race fault. (b) Bearing rolling ball fault. (c) Bearing inner race fault.
(a)
(b)
(c)According to the proposed procedure shown in Figures13 and 14, fault types and degrees can be effectively identified under variable loads with the intelligent diagnosis method. The fault recognition rate is 100%, and the fault degree recognition rate also reaches 90%; the overall recognition rate is 95%, which proves the validity of the method. Meanwhile, VRVM is compared with other classification methods in terms of classification accuracy and running time. The results are shown in Table 5. The experiment results show that the VRVM method has a good performance in classification accuracy, but the running time is longer. In practice, equipment fault samples are often scarce, and it is even more difficult to get data of different fault types under the same load. It is a practical significance to identify fault types and fault levels under unknown loads by using fault data under known loads.Table 5
Performance comparison of different classification methods.
Algorithms
Classification accuracy (%)
Run times (s)
Testing precision
Training precision
LSSVM
71.0
73.1
27.1
FSVM
80.6
82.1
24.1
KFCM-FSVM
84.3
87.8
36.2
RVM
82.4
84.7
34.8
M-RVM
83.2
88.2
38.7
V-RVM
89.1
95.0
31.6


## 6. Conclusions

In this paper, the optimal variational mode decomposition (VMD) based on the quantum chaotic fruit fly optimization algorithm and variational relevance vector machine (VRVM) are combined as a hybrid method to diagnose the bearings’ fault under variable load conditions. The results of the experiment and application demonstrate the superiority of the proposed method. The conclusions are summarized as follows:(1)
The VMD method is a new adaptive signal processing method. In the process of bearing fault signal processing, its performance is influenced by the two parameters such as the number of its own components and the penalty factor. Therefore, using the quantum chaotic fruit fly optimization algorithm to filter the two key parameters of the optimal value can guarantee the effectiveness and reliability of VMD performance.(2)
FOA is a new swarm intelligence optimization algorithm. Its principle is simple and easy to implement and has strong local search capability; however, its global search is weak. If the initial value is not set properly, it will easily fall into local minimum, thus losing population diversity and premature convergence. Compared with the traditional chaotic system, the quantum chaotic system has the characteristics of better aperiodicity, ergodicity, and class randomness and more sensitivity to system parameters and initial conditions. Therefore, the location of the fruit fly population is initialized by the proposed quantum chaotic system, which can improve the diversity of the population, strengthen the ergodicity of the search, avoid the search process falling into the local optimal value prematurely, and improve the search efficiency.(3)
The posterior distribution of all parameters and superparameters is obtained by using RVM. Then, the probit model is used to replace the logistic model in the original variational RVM classification so that the classification and regression are organically combined. The approximate deduction of the logistic model from continuous output to discrete output mapping is avoided so that the reasoning algorithm of the RVM regression model can be directly applied to the classification model.(4)
Usually, the fault diagnosis model is trained by the traditional diagnosis method under a specific load, and it has a great limitation. In practice, most mechanical equipment work under variable load conditions. A single marginal spectrum entropy cannot effectively characterize the degree of the fault under variable load conditions. Therefore, by introducing the central frequency of VMD, the one-dimensional marginal spectrum entropy is extended to two-dimensional marginal spectrum entropy as the learning sample of VRVM.(5)
To develop an intelligent fault diagnosis model which is a systematic problem, it includes experimental condition design, feature extraction, feature selection, and model training, and each step will affect the validity of the final model. The method proposed in this paper provides an effective diagnostic strategy for multiclass faults and fault degree diagnosis under variable load conditions. However, identification of the accuracy and running time of bearing damage recognition needs to be improved in this paper, which is also a problem that needs further study in the later stage.



---

*Source: 8213056-2019-01-03.xml*