# Rolling Bearing Fault Diagnosis Based on MFDFA-SPS and ELM

Authors: Yunfan Yang; Caiping Xi
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

Rolling bearings, as important parts on supporting rotating shafts, frequently suffer from fatigue failures. If these rolling bearing failures are not found in time, it will have a huge impact on the whole mechanical system’s operating safety and operating life. To improve the diagnosis of different faults as well as different degrees of faults, a fault diagnosis method based on the multifractal detrended fluctuation analysis (MFDFA) method-singularity power spectrum (SPS) with extreme learning machine (ELM) is proposed. First, MFDFA and SPS analyses are performed on vibration acceleration signals with different faults and different degrees of damage under the same operating conditions, the spectral parameters of stability and quantitative description of differentiation are selected for feature extraction, and then the selected six feature parameters are put into the extreme learning machine for fault classification. The effectiveness of the MFDFA-SPS feature extraction method is demonstrated by analyzing and testing the measured bearing signals. The fault diagnosis accuracy of the bearing fault signals can reach 99.2% based on the MFDFA-SPS with ELM method by using the Case Western Reserve database. The improvements are 6.79% and 18.42% compared to the fault diagnosis methods based on MFDFA with ELM and SPS with ELM. Compared with the methods based on MFDFA-SPS with LSSVM classifier and SVM classifier, the accuracy improvements are 3.54% and 4.25%, respectively. The results show that the method proposed in this paper can achieve the diagnosis of bearing faults and the method based on MFDFA-SPS with ELM is more efficient than the methods based on MFDFA-SPS with LSSVM and SVM classifiers, which is suitable for practical engineering problem-solving.


---

## Body

## 1. Introduction

Rolling bearing is one of the important components of the mechanical transmission system, due to corrosion and fatigue factors, and easy to cause bearing abrasion, pitting, and cracking failures [1]. Different faults cause different vibration noise and increased rotational resistance, making the whole transmission system life decay or even fail [2]. Therefore, it is important to ensure the normal operation of the bearings in the gearbox or drive shaft and the timely diagnosis and analysis of the structural components after a failure occurs for the normal operation of the whole unit.Since the concept of “fractal” was first introduced by Mandelbrot in the 1980s [3], this theory of studying nonlinear systems has been widely applied in many fields and has achieved fruitful results. The object of study of fractal theory is not smooth and irregular geometry in nature, and the bearing fault diagnosis method based on fractal theory is one of the more suitable and effective methods for bearing fault identification. Commonly used fractal theory basis includes iterative function system (IFS), fractal Brownian random field, single fractal, extended fractal, multifractal, and multifractal correlation [4, 5]. The difference in internal signal characteristics between faulty vibration signals and normal operation vibration signals is usually used in research for fault diagnosis by extracting fractal geometric change features.For decades, scholars worldwide have proposed many methods for fault diagnosis of bearings, and the current mainstream diagnosis methods are the temperature analysis method [6] and vibration analysis method [7]. The enriched fault features such as cyclic power spectra [8] and vibration cyclostationarity [9] parameters can be extracted when using bearing vibration signals for fault monitoring, which enables excellent diagnostic performance. The vibration analysis method mainly performs digital signal acquisition of the faulty bearing vibration signal, and its acquisition method is low cost and contains rich fault information in the vibration signal, which is convenient for real-time diagnosis and analysis.For the analysis of nonlinear signal characteristics from the perspective of multifractals, Shao et al. [10] used the fractal box dimension method to quantitatively characterize the vibration signals of faulty bearings and found that the fractal dimensions of the faults differed from each other in their manifestations. Liu et al. [11] used detrended fluctuation analysis (DFA) to extract features of wind turbine bearing faults and proved that the fault characteristic parameters can be obtained by multifractal detrended fluctuation analysis (MFDFA). Lin and Chen [12] used the shape and location feature parameters of the multifractal spectrum obtained by the detrended fluctuation analysis and combined with the martingale repulsion to achieve the identification of the bearing damage degree. Chen et al. [13] proposed a fault prediction method based on fractal dimension using an extreme learning machine (ELM). Xiong et al. [14] used multifractal detrended fluctuation analysis methods with PSO optimized LSSVM classification method to achieve fault diagnosis of bearings. Numerous worldwide research papers have demonstrated that fractal theory can be applied to faulty bearing vibration signal processing, and the difference lies in the different effects of using different feature parameters’ extraction criteria and classifiers. At present, the classifier based on combining multifractal theory and deep learning is the current niche use diagnosis method, and the diagnosis effect is fast and has a good diagnosis rate with further development trends. ELM and deep learning are complementary, and some applications combine the two to get promising results, such as using CNNs for feature extraction [15] and ELMs for classifiers.In summary, to further realize the applicability of fractal theory to bearing fault analysis and diagnosis, this paper proposes a rolling bearing composite fault method based on MFDFA-singularity power spectrum (SPS) and ELM. The method is a quantitative description of nonlinear signals using multifractals for feature extraction. ELM is applied for further exploration of multiple classifications, and the validity of the method is verified by the actual measurement of many different types of faulty bearing signals.This paper is organized as follows: Section2 describes the specific model of the MFDFA algorithm; Section 3 describes the specific steps of the SPS algorithm; Section 4 presents a brief description of the ELM tool and the way to extract features using the algorithms in Sections 2 and 3; Section 5 describes the source of data, the specific feature parameter distribution, and the results of the experiments; the main findings and conclusions of the article are summarized in Section 6.


## 2. The Multifractal Detrended Fluctuation Analysis

The MFDFA algorithm implementation steps are described as follows [12]:Step 1: for a times seriesxt of length N, t=1,2,…,N, construct a new sequence yt:(1)yt=∑i=1txi−1N∑j=1Nxj,t=1,2,…,N.Step 2: partition the sequenceyt into NS=N/S mutually disjoint equal-length subintervals. There are s data points in each interval, maxm+2,10≤s≤N/10. Since the sequence length N is not necessarily an integer multiple of s, then yt will have segments remaining, and to make the sequence of the last small segment processed, the partitioning process can be repeated from the tail of the sequence yt, thus obtaining 2Ns equal-length subintervals. m represents the order of polynomial fit, generally taken as m=1,2,3,…. Here, MFDFA can be denoted as MFDFA-m, and when m takes the value of 1, the method is denoted as MFDFA-1.Themth-order polynomial yvmt in the vth subinterval v=1,2,…,2Ns is fitted by the least-squares method, and the residual corresponding to the label t in this interval is calculated as follows:(2)y˜vt=yt−yvmt,t=v−1s+i,i=1,2,...,s.Step 3: calculate the mean square error separately:(3)Fvs=1s∑i=1sy˜v2v−1s+i1/2,v=1,2,3,…,Ns,Fvs=1s∑i=1sy˜v2N−v−Nss+i1/2,v=Ns+1,…,2Ns.Step 4: calculate the global root mean square value of orderq, when q≠0, with(4)Fqs=12Ns∑v=12NqFvsq1/q.Whenq=0, the calculation formula is(5)lnFq=0s=12Ns∑v=12NqlnFvs.Step 5: change the value of scales to obtain the corresponding Fqs; if the signal sequence has a power-law relationship scale-invariant signal, then(6)Fqs∼sHq,whereHq is the q-order Hurst exponent.Find the scaling range with a linear relationship betweenlog10Fqs and log10s. Find the slope of its linear regression line, denoted as Hq, by calculating the first-order polynomial fit coefficient between the scale s and the logarithmic value of the function Fqs in the scaling range. If the relationship between log10Fqs and log10s is S-shaped curve, the value of Hq cannot be estimated.Step 6: the correlation betweenHq and the multifractal mass exponent τq can be expressed as(7)τq=qHq−Df,whereDf denotes the topological dimension of the multifractal signal, and for a one-dimensional time series of signals, Df=1.Step 7: the singularity intensity functionαq and the singularity dimension fα can be obtained by the Legendre transformation:(8)αq=dτqdq,fα=infqqαq−τ9.In the process of changingq value from negative to positive, if Hq remains almost constant, the relationship between q-order mass exponent τq and q is linear; then, the signal has single fractal characteristics and the width of the singularity spectrum is small and tends to 0. On the contrary, if Hq changes to a large extent, the relationship between q-order mass exponent τq and q is not linear; then, the signal has multifractal characteristics and the width of the singularity spectrum is large. In this paper, we use MFDFA-1 method with one-direction partition process, just do the partitioning process from head to tail of the sequence yt, thus obtaining Ns equal-length subintervals.


## 3. The Multifractal Singularity Power Spectrum Algorithm (SPS)

### 3.1. SPS Algorithm

Suppose a set of signals or time series is a continuous signal.Definingxαt=xt:αxt=α as a fractal subband signal or a fractal subset of signal xt, the time-indexed set corresponding to xαt constitutes the set tα=t,αxt=α; then, xαt is denoted as(9)xαt=xt,t∈tα,0,t∈0,Ttα,where αxt is the singularity exponent of xt at the time t. According to the theory of multifractals, xt is a thick subset supported by a tight set ∪αxαt. Meanwhile, xα1t,xα2t,xα3t,… form the segments of xt, where xαit are mutually disjoint, and we obtain xαit∩xαjt=φ,i≠j.All measurable subsets constitute the set systemI=xα1t,xα2t,xα3t,…, and I is a topological measurable set. There exists a measurable function f:I⟶ℜ+, and X,I,f constitutes a measurable space.Based on the singularity decomposition, we can obtain the singular subsetxTαt of xTt under the fractal truncation signal. The fractal energy measure of xTαt can be expressed as(10)PxT=∫αPxTαdα,where PxTα is the singularity measure of each fractal subset xTαt and the singularity power spectral density of xTα. It reflects the distribution of the fractal power of xTt with the singularity exponent. The fractal power analysis of the singularity fractal subset xTαt gives the singularity power spectrum of xTt as(11)pxTα=1TxTt∫−T/2+xTαt21+tg2θtdHxTαt.The fractal power measurement of the signal sequencexTt and the singularity power spectrum distribution satisfy the following relationship:(12)pxT=∫αpxTαdα=∫α∫−T/2+T/21TxTtxTαt21+tg2θtdHxTαtdα.Therefore, the SPS distribution function of the signalxt can be defined as a limit form:(13)Pxα=limT⟶∞pxT=∫−∞+∞limT⟶∞xTαt2TxTt1+tg2θtHxTαdt,where θi is the local orientation angle of the truncated signal xTt (equivalent to the tangent direction of a conventional curve) and dHxTαt denotes the differential of the Hausdorff measure of xTt. Based on the modulus squared of the signal, the local azimuth of the fractal time element, and the Hausdorff measure, equation (15) represents the power distribution of the fractal signal concerning the singularity exponent. Thus, the fractal energy measurement and the singularity power spectrum satisfy the relationship(14)Px=∫αPxαdα.For fractal signals with a periodT or finite time support xt, t∈−T/2,T/2. Singularity energy spectral density and singularity power spectral density functions can be defined as(15)Exα=∫−T/2T/2xαt21+tg2θtαHxdt,Pxα=1Txt∫−T/2T/2xαt21+tg2θtαHxdt.It can be obtained that the singularity power spectrum reflects the fractal energy per unit singularity exponent interval, and it reflects how the fractal energy is distributed with the variation of the singularity exponent.


### 3.2. Specific Steps of the SPS Algorithm

The multifractal singularity power spectrum (SPS) algorithm is a quantitative analysis method to study the correlation of a sequence [16]. The steps of the SPS algorithm are as follows:Step 1: calculate the instantaneous singularity exponent: for a given time seriesxt, estimate the instantaneous singularity exponent αxkn of the series.Step 2: construct the singularity subset signal: according to the instantaneous singularity exponentαxkn, set the singularity variable change interval αmin,αmax, where αmin≤αxkn≤αmax, and divide the interval into Nα discrete singularity exponent sets αq,q=1,2,…,Nα+1 according to the singularity exponent interval equivalence, where α1=αminandαNα+1=αmax. Nα is the number of discretized singularity indices.Step 3: according to the interval segment of the singularity exponent after discretization, the singularity decomposition of this time series is performed; the data points whose singularity indices fall in the intervalαi,αi+1,i=1,2,…,Nα−1, are specified as subsets belonging to the same singularity exponent αi, and the data points whose singularity indices fall in the interval αi,αi+1,i=Nα, are specified as subsets belonging to the same singularity exponent αi; thus, it is possible to obtain the fractal subsets with the same singularity exponents. We denote the subset of the same singularity exponent αq by xαq=xi,αxi=αq and the total sequence by x=xi=∪q=1Nαxαq, q=1,2,…,Nα, where xαq denotes the set of all moments in the time series and xi has the instantaneous singularity exponent value αq. The number of points in the subset xαq is related to the specific calculation, and the subset corresponding to each singularity exponent αq is denoted as xαqt.Step 4: calculate the average powerPxαq of the singularity subset with the same singularity exponent α=αq: based on doing the singularity subset decomposition on xt, for the singularity subset xαt, for a given singularity exponent α=αq, calculate the power Pxαq of the subset (time series) with the same singularity exponent; for the convenience of presentation, denote this time series with its corresponding singularity subset xαq as xαqt and its discrete form as xαqk,k=1,2,…,Nαq, where Nαq is, respectively, the number of discrete point sets contained in the time series xt, as the average power Pαq=1/Nαq∑kxαqk2.Step 5: calculate the singular power spectrum: iterate over all the singularity indices to obtain the singularity power spectrumPxα.


## 3.1. SPS Algorithm

Suppose a set of signals or time series is a continuous signal.Definingxαt=xt:αxt=α as a fractal subband signal or a fractal subset of signal xt, the time-indexed set corresponding to xαt constitutes the set tα=t,αxt=α; then, xαt is denoted as(9)xαt=xt,t∈tα,0,t∈0,Ttα,where αxt is the singularity exponent of xt at the time t. According to the theory of multifractals, xt is a thick subset supported by a tight set ∪αxαt. Meanwhile, xα1t,xα2t,xα3t,… form the segments of xt, where xαit are mutually disjoint, and we obtain xαit∩xαjt=φ,i≠j.All measurable subsets constitute the set systemI=xα1t,xα2t,xα3t,…, and I is a topological measurable set. There exists a measurable function f:I⟶ℜ+, and X,I,f constitutes a measurable space.Based on the singularity decomposition, we can obtain the singular subsetxTαt of xTt under the fractal truncation signal. The fractal energy measure of xTαt can be expressed as(10)PxT=∫αPxTαdα,where PxTα is the singularity measure of each fractal subset xTαt and the singularity power spectral density of xTα. It reflects the distribution of the fractal power of xTt with the singularity exponent. The fractal power analysis of the singularity fractal subset xTαt gives the singularity power spectrum of xTt as(11)pxTα=1TxTt∫−T/2+xTαt21+tg2θtdHxTαt.The fractal power measurement of the signal sequencexTt and the singularity power spectrum distribution satisfy the following relationship:(12)pxT=∫αpxTαdα=∫α∫−T/2+T/21TxTtxTαt21+tg2θtdHxTαtdα.Therefore, the SPS distribution function of the signalxt can be defined as a limit form:(13)Pxα=limT⟶∞pxT=∫−∞+∞limT⟶∞xTαt2TxTt1+tg2θtHxTαdt,where θi is the local orientation angle of the truncated signal xTt (equivalent to the tangent direction of a conventional curve) and dHxTαt denotes the differential of the Hausdorff measure of xTt. Based on the modulus squared of the signal, the local azimuth of the fractal time element, and the Hausdorff measure, equation (15) represents the power distribution of the fractal signal concerning the singularity exponent. Thus, the fractal energy measurement and the singularity power spectrum satisfy the relationship(14)Px=∫αPxαdα.For fractal signals with a periodT or finite time support xt, t∈−T/2,T/2. Singularity energy spectral density and singularity power spectral density functions can be defined as(15)Exα=∫−T/2T/2xαt21+tg2θtαHxdt,Pxα=1Txt∫−T/2T/2xαt21+tg2θtαHxdt.It can be obtained that the singularity power spectrum reflects the fractal energy per unit singularity exponent interval, and it reflects how the fractal energy is distributed with the variation of the singularity exponent.


## 3.2. Specific Steps of the SPS Algorithm

The multifractal singularity power spectrum (SPS) algorithm is a quantitative analysis method to study the correlation of a sequence [16]. The steps of the SPS algorithm are as follows:Step 1: calculate the instantaneous singularity exponent: for a given time seriesxt, estimate the instantaneous singularity exponent αxkn of the series.Step 2: construct the singularity subset signal: according to the instantaneous singularity exponentαxkn, set the singularity variable change interval αmin,αmax, where αmin≤αxkn≤αmax, and divide the interval into Nα discrete singularity exponent sets αq,q=1,2,…,Nα+1 according to the singularity exponent interval equivalence, where α1=αminandαNα+1=αmax. Nα is the number of discretized singularity indices.Step 3: according to the interval segment of the singularity exponent after discretization, the singularity decomposition of this time series is performed; the data points whose singularity indices fall in the intervalαi,αi+1,i=1,2,…,Nα−1, are specified as subsets belonging to the same singularity exponent αi, and the data points whose singularity indices fall in the interval αi,αi+1,i=Nα, are specified as subsets belonging to the same singularity exponent αi; thus, it is possible to obtain the fractal subsets with the same singularity exponents. We denote the subset of the same singularity exponent αq by xαq=xi,αxi=αq and the total sequence by x=xi=∪q=1Nαxαq, q=1,2,…,Nα, where xαq denotes the set of all moments in the time series and xi has the instantaneous singularity exponent value αq. The number of points in the subset xαq is related to the specific calculation, and the subset corresponding to each singularity exponent αq is denoted as xαqt.Step 4: calculate the average powerPxαq of the singularity subset with the same singularity exponent α=αq: based on doing the singularity subset decomposition on xt, for the singularity subset xαt, for a given singularity exponent α=αq, calculate the power Pxαq of the subset (time series) with the same singularity exponent; for the convenience of presentation, denote this time series with its corresponding singularity subset xαq as xαqt and its discrete form as xαqk,k=1,2,…,Nαq, where Nαq is, respectively, the number of discrete point sets contained in the time series xt, as the average power Pαq=1/Nαq∑kxαqk2.Step 5: calculate the singular power spectrum: iterate over all the singularity indices to obtain the singularity power spectrumPxα.


## 4. Extreme Learning Machine

The ELM algorithm [17] has the advantages of strong generalization ability and fast learning speed, which overcomes the disadvantages of a traditional neural network learning algorithm (such as BP neural network algorithm) that requires continuous iteration, many training times, low learning efficiency, and slow convergence speed, and is a new type of single implicit layer feedforward neural network. The ELM consists of three parts: input layer, hidden layer, and output layer, and the specific structure is shown in Figure 1.Figure 1
Structure of basic ELM.ωij is the value of the connection between the input layer and the hidden layer, which is denoted as(16)ωij=ω11ω12⋯ω1nω21ω22⋯ω2n⋮⋮⋮ωl1ωl2⋯ωln.βjk is the value of the connection between the hidden layer and the output layer, which is denoted as(17)βjk=β11β12⋯β1mβ21β22⋯β2m⋮⋮⋮βl1βl2⋯βlm,where n is the number of neurons in the input layer, l is the number of neurons in the hidden layer, and m is the number of neurons in the output layer.The threshold of the implicit layer neuronb is expressed as(18)b=b1b2⋯blT.Letgx be the hidden layer neuron excitation function, and the output T can be expressed as the following equation:(19)T=t1t2⋯tmT,(20)tj=t1jt2j⋮tmj=∑i=1lβi1ωixj+bi∑i=1lβi2ωixj+bi⋮∑i=1lβimωixj+bi,(21)ωi=ωi1ωi2⋯ωin,xj=x1jx2j⋯xnjT.Equation (22) can also be expressed as Hβ=T′, where T′ is the transpose matrix of T and H is the hidden output array, and the expression is(22)H=gω1x1+b1gω2x1+b2⋯gωlx1+blgω1x2+b1gω2x2+b2⋯gωlx2+bl⋮⋮⋮gω1xQ+b1gω2xQ+b2⋯gωlxQ+bl.The system of equations’ least squaresβ is(23)minβHβ−T′.The solution isβ′=H†T′, where H† is the Moore–Penrose generalized inverse of the output matrix of the hidden layer.


### 4.1. Extraction of Fractal Characteristic Features in Bearing Signals

Statistical physics can be used to obtain a partition function for the set of probabilities of regular and irregular fractals, which uses theq−th power of the probability as a weight to distinguish subsets with different sizes of probabilities. In some cases, the collocation function has scale invariance of the form ετq over a certain range (and over an infinite range in the case of regular fractals). Figure 2 shows the sketch of the multifractal spectrum used to extract features, Figure 3 shows the sketch of the Hurst exponents, and Figure 4 shows the sketch of the singularity power spectrum.Figure 2
Multifractal spectrum.Figure 3
Hurst exponents.Figure 4
Singularity power spectrum.The generalized fractal dimensionDf can be obtained from τq, and the multifractal spectrum fα is obtained by the Legendre transform. αmax and fαmax of the multifractal spectrum reflect the properties of the smallest subset of probabilities. αmin and fαmin reflect the properties of the largest subset of probabilities. fα0 reflects the properties of the most contingent subset. The width Δα=αmax−αmin of the multifractal spectrum reflects the size of the probability distribution range.The more inhomogeneous the probability distribution is, the wider the correspondingfα curve is and therefore the larger the Δα value becomes. In the practical problem of bearing fault signals such as having irregular distribution in space, the quantitative characterization of the inhomogeneous distribution of physical quantities can be achieved with the specific parameters of the multiple fractal spectrum, so it is suitable for the feature extraction of different kinds of bearing faults and different degrees of damage to achieve the diagnosis and identification of faults.The value of the weighting factorq is also worth discussing when calculating the multifractal spectrum of the measured signal. When calculating the multifractal spectrum according to the multifractal detrended fluctuation analysis, different values of q are used to divide the signal into regions with different levels for the study. Theoretically, the larger the range of q is −∞<q<∞, the more complete the representation of the multifractal spectrum, but in practice, as q increases, the computational effort increases exponentially, and the increase of q to a certain level will inevitably cause computer overflow errors. However, if the range of q is too small, the change of fα is still large for every increasement of Δq=1. The obtained is only a part of the multifractal spectrum, which cannot fully reflect the probability distribution of the object of study. In general, when q increases to a certain degree, regardless of the degree of a signal probability distribution, αq and fα basically do not change with the increase of q, and the range of q can be cut off. In the calculation of the stochastic fractal, a criterion for determining qmax can be chosen based on this property, and a coefficient η for the variation of f with the value of q can be defined as follows:(24)η=fq−fq−1fq−fmax.Whenq is infinite, fq=fq−1, η=0. In the actual calculation, a value slightly larger than 0 is chosen as a criterion, for example, η<0.1%.In the fractal analysis of the actual faulty bearing signal, the data are divided into segments of equal length, and each segment contains 2048 points. The efficiency can be ensured whenqmax=5, and without losing the multiple fractal spectrum features, the spectral parameters are characterized comprehensively.Mandelbrot in his book “Fractals: Form, Chance, and Dimension” [15] gives a portrayal of the Hurst phenomenon, that is, river flow in which a larger-than-average flood may follow a larger flood and then, suddenly, a below-average flood, followed by more below-average floods; the whole process like there is a cycle, but not for the periodicity, and the standard statistical analysis also reveals significant correlations between observations. To separate the fractal time series from other time series, the Hurst exponent was therefore constructed to estimate the fractal dimension of the time series and the degree of series persistence or inverse persistence. When calculating the Hurst indices for measured bearing failure signals, the Hurst indices for different failures and signals with different fault levels also reflect different distributions, and the Hurst indices corresponding to different q values are taken for feature extraction when the value of qmax is given. The features taken in this paper are included in Table 1.Table 1
The six characteristic parameters extracted by MFDFA and SPS.
MFDFAα0The abscissa of the multifractal spectrum verticesΔαThe width of the multifractal spectrumHq=−5The value of the hurst exponent whenq=−5Hq=0The value of the hurst exponent whenq=0Hq=5The value of the hurst exponent whenq=5SPS∑PαSummation of the singularity power spectrum


## 4.1. Extraction of Fractal Characteristic Features in Bearing Signals

Statistical physics can be used to obtain a partition function for the set of probabilities of regular and irregular fractals, which uses theq−th power of the probability as a weight to distinguish subsets with different sizes of probabilities. In some cases, the collocation function has scale invariance of the form ετq over a certain range (and over an infinite range in the case of regular fractals). Figure 2 shows the sketch of the multifractal spectrum used to extract features, Figure 3 shows the sketch of the Hurst exponents, and Figure 4 shows the sketch of the singularity power spectrum.Figure 2
Multifractal spectrum.Figure 3
Hurst exponents.Figure 4
Singularity power spectrum.The generalized fractal dimensionDf can be obtained from τq, and the multifractal spectrum fα is obtained by the Legendre transform. αmax and fαmax of the multifractal spectrum reflect the properties of the smallest subset of probabilities. αmin and fαmin reflect the properties of the largest subset of probabilities. fα0 reflects the properties of the most contingent subset. The width Δα=αmax−αmin of the multifractal spectrum reflects the size of the probability distribution range.The more inhomogeneous the probability distribution is, the wider the correspondingfα curve is and therefore the larger the Δα value becomes. In the practical problem of bearing fault signals such as having irregular distribution in space, the quantitative characterization of the inhomogeneous distribution of physical quantities can be achieved with the specific parameters of the multiple fractal spectrum, so it is suitable for the feature extraction of different kinds of bearing faults and different degrees of damage to achieve the diagnosis and identification of faults.The value of the weighting factorq is also worth discussing when calculating the multifractal spectrum of the measured signal. When calculating the multifractal spectrum according to the multifractal detrended fluctuation analysis, different values of q are used to divide the signal into regions with different levels for the study. Theoretically, the larger the range of q is −∞<q<∞, the more complete the representation of the multifractal spectrum, but in practice, as q increases, the computational effort increases exponentially, and the increase of q to a certain level will inevitably cause computer overflow errors. However, if the range of q is too small, the change of fα is still large for every increasement of Δq=1. The obtained is only a part of the multifractal spectrum, which cannot fully reflect the probability distribution of the object of study. In general, when q increases to a certain degree, regardless of the degree of a signal probability distribution, αq and fα basically do not change with the increase of q, and the range of q can be cut off. In the calculation of the stochastic fractal, a criterion for determining qmax can be chosen based on this property, and a coefficient η for the variation of f with the value of q can be defined as follows:(24)η=fq−fq−1fq−fmax.Whenq is infinite, fq=fq−1, η=0. In the actual calculation, a value slightly larger than 0 is chosen as a criterion, for example, η<0.1%.In the fractal analysis of the actual faulty bearing signal, the data are divided into segments of equal length, and each segment contains 2048 points. The efficiency can be ensured whenqmax=5, and without losing the multiple fractal spectrum features, the spectral parameters are characterized comprehensively.Mandelbrot in his book “Fractals: Form, Chance, and Dimension” [15] gives a portrayal of the Hurst phenomenon, that is, river flow in which a larger-than-average flood may follow a larger flood and then, suddenly, a below-average flood, followed by more below-average floods; the whole process like there is a cycle, but not for the periodicity, and the standard statistical analysis also reveals significant correlations between observations. To separate the fractal time series from other time series, the Hurst exponent was therefore constructed to estimate the fractal dimension of the time series and the degree of series persistence or inverse persistence. When calculating the Hurst indices for measured bearing failure signals, the Hurst indices for different failures and signals with different fault levels also reflect different distributions, and the Hurst indices corresponding to different q values are taken for feature extraction when the value of qmax is given. The features taken in this paper are included in Table 1.Table 1
The six characteristic parameters extracted by MFDFA and SPS.
MFDFAα0The abscissa of the multifractal spectrum verticesΔαThe width of the multifractal spectrumHq=−5The value of the hurst exponent whenq=−5Hq=0The value of the hurst exponent whenq=0Hq=5The value of the hurst exponent whenq=5SPS∑PαSummation of the singularity power spectrum


## 5. Applications

### 5.1. Dataset Introduction

The bearing signals used in the experiments were obtained from the Bearing Data Center at the Electrical Laboratory of Case Western Reserve University, USA [18, 19]. The experimental platform consists of a 1.5 kW (2 hp) electric motor (electric motor in the left side of Figure 5), a torque sensor/translator (torque transducer and encoder in the middle of Figure 5), and a power test meter (the dynamometer in Figure 5). The rotational speed is 1772 r/min, the sampling frequency is 12 kHz, and the sampling time is 10 s. Single point faults were introduced to the test bearings using electrodischarge machining with fault diameters of 7 mils, 14 mils, 21 mils, 28 mils, and 40 mils (1 mil = 0.001 inches). SKF bearings were used for the 7, 14, and 21 mils diameter faults, and NTN equivalent bearings were used for the 28 mils and 40 mils faults. All experimental data use drive-end acceleration data, which are called vibration acceleration signals. The validity of the proposed method in this paper is verified by selecting the data of the normal state shown in Figure 6(a) and 6 fault states. The fault state data are inner-race 7 mils (slight inner-race fault) shown in Figure 6(b), inner-race 28 mils (serious inner-race fault) shown in Figure 6(c), outer-race 7 mils (slight outer-race fault) shown in Figure 6(d), outer-race 28 mils (serious outer-race fault) shown in Figure 6(e), rolling elements 7 mils (slight ball fault) shown in Figure 6(f), and rolling elements 28 mils fault (serious ball fault) shown in Figure 6(g). Because of the inconsistent length of the bearing signals, each of the data shown in Figure 6 is only for a sampling time of 5 seconds. Each type of data is divided into 50 equal segments with 2048 points per segment. We have 7 kinds of signals in different states. The validation acceleration data of each state are divided into 50 sets, and we can get 350 sets. During the analysis of the MFDFA-1 algorithm for all the selected bearing signals, the subinterval segments are taken as 12 values equally spaced and rounded between 10 and 128, and the values of q are taken as 31 values equally spaced between −5,5 according to the rule of equation (24). The data storage format, the feature extraction algorithms utilized, and the classifier application process all run on the MATLAB2020a platform. In the singularity power spectrum algorithm analysis process, the selection of subintervals is to distinguish the singularity exponent α into 50 equally spaced segments. The time-domain waveform of the vibration signal of each state bearing is shown in Figure 6.Figure 5
Experimental setup.Figure 6
The waveforms of seven vibration acceleration signals in each state. (a) Normal. (b) 7 mils, slight inner-race fault. (c) 28 mils, serious inner-race fault. (d) 7 mils, slight outer-race fault. (e) 28 mils, serious outer-race fault. (f) 7 mils, slight ball fault. (g) 28 mils fault, serious ball fault.
(a)(b)(c)(d)(e)(f)(g)


### 5.2. Steps of Bearing Fault Diagnosis

To ensure the reliability of the experimental results, the bearing signals in the database are tested with the same type of bearings and the same operating conditions except for the degree of damage. The flowchart of the proposed algorithm is shown in Figure7.Step 1: the MFDFA and SPS algorithms are used to process the normal bearing signals and the bearing fault signals with different damage faults and different damage degrees to obtain the corresponding multiple fractal spectrum and SPS spectrumStep 2: feature extraction using the six specific features mentioned in the above sectionStep 3: the ELM classifier is used to classify the bearings in seven states, and multiple classification tools are taken to compare the classification effectsFigure 7
Flowchart of algorithm procedures.


### 5.3. Experimental Results

First, the MFDFA and SPS algorithms were used to analyze the bearing signals in seven states, with ten samples for each state, and the corresponding multifractal spectrum and SPS spectrum were obtained, as shown in Figures8–10.Figure 8
The multifractal spectra of seven signals based on MFDFA-1.Figure 9
The Hurst exponents of seven signals based on MFDFA-1.Figure 10
The singularity power spectra of seven signals based on the SPS method.Figure8 shows multifractal spectra of seven signals based on MFDFA-1, where MFDFA-1 means the order of the polynomial of MFDFA is 1. It can be seen that the multifractal spectral lines in different states are hump-shaped and the spectrum widths and the location of the spectrum vertices of the corresponding multifractal spectra are different from each other. The location of the top of the multifractal spectrum of the normal bearing signal is in the rightmost part of Figure 8, and the length of the spectrum width is the shortest. The locations of the top of the spectra of the rest of faulty bearing signals are shifted to the left side of the spectrum of the normal signal. The spectrum widths of the rest of the faulty bearing signals are wider than the one of the normal signals, which indicates that the internal probability density of the faulty bearing signal has changed to different degrees.The multifractal spectrum width of the normal bearing vibration acceleration signal is narrower than that of the faulty bearing signal, which indicates that the distribution of the probability of the amplitude of the normal operation signal is more even. An increase in the multifractal spectrum width usually implies an increase in the irregular properties within the signal. The purple solid line corresponds to a larger value offαmax for slight outer-race fault signal, indicating that the signal has a larger number of subsets with small probability, and the αmin value is the smallest, which corresponds to the subsets with maximum probability. Compared with the serious outer-race fault (olive green dotted line), the two multifractal spectra are in the shape of right hook, which indicates that the signal of the outer-race fault signal is predominated by the subset with small probability. The α0 positions of the two signals are different, indicating that the most probable subsets are different. The spectrum widths of the two signals are similar to each other, which means that the probability density distributions of the corresponding signals are similar to each other, and the difference lies in the range of probability distributions of the two kinds of signals.The values offαmin corresponding to slight ball fault (blue-green dotted line) and serious ball fault (red solid line) are at a large level. The spectrum widths Δα of the spectra of the two kinds of signals are smaller than those of the other faults, indicating that the probability distribution of the ball fault signal is more uniform than the probability distributions of the other faults. The two spectrum lines are in the shape of left hook, which means that the large probability subsets dominate the signal.The spectrum lines of slight inner-race fault (red solid line) and serious inner-race fault (blue solid line) are both left hook-shaped, and theα0 positions of two curves are also close to each other. The values of fαmax and fαmin of serious inner-race fault signal are smaller than those of slight inner-race fault signal. This means that the number of both the small and large probability subsets of the serious inner-race fault signal is smaller. The value of Δα for serious inner-race fault signal is larger than that of slight inner-race fault signal, so the probability density distribution of serious inner-race fault signal is more uneven than that of slight inner-race fault signal.By comparing the multifractal spectra of signals with the same fault and different degrees of damage, it can be found that the multifractal spectra of signals with different degrees of damage for the same fault have similar shapes (left hook or right hook), but the difference lies in the difference in the uniformity of the probability density distribution (the size ofΔα value).Figure9 shows the q-order Hurst exponent Hq calculated by MFDFA-1. It can be seen in Figure 9 that when q is changing from −5 to 5 by increasement of Δq=1/3, the curve of Hq for the bearing vibration acceleration signal in the normal state is different from the line of Hq for the overall faulty bearing vibration acceleration signal distribution, indicating that once the bearing fails, the internal tortuous state of the signal will change. The q-order Hurst exponent Hq of serious bearing fault signal is generally larger than that of slight bearing faults. Hq for the normal state is calculated in the range of q∈−5,5; it is uniformly distributed. The distinction is evident when comparing with the q-order Hurst exponent Hq values of faulty signals. When q<0, the Hq values for the signals of serious inner-race fault, serious outer-race fault, and serious rolling element fault are greater than the Hq values of the corresponding slight fault signals. When q>0, as q increases, there appear noticeable differences of the Hq values for signals of slight inner-race fault, slight outer-race fault, and slight ball fault. Therefore, the two points at the endpoints Hq=−5andHq=5 with the points in the middle Hq=0 can be chosen as the feature values to distinguish the degree of fault for a given value of q. In the range q∈−5,−4, the Hq values for the signals of the serious outer-race fault and the serious inner-race fault are more distinguishable from the rest of the fault signals, and the curves of the Hq values for the signals of the two ball fault signals are difficult to be distinguished from the curves of Hq values for the slight outer-race fault and the slight inner-race fault in this interval. However, with the increase of q value, the Hq values for the signal of slight ball fault do not drop as rapidly as the Hq values for the rest of the fault signals. The Hq curve of slight ball fault signal is distinguishable from that of the serious ball fault in the interval q∈3,5. Similarly, as the value of q increases, the difference between Hq curves of the slight outer-race fault signal and serious outer-race fault signal is almost the same except in the range q∈−5,−2. Conversely, the Hq curves of signals of the two inner-race faults and the serious outer-race faults in the range q∈2,5 are difficult to be distinguished from each other. The difference between the Hq curves of the signals of the two kinds of inner-race faults is greater in the range q∈−5,0 than that in the interval q∈1,5. When q>0, Hq reflects the degree of fluctuation of the major trend within the signal, and when q<0, Hq reflects the degree of fluctuation of the minor trend within the signal. As a result, it can be obtained that the degree of fluctuation of minor trends within the signal of serious degree of bearing fault is usually greater than that of the signal of slight degree of bearing fault when q∈−5,0.Figure10 shows the singularity power spectra of seven signals based on the SPS method. In Figure 10, the characteristics of the SPS spectra of the faulty bearing signals in different states are not the same and the power spectra at specific different singularity values changed to a larger extent due to the shift in the distribution range of the singularity values. The normal bearing signal has a large singularity exponent distribution and a uniform power distribution. The distribution range of the singularity exponents of the faulty bearing signal changes, and the power distribution becomes uneven compared to the normal state. According to the occurrence of such a phenomenon, extracting the features with a large degree of variation can effectively distinguish the bearing fault and the degree of failure. The specific distribution of the 5 feature parameters α0,Δα,Hq=−5,Hq=0,andHq=5 of the 70 samples (all training sets consisted of 10 sets of each state, so we get 70 sets for training) calculated by MFDFA-1 and ∑Pα calculated by using the SPS method are shown in Figures 10–15. The horizontal coordination denotes the number index of the sample signal. The vertical coordination denotes the feature parameter chosen for classification. The values of the six specific characteristic parameters are shown in Table 2.Figure 11
The distribution range of featureα0 calculated by MFDFA-1.Figure 12
The distribution range of featureΔα calculated by MFDFA-1.Figure 13
The distribution range of featureHq=−5 calculated by MFDFA-1.Figure 14
The distribution range of featureHq=0 calculated by MFDFA-1.Figure 15
The distribution range of featureHq=5 calculated by MFDFA-1.Table 2
The detailed descriptions of seven vibration signals.
Operational statusα0ΔαHq=−5Hq=0Hq=5∑PαNormal0.94280.19200.97091.01800.97190.1640.90150.05050.88900.90160.92200.18890.86040.08100.86340.88100.87930.20010.94440.12050.94940.97770.96010.20480.8790.24750.90311.00360.90640.1870.84120.08750.85210.86330.86460.19410.89530.03580.87280.88390.90020.21250.88450.14720.92520.97040.92960.17950.89430.23200.93181.00880.93480.21520.88010.27730.95921.04470.96590.1808Slight inner-race fault0.02120.71660.24690.47840.26803.81630.0170.76090.24100.49530.26263.46120.02820.74460.23990.49680.26123.62630.01940.76370.24460.49340.26553.91880.01990.75140.24230.49260.26353.68430.03380.79900.23490.52160.25583.93120.02320.74520.24080.50130.26293.83480.0260.83900.24420.53250.26473.84260.02060.76320.23940.49590.26063.72960.01650.69420.23820.47220.26014.4186Serious inner-race fault0.01791.04540.29570.69710.32476.55670.01381.00070.30160.64990.32959.7837−0.00380.99710.29930.65810.32958.4325−0.02531.03160.28640.65470.31649.27390.00351.05470.29260.66720.32116.150.00270.99440.30190.67170.33297.4287−0.00391.17390.30560.74620.33615.56310.01781.01910.28990.66040.31588.1085−0.01721.18450.30050.74660.33177.18770.00321.16860.29950.75430.32977.1261Slight outer-race fault−0.1130.82480.24080.48750.278015.47−0.12010.83910.25840.51130.300714.3396−0.11660.82490.23920.48470.276614.1815−0.10820.83110.23010.47990.263913.2106−0.11680.86550.25040.51290.290012.5522−0.11690.85590.24260.50090.279715.8954−0.11790.94780.24210.52990.278814.9106−0.10990.84900.25010.51780.289517.2239−0.11110.81750.24460.49510.283421.1568−0.10811.05470.23370.58150.269918.6271Serious outer-race fault−0.00350.91190.39160.66940.432115.9467−0.01351.18380.39720.82840.444213.4769−0.04681.01790.38760.69460.430419.988−0.02311.00090.38870.71310.434217.8346−0.031.06570.39650.75330.443317.2280.00410.99670.39210.75760.438711.4255−0.01690.99720.37270.70980.414417.3455−0.03111.09280.38040.76040.426122.7398−0.02261.04150.38860.73330.434223.2125−0.01110.99410.38750.73210.434022.9646Slight ball fault0.21360.52860.30910.49420.31750.93910.20060.42390.29350.42540.30080.95050.22550.39320.30930.43950.31660.84770.23180.28920.29610.38940.30160.98750.23530.43360.30750.45800.31440.96260.2450.28820.31010.40530.31580.93580.22340.41320.31500.45920.32370.84140.23080.40950.29420.44480.30090.9340.24470.43840.31590.46890.32270.89980.23110.53420.31530.51910.32430.903Serious ball fault0.30110.55760.41280.60610.42290.67630.29250.46270.38970.53520.39750.39660.25330.64870.39980.59210.41050.59590.27510.49470.38350.54750.39300.54890.29310.41890.39270.53010.40130.61710.26970.40430.36230.48610.36980.50540.28230.45700.38610.53810.39520.51890.28430.41490.37820.50470.38560.63140.28210.39310.38730.49160.39440.82070.28310.45110.38730.52280.39530.657In Figure11, the normal state signal has the largest value of α0 and the α0 value of the signal for the occurrence of faults is smaller than the α0 value of the normal state. Similarly, the α0 values of the same faults at different degrees present a larger α0 value for the signal with serious degree faults than for the signal with slight degree faults. This indicates that the two parameters Δα and α0 can be used as characteristic parameters to distinguish the normal state from fault states and also provides some possibilities to distinguish the degree of the fault as well.In Figure12, the distribution of Δα values for normal signals is at a low level, accompanied by an increase in Δα values in the case of faults, while the Δα values for faults of the same type but of different degrees are larger for serious faults than for slight faults.In Figure13, when q=−5, the values of the q-order Hurst exponent Hq=−5 for different degrees of the same fault start to be distinguishable, especially for serious and slight faults in the inner-race as well as for serious and slight faults in the outer-race. Therefore, the q-order Hurst exponent value at q=−5 is sensitive to the degree of inner- and outer-race faults and is suitable for extraction as a feature parameter.In Figure14, when q=0, Hq=0 corresponding to serious rolling element faults can be distinguished as a whole from those corresponding to slight rolling element faults. It also shows good differentiation between serious outer-race faults and slight outer-race faults.In Figure15, when q=5, Hq=5 for the signal of rolling element fault is relatively larger than those for the inner-race fault and outer-race fault values, but it is less sensitive to the degree of damage of the fault.The singularity power spectra of the seven signals also exhibit distinguishability between normal and fault states. In Figure10, the distribution of singularity exponents α for the normal state lies in the rightmost of the figure. The distribution of the singularity power of the signal in the normal state is flat compared to the singularity power of the signal in the faulty state. Among them, in the distribution of α, the serious inner-race fault signal lies in the leftmost of the figure, and the serious outer-race fault signal has the largest power extrema in its distribution range of α. The values of the singularity power of slight inner-race fault signal are overall smaller than those of severe degree inner-race fault signal. However, the singularity power spectra cannot be separated from each other evidently.In Figure16, the distribution of ∑Pα for the seven states of the bearing signal is depicted. The distribution of ∑Pα for serious outer-race fault signal is the widest; comparing to slight outer-race fault signal, there are crossovers existed in some samples, but the majority of samples have greater ∑Pα values than slight degrees of outer-race faults. ∑Pα distinguishes the degree of two inner-race fault signals, and the distributions of both feature lines are distinguished evidently in Figure 16. Since the singularity power spectrum distribution of ball fault signal is stable, there is only a minor difference between the ∑Pα feature lines of serious ball fault signal and slight ball fault signal, but ∑Pα is more stable for the samples.Figure 16
The distribution range of feature∑Pα calculated by the SPS method.Some valuable patterns can be drawn from these plots; the magnitude of the spectrum vertex positionα0 usually increases as the degree of damage increases for the same fault type. The multifractal spectrum width Δα is more likely to increase for severely damaged faults than for slightly damaged faults. Except for rolling element faults, the SPS summation ∑Pα feature values for severely damaged inner and outer-race faults are usually larger than those for slightly damaged inner- and outer-race faults. The above analysis proves that the six characteristics selected in this paper can effectively quantify the types of bearing faults and the degrees of damage.According to the representation law of these six features, the six parameters are composed into a vector to represent the different bearing state characteristics.For classification of types using ELM, the number of hidden neurons is 50 and the transfer function is sigmoidal function. The training set consisting of 70 samples is fed into the ELM classifier, and the remaining samples after segmenting the data are used as the test set. The ELM classification results are shown in Figure17. Similarly, the training set is trained using SVM and LSSVM tools, respectively, and the same is verified for the test set. The classification results obtained are shown in Figures 18 and 19. In the multiclassification problem of bearing faults using SVM, Support Vector Machine for Regression (SVR) is chosen. The kernel function is chosen as radial basis function with a gamma parameter coefficient of 0.05. For classification of fault signals using LSSVM, the optimization function is chosen as simplex, and the training and test sets are classified using a cross-validation method.Figure 17
Classification results of the ELM model based on MFDFA-SPS.Figure 18
Classification results of the LSSVM model based on MFDFA-SPS.Figure 19
Classification results of the SVM model based on MFDFA-SPS.As can be seen in Figures20 and 21, we have 7 kinds of signals in different states; the validation acceleration data of each state are divided into 50 sets, so we can get 350 sets. All training sets consisted of 10 sets of each state, so we get 70 sets for training, and the number of total test sets is 350. The average runtime and correct recognition rates of five methods are shown in Table 3. The diagnostic results of the five classification tools illustrate that this feature extraction method can distinguish different kinds of bearing fault signals with different faults. Among the ELM classifier results, there are no errors in the test set except for only one data belonging to the serious rolling element fault in the test set is classified into the slight rolling element fault with a correct rate of 99.25%. In the results obtained by the MFDFA-SPS-based detection method with SVM classification, the highest probability of errors was found for ball faults, with an overall correct rate of 95%. In the LSSVM classification results, the correct rate also reached 95.71% due to the excessive number of parameters that needed to be preset and did not have optimization.Figure 20
Classification results of the ELM model based on MFDFA.Figure 21
Classification results of the ELM model based on SPS.Table 3
Average runtime and recognition accuracy of five methods.
AlgorithmAverage runtime (s)Average recognition accuracy (%)MFDFA + ELM54.4192.46SPS + ELM2.0580.83MFDFA-SPS + ELM54.5699.25MFDFA-SPS + LSSVM55.5295MFDFA-SPS + SVM54.2295.71Therefore, based on the overall classification correctness, it can be obtained that the MFDFA-SPS-based feature extraction method can be applied to the fault diagnosis of bearings, and the classification method can be further optimized to achieve better classification results.


## 5.1. Dataset Introduction

The bearing signals used in the experiments were obtained from the Bearing Data Center at the Electrical Laboratory of Case Western Reserve University, USA [18, 19]. The experimental platform consists of a 1.5 kW (2 hp) electric motor (electric motor in the left side of Figure 5), a torque sensor/translator (torque transducer and encoder in the middle of Figure 5), and a power test meter (the dynamometer in Figure 5). The rotational speed is 1772 r/min, the sampling frequency is 12 kHz, and the sampling time is 10 s. Single point faults were introduced to the test bearings using electrodischarge machining with fault diameters of 7 mils, 14 mils, 21 mils, 28 mils, and 40 mils (1 mil = 0.001 inches). SKF bearings were used for the 7, 14, and 21 mils diameter faults, and NTN equivalent bearings were used for the 28 mils and 40 mils faults. All experimental data use drive-end acceleration data, which are called vibration acceleration signals. The validity of the proposed method in this paper is verified by selecting the data of the normal state shown in Figure 6(a) and 6 fault states. The fault state data are inner-race 7 mils (slight inner-race fault) shown in Figure 6(b), inner-race 28 mils (serious inner-race fault) shown in Figure 6(c), outer-race 7 mils (slight outer-race fault) shown in Figure 6(d), outer-race 28 mils (serious outer-race fault) shown in Figure 6(e), rolling elements 7 mils (slight ball fault) shown in Figure 6(f), and rolling elements 28 mils fault (serious ball fault) shown in Figure 6(g). Because of the inconsistent length of the bearing signals, each of the data shown in Figure 6 is only for a sampling time of 5 seconds. Each type of data is divided into 50 equal segments with 2048 points per segment. We have 7 kinds of signals in different states. The validation acceleration data of each state are divided into 50 sets, and we can get 350 sets. During the analysis of the MFDFA-1 algorithm for all the selected bearing signals, the subinterval segments are taken as 12 values equally spaced and rounded between 10 and 128, and the values of q are taken as 31 values equally spaced between −5,5 according to the rule of equation (24). The data storage format, the feature extraction algorithms utilized, and the classifier application process all run on the MATLAB2020a platform. In the singularity power spectrum algorithm analysis process, the selection of subintervals is to distinguish the singularity exponent α into 50 equally spaced segments. The time-domain waveform of the vibration signal of each state bearing is shown in Figure 6.Figure 5
Experimental setup.Figure 6
The waveforms of seven vibration acceleration signals in each state. (a) Normal. (b) 7 mils, slight inner-race fault. (c) 28 mils, serious inner-race fault. (d) 7 mils, slight outer-race fault. (e) 28 mils, serious outer-race fault. (f) 7 mils, slight ball fault. (g) 28 mils fault, serious ball fault.
(a)(b)(c)(d)(e)(f)(g)


## 5.2. Steps of Bearing Fault Diagnosis

To ensure the reliability of the experimental results, the bearing signals in the database are tested with the same type of bearings and the same operating conditions except for the degree of damage. The flowchart of the proposed algorithm is shown in Figure7.Step 1: the MFDFA and SPS algorithms are used to process the normal bearing signals and the bearing fault signals with different damage faults and different damage degrees to obtain the corresponding multiple fractal spectrum and SPS spectrumStep 2: feature extraction using the six specific features mentioned in the above sectionStep 3: the ELM classifier is used to classify the bearings in seven states, and multiple classification tools are taken to compare the classification effectsFigure 7
Flowchart of algorithm procedures.


## 5.3. Experimental Results

First, the MFDFA and SPS algorithms were used to analyze the bearing signals in seven states, with ten samples for each state, and the corresponding multifractal spectrum and SPS spectrum were obtained, as shown in Figures8–10.Figure 8
The multifractal spectra of seven signals based on MFDFA-1.Figure 9
The Hurst exponents of seven signals based on MFDFA-1.Figure 10
The singularity power spectra of seven signals based on the SPS method.Figure8 shows multifractal spectra of seven signals based on MFDFA-1, where MFDFA-1 means the order of the polynomial of MFDFA is 1. It can be seen that the multifractal spectral lines in different states are hump-shaped and the spectrum widths and the location of the spectrum vertices of the corresponding multifractal spectra are different from each other. The location of the top of the multifractal spectrum of the normal bearing signal is in the rightmost part of Figure 8, and the length of the spectrum width is the shortest. The locations of the top of the spectra of the rest of faulty bearing signals are shifted to the left side of the spectrum of the normal signal. The spectrum widths of the rest of the faulty bearing signals are wider than the one of the normal signals, which indicates that the internal probability density of the faulty bearing signal has changed to different degrees.The multifractal spectrum width of the normal bearing vibration acceleration signal is narrower than that of the faulty bearing signal, which indicates that the distribution of the probability of the amplitude of the normal operation signal is more even. An increase in the multifractal spectrum width usually implies an increase in the irregular properties within the signal. The purple solid line corresponds to a larger value offαmax for slight outer-race fault signal, indicating that the signal has a larger number of subsets with small probability, and the αmin value is the smallest, which corresponds to the subsets with maximum probability. Compared with the serious outer-race fault (olive green dotted line), the two multifractal spectra are in the shape of right hook, which indicates that the signal of the outer-race fault signal is predominated by the subset with small probability. The α0 positions of the two signals are different, indicating that the most probable subsets are different. The spectrum widths of the two signals are similar to each other, which means that the probability density distributions of the corresponding signals are similar to each other, and the difference lies in the range of probability distributions of the two kinds of signals.The values offαmin corresponding to slight ball fault (blue-green dotted line) and serious ball fault (red solid line) are at a large level. The spectrum widths Δα of the spectra of the two kinds of signals are smaller than those of the other faults, indicating that the probability distribution of the ball fault signal is more uniform than the probability distributions of the other faults. The two spectrum lines are in the shape of left hook, which means that the large probability subsets dominate the signal.The spectrum lines of slight inner-race fault (red solid line) and serious inner-race fault (blue solid line) are both left hook-shaped, and theα0 positions of two curves are also close to each other. The values of fαmax and fαmin of serious inner-race fault signal are smaller than those of slight inner-race fault signal. This means that the number of both the small and large probability subsets of the serious inner-race fault signal is smaller. The value of Δα for serious inner-race fault signal is larger than that of slight inner-race fault signal, so the probability density distribution of serious inner-race fault signal is more uneven than that of slight inner-race fault signal.By comparing the multifractal spectra of signals with the same fault and different degrees of damage, it can be found that the multifractal spectra of signals with different degrees of damage for the same fault have similar shapes (left hook or right hook), but the difference lies in the difference in the uniformity of the probability density distribution (the size ofΔα value).Figure9 shows the q-order Hurst exponent Hq calculated by MFDFA-1. It can be seen in Figure 9 that when q is changing from −5 to 5 by increasement of Δq=1/3, the curve of Hq for the bearing vibration acceleration signal in the normal state is different from the line of Hq for the overall faulty bearing vibration acceleration signal distribution, indicating that once the bearing fails, the internal tortuous state of the signal will change. The q-order Hurst exponent Hq of serious bearing fault signal is generally larger than that of slight bearing faults. Hq for the normal state is calculated in the range of q∈−5,5; it is uniformly distributed. The distinction is evident when comparing with the q-order Hurst exponent Hq values of faulty signals. When q<0, the Hq values for the signals of serious inner-race fault, serious outer-race fault, and serious rolling element fault are greater than the Hq values of the corresponding slight fault signals. When q>0, as q increases, there appear noticeable differences of the Hq values for signals of slight inner-race fault, slight outer-race fault, and slight ball fault. Therefore, the two points at the endpoints Hq=−5andHq=5 with the points in the middle Hq=0 can be chosen as the feature values to distinguish the degree of fault for a given value of q. In the range q∈−5,−4, the Hq values for the signals of the serious outer-race fault and the serious inner-race fault are more distinguishable from the rest of the fault signals, and the curves of the Hq values for the signals of the two ball fault signals are difficult to be distinguished from the curves of Hq values for the slight outer-race fault and the slight inner-race fault in this interval. However, with the increase of q value, the Hq values for the signal of slight ball fault do not drop as rapidly as the Hq values for the rest of the fault signals. The Hq curve of slight ball fault signal is distinguishable from that of the serious ball fault in the interval q∈3,5. Similarly, as the value of q increases, the difference between Hq curves of the slight outer-race fault signal and serious outer-race fault signal is almost the same except in the range q∈−5,−2. Conversely, the Hq curves of signals of the two inner-race faults and the serious outer-race faults in the range q∈2,5 are difficult to be distinguished from each other. The difference between the Hq curves of the signals of the two kinds of inner-race faults is greater in the range q∈−5,0 than that in the interval q∈1,5. When q>0, Hq reflects the degree of fluctuation of the major trend within the signal, and when q<0, Hq reflects the degree of fluctuation of the minor trend within the signal. As a result, it can be obtained that the degree of fluctuation of minor trends within the signal of serious degree of bearing fault is usually greater than that of the signal of slight degree of bearing fault when q∈−5,0.Figure10 shows the singularity power spectra of seven signals based on the SPS method. In Figure 10, the characteristics of the SPS spectra of the faulty bearing signals in different states are not the same and the power spectra at specific different singularity values changed to a larger extent due to the shift in the distribution range of the singularity values. The normal bearing signal has a large singularity exponent distribution and a uniform power distribution. The distribution range of the singularity exponents of the faulty bearing signal changes, and the power distribution becomes uneven compared to the normal state. According to the occurrence of such a phenomenon, extracting the features with a large degree of variation can effectively distinguish the bearing fault and the degree of failure. The specific distribution of the 5 feature parameters α0,Δα,Hq=−5,Hq=0,andHq=5 of the 70 samples (all training sets consisted of 10 sets of each state, so we get 70 sets for training) calculated by MFDFA-1 and ∑Pα calculated by using the SPS method are shown in Figures 10–15. The horizontal coordination denotes the number index of the sample signal. The vertical coordination denotes the feature parameter chosen for classification. The values of the six specific characteristic parameters are shown in Table 2.Figure 11
The distribution range of featureα0 calculated by MFDFA-1.Figure 12
The distribution range of featureΔα calculated by MFDFA-1.Figure 13
The distribution range of featureHq=−5 calculated by MFDFA-1.Figure 14
The distribution range of featureHq=0 calculated by MFDFA-1.Figure 15
The distribution range of featureHq=5 calculated by MFDFA-1.Table 2
The detailed descriptions of seven vibration signals.
Operational statusα0ΔαHq=−5Hq=0Hq=5∑PαNormal0.94280.19200.97091.01800.97190.1640.90150.05050.88900.90160.92200.18890.86040.08100.86340.88100.87930.20010.94440.12050.94940.97770.96010.20480.8790.24750.90311.00360.90640.1870.84120.08750.85210.86330.86460.19410.89530.03580.87280.88390.90020.21250.88450.14720.92520.97040.92960.17950.89430.23200.93181.00880.93480.21520.88010.27730.95921.04470.96590.1808Slight inner-race fault0.02120.71660.24690.47840.26803.81630.0170.76090.24100.49530.26263.46120.02820.74460.23990.49680.26123.62630.01940.76370.24460.49340.26553.91880.01990.75140.24230.49260.26353.68430.03380.79900.23490.52160.25583.93120.02320.74520.24080.50130.26293.83480.0260.83900.24420.53250.26473.84260.02060.76320.23940.49590.26063.72960.01650.69420.23820.47220.26014.4186Serious inner-race fault0.01791.04540.29570.69710.32476.55670.01381.00070.30160.64990.32959.7837−0.00380.99710.29930.65810.32958.4325−0.02531.03160.28640.65470.31649.27390.00351.05470.29260.66720.32116.150.00270.99440.30190.67170.33297.4287−0.00391.17390.30560.74620.33615.56310.01781.01910.28990.66040.31588.1085−0.01721.18450.30050.74660.33177.18770.00321.16860.29950.75430.32977.1261Slight outer-race fault−0.1130.82480.24080.48750.278015.47−0.12010.83910.25840.51130.300714.3396−0.11660.82490.23920.48470.276614.1815−0.10820.83110.23010.47990.263913.2106−0.11680.86550.25040.51290.290012.5522−0.11690.85590.24260.50090.279715.8954−0.11790.94780.24210.52990.278814.9106−0.10990.84900.25010.51780.289517.2239−0.11110.81750.24460.49510.283421.1568−0.10811.05470.23370.58150.269918.6271Serious outer-race fault−0.00350.91190.39160.66940.432115.9467−0.01351.18380.39720.82840.444213.4769−0.04681.01790.38760.69460.430419.988−0.02311.00090.38870.71310.434217.8346−0.031.06570.39650.75330.443317.2280.00410.99670.39210.75760.438711.4255−0.01690.99720.37270.70980.414417.3455−0.03111.09280.38040.76040.426122.7398−0.02261.04150.38860.73330.434223.2125−0.01110.99410.38750.73210.434022.9646Slight ball fault0.21360.52860.30910.49420.31750.93910.20060.42390.29350.42540.30080.95050.22550.39320.30930.43950.31660.84770.23180.28920.29610.38940.30160.98750.23530.43360.30750.45800.31440.96260.2450.28820.31010.40530.31580.93580.22340.41320.31500.45920.32370.84140.23080.40950.29420.44480.30090.9340.24470.43840.31590.46890.32270.89980.23110.53420.31530.51910.32430.903Serious ball fault0.30110.55760.41280.60610.42290.67630.29250.46270.38970.53520.39750.39660.25330.64870.39980.59210.41050.59590.27510.49470.38350.54750.39300.54890.29310.41890.39270.53010.40130.61710.26970.40430.36230.48610.36980.50540.28230.45700.38610.53810.39520.51890.28430.41490.37820.50470.38560.63140.28210.39310.38730.49160.39440.82070.28310.45110.38730.52280.39530.657In Figure11, the normal state signal has the largest value of α0 and the α0 value of the signal for the occurrence of faults is smaller than the α0 value of the normal state. Similarly, the α0 values of the same faults at different degrees present a larger α0 value for the signal with serious degree faults than for the signal with slight degree faults. This indicates that the two parameters Δα and α0 can be used as characteristic parameters to distinguish the normal state from fault states and also provides some possibilities to distinguish the degree of the fault as well.In Figure12, the distribution of Δα values for normal signals is at a low level, accompanied by an increase in Δα values in the case of faults, while the Δα values for faults of the same type but of different degrees are larger for serious faults than for slight faults.In Figure13, when q=−5, the values of the q-order Hurst exponent Hq=−5 for different degrees of the same fault start to be distinguishable, especially for serious and slight faults in the inner-race as well as for serious and slight faults in the outer-race. Therefore, the q-order Hurst exponent value at q=−5 is sensitive to the degree of inner- and outer-race faults and is suitable for extraction as a feature parameter.In Figure14, when q=0, Hq=0 corresponding to serious rolling element faults can be distinguished as a whole from those corresponding to slight rolling element faults. It also shows good differentiation between serious outer-race faults and slight outer-race faults.In Figure15, when q=5, Hq=5 for the signal of rolling element fault is relatively larger than those for the inner-race fault and outer-race fault values, but it is less sensitive to the degree of damage of the fault.The singularity power spectra of the seven signals also exhibit distinguishability between normal and fault states. In Figure10, the distribution of singularity exponents α for the normal state lies in the rightmost of the figure. The distribution of the singularity power of the signal in the normal state is flat compared to the singularity power of the signal in the faulty state. Among them, in the distribution of α, the serious inner-race fault signal lies in the leftmost of the figure, and the serious outer-race fault signal has the largest power extrema in its distribution range of α. The values of the singularity power of slight inner-race fault signal are overall smaller than those of severe degree inner-race fault signal. However, the singularity power spectra cannot be separated from each other evidently.In Figure16, the distribution of ∑Pα for the seven states of the bearing signal is depicted. The distribution of ∑Pα for serious outer-race fault signal is the widest; comparing to slight outer-race fault signal, there are crossovers existed in some samples, but the majority of samples have greater ∑Pα values than slight degrees of outer-race faults. ∑Pα distinguishes the degree of two inner-race fault signals, and the distributions of both feature lines are distinguished evidently in Figure 16. Since the singularity power spectrum distribution of ball fault signal is stable, there is only a minor difference between the ∑Pα feature lines of serious ball fault signal and slight ball fault signal, but ∑Pα is more stable for the samples.Figure 16
The distribution range of feature∑Pα calculated by the SPS method.Some valuable patterns can be drawn from these plots; the magnitude of the spectrum vertex positionα0 usually increases as the degree of damage increases for the same fault type. The multifractal spectrum width Δα is more likely to increase for severely damaged faults than for slightly damaged faults. Except for rolling element faults, the SPS summation ∑Pα feature values for severely damaged inner and outer-race faults are usually larger than those for slightly damaged inner- and outer-race faults. The above analysis proves that the six characteristics selected in this paper can effectively quantify the types of bearing faults and the degrees of damage.According to the representation law of these six features, the six parameters are composed into a vector to represent the different bearing state characteristics.For classification of types using ELM, the number of hidden neurons is 50 and the transfer function is sigmoidal function. The training set consisting of 70 samples is fed into the ELM classifier, and the remaining samples after segmenting the data are used as the test set. The ELM classification results are shown in Figure17. Similarly, the training set is trained using SVM and LSSVM tools, respectively, and the same is verified for the test set. The classification results obtained are shown in Figures 18 and 19. In the multiclassification problem of bearing faults using SVM, Support Vector Machine for Regression (SVR) is chosen. The kernel function is chosen as radial basis function with a gamma parameter coefficient of 0.05. For classification of fault signals using LSSVM, the optimization function is chosen as simplex, and the training and test sets are classified using a cross-validation method.Figure 17
Classification results of the ELM model based on MFDFA-SPS.Figure 18
Classification results of the LSSVM model based on MFDFA-SPS.Figure 19
Classification results of the SVM model based on MFDFA-SPS.As can be seen in Figures20 and 21, we have 7 kinds of signals in different states; the validation acceleration data of each state are divided into 50 sets, so we can get 350 sets. All training sets consisted of 10 sets of each state, so we get 70 sets for training, and the number of total test sets is 350. The average runtime and correct recognition rates of five methods are shown in Table 3. The diagnostic results of the five classification tools illustrate that this feature extraction method can distinguish different kinds of bearing fault signals with different faults. Among the ELM classifier results, there are no errors in the test set except for only one data belonging to the serious rolling element fault in the test set is classified into the slight rolling element fault with a correct rate of 99.25%. In the results obtained by the MFDFA-SPS-based detection method with SVM classification, the highest probability of errors was found for ball faults, with an overall correct rate of 95%. In the LSSVM classification results, the correct rate also reached 95.71% due to the excessive number of parameters that needed to be preset and did not have optimization.Figure 20
Classification results of the ELM model based on MFDFA.Figure 21
Classification results of the ELM model based on SPS.Table 3
Average runtime and recognition accuracy of five methods.
AlgorithmAverage runtime (s)Average recognition accuracy (%)MFDFA + ELM54.4192.46SPS + ELM2.0580.83MFDFA-SPS + ELM54.5699.25MFDFA-SPS + LSSVM55.5295MFDFA-SPS + SVM54.2295.71Therefore, based on the overall classification correctness, it can be obtained that the MFDFA-SPS-based feature extraction method can be applied to the fault diagnosis of bearings, and the classification method can be further optimized to achieve better classification results.


## 6. Conclusions

Based on the above study, the following conclusions can be drawn:(1)
By applying the bearing faults from the Case Western Reserve database, the bearing fault diagnosis method combining MFDFA with SPS and ELM can be applied to diagnose different faults as well as faults with different damage levels(2)
The MFDFA-SPS fault feature extraction method based on ELM does not require preprocessing of fault data and is suitable for real-time monitoring of the operation status of bearing faultsIf only SPS is used to extract the specific feature parameters of the faulty bearing because the distribution range of the singularity values of different degrees of the faulty signal is different (except for some serious faulty signals), it will cause the stacking of different signal spectrum lines in different intervals of distribution, and the intuitive spectrum line discrimination is not strong only from the graph, so we need to combine other feature extraction methods and choose a more effective classifier for classification to achieve a better classification effect on the rolling bearing signal.The proposed method has 3 advantages. First, the proposed method has no preprocessing procedure for the sampled data, which substantially reduces the complexity of the overall diagnostic process. Second, the feature extraction time using the MFDFA-SPS algorithm is short, which improves the efficiency of online fault diagnosis. The various bearing signal data used in this paper have only 2048 points, and the need for fault diagnosis can be achieved in a short time. Finally, the effect of classification can reach more than 99%, which ensures the reliability of the diagnosis method. Therefore, the proposed method is suitable for application to practical engineering problems.With the approach proposed in the above section, we can still achieve good diagnostic results based on the Case Western Reserve bearing database using a smaller number of training set. The future direction of research is to test the sensitivity of this method to the feature extraction and fault diagnosis effect of different bearing vibration signals, as well as to realize the composite fault diagnosis of different faults superimposed.



---

*Source: 4034477-2022-07-18.xml*