Wavelet decomposition and nonlinear prediction of nonstationary vibration signals

Abstract

In order to overcome the shortcomings (such as the time–frequency localization and the nonstationary signal analysis ability) of the Fourier transform, time–frequency analysis has been carried out by wavelet packet decomposition and reconstruction according to the actual nonstationary vibration signal from a large equipment located in a large Steel Corporation in this article. The effect of wavelet decomposition on signal denoising and the selection of high-frequency weight coefficients for each layer on signal denoising were analyzed. The nonlinear prediction of the chaotic time series was made by global method, local method, weighted first-order local method, and maximum Lyapunov exponent prediction method correspondingly. It was found the multi-step prediction method is better than other prediction methods.

Keywords

Wavelet decomposition nonstationary signal nonlinear prediction vibration time series

Introduction

The most common and mature method for signal analysis is Fourier transform. Unfortunately, the Fourier transform is a global transformation, that is, it is either completely in time domain or completely in frequency domain. The signal cannot be analyzed in time–frequency localization simultaneously, and it cannot reveal information about nonstationary signals. But there are lots of nonsteady dynamic signals in engineering vibration, such as the speed down or up process of rotating machinery, friction, equipment foundation loosening, misalignment, crack, rotating stall, oil whirl, and oil whip fault in the operation process of machinery and equipment. The vibration signals are nonstationary; this feature can explain some fault characterization.^1–3 Because the statistical characteristics of nonstationary dynamic signals are related to time, the processing of nonstationary signals needs time–frequency analysis. And the whole picture and localization results of nonstationary signals in time and frequency domains should be obtained. Nonstationary signal processing method for engineering vibration generally has the following five types:^4–6 piecewise Fourier transform, rotational speed tracking analysis with Hanning window, short-time Fourier transform, Wigner–Ville distribution, wavelet analysis, and Hilbert–Huang transform.

As the vibration of the rolling mill is a very complex problem, it is usually difficult to find the vibration source and the vibration law. Until now, there is no perfect understanding of its vibration mechanism or a satisfactory theoretical system. In production, the vibration is often controlled by experience and test. In the thin slab continuous casting and rolling line of a steel company in China, the vibration of F2, F3, and F4 mills appears in varying degrees, especially in the F3 mill. Due to the severe vibration of the rolling mill, fatigue damage of equipment components has occurred, such as unstable roll transmission or even work roll neck fracture, electrical component damage, F3 top programmable controller damage, F3 main reducer high-speed gear shaft tooth surface peeling, strip steel, and roll surface vibration, which has brought great harm to production. The project team has carried out a large-scale finishing mill. Based on the comprehensive test and theoretical research, this article aims to explore the shape of roller vibration from the signal point of view and then the vibration source.

Wavelet packet decomposition and reconstruction

Given a basic function $ψ (t)$ , and $ψ (t) \in L^{2} (R)$ ( $L^{2} (R)$ is a square integrable real number space), its Fourier transform is $\hat{ψ} (ω)$ , and $\hat{ψ} (ω)$ meets the admissibility conditions as $C_{ψ} = \int_{R} | \hat{ψ} (ω) |^{2} / | ω | d ω < \infty$ . If there are no special instructions, the integration limits are from –∞ to +∞.^7–9 $ψ (t)$ is the base or mother wavelet. The base wavelet $ψ (t)$ meets $\hat{ψ} (ω = 0) = 0$ , that is, $\int ψ (t) d t = 0$ . The wavelet sequence is obtained from the mother wavelet dilation and translation, called sub-wavelet, that is

ψ_{a, b} (t) = \frac{1}{\sqrt{| a |}} ψ (\frac{t - b}{a}), a, b \in R, a \neq 0

(1)

where a is the scale factor, b is the translation factor, and $ψ_{a, b} (t)$ is derived by shifting and expanding successively the basic function $ψ (t)$ . If a and b are constantly changing, a functions family can be obtained, which is named as wavelet basis functions or wavelet bases. Given the square integrable signal $x (t)$ , that is, $x (t) \in L^{2} (R)$ , the wavelet transform of $x (t)$ is defined as follows

\begin{array}{l} W T_{x} (a, b) = \frac{1}{\sqrt{a}} \int x (t) ψ^{*} (\frac{t - b}{a}) d t \\ = \int x (t) ψ_{a, b}^{*} (t) d t = 〈 x (t), ψ_{a, b} (t) 〉 \end{array}

(2)

where a, b, and t are continuous variables, and then equation (2) is also called continuous wavelet transform. $ψ_{a, b}^{*} (t)$ is the complex conjugate function of $ψ_{a, b} (t)$ . The wavelet transform $W T_{x} (a, b)$ of $x (t)$ is a function of a and b; it can be interpreted as an inner product of $x (t)$ and a wavelet bases family. The role of b is used to determine the time position, that is, the time center. The role of scale factor a is used to expand basic wavelet $ψ (t)$ . When $a > 1$ , $ψ (t)$ becomes $ψ (t / a)$ . If a is larger, the time domain support of $ψ (t / a)$ is more than that of $ψ (t)$ . On the contrary, when $a < 1$ , a is smaller and the width of $ψ (t / a)$ is narrower. In this way, the central position of the $x (t)$ analysis and the time width are determined by a and b. This change justly adapts to meet the basic requirement of different resolutions in different frequency ranges for signal analysis, which is suitable for dealing with mutation signal and time-varying signal, frequency modulation, or amplitude modulation signal.

Taking $a = a_{0}^{m}, b = n b_{0} a_{0}^{m}, a_{0} > 1, and b_{0} \in R$ , the discrete wavelet transform of the signal $f (t)$ is as follows

W_{ψ} f (m, n) = a^{\frac{m}{2}} \int_{- \infty}^{\infty} f (t) \bar{ψ} (a_{0}^{- m} t - n b_{0}) d t

(3)

When $a_{0} = 2$ and $b_{0} = 1$ , equation (3) becomes a two progressive wavelet transform as follows

W_{ψ} f (m, n) = 2^{- \frac{m}{2}} \int_{- \infty}^{\infty} f (t) \bar{ψ} (2^{- m} t - n) d t

(4)

The reconstruction function of $f (t)$ is expressed as follows

f (t) = C_{ψ}^{- 1} \iint_{R^{2}} W_{ψ} f (a, b) ψ_{a, b} (t) \frac{d a d b}{a^{2}}

(5)

Currently, Haar wavelet, Mexican hat wavelet, Morlet wavelet, Daubechies wavelet, and spline wavelet are widely used. Based on practical experience, the Morlet wavelet is widely used in signal representation and classification, and image recognition feature extraction; the Haar wavelet is widely used in image processing and digital watermarking; the Mexican hat wavelet is used for system identification; the spline wavelet is used for material flaw detection; and the Shannon orthogonal basis is used for difference equation solution. Therefore, the Morlet wavelet is used here. Wavelet decomposition only implements decomposition for the low-frequency space step-by-step, and the high-frequency space W_j is no longer decomposed. Comparing $W_{1}$ and $W_{2}$ , $W_{1}$ has the best time domain resolution, but the worst frequency domain resolution. This does not satisfy obtaining good time resolution and good frequency resolution simultaneously, but wavelet packet decomposition can meet this need. Wavelet packet decomposition is carried out by the orthogonal image filter, and by setting the sampling signal as $y (t)$ , the following recursive equations are obtained

{\begin{matrix} y_{2 n} (t) = \sqrt{2} \sum_{k} h (k) y_{n} (2 t - k) \\ y_{2 n + 1} (t) = \sqrt{2} \sum_{k} g (k) y_{n} (2 t - k) \end{matrix}

(6)

The function ${y_{n} (t)}$ is called orthogonal wavelet packet; it is the decomposition result of the original signal at different scales for all bands. The time–frequency characteristics of signals can be well reflected by selecting the best base, which reflects self-adaption of wavelet packet decomposition.

In view of the strong vibration phenomenon of the rolling mill during the rolling of the thin strip, a comprehensive test was carried out on the spot, mainly collecting the vibration acceleration signals of the rollers in the transverse and vertical directions. The sampling frequency of all the measured signals is 512 Hz, which will not be commented on later in this article. Figure 1 is the details signal and power spectrum diagram of the five-layer db3 decomposition of an equipment measured vibration waveform. Figure 2 is the approximation signal. Figure 3 is the time domain waveform and its power spectrum after reconstruction. The amplitude or spectral values in Figure 3 are the regular normalization result. From the details signal in Figure 1 and approximation signal in Figure 2, every 450 sampling points (time interval is about 0.8 s) produces an impact with increasing amplitude. While the details signal d3 (64–32 Hz) showed the gradual process of the signal, the details signal d5 (16–8 Hz) displays the general features of the signal. From the reconstructed signal in Figure 3 and the measured signal, the two are almost identical. So, the wavelet decomposition method is credible.

Figure 1.

Five-layer db3 decomposition of the details signal and its frequency envelope spectrum of an equipment vibration signal.

Figure 2.

Five-layer db3 decomposition of the vibration signal.

Figure 3.

Time domain and power spectrum of the vibration signal after reconstruction.

But it is important to realize that due to the influence of background signal and noise signal, we cannot diagnose rubbing fault directly using wavelet transform processing. And because the rubbing fault signal and the random noise have the broadband frequency domain, we do not know the energy concentrated frequency range of rubbing avoidance signal. So by using the traditional filtering method, we cannot eliminate the noise signal and retain the fault information. But we can use the singular value decomposition method by attractor track matrix of time series reconstruction for eliminating the background signal and extracting the residual signal (noise and rubbing fault signal) and strengthening the fault signal and wavelet transforming the residual signal. The default threshold and given soft threshold denoising methods are used to filter the signal background noise, respectively, and it is found that the given soft threshold denoising method is better.

Nonlinear time series prediction

Prediction is an estimation of a thing or event behavior features or trend; the prediction of time series is defined by its future value and the number of properties according to the past and present value of the time series. Then we can carry out the analysis, prediction, and control of system behavior. Time series prediction of dynamic system can be divided into two stages: the first is to find the initial reference state data using the known data information, identification and optimization of various parameters in the prediction model, and the optimal approximation of the dynamic system evolution, and this process is called model learning stage. The second is to input the analyzed data into prediction model and obtain the prediction value of dynamic systems’ future evolution results, and this process is called data prediction stage. The existing time series prediction methods, according to the initial reference state data prediction model, can be divided into two categories: one is to find the initial reference state from the one-dimensional time series directly, and the other is to reconstruct the phase space of one-dimensional time series and select the initial reference state from phase space prediction model. At present, the main methods used to predict chaotic time series are the global domain method, local method, adding-weight one-rank local-region method, and maximum Lyapunov exponent prediction method.^10–14

Local prediction method based on phase space reconstruction

According to phase space attractors, the phase space reconstruction prediction method can be divided into global method and local method. The computation of the global prediction method is complicated, and the prediction accuracy will decline rapidly when choosing high embedding dimension, so the local method is used in this article. And the calculation principle is as follows. Suppose $x (t_{N})$ is a point in time series, and $x (t_{N} + 1)$ is the one-step prediction for one-dimensional time series. Suppose $Y (t)$ is the corresponding point $x (t_{N})$ of system phase space trajectory after the reconstruction of one-dimensional time series. Setting $Y (t)$ as the center point and some trace points $Y (t_{1}), Y (t_{2}), \dots, Y (t_{k})$ , which are the nearest center points, as the related points, further evolution points $Y (t_{1} + 1), Y (t_{2} + 1), \dots, Y (t_{k} + 1)$ of the related points are found. Then, according to the mapping relation between these related points and their evolution points, the mapping relationship $\hat{f}$ between $Y (t)$ and its further evolution point $Y (t + 1)$ is fit out, and then $Y (t + 1)$ is calculated. The final one-dimensional component $x (t + (m - 1) τ + 1)$ of $Y (t + 1)$ is the predicted value $x (t_{N} + 1)$ of the one-dimensional time series.

Compared with the global prediction method, the local method is applicable in most cases. The local prediction method has a very obvious physical interpretation, that is, the evolution law of the adjacent points is consistent in the predictable time step range, so we can describe the center point evolution according to the evolution of neighboring points. According to the different fitting methods of $\hat{f}$ , the local prediction method can be divided into the weighted zero-order local method, the weighted first-order local method, the local multiple linear regression method, and so on. In this article, the prediction results of the strong vibration signal and the slight vibration signal by the weighted zero-order local method are shown in Figures 4(c), (f) and 5(c), (f).

Figure 4.

Prediction of strong vibration signal: (a)–(c) horizontal vibration and (d)–(f) vertical vibration.

Figure 5.

Prediction of slight vibration signal: (a)–(c) horizontal vibration and (d)–(f) vertical vibration.

Prediction method based on maximum Lyapunov exponent

According to the measured vibration signal of rolling mill, it has been determined that it has nonstationary characteristics. In order to develop the rolling mill vibration prediction system in industrial production, it is necessary to study the advantages and disadvantages of various time series methods in rolling mill vibration prediction and make a comparative analysis.

The basic characteristic of chaotic motion is that the motion is very sensitive to initial condition, and the orbits produced by two very close initial values are separated by exponential form. Lyapunov index is a quantitative description of the characteristics of this phenomenon, and it is a very good prediction parameter for the nonlinear system.¹⁵

The basic idea of the maximum Lyapunov exponent prediction method is as follows: suppose $Y (t_{M})$ is the reference center points for predicting in phase space and $Y (t_{n b t})$ is the nearest neighboring point of $Y (t_{M})$ , after the one-step time evolution, the corresponding evolution points $Y (t_{M} + 1)$ and $Y (t_{n b t} + 1)$ are obtained, and they are as follows

e^{λ_{1}} = \frac{‖ Y (t_{M} + 1) - Y (t_{n b t} + 1) ‖}{‖ Y (t_{M}) - Y (t_{n b t}) ‖}

(7)

where $λ_{1}$ is the maximum Lyapunov index of nonlinear time series, ‖.‖ indicates the distance between two points in the phase space, and $Y (t_{M} + 1)$ is the forecast object in phase space. The coordinates of the last dimension $x (t_{M} + 1)$ are unknown; therefore, equation (7) can be solved. The solution $x (t_{M} + 1)$ is a one-dimensional time series prediction value. In this article, the prediction results of the strong vibration signal and the slight vibration signal by maximum Lyapunov exponent method are shown in Figures 4(a), (d) and 5(a), (d).

One-step prediction model of weighted first-order local method

For time series $x_{1}, x_{2}, \dots, x_{N}$ , N is the length of the sequence, and the correlation dimension d of time series is calculated according to G-P (Grassberger–Procaccia) algorithm. Then we select the embedding dimension $m = 2 d + 1$ by Takens’ theorem and acquire the averaged orbital period mt by analyzing the signal using the fast Fourier transform (FFT) method. The time delay $τ$ is obtained by $m t = (m - 1) τ$ , and then the reconstruction phase space is as follows: $Y_{i} = (x_{i}, x_{i + τ}, \dots, x_{i + (m - 1) τ})$ where i = 1, 2, . . ., M. M is the number of phase points in the reconstructed phase space, and $M = N - (m - 1) τ$ . Suppose the adjacent point of the center $Y_{M}$ is $Y_{M i}$ , i = 1, 2, . . ., q, and the distance from $Y_{M}$ is $d_{i}$ . Suppose $d_{\min}$ is the minimum value in $d_{i}$ , the weight definition of the point $Y_{k i}$ is as follows

P_{i} = \frac{\exp (- (d_{i} - d_{\min}))}{\sum_{i = 1}^{q} \exp (- (d_{i} - d_{\min})}

(8)

The first-order local linear fitting is as follows

Y_{M i + 1} = a e + b Y_{M i} ， i = 1, 2 ， \dots, q

(9)

where a and b are the real fitting coefficients; $e$ is a q-dimension vector, $e = {(1, 1, \dots, 1)}^{T}$ ; and $Y_{M i + 1}$ is a phase point after an evolutionary step of $Y_{M i}$ .

If the embedding dimension is $m = 1$ , equation (10) after using the weighted least square method is as follows

\sum_{i = 1}^{q} P_{i} {(x_{M i + 1} - a - b x_{M i})}^{2} = \min

(10)

The coefficients a and b can be obtained after solving the equations, and by substituting a and b into one-step prediction formula $Y_{M + 1} = a e + b Y_{M}$ , we can obtain the phase point prediction $Y_{M + 1}$ after one-step evolution

Y_{M + 1} = (x_{M + 1}, x_{M + 1 + τ}, \dots, x_{M + 1 + (m - 1) τ})

(11)

where the former (m–1) elements in $Y_{M + 1}$ are known values of the original sequence, and the mth element $x_{M + 1 + (m - 1) τ}$ is the one-step predictive value of the original sequence ${\hat{x}}_{N + 1}$ . This is the one-step prediction model of the weighted first-order local method.

For multi-step prediction, we can use the new information to add the original time sequence and repeat the steps. But the calculation is tremendous, and there is an error accumulation effect. Therefore, the weighted first-order local method of the multi-step prediction model can be used; the essence of this model is that the (m + 1) phase points are mostly similar to the reference point found in the reconstructed phase space. And we can carry out one-step prediction according to the evolution step of these (m + 1) phase points. For the chaotic time series, the evolution law over time of a pair of nearest neighbors in the phase space follows an exponential law as $e^{- λ t}$ . The parameter $λ$ is the largest Lyapunov exponent, which reflects the exponential divergence rate of the initial closed orbit. Therefore, when the embedding dimension $m > 1$ and if $k$ $(k > 1)$ step prediction is carried out, it is similar to the weighted one-order local law of the one-step prediction model. The k prediction can be implemented according to the evolution k step law of the (m + 1) phase points, and its specific derivation is as follows.

Suppose the reference vector set of the center point is ${Y_{M i}}$ , i = 1, 2, . . ., q, and the phase point set of evolution k step is ${Y_{M i + k}}$ , the one-order local linear fitting is as follows

Y_{M i + k} = a_{k} e + b_{k} Y_{M i}, i = 1, 2 ， \dots, q

(12)

According to the weighted least square method

\sum_{i = 1}^{q} P_{i} [\sum_{j = 1}^{m} {(x_{M i + k}^{j} - a_{k} - b_{k} x_{M i}^{j})}^{2}] = \min

The above formula is regarded as the two-variable function of the unknown $a_{k}$ and $b_{k}$ , and we can implement the partial derivative on their both sides

{\begin{cases} \sum_{i = 1}^{q} P_{i} \sum_{j = 1}^{m} (x_{M i + k}^{j} - a_{k} - b_{k} x_{M i}^{j}) = 0 \\ \sum_{i = 1}^{q} P_{i} \sum_{j = 1}^{m} (x_{M i + k}^{j} - a_{k} - b_{k} x_{M i}^{j}) x_{M i}^{j} = 0 \end{cases}

Simplify them as follows

{\begin{cases} a_{k} \sum_{i = 1}^{q} P_{i} \sum_{j = 1}^{m} x_{M i}^{j} + b_{k} \sum_{i = 1}^{q} P_{i} \sum_{j = 1}^{m} {(x_{M i}^{j})}^{2} \\ = \sum_{i = 1}^{q} P_{i} \sum_{j = 1}^{m} x_{M i + k}^{j} x_{M i}^{j} \\ a_{k} m + b_{k} \sum_{i = 1}^{q} P_{i} \sum_{j = 1}^{m} x_{M i}^{j} = \sum_{i = 1}^{q} P_{i} \sum_{j = 1}^{m} x_{M i + k}^{j} \end{cases}

Write it in matrix form as follows

(\begin{matrix} \begin{array}{l} c o e 1 \\ m \end{array} & \begin{array}{l} c o e 2 \\ c o e 1 \end{array} \end{matrix}) (\begin{array}{l} a_{k} \\ b_{k} \end{array}) = (\begin{array}{l} e_{k} \\ f_{k} \end{array})

where $c o e 1 = \sum_{i = 1}^{q} P_{i} \sum_{j = 1}^{m} x_{M i}^{j}$ $c o e 2 = \sum_{i = 1}^{q} P_{i} \sum_{j = 1}^{m} {(x_{M i}^{j})}^{2}$ , $e_{k} = \sum_{i = 1}^{q} P_{i} \sum_{j = 1}^{m} x_{M i + k}^{j} x_{M i}^{j}$ , and $f_{k} = \sum_{i = 1}^{q} P_{i} \sum_{j = 1}^{m} x_{M i + k}^{j}$ .

Then

(\begin{array}{l} a_{k} \\ b_{k} \end{array}) = {(\begin{matrix} \begin{array}{l} c o e 1 \\ m \end{array} & \begin{array}{l} c o e 2 \\ c o e 1 \end{array} \end{matrix})}^{- 1} (\begin{array}{l} e_{k} \\ f_{k} \end{array})

In programming calculation, the item $\sum_{j = 1}^{m} x_{M i}^{j}$ of the coefficient $c o e 1$ is element accumulation of the reference vector $Y_{M i}$ . The item $\sum_{j = 1}^{m} {(x_{M i}^{j})}^{2} = Y_{M i} Y_{M i}^{T}$ of coefficient $c o e 2$ , the item $\sum_{j = 1}^{m} x_{M i + k}^{j} x_{M i}^{j} = Y_{M i + k} Y_{M i}^{T}$ of $e_{k}$ , and the item $\sum_{j = 1}^{m} x_{M i + k}^{j}$ of $f_{k}$ are element accumulation of the reference vector $Y_{M i + k}$ . Then, the formula of $a_{k}$ and $b_{k}$ is obtained as follows

\begin{array}{l} (\begin{array}{l} a_{k} \\ b_{k} \end{array}) = {(\begin{matrix} \begin{array}{l} \sum_{i = 1}^{q} P_{i} \sum_{j = 1}^{m} x_{M i}^{j} \\ m \end{array} & \begin{array}{l} \sum_{i = 1}^{q} P_{i} Y_{M i} Y_{M i}^{T} \\ \sum_{i = 1}^{q} P_{i} \sum_{j = 1}^{m} x_{M i}^{j} \end{array} \end{matrix})}^{- 1} \\ (\begin{array}{l} \sum_{i = 1}^{q} P_{i} Y_{M i + k} Y_{M i}^{T} \\ \sum_{i = 1}^{q} P_{i} \sum_{j = 1}^{m} x_{M i + k}^{j} \end{array}) \end{array}

(13)

According to $a_{k}$ and $b_{k}$ , and substituting them into k step prediction formula $Y_{M + k} = a_{k} e + b_{k} Y_{M}$ , we can obtain the phase point prediction value after k step evolution $Y_{M + k}$

Y_{M + k} = (x_{M + k}, x_{M + k + τ}, \dots, x_{M + k + (m - 1) τ})

(14)

where the mth element $x_{M + k + (m - 1) τ}$ of $Y_{M + k}$ is the k step predicted value ${\hat{x}}_{N + k}$ of the original sequence.

Chaos prediction algorithm

When the equipment is running, for every time $τ$ , the state parameters of the mechanical system are measured, and the state time series are obtained as follows^16,17

x (t) = t = i τ, \dots, i = 0, 1, 2, \dots

In addition, simplify it as $x_{0}, x_{1}, \dots, x_{j}$ .

2. Build the phase space. For the attractor characterized by the time series ${x (t_{i})} (i = 1, 2, \dots, n)$ , suppose the embedding dimension is d and the delay time is τ, the phase space reconstruction of the time series is performed as follows

\begin{matrix} x (t_{n}) & x (t_{n - 1}) & \dots & x (t_{1} + (d - 1) τ) \\ x (t_{n} - τ & x (t_{n - 1} - τ) & \dots & x (t_{1} + (d - 2) τ) \\ ⋮ & ⋮ & ⋮ & ⋮ \\ x (t_{n} - (d - 1) τ) & x (t_{n - 1} - (d - 1) τ) & \dots & x (t_{1}) \end{matrix}

(15)

So, $Y_{n} = (x (t_{n}), x (t_{n} + τ), \dots, x (t_{n} + (m - 1) τ))$ .

3. Calculate the fractal dimension, Lyapunov exponent, and Kolmogorov entropy of the phase space and determine whether the system is chaotic. We can use chaotic prediction if it is chaotic. If not, we can use another method.

4. The maximum predictable time scale T is calculated.

5. In the phase space, the Euclidean distance between each point and the center point $Y_{n}$ is calculated, and the nearest $p$ points are found. In the time scale T, $Y_{n + 1} = \sum_{i = 1}^{N} Y_{i_{k}} e^{- l (d_{i} - d_{m})} / \sum_{i = 1}^{N} e^{- l (d_{i} - d_{m})}$ is obtained by the weighted local prediction method.

6. According to the phase space reconstruction theory, there are $Y_{n} = (x (t_{n}), x (t_{n} + τ), \dots, x (t_{n} + (m - 1) τ))$ and $n = N - (m - 1) τ$

And $x (t_{n + 1} + (m - 1) τ)$ is corresponding to $x (T_{N + 1})$ , so we take out the final one-dimensional sequence of $Y_{n + 1}$ , and then the predicted value of $x (T_{N + 1})$ can be obtained.

In this article, the prediction results of the strong vibration signal and the slight vibration signal are shown in Figures 4(b), (e) and 5(b), (e).

As seen from Figures 4 and 5, by comparing the deviation between the predicted value and the measured value, it can be seen that one-time multi-step prediction method is better than the other two methods.

Conclusion

In order to overcome the shortcomings of the Fourier transform (which cannot analyze the time–frequency localization and nonstationary signal), in this article, according to the actual nonstationary vibration signal of a large equipment located in a large Steel Corporation, the research group used wavelet packet decomposition and reconstruction, time–frequency analysis, and the application of wavelet decomposition in signal denoising. The effect of the high-frequency weight coefficient selection of each layer on signal denoising was analyzed. The chaotic time series nonlinear prediction was carried out using the global method, local method, weighted first-order local method, and maximum Lyapunov exponent prediction method, and it was found that the multi-step prediction method is better than other prediction methods.

Footnotes

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) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was supported by Key Scientific and Technological project of Henan Province (No. 172102210022, No.182102210086), Natural Science Foundation of Henan Province of China (No.182300410265), and National Undergraduate Training Program for Innovation and Entrepreneurship (No. 201810460069).

ORCID iD

Xiao-bin Fan

References

Bucher

Estimating the ratio between travelling and standing vibration waves under non-stationary conditions. J Sound Vibrat 2004; 270(1): 341–359.

Law

Kim

Liew

WYH

, et al. An approach based on wavelet packet decomposition and Hilbert–Huang transform (WPD–HHT) for spindle bearings condition monitoring. Mech Syst Signal Process 2012; 33(2): 197–211.

Petsounis

Fassois

SD.

Non-stationary functional series tarma vibration modelling and analysis in a planar manipulator. J Sound Vibrat 2000; 231(5): 1355–1376.

Heo

Kim

KJ.

Definitions of non-stationary vibration power for time–frequency analysis and computational algorithms based upon harmonic wavelet transform. J Sound Vibrat 2015; 336(336): 275–292.

Kougioumtzoglou

Spanos

PD.

Harmonic wavelets based response evolutionary power spectrum determination of linear and non-linear oscillators with fractional derivative elements. Int J Non-Linear Mech 2015; 80: 66–75.

Spanos

Kong

, et al. Harmonic wavelets based excitation-response relationships for linear systems: a critical perspective. Probabil Eng Mech 2015; 44: 163–173.

Singh

Kumar

Thrust bearing groove race defect measurement by wavelet decomposition of pre-processed vibration signal. Measurement 2013; 46(9): 3508–3515.

Lesecq

Gentil

Taleb

Fault detection based on wavelets transform. Application to a Roughing Mill. IFAC Proc 2007; 39(13): 1115–1120.

Risse

Combining wavelet decomposition with machine learning to forecast gold returns. Int J Forecast 2019; 35(2): 601–615.

10.

Inoussa

Peng

Nonlinear time series modeling and prediction using functional weights wavelet neural network-based state-dependent AR model. Neurocomputing 2012; 86(1): 59–74.

11.

Bertrand

A class of affine Wigner functions with extends covariance properties. J Math Phys 1992; 33: 2515–2527.

12.

Qian

Chen

Decomposition of the Wigner–Ville distribution and time-frequency distribution series. IEEE Trans Signal Process 1994; 42(10): 2836–2842.

13.

Cuomo

Lapenna

Macchiato

, et al. Linear and nonlinear dynamics in electrical precursory time series: implications for earthquake prediction. Tectonophysics 1998; 287(1–4): 279–298.

14.

Burdzik

Multidimensional identification of resonances analysis of strongly nonstationary signals, case study: diagnostic and condition monitoring of vehicle’s suspension system. Appl Acoust 2019; 144: 51–63.

15.

Zhang

Lam

Yan

, et al. Time series prediction using Lyapunov exponents in embedding phase space. Comp Electr Eng 2004; 30(1): 1–15.

16.

Song

Multi-step prediction of chaotic time-series with intermittent failures based on the generalized nonlinear filtering methods. Appl Mathemat Comput 2013; 219(16): 8584–8594.

17.

Kim

Kwon

Song

Response prediction of nonlinear hysteretic systems by deep neural networks. Neural Netw 2019; 111: 1–10.