Abstract
In this paper nonstationary vibrations are studied by means singular spectrum analysis (SSA) – a modelfree method of time series analysis and forecasting. SSA allows decomposing the nonstationary time series into trend, periodic components and noise and forecasting subsequent behavior of system. The method can be successfully used for processing the signals from the vibrating constructional elements and machine parts. This paper shows application of this method for random and nonlinear vibrations study on the examples of construction elements vibration under seismic action.
1. Introduction
Singular spectrum analysis (SSA) is new powerful method of time series analysis and forecasting, which was independently developed in Russia (called “Caterpillar” SSA) and also in UK and USA (under the name of SSA) [14]. The “Caterpillar” – SSA is a modelfree technique of time series analysis; it combines advantages of other methods, such as Fourier and regression analyses, with simplicity of visual control aids.
The “Caterpillar” – SSA approach is used to study the time series in many various spheres wherever the trend or periodic behaviour may present: in meteorology, climatology, hydrology, geophysics, medicine, biologyy, sociology. A lot of problems can be solved by means of “Caterpillar” – SSA technique: finding trends of different resolution, smoothing, extraction of the periodic components, changepoint detection, simultaneous extraction of cycles with small and large period; extraction of periodicities with varying amplitudes. Method is extensively developing: the modification of SSA for time series with missing data and the variants for analysis of multidimensional time series (MSSA) are elaborated.
Last decade the methods for nonstationary and nonlinear time series based on Empirical Mode Decomposition (EMD) are developed. EMD method used to decompose series as sums of zeromean components. The technique consists basically of computing an upper and lower spline envelope of the timeseries and extracting the average of this envelope up to when the residuals are random noise or stopping criterion is met. SSA method is more mathematically grounded, while EMD is more empirical. The aim of SSA is to make a decomposition of the original series into the sum of independent and interpretable components: trend – sum of slowly varying components, oscillatory components – sum of elementary harmonic components and a structureless noise – the random component formed by a set of random factors.
The “Caterpillar” – SSA method is useful tool for short and long, stationary and nonstationary, almost deterministic and noisy time series analyzing. The method can be successfully used for processing the signals from the vibrating constructional elements and machine parts. However, in mechanics application this method does not become widely used, but it may be successfully used for periodic motion detection, amplitude and frequency definition, as well as chaotic motion identification. The purpose of this work is to show the principle of the method, its ease of use and accessibility for mechanical applications. For this purpose the examples of vibration of construction element is considered: nonstationary vibration of pier abutment under seismic action.
2. Mathematical background
The main steps of “Caterpillar” – SSA algorithm described here, follow the methodology given in work [3]. The basic version of method consists in transformation of onedimensional series into multidimensional by oneparametric translation procedure, research of got multidimensional trajectory by means of principal components analysis (Singular Value Decomposition) and reconstruction of the series in accordance with the chosen principal components. The result of this method application is expansion of the time series into sum of simple components: slow trends, periodic or oscillating components and components of noises. The obtained decomposition may serve as basis of correct forecasting. Both of two complementary stages of SSA technique (decomposition and reconstruction) include two separate steps, therefore the analysis is performed by four steps.
At the first step (called the embedding step) a onedimensional series $F=({f}_{0},...,{f}_{N1})$, ${f}_{n}\in \mathbb{R}$, with length $N\ge 3$ is transferred to the $L$ – dimensional vector:
where $K=NL+1$ and $1<L<N\text{,}$$L\in \mathbb{N}$. $L$ is called by window length, it is very important parameter of the embedding step; it should be big enough but not greater than a half of series length. This delay procedure gives the first name to the whole technique.
Vectors ${\mathbf{X}}_{i}$ form columns of the trajectory matrix $\mathbf{X}=\left[{X}_{1}:\dots :{X}_{K}\right]$, or:
Trajectory matrix $\mathbf{X}$ is a Hankel matrix, which mean that all the elements along the diagonal $i+j=const$ are equal: ${x}_{ij}={f}_{i+j2}$.
The second step is the singular value decomposition (SVD) of the trajectory matrix into a sum of rankone biorthogonal elementary matrices: $\mathbf{X}={\mathbf{X}}_{1}+\dots +{\mathbf{X}}_{d}.$
Elementary matrix ${\mathbf{X}}_{i}$ is determined by the equality:
where $\sqrt{{\lambda}_{i}}$ ($i$th singular value) is the square root of the $i$th eigenvalue of the matrix $\mathbf{S}$, ${\mathbf{U}}_{1}\text{,...,}\text{}{\mathbf{U}}_{L}$ – orthonormal set of eigenvectors of matrix $\mathbf{S}$, $d=\mathrm{m}\mathrm{a}\mathrm{x}\{i:\lambda i>0\}$, ${\mathbf{V}}_{i}={\mathbf{X}}^{T}{\mathbf{U}}_{i}/\sqrt{{\lambda}_{i}}$. ${\mathbf{U}}_{i}$ and ${\mathbf{V}}_{i}$ are called as left and right singular vectors of the trajectory matrix. It is assumed that eigenvalues ${\lambda}_{i}$ are arranged in decreasing order of their magnitude$\left({\lambda}_{1}\ge \dots \ge {\lambda}_{L}\ge 0\right)\text{.}$ The collection $\left(\sqrt{{\lambda}_{i},}{\mathbf{U}}_{i},{\mathbf{V}}_{i}^{T}\right)$ is called the $i$th eigentriple of matrix $\mathbf{X}$.
The third step (grouping step) corresponds to splitting of the elementary matrices ${\mathbf{X}}_{i}$ into several groups and sunning the matrices within each group. Grouping procedure splits the set of indexes $\{1,\dots ,d\}$ into the disjoint subsets ${I}_{1}\text{,}\text{}{I}_{2}\text{,...,}\text{}{I}_{m}\text{.}$ Let ${I}_{k}=\{{i}_{1},...,{i}_{p}\}$ then ${\mathbf{X}}_{{I}_{k}}={\mathbf{X}}_{{i}_{1}}+...+{\mathbf{X}}_{{i}_{p}}$ under condition $I={\bigcup}_{k=1}^{m}{I}_{k}$. For given group ${I}_{k}$ the contribution of the component ${{\mathbf{X}}_{I}}_{k}$ into the expansion is measured by the share of the corresponding eigenvalues.
The last step (diagonal averaging) transfers each resultant matrix into a Hankel matrix, from which then the time series will be reconstructed (the additive components of the initial series $F$). Let ${y}_{ij}$ stands for an element of $L\times K$ matrix $\mathbf{Y}\text{,}$ where $1\le i\le L\text{,}$$1\le j\le K\text{.}$ Let ${L}^{*}=\mathrm{m}\mathrm{i}\mathrm{n}\left(L,K\right)\text{,}$${K}^{*}=\mathrm{m}\mathrm{a}\mathrm{x}\left(L,K\right)\text{,}$${y}_{ij}^{\mathrm{*}}={y}_{ij}$ if $L<K$ and ${y}_{ij}^{\mathrm{*}}={y}_{ji}$ otherwise. Diagonal averaging transforms matrix $\mathbf{Y}$ into series ${g}_{0}\text{,...,}\text{}{g}_{N1}$ ($N=L+K1$) by formulas:
${g}_{k}=\frac{1}{Nk}{\sum}_{m={K}^{\mathrm{*}}+2}^{N{K}^{\mathrm{*}}+1}{y}_{m,km+2}^{\mathrm{*}},{K}^{*}\le k<N.$
Elements of Hankel matrix is obtained by averaging of ${y}_{ij}$ over all $i\text{,}$$j$ such that $i+j=k+2$. In this way we obtain a decomposition of the initial series into several additive components. The result is the expansion of matrix ${\mathbf{X}}_{{I}_{k}}$ into simple components ${F}^{\left(k\right)}=\left({f}_{0}^{\left(k\right)},\dots ,{f}_{N1}^{\left(k\right)}\right)$ and initial series $F=\left({f}_{0},\dots ,{f}_{N1}\right)$ is decomposed into the sum of $m$ series:
Forecasting time series, unlike the analysis, is only possible involving the mathematical model of the series. The model must either originate from the data itself, or be tested on the available data. When forecasting with SSA method, we consider a set of time series described by linear recurrent formula (LRF). The class of time series governed by linear recurrent formulae is rather wide; it includes harmonics, polynomial and exponential time series. The time series governed by LRF generates recurrent extension as each its term is a linear combination of some of the previous ones. SSA forecasting method can be applied to the time series that approximately satisfy linear recurrent formulas.
For SSA procedure the correct choosing of window length $L$ and way of grouping of of the eigentriples (or grouping of elementary matrices, since each matrix component of the SVD is completely determined by the corresponding eigentriple) are very important. It is preferable to take the window length proportional to supposed period in order to achieve the better separability of the periodic components. As for the way of grouping, it is useful to mention that under the proper choice of window length $L$ singular vectors ‘repeat’ the behavior of corresponding time series components. In particular, trend of the series corresponds to slowly varying singular vectors; harmonic component produces a pair of left (and right) harmonic singular vectors with the same frequency, etc.
The main principles for identification of eigentriples corresponding to the time series component of interest are as follows: to extract the trend or to smooth the time series, slowly changing eigenfunctions are chosen; to extract oscillations, oscillating eigenfunctions are chosen; to extract harmonic components with period greater than 2 a pair of eigentriples with the same period should be found and chosen; to approximate the time series the leading eigentriples are chosen.
3. Numerical examples of SSA method application
In recent years corresponding computer software are elaborated for time series analysis: “Caterpillar” – SSA and Rssa package in St. Petersburg University [1, 3], SSAMTM Toolkit for Spectral Analysis, kSpectra Toolkit [5]. MATLAB program is also widely used. In this work MathCAD program are applied as easytouse environment for engineers and scientists research. MathCAD has a large number of builtin functions, including matrix, and programming tools. The analysis of horizontal and vertical displacement of the top of pier after seismic action is performed by means of “Caterpillar” – SSA method using MathCAD. Recording of the amplitudes (mm) was made with regular interval 0.04 s, hereinafter this data are called as time series.
3.1. Example 1
The analysis of horizontal vibration of the pier top is presented below, time series length is equal 248 ($N=$ 248). Initial time series and its trend (reconstructed with the first 4 eigentriples) are presented in Fig. 1, window length is set $L=$ 124. Plot in Fig. 2 shows the plot logarithms of eigenvalues in decreasing order; logarithm was introduced because of very large figures.
It is clear that the first four eigenvalues have the largest value and are able to explain the main part of variance. To determine the size of informative basis the next features of eigenvalue is used: the larger share of eqenvalue in the sum of elements of the eigenvalues vector, the more background information includes a projection on the corresponding eigenvector. Because a drop in values is observed around 14th eigenvalue, it may be supposed that there noise floor starts.
Plot in Fig. 2 may be used for periodical components finding. Since in ideal situation two main components with the same period (“sin” and “cos”), relevant to one eigenvalue, correspond to one sinusoidal component, in real situation these components correspond to close eigenvalues. Thus the “tread” may be observed in the plot: 56, 78, 910, 1112, i.e. four evident pairs, with almost equal leading singular values.
In Fig. 3 onedimensional plots of the first nine eigenvectors are given as the example. Analysis of eigenvectors shows that the first four eigenvectors are a trend components (slowly changing component), the 56 and 78 may be periodical components.
Fig. 1Initial series () and its trend (), reconstructed with first four eigentriples
Fig. 2Singular value spectrum: plot of logarithms of the first 60 from 124 eigenvalues
Fig. 3Onedimensional plots of the first nine eigenvectors
a)
b)
c)
d)
e)
f)
g)
h)
i)
Fig. 4 shows the examples of twodimensional projection of the eigenvectors; vectors number is indicated in superscript of function and argument. Analysis of plots in Fig. 4 confirms the conclusion that the 5th and 6th eigenvectors (and also 7th and 8th) refer to the periodical components.
In Fig. 5(a) the periodical components consisting of the 5th and 6th eigentriples is presented, in Fig. 5(b) – sum of periodical components from 5th to 12th eigentriples.
By the same way the rest eigenvectors may be investigated.
Fig. 4Twodimensional plots of eigenvectors
Fig. 5Plots of the initial series () and periodical components (): a) reconstructed with 5th and 6th eigentriples, b) reconstructed with eigentriples from 5th to 12th
a)
b)
3.2. Example 2
The analysis of vertical vibration of the pier top is presented in Figs. 611; time series length $N=$242, window length is set $L=$ 120. Initial time series is presented in Fig. 6, plot in Fig. 7 shows the logarithms of eigenvalues in decreasing order; after 14th component the drop in value occurs, it could be considered as start of the noise floor. Six evident pair with almost equal leading singular value correspond to six (almost) harmonic components of series: eigentriples pair 12, 34, 56, 78, 910 and 1314 are related to the harmonic components. Reconstructed series with first 14 eigentriples is given in Fig. 8, reconstructed series of noise (last 106 eigentriples) – in Fig. 9. The principal components of first six eigentriples are shown in Fig. 10, twodimensional plots of eigenvectors is given in Fig. 11.
This scatterplots of paired principal components tend to a center, it means that vibrations are damped. For nondamped vibration twodimensional plot appears as ring or circle. From this plots it is seen that eigenvectors 1112 are not periodical components.
Fig. 6Initial series of amplitudes
Fig. 7Singular value spectrum: plot of logarithms of the first 40 from 120 eigenvalues
Fig. 8Reconstructed series with first 14 eigentriples (periodical components)
Fig. 9Noise reconstructed series (15120 eigentriples)
Fig. 10Onedimensional plots of the first six eigenvectors
a)
b)
c)
d)
e)
f)
Fig. 11Twodimensional plots of eigenvectors
a)
b)
c)
d)
e)
f)
g)
h)
4. Conclusions
Examples of application of the SSA technique for analysis of onedimensional time series is presented in this work: task of analysis of nonlinear vibration of the elements of constructions under seismic action. Simple solution shows the convenience of this method, its usability to nonstationary timeseries investigation.
References

Danilov D., Zhigljavsky A. Principal Components of Time Series: “Caterpillar” Method. St. Petersburg University, St. Petersburg, Russia, http://www.gistatgroup.com/gus/, 1997.

Ghil M., et al. Advanced spectral methods for climatic time series. Reviews of Geophysics, Vol. 40, 2002, p. 11–141.

Golyandina N., Zhigljavsky A. Singular Spectrum Analysis for Time Series. Springer Briefs in Statistics, Springer, 2013.

Jolliffe I. T. Principal Component Analysis. SpringerVerlag, 1986.

SSA –MTM Toolkit for Spectral Analysis, http://www.atmos.ucla.edu/tcd/ssa/.