Encoder-based weak fault detection for rotating machinery using improved Gaussian process regression

Abstract

Encoder signal analysis has proven to be a novel and cost-effective tool for the health monitoring of rotating machinery. Nevertheless, how to effectively detect the potential fault utilizing encoder information, especially at an early stage, remains a challenging issue. In light of this limitation, an improved Gaussian process regression analysis is proposed for the weak fault detection of rotating machinery via encoder signal. In this article, the Gaussian process regression model is first introduced to estimate the instantaneous angular speed and its confidence interval. Subsequently, to improve the robustness of Gaussian process regression under weak fault conditions, a spectral density complex kernel is constructed through modeling the spectral density with a mixture of Gaussians. Finally, built upon the eigenvalue decomposition, the optimal inference approach of improved Gaussian process regression is proposed. Compared with other regression methods, the major contribution is that the new method not only enhances the weak fault-related features but also sets their confidence interval adaptively. Using the proposed improved Gaussian process regression, the interference components are suppressed, while the fault-related instantaneous angular speed outliers are accurately detected. In addition, the significance of fault can be quantitatively evaluated according to the confidence level of the improved Gaussian process regression. The simulated and experimental analyses manifest that the proposed improved Gaussian process regression method can effectively identify the early weak fault. It may offer an effective tool for early fault detection of rotating machinery in industrial applications.

Keywords

Encoder signal analysis weak fault detection Gaussian process regression spectral density complex kernel

Introduction

With the rapid advancements in the modern industry and scientific technology, rotating machines have been widely used in many crucial fields such as energy, metallurgy, and national defense.¹ Nevertheless, due to the long-term running under harsh working conditions, these machines are prone to local damages.² If they are not detected in time, the rotating machinery may eventually lead to fatal breakdowns or even catastrophic accidents.³ Accordingly, considerable attention has been paid to early weak fault detection of the rotating machinery both in industry and academia to guarantee the safety, reliability, and stability of the machines.⁴

During the past decades, condition monitoring of rotating machinery with the vibration signal has achieved considerable advancement and has thus become the most widely used tool.⁵ Despite its excellent performance, the vibration signal is, however, not always effective in practical applications.⁶ For instance, the fault signature is easily attenuated due to the long transfer path between the fault and the transducer. In the health assessment of a complex mechanical device, it is not easy to extract weak fault information because of the interference from other components. Furthermore, in some special circumstances, it is even impossible to install a vibration sensor due to environment or space limitations. Therefore, there is an urgent need to discover new dynamic information to guarantee the reliability of rotating machinery.

Recently, with the rapid advent of microelectronics, more and more rotary encoders are equipped in rotating machinery.⁵ These encoders were originally used for servo control and precise positioning of electromechanical equipment. Nevertheless, research shows that the rotary encoder signals, after proper pretreatment, can be successfully used for health assessment of electromechanical equipment including robots, wind turbines, aerospace actuators, and so on.⁷ Compared with the traditional vibration signal, the encoder signal primarily reflects the torsional vibration of the monitored machine, which is more sensitive to weak faults, transient shock, and sudden change of the stiffness. However, unlike vibration signal which can be directly used for health monitoring, the encoder signal first needs to be converted into kinetic information, like instantaneous angular speed (IAS) or instantaneous angular acceleration (IAA), for further diagnosis.⁸ As investigated in Li et al.,⁹ this conversion can be achieved in different ways, such as the numerical difference method,¹⁰ the central difference method (CDM),¹¹ and the local polynomial differentiator (LPD).⁷ Of course, the IAS can also be estimated by the Hilbert transform, but its resolution is related to the sampling frequency and is more sensitive to signal noise.^12,13 Despite the encouraging results obtained, their applications are greatly limited. This is not only because there are many tuning parameters to be selected manually by the users, but also that the converted IAS signal is easily corrupted by different sources of noise such that the fault transients cannot be easily identified.

In view of the above limitations of current methods, the development of new data-driven methods using artificial intelligence or statistical learning is the emerging and promising trend of encoder signal analysis (ESA). For example, the work of X Xu et al.¹⁴ proposed a new approach called singular spectrum analysis (SSA), which detects weak position fluctuations from encoder signals to examine the operating condition of the machine. Jing et al.¹⁵ proposed a new multi-information fusion method building on deep convolutional neural networks (DCNN) for fault detection. Although these methods have achieved desirable estimation results for the IAS, they still lack an adaptive detection mechanism for IAS outliers. Among the existed models, GPR has aroused widespread attention as a valid statistical method for data-driven modeling.¹⁶ Compared with other regression models,^17,18 the GPR is more robust to other interference and is able to set a suitable confidence interval (CI) for outliers detection. For this reason, the GPR model may provide an effective way for fault diagnosis of rotating machinery via ESA.

In theoretical investigations, the GPR integrates Bayesian nonparametric statistics,¹⁹ infinite neural networks,²⁰ kernel methods in machine learning, and spatial statistics.²¹ Research on GPR was triggered by Neal,²² who demonstrated that Bayesian neural networks became Gaussian process (GP) as the number of hidden units approached infinity, and conjectured simple inference ways in this case. These simple inference techniques became the cornerstone of subsequent GP models for machine learning.²³ These models directly assume a priority over functions, rather than parameters. By further assuming homoscedastic Gaussian noise, one can analytically conclude a posterior distribution over these functions when the training dataset are given.²⁴

Despite the GPR producing extremely satisfactory results for miscellaneous modeling tasks, there usually exist a series of noticeable problems remain to be solved. First, the choice of covariance kernel that controls the modeling power affects the performance of GPR,²¹ whereas the GPR with popular kernels, such as squared exponential kernel (SK), matérn kernel (MK), periodic kernel (PK), and rational quadratic kernel (RK), was hard to accurately estimate IAS and to set suitable CI for weak fault detection. Moreover, the GPR is typically applied to high-dimensional, nonlinear, and small-sample tasks.²⁵ As the encoder is generally used in complex electromechanical systems, a huge amount of monitoring data are collected from them every hour. Hence, it is still an open question to seek more efficient methods of covariance inversion and to select training data under massive datasets conditions in the GPR model.²⁶ Finally, poor data quality and unknown statistical properties of datasets bring serious challenges to the GPR model.²⁷ Therefore, improving the anti-interference performance and robustness of the GPR model is critical in industrial applications, especially for the health monitoring of rotating machinery.

In view of all this, we present an improved Gaussian process regression (IGPR), for weak fault detection of rotating machinery using encoder signal. In this article, a Gaussian process regression (GPR) model with the popular kernel is first introduced for IAS estimation and setting its outliers detection method called CI. With this index, the fault-related outliers are adaptively identified. Successively, a novel covariance kernel that models the spectral density with a mixture of Gaussians is proposed to enhance the weak fault feature in a data-driven manner. Finally, the expressive interpretable fast algorithm is presented to improve the generalization ability of the GPR in the big data context.

The rest of this article is organized as follows. In section “Overview of GPR model,” the advantages and limitations of the GPR model are reviewed. In section “IGPR,” aiming at those limitations, an IGPR analysis is proposed and elaborated for the weak fault detection with encoder signal. We analyze the proposed method and evaluate its results by both simulated signals and experimental data in section “Simulation analysis” and “Experimental validation,” respectively. Finally, some conclusions are presented in section “Conclusion.”

Overview of GPR model

The GPR belongs to the nonparametric statistical learning model, which can explain the probability of output, and applies Bayesian techniques for calculations.²⁸ Different from conventional machine learning models which are primarily used to model the distribution,²⁹ the GPR can not only model the statistical distribution of a process but also quantify the uncertainty of estimation and thus providing an adaptive way for ESA of rotating machinery.

The GPR has a strict statistical theoretical basis with good adaptability for dealing with complex problems such as small samples, nonlinearities, and high dimensionality. In a limited set of $D = {(x_{i}, y_{i})}_{i = 1}^{n}$ , where $x_{i} \in R^{p}$ is the observation point and $y_{i} \in R$ is the observation value, we can define a GP over function $f (x)$ . All statistical feature of the GP consists of mean function $E [f (x)] = μ (x)$ and covariance kernel function $cov (f (x), f (x')) = k (x, x')$ . Moreover, to simplify the model, the data are generally preprocessed so that its mean function is equal to 0

f (x) ~ GP (0, k (x, x'))

(1)

According to the definition of GP,³⁰ one can know that any collection of function values $f (X) = {[f (x_{1}), f (x_{2}), . . ., f (x_{n})]}^{T}$ has a joint Gaussian distribution

f (X) ~ N (0, K (X, X))

(2)

where $K$ is the $n \times n$ covariance matrix whose $(i, j) - th$ element is $K_{i j} = k (x_{i}, x_{j})$ .

Considering the noise contamination in the observation value y, a general regression model can be established by

y = f (x) + ξ

(3)

where $ξ$ indicates a Gaussian noise with mean 0 and variance $σ_{n}^{2}$ .

Note that the noise $ξ$ and function values $f (x_{i})$ are independent of each other. If we substitute a series of observation points $X = {[x_{1}, x_{2}, . . ., x_{n}]}^{T}$ into the regression model, we can get the prior distribution of observation values as $y = {[y_{1}, y_{2}, . . ., y_{n}]}^{T}$ .

y ~ N (0, K (X, X) + σ_{n}^{2} I)

(4)

To obtain the smoothed and interpolated target values $g_{*} = {[g_{* 1}, . . ., g_{* m}]}^{T}$ at observation points $X_{*} = {[x_{n + 1}, . . ., x_{n + m}]}^{T}$ , the joint prior distribution of observation values $y$ and target values $g_{*}$ are presented as below

[\begin{matrix} y \\ g_{*} \end{matrix}] ~ N (0, [\begin{matrix} K (X, X) + σ_{n}^{2} I & K {(X_{*}, X)}^{T} \\ K (X_{*}, X) & K (X_{*}, X_{*}) \end{matrix}])

(5)

The posterior distribution of target values $g_{*}$ introduced above can be represented as the following format

g_{*} | X, y, X_{*} ~ N (μ, Σ)

(6)

μ = K (X, X_{*})^{T} {(K (X, X) + σ_{n}^{2} I)}^{- 1} y

(7)

Σ = K (X_{*}, X_{*}) - K (X, X_{*})^{T} {(K (X, X) + σ_{n}^{2} I)}^{- 1} K (X, X_{*})

(8)

With all the equations introduced above, the systematic framework of the GPR model is as presented in Figure 1.

In the GPR model, the covariance kernel plays a crucial role, which actually measures the closeness between data points.³¹ Similar to the activation function it can impact the performance of the neural network.²⁰ The squared exponential kernel (SK) is a widely useful covariance kernel in GPR,³² which is defined as

k_{SK} (x, x') = s_{f}^{2} \exp (- \frac{{‖ x - x' ‖}^{2}}{2 l^{2}})

(9)

where $s_{f}$ indicates the variance of observation values and the parameter l is defined as the kernel-scale.

Clearly, if the hyperparameters, $θ = {l^{2}, s_{f}^{2}, σ_{n}^{2}}$ , are not chosen properly, the model may perform poorly and cannot represent the real distribution of the data. Therefore, the hyperparameters included in the GPR should be estimated using the training datasets. Although the Markov chain Monte Carlo (MCMC) method can carry out on GPR without the need for hyperparameters estimation,³³ it is not widely accepted due to its high computational cost. And the general method is minimizing the negative log marginal likelihood (NLML) as given below²⁶

\hat{θ} = \underset{θ}{\arg \min} L (θ) = \underset{θ}{\arg \min} (- \log p (y | X, θ))

(10)

It is noteworthy that the NLML pleasingly separates into automatically calibrated model fit and complexity penalty terms;³⁴ the equation introduced above can be represented in the following format

\begin{matrix} L (θ) = - \log p (y | X, θ) \\ = \underset{model fix}{\frac{1}{2} y^{T} Σ_{θ}^{- 1} y} + \underset{complexity penalty}{\frac{1}{2} \log det Σ_{θ}} + \frac{n}{2} \log (2 π) \end{matrix}

(11)

where $Σ_{θ} = K_{θ} + σ_{n}^{2} I$ is the covariance kernel matrix with unknown hyperparameters.

For solving the extrema of equation (10), the iterative method is derived via taking the partial derivative with respect to the hyperparameters, equating it to 0

\frac{\partial}{\partial θ_{i}} L (θ) = \frac{1}{2} tr (\sum_{θ}^{- 1} \frac{\partial Σ_{θ}}{\partial θ_{i}}) - \frac{1}{2} y^{τ} \sum_{θ}^{- 1} \frac{\partial Σ_{θ}}{\partial θ_{i}} \sum_{θ}^{- 1} y = 0

(12)

where tr() represents the trace of a matrix.

In light of the above discussions, the hyperparameters estimation can be modeled as an optimization problem. First, the hyperparameters are initialized by random values from the prior distribution $p_{0} (θ)$ . Subsequently, built upon the Bayesian framework, the optimal value is searched by minimizing the NLML. Finally, the iteration finishes when the stop condition is met. In addition, a condition was set to abort the iteration if the value of $L (θ)$ was found to be diverging. After obtaining the optimal hyperparameters, the GPR model can be constructed by means of equation (6). In sum, the flowchart of the hyperparameters estimation can be presented as follows (Figure 2).

While promising, the standard GPR model is task-specific and requires sophisticated Bayesian inference, which is more demanding than neural network models or deep learning models.³⁵ Moreover, due to the fact that weak fault can result in low kinetic energy, it is easily corrupted by strong background noise. Consequently, it is not an easy task to extract the fault-related feature from the raw signal by means of the standard GPR model. To overcome the above defects, an IGPR model is proposed for denoising and weak fault feature enhancement, which could improve the IAS estimation accuracy in a data-driven manner, but the model inference remains simple and easy to analyze.

IGPR

Spectral density complex kernel

The covariance kernel measures the similarity between data points and determines how the associated random functions change with different data structures. Any stable covariance kernel function with translation-invariant can be transformed into integral form by Bochner’s theorem³⁶

k (τ) = \int_{R^{P}} e^{2 π i s^{T} τ} ψ (d s)

(13)

where $τ = x - x'$ is the stabilization factor and $ψ$ represents the positive finite measure.

If $ψ$ has a density S, then k and S are Fourier duals, where S is the spectral density of the covariance kernel

S (s) = \int_{R^{P}} k (τ) e^{- 2 π i s^{T} τ} d τ

(14)

The above equation indicates that the properties of a stationary kernel k are entirely determined by its spectral density S. Therefore, substituting the $k_{SK}$ of equation (9) into (14), we can find that the spectral density of $k_{SK}$ is a Gaussian distribution centered at the origin

S_{S K} (s) = (2 π l^{2})^{P / 2} \exp (- 2 π^{2} l^{2} s^{2})

(15)

This phenomenon indicates that it is quite limited in the accurate approximation of distribution even with arbitrarily mixtures of SK, as they correspond only to Gaussian spectral densities centered on the origin. In view of all this, a natural approach is therefore to use a mixture of Gaussians that have nonzero means. This approach models the spectral density S as a Gaussian distribution of scale-location mixture as

S (s) = \sum_{i = 1}^{Q} ω_{i}^{2} [N (s; μ_{i}, σ_{i}^{2}) + N (- s; μ_{i}, σ_{i}^{2})] / 2

(16)

where Q denotes the number of mixtures and the weight coefficient $ω_{i}$ is used to measure the proportion of each component. Note that the spectral density S is even symmetric.³⁷

Successively, by means of performing the inverse Fourier transform of equation (15), a corresponding spectral density complex kernel (SDCK) can be constructed as follows

k_{SDCK} (τ) = \sum_{i = 1}^{Q} ω_{i}^{2} \exp (- 2 π^{2} τ^{2} σ_{i}^{2}) \cos (2 π τ μ_{i})

(17)

Furthermore, equation (16) can easily be extended to higher dimensional inputs, $x \in R^{P}$ , which is easy to deal with large datasets and multidimensional inputs

k_{SDCKP} (τ) = Π_{k = 1}^{P} k_{SDCK} (τ_{k} | θ_{k})

(18)

where $τ_{k}$ is the $k th$ component of $τ = x - x'$ and $θ_{k} = {μ_{i}, σ_{i}^{2}, ω_{i}^{2}}_{i = 1}^{Q}$ are the $k th$ hyperparameters.

Compared with other popular covariance kernels (including the SK, MK, and RK), the SDCK with enough components Q enjoys many nice properties. First, the SDCK can achieve a wider range of kernel approximation, which is flexible even with a small number of components Q. Furthermore, the SDCK is highly expressive, which can be used to extract the representative feature information from the complex data structure. Finally, the simplicity of SDCK is one of its strongest qualities, which means that the GPR inference remains simple and appropriate for large multidimensional datasets.

Optimal inference approach of IGPR

Since the GPR model is suitable for analyzing high-dimensional, small-sample, and nonlinear problems, therefore, with the increase of encoder data, the computation of the GPR will start to rise rapidly. In consideration of this limitation, an optimal inference approach of the GPR model, which exploits the eigenvalue decomposition (EVD) in the model,³⁸ is proposed for reducing the computations and improving its generalization ability on large multidimensional datasets.

It is remarkable to note that the covariance kernel matrix K is widely used to construct the GPR model and estimate the hyperparameters such as ${(K + σ_{n}^{2} I)}^{- 1}$ , $\log | K + σ_{n}^{2} I |$ . A natural approach is, therefore, to decompose the covariance kernel matrix K so as to improve the computational efficiency. This approach can be achieved by the EVD as

K = U \sum U^{T}

(19)

where $\sum = diag (λ_{1}, λ_{2}, . . ., λ_{n})$ is a diagonal matrix storing the eigenvalues and $U = [u_{1}, u_{2}, . . ., u_{n}]$ is an orthogonal matrix. The column vectors of U are the orthonormal eigenvectors of K.

Substituting the above formula into the GPR model, the ${(K + σ_{n}^{2} I)}^{- 1}$ and $\log | K + σ_{n}^{2} I |$ can be reformulated as follows

{(K + σ_{n}^{2} I)}^{- 1} = {(U \sum U^{T} + σ_{n}^{2} I)}^{- 1} = U {(\sum + σ_{n}^{2} I)}^{- 1} U^{T}

(20)

\log | K + σ_{n}^{2} I | = \log | U \sum U^{T} + σ_{n}^{2} I | = \sum_{i = 1}^{n} \log (λ_{i} + σ_{n}^{2})

(21)

where ${(\sum + σ_{n}^{2} I)}^{- 1} = diag (\frac{1}{λ_{1} + σ_{n}^{2}}, . . ., \frac{1}{λ_{n} + σ_{n}^{2}})$ .

Thus, the decomposition of covariance kernel matrix K will greatly reduce the amount of calculation, with $O (P n^{(P + 1) / P})$ operations and $O (P n^{2 / P})$ storage. We can easily evaluate the model distribution and marginal likelihood in equations (6) and (10) to perform exact inference and hyperparameters estimation. Compared with other machine learning models, the advantages of constructed IGPR using encoder signals for weak fault detection are mainly reflected in the following aspects.

Accurate IAS estimation

As previously mentioned, the Fourier transform of a single covariance kernel indicates that it is a Gaussian distribution centered at the origin, while the adoption of scale–location mixtures of Gaussians can be employed to recognize the weak fault information from the complex data structure. Besides, the computation of the IGPR model is significantly reduced and therefore easy to estimate IAS on the big data context. The integration of these merits may provide an accurate way for IAS estimation.

CI estimation

The IGPR model can denoise to enhance weak fault feature, but it can also set a suitable CI for outliers’ detection. Therefore, it may provide an adaptive tool for quantitative analysis the health condition of rotating machinery.

Effective generalization ability

As a flexible nonparametric method, the complexity of IGPR models will not increase sharply with the increase in the amount of available data, while its hyperparameters are automatically estimated by means of the minimum NLML, without the need for cross-validation. In addition, this model can also be flexibly adjusted by constructing different complex covariance kernels. As a consequence, the IGPR may improve the generalization ability of the standard GPR model for weak fault detection of rotating machinery.

In sum, the flowchart of the proposed method is presented as below (Figure 3).

Simulation analysis

To demonstrate the efficiency of the proposed IGPR method in weak fault diagnosis, some simulated studies are conducted in this subsection. In industrial applications, the measured encoder signals can generally be described in the following form³⁹

\begin{matrix} s (t) = p (t) + g (t) + m (t) + n (t) \\ = 360 f_{0} t + \sum_{j} B_{j} sinc (c_{0} (t - j T_{0})) \\ + \sum_{i} A_{i} \sin (2 π i f_{0} t + β_{i}) + n (t) \end{matrix}

(22)

In the above signal model, four influencing ingredients are considered. The first item p(t) denotes the accumulated angular position of a rotating machine with frequency f₀, as shown in Figure 4(a). In this work, a sinc(t) function is applied to simulate the weak fault effect, which is displayed in Figure 4(b), where B_j, c₀, and T₀ are the magnitude, oscillation frequency and fault period, respectively. As exemplified in Figure 4(c), the third item m(t) refers to the angular position fluctuations due to driving/load variation or gear meshing, which is commonly described by a family of sine waves with amplitude A_i and phase β_i for the ith harmonic. The last item n(t) in Figure 4(d) illustrates the noise disturbances resulting from measurement noise and quantization error. The detailed parameters of encoder signals are listed in Table 1. Adding the above four components together, the final simulated signal is presented in Figure 5. It is noteworthy that the signal to be analyzed manifests as a straight line, and the interested fault feature can hardly be detected due to the trending component and the noise disturbances.

Table 1.

Parameters for the simulated signal.

$f_{0}$ /Hz	$B_{j}$ /deg	$c_{0}$ /Hz	$T_{0}$ /s	A ₁/deg	β ₁/rad	A ₂ /deg	β ₂/rad	SNR/dB	$f_{s}$ /Hz
10	0.5	1280	0.19	0.3	0	0.2	π/3	−10	5000

To detect the weak fault embedded in the raw signal, the traditional central difference method (CDM) is first applied to transform those position series into IAS. Figure 6 illustrates the IAS obtained with CDM. As can be seen, the IAS signal is so noisy that the weak fault transients can hardly be detected. The primary reason is that the CDM is very sensitive to measurement noise and other interference.

For solving this issue, the local polynomial regression (LPR) method is then utilized in this work.⁷ It adopts a smoothing mechanism to suppress the influence of noise. However, the performance of LPR is highly dependent on the choice of different degree d and fitting points p. In order to seek the optimal parameters to suppress noise interference and improve fault resolution, a two dimensional (2-D) grid search strategy was presented as shown in Figure 7. To be specific, we randomly select an array as the value of the (d, p) from the 2-D grid, and then calculate Pearson’s correlation coefficient ρ between the regression curve and the noise-free signal. Finally, the parameters with maximum ρ are selected for LPF. It can be seen that the fault-related transients cannot be easily recognized even by 2-D grid searches. The underlying reason for this failure is that the polynomial basis functions lack a flexible signal approximation mechanism. Moreover, the LPR cannot provide a CI for the IAS estimation, which makes it difficult to detect the fault-related IAS outliers in an adaptive fashion.

Successively, to determine the CI of the IAS estimation, the traditional quantile regression (QR) is also employed for this work. The QR with different degree d, quantile τ, and bootstrap number n, denoting QR (d, τ, n), are adopted. Figure 8 illustrates the IAS CI estimation result. It can be seen that the interval estimation of QR, which fails to detect the weak fault, can only be used for statistical analysis. Just like LPF, it is difficult to adaptively capture the local transient characteristics of the signal.

For conducting a comprehensive comparison, the GPR with conventional kernels is first employed to signal denoising and enhance weak fault feature as shown in Figure 9. It can be seen that although the GPR with conventional kernels has a better effect on interval estimation, the weak fault information cannot be easily recognized according to the CI. Moreover, by checking the feature enhancement results, we also observe that the defect impulses are still noisy.

Finally, the proposed IGPR is performed on the same task. As displayed in Figure 10, the IGPR can not only achieve accurate IAS estimation but also provide a CI for fault identification. It can be seen from this figure that 11 defect impulses spaced by 0.19 s are clearly observed, which is well agreed with the simulated fault impulse location. Furthermore, the weak fault can be clearly detected by means of the CI.

For quantitative evaluation of the performance of the proposed method in IAS estimation, the Pearson’s correlation coefficient ρ, and correlated Kurtosis (CK),⁴⁰ of different kernels are presented in Table 2. It can be seen that the IGPR reaches the highest values whether ρ or CK, which indicates whose IAS estimation accuracy is far larger than other kernels. Moreover, in order to testify the suitability of interval obtained by IGPR, the coverage probability CP_α,⁴¹ and mean width percentage MWP_α,⁴² of different methods are shown in Table 3. In this table, IGPR and GPR have close MWP_α, but IGPR has a larger CP_α. The reason is that the interval estimation mechanism of IGPR and GPR is the same but IGPR achieves a higher IAS estimation accuracy. Compared with GPR, QR has a lower CP_α and higher MWP_α; such an estimation accuracy can hardly provide any useful information for diagnosis. These metrics manifest that the CI obtained by IGPR is more effective.

Table 2.

IAS estimation accuracy of different kernels.

Kernel	Period	RQ	SE	Matérn	SDCK
ρ	0.132	0.218	0.264	0.406	0.891
CK	0.006	0.014	0.058	0.135	0.541

SDCK: spectral density complex kernel.

Table 3.

Interval estimation metrics of different methods.

Method	QR	GPR (period)	GPR (RQ)	GPR (SE)	GPR (Matérn)	GPR (SDCK)
CP 95%	0.721	0.737	0.796	0.813	0.841	0.95
MWP 95%	0.647	0.427	0.378	0.430	0.364	0.324

QR: quantile regression; GPR: Gaussian Process Regression; SDCK: spectral density complex kernel; MWP: mean width percentage, which is defined as the mean percentage of the interval width to the observation to ensure the validity of the interval.

Experimental validation

Planetary gearboxes (PG) usually work in the harsh environment and are prone to damage under instantaneous impact and alternating load excitation. If the early weak fault can not be detected and eliminated in time, they will gradually expand, increase, and eventually lead to the failure of the gearbox. To detect the early weak fault, an ESA-based health monitoring of PG has achieved relevant research. However, due to weak fault response, the fault feature is usually susceptible to noise interference. Furthermore, because of the kinematics of PG it is far more complex than that of fixed-axis ones, thus posing difficulty for accurate detection. To overcome this limitation, the IGPR is introduced in this segment, and its effectiveness in early weak fault diagnosis is also demonstrated.

As shown in Figure 11, the experimental device is powered by a servo motor that transmits torque from the PG and ends with a magnetic brake for load control. The specific parameters of the experimental device are shown in Table 4.

Table 4.

The specific parameters of the experimental device.

Parameter item	Parameter value
Rated power of servo moto	1.2 kW
Planetary gearbox ratio	5.1
Encoder model	Heidenhain ERN180

Among them, the system of the PG is shown in Figure 12, which includes a sun gear, an inner ring gear, and three evenly distributed planet gears. Table 5 lists the detailed parameters of the PG, where the sun gear is connected with input shaft, the gear meshes to output power from inner ring gear, and rotary encoders are installed on output and input shafts to capture operating information of the gearbox. In this experiment, the rotating speed of driving motor was 600 r/min, and the output encoder signals were collected with a sampling frequency of 5000 Hz using an IK220 counter card.

Table 5.

The detailed parameters of the PG.

Parameter item	Number of teeth	Fault characteristic frequency
Sun gear	$z_{s} = 20$	$f_{s} = f_{m} / z_{s} \cdot P$
Planet gear	$z_{p} = 31$	$f_{p} = f_{m} / z_{p}$
Inner ring gear	$z_{r} = 82$	$f_{r} = f_{m} / z_{r} \cdot P$

$f_{m} = f_{i} \cdot z_{r} z_{s} / (z_{r} + z_{s})$ is the PG fault characteristic frequency, and $f_{i}$ represents the rotation frequency of the servo motor.

Detection of tooth surface wear

Due to the harsh service environment of PG, tooth surface damages frequently occur in real applications, such as tooth surface wear, root crack, and tooth breakage. Among them, tooth surface wear is one of the most common early weak faults. To simulate this type of damage, two mild surface wears spaced by 13 teeth were seeded on the tooth surface of a planet gear as shown in Figure 13.

Figure 14 illustrated the encoder signal measured at the output shaft of PG. Similar to the simulated signal, the raw encoder signal is a straight line from which no valid information can be obtained.

To extract fault information from the raw signal, the CDM is employed for obtaining IAS and IAA as shown in Figure 15. It can be seen that the weak fault information can hardly be detected due to measurement noise interference.

For the purpose of signal denoising, the widely used fast-Kurtogram analysis is applied to this work.⁴³ By extracting the suitable resonant bandwidth, the interferences coming from noise interference could be effectively suppressed. Figure 16 depicts the fast-Kurtogram analysis method for IAA signal. To be specific, the optimal resonance band can be recognized as 1953–2265 Hz. Building on this band, the signal components of this frequency band are then extracted for fault diagnosis. It can be seen from Figure 16(b) that this signal contains only a series of random transients, rather than equally spaced impulses. Finally, a narrowband envelope method is employed for this filtered signal. However, it is hard to see effective fault feature information from Figure 16(c) and (d). The main reason is that fast-Kurtogram analysis is susceptible to random shocks. And the resonance frequency caused by the two weak faults is the same and, therefore, difficult to identify the different faults on the same gear.

Hence, to deal with the aforementioned issue, the wavelet packet transform (WPT) is used for the same task.^44,45 Figure 17 illustrates the IAA decomposition results obtained by WPT at depth 4 with db10 wavelet packets. In order to detect the fault transient features, we calculate the kurtosis of the signal at each node by WPT. The decomposed signals of the first three nodes with the highest kurtosis are shown in Figure 17(b) to (d), from which no obvious equally spaced impulses can be easily recognized.

In order to effectively identify the transient impact caused by fault, the GPR model based on the SDCK is used for the encoder signal. It is worth pointing out that there exist two transients per revolution of planet gear, as presented in Figure 18. Furthermore, the interval of the in-between transients is 0.08 s, which exactly corresponds to 13 tooth meshes. These results again demonstrate the effectiveness of the IGPR in weak fault detection.

Detection of tooth surface pitting

Tooth surface pitting is always inevitable in the long-term of PG. As the pitting extends into pieces, it will cause the metal block on the tooth surface to peel off, which may cause equipment failure. Thus, the early detection of pitting is critical to avoid fatal accidents. To simulate this kind of fault, a slight pitting was seeded on the tooth surface of a planet gear as given in Figure 19.

As exemplified in Figure 20, the encoder signal collected at the output shaft of the PG is a straight line, so it is difficult to observe any distinctive features. Furthermore, due to the low kinetic energy of slight pitting, the fault signature obtained with CDM is easily overwhelmed by measurement noise as illustrated in Figure 21.

With respect to denoise, the fast-Kurtogram is employed to detect the potential signatures from the IAA signal of slight pitting damage, and the consequences are depicted in Figure 22. It can be seen that little fault information could be acquired. The main reason for this loss is that kurtosis is merely used to evaluate the peakedness of the signal and therefore tends to emphasize the random shocks that are fault-unrelated components.

As a reference, the WPT is used for the same work. Figure 23 illustrates the decomposition results by WPT at depth 4 with db10 wavelet packets. Similar to the above experimental analysis, due to the poor IAA estimation, it is not easy to extract effective weak fault feature by means of WPT.

Figure 24 presents the result of IAS estimation and its interval confidence by IGPR. It can be seen that a series of periodic transients are clearly identified. Besides, those transients are spaced by 0.19 s, which exactly corresponds to the fault period of the planet gear. All of this evidence indicates that the proposed method is effective in detecting the weak fault.

Conclusion

A novel approach IGPR is proposed for weak fault detection of rotating machinery in this article. To extract weak fault feature information, a GPR model is first introduced to estimate the IAS and its CI. It has shown that the GPR model could not only achieve the IAS estimation like general regression methods, but also quantitatively analyze the accuracy of estimated values. In addition, the nonparametric mechanism can avoid the process of selecting hyperparameters. In order to improve the estimated accuracy of IAS, a new kernel termed SDCK is then constructed, which is to model the spectral density with a mixture of Gaussians. Built upon enough mixture components Q, the SDCK could achieve a broad class of covariance kernel approximation, which considerably improves the estimated flexibility. Finally, to reduce the computations of IGPR and improve its generalization ability on large datasets, an optimal algorithm scheme based on the EVD is constructed. The performance of the proposed IGPR has been validated by both simulated signals and experimental analyses. The outcomes indicate that the proposed method could not only achieve IAS estimation in accuracy and efficiency, but also set suitable CI for fault feature detection. Therefore, it may provide a useful tool for the health monitoring of rotating machinery in industrial applications.

Footnotes

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research is supported by the National Natural Science Foundation of China (Grant No. 51875434), Natural Science Foundation of Shaanxi Province (Grant No. 2019JM-278), the National Key Laboratory of Science and Technology on Reliability and Environmental Engineering, and the Key Laboratory of Advanced Manufacture Technology for Automobile Part (Chongqing University of Technology), Ministry of Education, which are highly appreciated by the authors.

ORCID iD

Ming Zhao

References

McDonald

Zhao

Multipoint optimal minimum entropy deconvolution and convolution fix: application to vibration fault detection. Mech Syst Signal Pr 2017; 82: 461–477.

Gao

Karimi

Kudreyko

, et al. Spare optimistic based on improved ADMM and the minimum entropy de-convolution for the early weak fault diagnosis of bearings in marine systems. ISA Trans 2018; 78: 98–104.

Lee

Zhao

, et al. Prognostics and health management design for rotary machinery systems—reviews, methodology and applications. Mech Syst Signal Pr 2014; 42: 314–334.

Wang

Zhu

, et al. Fault diagnosis of rotating machines based on the EMD manifold. Mech Syst Signal Pr 2020; 135: 106443.

Zhao

Lin

Health assessment of rotating machinery using a rotary encoder. IEEE T Ind Electron 2018; 65: 2548–2556.

Zhang

A new strategy of instantaneous angular speed extraction and its application to multistage gearbox fault diagnosis. J Sound Vib 2017; 396: 340–355.

Zhao

Jia

Lin

, et al. Instantaneous speed jitter detection via encoder signal and its application for the diagnosis of planetary gearbox. Mech Syst Signal Pr 2018; 98: 16–31.

Guo

, et al. A novel contactless angular resampling method for motor bearing fault diagnosis under variable speed. IEEE T Instrum Meas 2016; 65: 2538–2550.

Harris

, et al. The measurement of instantaneous angular speed. Mech Syst Signal Pr 2005; 19: 786–805.

10.

Jia

Zhao

, et al. Sparse filtering with the generalized l_p /l_q norm and its applications to the condition monitoring of rotating machinery. Mech Syst Signal Pr 2018; 102: 198–213.

11.

Diamond

Heyns

Oberholster

AJ.

Online shaft encoder geometry compensation for arbitrary shaft speed profiles using Bayesian regression. Mech Syst Signal Pr 2016; 81: 402–418.

12.

Luo

Fang

Ertas

Hilbert transform and its engineering applications. AIAA J 2009; 47: 923–932.

13.

Sandoval

De Leon

PL.

Theory of the Hilbert spectrum, 2015, pp. 1–61, https://arxiv.org/abs/1504.07554

14.

Zhao

Lin

Detecting weak position fluctuations from encoder signal using singular spectrum analysis. ISA Trans 2017; 71: 440–447.

15.

Jing

Wang

Zhao

, et al. An adaptive multi-sensor data fusion method based on deep convolutional neural networks for fault diagnosis of planetary gearbox. Sensors 2017; 17: 414.

16.

Seeger

Gaussian processes for machine learning. Int J Neural Syst 2004; 14: 69–106.

17.

Cai

Jia

Feng

, et al. A combined filtering strategy for short term and long term wind speed prediction with improved accuracy. Renew Energ 2019; 136: 1082–1090.

18.

Exterkate

Model selection in kernel ridge regression. Comput Stat Data Anal 2013; 68: 1–16.

19.

Pe’er

Bayesian network analysis of signaling networks: a primer. Sci Signal 2005; 2005: pl4.

20.

Wilson

Salakhutdinov

, et al. Deep kernel learning, 2015, http://proceedings.mlr.press/v51/wilson16.pdf.

21.

Wilson

Adams

. Gaussian process kernels for pattern discovery and extrapolation, 2013, p. 28, https://arxiv.org/abs/1302.4245

22.

Neal

. Bayesian learning for neural networks (lecture notes in statistics). New York: Springer Verlag, 1996, p.94724.

23.

Maddox

Garipov

Izmailov

, et al. A simple baseline for Bayesian uncertainty in deep learning, 2019, https://papers.nips.cc/paper/9472-a-simple-baseline-for-bayesian-uncertainty-in-deep-learning

24.

Wilson

AG.

Covariance kernels for fast automatic pattern discovery and extrapolation with Gaussian processes. PhD thesis, University of Cambridge, Cambridge, 2014, p. 226.

25.

Wang

Pleiss

Gardner

, et al. Exact Gaussian processes on a million data points, 2019, https://arxiv.org/abs/1903.08114

26.

Al-Shedivat

Wilson

Saatchi

, et al. Learning scalable deep kernels with recurrent structure. J Mach Learn Res 2016; 18: 1–37.

27.

Athiwaratkun

Wilson

AG.

Multimodal word distributions, 2017, https://arxiv.org/abs/1704.08424

28.

Oliva

Dubey

Wilson

, et al. Bayesian Nonparametric Kernel-Learning, 2015, p. 41, http://proceedings.mlr.press/v51/oliva16.pdf

29.

Poloczek

Wilson

, et al. Bayesian optimization with gradients, 2017, pp. 1–17, https://papers.nips.cc/paper/7111-bayesian-optimization-with-gradients.pdf

30.

Plataniotis

Hatzinakos

Gaussian mixtures and their applications to signal processing. In: Stergiopoulos

(ed.) Advanced signal processing handbook: theory and implementation for radar, sonar, and medical imaging real-time systems. Boca Raton, FL: CRC Press, 2017, pp. 3-1–3-36, https://doi.org/10.1201/9781315149790

31.

Wilson

AG.

A process over all stationary covariance kernels, 2012, pp. 1–6, https://www.cs.cmu.edu/∼andrewgw/spectralkernel.pdf

32.

Wilson

AG.

The change point kernel, 2013, https://www.cs.cmu.edu/∼andrewgw/changepoints.pdf

33.

Wilson

Ghahramani

Modelling input varying correlations between multiple responses. Lecture Notes in Computer Science (Including Subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics) 2012; 7524: 858–861.

34.

Rasmussen

Ghahramani

. Occam’s razor advances in neural information systems. Network-Comp Neural 2001; 5: 517–548.

35.

Wilson

Holland

, et al. Bayesian inference for NMR spectroscopy with applications to chemical quantification, 2014, pp. 1–26, https://arxiv.org/abs/1402.3580

36.

Lázaro-Gredilla

Com

Rasmussen

, et al. Sparse spectrum Gaussian process regression Joaquin quiñonero-candela aníbal R. figueiras-vidal. J Mach Learn Res 2010; 11: 1865–1881.

37.

Rea

Book reviews. J Am Acad Dermatol 1985; 12: 749–750, https://doi.org/10.1016/S0190-9622(85)80183-3

38.

De Lathauwer

De Moor

Vandewalle

. A multilinear singular value decomposition. SIAM J Matrix Anal Appl 2000; 21: 1253–1278.

39.

Feng

Zuo

MJ.

Vibration signal models for fault diagnosis of planetary gearboxes. J Sound Vib 2012; 331: 4919–4939.

40.

Cheng

Chen

Mei

, et al. A novel blind deconvolution method and its application to fault identification. J Sound Vib 2019; 460: 114900.

41.

Zhang

Qin

, et al. Wind speed prediction method using shared weight long short-term memory network and Gaussian process regression. Appl Energ 2019; 247: 270–284.

42.

Jin

A wind speed interval prediction system based on multi-objective optimization for machine learning method. Appl Energ 2018; 228: 2207–2220.

43.

Zhang

, et al. Adaptive Kurtogram and its applications in rolling bearing fault diagnosis. Mech Syst Signal Pr 2019; 130: 87–107.

44.

Haidong

Junsheng

Hongkai

, et al. Enhanced deep gated recurrent unit and complex wavelet packet energy moment entropy for early fault prognosis of bearing ✩. Knowl Based Syst 2020; 188: 105022.

45.

Fan

Zuo

MJ.

Gearbox fault detection using Hilbert and wavelet packet transform. Mech Syst Signal Pr 2006; 20: 966–982.