LOW SPEED ROLLING BEARING DIAGNOSTICS USING ACOUSTIC EMISSION AND HIGHER ORDER STATISTICS TECHNIQUES

Diagnostics in low speed rolling element bearings is difficult. Not only are normal frequency domain diagnostics methods not appropriate for this application, but the bearing response signals are usually immersed in background noise which make it difficult to detect these faults. Higher order statistics (HOS) techniques have been available for decades but have not been widely applied to machine condition monitoring with the exceptions of skewness and kurtosis. There is however reason to believe that these HOS techniques could play an important role in acoustic emission (AE) based condition monitoring of rolling element bearings at low speeds provided appropriate care is taken. To explore this hypothesis, AE signals at low bearing rotational speeds of 70, 80, 90 and 100 rpm respectively were used in this work for the monitoring of tapered roller bearings. In addition to the well-established statistical parameters such as mean, standard deviation, skewness and kurtosis, higher moments such as hyper flatness are considered in this study. A novel diagnostic method is proposed for fault extraction based on hyperflatness, combined with Kullback-Leibler divergence, and an indicator formula derived with the use of Lempel-Ziv Complexity is given. The Kullback-Leibler divergence is used together with the skewness and hyperflatness to obtain the Kullback-Leibler information Wave (KLW)with which the analysis is performed, and better results obtained as compared to conventional frequency domain analysis.


INTRODUCTION
Condition monitoring (CM) of rolling element bearings (REB) is commonly used in determining the operational state and health of machines to detect early stages of component degradation.Based on a study, at low speeds however, little energy is generated, and conventional spectral based vibration monitoring becomes unsuitable [1][2][3].
According to research, acoustic emission (AE) is a high frequency phenomenon whereby transient elastic waves are generated by rapid (and spontaneous) release of energy from a localized source or sources within a material [4,5].AE-based condition monitoring in rolling element bearings is not new and it has been studied widely for condition monitoring at higher speeds (> 600 rpm).However not very much has been done at low speeds between 10-600 rpm.Study showed low speeds such as these are also often associated with widely varying operating conditions which makes spectral based analysis techniques less appropriate [6].Time domain parameters like crest factor, skewness and kurtosis have been used for many years in vibration monitoring.These parameters are however not sensitive enough to detect early faults, even with AE.It is therefore tempting to use higher statistical orders which would emphasize the effects of slight irregularities in the race even more and thus solve the problem of multiple events which sequentially occur at localized measurement points during testing in AE signal processing.
Statistical moments tend to describe the shape of the amplitude distribution of vibration data collected from a bearing and are sensitive to the impact impulses.From previous research on the application of statistical moments to condition monitoring in rolling element bearings, it is known that the third and fourth normalized central moments restrain the selective range of statistical parameters.
A scholar suggests that the ability of skewness and kurtosis for fault detection decreases at low speeds.Since these indicators do not perform well at low speeds, there is a need to formulate better indicators using other higher order statistical parameters to be able to propose methods other than the traditional ones that could detect faults at lower speeds [7][8][9][10].
When calculating the statistical moments of signals measured in machines, the presence of outliers which result from noise generated by moving parts within and around the machine renders the estimated moments unstable and this could become severe for higher order statistical moments like skewness, kurtosis, hyperskewness and hyperflatness.For heavy tailed distributions this implies that these estimates have high variance and are generally too unstable to capture the properties of the distribution.Higher orders expose one to the effects of outliers (measurement errors) which cause significant variance and no longer captures the correlation with damage.
Being nonlinear functions, the utilization of higher order statistical (HOS) techniques allow the analysis of the systems operating under the influence of random inputs, where the processes deviate from Gaussianity to indicate that there is a fault developing.This stems from the property of Gaussian processes to have zero higher-order spectra.HOS techniques provide high signal to noise ratio domains in which one can perform detection, signal reconstruction, if the time domain noise is spatially correlated.HOS parameters which are defined in terms of higher-order moments of the data (orders greater than 2), contain much information if one were to look beyond the power spectrum domain.
A scholar recently proposed a feature extraction method of the AE time domain waveform signal using the largest Lyapunov exponent (LLE) algorithm thereby demonstrating that the LLE feature can detect indications of failure from AE hit parameters such as the amplitude, RMS, counts and AE Burst, energy with the Average signal level (ASL) etc.However, the process is too cumbersome and unfit for online fault detection.
Based on previous research, the use of Kullback-Leibler (KL) divergence and Lempel-Ziv complexity may offer a solution to this issue; since by introducing KL divergence one reduces the sensitivity to outliers but essentially retain the advantage of higher emphasis on irregularity [11].The Lempel-Ziv complexity is an analysis tool used for non-linear dynamic systems.Studies showed it specifically measures the generational rate of new patterns along a digital sequence and is closely related to important source properties such as entropy and compression ratio [12][13][14].KL divergence is defined as the mean of the log-likelihood ratio which is the exponent in large deviation theory.It is also known as information divergence and relative entropy that measures the distance between two density distributions.
The problem associated to conditions found in many industrial applications which include draglines in the mining industry and large rolling mills in materials processing environments, is the difficulty in detecting faults due to immersed background noise and weak signal acquired from transducers, especially at low speed and under varying load and speed conditions.According to research, higher order statistical (HOS) techniques have not found wide application in this condition monitoring area, and this is due to the fact that if a process is Gaussian then HOS provide no additional information that can be obtained from the second or higher order statistics [15][16][17][18][19][20].Most of the conventional methods of analyzing faults in literatures are limited to stationary and linear situation at high speed (above 600rpm) [21].However, there is the need to believe that these HOS techniques could play an important role in condition monitoring, provided appropriate care is taken.In problems that are non-Gaussian as in the case of consideration here, HOS techniques like the 6 th order statistical moment (hyperflatness) could play an important role in its condition monitoring (CM) and that is what this work contribution is about.By applying this analytical method, processing of acoustic emission signal at low speed and varying load condition could be very useful thereby providing details about the signal which the conventional second order statistics cannot [22,23].
The main objective of this paper is proposing a new indicator called the energy coherent factor for detecting anomalies in faulty bearings through the use of Kullback-Leibler (KL) divergence and HOS elements as it utilizes the AE signal in detecting the state of health of bearings.A new indicator is proposed to this effect, which responds to the amplitude of the energy of the emission.Sections 2 describe the proposed methods used to detect the parameters and features extracted based on statistical higher moments for KL divergence formulation with Lempel-Ziv complexity, together with the proposed indicator derived.While section 3 describes the experimental investigation.The experimental results are reported in section 4 with conclusions in section 5.

Parameter extraction based on statistical higher moments for kl formulation
A bearing in a good condition has a Gaussian acceleration probability density distribution, whereas damage in a bearing results in a non-Gaussian distribution with dominant tails, because of the relative increase in the number of high level acceleration events.AE signals are considered in this work to test the computational robustness of the proposed indicator energy coherent factor and to validate it, which helps in the analysis to capture the impact within a narrow broad band spectrum being that the impact response in the damage bearing is very short, i.e.∆ = 1  .
(1) Skewness (3rd Order Statistical Moment) Skewness characterizes the degree of asymmetry of distribution around its mean and a measure of the lopsidedness of the distribution.A distribution that is skewed to the left (i.e. the tail of the distribution is longer on the left) will have a negative skewness while a distribution that is skewed to the right will have a positive skewness.Hence it is reasonable that measurements must be conducted over a sufficiently long period to encompass at least one complete cycle of the modulation in order to be certain to have measured the maximum skewness.
In order to extract the feature parameter from a faulty signal contaminated by noise and to accurately identify the fault type, a skewness wave (SKWskewness) of the faulty damage bearing is used to compute Kullback-Liebler information wave (KLW) which is presented later.
where M is the number of data points in a small region, and  and   respectively are the mean and standard deviation of the signal series   in the  ℎ small region.The skewness formula of ( 2) is now used to first generate a skewness wave to the fault-based signal from the bearing housing.Table 1 shows the statistical moment order relation giving its moment number in increasing order in relation to the central moment and standardized moment.
where   and   are, as above the mean and standard deviation of the signal series   in the  ℎ small region. is the number of data points in a small region and = 1 − ,  =   ⁄ ≤     ⁄ where   and   are the sampling and the analysis frequencies respectively.The signal data   ( = 1 − ) is divided into smaller regions.The points of the kurtosis are connected in order to derive the kurtosis wave KW(j) [16] (3) Hyperflatness (6th Order Statistical Moment) Hyperflatness is a statistical power of the sixth order.It is a powerful means of estimating fault diagnosis in bearing signal since it belongs to the kurtosis family.

Feature extraction using kl and lempel-ziv complexity
Kullback-Leibler (KL) divergence is used to compare the current estimate of the feature variable of  1 () to the reference feature variable  2 () as expressed in (5).It is defined as where  1 and  2 are two probability distributions which have probability density  1 () and  2 () respectively and ( 1  2 ) can be used to compare the two PDFs.
Kullback-Leibler divergence captures the expected log-likelihood ratio which gives a statistical interpretation of power loss, when the wrong distribution is used for one of the hypotheses.More, technically it is a vital part of probability theory with a deep connection to large deviations theory and to statistical inference to ergodic theory i.e. outliers associated to data measurement.This could be important in the analysis if the damage detection is compensated so that in the use of HOS, these deviations which could have irregular results in the analysis are well taking care of.According to a study, KL divergence has interesting statistical properties; in particular finding parameters of a statistical model by maximizing the likelihood is analogous to finding the parameters minimizing the divergence [24,25].The KL information quantity may be calculated from the expected value of the reference feature variables given by Hence to extract the Kullback feature of the fault signals, the "Kullback-Leibler Wave (KLW)" is expressed based on the KL information quantity expressed above in (6). Therefore It is on the KLW(t) signal that the spectral analysis will be performed.
Getting the envelope wave from the absolute values of the Kullback-Leibler information Wave (KLW) which is presented below as Kullback-Leibler (KL) divergence was used together with the hyperflatness wave to obtain the final reduced signal on which the Lempel-Ziv complexity was used for fault diagnosis.To be able to detect a defect (change) in the vibrating signal, there is the need to be able to interpret the contribution of each variable to the divergence, then this contribution is normalized.This contribution which is based on the idea of can be evaluated in (9).

Complexity Measurement
The Lempel-Ziv complexity is an alternative tool for signal analysis involving nonlinear dynamics and has been used for monitoring the effect of anesthesia in patients.It is used in transforming signal to be analyzed into a data sequence whose elements are given in symbols of which such transformation is often referred to as "coarse-graining" operation.
Complexity analysis is known to help in focusing on the intrinsic characteristics of the overall dynamics of a signal and neglecting details contained in lower hierarchical components especially as it relates to AE signal whose analysis cut across a narrow broad band spectrum where the impact fault occurs in the damage bearing.Thus, it specifically measures the generational rate of new patterns caused by the subsequent change in speed along a digital sequence which is closely related to such important source as entropy and compression ratio [26,27].
The Lempel-Ziv complexity reflects on the number of all different subsequences contained in the original signal data sequence.Considering a generalized normalized complexity, where A* denotes the length of all finite length sequences over the finite symbol set A, and ()denotes the length of a sequence  ∈  * with (10).
For every ∈   , the Lempel-Ziv complexity can be expressed as where   → 0 if  → ∞, and is the number of different symbols in the symbol set A.

Proposed indicator
The reason for proposing the damage indicator is based on the contribution in the formulation from the fundamental equation at (9) which is based on the formulation by Gonzalez De la Rosa et al.Here it is further expatiated upon and a complexity measurement is introduced.The indicator is computed from the sixth statistical power, the Kullback-Leibler information Wave (KLW), kurtosis wave (KWkurtosis), skewness wave (SKWskewness) and the skewness series.The motivation for proposing this indicator as a damage analyzing tool is that it tends to describe the shape of the amplitude distribution of the acoustic signal by trying to explain its sensitivity phenomenon at low speeds as indicators in the past is needed to determine the damage in a bearing (which from the literature has some drawback such as been poor indicators when dealing with very low speeds in bearings and hence need to be improved upon).The following describes the calculation procedure: (1) Compute the Lempel-Ziv complexity C(n) of the original Kullback-Leibler wave C(KLW), and the sixth statistical moment (hyperflatness C(HFIW)), (and complexityof the series skewness, i.e.C(skewness)).
(2) Compute the original information wave of the kurtosis (KWkurtosis)and of the skewness (SKWskewness).
(3) Finally, in (12) we have the energy coherent factor formulation

EXPERIMENTAL INVESTIGATION
The test rig shown in Figure 1 was used to simulate the progression of damage in a non-linear system.Data acquired from this set-up were subsequently used to extract the features required for the computation of damage indicators based on HOS.The progression of the damage was traced with an acoustic emission signal from an acoustic emission transducer.Only the radial response was monitored in this case.

Figure 1: Experimental test rig
The test rig depicted in Figure 1 features the use of two servo-hydraulic actuators to introduce constant amplitude axial and radial loads on a test bearing.The purpose of introducing the two actuators is to allow simulating a scenario where coupling between axial and radial loads could be considered.The rotating speed for the slow rotating bearing ranges from 70 to 100 rpm.The actuator forces applied is however effectively sinusoidal.Figure 2 shows the schematic diagram of the built test rig.Three test bearing scenarios were considered.Firstly, an undamaged bearing was considered.In a second bearing grounded metallic debris was mixed with grease and introduced into the bearing.In the last bearing a simulated crack was introduced.A taper roller bearing (Timken HR 30307 J) was used to be able to artificially introduce the localized-defect, since it can be dismantled from the outer raceway.Surface damage was seeded on the outer raceway of the bearing (as shown in Figure 3) with the use of a small hand drilling machine to which a small disk was mounted, which was then used to introduce a groove on the outer raceway of the taper bearing.Details of the selected bearing are reported in Table 2.A brushless AC motor (Rockwell Automation MPL-B680B), mounted on a NSK 6309 single row bearing was used to drive the system.The angular velocity of the motor was retrieved from one of the analogue outputs available on the motor drive, (Rockwell Automation Kinetix 6000 series BM-01).This system allows continuous speed variation from 0 to 3600 rpm.A Soundwel AE sensor with model number SR 150M with a frequency range of 25-530kH was used in the experiment.
The first bearing (undamaged) was loaded sinusoidally with forces of amplitude of 300N at a frequency of 2Hz on the axial load and an amplitude of 700N at a frequency of 1Hz in the radial, while bearing two which was debris induced was loaded sinusoidally with forces of amplitude of 400N at a frequency of 2Hz in the axial direction and an amplitude of 800N at a frequency of 1Hz in the radial direction.Bearing three with crack at the outer race was loaded sinusoidally with forces of amplitude of 500N on the axial at a frequency of 2Hz and 900N at a frequency of 1Hz in the radial direction.The reason for applying different loads is to simulate real life scenarios of how the application of load affects the bearing in its life cycle.The speed of the servo motor was set at 70 rpm, 80 rpm, 90 rpm and 100rpm for bearings 1, 2, and 3 respectively.The vibration signatures for the three test bearings were collected for the four speeds, using an FFT analyzer, a National Instruments data acquisition card (BNC-2110) with a shielded BNC connector block.
Figure 4 shows the flow chart for arriving at the indicator proposed for use in this paper.First the loads are applied in both the axial and radial direction and the acoustic transducer is used to obtain the signatures from which the skewness, kurtosis and hyperflatness values are obtained which are then used to find the KLW before been applied to the indicator function.Sinusoidal axial and radial load where applied to the test bearings and acoustic emission signals were collected from the bearing.These signals were grouped into smaller groups and analyzed further obtaining skewness and the kurtosis wave to which the Kullback-Leibler principle was applied to obtain Kullback-Leibler information wave (KLW).

EXPERIMENTAL RESULTS
Figure 5 shows a typical acoustic signal obtained from the bearing test rig; the length of the data taken was 60000 data samples.It cannot simply be visually observed that the amplitude of the wave formed as shown in the figure varies within some given time interval.Frequencies of the order of 100 kHz are involved with AE and the fault characteristic frequencies caused by the defective bearing and its harmonics are difficult to detect in the corresponding spectrum by conventional FFT-based envelope analysis especially at low speeds and very low speeds (< 10 rpm), as it occurs within a narrow band spectrum their harmonics as depicted in fig.6 below are also difficult to obtain especially for the outer race cracked bearings.In order to be able to cluster at the different speeds, it is necessary to obtain an absolute value of the difference in the subgroup division of the indicator equation formulated as sometimes negative values are obtained which could cause the signal at the various speed level to cross each other when plotted causing clustering/separation of the signal to be difficult and also since the complexity of the signal grouped into groups are always positive.Skewness could sometimes give a negative result hence the denominator should always be positive and smaller than the numerator so as to obtain favorable answers from the indicators.The high energy and amplitude exhibited by the ringing pulses generated as the roller passes over the groove induced on the outer ring cannot be observed on the time signal given in Figure 5, that is why the indicator is formulated to see if they can reflect this energy exhibited and this is shown in Figure 7 by the plot shown from energy coherent factor formula.By observing figure 7, we find that at a low speed of 70 rpm the energy level of the three scenarios is clearly spelt out and as these speed increases (after a longer life cycle of use) there is a reduction in the energy level which further reduces as the speed and life cycle increases which explains the description provided by literatures on how the kurtosis value relating to that of a good bearing stays below 3.These increases as fault is introduced but which reduces over time as the fault is prolonged.
The good bearing indicator line as shown in the plot above did not show much deviation from the zero line, while the debris induced, and the outer race crack bearing indicator lines show significant deflection away and from the zero line respectively.

CONCLUSION
HOS seem like an intuitive thing to do.Examples are the use of skewness and kurtosis.However, the direct outliers (which originate from nonuniform load) that make it worse for signal analysis as a peaked-ness fault in the measured data become difficult to measure.Also, the problem that outliers create in measured data being analyzed, that is common with higher order moment is a problem it has.Lower order statistical moments such as skewness and kurtosis which are commonly used could not be easily used to identify the faults at incipient stage especially as it relates to very low speeds, but in this work the condition of bearings at very low speeds were successfully classified into clusters of good bearing, debris induced, and outer race cracked.This is because HOS are much more sensitive than lower order statistical moments though it has the problem of outliers affecting the analysis.This has been taking care of here through the use of KL divergence and Lempel-Ziv complexity.HOS is effectively combined with KL divergence and Lempel-Ziv complexity to formulate an indicator that made use of kurtosis which is said to be a poor indicator of fault especially on outer race of a bearing.With HOS, information was successfully extracted which deviated from Gaussianity, making it easier to detect and quantify non-linearities in time series especially at very low speed.For this, it becomes necessary to introduce statistical techniques to "smooth" the signal, by introducing even higher orders parameters like (HFIW in energy coherent factor), and this leads to results where the energy coherent factor is less sensitive to speed but separate well on the basis of damage.For heavily tailed distributions which imply that their estimates have high variance and are generally too unstable to capture the properties of their distribution, the energy coherent factor derived in this work was introduced to help absorb the effects of useful outliers in measurements which would have caused significant variance as found with other application analysis.With the approach developed here in this work, simple representation and interpretation of the online extracted information could be made possible.For example, different colors of light emitting diode (LED) could be used to indicate the types of faults, such that a non-expert or a simple classification algorithm may interpret the result.
[22] Caesarendra, W., Tjahjowidodo, T. 2017.A review of feature extraction methods in vibration-based condition monitoring and its application for degradation trend estimation of low-speed slew bearing.MDPI Machines, 5, 21.

Figure 2 :
Figure 2: Schematic diagram of test rig setup

Figure 3 :
Figure 3: Seeded damage on outer race of a bearing.

Figure 4 :
Figure 4: A diagram for arriving at the indicator values.

Figure 5 :
Figure 5: Acoustic emission obtained from test rig at 100rpm

Figure 6 :Figure 7
Figure 6: The FFT of the signal at 70 and 100 rpmFigure7show the evolution of the descriptor established for the acoustic signal of the different defect for the indicator.With regard to the energy coherent factor, the value obtains did consider the energy and the amplitude of the signals for the indicator.The rotational speeds used for the test are 70rpm, 80rpm, 90rpm and 100rpm.In order to be able to cluster at the different speeds, it is necessary to obtain an absolute value of the difference in the subgroup division of the indicator equation formulated as sometimes negative values are obtained which could cause the signal at the various speed level to cross each other when plotted causing clustering/separation of the signal to be difficult and also since the complexity of the signal grouped into groups are always positive.Skewness could sometimes give a negative result hence the denominator should always be positive and smaller than the numerator so as to obtain favorable answers from the indicators.The high energy and amplitude exhibited by the ringing pulses generated as the roller passes over the groove induced on the outer ring cannot be observed on the time signal given in Figure5, that is why the indicator is formulated to see if they can reflect this energy exhibited and this is shown in Figure7by the plot shown from energy coherent factor formula.
bearing signal at 100rpm

Table 1 :
Statistical moments of higher order relation.Kurtosis as a fitness parameter also offers the advantage of having high values in the presence of a faulty signal while it is usually zero when only background noise is present.The kurtosis spectral wave (KWkurtosis) was generated on the fault diagnostic signal (i.e the signal generated from the faulty bearing housing).