Estimation of Evolutionary Spectra of Monitored Seismic Ground Motions by Transformation of Correlation Functions

Spatially varying seismic ground motions recorded by distributed structural health monitoring systems (SHMS) can be used to improve the performances of civil engineering structures, necessitating estimation of the evolutionary power spectral density as an indispensable procedure for utilizing records of SHMS. This paper proposes a method for the estimation of evolutionary power spectral density of a nonstationary process by transforming the correlation functions of its sample time histories. First, the background of the theory of evolutionary power spectral density is reviewed in detail. Relationship between the EPSD and the correlation function of a reference stationary process is then established. Formulas are derived for estimating this correlation function directly from nonstationary sample time histories so that the EPSD can be obtained.The implementation procedure of the proposed method is also detailed. Finally, a numerical example is presented, which validates the proposed method, demonstrates its application for SHMS, and displays its capabilities by comparison with the traditional method.


Introduction
With the emergence of large-scale infrastructures, such as long-span bridges, skyscrapers, and space structures, structural health monitoring system (SHMS), which aims at enhancing the safety and reliability and reducing the operation and maintenance costs, has been an active research regime [1][2][3][4][5].One necessary application of the SHMS is to record and analyze structural responses under extreme natural effects such as earthquake and wind storms.Estimating the statistical characteristics of these natural effects is therefore an indispensable procedure for exploiting SHMS.Because these natural effects, especially the seismic ground motion records, usually display nonstationary features, methods for accurately estimating some useful statistical characteristics of nonstationary seismic ground motions are of considerable importance for the application of seismic records given by SHMS equipped on long-span civil engineering structures.
For nonstationary stochastic processes such as seismic ground motions or downburst winds, their time-frequency properties or time-dependent spectra attracted lots of researches in the past fifty years.Some existing approaches for depicting these properties are listed as follows: the evolutionary power spectral density proposed by Priestley [6], which is interpreted as energy distribution over frequency at different time for those can be regarded as quasi-stationary process; the short-time Fourier transform (STFT), obtained by using a sliding window for filtering [7]; the Wigner-Ville distribution (WVD), which is defined as the Fourier transform of an instantaneous correlation function [8]; the wavelet transform (WT) scalogram, which is represented by a group of varied sized wavelets [9]; the Hilbert spectrum acquired by Hilbert Huang transform (HHT), which is represented by Hilbert transform of intrinsic mode function (IMF) components [10]; the Chirplet transform (CT), which can be treated as the extension and combination of STFT and WT [11].Among them, due to its physical comparability with the power spectral density of stationary processes and its mathematical completeness, the approach of evolutionary power spectral density (EPSD) is most widely used in earthquake engineering for characterizing seismic ground motions.

Shock and Vibration
Although the concept of EPSD has been developed for a relatively long time period, the methods for estimating EPSD from sample time histories of seismic ground motion still remain only partially resolved.A variety of methods have been proposed for estimation of EPSD.Priestley proposed an original EPSD estimating method through windowed filtering [6].Conte and Peng introduced multiple-windowed STFT method in the EPSD estimation [12].Spanos and Failla proposed an innovative approach to estimate the EPSD of univariate nonstationary stochastic process via wavelet transform (WT) [13], which was later extended to multivariate cases by Huang and Chen [14].Chen estimated the EPSD of downburst wind through time-varying autoregressive mode [15].Kayhan proposed a data-adaptive method for EPSD estimation [16,17].Nevertheless, due to the basic difficulty of estimation of EPSD posed by the limited amount of sample time histories and by the Heisenberg uncertainty principle, none of these methods could perfectly be applicable for most situations.In addition, these methods only take advantage of the time series directly but not of the transformation of correlation function, because of the blurred relationship from nonstationary time-varying correlation function to EPSD.
This paper aims to propose an alternative method for EPSD estimation of a set of multivariate sample time histories, which yields the EPSD by transforming their estimated correlation functions.First, the background information of the theory of EPSD and its traditional estimation method is briefly introduced.On top of that, the relationship between the EPSD and the correlation function of a reference stationary process is established.A method is then derived for estimating this correlation function directly from nonstationary sample time histories, on the basis of which the EPSD can be obtained by transformation.Smoothing procedure necessary for eliminating random fluctuations is also presented.Then the implementation procedure of the proposed method of estimation of EPSD is detailed.Finally, a numerical example is presented in order to demonstrate the validity and the application of the proposed method.

Background of Evolutionary Power Spectral Density (EPSD)
Nonstationary seismic ground motions recorded by a SHMS installed on long-span structures usually are modelled as a one-dimensional N-variates (1D-Nv) nonstationary stochastic vector process X () () consisting of N components  () 1 (),  () 2 (), . . .,  ()  (), where the superscript () denotes the sample of a nonstationary stochastic vector process.For seismic ground motions, the mean values of X () () are set equal to zero without loss of generality.
The component process of X () () has the following representation [6]: where   (, ) is the time-frequency modulating function and where both d ()  () and d ()  (, ) denote Gaussian stochastic processes with orthogonal increments, which satisfy where ( * ) denotes complex conjugation and    () is the stationary power spectral density of the stationary process to be modulated by (1); and therefore where  0  (, ) is the auto-EPSD of  ()  ().Accordingly, the evolutionary power spectral density (EPSD) matrix of X() is given by The elements of the evolutionary power spectral density matrix are defined as  0  (, ) =   (, )   (, )    () , ,  = 1, 2, . . . ,, (6) where   (, ) is the modulating function and    () is the cross stationary power spectral density to be modulated.
The time-varying correlation matrix of X() is given by The relationship between S 0 (, ) and R 0 (,  + ) is shown by the following transformation: The theory of evolutionary processes proposed by Priestley [6] can be applied only for semistationary processes; that is, the nonstationary process can be regarded as "stationary" in a certain period of time Γ.Generally, this time period Γ is defined as follows [6]: Apparently, the value of Γ cannot be determined, which is a priori.

Priestley's Method for Estimation of EPSD.
For estimating the EPSD from sample time histories, traditionally the method proposed by Priestley [6] is invoked, which is briefly introduced here.This approach is showed by the following transforms: where (, ) is the estimation of EPSD, () is a window function, and () is a weight function, which are usually selected as ) . (12)

EPSD Estimation via Estimated Correlation Function
To utilize the correlation functions to estimate EPSD, firstly, a distinct relationship between EPSD and the correlation function of a reference stationary process of the original nonstationary process is derived, by which the problem is transformed to estimate the correlation function of that specified reference stationary process.Consequently, a set of formulas is derived for a specially designed method to estimate that correlation function from sample time histories and thus enables the estimation of EPSD.The estimated EPSD from a limited number of samples may possess large random fluctuations, and a smoothing procedure is presented to address this.For facilitating the application of the proposed method, its implementation procedure is also presented.

Transforming Correlation Function to EPSD.
With reference to the spectral representation of an evolutionary process in (1), we can define a family of stochastic processes  ()  ()  by fixing the time instant  as follows: that is, for a specified time instant ,  ()  ()  is a stochastic process with respect to the pseudo time  but not the real time , and  ()   ()  can be regarded as stationary since d ()   (, ) does not vary with the pseudo time .Consequently, defining   ()  as the cross power density spectra between the stationary process  ()   ()  and the stationary process  ()   ()  , with   ()  defined as the correlation function between  ()  ()  and  ()  ()  , the relationship between   ()  and   ()  can be depicted as a Fourier transformation pair: In this context, when the correlation function is calculated by estimating over a period T, for example, by then   ()  can be calculated by (14).On the other hand, using (13)   ()  has the following representation: By comparing (7) with (17), it can derived that once the pseudo time  is replaced with real time , the power spectral density of  ()  ()  can be equal to the EPSD of  ()  (): Jointly invoking ( 14)- (18), by estimation of the correlation function   ()  of the reference stationary process  ()  ()  , the EPSD of  ()  () can be readily obtained by Fourier transformation.However, since only the sample time histories of  ()   () but not of  ()  ()  are available, formulas for estimating  ()   ()  from   ()  need to be developed.

Formulas for Estimating Correlation Function.
Recalling that the evolutionary processes are supposed to be semistationary over a time period of Γ such that  − Γ/2 ≤  +  ≤  + Γ/2, the following approximate relationship can be presumed for EPSD: Utilizing this relationship and taking into consideration the relationship between EPSD and time-varying correlation function of  ()  () (9), we can derive that Equation ( 20) indicates that, under the assumption of local stationarity, the correlation function of the reference process   ()  can be approximated by the nonstationary timevarying correlation functions  0  (,  + ) in the time period −Γ/2 ≤  ≤ Γ/2.Therefore, the correlation function   ()  can be calculated approximately through the following formula: ()  ()  ()  ( + ) d. (21) Consequently, Ŝ  (, ) can now be estimated by the Fourier transform of R  ()  : 3.3.Smoothing.For estimating EPSD of nonstationary processes from the ground motions recorded by SHMS, an inherent difficulty is that the amount sample of time histories is very limited: usually one set of sample time histories is available.This entails considerable part of the stochastic error in the estimated EPSD Ŝ  (, ), which indicates that Ŝ  (, ) is highly unreliable and thus not applicable in practice.In order to reduce the fluctuations of Ŝ  (, ), two weight functions (V) and () must be included to smooth Ŝ  (, ) in the time domain and in the frequency domain so that where   (, ) represents a smoothed result of Ŝ  (, ), which significantly decreases stochastic errors although it may increase the bias of estimated EPSD.In consideration of the constraints due to limited samples,   (, ) can be used as an applicable estimation of EPSD.3.4.Implementation Procedure.In practice, suppose  sets of sample time histories X ()   (Δ) of the seismic ground motion recorded by SHMS are presented.For using the proposed method to estimate the EPSD, the following three steps are included.
(1) Estimate the correlation function R  ()  of the reference stationary processes from the given sample time histories: where  = Γ/2Δ, Δ is the time step, and P is the total amount of samples.
( where  = Γ/2Δ and Δ is the frequency step. (3) Smoothing the unsmoothed estimation of evolutionary spectra Ŝ  (, ) in both the time domain and the frequency domain, we obtain the final result:

Numerical Example
In this part, a numerical example is presented to demonstrate the capabilities of the proposed approach in estimating the EPSD of seismic ground motions.Due to the lack of multipoint measurements, the SHMS-recorded ground motions are represented by an ensemble of multivariate nonstationary time histories generated conforming to EPSD of a well-known measured seismic ground motion time history and the coherence function given by statistical results of measurements.Then the method proposed in this paper was adopted to estimate the ESPD of sample time histories.The estimated results are compared with the target EPSD of the generated samples as well as the corresponding results yielded by traditional estimation method.

Sample Time Histories of Seismic Ground Motions.
In this example, the ground motions were modeled as multivariate nonstationary stochastic vector processes according to the EPSD of the seismic ground motion record of El-Centro.
Figure 1 shows the time history of El-Centro wave.The evolutionary spectrum of the El-Centro wave is denoted by  El (, ) as showed in Figure 2. The ground motions at three points are treated as trivariate nonstationary stochastic vector processes.The configuration of the three points is displayed in Figure 3.For simplicity, the effect of wave propagation is not considered in this example.The simulated ground motion at point  is denoted by  ()   (), and the cross-ESPD of X () () are given by  0  (, ) = √  (, )   (, )  () , ,  = 1, 2, 3, ( where   () is the coherence function between  ()  () and  1 (, ) = 1.0El (, ),  2 (, ) = 1.2El (, ), and  3 (, ) = 0.8 El (, ).The Harichandran-Vanmarcke model [18] is chosen to describe the coherence functions   (): where   is the distance between point  and point  and () is frequency-dependent correlation distance: The spectral representation method (SRM) is employed to generate the sample time histories of X () () [19].Figure 4 shows a set of time histories of simulated ground motions at the three points.

Estimation by the Proposed Method.
The proposed method is used to estimate the EPSD and the correlation functions of the generated sample time histories of seismic ground motions.Parameters used for the estimation are as follows: the time step Δ = 0.02 s, the length of time period  = 20 s, the frequency step Δ = 1 rad/s, and the upper cut-off frequency   = 126 rad/s.The two window functions () and (V) are chosen as the following forms: where Ω = 12 rad/s and   = 1.6 s.
Figure 5 shows the auto-EPSD estimated (a) from one sample time history and (b) from the ensemble of ten sample time histories for Point 1.In each subfigure, the estimated EPSD with and without smoothing are plotted against the smoothed and unsmoothed target EPSD (i.e., the EPSD of Point 1 used for generated sample time histories).In Figure 5(a), it can be observed that, for only one sample, the smoothed estimated EPSD   (, ) is close to the smoothed target EPSD and is also close to the target EPSD  0  (, ), which indicates that the proposed method can estimate the EPSD from just one sample time history with satisfying accuracy.If the number of the samples can be increased, for example, a set of 10 samples are available, comparison between Figures 5(b) and 5(a) shows that even the unsmoothed EPSD Ŝ  (, ) can also be a satisfactory estimation of EPSD.That is, increasing the number of samples would increase the accuracy of the estimation of EPSD.Corresponding to Figures 5, 6  compare the mean of the EPSD estimated by the proposed method and the traditional method and the standard deviations of EPSD estimated by the two methods, respectively; each figure contains three subfigures associated with different time instants.From Figure 7 it can be observed that the biases between results of these two methods and the target EPSD are comparable, indicating that the proposed method can lead to accuracy similar to the traditional method.Meanwhile, Figure 8 displays the standard deviation of the EPSD estimated by the proposed method approach which is much lower than those of the traditional Priestley method, which demonstrates that results of the proposed method have much less fluctuations so that it is more reliable than those of the traditional method.

Conclusion
In this paper, a new approach has been proposed to estimate the EPSD of multivariate nonstationary seismic ground motions recorded by SHMS.The method is based on transformation of correlation functions of the reference stationary process, which can be estimated from sample time histories directly by using formulas derived.The smoothing procedure has also been included in the procedure method.The whole proposed method can be readily applied by using the implementation steps detailed.The validity and the accuracy of the proposed method have been demonstrated by a numerical example and the comparison between results yielded by this method and by the traditional Priestley method.This method could approximately estimate the EPSD from only one sample time history, while increasing the number of samples could efficiently increase the accuracy of the estimation dramatically.The biases of the EPSD estimated by these two methods are comparable, while the stochastic fluctuations of the proposed method are much lower than those of Priestley's method.It is worthy to point out that the proposed method is of suitable accuracy and more reliable with less random fluctuation than the traditional method and thus has potential application for exploiting data recorded by SHMS.

2 )Figure 4 :
Figure 4: A set of sample time histories of seismic ground motions.
(a) and 6(b) show the estimated cross EPSD between Point 1 and Point 2 from one and ten samples respectively, and observations similar to those from

Figure 7 :
Figure 7: Comparison between the bias errors of the proposed method  12 (, ) and Priestley's method  12 (, ): (a) at time instant t = 2 s; (b) at time instant t = 4 s; (c) at time instant t = 12 s.

Figure 8 :
Figure 8: Comparison between the stochastic errors of the proposed method  12 (, ) and Priestley's method  12 (, ): (a) at time instant t = 2 s; (b) at time instant t = 4 s; (c) at time instant t = 12 s.