Constructing a Reliable Health Indicator for Bearings Using Convolutional Autoencoder and Continuous Wavelet Transform

: Estimating the remaining useful life (RUL) of components is a crucial task to enhance reliability, safety, productivity, and to reduce maintenance cost. In general, predicting the RUL of a component includes constructing a health indicator ( 𝐻𝐼 ) to infer the current condition of the component, and modelling the degradation process in order to estimate the future behavior. Although many signal processing and data ‐ driven methods have been proposed to construct the 𝐻𝐼 , most of the existing methods are based on manual feature extraction techniques and require the prior knowledge of experts, or rely on a large amount of failure data. Therefore, in this study, a new data ‐ driven method based on the convolutional autoencoder (CAE) is presented to construct the 𝐻𝐼 . For this purpose, the continuous wavelet transform (CWT) technique was used to convert the raw acquired vibrational signals into a two ‐ dimensional image; then, the CAE model was trained by the healthy operation dataset. Finally, the Mahalanobis distance ( MD ) between the healthy and failure stages was measured as the 𝐻𝐼 . The proposed method was tested on a benchmark bearing dataset and compared with several other traditional 𝐻𝐼 construction models. Experimental results indicate that the constructed 𝐻𝐼 exhibited a monotonically increasing degradation trend and had good performance in terms of detecting incipient faults. The exponential function is an increasing function that uses the mean from the starting time to the current time; thus, it can effectively eliminate oscillation and enhance monotonicity. It can be seen by the naked eye that the modified 𝐻𝐼 s were smoother, and gradually increasing, while the degradation trends were effectively captured as well. The results indicate that the proposed method has a good performance in detecting the early bearing defects and abnormal bearing health conditions. Moreover, this method provides an 𝐻𝐼 that is well correlated with progressively increasing bearing degradations, and it can lead to better RUL prediction.


Introduction
Performance degradation, which is almost inevitable for mechanical equipment, results in machinery damage, severe financial losses due to replacement or repair work and machine downtimes, or even personnel injury. Thus, prognostics and health management (PHM) has emerged as an engineering discipline to improve availability, reliability, and safety of equipment. As a crucial task in the lifecycle monitoring of complex equipment, PHM is used to monitor the equipment condition and to design robust and accurate models in order to assess the health state of equipment, as well as to define appropriate maintenance strategies [1]. In recent years, improving PHM methods by the Industry 4.0 paradigm, such as digital twin and predictive maintenance, attracts the attention of researchers [2][3][4][5].
In a digital twin, a virtual counterpart of the physical system during its whole life is created, with abilities such as analyzing, evaluating, optimizing, and predicting [6]. Jinjian et al. [5] presented a digital twin model for rotating machinery to diagnose the unbalance faults, on the basis of the dynamic behavior of the rotor system and vibrational status monitoring. Fei et al. [4] proposed a new approach for PHM, driven by digital twin for complex equipment. In this approach, a fivedimensional digital twin model is constructed to identify the health conditions of wind turbine gearboxes. Dinardo et al. [7] proposed a prognostic approach to detect the incipient faults of rotating machines by means of their vibrational status monitoring. Yan et al. [8] presented a two-phase digital twin to diagnose the fault using a deep transfer learning method. In this approach, the trained knowledge of the deep neural network is transferred from the virtual space to the physical space for real-time monitoring and predictive maintenance.
Traditionally, a digital twin uses physical-based simulation tools to describe the current behavior of a system [3]. However, due to manufacturing tolerances and material variances, describing a complex system in a simulation environment usually contains a strong deviation from reality [9]. One solution is to obtain a digital representation of the expected behavior of the physical system directly from measured data [3]. For this purpose, the first step is to construct a (multi) digital health indicator ( ) that describes different aspects of the physical component state during the whole life of the component. This should represent the deviation between the initial conditions of the component and its actual conditions during lifetime [1]. This can be further used for remaining useful life (RUL) estimation by implementing statistical estimation techniques, such as exponential degradation model [10], particle filter [11], or Kalman filter [12]. Therefore, defining an appropriate and sensitive that reflects the deviation degree from normal health conditions is now a hot research topic in the RUL estimation field.
In general, constructing a can be performed in three steps: (1) signal acquisition, (2) signal processing, and (3) feature extraction [13]. Vibration measurement provides a very efficient way of monitoring the dynamic conditions of a machine, such as unbalancedness, misalignment, mechanical looseness, structural resonance, wear, and shaft bow. [14]. Developing each failure mode leads to varying system dynamic behavior, resulting in significant deviation in vibration patterns [15]. Vibration signals generated by the faulty component can be analyzed in the time domain [16], frequency domain [17], or time-frequency domain [18]. Using the time-domain techniques for feature extraction requires the recording of the time-series vibrations over a long period of time to obtain suitable parameters to reveal fault evolution. However, obtaining the necessary data for a complex equipment may be expensive or even impossible. Using frequency-domain techniques such as fast Fourier transform (FFT) are powerful diagnostic tools in stationary conditions [18]. Since the FFT is essentially an integral over time, it fails to do so for non-stationary data, which could result from intermittent defect or evolutionary faults [18]. To address the FFT limitation, time-frequency signalprocessing tools such as the short-time Fourier transform (STFT) [19], Hilbert-Huang transform (HHT) [20], Wigner-Ville distribution (WVD) [21], and wavelet transform [22] are introduced. The wavelet transform is a relatively new and powerful tool, able to perform a local analysis of a signal and revealing some hidden aspects of the data that the other signal analysis fails to detect [23]. In this work, the wavelet transform was selected for signal processing to detect changes in vibration signatures that are caused by the faulty components.
Once the raw signal is acquired and processed, feature extraction techniques should be employed to extract the representative features that are used for construction. Feature extraction methods could be roughly classified into model-based methods and data-driven methods [24]. Rodney et al. [12] obtained a bearing by fusion of vibrational signal variance from the time domain and Choi-Williams distribution from the time-frequency domain. Yaguo et al. [25] presented a method to extract multiple features from the vibrational signal with multiple signal processing techniques, and then these features are selected and weighted to form the new . In [26], the authors implemented the discrete wavelet packet transform to decompose the raw signal into different subbands, and the was extracted from each signal. Although model-based methods do work and achieve an extraction of an accurate , they still have two deficiencies: (1) Feature selection is heavily dependent on prior knowledge and diagnostic expertise. Moreover, it often focuses on a specific fault type, and thus it may be unsuitable for other faults [27,28]. (2) In real industries, acquired signals are usually exposed to environmental noises, and are transient and non-stationary. Therefore, signal processing technologies need to be employed to filter the collected signals, which can result in a loss of information [27,29].
Data-driven methods attempt to extract features from measured data using machine learning techniques. In recent years, deep learning has emerged as a powerful tool to extract the representative feature from the collected signals [30,31]. Different deep learning architecture, including convolutional neural network (CNN) [29,32], recurrent neural network (RNN) [33], autoencoder [27], and generative adversarial network (GAN) [34], are successfully used to extract features automatically. The greatest advantage of deep learning is that it needs no prior expert knowledge and represents more accurate features [35]. Until now, most studies employ deep learning methods in a supervised setting to extract features for classification problems. For this purpose, performance data of different degradation levels of a component are prerequisites to creating labeled healthy and unhealthy datasets. However, gathering different degradation-level data requires large failure data, which is not available in practice, especially for high-reliability component [36]. On the other hand, a recent review on the state of deep learning on PHM [37] revealed that studies from a healthmanagement point of view have been rather limited, largely due to the unavailability of fault data. Moreover, many implementations of deep learning models in the literature are still constrained to specific equipment or applications and are not reusable when the predefined conditions change. In order to address the aforementioned restrictions, developing a single framework that can systematically be extended to all aspects of system health management is necessary [3,36]. This framework has to be able to be trained on-line without requiring historical data, and must use only healthy operational data for training [3]. In addition, it should be applicable to any equipment that operates under stationary and non-stationary conditions; it should also be extendable for different components [3].
As a step toward the development of a single framework for system health management, this paper proposes a method to construct an from the vibrational signal, on the basis of unsupervised deep learning. This method establishes an online construction of in the sense that the input data can be acquisitioned while the equipment is being exploited. The proposed method mainly includes three steps: First, healthy raw vibrational signals of the equipment are processed with the continuous wavelet transform (CWT) technique. These 2D images are considered as input of the deep learning model. In the second step, a convolutional autoencoder (CAE) model is developed that is solely trained by the healthy data. Lastly, during online monitoring, in each assessment interval, throughout the entire lifetime of the equipment, the CWT image of the vibrational signal is fed to the trained CAE model. Similar to the training data, the trained autoencoder can reconstruct images with small reconstruction errors. The distance between the normal condition data and the failure stage is measured by the MD formula and, thus, the is created. In this study, to experimentally evaluate the effectiveness of the proposed methods, we chose the ball bearing. Ball bearings are known as the most widely used rotating machine components, playing an important role in successful and reliable operation of rotary machines. Health prognostic of the ball bearing has great practical significance in reducing the failures of rotating machinery and enhancing machine availability. Thus far, fault detection techniques exist for rolling bearing monitor vibration, acoustic emission, motor current consumption, temperature, and oil debris. Among these techniques, vibration monitoring has proved to be a reliable and effective technique for fault detection in bearings [28]. Therefore, the vibrational analysis technique is selected in this work, and the CAE model is used to extract features from the vibration data. Overall, this study proposes a method to construct a on the basis of an unsupervised deep learning method that describes every instant condition of the bearing and can be regarded as an indicator for a digital twin. In brief, the main contributions of the current work are (1) The CAE model is only trained by using healthy operation data at the beginning of an asset's life cycle. Therefore, unlike most methods to construct a , this model can be trained online without requiring historical failure data from similar assets or fleets. In addition, since the CAE model is trained by the CWT image, it is applicable for equipment that operates under stationary and non-stationary conditions.
(2) The values of the bottleneck nodes of the CAE model are used as extracted features. Using these values reduces any dependencies on the prior knowledge, and thus the is constructed automatically.
The further course of the paper is organized as follows: Section 2 briefly introduces the theoretical background of CWT and convolutional networks. Section 3 presents the proposed methodology to construct the in detail. In Section 4, the results of the experimental evaluation are presented and discussed. Section 5 provides the conclusions and future work guidelines.

Continuous Wavelet Transform
The purpose of CWT of the raw vibrational signal is to preprocess raw vibration in the timefrequency domain and convert a 1D signal to a 2D image, as the input of the CAE model. The wavelet transform is widely used to process non-stationary signals over many different frequencies. The wavelet transform can analyze a localized area of a large signal without losing the spectral information contained therein. Therefore, the wavelet transform can reveal some hidden aspects of the signal that other techniques fail to detect. This property can particularly be employed to identify the damage (crack) or fault of a component that evolves during the time. There are two main trends in how wavelet transforms are used, the CWT and the discrete wavelet transform. Both Fourier transform and CWT use inner products to measure the similarity between a signal and an analytic function. In the Fourier transform, the analytic function is complex exponentials ( ) and in the CWT, the analytic function is a mother wavelet function, . The mother wavelet, ∈ , is a function of finite length and zero average; is the space of square-integrable complex functions [32]. The CWT compares the signal to shifted and compressed or stretched versions of the mother wavelet function. Stretching or compressing a function is collectively referred to as dilation or scaling and corresponds to the physical notion of scale. The family of time-scale waveform is obtained by shifting and scaling the mother wavelet, which can be expressed as By comparing the signal with the mother wavelet function at various scales and positions, we obtained two continuous variables, and ; is the dilation and is the translational parameter variable. For the given signal, , wavelet coefficient , can be represented as where * denotes the complex conjunction of . Since selecting of a mother wavelet function is application-dependent, the selection of the appropriate function is the first and most important step in the wavelet analysis. As a rule of thumb, the most appropriate mother wavelet is a function that has more similarity with the signal. Although there is no standard or general method to select mother wavelet, the shape matching by visual inspection is commonly used to select the appropriate mother wavelet function for the signal. For this study, on the basis of the visual inspection and the result of the previous studies [32,38], we selected Morlet wavelet in order to extract image features from the raw vibration signal. The Morlet function is a Gaussian function modulated by complex exponential, defined as where depends on the sampling frequency and usually is taken as 5 [39]. For wavelet transform of a real signal, the real part of the Morlet function is employed as the mother wavelet:

Convolutional Networks
Convolutional neural network (CNN) is a type of deep network that uses convolutional and pooling operation to extract the topological features of the input data. CNN is primarily used to solve difficult image-driven pattern recognition tasks, and if trained well, it will learn the features of the image completely. Therefore, in recent years, CNN has been widely used in image pattern recognition and image classification. CNN architectures come in several variations; however, a typical CNN includes convolutional layers, pooling layers, and fully connected layers. In the convolutional layer, a features map of the previous layer is convolved with multiple filters (also called Kernel) and is sent to the activation function to construct the output features map [29]: where is the th input feature map, * stands for the convolutional operator, denotes a convolutional filter, is an additive bias, is a set of input feature map, is the th layer in the network, and • is a nonlinear activation function. The mathematical inverse of the convolutional layer used in the decoder is known as the deconvolution layer. Different nonlinear functions such as a rectified linear unit (ReLU), sigmoid function, and scaled exponential linear unit (SELU) function can be used in convolutional layers. SELU is a variant of the ReLU activation function that, due to its self-normalizing properties, makes learning highly robust and allows training of networks which have many layers. SELU, ReLU, and sigmoid activation functions are defined as For standard scale inputs (zero mean and standard deviation), the selected values for the parameters are 1.6732 and λ 1.0507 [40].
The pooling layer usually follows the convolutional layer and is used to reduce the computational load by reducing the size of the features map. Two common pooling methods are max pooling and average pooling, which perform local max and average operations over the input features, respectively. The calculation process of the pooling layer is given as [29] • (9) where is the weight of pooling and is the additive bias, and denotes the downsampling function, e.g., max pooling. In contrast to the pooling layer, an upsampling layer is a simple layer in the decoder with no weights that is used to increase the dimensions of input.
In a fully connected layer, the features maps are converted into a one-dimensional feature vector, and all neurons of both layers are connected, like a traditional multilayer neural network. The output of the fully connected layer can be obtained as [29] (10) where is the output value; is the th neuron in the fully connected layer; and are the weight and the additive biases corresponding to and , respectively; and • is an activation function.

Methodology
In this work, a method is proposed to construct the automatically from the image of CWT of the vibrational signal of the ball bearing by using a deep learning model. The method consists of three main stages. The first stage involves acquiring and analyzing a healthy vibrational signal from the ball bearing and establishing the training repository for deep learning model. In this stage, it is assumed that the bearing is in a healthy condition and it is free from the defects at the beginning of its life cycle. In the second stage, the deep learning model is trained by the established healthy dataset. Finally, in the last stage, the is constructed to capture the bearing degradation throughout its failure phase. For this purpose, the difference in values of the bottleneck nodes between the failure stage and the normal stage is measured by using the MD formula for each assessment interval. The proposed method is shown in Figure 1. More details on these three stages are given below.

Data Acquisition and Analyzing
In the Industry 4.0 era, by emerging technologies such as data mining, internet of things, and cloud computing, the online data acquisition and processing becomes more pervasive than ever. The goals of data acquisition and data processing in this work are to create a dataset from the healthy conditions and to construct the throughout the failure stage. In the absence of sudden and unexpected failures, the degradation process of the bearing generally includes two stages: the normal operation stage and the failure stage, shown in Figure 2. In practical experiences, most of the bearing lifetime is passed in a stable and healthy stage; therefore, it is possible to acquire sufficient healthy data to train the deep learning model. During a normal operation stage of the bearing, a sliding window is used to capture the vibrational signal, and for each window, the power spectrum of the CWT is used to convert a 1D vibrational signal into a 2D image. The transformed image contains both time and frequency domain information and can represent the non-stationary and transient evolution of the signal. In practical application, a component is evolved from the normal stage to the failure stage gradually (excluding sudden and unexpected failures) through a series of degradation states. In addition, there are high uncertainties about the ambient conditions and component properties. Therefore, defining a fixed failure threshold that clearly separates the normal stage from the failure stage is not feasible. To address this issue, in this paper, we introduced an adaptive failure threshold approach, as depicted in Figure 3. According to the theory of statistical process control (SPC), measured vibration signals under the normal operation follow a normal distribution (with the mean μ and the variance σ) [41]. With the transfer from the normal operation to the failure operation, the distribution pattern of the vibration signals might vary from normal distribution to unknown distribution and, consequently, the mean and variance change. In this paper, the Pauta criterion [42] was employed to determine the failure threshold. If the mean of the acquired vibrational signal amplitude for each assessment interval is within the 3 , 3 range, the stage is recognized as a normal stage. Otherwise, the measured data are recognized as abnormal. In order to obtain an adaptive failure threshold, we established a healthy dataset ̅ , ̅ , … , ̅ by collecting the mean amplitude of the vibrational data points in the normal stage. Then, a reference range 3 , 3 is computed by using the data. When a new datum ̅ is observed, the Pauta criterion is used to determine whether the new data belongs to the healthy dataset. If the ̅ falls within the 3 , 3 range, it is added to the healthy dataset, i.e., ̅ , ̅ , … , ̅ , ̅ , and the reference range is updated. Otherwise, the datum is recognized as an abnormal one and the original healthy dataset remains unchanged. If for more than, e.g., 500 consecutive assessment intervals the means of the vibrational amplitudes are not in the reference range, the failure threshold is recognized, and the collected healthy dataset is not changed anymore.

Convolutional Autoencoder Model
The main purpose of using a deep learning model in this work is a dimensionality reduction through the feature extraction, and thus the extracted features represent the conditions of every moment of the ball bearing. Among the developed deep learning models, the CAE was selected in this study. CAE is a type of autoencoder (AE) neural network that is used to extract hidden features from unlabeled images. CAE is characterized by having identical input and output sizes and is trained to predict the input in the output (ℎ , ). One CAE algorithm consists of three layers: the input layer, the hidden layer(s), and the output layer. The idea is that one or several hidden layers have lower dimensions than the visible layers (input and output), and thus the input information is reconstructed and compressed in the hidden layer(s). The hidden layer that contains the fewest nodes is known as a bottleneck. The bottleneck layer represents the maximum point of compression of the input data, which contains all necessary data to reconstruct the input data again. Therefore, the CAE is constituted by two main parts: an encoder that maps the input into the code, and a decoder that maps the code to a reconstruction of the original input in the output layer.
In the encoding section, some convolutional layers and pooling layers are stacked on the input image to extract hierarchical features. Then, all units in the last convolutional layer have been flattened to form a vector followed by a fully connected layer(s). The bottleneck layer usually has 2 neurons. Accordingly, the input 2D image (for this study, the input image was 60 × 60 RGB pixels) is transformed into a two-dimensional vector space (ℝ → ℝ ). To train the CAE model in an unsupervised manner, we used the mirror architecture of the encoder in the decoder section. Thus, fully connected layer(s) followed by some deconvolutional and upsampling layers are used to transform the embedded features back to the original image. The structure of the proposed CAE model is shown in Figure 4. After setting up the CAE, it is essential to optimize weights and biases by minimizing the reconstruction error, i.e., the loss function. Backpropagation algorithm [43] is used to compute the gradient of the loss function with respect to any weight and bias in the network. To make the output of the decoder as equivalent as possible to the input, the binary_crossentropy function is employed as the standard CAE loss function, and the Adam optimizer is used to optimize the loss function [43].

Construction of
Once the CAE model is trained by the wavelet power spectrum images in normal operation, should be constructed for the failure stage automatically. The constructed is expected to exhibit a monotonically increasing trend and should be robust to noise and stochastic fluctuations. The CAE model learns to extract distribution characteristics of the normal data through its deep structure, and to reproduce images similar to the training dataset with a small reconstruction error. Over the time in the failure stage, damage evolution of bearing leads to a more turbulent vibration pattern. Consequently, with the development of the degradation, high energies appear in the low-frequency bands and produce a different wavelet power spectrum image from the normal stage. When the wavelet power spectrum image of the failure stage is input to the trained CAE model, the dissimilarity between the extracted vector of features in normal stage images and the faulty sample image is estimated to achieve the corresponding degradation indicator. In this regard, the mapped features in the bottleneck layer are used to measure the distance between the normal stage and the failure stage by the MD formula. This strategy is used to construct during the failure stage in this work. Since the feature extraction and construction are performed automatically by the CAE model, prior knowledge and diagnosis expertise are not required.
The , which is an effective multivariate distance metric that measures the distance between a vector and a distribution, still remains to be defined. An anomaly vector is an observation that has a different distribution from the rest of the data, and MD has an excellent ability to identify this deviation [44]. The MD is defined as (11) where A = ( , , , . . , gives the values of the bottleneck nodes in each assessment interval, , , , … , is the mean values vector in healthy condition, and is the reverse covariance matrix of the healthy condition.

Dataset Description
To evaluate the effectiveness of the proposed method, we used the bearing dataset, generated by the University of Cincinnati Center for Intelligent Maintenance Systems (IMS). The IMS bearing test rig is illustrated in Figure 5. It consists of an AC motor coupled to the shaft via a rub belt and four double row bearings installed on the shaft. The sampling rate of the record data was 20 kHz under the radial load of 6000 lbs for each of the four bearings at the constant rotation speed of 2000 rpm. Data acquisition was made every 10 minutes and for each assessment interval; a 1-second vibrational signal snapshot, which included 20,480 points, was recorded. The IMS bearing datasets contain three run-to-failure tests; both normal stage and failure stage data exist in each test. For each test, data collection continued until any failure in inner race, outer race, or roller elements occurred at least for one bearing. At the end of the test, bearings 3 and 4 from the first dataset, bearing 2 from the second dataset, and bearing 3 from the third dataset showed signs of failure. In the current work, the datasets for all three cases were considered. The time-domain vibrational signals for the four bearings are shown in Figure 6. For each bearing datum, the failure threshold was identified by the adaptive failure threshold method, and healthy samples and faulty samples were established. Details of IMS bearing datasets are described in Table 1. Since different bearing defect frequencies are proportional to the revolutions per minute, all defects occurring in the bearing are revealed in every revolution. Therefore, given the acquisition frequency of 20 kHz and rotation speed of 2000 rpm, the vibrational signals were split into equal chunks of length 600 points. For each chunk, the CWT was performed, and the power spectrum image was obtained.

Analysis of Wavelet Power Spectrum Images
To demonstrate the advantage of the wavelet transform technique, we depicted the time-domain vibrational signals and their corresponding wavelet power spectrums for a normal and a failure stage in Figure 7. While periodic vibrations with low amplitude are observed in the normal stage, more severe vibrations appear in the failure stage. The wavelet power spectrum represents the variations of the energy distribution of the vibration signal for different frequencies over time. As shown in Figure 7, the wavelet transform can clearly distinguish between the normal stage and the failure stage signals. For normal operations, most of the energy is concentrated in high frequencies, but for a failure stage, due to the evolution of defects, the burst of energy is observed in a broader range.

Training the CAE Model
As mentioned above, power spectrum images of the normal health stage were used to train the CAE model. Before training the CAE model, some parameters needed to be configured such as the number of layers, number of nodes in each layer, activation function for each layer, iteration number, and the learning rate. Regarding the number of layers and the number of nodes in hidden layers, the parameters are related to the dimension of the input data. These parameters are usually obtained experimentally and provide a good visual conformity. The convolutional autoencoder model starts with propagating from the input layer to the convolution layer. The convolutional encoder consists of three convolutional layers, each followed by a max pooling layer. The number of filters for the three convolutional layers were (1,64,128), and the filter size for all layers was 3 3. Features map from the convolutional encoder section were flattened and passed through the autoencoder layers. The autoencoder section consisted of three fully connected layers with 64, 2, and 64 neurons in each layer. Therefore, the input 2D RGB image was transformed into a two-dimensional feature space in the bottleneck layer that was used to construct the . To transform the extracted features back to the original image, we added three deconvolutional layers with three upsampling layers in between to the autoencoder section. The number of filters for the three deconvolutional layers were [1,64,128], with the filter size of 3 3.
To attain better visual conformity, after several trials, we chose the SELU as the activation function for the convolutional and deconvolutional layers, ReLU for the autoencoder layers, and the sigmoid for the last layer. To illustrate the influence of the learning rate on training the model, we show the reconstruction error under various learning rates in. If the learning rate was too low, the convergence was too slow or overfitting; if the learning rate was too high, it would hinder the convergence. Therefore, the learning rate was set to 0.005 for this study. As illustrated in Figure 8, the reconstruction error did not significantly decrease after 30 cycles; therefore, 30 iterations were performed in this study. The experimental model was developed using the Python-based Keras library [46] with a TensorFlow backend [47]. Since the CAE model is trained in an unsupervised manner, the successful trained CAE should learn to extract meaningful features from the images of the normal operation, and to reconstruct a closely similar image to the original image in the output. However, perfect reconstruction is usually a sign of overfitting where it only learns to copy its input to the output without learning to extract intelligent features and generalize to a new instance. Indeed, reasonably close reconstruction with a small error demonstrates that the CAE learned the meaningful features of the training dataset and had an acceptable generalization. Figure 9 represents the comparison of the original images with the reconstructed images by the developed model. Although the CAE model was solely trained by the normal operation dataset, it also reconstructed faulty images very well, indicating that the model could extract subtle features from images.

Smoothing the
The preliminary designed almost always exhibits local random fluctuations. To reveal an underlying long-term trend in the designed , we should smooth local spurious fluctuation in the curve. In this study, an exponential function was used to remove any sharp changes in the curve and to improve the monotonicity of the designed [48]. This function was given by where denotes the historical measured , is the total number of values, and represents the modified value of the current . In Equation (12), the mean value of the historical measurement from the first value to the current th value is used to calculate the . Therefore, if the curve exhibits a significant oscillation at point, it will be weakened. Furthermore, the exponential function is a monotonically increasing function and it can reveal a monotonically increasing trend in the curve.

Results
Since the normal stage images are used to train the deep learning model, the is solely constructed for the failure stage. The trained deep learning model was used to construct on-line by applying the MD formula to measure the distance of the values of the bottleneck nodes between the normal stage and failure stage. The constructed s for the four IMS bearings are shown in Figure  10a-d. Intuitively, it can be observed that s evolved gradually at the beginning and dramatically at the end; they revealed the true degradation in bearings. Although the initial curves represented global monotonicity, there still were severe local spurious fluctuations that may have been the result of highly inaccurate and unreliable data. Therefore, the smoothness and monotonicity of the constructed s were improved by using the exponential function defined in Equation (12); the modified s are shown in Figure 10e-h. The exponential function is an increasing function that uses the mean from the starting time to the current time; thus, it can effectively eliminate oscillation and enhance monotonicity. It can be seen by the naked eye that the modified s were smoother, and gradually increasing, while the degradation trends were effectively captured as well. The results indicate that the proposed method has a good performance in detecting the early bearing defects and abnormal bearing health conditions. Moreover, this method provides an that is well correlated with progressively increasing bearing degradations, and it can lead to better RUL prediction.

Comparison with Other Traditional Methods
To verify the effectiveness and superiority of the proposed method, we considered a comparison among the proposed method and several traditional methods. These methods were categorized in the two separate groups: five methods that are based on the time-domain features and a method that is based on the time-frequency domain feature. To compare, the vibration signal of the "subset 2 bearing 1" was used to construct the in this section. (1) In the first group, several commonly used features from the time domain were extracted to construct the . The selected features include:  Root mean square (RMS), is the vibrational signal series; ̅ and are the mean value and the variance of the series, respectively.
Approximate entropy (ApEn). ApEn expresses the regularity of fluctuations over time-series data. The detailed steps to construct an using the ApEn are introduced in [49]. During the failure stage, these equations were applied for each assessment interval to construct . Figure 11a-e illustrates the extracted from these methods. The RMS and variance curves of the are shown in Figure 11a,b, respectively. It can be noticed that these curves were insensitive during the early stage degradation, making it difficult for RUL prediction. The ApEn and kurtosis curves depicted in Figure 11c,d overcame the weakness of RMS and variance curves and recognized the infant mortality period. However, in these curves, oscillations and sudden changes near the end of the life of the bearing were obvious. This sudden change in the curve may cause problems in predicting the RUL accurately. As can be seen from Figure 11e, the skewness curve contained severe noises, and no up-and-down trend was visible, especially during the end of failure period.
(2) In the second group, the empirical mode decomposition (EMD) process was applied to decompose the vibrational signal into a series of intrinsic mode functions (IMFs). Afterward, the concept of singular value decomposition (SVD) was used to compute singular values (SVs) from the first two IMFs, known as defect feature vectors. Finally, the extracted feature vectors were taken as the input of the K-medoids algorithm to clustering normal and abnormal conditions and constructing the . The details steps can be found in [50]. The constructed by the second method, EMS-SVD-K-medoids, is shown in Figure 11f. It is realized from Figure 11f that during the early stages, the curve had a smoothly increasing trend. However, after about 50 10 s, unpredictable stochastic fluctuations were obvious. To evaluate the performance of as defined in different methods, we employed three metrics, namely, correlation ( ), monotonicity ( ), and robustness ( ) for different s. It was expected that a good would exhibit a monotonically increasing or decreasing trend, and that it would be robust to noise and stochastic fluctuations. was used to assess consistently increasing or decreasing trend of the curve. In , the difference between the values of any two adjacent points of the curve is measured. For the rising monotonicity, the total number of positive values is more than the total number of negative values and is close to 1. On the other hand, for the turbulent and oscillation curves, the total number of positive values is close to the total number of negative values and the value is close to 0. is calculated as follows: . 0 1 where is the difference in the values of any two adjacent points in curve and is the total number of values. reflects the tolerance of the to random fluctuations, which may arise due to faulty sensors, variations in operating conditions, or unexpected events.
is defined as where is the mean trend value of the . Similarly, measures the degree of linear correlation between the and time. It is expected that a good gradually increases by time. In a strong positive correlation, the value is close to 1 and vice versa.
is defined as Here, is the mean value of all the values. Table 2 presents , , and values for the proposed and the traditional health indicators mentioned earlier. The results in Table 2 show that all , , and of the current model are higher than those in other models. The obtained results demonstrate that the proposed model is superior to other models, and it yields a better .

Discussion
Constructing a reliable is the first and most important step in order to estimate the accurate RUL for bearings; therefore, this has been the focus of many studies [24]. In general, these methods are classified into three categories: mechanical signal processing-based, model-based, and machine learning-based. In mechanical signal processing-based methods, after pre-processing of the vibration signal, statistical parameters are directly used to construct the . Due to the flexibility and simplicity of mechanical signal processing methods, these methods are widely used in industries. These methods also have an acceptable performance to detect early bearing defects and abnormal bearing health conditions. However, it has been experimentally shown that the indicator performance decreases in the presence of transient conditions caused by bearing's defects [1]. Compared to these methods, the proposed method is sensitive to initial degradation, and is consistent with the degradation process. Nevertheless, in this work, the data of the run to failure vibrations is divided into two parts: the first part is used to train the CAE model and the second part is used to construct the . However, nothing ensures that a sudden degradation or failure does not happen during the training phase. Therefore, the method proposed in this work is limited to those faults that cause particular vibration patterns. In the case of any sudden failure or extremely slow degradation, this method is not able to construct the .
In contrast to the time or frequency techniques, which only represent the information in time or frequency domain, time-frequency techniques provide more information in both domains. In the present work, the CWT technique was used to pre-process the vibrational signals. The CWT method is a joint time-frequency analysis method that can decompose a time series into time and frequency spaces simultaneously. Therefore, the outputs of the CWT analysis are images that contain information on both time and frequency domains. When a defect appears in the bearing, it generates an impulsive force and excites resonances in the bearing and surrounding elements. With the progress of the defect over time, the frequency spectrum changes drastically. Since the faulty signals are non-stationary and transient in nature, using the CWT for pre-processing the vibration signals has better performance than time or frequency techniques in constructing the . Furthermore, in the proposed method, the is constructed by comparing the images of normal and failure stages, which are acquired for an identical bearing. Therefore, the perpetual background noise will not affect the accuracy. In addition, in this work, a deep learning model is used in extracting features, as well as for dimensionality reduction from the pre-processed vibration signals. This provides a more powerful capability of learning complex nonlinear relationships, which is able to extract the bestsuited features automatically. Moreover, using the exponential function improves the smoothness and monotonicity of the preliminary designed , which leads to better RUL estimation.

Conclusions
A new data-driven approach to construct the is presented. This represents every moment conditions of the bearing and can be considered as a digital twin of the bearing during its failure stage. Furthermore, this can be used for RUL estimation. First, the Pauta criterion was employed to determine the failure threshold and a normal dataset. Since the CWT is suitable for analyzing the non-stationary signals, it was used to convert raw vibrational signals into twodimensional feature images. The wavelet power spectrum image clearly revealed the degradation process of the bearing and included information in both time and frequency domains. Subsequently, the CAE model was used for dimensionality reduction through the feature extraction, and it was trained by the normal operation dataset. The values of the bottleneck nodes of the trained CAE represented the conditions of every moment of the ball bearing life; they were used to construct . Finally, the wavelet power spectrum image of the failure stage was fed to the trained CAE. The distance between the values of bottleneck nodes in normal and failure stages was measured by MD formula, and then the was constructed. To improve the curve monotonicity, we used an exponential function to remove random fluctuations in the curve. Experiments were conducted on the run-to-failure IMS dataset to verify the performance of the proposed method.
The results indicate that the constructed was capable of representing the health status of the bearing and tracking the evolution of degradation over the whole lifetime of the bearing. Moreover, constructing the with the proposed method required no prior knowledge or failure history data. Therefore, it is suitable for industrial applications. Furthermore, to prove the effectiveness of the proposed method, we compared this method with several other methods, such as RMS, EMD-SVD-K-medoids, skewness, kurtosis, and ApEn, showing a considerable superiority. The method, at the current state, is limited to gradual degradation and excludes any sudden failure.
Future research directions are to use the proposed method for other mechanical components such as ball screws, gears, and cutting tools. Funding: This research was partially funded by Leading House South Asia and Iran seed money grant at ZHAW Zurich University of Applied Science.

Conflicts of Interest:
The authors declare no conflict of interest.