Time–frequency analysis of torsional vibration signals based on the improved complete ensemble empirical mode decomposition with adaptive noise,Robust independent component analysis,and Prony’s methods

Abstract

The torsional vibration signals of rotating shafts are multimodal non-stationary noisy signals. The harmonics and attenuation characteristics of these non-stationary signals cannot be obtained effectively by the ordinary short-time Fourier transform algorithm. Although Prony analysis can accurately fit and identify the characteristic coefficients of such non-stationary signals, it is still sensitive to noise. In this article, we propose a system for denoising and identification of torsional vibration signals. In particular, the improved complete ensemble empirical mode decomposition with adaptive noise (ICEEMDAN), wavelet transform, and robust independent component analysis methods are used to denoise the torsional vibration signals, and then, Prony analysis is used to obtain the characteristic parameters of these signals. The proposed algorithm has good denoising performance and it can improve the identification accuracy and reduce the order of the Prony analysis.

Keywords

Improved complete ensemble empirical mode decomposition with adaptive noise (ICEEMDAN)Prony analysis robust independent component analysis algorithm torsional vibration analysis wavelet transform

1. Introduction

In mechanical systems, the torsional vibration signals (TVS) of rotating shafts are multimodal non-stationary noisy signals. Taradai et al. (2018) have described three experimental methods to determine the lowest natural frequencies of torsional vibration excited in the shaft system, Lv et al. (2016) used a T-like beam structure to measure the torsional vibration, Xue and Howard (2018) used the TVS as a diagnostic tool for planetary gearbox faults detection, and Janssens et al. (2011) developed a zebra tape butt joint algorithm for torsional vibrations. Although harmonic analysis of these signals based on the short-time Fourier transform (STFT) can be carried out, the higher frequency resolution can be obtained only through a large sampling window. Also, it is difficult to obtain the attenuation characteristics of these non-stationary signals. The wavelet transform (WT) (Daubechies and Heil 1998; Daubechies 2002) and the Hilbert–Huang transform (HHT) (Huang et al., 1998; Huang and Wu, 2008; Wu and Huang, 2009; Yeh et al., 2010) are good alternatives for the analysis of such non-stationary signals. However, the results of wavelet analysis depend on the selection of the wavelet basis and are affected by the frequency-domain aliasing, and hence are not accurate enough. For the HHT, a signal is decomposed into intrinsic mode functions (IMF) by empirical mode decomposition (EMD). The signal parameters are thus obtained by means of Hilbert transform (HT). This transform can process non-stationary signals, and the generated parameters are independent of the selection of basic functions, but there are problems of modal aliasing and limited accuracy. Prony analysis can accurately fit and identify the characteristic coefficients of non-stationary signals. However, for noisy signals, the fitting results of this analysis are typically distorted (Marple, 1989; Okamoto et al., 1996). The WT combined with EMD can be used for signal denoising, but the obtained signal-to-noise ratios (SNR) are usually limited. The robust independent component analysis (ICA) algorithm (Zarzoso and Comon, 2010) exploits the mutual statistical independence between mixed source signals to blindly separate such signals and uses the historical data similarities to match the separation results. Nejad and Rizi (2018) used the joint method of ICA, ICEEMDAN, Wavelet Shrinkage to extract and anlysis an adaptive FECG. Nevertheless, the signal transmission in the shaft torsional vibration testing is random and uncertain, thus it is difficult for this algorithm to get the ideal independent signal components. In this article, we developed a joint algorithm to analyze the signal. The signal can be divided by EMD into IMFs, which can be used as the input for the robust ICA algorithm, which can produce physically meaningful signals. Then, wavelet thresholding is used to denoise the robust ICA output and obtain a higher SNR. The denoised signals can be analyzed by the Prony method to obtain higher identification accuracies.

2. Signal parameter identification based on the Prony method

2.1. The Prony method

The Prony method can be used to identify the modal parameters of a signal as well as the signal frequency, amplitude, phase, and attenuation factor (Marple, 1989; Okamoto et al., 1996). In this work, we essentially use a linear combination of exponential functions to fit the equal-interval sampled data, for which the fitting coefficients can be found by the singular value decomposition and the total least-square regression (SVD-TLS). The steps of this approach are as follows:

Suppose the time series of the TVS can be represented as the following linear combination of complex exponentials

\hat{x (n)} = \sum_{m = 1}^{p} b_{m} z_{m}^{n}, (n = 0,1, \dots, N - 1)

(1)

where

p

is the fitting order and

N \geq 2 p

is the number of data samples. The coefficient

b_{m}

is the retention number and is given by

b_{m} = A_{m} \exp (j θ_{m})

, where

A_{m}

and

θ_{m}

are amplitude and phase parameters, respectively. The base

z_{m}

is given by

z_{m} = \exp [(α_{m} + j 2 π f_{m}) Δ t]

, where

α_{m}

is the attenuation factor,

f_{m}

is the oscillation frequency, and

Δ t

is the sampling interval.

To make a good fit to the actual signal, minimize the sum of the squared errors, that is

\min (ε = {\sum_{n = 0}^{N - 1} | x (n) - \hat{x (n)} |}^{2})

(2)

The solution of the minimization problem in (2) can be regarded as the homogeneous solution of a difference equation with constant coefficients. Hence, a linear prediction model can be constructed as

x (k) = a_{1} x (k - 1) + a_{2} x (k - 2) + \dots + a_{n} x (k - n)

(3)

For sampled data, the linear prediction model can be written in matrix form as

[\begin{matrix} x (n) & \dots & x (1) \\ \dots & \dots \\ x (N - 1) & \dots & x (N - n) \end{matrix}] [\begin{matrix} a_{1} \\ \dots \\ a_{n} \end{matrix}] = [\begin{matrix} x (n + 1) \\ \dots \\ x (N) \end{matrix}]

(4)

As for the given $\dot{x} (n)$ values, when $N = 2 n$ , the number of unknowns is equal to the number of equations, and the system of equations in (4) can be solved directly. If $N > 2 n$ , the SVD-TLS regression can be used to obtain the coefficients $a_{i}$ .

By computing the eigenvalues of the characteristic equation of the system in (4)

z^{n} - (a_{1} z^{n - 1} + a_{2} z^{n - 2} + \dots + a_{n} z^{0}) = 0

(5)

The eigenvalues can be replaced, and the coefficients $b_{i}$ can be obtained by solving the least-squares system

[\begin{matrix} 1 & \dots & 1 \\ z_{1} & \dots & z_{n} \\ \dots & \dots & \dots \\ z_{1}^{N - 1} & \dots & z_{n}^{N - 1} \end{matrix}] [\begin{matrix} b_{1} \\ b_{2} \\ \dots \\ b_{n} \end{matrix}] = [\begin{matrix} x (1) \\ x (2) \\ \dots \\ x (N) \end{matrix}]

(6)

By solving the amplitude $A_{i}$ , the phase $θ_{i}$ , the frequency $f_{i}$ , and the attenuation factor $α_{i}$ of each component as

A_{i} = | b_{i} |

(7)

θ_{i} = \arctan [\frac{Im (b_{i})}{Re (Z_{i})}]

(8)

a_{i} = \frac{\ln | Z_{i} |}{Δ t}

(9)

f (x) = \frac{\arctan [(Im (Z_{i}) / Re (Z_{i}))]}{2 π Δ t}

(10)

2.2. Prony analysis of TVS

The damped TVS of a rotating shaft can be given as

s (t) = A e^{- α t} \cos (ω t + θ_{0}) = A e^{- ξ ω_{0} t} \cos (ω_{0} \sqrt{1 - ξ^{2}} t + θ_{0})

(11)

In this expression, $ω$ and $ω_{0}$ are the damped and undamped (natural) frequencies, respectively, $A$ is the oscillatory amplitude, $ξ$ is the damping ratio, and $θ_{0}$ is the initial phase. The damping frequency is related to these parameters by

{\begin{array}{l} ξ ω_{0} = - α \\ ω_{0} \sqrt{1 - ξ^{2}} = ω \end{array}

(12)

ξ = \frac{- α}{\sqrt{ω^{2} + α^{2}}}

(13)

For example, consider the TVS

φ (t) = 1 e^{- 0.5 t} \sin (100 π t) + 0.5 e^{- 2 t} \sin (90 π t)

(14)

We show here the experimental results of applying the Prony analysis to the signal in (14). For the experimental setup, we used a ThinkPad X270 machine with an Intel Core i7-7500U CPU, a 2.7-GHz processor, an 8 GB memory, a 64-bit Windows 10 operating system, and a MATLAB R2018b environment.

The TVS (14) is analyzed with a sampling rate of 1 kHz and a sampling duration of 1 s. The fitting results are shown in Table 1.

Table 1.

Prony analysis results for the torsional vibration signals in (14).

Order	Time(s)	$f_{1}$	$f_{2}$	$A_{1}$	$A_{2}$	$α_{1}$	$α_{2}$	$θ_{1}$	$θ_{2}$
2	0.0102	45.046	50.05	0.5	1	−2.002	−0.5005	−1.5708	−1.5708
3	0.0112	45.046	50.05	0.5	1	−2.002	−0.5005	−1.5708	−1.5708
10	0.0226	45.046	50.05	0.5	1	−2.002	−0.5005	−1.5708	−1.5708

We also investigate the effect of noise by contaminating the signal in (14) with white Gaussian noise to produce SNR of 20 dB and 30 dB, respectively, where the sampling rate is 1 kHz and the sampling time is 1 s. The fitting results for the two noise levels are shown in Tables 2 and 3, respectively.

Table 2.

Prony analysis results for the signal in (14) with 20-dB signal-to-noise ratios.

Order	Time(s)	$f_{1}$	$f_{2}$	$A_{1}$	$A_{2}$	$α_{1}$	$α_{2}$	$θ_{1}$	$θ_{2}$
10	0.034	47.928	49.71	2.25	1.78	−47.07	−4.34	−2.861	−1.072
50	0.22	44.99	50.064	0.012	0.0167	−17.39	−9.556	−1.174	−2.114
100	0.3796	45.054	50.053	0.497	1.013	−1.975	−0.5244	−1.593	−1.5837

Table 3.

Prony analysis results for the signal in (14) with 30-dB signal-to-noise ratios.

Order	Time(s)	$f_{1}$	$f_{2}$	$A_{1}$	$A_{2}$	$α_{1}$	$α_{2}$	$θ_{1}$	$θ_{2}$
10	0.025	45.21	50.12	0.682	1.307	−6.934	−1.326	−1.858	−1.723
50	0.3528	45.06	50.04	0.4998	1.016	−1.998	−0.534	−1.598	−1.567
100	0.4633	45.05	50.05	0.5023	1.000	−2.026	−0.500	−1.567	−1.57

As can be seen from Tables 1, 2, and 3, when the signal is noisy, the Prony method with lower orders shows inferior identification performance. The performance generally improves as the Prony order increases. Also, for the signal with high SNR, Prony fitting accuracy is higher under the same sampling rate and sampling time.

3. Prony analysis with ICEEMDAN-based denoising

3.1. The ICEEMDAN method

The EMD is a time–frequency signal processing method (Huang et al., 1998, Huang and Wu, 2008; Wu and Huang, 2009; Yeh et al., 2010) that decomposes nonlinear non-stationary signals into IMF. However, there are endpoint effects and pattern aliasing problems with this approach. Wu and Huang (2009) developed the noise-assisted ensemble empirical mode decomposition (EEMD) to solve the above-mentioned EMD problems. The EEMD method performs EMD decomposition on multiple contaminations of the original signal with white noise and averages the results of the multiple decompositions to obtain the final IMF. The quality of the reconstructed signal is affected by the addition of noise in EEMD. Yeh et al. (2010) adopt the complementary empirical mode decomposition (CEEMD) method that mixes two complementary white noise signals with the input signal, decomposes the noisy signals, and performs EMD analysis. This method can significantly reduce the influence of noise on the reconstructed signal. With the addition of noise in the EEMD and CEEMD methods, the EMD computational complexity increases. Also, since a lot of pseudo components may be produced after EMD completion, the EMD-generated components may need to be reanalyzed and reprocessed. For the actual EEMD implementation, the number of iterations is limited. Consequently, it is likely that the decomposition iteration limit is reached while the IMF termination condition is not met. Torres et al. (2011) proposed an improved complete EEMD with adaptive noise (ICEEMDAN), which can reduce decomposition times by improving the noise added in the EEMD analysis, and reducing the noise impact without affecting the signal quality. Terrien et al. (2011) estimated the noise content in the first IMF by energy statistics and extracted effective signals for noise estimation correction.

The steps of the ICEEMDAN algorithm are as follows:

Step 1: Calculate the local mean value of the signal $x_{i} (t) = x (t) + ε_{0} E_{1} (w_{i} (t))$ by EMD to obtain the first residual

r_{1} (t) = M (x_{i} (t))

(15)

Step 2: Calculate the first modal component of the first order $(j = 1)$

im f_{1} (t) = x (t) - r_{1} (t)

(16)

Step 3: Estimate a second residual as a local mean of the signal $r_{1} (t) + ε_{1} E_{2} (w_{i} (t))$ and define a second modal component

Step 4: Calculate the residual of order $j$ for $(j = 3, \dots, J)$

r_{j} (t) = M (r_{j - 1} (t) + ε_{j - 1} E_{j} (w_{i} (t)))

(17)

Step 5

Calculate the $j$ modal component

i m f_{j} (t) = r_{j - 1} (t) - r_{j} (t)

(18)

Step 6

Repeat Steps 4–6 until the termination condition is met.

3.2. Prony analysis after ICEEMDAN denoising

We added white Gaussian noise to the signal in (14), consequently that the SNR is 20 dB. The ICEEMDAN method was used to decompose this noisy signal. The decomposition results are shown in Figure 1.

Figure 1.

The improved complete ensemble empirical mode decomposition with adaptive noise–based decomposition of the signal in (14) with 20-dB signal-to-noise ratios.

As shown in Figure 1, the signal and the noise can be separated from the modal components by the ICEEMDAN-based decomposition. However, the decomposed signal is a time-varying inter-harmonic signal, whose frequencies are comparatively close. Thus, the ICEEMDAN-based signal decomposition is mixed to a certain extent.

In order to detect and remove the abnormal components, we calculate the permutation entropy (Bandt and Pompe, 2002; Yan and Gao, 2007; Zheng, et al., 2013) of the ICEEMDAN components. When the permutation entropy of a component is greater than 0.8, this component is considered abnormal and is subsequently discarded. After removing the abnormal components, the remaining components are synthesized. Compared with the signal in (14), the reconstructed signal has a SNR of 24.345 dB, a correlation coefficient of 0.9992, and a root-mean-square error (RMSE) of 0.417. Compared with the 20-dB noisy signal, the SNR increases from 20 to 24.345, the correlation coefficient is not significantly altered, and the RMSE increases. This indicates that the ICEEMDAN method can reduce the noise to a certain extent and remove the abnormal components based on the permutation entropy. However, the first discarded abnormal component may have essential information, and the residual components still contain noise after eliminating the abnormal components. This leads to the RMSE increase for the reconstructed signal. The Prony analysis results of the reconstructed signal are shown in Table 4.

Table 4.

Prony analysis results after improved complete ensemble empirical mode decomposition with adaptive noise–based denoising for the signal in (14) with 20-dB signal-to-noise ratios.

Order	$f_{1}$	$f_{2}$	$A_{1}$	$A_{2}$	$α_{1}$	$α_{2}$	$θ_{1}$	$θ_{2}$
200	45.0547	50.0512	0.4927	0.9997	−1.9492	−0.4963	−1.5946	−1.5767

4. Prony analysis with wavelet-based denoising

Weaver et al. (1991) proposed the wavelet-threshold noise reduction method. Then, Daubechies (2002) made a more refined analysis of this denoising method, which is simple, convenient, and widely applicable. The main steps of the wavelet denoising method are as follows:

Step 1: Transform the target signal into the wavelet domain;

Step2: Carry out threshold processing for the wavelet coefficients, mainly to remove the random noise; and

Step 3: Apply the inverse WT to obtain the denoised signal.

The wavelet thresholding is used to extract the ground signal from the noise-contaminated signal. The discrete Meyer wavelets are selected as the wavelet basis functions, and the decomposition tree has a height of three layers. When the heuristic threshold is 4.29, the denoising performance is relatively good. The reconstructed waveform is shown in Figure 2.

After wavelet denoising, the SNR is 24.87 dB, the correlation coefficient is 0.9998, and the RMSE is 0.0336. Compared with the ICEEMDAN method, the SNR and the correlation coefficient are improved, and the RMSE is reduced. This indicates that the denoising performance with wavelet thresholding is better than that of the ICEEMDAN algorithm. The reconstructed signal is analyzed by the Prony method and the results are shown in Table 5.

Figure 2.

The wavelet-based denoising output for the signal in (14) with 20-dB signal-to-noise ratios.

Table 5.

Prony analysis results after wavelet-threshold denoising for the signal in (14) with 20-dB signal-to-noise ratios.

Order	$f_{1}$	$f_{2}$	$A_{1}$	$A_{2}$	$α_{1}$	$α_{2}$	$θ_{1}$	$θ_{2}$
200	45.0548	50.0505	0.4958	0.9965	−1.9648	−0.4949	−1.5949	−1.5749

5. Time–frequency analysis based on the ICEEMDAN, WT, Prony, and robust ICA methods

5.1. Prony analysis with joint denoising based on ICEEMDAN and WT

For signal analysis, the modal components are first decomposed by ICEEMDAN and then denoised by wavelet thresholding. Again, a discrete Meyer wavelet basis and a wavelet tree of three layers are selected. The effect of the denoising of each modal component by using the heuristic threshold is shown in Figure 3.

Figure 3.

Signal decomposition and denoising with Improved complete ensemble empirical mode decomposition with adaptive noise and wavelet thresholding.

After the wavelet-threshold denoising, the SNR is 24 dB, the correlation coefficient is 0.998, and the RMSE is 0.037. Compared with the ICEEMDAN-based and wavelet-based methods, the correlation denoising index is slightly reduced. When the wavelet thresholding is used, many influencing factors are involved. When denoising the ICEEMDAN components, there is a modal aliasing problem in the ICEEMDAN decomposition. In the denoising of each modal component, the effective components are also removed. This leads to the unsatisfactory results.

5.2. Prony analysis with joint denoising based on ICEEMDAN, robust ICA, and WT

5.2.1. Robust ICA

Zarzoso et al. (2010) proposed a robust independent component analysis (robust ICA) algorithm based on the kurtosis and the optimal step size. In this algorithm, the optimal step size search is carried out by a kurtosis control function, the unmixing matrix is obtained, and the source signal is estimated.

Let the mixed noisy signal observations be $x$ . Without whitening, only the mean value of $x$ is zero, and the signal transformed by the separation matrix is $y = w_{i} x (i = 1,2, \dots, m)$ . The kurtosis is

K (w_{i}) = \frac{E {{| y |}^{4}} - 2 E^{2} {{| y |}^{2}} - {| E {{| y |}^{2}} |}^{2}}{E^{2} {{| y |}^{2}}}

(19)

In (19), $E {•}$ represents mathematical expectations. It can be seen from (19) that the kurtosis index is not sensitive to the input amplitude, that is, $K (λ w_{i}) = K (w_{i})$ , $\forall λ \neq 0$ . Therefore, without loss of generality, $w_{i}$ is standardized to a unit-norm input, that is, $w_{i} = 1$ . In the general ICA algorithm, it is necessary to whiten and center the data, but these steps are not needed in robust ICA. The sentinel control function is optimized by finding the optimal linear step size

μ_{opt} = \underset{μ}{argmax} | K (w_{i} + μ g) |

(20)

The typical search direction is the direction of the gradient $g$

g = \nabla_{w_{i}} K (w_{i})

(21)

This gradient can be computed as

\nabla_{w_{i}} K (w_{i}) = \frac{4}{E^{2} {{| y |}^{2}}} {E {{| y |}^{2} y x} - E {y x} E {y^{2}} - \frac{(E {{| y |}^{4}} - {| E {y^{2}} |}^{2}) E {y x}}{E {{| y |}^{2}}}}

(22)

Some other restrictive conditions are introduced into the exact linear search, which is very computationally intensive. Hence, this optimization approach is difficult to be applied in practice. For more practicality and accuracy, kurtosis can be converted into a functional relation in the step size where we seek to maximize the absolute value of kurtosis. The calculated step size value can be optimized by iterating to obtain the global optimal step size OS. The optimization steps are as follows:

Step 1: The coefficients of OS are calculated. The relationship between the step size and kurtosis is given by

p (μ) = \sum_{k = 0}^{4} a_{k} μ^{k}

(23)

In each iteration, the coefficients ${a_{k}}_{k = 0}^{4}$ can be obtained from the observed signal $x$ and the current values of $w_{i}$ and $g$ . The optimal value $μ_{opt}$ can be adjusted by standardizing the value of the gradient vector in advance.

Step 2: The roots ${μ^{k}}_{k = 1}^{4}$ of the extraction polynomial are computed. The four polynomial roots can be obtained by any standard algebraic method, such as Ferrari’s method.

Step 3: The maximum root of the polynomial in (23) can be found iteratively, and hence can be used to iterate through the OS coefficients calculated in Step 1

μ_{opt} = argmax | K (w_{i - 1}) + μ_{i} g |

(24)

Step 4: Update $w_{i}$ with the value of the optimized step size $μ_{opt}$ from Step 3

w_{i} = w_{i - 1} + μ_{opt} g

(25)

Step 5: Normalize $w_{i}$ by its norm

w_{i}^{*} = \frac{w_{i}}{w_{i}}

(26)

5.2.2. Prony analysis with joint denoising based on ICEEMDAN, robust ICA, and WT

First, the noisy signal is decomposed by the ICEEMDAN method. Then, the obtained components are then used as a multi-sensor input to the robust ICA algorithm. The independent components obtained by the robust ICA analysis are shown in Figure 4. After this analysis, the obtained signal is still irregular because of its time dependence, and there is modal aliasing in the ICEEMDAN outcome, but the shape of the output signal is obviously different from that of the ICEEMDAN input signal. The independent components are thus denoised by wavelet thresholding with a discrete Meyer wavelet basis, a wavelet tree of three layers, and a heuristic threshold. The denoised signal is then mixed with the mixed matrix to reconstruct the updated IMFs, as shown in Figure 5. The SNR is 25.22 dB, the RMSE is 0.0322, and the correlation coefficient is 0.9995. The denoising performance index is higher than that of the methods described above.

Figure 4.

Independent component analysis components for the signal in (14).

Figure 5.

Independent component analysis components after denoising.

The reconstructed signal is analyzed by the Prony analysis. The results are shown in Table 6. Indeed, the accuracy of the Prony analysis results has been slightly improved.

Table 6.

Prony analysis results after denoising based on the improved complete ensemble empirical mode decomposition with adaptive noise, wavelet thresholding, and robust independent component analysis methods. The tested signal is the one in (14) with 20-dB signal-to-noise ratios.

Order	$f_{1}$	$f_{2}$	$A_{1}$	$A_{2}$	$α_{1}$	$α_{2}$	$θ_{1}$	$θ_{2}$
200	45.0554	50.0520	0.4904	0.9998	−1.9445	−0.494	−1.5972	−1.581

The signal is segmented into 300, 300, and 400 ms. Then, the signal is analyzed by the method of the fourth experiment. The ICEEMDAN outcome for the former 300 ms is divided into noise and signal components as shown in Figure 6. (a). Because the two frequency segments of the test signal are in the frequency-doubling adjacent ranges, the EMD method is generally not suitable, and the ICEEMDAN components can be identified by the robust ICA algorithm. As a result, as shown in Figure 6(b), the obtained independent components are subjected to wavelet threshold denoising, whose results are shown in Figure 6(c). The denoised signal is multiplied with the mixed matrix to obtain the updated IMFs components, which are then added to complete the signal reconstruction. The residual error of the reconstructed and original signals is shown in Figure 6(d). The horizontal axis in Figure 6 is the time and the vertical axis is the signal amplitude. The SNR is 26.67 dB, the RMSE is 0.0315, and the correlation coefficient is 0.999. Compared with the results of the previous methods, the denoising index is obviously improved.

Figure 6.

Joint algorithm results after segmentation.

The reconstructed signal is analyzed by the Prony method. The results are shown in Table 7. The accuracy is slightly improved, and the Prony analysis order is reduced. The analyzed signal can be reasonably truncated after appropriate segmentation, and the performance of the denoising algorithm can be improved.

Table 7.

Prony analysis results after segmentation.

Order	$f_{1}$	$f_{2}$	$A_{1}$	$A_{2}$	$α_{1}$	$α_{2}$	$θ_{1}$	$θ_{2}$
50	44.999	50.0691	0.4985	1.0028	−2.035	−0.5047	−1.5425	−1.5947

6. Analysis of real TVS

We investigate here the analysis of real TVS in a power plant. As shown in Figure 7, the intermediate shaft of an induced draft fan is broken after its setup with a frequency converter.

Figure 7.

Cracks in the intermediate shaft of an induced draft fan.

After the shaft is replaced or repaired, the frequency converter is put back into operation. When the operating frequency reaches 22 Hz, the system vibration is significantly enhanced. As shown in Figure 8, the torque signal is collected through a wireless torque sensor with a sampling rate of 2 kHz.

Figure 8.

Torque signal of the induced draft fan unit.

As shown in Figure 9, the signal with a sampling period of 0.5s is captured and denoised by the proposed joint algorithm of this article. This denoising algorithm takes about 12.7s; Fig. 9(a) shows the ICEEMDAN decomposition results, while Fig.9(b) shows the robust ICA analysis results. Fig.9(c) exhibits the results after wavelet denoising, while Fig.9(d) shows the robust ICA reconstruction results.

Figure 9.

De-noising of real torsional vibration signals with the proposed joint algorithm.

After denoising, the denoised signal is analyzed by 200-order Prony and an analysis time of 0.63 s. The Prony components are ranked according to their correlation with the original signal. The top five components (in terms of correlation) are selected as shown in Table 8. The last column in this table shows the correlation coefficient between the original signal and the respective Prony component.

Table 8.

Prony analysis results.

NO.	$Frequency (f)$	$Amplititude (A)$	$Attenuation factor (α)$	$phase (θ)$	correlation coefficient (corr)
1	17.69	28.52	0.038	1.637	0.988
2	51.09	1.341	−2.498	−0.2162	0.0257
3	46.98	2.66	−9.736	−0.7	0.0251
4	177.3	0.7523	−1.075	−1.377	0.01999
5	111.66	1.1597	−1.384	−0.8451	0.01798

The inertia constant of the induced draft fan $(J_{i})$ is 5400 kg.m², while that of the motor $(J_{e})$ is 1900 kg.m². The elastic constant of the intermediate shaft $(e_{i e})$ is 1.68×10⁷ N.m/rad. If the system under consideration is regarded as a double-mass system, the natural torsional vibration frequency of the system $(f_{n})$ is given by

f_{n} = \frac{1}{2 π} \sqrt{\frac{e_{i e} (J_{i} + J_{e})}{J_{i} J_{e}}}

(27)

where

f_{n}

= 17.4 Hz. The first Prony component has a frequency of 17.69 Hz, which is quite close to the estimated

f_{n}

. As well, the attenuation factor of 0.038 is greater than 0. Consequently, the system is unstable.

7. Summary

For a multimodal noisy non-stationary signal, the signal can be decomposed by ICEEMDAN, and then, the abnormal component can be judged by permutation entropy and the signal is reconstructed after removing the abnormal component. Thus, signal denoising can be realized. While SNR is improved, the RMSE will increase. Also, while the noise modal component is deleted by the EMD method, the noise modal component still contains the effective signal that has not been decomposed, and the wavelet threshold is used to denoise the signal. The selection of wavelet basis function has a large workload and SNR is improved, and at the same time, the RMSE will not increase significantly. The selection of wavelet function, the number of decomposition layers, and the selection of threshold have influence on the denoising effect. The signal is first decomposed by ICEEMDAN and then denoised by wavelet threshold for each modal component. The method can improve the SNR and the RMSE does not increase significantly. But compared with wavelet threshold denoising, this method has no significant advantage. The decomposed signal is decomposed first by ICEEMDAN, the independent components are analyzed by robust ICA method, and the wavelet threshold denoising is carried out on the obtained independent components. Under the same wavelet conditions, decomposition layer number and threshold selection method, this method provides the possibility of further improving the SNR. The SNR improvement effect is related to the signal itself and the signal interception length. After the signal is segmented and then analyzed by the above methods, it is possible to further improve the SNR. The Prony analysis of the denoised signal can reduce the order of Prony analysis and obtain more accurate results.

We try to use this combined method to analyze the torsional vibration of shafting and the subsynchronous oscillation harmonics of power system. We get more accurate resonant frequency and related parameters. This combined method is suitable for analyzing signals with low SNR and complex inter-harmonics. In our experiments, we found that the segmentation of time-varying signals can get better results, but with the change of signals, we will make further efforts to study how to segment more reasonably.

Footnotes

Acknowledgments

This research was supported by Department of Electrical Engineering, College of Engineering, Nanjing Agricultural University, Nanjing, China, SCHOOL of Electrical Engineering, Southeast University, Nanjing, China, School of Automation, Nanjing Institute of Technology, Nanjing, China.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

Data Availability Statement

The data used to support the findings of this study are available from the corresponding author upon request.

ORCID iD

Qingjie Zhang

References

Bandt

Pompe

(2002) Permutation entropy: a natural complexity measure for time series. Physical Review Letters 88(17): 174102.

Daubechies

Heil

(1998) Ten lectures on wavelets. Computers in Physics 6(3): 1671.

Daubechies

(2002) The wavelet transform, time-frequency localization and signal analysis. IEEE Transactions on Information Theory. 36(5): 961–1005.

Huang

Shen

Long

, et al. (1998) The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proceedings Mathematical Physical & Engineering Sciences, 454(1971): 903–995.

Huang

(2008) A review on Hilbert‐transform: method and its applications to geophysical studies. Reviews of Geophysics 46(2): RG2006.

Janssens

Vlierberghe

D’Hondt

, et al. (2011) Zebra Tape Butt Joint Algorithm for Torsional Vibrations. New York, NY: Springer

Ouyang

Shuai

, et al. (2016) An indirect torsional vibration receptance measurement method for shaft structures. Journal of Sound and Vibration 372: 11–30

Marple

Carey

(1989) Digital spectral analysis with applications. The Journal of the Acoustical Society of America 86(5): 2043.

Nejad

Rizi

(2018) An Adaptive FECG Extraction and Analysis Method Using ICA, ICEEMDAN and Wavelet Shrinkage. Electrical engineering (ICEE), Iranian conference on, Mashhad, IR, 8–10 May 2018.

10.

Okamoto

Kurita

Sanchez-Gasca

, et al. (1996) Identification of equivalent linear power system models from electromagnetic transient time domain simulations using Prony's method. Decision and Control, 1996. Proceedings of the 35th IEEE, Kobe, Japan, 13 December 1996.

11.

Taradai

Deomidova

Zile

, et al. (2018) Results from investigations of torsional vibration in turbine set shaft systems. Thermal Engineering 65(1): 17–26.

12.

Terrien

Marque

Karlsson

(2011) Automatic detection of mode mixing in empirical mode decomposition using non-stationarity detection: application to selecting IMFs of interest and denoising. EURASIP Journal on Advances in Signal Processing 2011(1): 1–8.

13.

Torres

Colominas

Schlotthauer

(2011) A complete ensemble empirical mode decomposition with adaptive noise, IEEE International Conference on Acoustics, Prague, CR, 22 July 2011

14.

Xue

Howard

(2018) Torsional vibration signal analysis as a diagnostic tool for planetary gear fault detection. Mechanical Systems and Signal Processing 100: 706–728.

15.

Weaver

Healy

, et al. (1991) Filtering noise from images with wavelet transforms. Magnetic Resonance in Medicine 21(2): 288–295.

16.

Huang

(2009) Ensemble empirical mode decomposition: a noise-assisted data analysis method. Advances in Adaptive Data Analysis 01(01): 1–41.

17.

Yeh

Shieh

Huang

(2010) Complementary Ensemble empirical mode decomposition: a novel noise enhanced data analysis method. Advances in Adaptive Data Analysis. 2(2): 135–156.

18.

Yan

Gao

(2007) Approximate entropy as a diagnostic tool for machine health monitoring. Mechanical Systems and Signal Processing 21(2): 824–839.

19.

Zarzoso

Comon

(2010) Robust independent component analysis by iterative maximization of the kurtosis contrast with algebraic optimal step size. IEEE Transactions on Neural Networks 21(2):248–261.

20.

Zheng

Cheng

Yang

(2013) Modified EEMD algorithm and its applications. Journal of Vibration and Shock. 32(21):21–26.