Residual life prediction of progressive failure bearings based on NGO-AVMD hybrid domain features

1 Introduction

Rolling bearings are vital components in rail transit vehicles, crucial for running gear and equipment like motors. Their condition directly impacts system safety and reliability. Unexpected failures cause increased downtime and maintenance costs, leading to significant economic losses (Sun et al., 2025). Thus, predicting their Remaining Useful Life (RUL) is paramount for equipment health assessment.

While bearing state analysis involves methods like vibration signals, temperature monitoring, and oil analysis, vibration analysis is the most suitable and widely used. This is because damaged bearings emit distinct high-frequency pulse excitations, and vibration signals are stable and easily collected. However, collected vibration signals are often complex, nonlinear, and non-stationary due to track vibrations and environmental noise. Sensors measure these acceleration data for fault information (Heng et al., 2009). While features extracted from these signals are crucial for RUL prediction, the presence of significant noise in initial signals makes discerning the bearing’s operational state challenging. Therefore, various signal processing techniques across time, frequency, and time-frequency domains are essential for extracting robust features.

Vibration signal processing methods typically include time-domain (Lei, 2016), frequency-domain (Kong et al., 2022), and time-frequency domain analysis (Wang and Du, 2020), with common techniques such as EMD (Meng et al., 2022), EEMD (Yu et al., 2021), CEEMDAN (Civera and Surace, 2021), and VMD (Kumar et al., 2021). While EMD improves fault detection for weak signals (Shi et al., 2020), it often encounters mode mixing and endpoint effects. EEMD (Zhijian, 2018) partially addresses mode mixing but remains complex. Jin et al. (2022) introduced VMD, which adaptively determines frequencies and modes for effective signal component separation. Recent studies have further substantiated the superiority of VMD in structural health monitoring. For instance, comparison experiments in demonstrated (Civera and Surace, 2021) that VMD significantly outperforms other common signal decomposition techniques, such as CEEMDAN and HVD, in terms of defect identification and feature extraction robustness. Liu et al. (2018) confirmed VMD’s enhanced noise removal and fault feature extraction for train wheelset bearings. Crucially, VMD’s decomposition effectiveness is significantly influenced by parameters like the number of modal components (K) and the penalty factor (α). Recent VMD parameter optimization efforts have utilized algorithms such as PSO (Tang and Wang, 2015), WOA (Mirjalili and Lewis, 2016), Genetic Optimization (Gu et al., 2022), Sparrow Search (Liu et al., 2022), and Multiverse Optimization (Sayed et al., 2019). Despite their contributions, these methods often exhibit challenges with weak convergence and limited search capabilities. To overcome these limitations, this study employs the Northern Goshawk Optimization Algorithm (NGO) to fine-tune VMD’s K and α. NGO is a novel optimization technique demonstrating superior convergence and search capabilities compared to other methods (Dehghani et al., 2021). Therefore, NGO is utilized for adaptive VMD parameter optimization, termed NGO-AVMD, to effectively determine both K and α.

In some of the existing studies, Su et al. (2012) obtained time-frequency domain features from bearings and employed principal component analysis (PCA) to reduce dimensionality. They then used the Least Squares Support Vector Machine (LSSVM) regression method to predict the remaining useful life (RUL). Guo et al. (2017) extracted six relevant similarity features and eight classic time-frequency features. The optimal set of ten features based on scoring metrics are selected. In the RUL prediction phase, these features were mapped to RUL labels using Recurrent Neural Networks (RNN). Ren et al. (2019) extracted time-domain, frequency-domain, and time-frequency domain features. They selected a 34-dimensional feature set as input for Gated Recurrent Unit (GRU) networks to predict RUL labels. Wang W. et al. (2022) extracted six time-frequency domain features from denoised empirical mode functions to represent the degradation of bearings. Some methods for estimating the RUL and predicting degradation trends have achieved good results by establishing a Health Index (HI). Wang H. et al. (2022) extracted and integrated degradation information from a subset of time-domain features to create an initial Health Index (HI) for predicting the Remaining Useful Life (RUL) of bearings. Xia and Lu Ql (2018) employed MD-CUSUM to reduce feature dimensionality and derive a monotonically increasing HI. Additionally, Li et al. (2019) computed a novel health index (HI) using convolutional integration techniques. In comparison to raw data, this new HI proved instrumental in precisely recognizing and quantifying the degree of bearing degradation. However, HI models constructed using single features or signal processing methods often contain limited information and may inadequately represent the bearing’s degradation process. Deep learning has found extensive application in HI construction. Ding et al. (2021) introduced a deep adaptive network for cross-domain regression, effectively aligning relevant subdomains between source and target domains, showcasing robust cross-domain generalization. Deep learning models are typically comprised of multi-layer neural networks, which possess powerful nonlinear modeling capabilities. However, the challenge with deep learning lies in its ability to generate complex and opaque features, which can hinder model interpretability. Consequently, explaining the feature extraction process within the model and establishing the foundation for health indicators can become a challenging endeavor. In addition, single-domain feature extraction typically concentrates solely on the information within a specific domain, disregarding other relevant information. This approach often leads to insufficient feature information, which fails to fully characterize the bearing state and accurately predict its lifetime. Multi-domain feature extraction methods primarily derive parameters from raw signals, including those related to time, frequency, and time-frequency domains. However, these methods typically capture information from only one scale space, failing to fully exploit the feature information needed to comprehensively describe the state of the bearing.

For predicting the Remaining Useful Life (RUL) of bearings, various neural networks, including CNNs and RNNs like Long Short-Term Memory (LSTM) and Gated Recurrent Units (GRU), have been widely used. Ren et al. (2018) introduced an innovative approach to predict bearing RUL using Deep Convolutional Neural Networks (DCNN). Yang et al. (2022) introduced a new approach for bearing RUL prediction using Uncertainty Quantification-based LSTM, which adaptively extracts features. Eknath and Diwakar (2023) introduced an innovative prediction approach that combines DCNN and GRU. In this study, a bidirectional recurrent layer was implemented to capture both historical data and future information, thereby enhancing the GRU model’s capacity for data assimilation. Wang et al. (2023) presented a Time Convolutional Network (TCN) incorporating a convolution attention mechanism. TCN is known for its superior prediction accuracy and rapid computational capabilities, rendering it well-suited for time-series challenges. It overcomes the constraints of non-parallelizable RNNs and similar networks, allowing for end-to-end Remaining Useful Life (RUL) prediction. Nevertheless, the approaches mentioned earlier primarily emphasize the temporal connections within vibration data, often neglecting the increased vulnerability of bearings to external interferences and the impact of mounting signal instability as bearings degrade progressively. The intricate variations in the amplitude and severity of vibration signals leads to significant errors when networks learn the characteristics of bearing degradation stages. Furthermore, most current predictions of bearing RUL are made without determining the failure mode. For instance, Ding and Jia (2022) proposed an end-to-end deep framework using cross-validation. However, since different bearings take unique failure modes throughout their lifespan, treating all bearings equally for RUL prediction without distinguishing their failure modes may induce errors. Bearings with different failure modes exhibit inconsistent distribution patterns, which impacts the accuracy of RUL prediction. Therefore, Kundu et al. (2019) categorized various failure behaviors of bearings into two categories: gradual degradation and sudden failure behaviors. Subsequently, Dong et al. (2023) conducted RUL predictions for bearings with various failure behaviors. They used transfer learning to classify the failure behaviors of the bearings, leading to improved accuracy in predictions.

This study addresses incomplete bearing degradation information and imprecise RUL prediction stemming from undifferentiated failure modes. We propose a hybrid-domain feature extraction method, combining frequency and time-domain features from unprocessed vibration signals with frequency, time, and time-frequency domain features from NGO-AVMD reconstructed signals. This hybrid set, used for RUL prediction, enriches feature diversity and comprehensively unearths degradation information by leveraging the distinct physical meanings of original and reconstructed signals. A robust health index (HI) is then constructed by integrating deep features from deep learning networks with multi-scale spatial information from the hybrid-domain vibration signals. This combination merges the complex pattern-capturing ability of deep features with the comprehensive statistical insights of hybrid-domain features, providing a more holistic description of bearing degradation by integrating global and local information. This HI is crucial for distinguishing bearing failure modes. Finally, an Multi-Head Attention Temporal Convolutional Network (MA-TCN) network predicts the RUL of progressively degrading bearings. Its multi-head attention mechanism enhances focus on time-step relationships, improves long-term dependency capture, offers robust representation learning, and enables parallel processing for computational efficiency. To further clarify the positioning of this study, we compare the proposed framework with existing state-of-the-art methods, such as Wang H. et al. (2022) and Wang et al. (2023). As summarized in Table 1, the primary distinctions lie in the comprehensive feature fusion strategy (combining physical and deep features), the advanced adaptive optimization (NGO), and the multi-head attention mechanism tailored for progressive degradation.

TABLE 1. Comparison between the proposed method and related works.

This paper offers the following key contributions:

• Due to the incomplete nature of degradation information in rolling bearings, the proposal to construct a multi-scale hybrid domain feature set involves extracting statistical features from both the initial vibration signals and the reconstructed signals after NGO-AVMD processing.
• The proposal involves combining deep features with hybrid domain features to calculate a health index that can differentiate between bearing failure modes. By incorporating them into the health index, the condition and degradation of the bearing can be more accurately represented, thereby improving the ability to evaluate the condition of the bearing. Based on the health index, bearing failure modes can be distinguished as either progressive degradation or sudden failure, which leads to improved accuracy in bearing RUL prediction.
• A Time Convolutional Network (TCN) with the integration of a multi-head attention mechanism is applied for predicting bearing life. The TCN, with its deep structure, is capable of capturing prolonged dependencies within time-series data.

The remainder of this paper is structured as outlined below: Section 2 presents the principles of the methods applied. Section 3 outlines the overall method flow employed in this study. The effectiveness and practicality of the proposed method are assessed with bearing datasets in Section 4. Finally, Section 5 offers a conclusion and future prospects.

2 Theoretical background

2.1 VMD

Obtaining vibration signals through accelerometers can introduce noise from rolling bearings, gears, and other mechanical components. Additionally, it can lead to noise in the measurement system. To eliminate this noise and better track degradation, signal processing techniques are employed. VMD is an adaptive decomposition method. During the process of obtaining decomposition components, this method determines the frequency center and bandwidth of each component through iterative search for the optimal solution of the variational model, thus achieving adaptive frequency domain segmentation of the signal and effective separation of each component’s Intrinsic Mode Function (IMF). The key to signal decomposition is the construction and solution of the variational model. The decomposition effect is influenced by the number of modal components and the penalty factor. Therefore, it is necessary to optimize both the modal components and the penalty factor. The core of VMD is to construct the variational problem and solve it through the alternating direction multiplier method. In this paper, we will briefly introduce the theory of VMD from two aspects.

First, VMD defines Intrinsic Mode Functions (IMFs) as amplitude-frequency modulation signals, which have limited bandwidth and center frequency. If the pre-decomposition number is K, the expression for the k-th IMF u k t is as follows in Equation 1:

u k t = A k t cos φ k t (1)

Here, k is a positive integer between 1 and K, φ k t represents a non-decreasing phase function, and A k t represents the instantaneous amplitude.

Next, the Hilbert transform is used on each IMF to extract the demodulated signal. Subsequently, the spectrum is shifted to the respective baseband by utilizing the estimated center frequency, as formulated in Equation 2.

δ t + j π t * u k t e − j ω k t (2)

Lastly, compute the square of the time-gradient L2 parameterization of the demodulated signal and estimate the bandwidth of each mode signal. Furthermore, a constrained variational model can be formulated as follows:

min u k , ω k ∑ k = 1 K ∂ t σ t + j π t * u k t e − j ω k t 2 2 s . t . ∑ k = 1 K u k = f (3)

In the equation: f denotes the initial time-domain signal; K stands for the count of Intrinsic Mode Functions (IMFs); ∂ t represents the time derivative at time t; u k t and ω k are the central frequencies of the k-th IMF obtained through signal decomposition; * denotes convolution calculation; σ t represents the unit impulse function.

In order to tackle the constraint variation issue in Equation 3, we introduce Lagrange multipliers λ(t) and a penalty factor α, leading to the following augmented Lagrangian function in Equation 4:

u k , ω k , λ = α ∑ k = 1 K ∂ t σ t + j π t * u k t e − j ω k t 2 2 + f t − ∑ k = 1 K u k t 2 2 + λ t , f t − ∑ k = 1 K u k t (4)

The Lagrangian function is resolved using a variable iterative optimization method with alternations, ultimately yielding the modulus u k t and the frequency band center ω k for the corresponding IMF.

2.2 Northern Goshawk Optimization (NGO) algorithm

Considering the earlier analysis, the decomposition layers (K) and the secondary penalty factor (α) exert a substantial influence on the outcomes of VMD decomposition. Incorrect parameter selection may lead to modal information loss or modal aliasing. To address this concern, we integrate the NGO algorithm into the VMD parameter optimization procedure to dynamically determine the optimal parameter combination. This eliminates the need for manual parameter setting based on empirical or prior knowledge, and we refer to this approach as Adaptive Variational Mode Decomposition (AVMD).

The Northern Goshawk Optimization (NGO) algorithm, proposed by Dehghani et al. (2021), models the hunting behavior of the northern goshawk, encompassing prey identification, attack, pursuit, and evasion activities.

In the Northern Goshawk Optimization algorithm, the hunting process of the northern goshawk can be categorized into two distinct phases: prey identification and assault (exploratory phase), followed by chase and evasion (exploitation phase).

2.2.1 Initialization

In the Northern Goshawk Optimization algorithm, the population of Northern Goshawks can be described using the following population matrix in Equation 5:

X = X 1 ⋮ X i ⋮ X N N × m = x 1 , 1 ⋯ x 1 , j ⋯ x 1 , m ⋮ ⋱ ⋮ ⋮ x i , 1 ⋯ x i , j ⋯ x i , m ⋮ ⋮ ⋱ ⋮ x N , 1 ⋯ x N , j ⋯ x N , m N × m (5)

Where: X denotes the population matrix of Northern Goshawks; X i represents the location of the i-th Northern Goshawk; x i , j stands for the position of the i-th Northern Goshawk in the j-th dimension.; N is the population size of Northern Goshawks, and m is the dimensionality of the problem being solved.

In the Northern Goshawk Optimization algorithm, the objective function of the problem determines the objective function value for each Northern Goshawk in the population. The objective function values of the Northern Goshawk population can be represented as a vector of objective function values in Equation 6:

F = F 1 ⋮ F i ⋮ F N N × 1 = F X 1 ⋮ F X i ⋮ F X N N × 1 (6)

Where: F represents the objective function vector for the Northern Goshawk population; F i denotes the objective function value of the i-th Northern Goshawk.

2.2.2 Phase 1: prey detection (exploration phase)

In the initial hunting phase, the Northern Goshawk randomly chooses a prey target and swiftly initiates an attack. As prey selection in the search space is randomized during this phase, it strengthens the exploratory capability of the NGO algorithm. In this phase, Northern Goshawks conduct a thorough search of the space to identify the optimal region, employing Equations 7–9 to describe their selection and attack behavior on prey:

P i = X k , i = 1 , 2 , … , N , k = 1 , 2 , … , i − 1 , i + 1 , … , N (7)

x i , j n e w , P 1 = x i , j + r p i , j − I x i , j , F P i < F i x i , j + r x i , j − p i , j , F P i ≥ F i (8)

X i = X i n e w , P 1 , F i n e w , P 1 < F i X i , F i n e w , P 1 ≥ F i (9)

Where: P i represents the prey position for the i-th Northern Goshawk; F P i is the prey’s objective function value for the i-th Northern Goshawk, K is a positive integer within the range [1, n]; X i n e w , P 1 is the updated position for the i-th Northern Goshawk; x i , j n e w , P 1 is the updated position in the j-th dimension for the i-th Northern Goshawk; F i n e w , P 1 represents the updated objective function value of the i-th Northern Goshawk after the first phase; r is a random number within the interval [0,1]; I is a random integer, taking a value of either 1 or 2.

2.2.3 Second phase: chase and escape (development phase)

In the aftermath of the Northern Goshawk’s attack on prey, the prey attempts to escape. Therefore, during the tail-end of the chase, the Northern Goshawks continue to pursue their prey. Because of their rapid pursuit speed, their exceptional speed enables them to consistently track and capture prey under various conditions. This behavior, when simulated, improves the algorithm’s local search capabilities in the search space. We assume that this hunting activity is approximated by an attack radius of Rdescribed by Equations 10–12:

x i , j n e w , P 2 = x i , j + R 2 r − 1 x i , j (10)

R = 0.02 1 − t T (11)

X i = X i n e w , P 2 , F i n e w , P 2 < F i X i , F i n e w , P 2 ≥ F i (12)

Where: t represents the current iteration number; X i n e w , P 2 represents the new position of the i-th Northern Goshawk; x i , j n e w , P 2 denotes the updated position of the j-th dimension of the i-th Northern Goshawk; T represents the maximum iteration number; F i n e w , P 2 represents the target function value of the i-th Northern Goshawk updated in the second phase.

2.3 Temporal Convolutional Network (TCN)

TCN is a novel time series model based on an improved 1-D Convolutional Neural Network (CNN). The core idea of this approach is to utilize causal convolutions to process time data and incorporate dilated convolutions to effectively handle long-range dependencies in time series models. Compared to traditional CNNs, TCN offers the following improvements: (1) It incorporates causal relationships between convolutional layers, enabling predictions of future information based on historical data and effectively reducing information loss. (2) Residual connection layers and dilated convolutions are used to increase the depth of the network, enabling it to have long-term memory capabilities, and providing advantages in handling time series data. The causal dilated convolution structure of the TCN is depicted in Figure 1 in this paper.

In TCN, the relationship between input and output is as follows in Equation 13:

Y = F T X (13)

Where: F T · represents the TCN model, X = {x 0, x 1, x 2,…, x t − 1, x t} is the input mixed-domain feature sequence at time t, and Y = {y 0, y 1, y 2,…, y t − 1, y t} is the corresponding output sequence.

In TCN, the performance gradually deteriorates as the network depth increases due to the introduction of causal dilated convolutions. This could be attributed to the following factors: Firstly, due to the deeper network structure, the training time cost is higher, and the network is more difficult to converge during backpropagation. Secondly, As the network structure becomes deeper, achieving constant mapping becomes more challenging, making model optimization more difficult. To address these issues, TCN introduces residual connection modules. Compared to constant mappings formed by fitting a linearly superposed network of layers, residual mappings offer enhanced optimization and expedite network convergence effectively. The structure of the residual connection is shown in Figure 2. The structure of the residual connection module is specified as follows in Equation 14:

x n + 1 = φ x n + Ϝ x n , ς n (14)

Where: x n represents the input sequence at the n-th level; x n + 1 represents the output of the residual block; Ϝ x n , ς n represents the residual mapping; usually consisting of two to three causal dilated convolution operations denoted as ς n, and φ · represents the activation function.

2.4 Multihead attention mechanism

The degradation information contained in the bearing vibration signal typically changes throughout its operational lifetime, and it often becomes more pronounced as degradation increases. The purpose of the attention module is to encourage the network to closely examine hidden knowledge in the time-series data that is closely related to the degradation process, while suppressing redundant information during the feature learning process. Reflected in the internal computations of the structure, an attention matrix is constructed by assigning different weights to each incoming degradation data, thereby achieving matrix reconstruction.

In general, the self-attention mechanism consists of a query matrix Q, a key matrix K, and a value matrix V. The definition of self-attention is as shown in Equation 15:

A t t e n t i o n Q , K , V = s o f t m a x Q K T d k V (15)

Where: Q, K, and V are linearly mapped formats, defined as Q = W q e i, K = W k e i, and V = W v e i; d k represents the dimension to project.

Multihead attention (Vaswani et al., 2017) projects input features into different subspaces through multiple self-attention operations defined in Equation 15, obtaining attention vectors for multiple subspaces, and then collecting these attention vectors together. This allows exploration of the correlations between different embeddings from multiple perspectives, enhancing the performance of the self-attention model. The definition of multi-head attention is shown in Equations 16–18. The structure has been shown in Figure 3.

M u l t i h e a d Q , K , V = C o n c a t h e a d 1 , … h e a d h (16)

h e a d i = A t t e n t i o n Q W i Q , K W i K , V W i V (17)

The result of the multi-head attention is concatenated with a single-layer feed-forward neural network (FFNN). The attention carrier is input into the single-layer feed-forward neural network, resulting in the output, which is defined as follows:

h a t t e n t i o n = s i g m o i d W f f n n a i + b f f n n + a i (18)

3 Our proposed approach

The detailed procedure is illustrated in Figure 4:

• Step 1: Gather vibration data from bearings.
• Step 2: Extract time and frequency domain characteristics directly from the raw vibration signals acquired by the sensors.
• Step 3: Optimize VMD using the NGO optimization algorithm to determine the optimal values for the decomposition levels, K, and the penalty factor, α. This approach is referred to as adaptive NGO-AVMD.
• Step 4: Apply NGO-AVMD to the raw vibration signals, filter out effective IMF components using Effective Weighted Sparse Kurtosis (EWSK), obtain the reconstructed vibration signals, and extract characteristics in the time, frequency, and time-frequency domain from the reconstructed signals.
• Step 5: Construct a mixed-domain feature set by combining the time-domain and frequency-domain characteristics extracted from the raw vibration signals with those extracted from the reconstructed signals.
• Step 6: Use the unprocessed vibration data collected by the bearing sensor as input to a Convolutional Autoencoder (CAE), extracting deep features from the bearing vibration signals. Combine these features with the mixed-domain features to calculate the health index.
• Step 7: Use the calculated health index to assess and classify the bearing’s failure behavior.
• Step 8: Use the mixed-domain feature set as input for the Multihead-TCN (MA-TCN) model for remaining useful life prediction.
• Step 9: Split the dataset into training and testing sets. Then, train the MA-TCN network using the training set. Optimize and adjust the network structure and parameters based on the training results.
• Step 10: Validate the effectiveness of the network by using the testing set, predicting the remaining useful life of the bearing, obtaining prediction results, and analyzing conclusions.

3.1 Construction of mixed-domain features

3.1.1 IMF component selection

The VMD method obtains various IMF components, each of which encapsulates local features at different time scales of the original signal. The initial IMF components effectively capture the primary characteristics of the original signal. To ensure that the constructed hybrid-domain features effectively retain the original signal’s characteristics while avoiding interference from noise and other components, some of the modal components contain several bearing vibration information, while others are noise interference. Therefore, in modal selection, kurtosis or correlation coefficients are usually used to identify whether the modal components are valid. However, a single indicator cannot fully explain the vibration properties of the signal. Hence, taking into account the signal’s sparsity, kurtosis, and correlation sparsity, the Effective Weighted Sparse Kurtosis (EWSK) index (Lv et al., 2022) is employed for modal component filtering. The advantage of the EWSK index is that it leverages the strengths of sparsity, kurtosis, and correlation coefficients. To judge the validity of the modal component, the specific criterion is to determine whether the calculated EWSK index value is positive. All IMF components satisfying EWSK >0 are considered effective and are retained. The final reconstructed signal is obtained by summing all these effective components. This selection strategy preserves the effective failure information from the original signal while efficiently removing noise components. The formula for calculating the EWSK index is as follows in Equations 14–23:

S p a = 1 N ∑ i = 1 N x 2 t 1 N ∑ i = 1 N x t (19)

K u r = 1 N ∑ i = 0 N − 1 x 4 t 1 N ∑ i = 1 N x 2 t 2 (20)

C o r = E x − x ¯ y − y ¯ E x − x ¯ 2 E y − y ¯ 2 (21)

W S K i = S p a · K u r · C o r (22)

E W S K i = W S K i − 1 K ∑ i = 1 K W S K i (23)

Where: Spa, Kur represent the sparsity and kurtosis of the signal x(t); N is the signal length; Cor is the correlation coefficient between the signals x(t); E (∙) represents the expectation operator, and i is the ith IMF component.

3.1.2 Feature extraction

In the time domain, 15 common time domain features are extracted. These features can provide insights into the changes and characteristics of vibration signals over time, such as the mean, variance, and peak value. In the frequency domain, extract 10 common frequency domain features, including spectral energy distribution and frequency peaks. In the frequency domain, these features are used to analyze the characteristics of vibration signals. In the time-frequency domain, 4 common time frequency domain features are extracted, such as energy and entropy. Energy reflects the magnitude of the signal’s energy and the trend of energy distribution over time. Entropy features represent the level of disorder in the signal, and common entropy features include permutation entropy, sample entropy, and fuzzy entropy. By extracting and analyzing these features, a more comprehensive understanding of the characteristics and properties of vibration signals can be obtained. The specific mixed-domain feature set was selected to leverage complementary physical interpretations across domains. Time-domain features, such as peak-to-peak and kurtosis, were chosen for their sensitivity to early-stage impulsive shocks, while frequency-domain metrics track energy shifts in resonance bands. Additionally, entropy-based time-frequency features quantify the increasing signal complexity associated with severe degradation. Although this high-dimensional set may introduce redundancy, explicit dimensionality reduction is avoided to preserve physical interpretability. Instead, the Multi-Head Attention mechanism within the MA-TCN model dynamically assigns weights to inputs, automatically focusing on discriminative indicators while suppressing redundant information during training. The selected features are shown in Table 2.

TABLE 2. Mixed-domain features.

The specific steps are as follows:

• Perform NGO-AVMD on the original vibration signal.
• Use the EWSK indicator to select the most effective IMFs for signal reconstruction.
• Calculate 15-dimensional time domain characteristics and 10-dimensional frequency domain characteristics for the original vibration signal. Additionally, calculate 15-dimensional time domain characteristics, 10-dimensional frequency domain characteristics, and 4-dimensional time frequency domain characteristics for the reconstructed signal.

Concatenate the extracted features sequentially to construct a multi-channel fused mixed-domain feature set.

In Table 2, t f 1-t f 11 are dimensional feature parameters, among which peak and peak-to-peak values reflect the magnitude of impact amplitudes on the rail surface in different states. RMS reflects the magnitude of axle box vibration energy, and its value changes with different defect levels, but it lacks sensitivity to small defect signals. The last four features, t f 12- t f 15 are dimensionless feature parameters, which exhibit good sensitivity to impulsive signals. f f 1-f f 10 are extracted through Fourier Transform (FFT), where y(k) represents the spectrum of the given time series, k = 1,2,3,…,K, and K is the number of spectral lines in the FFT spectrum, while f k is the frequency value of the k-th spectral line.

The time frequency domain characteristics introduce energy and entropy characteristics. Energy reflects the distribution trend of signal energy over time, while entropy features represent the degree of signal randomness. Commonly used entropy features include permutation entropy, sample entropy, and fuzzy entropy.

In the Northern Goshawk Optimization algorithm, the hunting process of the northern goshawk can be categorized into two distinct phases: prey identification and assault (exploratory phase), followed by chase and evasion (exploitation phase).

3.2 The calculation of health index

3.2.1 Convolutional Autoencoder (CAE)

A Convolutional Autoencoder (CAE) is an autoencoder that utilizes convolutional operations for feature extraction and data reconstruction. Unlike a regular autoencoder, which typically consists of a fully connected layer, the encoding part of a CAE includes convolutional layers. This enables it to capture local correlations and spatial structures within the input data. As a result, CAEs are known to perform better in feature extraction and can automatically learn important patterns, textures, edges, and other local features present in vibration signals from bearings. They are capable of extracting deep-level features from the data more effectively.

3.2.2 Health index

The deep features from the convolutional autoencoder are abstract representations learned from vibration signals via deep learning, lacking direct physical interpretation. To address this interpretability limitation, this paper draws inspiration from Dong et al. (2023). It integrates these deep features with the constructed mixed-domain features to calculate the health index. This integration effectively combines global and local vibration characteristics. The mixed-domain features provide overall trend and distribution information of the vibration signals. When merged with deep features, they offer a comprehensive reflection of vibration data characteristics in health index calculation. Deep features capture complex patterns within the signals, while mixed-domain features supply global statistical properties. This combined approach enables a more accurate representation of bearing degradation status, enhancing both health index accuracy and state assessment reliability.

The correlation coefficient can reflect the degree of correlation between data, serving as an indicator to measure the similarity between two sets of sample data. The higher the similarity, the larger the correlation coefficient. By calculating the correlation coefficient, it can be used to establish the health index (HI), which assesses the degradation process. The health index HI is defined as follows in Equation 24:

H I j = C o v X , Y S X S Y = ∑ i = 1 n x i − x ¯ y i − y ¯ ∑ i = 1 n x i − x ¯ 2 ∑ i = 1 n y i − y ¯ 2 (24)

Where X represents the feature vector of the baseline sample, which is extracted from the first vibration data file collected at the start of the operation (t = 0). The length of this baseline segment corresponds to the single sampling duration (0.1 s containing 2560 points for the PHM2012 dataset) with no overlap. Y is the feature vector of other data (Y 1, Y 2, …, Y j), Cov(∙) represents sample covariance, S X is the standard deviation of sample X, S Y is the standard deviation of sample Y. “n” denotes the number of data points in the sample, and “j” represents the jth moment of the bearing operation.

In this paper, Concatenating the deep features extracted by CAE with the mixed-domain features established as sample data, we first select the normal state of the sample data and the test sample as two feature sets. Then the correlation between the two sample sets is calculated via using the correlation coefficient for forming the health index.

Before extracting deep features, the CAE network is trained. The aim in this training stage is to minimize the disparity between the ultimate output and the initial input data, allowing the encoder’s output to closely mirror the characteristics of the original input data. During the encoding process, two pooling layers and convolutional layers are applied, while the decoding process involves three deconvolution layers. In addition, for optimization during the error backpropagation process, the Adam optimizer is employed. The ReLU activation function is utilized to introduce non-linearity, enabling the model to learn complex vibration patterns. The network is trained by minimizing the reconstruction error between the input and the decoder output using the Adam optimizer. The specific structure is shown in Figure 5.

4 Experimental validation

4.1 Introduction to the PHM2012 dataset

The experiment utilizes data from the comprehensive dataset released at the 2012 IEEE International PHM Conference, collected during the PRONOTIA experiment. A detailed description can be found in reference (Nectoux et al., 2012), and the experimental platform is depicted in Figure 6. The platform mainly consists of a rotating module, loading module, and measurement module.

Three different datasets were collected: the first with a motor speed of 1800 rpm and a load of 4000 N, and the other two with motor speeds of 1625 rpm and 4200 N, 1500 rpm and 5000 N, individually. No artificial faults were introduced during the experiment, and all faults were a result of the natural degradation of the bearings, mirroring real industrial scenarios. The data was collected at a 25.6 KHz frequency with a 0.1 s sampling interval, and measurements were recorded every 10 s. In this paper, four bearings under Operating Condition 1 were selected for validation. These specific bearings were chosen because they exhibit clear progressive degradation patterns, which aligns with the scope of this study. Sudden failure cases were excluded from the prediction task.

To ensure a robust evaluation, a Leave-One-Out Cross-Validation strategy was employed. As detailed in Table 3, four independent experiments were conducted. In each experiment, one bearing was used as the testing set, while the remaining three served as the training set.

TABLE 3. Experimental data and operating conditions.

4.2 Construction of hybrid domain features

4.2.1 Validation of the NGO-AVMD method

4.2.1.1 Simulation signal verification

Constructed a synthetic signal (Lv et al., 2022) to validate the performance of the proposed NGO-AVMD approach. This synthetic signal was sampled at a rate of 1200Hz, comprising 1200 sampling points, and lasting for 1 s. It involves the following signal details in Equations 25–28:

s 1 t = 0.5 ⁡ sin 2 π · 42 t (25)

s 2 t = 1.5 ⁡ sin 2 π · 113 t t ∈ 0.3 n − 0.3 , 0.3 n − 0.2 , n ∈ N 0 t ∈ 0.3 n − 0.2 , 0.3 n , n ∈ N (26)

s 3 t = 0.3 ⁡ cos 2 π · 250 t (27)

S t = s 1 t + s 2 t + s 3 t + s 4 t (28)

In the equations: s 1 t represents the low-frequency signal; s 2 t represents the intermittent signal used to simulate fault pulses; s 3 t is the high-frequency signal, and s 4 t is random noise with a mean of 0 and a variance of 0.1.

The frequencies of the three fault pulse signals are 42 Hz, 113 Hz, and 250 Hz. Refer to Figure 7 for the illustration of the simulated signal construction process.

Initially, the EMD, EEMD, and CEEMDAN methods are applied to decompose the original simulated signal. The Signal-to-Noise Ratio (SNR) and Mean Squared Error (MSE) for each component are presented in Table 3. These metrics serve as common indicators for assessing the reconstructed signal quality.

Next, the proposed adaptive NGO-AVMD method in this paper is employed to address the issue of manually determining parameters such as the decomposition level in traditional VMD. The adaptive AVMD method automatically determines these parameters through optimization algorithms, eliminating the cumbersome process of manually tuning the parameters. Traditional VMD is sensitive to parameter choices, and different parameter combinations can yield different decomposition results. The adaptive AVMD method, by adaptively adjusting parameters, enhances the algorithm’s stability and produces more robust decomposition results.

According to the NGO-AVMD method used in this paper, the best combination obtained is the optimal [K,α] combination: K = 4, α = 503. Table 4 lists the Signal-to-Noise Ratio (SNR) and Mean Squared Error (MSE) for the first four components of EMD, EEMD, and CEEMDAN. These metrics are commonly used to evaluate the quality of reconstructed signals. It can be observed that the original signal is decomposed into four components with central frequencies of 42 Hz, 113 Hz, 250 Hz, and 411 Hz. Through decomposition, it was found that the original simulated signal was decomposed into three valid components and one high-frequency noise component. The method used in this paper achieved the highest SNR and the lowest MSE values. Figure 8 shows the results of the decomposition.

TABLE 4. SNR and MSE of the simulated signal.

4.2.1.2 PHM2012 data validation

Using NGO algorithm to optimize the parameter combination [K 0, α 0] for VMD, the objective of the optimization is discussed by Gai et al. in relation to the fitness function (Gai et al., 2020). The envelope entropy serves as an indicator of variation in the envelope signal, and the optimization objective is to minimize the envelope entropy, selecting the minimum value. The range for K is set to (Heng et al., 2009; Mirjalili and Lewis, 2016), and the range for α is set to [500, 4000] (Ni et al., 2022). Specifically, the population size of the NGO algorithm was set to 30, and the maximum number of iterations was set to 50. The search dimension was set to 2 to optimize the pair [K, α]. Taking Bearing1_1 as an example, a sample from the full-life experimental data of Bearing1_1 is selected as a test sample. During the optimization of the VMD model, six different optimization algorithms were used to obtain fitness convergence curves for comparison, as shown in Figure 8. Among them, GA stands for Genetic Algorithm, SSA stands for Sparrow Search Algorithm, PSO stands for Particle Swarm Optimization, MOV stands for Multiverse Optimization, and WOA stands for Whale Optimization Algorithm.

The curves in Figure 9 indicate that the Sparrow Search Algorithm and Genetic Algorithm have slower convergence rates, while the PSO algorithm is relatively faster. However, the Particle Swarm Optimization algorithm may suffer from the issue of getting stuck in local optima. The Whale Optimization Algorithm exhibits local optima, while the Multiverse Optimization Algorithm converges slowly and faces significant local optima problems. Compared to the aforementioned optimization algorithms, the NGO algorithm not only possesses the ability to effectively avoid local optima but also demonstrates faster convergence performance.

As depicted in Figure 9, the search results exhibit convergence after the 5th iteration. The optimal [K, α] combination obtained through the Northern Goshawk Optimization algorithm is K = 5, α = 638, which breaks down the signal into 5 IMF components. The marginal spectra of each mode resulting from the signal decomposition are presented in Figure 10.

The results of decomposing the sample signal using EMD, EEMD, and CEEMDAN are compared in Table 4.

According to Table 5, it can be observed that the SNR of IMF1 components for EMD, EEMD, and CEEMDAN are all lower than those for NGO-AVMD, and their Mean Square Errors (MSE) are higher than those for NGO-AVMD. As depicted in Figure 10, the NGO-AVMD algorithm effectively extracts the low-frequency components from the original signal, successfully addressing mode mixing and facilitating the separation of each component.

TABLE 5. Sample SNR and MSE.

4.2.2 Selecting effective components

Calculate the Weighted Sparse Kurtosis Index (EWSK) for each modal component, as indicated in Table 5. The initial two rows represent commonly used metrics: Signal-to-Noise Ratio (SNR) and Mean Squared Error (MSE). These metrics serve to assess the reconstructed signal’s quality. The criteria for judging the components are strictly based on the EWSK index. As shown in Table 6, only IMF1 exhibits a positive EWSK value, whereas the EWSK values for IMF2 to IMF5 are all negative. Therefore, according to the selection rule (EWSK >0), IMF1 is identified as the sole effective modal component in this specific case, while the others are considered spurious and discarded. The effective modal component is used for signal reconstruction, as shown in Figure 11.

TABLE 6. Selection of effective components.

4.2.3 Extracting mixed-domain features

Extract 15 time domain characteristics and 10 frequency domain traits from the primary vibration signal. Then, compute 15-dimensional time domain characteristics, 10-dimensional frequency domain characteristics, and 4-dimensional time-frequency domain characteristics from the reconstructed signal obtained after selecting the IMF components. These characteristics collectively constitute the mixed-domain feature set.

4.3 Health index results and analysis

Select the complete lifecycle original vibration signal data of Bearing 1 for CAE training. After extracting deep features from the original vibration signals using CAE, concatenate the extracted hybrid domain feature set and compute the correlation coefficient as the health index. To validate the discriminability of the extracted feature set, we employed the t-SNE algorithm for dimensionality reduction and visualization. As illustrated in Figure 12, the features from different degradation stages map to distinct clusters in the two-dimensional space. The clear separation between clusters indicates that the fused feature set, enriched by CAE-extracted deep features, effectively captures the evolving characteristics of the bearing degradation process, providing a robust foundation for the subsequent Health Index construction. The health indices extracted for Bearing1-1, 1-3, 1-4, 1-5, 1-6, and 1-7 are shown in Figure 13. Based on the points where the health index curves decline, the operating states of the bearing are divided into three stages: normal operation, degradation, and failure. Additionally, based on the health index curves, bearings can be classified into two failure modes. Bearings 1_1, 1_3, 1_4, and 1_7 fall under the gradual slow failure mode, whereas bearings 1_5 and 1_6 belong to the sudden failure mode.

The comparison of some bearings combining deep features with hybrid domain features is shown in Figure 14.

Based on the constructed health index graphs, it can be observed that the nonlinear health index, which includes hybrid domain features, exhibits more instability and spikes compared to the linear health index. Therefore, the linear health index proposed in this paper has been chosen. The trend of the linear health index is more pronounced than the health index constructed solely using deep features. Taking Bearing1-7 as an example, from the time-domain vibration image in Figure 15, it can be observed that the impact of the vibration signal exhibits a pattern of fluctuation. The health index constructed solely from deep features cannot capture the subtle nuances of this variation. However, when hybrid domain features are incorporated, the health index constructed in this paper can more effectively reflect the subtle changes in the bearing.

Therefore, the health index reveals that bearing degradation can be categorized into two failure modes: gradual degradation and sudden failure. Sudden failure bearings typically fail abruptly, while gradually degrading bearings deteriorate over time. Treating them indiscriminately as the same type of bearing may lead to inaccurate predictions of failure time. Gradual degradation of bearings may be predicted to cause earlier failure, while sudden failure of bearings may be predicted to occur later. Furthermore, there may be differences in the data between sudden failure and gradual degradation of bearings. Using mismatched data to build predictive models together could introduce compatibility issues. Sudden bearing failure can be caused by external factors. Since sudden failure events are often influenced by multiple interacting factors, including workload, vibration, temperature, lubrication conditions, external damage, and more, the complex combination of these factors makes accurate failure prediction challenging. Therefore, in this paper, we solely focus on predicting the remaining useful life of bearings that gradually degrade.

4.3.1 Health index evaluation metrics

To assess the effectiveness of the health index, it is important for a good health index to possess monotonicity and consistency. Two metrics were selected for evaluation, namely, monotonicity (Mon, Equation 29) and consistency (Tred, Equation 30).

M o n = d F > 0 T − 1 − d F < 0 T − 1 (29)

T r e d H , T = ∑ k = 1 K h k − H ¯ t k − T ¯ ∑ k = 1 K h k − H ¯ 2 ∑ k = 1 K t k − T ¯ 2 (30)

The experiments incorporating a linear health index with mixed-domain features are denoted as A. Experiments with only a linear health index using depth features are denoted as B. Experiments incorporating a non-linear health index (Maximum Information Coefficient, MIC) with mixed-domain features are denoted as C. The results are shown in Table 7.

TABLE 7. Health index evaluation metrics.

4.4 Life prediction results and analysis

4.4.1 Evaluation metrics

To assess the model’s performance, five performance metrics were chosen, encompassing Mean Squared Error (MSE), Mean Absolute Error (MAE), Coefficient of Determination (R2), Root Mean Squared Error (RMSE) and Performance Score, defined as follows in Equations 31–35:

M S E = 1 n ∑ i = 1 n y i ^ − y i 2 (31)

R M S E = 1 n ∑ i = 1 n y i ^ − y i 2 (32)

M A E = 1 n ∑ i = 1 n y i ^ − y i (33)

R 2 = 1 − ∑ i = 1 n y i ^ − y i 2 ∑ i = 1 n y i ¯ − y i 2 (34)

S c o r e = 1 − M S E ∑ i = 1 n y i ¯ − y i 2 / n − M A E ∑ i = 1 n y i / n (35)

In the equations: y i represents the true values. y i ^ represents the predicted values. y i ¯ represents the mean of the true values.

The reason for selecting these indicators is to comprehensively assess the performance of the model. Metrics such as MSE, RMSE, and MAE provide information about the model’s prediction errors. The closer their values are to 0, the better the prediction performance. R^2 indicates the degree of fit between the model and the data, with values closer to 1 indicating a higher level of performance. The score provides an intuitive measure of the model’s overall performance, with values closer to 1 indicating better performance. By considering these metrics collectively, we can enhance our comprehension of the model’s performance in the task of lifespan prediction and implement necessary improvements and optimizations.

4.4.2 Temporal Convolutional Network (TCN) parameter configuration

In order to achieve accurate lifetime predictions, it is necessary to design the structure of the TCN network. Experiments were conducted on TCN using the method of controlled variables to determine the network’s structure and parameters. The experimental categories are displayed in Table 8.

TABLE 8. Experimental grouping conditions.

The various network designs in Table 8 were individually trained, tested, and compared using the dataset, ultimately determining the optimal network structure. The specific hyperparameters were configured as follows: the convolutional kernel size is 3 × 3, and the dilation factors are set to (1, 2, 4, 8). The LeakyReLU activation function is applied to prevent the dying ReLU problem. Based on the optimization experiments, the number of attention heads is set to 8 to capture multi-subspace correlations.

To facilitate reproducibility, the detailed hyperparameters for the MA-TCN model architecture and the training process are summarized in Table 9.

TABLE 9. Detailed hyperparameters of the MA-TCN model.

4.4.3 Setting the number of heads for multi-head TCN (MA-TCN)

After determining the structure of the Temporal Convolutional Network (TCN), we conducted experiments to incorporate a multi-head attention mechanism and explored the optimal number of attention heads. The results are presented in Table 10. The experimental results are shown in Figures 16, 17.

TABLE 10. Number of multi-head attention.

For comprehensive cross-validation (Kuo et al., 2022), three out of the four progressively degrading bearings were selected sequentially to train the remaining life prediction model. Then, another bearing with progressively degrading conditions was used to test the model.

The experiments aimed to investigate the impact of varying numbers of multi-head attention on the prediction results for bearing remaining life. The evaluation metrics for the four bearings in the test set included MSE, RMSE, MAE, R^2, and Score. Lower values for the first three metrics and higher values for the last two indicate better performance. The results represent the averages of five experiments.

From Figures 16, 17, observing the scenario where the number of multi-head attention is configured as 8, the MSE, RMSE, MAE metrics are the lowest, while the R² and Score metrics are the highest. Therefore, the optimal hyperparameter for the number of multi-head attention is chosen to be 8.

Its prediction results are shown in Figure 18.

Based on the prediction results, it can be seen that Bearing1-1 and 1-3 have the best performance. This may be because in the gradual failure mode, the data patterns of Bearing1-1 and Bearing1-3 are most similar, while the data patterns of Bearing1-4 and Bearing1-7 may have slight differences compared to Bearing1-1 and Bearing1-3. However, the overall performance is good, and as summarized in Table 11, the coefficients of determination (R²) are all 0.98 or higher, indicating a strong fit.

TABLE 11. Metrics results.

4.4.4 Generalization validation on additional datasets

To verify the generalization capability under different operating conditions, we extended the validation to the XJTU-SY dataset (Condition: 2100 rpm, 12 kN) and a self-designed test bench (Condition: 2800 rpm, 10 kN). As shown in Figure 19 and Table 12, the model achieved an R 2 exceeding 0.98 on the XJTU-SY dataset. Furthermore, the test bench results (Figure 20; Table 13) maintained an R 2 above 0.90 across all tested bearings. These results confirm that the proposed NGO-AVMD and MA-TCN framework possesses robust transferability across varying speeds and loads.

TABLE 12. Evaluation of indicator results.

TABLE 13. Results of lifetime pediction indicators.

4.4.5 Ablation experiments

In order to better compare our method, we set up sensitivity analysis and conducted comparative experiments with other methods, using Bearing1-7 as an example. The results are shown in Figure 21 and Table 14.

TABLE 14. Metric comparison.

Comparing different methods as shown in Figure 18, it can be observed that TCN, without any data preprocessing, outperforms LSTM, CNN, and GRU in terms of accuracy. After adding the multi-head attention mechanism, there is a decrease in MSE, RMSE, and MAE, with an improvement of 0.02 in and performance score. When applying VMD decomposition to construct the hybrid domain, the prediction accuracy significantly improves compared to the case without data preprocessing. In this case, the proposed method yields the lowest MSE, RMSE, and MAE, achieving an R 2 of 0.98. The performance score also increases by 0.14 compared to only adding the multi-head attention mechanism and by 0.09 compared to only performing data preprocessing to construct the hybrid domain, indicating higher prediction accuracy.

5 Conclusion

In this study, a bearing life prediction approach is introduced, which leverages NGO-AVMD mixed domain features in conjunction with MA-TCN.

In this method, information from the initial vibration signals and the different scale information from the reconstructed signals after NGO-AVMD processing jointly constitute a more comprehensive mixed-domain feature. This feature is then used as the input for the MA-TCN prediction model. This allows for a more comprehensive understanding of the bearing’s overall condition.

In the calculation of the health index, a convolutional autoencoder is utilized to extract deep features from the original vibration signals at a different scale space. These features are then combined with the constructed mixed feature set to determine the health index. The inclusion of mixed-domain features in the health index provides a more comprehensive reflection of the bearing’s operating state, enabling a more thorough understanding of the bearing’s degradation process. This enhances the effectiveness of equipment monitoring. A comparison with the health index that does not include mixed-domain features validates the effectiveness of this method. Bearing failure modes are categorized as progressive degradation and sudden failure. Since the reasons for sudden failures in bearings are unknown and difficult to accurately predict, the primary objective of this study is to predict the remaining useful life (RUL) of bearings experiencing progressive degradation.

In the life prediction model, this paper utilizes a time convolutional network (TCN) with incorporated multi-head attention mechanisms. TCNs have a deep structure that can capture long-term dependencies in time series data. Adding multi-head attention mechanisms enables the network to better focus on correlations between different time steps, facilitating the more effective learning of long-term dependencies. The multi-head attention mechanism enables the network to learn different time-step correlations simultaneously, thereby combining features from various time steps to create more comprehensive and meaningful representations. TCNs with multi-head attention mechanisms achieve higher prediction accuracy and effectively capture information about bearing degradation. According to the prediction results, the coefficient of determination (R2) for all cases is above 0.98, which enhances the accuracy of bearing RUL.

Future research directions will include the development of more effective health indices and the implementation of bearing life prediction under various operating conditions. These improvements will hold more practical significance.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Shalabh Gupta, University of Connecticut, United States

Reviewed by: Marco Civera, Polytechnic University of Turin, Italy

Chandrabhanu Malla, Radhakrishna Institute of Technology and Engineering, India