Abstract: The purpose of this study is to propose and investigate a new approach for extracting spectral information of motion response of offshore structures. The approach is based on applying Time-varying Autoregressive (TVAR) model. This study is virtually unexplored in offshore engineering field. In the literatures, a number of works have shown that spectral content are extracted using Discrete Fourier Transform (DFT) for the frequency-domain analysis. Here, we outline a practical algorithm for TVAR model which uses Expectation-maximization (EM) algorithm based Kalman smoother. Short time Fourier transformation and Hilbert transformation are used as benchmark. The method is then applied to sampled discrete displacements of a fixed platform as a time series generated from field measurements. All the methods reveal that the spectrum characteristics of sampled platform displacement are time- varying frequency and time- varying gain distribution. The results indicate that TVAR model using KS with EM algorithm is superior to other methods in tackling frequency or amplitude modulation and systems that have low frequency dynamics. It is also found out that the mean frequency derived from the Hilbert transform is lower 8.2%, around 4.8% for short time Fourier transformation and 6.2% for TVAR model than the FFT spectrum.
INTRODUCTION
An important characteristic of offshore structures either fixed or compliant is their motion responses. The Discrete Fourier Transform (DFT) has been the most widely used technique to extract the spectral contents on those motion responses. One of the drawbacks of the DFT is that it does not provide any information about the time at which a frequency component occurs. Nonetheless, when the signals are nonstationary, then the DFT is not applicable anymore. Leakage is also big problem for the DFT. Prior researches proved that application of the DFT is not recommended to process ocean wave elevation because of non-stationary and nonlinearity of sea states (Huang et al., 1998; Schlurmann, 2002; Liu, 2000) as well as structural motion responses. The application of the DFT also may affect the frequency response of offshore structures (Hwang et al., 2003).
Representation of spectral information of measured motion responses in time-frequency plane undoubtedly provides additional insights of the analyzed system that may have never been revealed using the DFT. Some tools that can be used for time-frequency analysis such as Short Time Fourier Transform (STFT), Wavelet Transform (WT) together with their variants are popular and have been implemented successfully. However, those non parametric approaches suffer of resolution conflict in both frequency and time domain due to Heisenberg uncertainty principle. The best solution is to employ time-varying spectral analysis which is not affected by resolution conflict. Generally, parametric approach based spectral contents extraction is a promising technique to solve such a problem.
This study proposes the application of TVAR model as an alternative method in extracting frequency content of motion responses for offshore structures. That is because spectral contents extraction is the basic stage for Response Amplitude Operators (RAO) estimation, response Transfer Function (TF) and coherence analysis. It is also then can be used for modal property analysis of existing offshore structures, damage identification and other purposes.
APPROACH AND METHODS
TVAR model is an extended Autoregressive (AR) model, but its coefficients are time-variant. As a parametric approach, TVAR model estimates time-varying spectrum by modeling the signals as a time-series. This realization enables to produce the poles of the system through TVAR coefficients. TVAR model in discrete time index k is given by:
| (1) |
Notation p represents the order of TVAR model and y(k) is the signal. Term ai(k) is the TVAR coefficients and e(k) is the prediction error term which is Gaussian with zero mean and variance σ2e.
Estimation of ai(k) can be computed through adaptive method and basis function approach (Sodsri, 2003). Until now, criteria for selecting the proper basis function is not available yet and still open research (Sodsri, 2003; Zhang et al., 2010; Nguyen et al., 2009; Khan and Dutt, 2007), while adaptive method is very popular due to its simplicity and generality. Hence, adaptive method is addressed in this study. To accommodate the use of adaptive method in estimation of ai(k), Eq. 1 must be converted into a measurement equation in vector notation as:
| (2) |
Notation C(k) = [y(k-1), …, y(k-p)] is the vector of the past measurements; vector x(k) = [a1(k), …, ap(k)]T is the array of TVAR coefficients and v(k) is the measurement noise with covariance matrix R. By simplifying the TVAR coefficients evolve over a time linearly and first-order Gauss-Markov process and then x(k) can be expressed as state equation as follows:
| (3) |
The term A is the state transition matrix and w(k) is the state noise with covariance matrix Q. Equation 2 and 3 represent a state-space model which enables the application of Kalman smoother in order to estimate TVAR coefficients. Both equations contain model parameters which are assumed before the application of the adaptive method. These parameters are initial conditions x0~N(μ0, Σ0), A, Q, R and denoted by θ = {A, Q, R, μ0, Σ0}. If some simplifications are introduced in Eq. 3, then the equation calls for two remarks:
| • | If there is no state noise in the state equation and state transition matrix A is constrained to a scaled identity matrix, then Eq. 3 can be written in Eq. 4 where state variable depends on the choice of the forgetting factor λ, expressed in Eq. 4: |
| (4) |
Equation 4 is called Adaptive Autoregressive (AAR) model and the only tuning parameter is λ and can be estimated with Least Mean Square (LMS) or Recursive Least Square (RLS) algorithm.
| • | If state equation in Eq. 3 is modeled as a random-walk model, then it can be expressed: |
| (5) |
Noise covariance matrix is constrained to an identity matrix: Q = 1pxpσ2w where σ2w is a noise state variance. The unknown parameter in random-walk model is σ2w.
Estimation of Eq. 4 using LMS and RLS algorithm had been investigated by Sodsri (2003). He revealed that the adaptive method under this class is sensitive to the noise and fail to track the systems with fast or broad frequency. Estimation of Eq. 5 was successfully carried out by Nguyen et al. (2009) using amplitude demodulation-Kalman smoother (AD-KS). However, both covariance matrix (Q and R) are set up manually. Further, the drawback of both models had been investigated by Khan and Dutt (2007) and they found out that it might deteriorate the performance of time-varying spectrum estimation and proposed the use of Kalman smoother with EM algorithm because of its superiority. This study addresses the use of Kalman smoother with EM algorithm which allows the model parameters θ to be estimated and employed in spectral estimation. As benchmark, STFT and Hilbert Transform (HT) are used.
KALMAN SMOOTHER WITH EM ALGORITHM
Kalman smoother with EM algorithm (KS with EM) is smoothed Kalman filter which is optimized with EM algorithm. Basic theory covers Kalman filter, smoothing equations and expectation-maximization of log-likelihood function. Kalman filter calculates the states x(k|k) and states covariance matrix P(k|k) in Eq. 5 in two stages: time update equations (predictor) and measurement update equations (corrector). The time update equations project the states and state covariance matrix estimates forward from time index k-1 to k, written in Eq. 6 and 7:
| (6) |
| (7) |
The measurement update equations incorporate a new measurement into the a priori estimate to obtain an improved a posteriori estimate. The first step in measurement update equations is to compute Kalman gain:
| (8) |
The next step is to compute an a posteriori state estimate x(k|k) as a linear combination of an a priori state estimate x(k|k-1), Kalman gain and weight difference between an actual measurement y(k) and a measurement prediction (C(k)x(k|k-1) as shown in Eq. 9. The difference (y(k)-C(k)x(k|k-1) is called residual. The residual reflects the discrepancy between the predicted measurement (C(k)x(k|k-1) and the actual measurement y(k):
| (9) |
The final step is to obtain an a posteriori error covariance estimate via Eq. 10:
| (10) |
The use of procedures above will generate lagged response of x(k|k)estimate. Smoothing equations can be used by reducing delay and decreasing the variance of states estimate (Nguyen et al., 2009). Combination between Kalman filter and smoothing equations is called Kalman smoother. Because the signal is processed offline, then fixed-interval smoother is applied in this paper. This method has performance to improve the accuracy of the states estimate (Khan and Dutt, 2007) and derived in Eq. 11 and 13:
| (11) |
| (12) |
| (13) |
EM algorithm is utilized to tune the model parameters θ, based on maximum likelihood of y(1:k) in the presence of hidden variables x(k|K), k = 1, …, K. EM algorithm consists of two steps. First step is calculation of the expected complete log-likelihood as a function of θ. The expected complete log-likelihood is expressed as follows:
| (14) |
The expected likelihood depends on three quantities below:
| (15) |
| (16) |
| (17) |
One quantity in Eq. 17, P(k,k-1|K) must be calculated through Eq. 18 while all the quantities in Eq. 15 and 16 are calculated using the Kalman smoother equations:
| (18) |
Second step is maximization by direct differentiation of F with respect to the θ. These two steps are applied iteratively until convergence achieved. The estimates for the model parameters θ are given by as follows:
| (19) |
| (20) |
| (21) |
| (22) |
| (23) |
It is noted that differentiation of F with respect to the θ can be referred to Khan and Dutt (2007) because of the lengthy of the expressions. Finally, TVAR coefficients are obtained in term x(k|K) and time-varying spectrum can be estimated as:
| (24) |
Term aj(k|K) is the jth element of the TVAR model coefficients, σ2e is the prediction error variance and f is the observed frequency.
NUMERICAL EXAMPLES
Here, numerical examples of synthetic signals as a verification of the method in estimating time-varying spectrum is presented.
| Fig. 1(a-b): | Numerical examples of time-varying spectra of (a) Linear chirp signal and (b) Jump signal |
Two kinds of synthetic signals are linear chirp signal and jump signal. Both synthetic signals have amplitude and frequency modulation. The frequency of chirp signal changes linearly from 0.05 to 2 Hz as shown in the top panel of Fig. 1a, while there is a jump from 0.05 Hz into 2 Hz at the time instant 50 sec in the top panel of Fig. 1b. These examples are a kind of an attempt in capturing systems that have slow varying dynamic parameters. Simulations are carried out in 100 sec with a 200 Hz sampling frequency and no measurement noise is injected.
In Fig. 1, HT and STFT as non-parametric approach clearly show their drawbacks in tracking nonstationary signal for the underlying system. In linear chirp signal, frequency estimate from using HT has tail effect as shown in the bottom panel of Fig. 1a. The bottom panel of Fig. 1b displays that frequency estimate from using HT bounces roughly after the jump, while STFT behave similarly before the jump. Compared to the others, TVAR model using KS with EM algorithm has better frequency and temporal resolution in tracking the sudden change of the signals.
APPLICATION TO FIELD MEASUREMENT DATA
Here, the application of the proposed method to analyze the field measurement data is demonstrated.
| Fig. 2: | Field measurement data of a fixed platform displacement |
Displacement time history of an offshore structure was recorded from a fixed platform in Malaysian water. The time history is depicted in Fig. 2. The data was recorded continuously with sampling period 0.29 sec. It is noted that statistical properties of the platform displacement record are not discussed here.
| Fig. 3: | Platform displacement spectrum estimation by FFT and AR model |
| Fig. 4: | Time-varying spectrum by TVAR model using KS with EM algorithm |
Time-series of Fig. 2 reveal that the platform displacements record seems to be almost linear random time-series because its amplitudes appear to be almost symmetric with respect to its mean value. The record also shows the nonstationary behaviour visually. It implies that the platform displacements in the rough or extreme sea state contains more nonlinearity and nonstationarity.
Spectral contents are calculated by FFT and AR model as basic information. Solution of AR model is done using Yule-Walker equation and solved by Levinson-Durbin algorithm. It is found that AR model has model order of 77 to fully capture the dynamics of the platform displacements record as well as from those obtained by FFT.
| Fig. 5: | Time-varying spectrum by HT |
This optimum model order is determined by Akaike’s Information Criteria (AIC). Figure 3 depicts the spectra obtained from the two methods.
The similarities and differences of the both methods in calculating the spectral peak frequency can be observed clearly. It can be seen that the dominant peak frequency between the FFT and the AR model is almost similar which is close to 0.12 Hz. The spectral magnitude of the FFT spectrum is higher around 1.28 times than the AR model magnitude. However, the AR model produces sharper and smoother spectrum than FFT. The result shows that estimation of the spectral peak frequency using AR model enables the frequency components to be determined more exactly than FFT in situation where the signal is contaminated with noise. However, both methods only give averaged spectrum and cannot give information on time localization.
This is the reason why the time-varying spectrum is presented. Time-varying spectrum of the platform displacement record in time-frequency plane using the HT, STFT and TVAR model via KS with EM algorithm are displayed in Fig. 4-6, respectively. One similarity is observed clearly in those Fig. 4-6, that the measured platform displacement has nonstationary characteristic. Time-varying spectrum obtained from the three methods has the mean frequency as shown in Fig. 4-6, denoted by dashed line (- - -). Model order of 2 is adequate for the TVAR model.
As benchmark for the TVAR model, the HT and STFT are performed. It should be noted that in STFT, every window is overlapped by 50% and multiple window procedure is carried out under Kaiser window. Compared to the TVAR model spectrum, the STFT and HT spectrum produce a rather broad time-frequency band.
| Fig. 6: | Time-varying spectrum by STFT |
| Fig. 7: | Temporally averaged time-varying spectrum |
However, the time-varying frequency obtained from the HT method is more transient than the others. It might be the effect of frequency modulation as stated by states (Huang et al., 1998; Schlurmann, 2002).
The marginal spectrum averaged over the whole time sequence (mean spectrum) of all methods is depicted in Fig. 7. This figure underlines results depicted in Fig. 4-6. There is difference in shape and magnitude among the methods. The mean spectrum obtained from TVAR model is much sharper than the others, while the HT produces the flattest spectrum. Each method produces different magnitude, where the HT has the lowest and the STFT has the highest magnitude, followed by TVAR model at the middle, respectively.
| Fig. 8: | Time-varying spectrum by TVAR model in 3D distribution using KS with EM algorithm |
If the FFT spectrum is taken as reference, the mean frequency derived from the HT is lower 8.2%, around 4.8% for STFT and 6.2% for TVAR model. From Fig. 7 also can be observed that even STFT is Fourier-based method, its spectrum is slightly different with the FFT spectrum due to multiple windowing procedure. It is strongly believed that the more stationary or nonlinear the platform displacement is (i.e., in rough or extreme sea state), the more the difference between the non Fourier-based methods and Fourier-based methods will be. These finding results certainly have impacts on the ocean engineering designs such as RAOs calculation.
3D distribution of time-varying spectrum provides good description as shown in Fig. 8. As observed in that figure, non-stationary is evident. The time evolution of the spectral peak with different intensity appears within the considered frequency range. In 3D distributions, it can be seen that the displacement record characteristics are time-varying frequency and time-varying magnitude distributions. TVAR model using KS with EM algorithm can extract such informations accurately and produces high resolution.
The validation of TVAR model output is displayed in the top panel of Fig. 9. The true means the original displacement record and the estimate is the prediction result of TVAR model. The modeling error of TVAR models is relatively small as shown in the bottom panel of Fig. 9, meaning that TVAR models can fit the displacement record accurately.
| Fig. 9(a-b): | TVAR model (a) Output and (b) its modeling error |
Accurate fitting will produce the accurate poles of the system through TVAR coefficients. It implies to the accurate time-varying spectrum of the displacement record.
CONCLUSION
Field measurements based time-varying spectrum has been estimated in this study using TVAR model. Performance of the TVAR model using KS with EM algorithm compared to the Hilbert transformation and STFT has been carried out. The results show that TVAR model using KS with EM algorithm is superior to other methods. All TVAR model-based methods, together with STFT and HT reveal that sampled offshore structure responses (displacement) is time-varying frequency and time- varying magnitude distributions. TVAR model via KS with EM algorithm can estimate the time-varying spectrum with high resolution. This study recommends that TVAR model is prospective for spectral analysis of offshore structures responses, especially for offshore structures vibration monitoring.
ACKNOWLEDGMENTS
The authors are thankful to Universiti Teknologi PETRONAS for providing the research facilities.