A Review of Signal Analysis Methods and Their Applications in the Reversible Pump Turbine

Abstract

This chapter summarizes the signal analysis methods and their applications in the field of reversible pump turbine, including time-domain analysis, time-frequency analysis, mode decomposition and de-noising. The time-domain analysis methods include the time-domain statistical parameters, correlation analysis and chaos evaluation. The time-frequency analysis methods include the short-time Fourier transform, Hilbert-Huang transform and VMD-Hilbert transform. Furthermore, the signal mode decomposition and signal de-noising will be introduced together with the related evaluation indexes. The applications of the aforementioned methods are demonstrated based on both the simulated ideal signals and the measured signals from the prototype reversible pump turbines.

1. Introduction

There are a large number of rotating machineries being intensively employed in the modern industry, such as hydro-turbines, pumps, steam turbines, compressors, and generators. They are widely employed in the power generation, heating, refrigeration, and chemical industries. The rotating machineries often involve extreme conditions such as extremely low load, high temperature, high pressure, high speed, and overload. Hence, the malfunctions may occur after the long-term operation with possibly catastrophic consequences. Therefore, it is very important to understand and monitor the operational states of rotating machineries.

With the developments of data analysis technology, signal analysis methods become gradually mature and have been widely applied in the field of condition monitoring and fault diagnosis of various kinds of rotating machineries. The signal analysis could extract useful information from the original signal, and judge the operational states of the equipment based on the obtained information. Signal analysis is of great significance to the rotating machineries both for condition monitoring and fault diagnosis.

A reversible pump turbine is the core part of the pumped hydro energy storage power station. It can switch between the pumping and the generating mode according to the actual demand. When the reversible pump turbine deviates from the design working condition, it is easy to produce serious pressure pulsation and vibration due to the influence of rotor-stator interaction in the vaneless space and the vortex rope in the draft tube, which will affect the normal operation of the pump turbine and the safety of the whole power station [1, 2, 3]. Therefore, to monitor the operational states of reversible pump turbine and ensure its safety and efficient operation, it is very necessary to employ the signal analysis methods extensively.

This chapter will review several typical signal analysis methods and introduce their applications in the field of the reversible pump turbine in detail through specific cases. This chapter will be divided into the following four parts. The first part will introduce the time-domain analysis methods of the signals. The second part will introduce the time-frequency analysis methods of the signals. The third part will introduce the signal decomposition and the signal de-noising. The fourth part will demonstrate the applications of the introduced signal analysis methods in the field of the reversible pump turbines with the aid of on-site measured signals.

2. Time-domain analysis methods

Signals can convey information by expressing the relationship between time and other physical quantities. The time-domain analysis methods of signals refer to a series of processing such as amplification, filtering, statistical feature calculation and correlation analysis in the time domain. Through time-domain analysis methods, the characteristic parameters reflecting the operational state of the mechanical equipment can be extracted from the signal, which could be further employed for the purpose of the evaluation of the operational state and fault diagnosis of the equipment. Time-domain analysis methods could also analyze the auto-correlation (AC) and cross-correlation characteristics of the signal together with the degree of chaos. This section will be divided into three parts including the time domain statistical parameters, the correlation analysis and the chaos evaluation.

2.1 Time-domain statistical parameters

The classical time-domain statistical parameters include the peak value, peak-to-peak value, mean value, mean square value, root mean square, variance and standard deviation. The meanings of these statistical parameters are explained in detail below. Here, the $x\left(t\right)$ is the signal varying with the time. The $t_1$ and $t_2$ are the starting and the ending time of the signal in the time domain.

Peak value

Peak value refers to the maximum amplitude of a signal.

$$x_p=max\left|x\left(t\right)\right|$$E1

Peak-to-peak value

Peak-to-peak value refers to the difference between the maximum value and the minimum value of a signal within a given period. In the practical application, to eliminate the influences of distortion and noise in the signal, the peak-to-peak value is usually calculated by the confidence interval method with employing the confidence coefficient (e.g., 97 or 95%). Then, the peak-to-peak value could be calculated by using the upper and the lower limits of the obtained signal.

Mean value

The mean value is the averaged value of the whole signal.

$$\overline{x}=\frac{1}{t_2-t_1}{\int }_{t_1}^{t_2}x\left(t\right)\mathit{dt}$$E2

Mean square value

Mean square value refers to the mean of the square of the signal, indicating the strength and average power of the signal.

$$x_{\mathit{ms}}=\frac{1}{t_2-t_1}{\int }_{t_1}^{t_2}{x}^2\left(t\right)\mathit{dt}$$E3

Root mean square

Root mean square refers to the arithmetic square root of the mean square value.

$$x_{\mathit{rms}}=\sqrt{\frac{1}{t_2-t_1}{\int }_{t_1}^{t_2}{x}^2\left(t\right)\mathit{dt}}$$E4

Variance

Variance is the average square value of the difference between the amplitude of the original signal and the mean value. It describes the fluctuation range of the signal and represents the strength of the fluctuating component in the signal.

$$s^2=\frac{1}{t_2-t_1}{\int }_{t_1}^{t_2}{\left(x-\overline{x}\right)}^2\mathit{dt}$$E5

Standard deviation

Standard deviation is the arithmetic square root of variance.

$$s=\sqrt{\frac{1}{t_2-t_1}{\int }_{t_1}^{t_2}{\left(x-\overline{x}\right)}^2\mathit{dt}}$$E6

2.2 Correlation analysis

When there are multiple signals, sometimes it is necessary to study the relationship between them. Correlation analysis (including auto-correlation function and cross-correlation analysis) could show the similarity or dependence between signals.

Auto-correlation function

The auto-correlation (AC) function describes the correlation of a signal between its values at a certain time and after a certain time delay $\tau$. It can quantitatively describe the degree of similarity between the time-domain waveform of a signal after shifting $\tau$ on the time and the original time-domain waveform of a signal.

$$R_{\mathit{xx}}\left(\tau \right)={\int }_{-\infty }^{\infty }x\left(t\right)x\left(t+\tau \right)\mathit{dt}={\int }_{-\infty }^{\infty }x\left(t-\tau \right)x\left(t\right)\mathit{dt}$$E7

$R_{\mathit{xx}}\left(\tau \right)$ is an even function. When $\tau =0$, $R_{\mathit{xx}}\left(\tau \right)$ reaches its maximum value and is equal to the energy of the signal.

Cross-correlation function

The cross-correlation function describes the degree of similarity of the time-domain waveforms of two signals $x\left(t\right)$ and $y\left(t\right)$ at $t$ and $t+\tau$.

$$R_{\mathit{xy}}\left(\tau \right)={\int }_{-\infty }^{\infty }x\left(t\right)y\left(t+\tau \right)\mathit{dt}={\int }_{-\infty }^{\infty }x\left(t-\tau \right)y\left(t\right)\mathit{dt}$$E8

The higher the $R_{\mathit{xy}}\left(t\right)$ is, the higher the similarity between $x\left(t\right)$ and $y\left(t+\tau \right)$ is. It should be noted that $R_{\mathit{xy}}\left(\tau \right)$ is not equal to $R_{\mathit{yx}}\left(\tau \right)$, and the relationship between them is as follows:

$$R_{\mathit{xy}}\left(\tau \right)=R_{\mathit{yx}}\left(-\tau \right)$$E9

Normalization of the correlation function

In the practical application of signal correlation analysis, to compare the magnitude of the correlation, we often use the normalized form of the correlation function. The normalized form of the auto-correlation function is defined as follows:

$${\rho }_{\mathit{xx}}\left(\tau \right)=\frac{R_{\mathit{xx}}\left(\tau \right)}{\sqrt{R_{\mathit{xx}}\left(0\right)}\sqrt{R_{\mathit{xx}}\left(0\right)}}=\frac{R_{\mathit{xx}}\left(\tau \right)}{R_{\mathit{xx}}\left(0\right)}$$E10

The normalized form of the cross-correlation function is defined as follows:

$${\rho }_{\mathit{xy}}\left(\tau \right)=\frac{R_{\mathit{xy}}\left(\tau \right)}{\sqrt{R_{\mathit{xx}}\left(0\right)}\sqrt{R_{\mathit{yy}}\left(0\right)}}$$E11

where $\left|{\rho }_{\mathit{xx}}\left(\tau \right)\right|\le 1$,$\left|{\rho }_{\mathit{xy}}\left(\tau \right)\right|\le 1$. When $\rho \left(\tau \right)=1$, it means that the correlation between the two signals reaches the maximum; when $\rho \left(\tau \right)=0$, it means that two signals are not correlated; when $\rho \left(\tau \right)=-1$, it means that two signals are linearly negatively correlated.

2.3 Measurement of chaos

Chaos is an inherent characteristic of a nonlinear dynamic system. When a signal is generated by a nonlinear system, it may show the characteristic of chaos. At this time, the signal waveform is very irregular and is very similar to the characteristic of random noise. Entropy can be adopted to measure the chaos of a nonlinear dynamic system. The more chaotic the system, the higher its entropy. Therefore, we can use entropy to evaluate the degree of chaos in a signal. As a typical method, the permutation entropy for evaluating the degree of chaos in a signal is introduced below [4].

Permutation entropy

Permutation entropy is an index to measure the complexity of the time series [4]. The more regular the time series is, the smaller the permutation entropy is. The corresponding calculation process is given as follows [4].

Considering a time series $x=\left\{x_i,i=1,2,\dots ,N\right\}$, it can be reconstructed into the following form [4]:

$$\left\{\begin{array}{l}{X}_1=\left\{x_1,x_{1+\tau },\dots ,x_{1+\left(\mathit{em}-1\right)\tau }\right\}\\ ⋮\\ X_i=\left\{x_i,x_{i+\tau },\dots ,x_{i+\left(\mathit{em}-1\right)\tau }\right\}\\ ⋮\\ X_{N-\left(\mathit{em}-1\right)\tau }=\left\{x_{N-\left(\mathit{em}-1\right)\tau },x_{N-\left(\mathit{em}-2\right)\tau },\dots ,x_N\right\}\end{array}\right.$$E12

Here, the “N” is the total number of the data points within the whole signal. The “$\tau$” refers to the delay time and the “em” refers to the embedded dimension.

Now, each $X_i$ could be sorted in the ascending order internally [4]:

$$x_{i+\left(j_1-1\right)\tau }\le x_{i+\left(j_2-1\right)\tau }\le \dots \le x_{i+\left(j_{\mathit{em}}-1\right)\tau }$$E13

Here, the “$j_1,j_2,\dots ,j_{\mathit{em}}$” in the subscripts refer to the indexes for different columns. Thus, a sequence of symbols can be obtained for any vector $X_i$ [4]:

$$S\left(h\right)=\left[j_1,j_2,\dots ,j_{\mathit{em}}\right]$$E14

There are up to “$\mathit{em}!$” kinds of such index sequences. The occurrence probability of each index sequence is calculated and expressed as $P_1,P_2,\dots ,P_r$, respectively. Here, $h=1,2,\dots ,r$ and $r\le \mathit{em}!$. Then, the permutation entropy of the original signal can be expressed as follows [4]:

$$H_p\left(\mathit{em}\right)=-\sum _{h=1}^r{P}_{h}ln{P}_h$$E15

When $P_h=1/\mathit{em}!$, the above formula can reach the maximum value $ln\left(\mathit{em}!\right)$. Therefore, the normalized form of Eq. (15) can be obtained as follows [4]:

$$H\left(\mathit{em}\right)=\frac{H_P\left(\mathit{em}\right)}{ln\left(\mathit{em}!\right)}$$E16

The variation range of the normalized permutation entropy is between 0 and 1. In the practical application of calculating the permutation entropy of a signal, the whole data points within the signal can be intercepted by a window with a certain length. And the permutation entropy of the data points in each window is calculated by moving the window one data at a time. In the present chapter, the average permutation entropy of all the windows will be taken as the permutation entropy of the whole signal.

2.4 Examples of a signal time-domain analysis

In this section, the pressure pulsation signal (the pressure pulsation signal 1) in the vaneless space of a prototype pump turbine under the dimensionless load condition (the ratio of the actual operational load to the rated load) P^*=52% is taken in the analysis to demonstrate the applications of the time-domain analysis method introduced in Sections 2.1–2.3. The length of the signal is 1 s and the sampling frequency is 2000 Hz. In this section, p^* refers to the dimensionless pressure pulsation using the water head (326.01 m). Figure 1 shows the dimensionless pressure pulsation 1 in the time domain.

Figure 1.

The time-domain diagram of the dimensionless pressure pulsation signal 1.

Time-domain statistical parameters

The statistical parameters of the dimensionless pressure pulsation signal 1 in the time domain are shown in the Table 1 as follows:

Statistical parameter	Value
Peak value	0.9492
Peak-to-peak value (confidence coefficient 97%)	0.0911
Peak-to-peak value (confidence coefficient 95%)	0.0883
Mean	0.9081
Mean square value	0.8252
Root mean square	0.9084
Variance	0.0006
Standard deviation	0.0244

Table 1.

The time-domain statistical parameters of the dimensionless pressure pulsation signal 1.

Auto-correlation analysis

Figures 2 and 3 show the normalized AC functions of the dimensionless pressure pulsation signal 1 and a random noise signal. When $\tau =0$, the AC function of the dimensionless pressure pulsation signal 1 in the Figure 2 reaches its maximum value, but it gradually decreases to near 0 on both sides of $\tau =0$. For the random noise signal, its AC function reaches the maximum also at $\tau =0$. When $\tau \ne 0$, due to the randomness of noise, its AC function decreases rapidly. Therefore, the normalized AC function of the random noise in Figure 3 is close to 1 at $\tau =0$ and approximately equal to 0 on both sides of $\tau =0$.

Figure 2.

The normalized AC function of the dimensionless pressure pulsation signal 1.

Figure 3.

The normalized AC function of the random noise.

Cross-correlation analysis

In the experiment, another dimensionless pressure pulsation signal (the dimensionless pressure pulsation signal 2) was measured in the vaneless space at the same operational condition (P^*=52%) during the different time period. Figure 4 shows the time-domain diagram of the dimensionless pressure pulsation signal 2.

Figure 4.

The time-domain diagram of the dimensionless pressure pulsation signal 2.

When the delay time $\tau =0$, the normalized cross-correlation function of the dimensionless pressure pulsation signals 1 and 2 is 0.98, indicating that there is a strong correlation between them. On the contrary, when the delay time $\tau =0$, the normalized cross-correlation function of the dimensionless pressure pulsation signal 1 and random noise signal is 0.0088, indicating that there is almost no correlation between them.

Permutation entropy analysis

In the present chapter, the embedded dimension, delay time and window length required in the analysis procedure of the permutation entropy are 6, 3 and 128, respectively. The analysis results of the permutation entropy of the dimensionless pressure pulsation signal 1, the dimensionless pressure pulsation signal 2 and the random noise signal are as follows:

It can be seen from the Table 2 that the analysis results of the permutation entropy of the dimensionless pressure pulsation signals 1 and 2 are smaller than that of the random noise, because the dimensionless pressure pulsation signals contain certain characteristic frequencies with regularity. The permutation entropy of the random noise is large, which indicates its high degree of chaos.

Signal	Permutation entropy
Dimensionless pressure pulsation signal 1	0.4254
Dimensionless pressure pulsation signal 2	0.4271
Random noise	0.5666

Table 2.

Analysis results of permutation entropy of different types of signals.

3. Time-frequency analysis methods

The aim of the time-frequency analysis is to study the time-varying signals, which can reflect the relationship of the variations of the frequency and amplitude with time within the signal. Generally speaking, the typical time-frequency analysis methods include the short-time Fourier transform [5], the Hilbert-Huang transform [6] and the VMD-Hilbert transform with a detailed introduction in this section.

3.1 The traditional Fourier transform

Before introducing the time-frequency analysis methods, the traditional Fourier transform is:

$$F\left(\omega \right)={\int }_{-\infty }^{\infty }x\left(t\right)e^{-\mathit{j\omega t}}\mathit{dt}$$E17

where $x\left(t\right)$ is the original signal and $e^{-\mathit{j\omega t}}$ is the basis function of the Fourier transform. The traditional Fourier transform can accurately reflect the frequency components contained in the signal, but it cannot reflect the time information of these frequency components.

3.2 The time-frequency analysis methods

Short-time Fourier transform

$${\text{STFT}}_x\left(\tau ,\omega \right)={\int }_{-\infty }^{\infty }x\left(t\right)g\left(\tau -t\right)e^{-\mathit{j\omega t}}\mathit{dt}$$E18

The formula of the short-time Fourier Transform (STFT) is presented above [5]. Compared with the traditional Fourier transform, the window function is added to the formula. STFT needs to select a window function, which can be moved continuously to identify the frequencies at different times. The window function will affect the resolution of the spectrum. The smaller the window size, the higher the time resolution of the spectrum. Normally, the low frequency band requires high frequency resolution and low time resolution. And, the high frequency band requires low frequency resolution and high time resolution. However, after the window function is selected, the resolution of the entire spectrum is fixed and cannot be adjusted. Therefore, STFT is suitable for analyzing the segmented stationary signals or approximately stationary signals.

Hilbert-Huang transform

Hilbert-Huang transform (HHT) consists of two parts including the empirical mode decomposition (EMD) and the Hilbert spectrum analysis (HSA). These two parts will be introduced in detail below.

The EMD method can decompose the original signal into several intrinsic mode functions (IMFs) and the residual component $r\left(t\right)$. The specific steps are given as follow.

Firstly, calculating all the extreme points of the signal $x\left(t\right)$, and using the cubic spline function to interpolate and fit the maximum and minimum points respectively to obtain the upper envelope $a\left(t\right)$ and lower envelopes $b\left(t\right)$ of the signal;

Secondly, calculating the mean curve $c\left(t\right)$ of the upper and lower envelopes and the difference $h\left(t\right)$ between the signal $x\left(t\right)$ and $c\left(t\right)$ [6]:

$$c\left(t\right)=\frac{a\left(t\right)+b\left(t\right)}{2}$$E19

$$h\left(t\right)=x\left(t\right)-c\left(t\right)$$E20

Thirdly, the screening stop condition is given as [6]:

$$\sum _{t=0}^T\frac{{\left|h_{g-1}\left(t\right)-h_g\left(t\right)\right|}^2}{h_{g-1}^2\left(t\right)}\le \epsilon$$E21

Here, the “$g$“ refers to the number of screenings. The “$T$” refers to the signal time length and the “$\epsilon$” refers to a preset value. If $h\left(t\right)$ does not meet the condition, replacing$x\left(t\right)$ with $h\left(t\right)$ and repeating the above steps. And if $h\left(t\right)$ meets the condition, then $h\left(t\right)$ will be the first IMF to be extracted.

Fourthly, repeating the above steps to obtain all the remaining IMF until the residual component $r\left(t\right)$ becomes a constant value sequence with a frequency approximately 0.

Although, the EMD can realize signal decomposition, it has the problems of endpoint effect and mode mixing. Endpoint effect means that analysis results at the signal endpoints will produce errors, which will affect the internal data. Mode mixing means that the EMD method cannot effectively separate different mode components based on the time characteristic scales. These problems can lead to serious performance degradation of EMD methods.

HSA performs a Hilbert transform on each IMF obtained by EMD. The specific steps are given as follows:

Firstly, performing the Hilbert transform on each IMF $h_i\left(t\right)$ [6]:

$$H\left[h_i\left(t\right)\right]=h_i\left(t\right)\ast \frac{1}{\mathit{\pi t}}=\frac{1}{\pi }{\int }_{-\infty }^{\infty }\frac{h_i\left(\tau \right)}{t-\tau }\mathit{d\tau }$$E22

Here, i represents the i-th IMF component; “*” represents the convolution operation.

Secondly, using the above formula to construct the analytical signal $d_i\left(t\right)$ of each IMF $h_i\left(t\right)$ [6]:

$$d_i\left(t\right)=h_i\left(t\right)+\mathit{jH}\left[h_i\left(t\right)\right]=a_i\left(t\right)e^{j{\phi }_i\left(t\right)}$$E23

Here [6],

$$a_i\left(t\right)=\sqrt{h_i^2\left(t\right)+H^2\left[h_i\left(t\right)\right]}$$E24

$${\phi }_i\left(t\right)=arctan\left\{\frac{H\left[h_i\left(t\right)\right]}{h_i\left(t\right)}\right\}$$E25

Thirdly, the instantaneous frequency of each IMF can be obtained from the following formula [6]:

$$w_i\left(t\right)=\frac{d{\phi }_i\left(t\right)}{\mathit{dt}}$$E26

After expressing the results on the time-frequency plane, the Hilbert amplitude spectrum of each IMF $h_i\left(t\right)$ can be obtained as follows [6]:

$$H_i\left(w,t\right)=\left\{\begin{array}{l}0,w\ne w_i\left(t\right)\\ a_i\left(t\right),w=w_i\left(t\right)\end{array}\right.$$E27

Fourthly, by presenting the amplitude spectrum of all IMFs in a spectrum, the Hilbert spectrum of the original signal can be obtained as follows [6]:

$$H\left(w,t\right)=Re\sum _{i=1}^n{a}_i\left(t\right)e^{j\int w_i\left(t\right)\mathit{dt}}$$E28

The advantage of HHT is that it is not restricted by the linear and stationary characteristics of the signal and can analyze non-linear and non-stationary signals. At the same time, the HHT is self-adaptive through using the EMD to eliminate the choice of basis function. In addition, HHT is not restricted by Heisenberg’s uncertainty principle.

VMD-Hilbert transform

The principle of VMD-Hilbert transform is similar to HHT. It also uses the signal decomposition method to obtain the IMF component of the signal, and then obtains the Hilbert spectrum of the signal through HAS. The difference between these two is that the variational mode decomposition (VMD) [7] is adopted as the signal decomposition method.

The core idea of the VMD method is to construct and solve the variational optimization problem in the frequency domain. The specific steps include the following ones.

Firstly, decomposing the original signal into K IMFs with respective central angular frequencies and limited bandwidths [7]:

$${\mathrm{IMF}}_k\left(t\right)=A_k\left(t\right)cos\left[{\phi }_k\left(t\right)\right]$$E29

Here, ${\phi }_k\left(t\right)$ refers to phase of the signal and $A_k\left(t\right)$ refers to instantaneous amplitude.

Secondly, the constraints are added to ensure that the sum of the estimated bandwidths of each IMF is minimum and the sum of K IMFs is equal to the original signal [7].

$$\left\{\begin{array}{l}\underset{\left\{{\mathrm{IMF}}_k\right\},\left\{{\omega }_k\right\}}{min}\left\{\sum _{k=1}^K{‖{\partial }_t\left[\left(\delta \left(t\right)+\frac{j}{\mathit{\pi t}}\right)\ast {\mathrm{IMF}}_k\left(t\right)\right]e^{-j{\omega }_{k}t}‖}_2^2\right\}\\ s.t.\sum _{k=1}^K{\mathrm{IMF}}_k\left(t\right)=x\left(t\right)\end{array}\right.$$E30

Here, $\delta \left(t\right)$ refers to the Dirac function. ${\omega }_k$ is the central angle frequency of the k-th IMF_k. “*” refers to the convolution operation.

Thirdly, the Lagrange multiplication operator is introduced to transform the constrained variational problem into an unconstrained variational problem. And an augmented Lagrange expression is obtained as follows [7]:

$$\begin{array}{l}L\left(\left\{{\mathrm{IMF}}_k\right\},\left\{{\omega }_k\right\},\lambda \right)=\alpha \sum _{k=1}^K{‖{\partial }_t\left[\left(\delta \left(t\right)+\frac{j}{\mathit{\pi t}}\right)\ast {\mathrm{IMF}}_k\left(t\right)\right]e^{-j{\omega }_{k}t}‖}_2^2+\\ {‖\mathrm{IMF}\left(t\right)-\sum _{k=1}^K{\mathrm{IMF}}_k\left(t\right)‖}_2^2+⟨\lambda \left(t\right),\mathrm{IMF}\left(t\right)-\sum _{k=1}^K{\mathrm{IMF}}_k\left(t\right)⟩\end{array}$$E31

Here, $\alpha$ refers to the quadratic penalty factor.

Fourthly, the alternating direction multiplier method is used to iteratively solve IMF_k, ${\omega }_k$ and $\lambda$ to find the saddle point of the unconstrained variational problem. The update formulas of $\mathrm{I}\widehat{\mathrm{M}}{{\mathrm{F}}_k}_{n+1}^{}\left(\omega \right)$, ${\omega }_k^{n+1}$ and ${\widehat{\lambda }}^{n+1}\left(\omega \right)$ obtained by Parseval Fourier equidistant transformation could be expressed as [7]:

$$\mathrm{I}\widehat{\mathrm{M}}{{\mathrm{F}}_k}_{n+1}^{}\left(\omega \right)=\frac{\widehat{x}\left(\omega \right)-{\sum }_{i=1}^{k-1}\mathrm{I}\widehat{\mathrm{M}}{\mathrm{F}}_i^{n+1}\left(\omega \right)-{\sum }_{i=k+1}^K\mathrm{I}\widehat{\mathrm{M}}{\mathrm{F}}_i^n\left(\omega \right)+\frac{{\widehat{\lambda }}^n\left(\omega \right)}{2}}{1+2\alpha {\left(\omega -{\omega }_k^n\right)}^2}$$E32

$${{\omega }_k}_{n+1}^{}=\frac{{\int }_0^{\infty }\omega {\left|\mathrm{I}\widehat{\mathrm{M}}{\mathrm{F}}_k\left(\omega \right)\right|}^2\mathit{d\omega }}{{\int }_0^{\infty }{\left|\mathrm{I}\widehat{\mathrm{M}}{\mathrm{F}}_k\left(\omega \right)\right|}^2\mathit{d\omega }}$$E33

$${\widehat{\lambda }}^{n+1}\left(\omega \right)={\widehat{\lambda }}^n\left(\omega \right)+\gamma \left(\widehat{x}\left(\omega \right)-\sum _{k=1}^K\mathrm{I}\widehat{\mathrm{M}}{\mathrm{F}}_k^{n+1}\left(\omega \right)\right)$$E34

Here, “^” refers to the Fourier transform operation. The “$\gamma$” refers to the noise tolerance. And the “$n$” refers to the number of iterations in the calculation process.

When the calculation result meets the given solution accuracy $\epsilon$, the iterative operation is stopped and K IMFs are outputted [7].

$$\sum _{k=1}^K\left(\frac{{‖\mathrm{I}\widehat{\mathrm{M}}{\mathrm{F}}_k^{n+1}-\mathrm{I}\widehat{\mathrm{M}}{\mathrm{F}}_k^n‖}_2^2}{{‖\mathrm{I}\widehat{\mathrm{M}}{\mathrm{F}}_k^n‖}_2^2}\right)<\epsilon$$E35

Fifthly, the $\mathrm{I}\widehat{\mathrm{M}}{{\mathrm{F}}_k}_{n+1}^{}\left(\omega \right)$ is transformed from the frequency domain to the IMF in the time domain through the operation of the inverse Fourier transform.

Compared with the EMD, the VMD overcomes the problem of mode mixing and has a more solid mathematical theoretical foundation.

In the present chapter, the “K” is 7, the “$\alpha$” is 2000, the “$\gamma$” is 0 with “$\epsilon$” is 10⁻⁷.

3.3 Examples of signal time-frequency analysis

This section uses the vibration signal at the top cover of a prototype reversible pump turbine during the start-up transient process of generating mode as the analysis object to demonstrate the applications of the time-frequency analysis methods introduced in Sections 3.1 and 3.2. The length of the signal is 2 s and the sampling frequency is 1000 Hz. Figure 5 shows the time-domain diagram of the vibration signal at the top cover.

Figure 5.

The time-domain diagram of the vibration signal at the top cover.

Figure 6 shows the result of the fast Fourier transform (FFT) of the vibration signal. It can be seen from the Figure 6 that the most obvious frequency component in the vibration signal is 40 Hz with an amplitude of about 0.46 mm. In addition, there is also a relatively obvious peak near 80 Hz with an amplitude of about 0.17 mm. However, Figure 6 cannot provide information about the amplitude variation of different frequency components with time, which is a shortcoming of the traditional Fourier transform.

Figure 6.

The FFT spectrum of the vibration signal.

Figures 7 and 8 are the STFT results of the vibration signal with the window length of 64 and 256 respectively. From these two figures, not only the amplitudes of different frequency components in the vibration signals but also the variations of the amplitude of each frequency varying with the time can be seen. However, the STFT results are significantly affected by the window size. By comparing Figures 7 and 8, it can be concluded that there exist obvious amplitudes at and around 40 Hz in the Figure 7, while it almost has obvious amplitudes only at 40 Hz in the Figure 8, which is caused by the difference in the window size between them. In the actual application of the STFT, it is difficult to choose the appropriate window size, which is the main problem of this method.

Figure 7.

The STFT result of the vibration signal with the window length of 64 and the window function of Blackman.

Figure 8.

The STFT result of the vibration signal with the window length of 256 and the window function of Blackman.

Different with the fixed basis functions of traditional Fourier transform and STFT, HHT and VMD-Hilbert transform can obtain basis functions through adaptive mode decomposition. In addition, HHT and VMD-Hilbert transform do not need to select the window functions and they can achieve high resolutions both in time and frequency domains. Figures 9 and 10 are the results of the HHT and VMD-Hilbert transform of the signal, respectively. It can be found that the analysis result in the Figure 9 is unclear because the dominant frequencies of the signal cannot be clearly observed, while the analysis result in the Figure 10 is very clear. This can be attributed to that the VMD overcomes the mode mixing problem in the EMD adopted in the HHT.

Figure 9.

The analysis result of HHT of the vibration signal.

Figure 10.

The analysis result of VMD-Hilbert transform of the vibration signal.

4. Signal mode decomposition and signal de-noising

Mode refers to the natural vibration characteristics of the mechanical structure, which is composed of natural frequency, damping ratio and mode shape. In practical engineering, the vibration of rotating machinery is often composed of different complex modes. The process of decomposing the vibration signal into sub-signals with different frequency bands on the basis of different modes is called the mode decomposition of the signal. Through the mode decomposition, the modes containing noise can be found and will be further deleted through the signal de-noising. In our previous research [8], similar methods with those shown in this section have been adopted to analyze the vibration signal from the rotating machinery of nuclear power. However, the specific frequency components between the signals between the rotating machinery of nuclear power and the reversible pump turbine are a little different.

4.1 Mode decomposition of simulated swing signal

Typical signal analysis methods, such as EMD and VMD, have been fully described in Section 3.2. In this section, a simulated swing signal at the turbine guide bearing of a prototype pump turbine will be used to demonstrate the signal de-noising effects based on these two methods. The simulated rotational speed of the unit is 500 rpm. The characteristic frequencies contained in the simulated signal are 1 time, 2 times and 3 times, the rotational frequency of the impeller with the corresponding amplitudes of 10, 2 and 1 μm, respectively. The signal is dimensionless by dividing the corresponding amplitude of the rotational frequency of the impeller (10 μm). The length of the signal is 2 s and the sampling frequency is 1000 Hz. The white noise with the signal-noise ratio (SNR) of 20 dB is added to the simulated swing signal. Figure 11 shows the time-domain diagrams of the simulated swing signals before and after adding the white noise. Figures 12 and 13 are the EMD and VMD analysis results of the simulated swing signal, respectively.

Figure 11.

The time-domain diagrams of simulated swing signals. (a) The simulated swing signal without noise. (b) The added noise with the SNR of 20 dB. (c) The simulated swing signal with noise added.

Figure 12.

The EMD analysis results of the simulated swing signal.

Figure 13.

The VMD analysis results of the simulated swing signal.

4.2 Results of signal de-noising

Through mode decomposition, several IMFs are obtained to further constitute the given signal with low noise. The IMFs contain different frequency bands. By removing the IMFs dominated by the random noises and retaining the IMFs dominated by the useful signals, the signal de-noising will be proceeded. Generally speaking, the cross-correlation coefficient between the noise-dominated IMF and the original signal is small, while the cross-correlation coefficient between the IMF component dominated by the useful signal and the original signal is relatively large. Based on this principle, the noise-dominated IMFs can be identified according to the cross-correlation coefficients.

Table 3 shows the cross-correlation coefficients of the IMFs obtained by EMD and VMD and the simulated swing signal. Here, the threshold of the cross-correlation coefficient is set to 0.12 [8]. And the IMFs with the cross-correlation coefficients larger than this threshold are retained.

Mode component	EMD	VMD
IMF1	0.0757	0.9732
IMF2	0.0537	0.2378
IMF3	0.7777	0.0453
IMF4	0.3762	0.0430
IMF5	0.0067	0.0407
IMF6	0.0242	0.0412
IMF7	0.0349	0.0400
residual	0.0103	—

Table 3.

Cross-correlation coefficients of IMFs based on EMD and VMD and the simulated swing signal.

Based on the analysis results in Table 3, the IMF₃ and IMF₄ in the EMD analysis results and the IMF₁ and IMF₂ in the VMD analysis results are selected to reconstruct the de-noised swing signals, respectively. Figure 14 shows the de-noising results of the simulated swing signal based on EMD and VMD, respectively.

Figure 14.

The de-noising results of the simulated swing signal based on EMD and VMD. (a) The simulated swing signal without noise. (b) The de-noised signal is based on EMD. (c) The de-noised signal is based on VMD.

Figure 14 shows that both the de-noised signals based on EMD and VMD are relatively smooth oscillation curves. Compared with the simulated swing signal without noise, the de-noised signal based on EMD has some abnormal spikes or irregularities as indicated by the red circles in the Figure 14b, while the time-domain waveform of the de-noised signal based on VMD is much closer to the simulated swing signal without noise, which indicates that the signal de-noising based on VMD has a better effect than that based on EMD.

4.3 Evaluation indexes of signal de-noising

The evaluation indexes of signal de-noising can help us quantitatively analyze the de-noising effects. Some evaluation indexes will be introduced below. Here $x\left(t\right)$ refers to original simulated swing signal without noise; “i” refers to the i-th data point of signal $x\left(t\right)$ which contains M data points in total; ${\tilde{x}}_i$refers to de-noised swing signal.

Signal-noise ratio (SNR)

$$\mathit{SNR}=10lg\left(\frac{{\sum }_{i=1}^M{x}_i^2}{{\sum }_{i=1}^M{\left(x_i-{\tilde{x}}_i\right)}^2}\right)$$E36

The larger the SNR, the better the de-noising effect.

Correlation coefficient (R)

$$R=\frac{{\sum }_{i=1}^M{x}_i{\tilde{x}}_i}{\sqrt{{\sum }_{i=1}^M{x}_i^2{\sum }_{i=1}^M{\tilde{x}}_i^2}}$$E37

The larger the R, the better the de-noising effect.

Root-mean-square error (RMSE)

$$\mathit{\text{RMSE}}=\sqrt{\frac{1}{M}{\sum }_{i=1}^M{\left(x_i-{\tilde{x}}_i\right)}^2}$$E38

The smaller the RMSE, the better the de-noising effect.

The evaluation indexes of signal de-noising effects based on EMD and VMD are calculated respectively and the results are shown in Table 4.

Evaluation index	EMD	VMD
SNR	12.1328	29.5778
R	0.969	0.9995
RMSE	0.1856	0.024

Table 4.

The evaluation indexes of de-noising for the simulated swing signal based on EMD and VMD.

It can be seen from Table 4 that compared with the results of signal de-noising based on EMD, the results of signal de-noising based on VMD have a larger SNR, a larger R, and a smaller RMSE, which further shows that the signal de-noising based on VMD has a better effect than that based on EMD.

4.4 Decomposition and de-noising of measured swing signal

In Sections 4.1–4.3, the simulated swing signal is adopted to demonstrate the signal decomposition and the signal de-noising processes. In this section, the analysis results of mode decomposition and de-noising of measured swing signal of a prototype pump turbine are presented. The swing signal at the turbine guide bearing of a pump turbine was measured at P^*=15%. The rated rotational speed of the unit is 500 rpm. The length of the signal is 2 s and the sampling frequency is 1000 Hz. Figure 15 shows the time-domain diagram of the measured swing signal.

Figure 15.

The time-domain diagram of the measured swing signal.

The EMD and VMD analysis results of the measured swing signal are shown in Figures 16 and 17.

Figure 16.

The EMD analysis results of the measured swing signal.

Figure 17.

The VMD analysis results of the measured swing signal.

The cross-correlation coefficients of the IMFs obtained by EMD and VMD and the measured swing signal are shown in Table 5.

Mode component	EMD	VMD
IMF1	0.0635	0.9471
IMF2	0.0797	0.4739
IMF3	0.7677	0.1096
IMF4	0.2904	0.0323
IMF5	0.0500	0.0321
IMF6	0.0715	0.0200
IMF7	0.1455	0.0154
residual	0.0444	—

Table 5.

The cross-correlation coefficients of the IMFs based on EMD and VMD and the measured swing signal.

The threshold of the cross-correlation coefficient is also set to 0.12 [8]. Based on the analysis results in Table 5, the IMF₃, IMF₄, and IMF₇ in the EMD analysis results and the IMF₁ and IMF₂ in the VMD analysis result are selected to reconstruct the de-noised swing signals, respectively. Figure 18 shows the de-noising results of the measured swing signal based on EMD and VMD, respectively. As shown in Figure 18, the VMD shows a better de-noising effect by comparing the smoothness of the curves.

Figure 18.

The de-noising results of the measured swing signal is based on EMD and VMD. (a) The original measured swing signal. (b) The de-noised signal is based on the EMD. (c) The de-noised signal based on VMD.

5. Conclusions

In summary, signal analysis and processing can help us better understand the operational states of the reversible pump turbines. Through the time-domain analysis, the time-domain characteristic parameters such as peak value, peak-to-peak value and average value can be obtained. Through the time-frequency analysis, the time-frequency variation characteristics of signals can be grasped. Through the mode decomposition of the signal, various components within the signal can be distinguished and the signal de-noising can be further performed. Therefore, signal analysis plays a very important role for the condition monitoring and fault diagnosis of the reversible pump turbines.

Acknowledgments

This work was financially supported by the National Natural Science Foundation of China (Project Nos.: U1965106, 51976056 and 52076215).