Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Next Article in Journal
Introducing Robotized Stator Cable Winding to Rotating Electric Machines
Previous Article in Journal
Specific Features of Operation of Distributed Generation Facilities Based on Gas Reciprocating Units in Internal Power Systems of Industrial Entities
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Adaptive Band Extraction Based on Low Rank Approximated Nonnegative Tucker Decomposition for Anti-Friction Bearing Faults Diagnosis Using Measured Vibration Data

1
Dynamics Laboratory, School of Mechanical, Aerospace and Civil Engineering, The University of Manchester, Manchester M13 9PL, UK
2
School of Electrical and Electronic Engineering, The University of Manchester, Manchester M13 9PL, UK
*
Author to whom correspondence should be addressed.
Machines 2022, 10(8), 694; https://doi.org/10.3390/machines10080694
Submission received: 13 July 2022 / Revised: 3 August 2022 / Accepted: 12 August 2022 / Published: 15 August 2022
(This article belongs to the Section Machines Testing and Maintenance)

Abstract

:
Condition monitoring and fault diagnosis are topics of growing interest for improving the reliability of modern industrial systems. As critical structural components, anti-friction bearings often operate under harsh conditions and are contributing factors of system failures. Efforts have been cast on bearing diagnostics under the sensor fusion and machine learning framework, whilst challenges remain open on the identification of incipient faults. In this paper, exploiting multi-way representations and decompositions of measured vibration data, a novel band separation method based on the factorization of spectrogram tensors using the low rank approximated nonnegative Tucker decomposition (LRANTD) is proposed and applied to identify detailed fault signatures from the spectral, temporal, and spatial dimensions, flexible for extracting multi-sensor features and multi-dimensional correlations. With the proposed method, informative frequency bands of the latent vibrational components can be automatically extracted, in accordance with the inherent temporal patterns that can be conveniently fed for spectral analysis and fault discrimination. Furthermore, an improved cross-spectrum can be calculated from multi-channel vibrations via LRANTD with enhanced fault features. Based on the real-world vibration data of the accelerated bearing life tests, detailed experimental studies and thorough comparisons to the conventional benchmarks have verified the effectiveness of the reported diagnostic methodology. The proposed method significantly improves the presence of the bearing frequency peaks distinctly over the background noises in the spectrum and hence improves the bearing defect detection process.

1. Introduction

In recent years, condition-based maintenance, involving minimally invasive condition measurements, massive data analytics, and decision making for complex engineering processes, has been considered a critical part of modern intelligent manufacturing for system optimization purposes [1], such as sustainable production and quality improvement. To prevent fatal system failures and allow scheduled maintenance, advanced prognostics and diagnostics are essential for the timely detection of anomalies.
Anti-friction bearings (AFBs) are critical structural components in various mechanical systems such as hydraulic pumps, induction motors, wind turbines, etc., but they are also sensitive parts that account for 45–55% of motor failures [2]. Without maintenance, defects in the AFBs could develop into critical disasters with considerable economic losses and physical injuries. In this context, condition monitoring and fault diagnosis for AFBs can be highly cost-efficient and of great practical significance.
Current approaches of bearing diagnosis can be divided into signal processing methods and machine learning (ML) methods. On the one hand, signal processing methods rely on a series of analytic and synthetic manipulations for processing signals to reveal specific fault characteristics from the time domain, frequency domain, and time-frequency (TF) domain. In this context, typical mainstream techniques include noise cancelation (e.g., time-synchronous average and minimum entropy deconvolution [3]), spectral analysis (e.g., envelope spectrum, spectral kurtosis [4], and cyclostationarity analysis [5,6]), time-frequency analysis (e.g., short-time Fourier transform [7]), sparsity diagnostics (e.g., sparse representations [8] and compressive sensing [9]), and matrix decomposition models for source separation (e.g., singular value decomposition (SVD) [10]). Among them, envelope analysis is one of the most successful approaches and plays an important role in bearing signal enhancement. On the other hand, machine learning methods attempt to identify the latent fault features from a data-driven perspective and learn the diagnostic model from the condition monitoring data in a reverse manner for state inference, such as support vector machines [11,12] and artificial neural networks [13,14]. Despite the successful applications in many data analysis challenges, ML methods are likely confronted with some last mile issues in industrial diagnostics, such as poor generalization to new machine sets and the data imbalance biased towards normal conditions. For both method categories, due to the complicated components and interferences inside the measured signals, challenges remain for robust solutions to bearing fault diagnosis.
With the advancement of sensing technology, more sensors are deployed for monitoring asset conditions and the resulting data exhibits in higher dimensions in the sensor networks. Based on the counteraction assumption of the random noise, sensor fusion and dimension reduction schemes are developed to enhance the inherent fault signatures from multiple measurements [15]. For example, Wodecki et al. [16] proposed a novel adaptive band selection method based on non-negative matrix factorization (NMF) of TF spectrograms and then further applied NMF to multi-sensor fusion for compound fault separation [17]. The diagnostic scope of NMF from the aspect of sparse regularization was studied for bi-spectrum decomposition by Liang et al. [18,19]. He et al. [20] introduced manifold learning for representing fault features in the TF domain by nonlinear mapping.
However, owing to the dimensional limitation of matrix representations, vectorization must be performed to further process higher-order (i.e., order larger than three, or multi-way) data, inevitably neglecting the inherent localized features such as multi-dimensional correlations. Hence, based on multilinear algebra, researchers resort to tensors for representing multi-way arrays and develop tensor decomposition (TD) for multidimensional data mining. As an extension to matrix component analysis, TD fills the gap from two-way analytics to multi-way component analytics. Examples of well-known TD models include the canonical polyadic decomposition or parallel factor analysis (CANDECOMP/PARAFAC decomposition or CPD in short) [21], the Tucker decomposition (aka the higher-order singular value decomposition (HOSVD) with orthogonal constraints) [22], nonnegative CPD (NCPD), nonnegative Tucker decomposition (NTD), etc. In [23], Cichocki and Mandic et al. gave a thorough review of TD and its applications for multi-way component analysis. A detailed guide of tensor methods for multi-sensor signal processing was summarized by Miron et al. [24]. Whereas in condition-based maintenance, limited applications of TD were found within current literature. Zhao et al. represented multiple intrinsic mode functions (IMF) into a Hankel tensor and applied tensor rank decomposition for bearing fault extraction [25]. The improved support tensor machine was proposed to classify the high-order condition monitoring data for rotating machinery diagnosis with imbalanced data [26]. Multi-dimensional noise reduction was introduced by Hu et al. to signal enhancement based on HOSVD for rotating machine diagnosis [27]. Sparse NCPD was applied to the impact feature extraction from the three-way TF representation of a single vibration signal [28] and to the fault identification using multi-channel measurements [29]. These examples have preliminarily demonstrated the efficacy of the high-order tensor representations and the corresponding data mining techniques for fault diagnosis purposes. Following this line of research, some issues that limit the diagnostic performance of current TD models are notable and deserve necessary improvements:
(1)
For CPD, the rank-one decomposition property requires the same number of lower ranks for each tensor mode (the dimensions of a tensor are often referred to as modes), which might be incompatible to express the inherent correlations where each mode has various number of latent components;
(2)
Although different low-rank dimensions are allowed, the Tucker decomposition is often criticized for its poor uniqueness and the curse of dimensionality (i.e., the core tensor is dense whose entry number scales exponentially with the tensor order) [30];
(3)
Diverse ways of representing vibration data as tensors can be found while a physically interpretable way is needed. To this end, the multi-dimensional correlation across different modes (such as the spectral, temporal, and spatial modes) might need more intuitive illustrations for diagnostic purposes.
To address the above issues, a novel diagnostic methodology has been formulated on the basis of time-frequency representations and the low rank approximated nonnegative Tucker decomposition [30,31], or LRANTD, a two-step model combining low-rank approximation (LRA) and NTD, which adaptively extracts the informative bands of the latent temporal patterns from the vibration signals. In the proposed method, a three-way TF tensor is generated from the single-channel or multi-channel vibration measurements by the introduced segmentation operation; it can be efficiently factorized by LRANTD into three factor matrices (embedding the latent features) and a core tensor (preserving the multi-dimensional correlations). The interpretation from the decomposed factors to the physical dimensions (including the spectral, temporal, and spatial dimensions) is elaborated in detail and the multi-dimension correlations preserved by the core tensor are illustrated graphically via the Hinton diagrams. Furthermore, the LRANTD model can be solved efficiently based on the iterative algorithms that integrate hierarchical alternating least squares (HALS) and multiplicative update (MU), which allow the flexible choice of the lower rank numbers for each tensor mode as compared to the CPD model and reduces the computational cost against the standard NTD.
For performance verification, the proposed method was applied to fault diagnosis using the bearing vibration data measured from multiple accelerated run-to-failure experiments [32]. Three typical health states including the incipient fault state, the maximal kurtosis state, and the severe fault state (later marked by Case A, B, and C) were analysed, whose results indicated the effectiveness of the proposed method for informative band extraction and fault signature enhancement. For the analysis using the single-channel signal, the diagnostic performance of LRANTD was verified through the comparative studies against a few successful benchmarks including smoothed envelope analysis, fast Kurtogram (FK), and empirical mode decomposition (EMD). For multi-channel sensor fusion, the cross spectrum computed based on LRANTD also outweighed the conventional approach based on the cross-power spectral density (CPSD).
This paper first provides the details of experimental vibration data and the observations made on their trend of the measured vibration RMS and Kurtosis values in Section 2. The principle of the LRANTD model is presented in Section 3, where the diagnostic methodology is also proposed for the defect detection in the anti-frication bearings. In Section 4, the LRANTD model is applied to identify fault signatures from the single-channel acceleration signal. Comparative studies are subsequently provided in Section 5 to show the effectiveness of the proposed method. Then, the cross-spectral analysis in Section 6 demonstrates the usage of the proposed LRANTD method for multi-sensors fusion to further reduce the background noises in the spectrum and hence improve the bearing fault detection process. Finally, the concluding remarks are given in Section 7. Basic notations are listed in Table 1.

2. Experimental Rig and Measured Vibration

The bearing vibration data collected from the accelerated life tests available in the XJTU-SY bearing database [32] are beneficial for the diagnostic performance validation under deteriorating health conditions and thus used in this study. A simple schematic of the experimental bearing test rig is shown in Figure 1. It mainly consists of an alternating current (AC) induction motor, an intermediate shaft supported through two ball (LDK UER204) bearings, and a hydraulic loading system to apply the load on the bearing in the horizontal direction. The vibration database [32] are composed of a total of 15 experimental groups based on the three different sets of operating conditions of the running speeds and load levels. Further details on the rig, conducted experiments, and the vibration data can be found in the study by Wang et al. [32].
In this work, the measured vibration acceleration data for the ‘Bearing 1_1’ [32] are considered and analysed to demonstrate the usefulness of the proposed method. The vibration acceleration data for the ‘Bearing 1_1’ were collected by Wang et al. [32] with the shaft speed of 35 Hz and a hydraulic load of 12 kN being applied to the bearing. The ‘Bearing 1_1’ data set contains 123 data records, each of which is acquired at an interval of 1 min and encompasses 1.28 s of vibration acceleration measured on the bearing in the vertical and horizontal directions. These 1.28 s data were collected at the sampling frequency, f s of 25,600 Hz. The 123 data records means that the run-to-failure experiment lasted for two hours and 3 min until the stopping criteria was met (i.e., the maximum acceleration amplitude exceeding 10 times that in the normal state). The bearing specifications and the calculated faults characteristic frequencies are listed in Table 2.
To better observe the degradation process of the test bearing, the trends of overall bearing vibration in both directions are presented in Figure 2 based on two different indicators, i.e., Root Mean Square (RMS) and Kurtosis, for each acceleration record. From the normal state to severe-fault state, the RMS indicator follows an increasing trend and surges since the 78th minute, whereas Kurtosis is more fluctuant and rises much earlier from around the 66th minute, then peaking at the 78th. Based on this observation, as an impulsiveness indicator, Kurtosis is justified for recognizing the degradation of AFBs. By comparisons, both indicators have higher amplitudes in the horizontal direction than in the vertical direction, which is due to the hydraulic loading on the horizontal axis.
Based on the vibration trend, three typical timestamps indicating incipient-fault state, the maximum Kurtosis state, and the severe-fault state were selected for further analysis, i.e., the 66th, 78th, and 80th records, as denoted by the capital letters ‘A’, ‘B’, and ‘C’ in Figure 2.

3. Proposed Method

The concept of the LRANTD method is presented first and then the diagnostic approach using this concept is proposed.

3.1. Nonnegative Tucker Decomposition Based on Low-Rank Approximation

This section is divided into a number of steps to achieve the optimised LRANTD model to aid the understanding.

3.1.1. Standard Nonnegative Tucker Decomposition (NTD)

The tensor decomposition theory was first developed for factor analysis and data mining of multi-way psychometric data [22]. As an extension of matrix component analysis, the Tucker decomposition, or higher-order SVD of an N -way real-value tensor 𝓨 I 1 × I 2 × I N , can be defined as the mode- n products of the core tensor and the factor matrices [22],
𝓨 = 𝓖 × 1 U ( 1 ) × 2 U ( 2 ) × N U ( N )
where 𝓖 R 1 × R 2 × R N is the core tensor and U ( n ) I n × R n ,   n = 1 ,   2 , , N , represent N mode- n matrices, respectively. The multilinear rank of the Tucker decomposition can be denoted by the N -tuple ( R 1 , R 2 , , R N ) . Note that R n I n , a lot of multi-dimensional data, e.g., RBG images and multi-sensor signals, could be well represented by low-rank representations. To that end, the HOSVD method is the workhorse for computing such representations for unconstrained Tucker decomposition. Similar to the column-wise orthogonality of the singular matrices in SVD, the mode- n matrices of the HOSVD are also semi-orthogonal, i.e., U ( n ) T U ( n ) = I . Due to the rotational indeterminacy, the unconstrained Tucker decomposition can be restricted by the lack of uniqueness and the curse of dimensionality [30].
For the better uniqueness, Nonnegative Tucker Decomposition (NTD) is put forward by relaxing the orthogonality requirement in Tucker decomposition to approximate the original tensor via a set of nonnegative factors, i.e., 𝓖 and   U ( n ) are subject to elementwise nonnegativity constraint [33],
𝓨 ^ = 𝓖 × 1 U ( 1 ) × 2 U ( 2 ) × N U ( N ) ,     s . t .   𝓖 ,     U ( n ) 0
or equivalently in a matricization form,
Y ^ ( n ) = U ( n ) ( U ( N ) U ( N 1 ) U ( n + 1 ) U ( n 1 ) U ( 1 ) ) G ( n )   T = U ( n ) V ( n )   T
where Y ^ ( n ) I n × i n I i and G ( n ) R n × i n R i are the mode- n matricizations of tensor 𝓨 ^ and 𝓖 , and V ( n ) is an auxiliary notation defined by,
V ( n ) = ( U ( N ) U ( n + 1 ) U ( n 1 ) U ( 1 ) ) G ( n )   T
The above nonnegative approximation is achieved by minimizing the approximation error based on the squared Frobenius norm, as in Formula (5) that further unfolds into matrix form, which can be solved effectively by a few alternating algorithms in the coordinate descent family.
min   𝓖 ,     U ( n ) 0 1 2 𝓨 𝓨 ^ F 2 = 1 2 Y ( n ) U ( n ) V ( n )   T F 2 ,   n { 1 , , N }
Two main attributes are endowed by the imposed nonnegativity: (1) only additive operations are allowed without subtractions, thus enhancing the interpretability in feature engineering applications; (2) since many entries of the core tensor becomes zeros, the resulting factors are often sparse, exhibiting parts-based features and alleviating the curse of dimensionality.

3.1.2. The Low Rank Approximated NTD (LRANTD) Model

Many efforts have been cast to reduce the computational cost in NTD for the low-rank approximation (LRA) of large-scale data, especially the computation of the gradients for the tensor-based cost function. Given its flexibility to other well-established LRA techniques and noise robustness, the low rank approximated nonnegative Tucker decomposition, or LRANTD, developed by Zhou, Cichocki et al. [30,31] is introduced as the key tensor factorization model in this work, which consists of two main steps:
(1)
LRA. Perform unconstrained Tucker decomposition as in (1) using HOSVD to obtain an approximation of the original tensor, 𝓨 , which gives
𝓨 𝓨 ˜ = 𝓖 ˜ × 1 U ˜ ( 1 ) × 2 U ˜ ( 2 ) × N U ˜ ( N )
and the mode- n unfolding of the LRA tensor, 𝓨 ˜ ,
Y ˜ ( n ) = U ˜ ( n ) G ˜ ( n ) ( U ˜ ( N ) U ˜ ( n + 1 ) U ˜ ( n 1 ) U ˜ ( 1 ) ) T = U ˜ ( n ) V ˜ ( n )   T
(2)
NTD. Minimize the new cost functions with the nonnegativity constraints on the factors as in (8). Note that given the tensor order N , a total of N unfolding forms from 𝓨 ˜ are defined in (7) and the sequence of N minimization problems should be calculated for each factor matrix, U ( n ) .
min   𝓖 ,     U ( n ) 0 D N T D = 1 2 𝓨 ˜ 𝓨 ^ F 2 = 1 2 U ˜ ( n ) V ˜ ( n )   T U ( n ) V ( n )   T F 2 ,   n { 1 , , N } .
Without frequent visits to the original big tensor, the LRA step not only significantly reduces the computation complexity for the subsequent optimizations, but also suppresses the interference in the raw data tensor [30]. Rather than directly implementing NTD of a tensor, LRANTD combines HOSVD with NTD to provide a trade-off between efficiency and approximation error.

3.1.3. Optimization Algorithms

Both (5) and (8) are known as nonnegative least squares problems, which can be solved by a few well-known first-order methods from the area of NMF, such as multiplicative update (MU) [34], accelerated proximal gradient (APG) [35], the hierarchical alternating least squares (HALS) algorithm [36], active set methods [37], etc. By a crafty choice of learning step, these methods avoid the searching process for an appropriate step size in gradient descent and can be further extended to various tensor decomposition models.
For ease of comprehension, the following optimization algorithm is provided for three-way tensor decomposition, i.e., N = 3 . The lower rank dimensions of the decomposition are given as ( R 1 , R 2 , R 3 ) ideally in accordance with the number of latent components. And the desired variables U ( n ) ,   n { 1 ,   2 ,   3 } and 𝓖 are initialized randomly within the interval of 0 and 1. See Figure 3 for an intuitive illustration of the three-way tensor decomposition and consult [30] for further details of extensions to higher-order tensor ( N > 3 ).
  • The iterative update rule for mode- n (factor) matrices U ( n )
The LRANTD of the three-way tensor is computed based on the matricization operation undertaken for each tensor mode. For example, Y ( 1 ) I 1 × I 2 I 3 is the mode- 1 matricization form of the tensor, 𝓨 I 1 × I 2 × I 3 . Taking the LRA step in (7), the approximation term Y ˜ ( 1 ) can be obtained and the mode- 1 factor matrix U ( 1 ) I 1 × R 1 can be calculated via the NTD step in (8).
To generalize, for the mode- n matrix, U ( n ) I n × R n ,   n { 1 ,   2 ,   3 } , the low-rank approximation Y ˜ ( n ) I n × i n I i is first attained in the LRA step, and the NTD step applies the HALS algorithm to update only one column of U ( n ) at a time subsequently. By defining the residue term as
Y ˜ ( n )   r = Y ˜ ( n ) i r u i ( n ) v i ( n ) T ,   r { 1 , , R n } ,   n { 1 ,   2 ,   3 } ,
Equation (8) expands to a series of minimization problems that sequentially optimize u r ( n ) while fixing the remaining entries and variables, as given by
min   u r ( n ) 0 D N T D = 1 2 Y ˜ ( n )   r u r ( n ) v r ( n ) T F 2 ,   r { 1 , , R n } ,
where u r ( n ) I n × 1 and v r ( n ) ( i n I i ) × 1 are the r -th column vectors of U ( n ) and V ( n ) with r { 1 , 2 , , R n } . Note that v r ( n ) always serves as an auxiliary term and is decoupled from U ( n ) as in (4).
Then, the gradient of D N T D with respect to u r ( n ) can computed as
D N T D u r ( n ) = u r ( n ) v r ( n ) T v r ( n ) Y ˜ ( n )   r v r ( n )
Assigning zero to the gradient and imposing the nonnegativity constraints further gives
u r ( n ) = max { 0 , Y ˜ ( n )   r v r ( n ) ( v r ( n ) T v r ( n ) ) 1 }
Substituting the residue in (7) further gives the following update rule,
u r ( n ) u r ( n ) + 1 v r ( n ) T v r ( n ) max { 0 , Y ˜ ( n ) v r ( n ) U ( n ) V ( n )   T v r ( n ) }
Thus, by updating a total of R n columns, the mode- n matrix U ( n ) can be obtained. Consequently, the three factor matrices, U ( 1 ) , U ( 2 ) , and U ( 3 ) can be successively computed, as in Figure 3.
  • The iterative update rule for core tensor 𝓖
Next, the MU algorithm can be applied in the NTD step to update the core tensor, 𝓖 R 1 × R 2 × R 3 . From (6) and (8), the gradient of the cost function with respect to 𝓖 can be computed as in
D N T D 𝓖 = 𝓨 ^ × n I U ( n )   T 𝓨 ˜ × n I U ( n )   T
where I reduces to I = { 1 ,   2 ,   3 } as N = 3 . By taking the step size in the gradient descent as η = 𝓖 𝓨 ^ × n I U ( n )   T , the additive update rule reduces to pure multiplicative operations, providing an efficient computation algorithm for the core tensor,
𝓖 𝓖 max { 0 ,     𝓨 ˜ × n I U ( n )   T }   max   { 0 ,     𝓨 ^ × n I U ( n )   T } = 𝓖 max { 0 ,     𝓖 ˜ × n I U ˜ ( n ) × n I U ( n )   T }   max   { 0 ,     𝓖 × n I U ( n ) × n I U ( n )   T }
where denotes the element-wise multiplication.

3.2. Diagnostic Scope of the LRANTD Model

In this section, the diagnostic application of the LRANTD model is illustrated by analysing multi-way vibrational data. The fundamentals of the tensor-based diagnostics are explained with practical considerations. Then, the methodology flowchart for fault diagnosis of AFBs is presented.

3.2.1. Data Pre-Processing: Moving Segmentation

Given the input vibration data, either from single-channel or multi-channel sensor measurements, the pre-processing operation termed moving segmentation is introduced to the raw signals to generate multi-way arrays that allow the subsequent multilinear data mining. Given a signal from i -th channel measurement, y i L , I { 1 , , M } , where M denotes the total number of sensor channels and L represents the signal length in data points, a total of S segments can be resampled from y i by a moving-window fashion. By segmentations, the j -th segment, y i ,   j L s , of the source signal, y i , is defined by
y i ,     j = [ y i ( j ) ,   y i ( j + 1 ) ,   , y i ( L s + j 1 ) ] ,     j = 1 ,   2 ,   , S
where L s   ( L s L S + 1 ) controls the length of each segment and is suggested to be an integral multiple of the shaft-rotation period. The step size of the segmentation is set at 1 such that the adjacent segments are highly overlapped and self-similar. The above settings are important to enhance the inherent components in the signals as the random interference could be suppressed by the subsequent decomposition in ideal conditions.

3.2.2. Tensorization: Generating the Time-Frequency Tensor

The vibration signals from a localized fault in the AFBs are known to be non-stationary, i.e., the statistical behaviours of the signals are inconstant, and are indeed pseudo-cyclostationary as pointed out by Antoni and Randall [38]. As effective non-stationarity analysis tools, time-frequency (TF) representations provide a thorough depiction of signal variations by balancing the time-domain and frequency-domain features. On each one-way signal segment, y i ,   j , the short-time Fourier transform (STFT) is first performed, as given by
STFT ( i 1 , i 2 ) = τ = 0 L w 1 y i , j ( i 2 + τ ) w ( τ ) e j 2 π τ i 1 / I 1
where i 1 = 0 ,   1 ,   , I 1 1 denotes the frequency bin with I 1 bins in total, i 2 = 0 ,   1 , , I 2 1 is the time bin for I 2 bins in total, and w ( τ ) is the window function of length L w . In practice, STFT can be efficiently computed by implementing fast Fourier transform for each windowed frame. Then, by taking the absolute value of the STFT result, the corresponding two-way TF spectrogram, Y i , j , is attained with each entry defined as
Y i , j ( i 1 , i 2 ) = | STFT ( i 1 , i 2 ) |
Consequently, by arranging the adjacent TF spectrograms, one can construct a three-way TF tensor from single or multi-sensor measurements by absorbing a total of I 3 = M S spectrograms into the last mode. For illustration, Figure 3 provides an example of the TF tensor, 𝓨 I 1 × I 2 × I 3 , spanned by three physical modes, i.e., the frequency mode, the time mode, and the data-segment mode.
Extensive research utilized multiple sensors for vibration measurements under the assumption that comprehensive information, including the fault features, could be obtained, while it might be more often a case where a single measurement is analysed without detailed clarification. The motivations of tensorization the raw vibrational signals are two folds:
  • Tensors are natural representations for multi-way data beyond the dimensional restriction of matrix analysis, which makes possible the data mining for higher-order connections among multi-channel vibration measurements;
  • The time-frequency signatures of localized faults in AFBs are likely to be low-rank, as inherited by the repetitive occurrences of the impact-like structures in the TF spectrogram. It is therefore straight-forward to represent the signature-of-interest through low-rank learning models.

3.2.3. Tensor Decomposition by LRANTD

After selecting the low-rank dimensions as ( R 1 , R 2 , R 3 ) for decomposition, one can further apply the presented LRANTD algorithm to the TF tensor for multilinear data mining and dimensionality reduction, which returns four factors with interpretable physical meaning (recall Figure 3 for illustrations):
(1)
Based on the proposed method, the latent frequency bands of the TF spectrograms can be automatically extracted without subjective selection. The mode- 1 matrix U ( 1 ) I 1 × R 1 is called the spectral factor that preserves the R 1 dominant frequency bands or spectral subfactors (SS) in its columns, u r ( 1 ) I 1 . Each SS is normalized within the interval [ 0 ,   1 ] to indicate the informative bands, which in turn rescales the corresponding temporal pattern to have meaningful amplitude;
(2)
The mode- 2 matrix U ( 2 ) I 2 × R 2 serves as the temporal factor that adaptively identifies the R 2 principal encoding patterns or the temporal subfactors (TS) embedded in the corresponding frequency bands. Essentially, these temporal components are the envelope profiles that reveal the variation in TF energy within the informative bands as extracted in the corresponding SS;
(3)
Similarly, the mode- 3 matrix U ( 3 ) I 3 × R 3 clusters I 3 data segments into R 3 segment-wise subfactors (SWS). Each SWS gives a hint about its relevance to the segments and sensor channels (or spatial locations). For multi-sensor measurements, the multiple SWS represent the dominating spatial components;
(4)
Therefore, the entries of the core tensor 𝓖 R 1 × R 2 × R 3 reflects the correlations across the spectral, the temporal, and the spatial dimensions, which can be visualized intuitively. In other words, the core tensor provides direct information about the spectral locations of the linked temporal patterns.
Note that the multilinear rank is always smaller than the input dimension ( R n I n ) , the presented LRANTD also provides an efficient multi-sensor fusion scheme that greatly reduces the computational storage, from O ( n I n ) to O ( n I n R n + n R n ) , which is of practical value in modern manufacturing systems with sensor networks and cloud-based computing for condition monitoring.

3.2.4. Evaluation Based on Spectrum Inspection

Finally, based on the separated TS, Kurtosis values are used to assess their impulsiveness and the frequency spectra are calculated to identify the characteristic frequencies of bearing faults. Overall, the proposed the diagnostic methodology based on LRANTD can be summarized into the flowchart in Figure 4.

4. Data Analysis and Case Studies

In this section, the proposed method was first applied to single-channel vibration (i.e., M = 1 ) for AFB diagnosis under three typical health conditions. The diagnostic results based on LRANTD were elaborated. Further comparative studies against a few benchmark models were provided in Section 5.

4.1. Case A

It can be seen from the vibration trend in Figure 2 that the Kurtosis at timestamp A just begins to climb during the incipient stage of bearing degradation. Provided that the hydraulic load was applied from the horizontal axis, the horizontal acceleration was selected for further analysis. In Figure 5a, a few irregular impacts appear in the acceleration signals, where judgment can hardly be made about the health state of the test bearing. In the amplitude spectrum in Figure 5b, the spectral peak points at around 1.2 kHz with multiple bands rich in energy. Then, the signal was high-pass filtered with the cut-off frequency at 1 kHz. After performing STFT with the 63-point Hanning window (using consistent settings hereafter), the TF spectrogram of the filtered signal is still complex as shown in Figure 6a, where some impacts signatures are slightly clearer in the bands centred at 1.2 kHz, 5.5 kHz, 6.5 kHz, and 11 kHz. Hence, the early fault signal is hardly recognizable due to considerable background contaminations.
Setting the segmentation length to 20 times of rotation periods, i.e., 20 f s / f r 14640 (0.57 s), the frequency bins of STFT as 256, and the number of segments as S = 8   ( I 3 = 1 S = 8 ) , the TF tensor in Case A, 𝓨 256 × 14640 × 8 , could be constructed according to the procedure in Section 3, as shown in Figure 6b.
As the number of latent features cannot be priorly known and the source determination remains a major challenge in current research, the low-rank dimensions ( R 1 ,   R 2 ,   R 3 ) were selected based on trial and error as ( 4 ,   4 ,   2 ) . Then, LRANTD was performed on 𝓨 for multi-way data mining based on the algorithms combining the HALS and MU update rules, which returned a factor set including the spectral factor U ( 1 ) 256 × 4 , the temporal factor U ( 2 ) 14640 × 4 , the segment-wise factor U ( 3 ) 8 × 2 , and the core tensor 𝓖 4 × 4 × 2 , as in Figure 7.
(1)
Note that each SS in U ( 1 ) had been normalized to intuitively reveal the spectral weight of the correlated TS, U ( 1 ) thereby adaptively extracted multiple dominant bands from the TF tensor, as shown in Figure 7a;
(2)
For the temporal factor U ( 2 ) , it was found that the TS 3 and TS 4 exhibit apparent impulsive features and have relatively higher kurtosis, as shown in Figure 7b:
(3)
Each SWS in the segment-wise factor U ( 3 ) was found to have a correlation of over 97% to all data segments. This is because all data segments are sampled from the same channel measurement and hence resemble one another;
(4)
Unlike the NMF factors, the mode matrices of LRANTD are implicitly linked with one another due to the flexible choice of lower rank dimensions. Such interconnections are embedded in entries of the core tensor, 𝓖 . For visualization, the Hinton diagram of each frontal slice of the core tensor is shown in Figure 7d,e. For example, from Figure 7d, TS 4 was found to have a more distinctive correlation with SS 2 against SS 1, indicating TS 4 is more likely to contain high-frequency components as inferred from Figure 7a.
Then, the amplitude spectrum of each TS was computed and that with respect to the 4th TS is shown in Figure 8. Both BPFO and its 2nd harmonics are clear and dominant, as marked by the green dash-dot line. Tracing the spectral property of this temporal pattern (TS 4) according to the Hinton diagram in Figure 7, it was found that the dominant frequency band was extracted by the spectral factor in SS 2, indicating that the band from 11 to 12 kHz is likely to comprise the outer-race fault components.
Therefore, the above result has shown the efficacy of the proposed diagnostic methodology based on LRANTD for informative bands extraction, which is further verified by comparative studies in Section 5.

4.2. Case B

In Case B, the bearing vibration data with the maximum kurtosis (at the 78th minute) were analysed. From the horizontal acceleration signal in Figure 9a, it can be recognized that more impacts appear in repetitive but uneven patterns as compared to the Case A signal (Figure 5a). In the corresponding frequency spectrum, the dominant bands remain similar to the early stage as in Figure 5b, whereas the bands centred at 7.5 kHz and 12 kHz have gained more intense energy. Correspondingly, the TF spectrogram based on STFT is given in Figure 9c. The impact-like signatures with high amplitudes lie above 10 kHz in the TF domain. While the 1k Hz components, dominating the frequency spectrum, show more regular patterns but much lower energy in the spectrogram. Given the maximum kurtosis of the analysed signal, no pivotal diagnostic information can be learned from the above observations directly. By TF analysis, one can better recognize the spectral locations of the inherent temporal behaviours, while the complex TF features still limit its diagnostic performance, for which more advanced analysis should be introduced.
Presumably, the fault signatures in Case B (maximum kurtosis state) are more prominent than that in Case A (at the early stage). Therefore, the low-rank dimensions were reduced to attain more concentrated fault features, i.e., ( R 1 ,   R 2 ,   R 3 ) was set to ( 3 ,   3 ,   2 ) . Keeping other settings same to the Case A, LRANTD was performed for factorizing the TF tensor, 𝓨 256 × 14640 × 8 , obtained in the similar fashion according to Section 3.
The resulting factors of the LRANTD are shown in Figure 10. Three SS with different dominant frequency bands were extracted in the spectral factor, i.e., [ 11 ,   12.5 ] , [ 1 , 6 ] , and [ 7 , 8 ] kHz. The Kurtosis values of the three TS in Figure 10b are 19.5, 3.54, and 8.36, respectively. Furthermore, the segment-wise factor in Figure 10c shows an opposite trend between the two SWS, indicating apparent differences between the separated components.
Next, the amplitude spectrum of each TS was computed and those of the first two temporal factors were shown in Figure 11. As the test bearing had further degraded in Case B, both the amplitude spectra of the TS 1 and TS 2 have captured the characteristics of the outer-race fault, whereas the spectrum of TS 2 shows more prominent signs with higher peaks at the BPFO and its 2nd multiple. Recalling the core tensor in Figure 10d,e, the TS 1 are dominated by SS 1 and TS 2 by SS 2, indicating that the outer-race fault component is more likely to locate itself in the bands of [ 1 , 6 ] and [ 11 ,   12.5 ] kHz.

4.3. Case C

In Case C, the source separation capability of the diagnostic algorithm based on LRANTD was demonstrated using a severe fault case (at the 80th minute) as the outer-race fault developed. For consistency, the horizontal channel acceleration signal was analysed. As shown in Figure 12, the time-domain acceleration signal exhibits impulsive behaviour of high-amplitude (over 10 g) in a cyclic (or regular) manner, while the frequency spectrum shows that the most dominant band lies in [ 1 ,   2 ] kHz. Moreover, a closer look at the STFT spectrogram in Figure 12c provides that the dominant band centred at around 1.2 kHz contains a series of impact signatures, while other bands remain highly complex and may scarcely help diagnosis.
Next, the same settings as in the flowchart in Figure 4 were used to obtain the three-way TF tensor 𝓨 256 × 14640 × 8 and to perform LRANTD for tensor factorization. The resulting four factors are shown in Figure 13. The spectral factor in Figure 13a has extracted three SS whose dominant frequency bands are centred at 1.2 kHz, 11.5 kHz, and 7.5 kHz, respectively. Three TS in Figure 13b show obvious impulsive patterns whose kurtosis values from TS 1 to TS 3 are 2.56, 8.06, and 4.30. Note that the core tensor in Figure 13d,e have large values in the diagonal entries, indicating the one-to-one correspondence between each SS and TS (e.g., SS 1 is most relevant to TS 1, so is SS 2 to TS 2, etc.).
Further, the amplitude spectrum of each TS was calculated and those of the 1st and 3rd TS (or TS1 and TS3) were given in Figure 14. Both spectra contain apparent frequency component at BPFO and its high-order harmonics. Interestingly, it is notable that the impact sequences in TS1 and TS3 are indeed quite different, as they possess different phase properties. In Figure 13b, this difference is highlighted by marking the local peaks of TS1 by red dash lines and those of TS 3 by green dots. The observations suggest that although TS1 and TS3 are similarly modulated by BPFO, they are intrinsically two different sources of impact components with different spectral excitations. Such spectral differences are revealed by the extracted spectral factor in Figure 13a, which shows that the fault impacts of TS1 was excited at about 1.2 kHz while the fault impacts of TS 3 are dominated by a broad band that peaks at 7.5 kHz.
On the above basis, the source separation capability of LRANTD is demonstrated. As the test bearing had undergone a severe outer-race fault, it is justifiable to infer that these two different impact patterns were resulted by the rolling elements entering and exiting the defect zone.

5. Comparisons with Benchmark Methods

In this section, to verify the performance of the proposed diagnostic methodology, the result of the LRANTD in the early fault case (Case A) was compared against a few practical benchmarks including envelope analysis, FK, and EMD. For illustration, the spectral result based on the LRANTD in Case A is recalled and shown in Figure 15.

5.1. Envelope Analysis

Figure 16 shows the smoothed envelope spectrum of the filtered (high pass above 1 kHz) signal of the horizontal acceleration in Case A. For noise reduction, the averaging operation was performed based on 90% overlap between the adjacent Fourier frames. Apparently, the envelope spectrum peaks at the shaft rotation frequency, as marked by the red dot, whereas the BPFO component is less prominent due to the masking of ambient interferences.
Compared to the envelope analysis, the spectral result of LRANTD in Figure 15 effectively removes the cyclic interference from the shaft rotation and presents clear fault characteristics for diagnosis.

5.2. Fast Kurtogram

As an effective band selector, fast Kurtogram is another pragmatic benchmark for rotating machine diagnosis [4]. The Kurtogram was computed based on the decimated filter banks and shown in Figure 17a. The result indicates that the impulsive components are more likely to dominate over 10 kHz, which is consistent with the spectral extraction result of LRANTD, as the SS 2 in Figure 7a.
Based on the returned optimal filter band with the highest kurtosis value, i.e., [ 11200 ,   12800 ] Hz, the envelope spectrum of the filtered signal is further given in Figure 17b. Here, the square operation of the envelope has been omitted for comparisons. In this result, the BPFO component is also dominated by the rotation frequency (dotted in red) and is of much lower energy than the envelope spectrum in Figure 16 due to a small band width.
In contrast, the LRANTD method benefits from the source separation capability of tensor decomposition that facilitates adaptive band extraction and hence provides more prominent diagnostic signatures.

5.3. Empirical Mode Decomposition

Lastly, the results were compared for the Case A signal using the empirical mode decomposition and the LRANTD method. Figure 18 shows the envelope spectrum of the 1st intrinsic mode function (IMF) (for comparison, the best result among 10 IMFs was selected), where the similar masking behaviours to the BPFO can be found from the shaft speed (dotted in red) and the other low-frequency components. This can be attributed to the incipient degradation condition where the BPFO component is low in energy and heavily contaminated.
On the other hand, the proposed method effectively enhanced the inherent signatures within the original measurement by fusing multiple signal segments in an overlapping fashion and further mining the multi-way arrays by the LRANTD model.

6. Cross-Spectrum Methods Based on Sensor Fusion

Next, the proposed methodology was applied to multi-sensor fusion based on the two-channel vibration data. Since both horizontal and vertical acceleration signals were measured, we set M = 2 ,   S = 4 , and the total number of data segments I 3 = M S = 8 , to create a three-way TF tensor, 𝓨 256 × 14640 × 8 , in a consistent manner as in Section 3.
In the following, the three typical cases (from Case A to C) were revisited but analysed with two-channel fusion via LRANTD to present a novel cross-spectrum. The cross spectral results from the LRANTD model are shown in Section 6.1 and the comparative results based on the conventional cross power spectral density are given in Section 6.2.

6.1. The Cross-Spectrum Based on LRANTD

By performing the proposed method on the data in Cases A to C, two-directional acceleration signals were first fused and then factorized using the LRANTD model. Consequently, the cross spectra based on the fast Fourier transform of the temporal factors were obtained in succession, as shown in Figure 19. From Case A to Case C, it can be observed that the BPFO component gradually increases in spectral amplitude, and the harmonics of BPFO also become more prominent as the out-race fault developed.

6.2. The Conventional Method

For comparisons, the cross power spectral density (CPSD) of the envelope signals from two orthogonal acceleration measurements was calculated using the Welch’s averaged periodogram method with 90% overlap. The resulting spectra of Case A, B, and C are, respectively, shown in Figure 20 and adjusted to the same scale as those in Figure 19. For Case A and Case B, the BPFO component and its harmonics in the CPSD are interfered by other cyclic components, in particular the shaft rotation frequency, as marked by the red dots in Figure 20a,b.
In contrast, the cross spectra obtained by the LRANTD-based sensor-fusion in Figure 19 provide more prominent fault features in terms of spectral amplitude and concentration, which eventually confirms the effectiveness of the proposed diagnostic method.

7. Summary of Observations

Table 3 and Table 4 summarise the observations made from this study using one sensor and 2-sensors data fusion. The results of the proposed LRANTD method are also compared. The background noise is calculated based on the mean of spectral amplitude within the range 0-400 Hz. The ratio in the bracket is calculated using the ratio of the spectral amplitude at a frequency to the mean background noise. This study clearly indicates the predominance of the outer-race defect in the bearing. The “Remark” in Table 3 and Table 4 is given by comparing the BPFO amplitude ratio (BPFO amplitude by background noise) with respect to the BPFO Baseline ratio of the envelope analysis as the reference (underlined in the table) within each case study. The following advantages are clearly observed for the proposed LRANTD method that are expected to improve the bearing fault diagnosis process significantly:
(i)
Reducing the presence of the shaft speed related frequency;
(ii)
The signal-to-noise ratio is much better for the proposed LRANTD method, hence improving the appearance of the bearing defect frequencies in the spectrum;
(iii)
Significant improvements in the BPFO amplitude, generally in the order of 90% for the sensor fusion cases (B and C) compared to the envelope cross-spectrum.

8. Conclusions

In the presented work, a novel diagnostic methodology based on the LRANTD model was proposed to tensorize the single-channel or multi-channel vibration into multi-way time-frequency tensors and to identify the latent fault signatures with different spectral excitations, temporal patterns, and sensor locations. Detailed observations and diagnostic analysis were provided based on the XJTU-SY bearing run-to-failure data, including the illustrations on three typical cases during the bearing lifespan to demonstrate its superiority in adaptive frequency band extraction and fault feature separation.
The LRANTD model was first introduced to bearing fault diagnostics to overcome the drawbacks of unconstrained tucker decomposition while enhancing the performance of NTD by low-rank approximation. The higher-order correlations among the spectral, temporal, and segment-wise dimensions were intuitively visualized using the Hinton diagrams. Moreover, the cross-spectrum derived from LRANTD also shows improved capability for sensor-fusion. The experimental and comparative studies have confirmed the diagnostic performance of the proposed tensor-based method. The proposed LRANTD method for the bearing fault diagnosis is significantly increasing the signal-to-noise ratio and reducing the presence of the shaft speed frequency and its harmonics. This improves the appearance of the bearing defect frequencies in the spectrum and hence eases the bearing defect detection process. The multi-sensors fusion further helps in improving the signal-to-noise ratio and the amplification in the peak amplitude of the bearing defect frequency by nearly 90% as compared to the envelope analysis. The proposed method is showing encouraging potential for the reliable bearing fault diagnosis, and useful for the industrial application. Further studies are underway with more case studies to observe the robustness of the proposed method for the industrial application. Since the computational time is not significant compared to the conventional approaches, the method can further be explored for the online monitoring.

Author Contributions

Conceptualization, H.W., L.Z. and J.K.S.; methodology, software, writing—original draft preparation, H.W.; resources, supervision, project administration, L.Z. and J.K.S.; writing—review and editing, L.Z. and J.K.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The XJTU-SY bearing dataset is available online at https://biaowang.tech/xjtu-sy-bearing-datasets (accessed on 5 January 2022).

Acknowledgments

The first author would like to thank the China Scholarship Council and the University of Manchester for their joint sponsorship.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Wang, J.; Xu, C.; Zhang, J.; Zhong, R. Big Data Analytics for Intelligent Manufacturing Systems: A Review. J. Manuf. Syst. 2022, 62, 738–752. [Google Scholar] [CrossRef]
  2. Loutas, T.H.; Roulias, D.; Georgoulas, G. Remaining Useful Life Estimation in Rolling Bearings Utilizing Data-Driven Probabilistic E-Support Vectors Regression. IEEE Trans. Reliab. 2013, 62, 821–832. [Google Scholar] [CrossRef]
  3. McDonald, G.L.; Zhao, Q. Multipoint Optimal Minimum Entropy Deconvolution and Convolution Fix: Application to Vibration Fault Detection. Mech. Syst. Signal Process. 2017, 82, 461–477. [Google Scholar] [CrossRef]
  4. Antoni, J. Fast Computation of the Kurtogram for the Detection of Transient Faults. Mech. Syst. Signal Process. 2007, 21, 108–124. [Google Scholar] [CrossRef]
  5. Antoni, J.; Xin, G.; Hamzaoui, N. Fast Computation of the Spectral Correlation. Mech. Syst. Signal Process. 2017, 92, 248–277. [Google Scholar] [CrossRef]
  6. Randall, R.B.; Antoni, J.; Gryllias, K. Alternatives to Kurtosis as an Indicator of Rolling Element Bearing Faults. In Proceedings of the ISMA 2016—International Conference on Noise and Vibration Engineering and USD2016—International Conference on Uncertainty in Structural Dynamics, Leuven, Belgium, 19–21 September 2016. [Google Scholar]
  7. Ding, Y.; He, W.; Chen, B.; Zi, Y.; Selesnick, I.W. Detection of Faults in Rotating Machinery Using Periodic Time-Frequency Sparsity. J. Sound Vib. 2016, 382, 357–378. [Google Scholar] [CrossRef]
  8. Feng, Z.; Zhou, Y.; Zuo, M.J.; Chu, F.; Chen, X. Atomic Decomposition and Sparse Representation for Complex Signal Analysis in Machinery Fault Diagnosis: A Review with Examples. Measurement 2017, 103, 106–132. [Google Scholar] [CrossRef]
  9. Chen, Z. Rolling Bearing Fault Diagnosis with Compressed Signals Based on Hybrid Compressive Sensing. J. Vibroeng. 2022, 24, 18–29. [Google Scholar] [CrossRef]
  10. Yang, B.; Liu, R.; Chen, X. Fault Diagnosis for a Wind Turbine Generator Bearing via Sparse Representation and Shift-Invariant K-SVD. IEEE Trans. Ind. Inform. 2017, 13, 1321–1331. [Google Scholar] [CrossRef]
  11. Arcos Jiménez, A.; Zhang, L.; Gómez Muñoz, C.Q.; García Márquez, F.P. Maintenance Management Based on Machine Learning and Nonlinear Features in Wind Turbines. Renew. Energy 2020, 146, 316–328. [Google Scholar] [CrossRef]
  12. Senanayaka, J.S.L.; Van Khang, H.; Robbersmyr, K.G. Toward Self-Supervised Feature Learning for Online Diagnosis of Multiple Faults in Electric Powertrains. IEEE Trans. Ind. Inform. 2021, 17, 3772–3781. [Google Scholar] [CrossRef]
  13. Sepulveda, N.E.; Sinha, J. Parameter Optimisation in the Vibration-Based Machine Learning Model for Accurate and Reliable Faults Diagnosis in Rotating Machines. Machines 2020, 8, 66. [Google Scholar] [CrossRef]
  14. Espinoza-Sepulveda, N.; Sinha, J. Mathematical Validation of Experimentally Optimised Parameters Used in a Vibration-Based Machine-Learning Model for Fault Diagnosis in Rotating Machines. Machines 2021, 9, 155. [Google Scholar] [CrossRef]
  15. Luwei, K.C.; Sinha, J.K.; Yunusa-Kaltungo, A.; Elbhbah, K. Data Fusion of Acceleration and Velocity Features (DFAVF) Approach for Fault Diagnosis in Rotating Machines. MATEC Web Conf. 2018, 211, 21005. [Google Scholar] [CrossRef]
  16. Wodecki, J.; Kruczek, P.; Bartkowiak, A.; Zimroz, R.; Wyłomańska, A. Novel Method of Informative Frequency Band Selection for Vibration Signal Using Nonnegative Matrix Factorization of Spectrogram Matrix. Mech. Syst. Signal Process. 2019, 130, 585–596. [Google Scholar] [CrossRef]
  17. Wodecki, J.; Michalak, A.; Zimroz, R.; Wyłomańska, A. Separation of Multiple Local-Damage-Related Components from Vibration Data Using Nonnegative Matrix Factorization and Multichannel Data Fusion. Mech. Syst. Signal Process. 2020, 145, 106954. [Google Scholar] [CrossRef]
  18. Liang, L.; Ding, X.; Wen, H.; Liu, F. Impulsive Components Separation Using Minimum-Determinant KL-Divergence NMF of Bi-Variable Map for Bearing Diagnosis. Mech. Syst. Signal Process. 2022, 175, 109129. [Google Scholar] [CrossRef]
  19. Liang, L.; Shan, L.; Liu, F.; Niu, B.; Xu, G. Sparse Envelope Spectra for Feature Extraction of Bearing Faults Based on NMF. Appl. Sci. 2019, 9, 755. [Google Scholar] [CrossRef]
  20. He, Q.; Ding, X. Time-Frequency Manifold for Machinery Fault Diagnosis. In Smart Sensors, Measurement and Instrumentation; Structural Health Monitoring; Springer: Cham, Switzerland, 2017; Volume 26, pp. 131–154. [Google Scholar]
  21. Harshman, R. A Foundations of the PARAFAC Procedure: Models and Conditions for an “Explanatory” Multimodal Factor Analysis. UCLA Work. Pap. Phon. 1970, 16, 1–84. [Google Scholar]
  22. Tucker, L.R. Some Mathematical Notes on Three-Mode Factor Analysis. Psychometrika 1966, 31, 279–311. [Google Scholar] [CrossRef]
  23. Cichocki, A.; Mandic, D.; De Lathauwer, L.; Zhou, G.; Zhao, Q.; Caiafa, C.; PHAN, H.A. Tensor Decompositions for Signal Processing Applications: From Two-Way to Multiway Component Analysis. IEEE Signal Process. Mag. 2015, 32, 145–163. [Google Scholar] [CrossRef]
  24. Miron, S.; Zniyed, Y.; Boyer, R.; de Almeida, A.L.F.; Favier, G.; Brie, D.; Comon, P. Tensor Methods for Multisensor Signal Processing. IET Signal Process. 2020, 14, 693–709. [Google Scholar] [CrossRef]
  25. Zhao, H.; Zhang, W.; Wang, G. Fault Diagnosis Method for Wind Turbine Rolling Bearings Based on Hankel Tensor Decomposition. IET Renew. Power Gener. 2019, 13, 220–226. [Google Scholar] [CrossRef]
  26. He, Z.; Shao, H.; Cheng, J.; Zhao, X.; Yang, Y. Support Tensor Machine with Dynamic Penalty Factors and Its Application to the Fault Diagnosis of Rotating Machinery with Unbalanced Data. Mech. Syst. Signal Process. 2020, 141, 106441. [Google Scholar] [CrossRef]
  27. Hu, C.; Wang, Y. Multidimensional Denoising of Rotating Machine Based on Tensor Factorization. Mech. Syst. Signal Process. 2019, 122, 273–289. [Google Scholar] [CrossRef]
  28. Liang, L.; Wen, H.; Liu, F.; Li, G.; Li, M. Feature Extraction of Impulse Faults for Vibration Signals Based on Sparse Non-Negative Tensor Factorization. Appl. Sci. 2019, 9, 3642. [Google Scholar] [CrossRef]
  29. Wen, H.; Liang, L.; Niu, B.; Shan, L.; Li, M.; Li, G. The Separation of Vibration Components Based on Sparse Nonnegative Tensor Factorization. In Smart Innovation, Systems and Technologies: Advances in Asset Management and Condition Monitoring; Springer: Berlin/Heidelberg, Germany, 2020; pp. 1297–1306. ISBN 9783030577445. [Google Scholar]
  30. Zhou, G.; Cichocki, A.; Zhao, Q.; Xie, S. Efficient Nonnegative Tucker Decompositions: Algorithms and Uniqueness. IEEE Trans. Image Process. 2015, 24, 4990–5003. [Google Scholar] [CrossRef]
  31. Zhou, G.; Cichocki, A.; Xie, S. Fast Nonnegative Matrix/Tensor Factorization Based on Low-Rank Approximation. IEEE Trans. Signal Process. 2012, 60, 2928–2940. [Google Scholar] [CrossRef]
  32. Wang, B.; Lei, Y.; Li, N.; Li, N. A Hybrid Prognostics Approach for Estimating Remaining Useful Life of Rolling Element Bearings. IEEE Trans. Reliab. 2018, 69, 401–412. [Google Scholar] [CrossRef]
  33. Kim, Y.D.; Choi, S. Nonnegative Tucker Decomposition. In Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition 2007, Minneapolis, MN, USA, 17–22 June 2007. [Google Scholar]
  34. Lee, D.D.; Seung, H.S. Algorithms for Non-Negative Matrix Factorization. Adv. Neural Inf. Process. Syst. 2001, 13, 556–562. [Google Scholar]
  35. Phan, A.H.; Cichocki, A. Extended HALS Algorithm for Nonnegative Tucker Decomposition and Its Applications for Multiway Analysis and Classification. Neurocomputing 2011, 74, 1956–1979. [Google Scholar] [CrossRef]
  36. Cichocki, A.; Zdunek, R.; Amari, S.I. Hierarchical ALS Algorithms for Nonnegative Matrix and 3D Tensor Factorization. In Proceedings of the 7th International Conference, ICA 2007, London, UK, 9–12 September 2007; Volume 4666, pp. 169–176. [Google Scholar] [CrossRef]
  37. Kim, H.; Park, H. Nonnegative Matrix Factorization Based on Alternating Nonnegativity Constrained Least Squares and Active Set Method. SIAM J. Matrix Anal. Appl. 2008, 30, 713–730. [Google Scholar] [CrossRef]
  38. Antoni, J.; Randall, R.B. Differential Diagnosis of Gear and Bearing Faults. J. Vib. Acoust. Trans. ASME 2002, 124, 165–171. [Google Scholar] [CrossRef]
Figure 1. Schematic of the experimental rig with accelerometers locations.
Figure 1. Schematic of the experimental rig with accelerometers locations.
Machines 10 00694 g001
Figure 2. Trends of the overall bearing vibration by the acceleration RMS values and Kurtosis: (a) the horizontal-channel vibration and (b) the vertical-channel vibration.
Figure 2. Trends of the overall bearing vibration by the acceleration RMS values and Kurtosis: (a) the horizontal-channel vibration and (b) the vertical-channel vibration.
Machines 10 00694 g002
Figure 3. An illustration of the tensor decomposition model for TF tensor factorization.
Figure 3. An illustration of the tensor decomposition model for TF tensor factorization.
Machines 10 00694 g003
Figure 4. Flowchart of the diagnostic methodology.
Figure 4. Flowchart of the diagnostic methodology.
Machines 10 00694 g004
Figure 5. (a) The horizontal acceleration signal and (b) its amplitude spectrum in Case A.
Figure 5. (a) The horizontal acceleration signal and (b) its amplitude spectrum in Case A.
Machines 10 00694 g005
Figure 6. Time−frequency representations of the analysed data: (a) STFT spectrogram. (b) Time−frequency tensor reorganized from the spectrogram segments, 𝓨 256 × 14640 × 8 .
Figure 6. Time−frequency representations of the analysed data: (a) STFT spectrogram. (b) Time−frequency tensor reorganized from the spectrogram segments, 𝓨 256 × 14640 × 8 .
Machines 10 00694 g006
Figure 7. Factors of LRANTD: (a) the Spectral Factor, (b) the Temporal Factor, consisting of 4 temporal subfactors with kurtosis of 2.98, 4.1, 15.7, and 13.3, respectively, (c) the Segment-wise Factor with 2 subfactors, and (d) the Core Tensor represented via its frontal slice w.r.t. SWS1, and (e) that w.r.t. SWS2.
Figure 7. Factors of LRANTD: (a) the Spectral Factor, (b) the Temporal Factor, consisting of 4 temporal subfactors with kurtosis of 2.98, 4.1, 15.7, and 13.3, respectively, (c) the Segment-wise Factor with 2 subfactors, and (d) the Core Tensor represented via its frontal slice w.r.t. SWS1, and (e) that w.r.t. SWS2.
Machines 10 00694 g007
Figure 8. Case A result via LRANTD: The amplitude spectrum of TS 4 with clear outer-race fault component.
Figure 8. Case A result via LRANTD: The amplitude spectrum of TS 4 with clear outer-race fault component.
Machines 10 00694 g008
Figure 9. (a) The horizontal acceleration signal, (b) the amplitude spectrum, and (c) the corresponding STFT spectrogram in Case B.
Figure 9. (a) The horizontal acceleration signal, (b) the amplitude spectrum, and (c) the corresponding STFT spectrogram in Case B.
Machines 10 00694 g009
Figure 10. Factors of LRA-NTD: (a) Spectral Factor, (b) Temporal Factor, (c) Segment-wise Factor, and (d) the Core Tensor represented via the frontal slice 1 and (e) the frontal slice 2.
Figure 10. Factors of LRA-NTD: (a) Spectral Factor, (b) Temporal Factor, (c) Segment-wise Factor, and (d) the Core Tensor represented via the frontal slice 1 and (e) the frontal slice 2.
Machines 10 00694 g010
Figure 11. Case B results via LRANTD: (a) The amplitude spectrum of TS1 and (b) that of TS2.
Figure 11. Case B results via LRANTD: (a) The amplitude spectrum of TS1 and (b) that of TS2.
Machines 10 00694 g011
Figure 12. (a) The horizontal acceleration signal, (b) the amplitude spectrum, and (c) the corresponding STFT spectrogram in Case C.
Figure 12. (a) The horizontal acceleration signal, (b) the amplitude spectrum, and (c) the corresponding STFT spectrogram in Case C.
Machines 10 00694 g012
Figure 13. Factors of LRANTD: (a) Spectral Factor, (b) Temporal Factor, (c) Segment-wise Factor, and (d) the Core Tensor represented via the frontal slice 1 and (e) the frontal slice 2.
Figure 13. Factors of LRANTD: (a) Spectral Factor, (b) Temporal Factor, (c) Segment-wise Factor, and (d) the Core Tensor represented via the frontal slice 1 and (e) the frontal slice 2.
Machines 10 00694 g013
Figure 14. Case C results via LRANTD. (a) The amplitude spectrum of TS1 and (b) that of TS3.
Figure 14. Case C results via LRANTD. (a) The amplitude spectrum of TS1 and (b) that of TS3.
Machines 10 00694 g014
Figure 15. The amplitude spectrum attained by the LRANTD method in Case A (Proposed method).
Figure 15. The amplitude spectrum attained by the LRANTD method in Case A (Proposed method).
Machines 10 00694 g015
Figure 16. Smoothed envelope spectrum (with 90% overlap and resolution at 1.5625 Hz).
Figure 16. Smoothed envelope spectrum (with 90% overlap and resolution at 1.5625 Hz).
Machines 10 00694 g016
Figure 17. Results of Fast Kurtogram: (a) Kurtogram with the maximum kurtosis, 7.1, found at level 3, with band width 1600 Hz, and the central frequency at 12,000 Hz. (b) The envelope spectrum of the filtered signal using the aforementioned optimal filter band.
Figure 17. Results of Fast Kurtogram: (a) Kurtogram with the maximum kurtosis, 7.1, found at level 3, with band width 1600 Hz, and the central frequency at 12,000 Hz. (b) The envelope spectrum of the filtered signal using the aforementioned optimal filter band.
Machines 10 00694 g017
Figure 18. The envelope spectrum of the 1st IMF.
Figure 18. The envelope spectrum of the 1st IMF.
Machines 10 00694 g018
Figure 19. (ac): The cross-spectra based on LRANTD fusion for Cases A to C. From Case A to Case C, the low-rank dimensions are, respectively, ( 3 , 4 , 2 ) , ( 4 , 4 , 2 ) , and ( 3 , 3 , 2 ) .
Figure 19. (ac): The cross-spectra based on LRANTD fusion for Cases A to C. From Case A to Case C, the low-rank dimensions are, respectively, ( 3 , 4 , 2 ) , ( 4 , 4 , 2 ) , and ( 3 , 3 , 2 ) .
Machines 10 00694 g019aMachines 10 00694 g019b
Figure 20. (ac): Results based on the cross-power spectral density of the two-channel envelope signals, from Case A to Case C.
Figure 20. (ac): Results based on the cross-power spectral density of the two-channel envelope signals, from Case A to Case C.
Machines 10 00694 g020aMachines 10 00694 g020b
Table 1. Basic notations and properties.
Table 1. Basic notations and properties.
0The vector/matrix with all its entries being zero.
U   0 , u r Nonnegative matrix U and its r -th column vector.
U ( n ) I n × R n , u r ( n ) I n × 1 Mode- n matrix and its r -th column vector, r = { 1 , 2 , , R n } .
𝓨 I 1 × I 2 × I N ,   Y ( n ) I n × i n I i N-way arrays or tensor and its mode- n matricizaton.
𝓖 R 1 × R 2 × R N , G ( n ) R n × i n R i Core tensor and its mode- n matricizaton.
( R 1 , R 2 , , R N ) Low-rank dimensions for NTD, R n < I n
I The index set of positive integers, I = { 1 , 2 , , N } .
𝓨 = 𝓨 × n U ( n ) Mode- n product, 𝓨 I 1 × I n 1 R n I n + 1 × I N
K = P Q Kronecker product of P I 1 × I 2 and Q R 1 × R 2 yields K = [ p i 1 i 2 Q ] I 1 R 1 × I 2 R 2
,   U U Element-wise multiplication, division
𝓖 × n 𝕀 U ( n ) is equal to 𝓖 × 1 U ( 1 ) × N U ( N ) . If 𝓨 = 𝓖 × n 𝕀 U ( n ) , then Y ( n ) = U ( n ) G ( n ) ( U ( N ) U ( n + 1 ) U ( n 1 ) U ( 1 ) ) T .
y i L , i = 1 , , M The i -th channel acceleration with L points in length, within M sensors measurement.
y i , j L s , y i , j ( j ) The j -th segment resampled from the i -th channel signal, the j -th entry of   y i , j .
Y i , j I 1 × I 2 ,   Y ( i 1 , i 2 ) Spectrogram matrix with I 1 frequency bins and I 2 time points, ( i 1 , i 2 ) th -entry of Y
S ,   L s , I 3 The number of data segments resampled from one signal, the length of each segmentation, the total number of segments
Table 2. Test bearing specifications ( f r : shaft speed, f c : cage frequency, f b : ball spin frequency (BSF), f o : ball pass frequency outer (BPFO), f i : bass pass frequency inner (BPFI).).
Table 2. Test bearing specifications ( f r : shaft speed, f c : cage frequency, f b : ball spin frequency (BSF), f o : ball pass frequency outer (BPFO), f i : bass pass frequency inner (BPFI).).
ParameterValueCharacteristicsValue
Mean diameter34.55 mm f r 35 Hz
Ball diameter7.92 mm f c 13.49 Hz
Number of balls8 f b 72.33 Hz
Contact angle f o 107.91 Hz
Load rating12 kN f i 172.09 Hz
Table 3. Characteristic spectral amplitude in g (with its ratio to the mean background noise) for Case A.
Table 3. Characteristic spectral amplitude in g (with its ratio to the mean background noise) for Case A.
Case A (Single-Channel)Case A (Sensor-Fusion)
ComponentEnvelope SpectrumFast
Kurtogram
EMDLRANTD(Envelope)
CPSD
LRANTD
Shaft Speed fr0.1115 (5.21)0.0138 (6.27)0.1181 (5.15)0.0211 (2.63)0.043 (5.58)0.0214 (2.05)
Cage Freq. fc0.0819 (3.82)0.0109 (4.95)0.0866 (3.78)0.0292 (3.65)0.0362 (4.70)0.0416 (4.00)
BSF fb0.0202 (0.94)0.0019 (0.86)0.0244 (1.06))0.0090 (1.12)0.0103 (1.33)0.0069 (0.66)
BPFOfo0.0656 (3.06)0.0096 (4.36)0.0661 (2.88)0.0375 (4.68)0.0332 (4.31)0.0597 (5.74) *
BPFO 2x0.0388 (1.81)0.0047 (2.13)0.0348 (1.51)0.0211 (2.63)0.0226 (2.93)0.0294 (2.82) *
BPFO 3x0.0230 (1.07)0.0021 (0.95)0.018 (0.78)0.0060 (0.75)0.0080 (1.03)0.0118 (1.13) *
BPFI fi0.0209 (0.97)0.0043 (1.95)0.0218 (0.95)0.0123 (1.53)0.0117 (1.51)0.0115 (1.10)
Background Noise0.02140.00220.02290.0080.00770.0104
RemarkBPFO Baseline+1.30 (+42.48%)−0.18 (−5.88%)+1.62 (+52.94%)BPFO Baseline+1.43 (+33.18%)
* The best amplitude ratio for the specific spectral component in each case.
Table 4. Characteristic spectral amplitude in g (with its ratio to the mean background noise) for Case B and Case C.
Table 4. Characteristic spectral amplitude in g (with its ratio to the mean background noise) for Case B and Case C.
Case B (Sensor-Fusion)Case C (Sensor-Fusion)
ComponentCPSDLRANTDCPSDLRANTD
Shaft Speed fr0.1244 (7.36)0.0217 (2.08)0.0958 (6.02)0.0394 (4.37)
Cage Freq. fc0.112 (6.62)0.0205 (1.97)0.0562 (3.53)0.0157 (1.74)
BSF fb0.0405 (2.39)0.0184 (1.76)0.0404 (2.54)0.0137 (1.52)
BPFOfo0.0969 (5.73)0.1134 (10.90) *0.2335 (14.68)0.2525 (28.05) *
BPFO 2x0.0451 (2.66)0.0488 (4.69) *0.1021 (6.42) *0.0464 (5.15)
BPFO 3x0.0170 (1.01)0.0232 (2.23) *0.0690 (4.33)0.0405 (4.50) *
BPFI fi0.0128 (0.75)0.0146 (1.40)0.0126 (0.79)0.0054 (0.60)
Background Noise0.01690.01040.01590.009
RemarkBPFO Baseline+5.17 (+90.23%)BPFO Baseline+13.37 (+91.07%)
* The best amplitude ratio for the specific spectral component in each case.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Wen, H.; Zhang, L.; Sinha, J.K. Adaptive Band Extraction Based on Low Rank Approximated Nonnegative Tucker Decomposition for Anti-Friction Bearing Faults Diagnosis Using Measured Vibration Data. Machines 2022, 10, 694. https://doi.org/10.3390/machines10080694

AMA Style

Wen H, Zhang L, Sinha JK. Adaptive Band Extraction Based on Low Rank Approximated Nonnegative Tucker Decomposition for Anti-Friction Bearing Faults Diagnosis Using Measured Vibration Data. Machines. 2022; 10(8):694. https://doi.org/10.3390/machines10080694

Chicago/Turabian Style

Wen, Haobin, Long Zhang, and Jyoti K. Sinha. 2022. "Adaptive Band Extraction Based on Low Rank Approximated Nonnegative Tucker Decomposition for Anti-Friction Bearing Faults Diagnosis Using Measured Vibration Data" Machines 10, no. 8: 694. https://doi.org/10.3390/machines10080694

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop