A Grid-Based Method to Represent the Covariance Structure for Earthquake Ground Motion

Spatial variation of earthquake ground motion is an important phenomenon that cannot be ignored in the design and safety of strategic structures. However, almost all the procedures for the evaluation of variation assumed that the random field is homogeneous in space. It is obvious that reality does not fully conform to the assumption. How to investigate the inhomogeneous feature of groundmotion in space is a challenge for researcher. A body-fitted grid-coordinates-basedmethod is proposed to estimate and describe the local spatial variations for the earthquake ground motion; it need not to make the assumption that the random field of earthquake is homogeneous in space. An analysis of spatial variability of seismic motion in smart-1 array monitored in Lotung, Taiwan demonstrates this methodology.


Introduction
Spatial variation of earthquake ground motion is an important phenomenon that cannot be ignored in the design and safety of strategic structures 1-5 .The seismic motion, which could be considered as the result of the complex wave propagation through a heterogeneous soil, is affected by many factors such as propagation path, source mechanism, and the amplification effect of the subsurface layers Thus, the record from the seismograph array is not the simple duplication of the traveling wave with the shift of time, but the difference will be amplified when the distance between stations becomes larger.How to make the description and estimation of correlation structure of spatial ground motion in efficient and robust way is a fundamental prerequisite for problems of simulation of earthquake ground motion and design of structures with spatially extended foundations and lifeline systems.
The development of procedure for evaluating the spatial instationarity has lagged well behind the development of procedure for evaluating the spatial stationary or isotropy/ homogeneous in space .The classical way to describe the spatial variation of earthquake ground motion problem is performed by using the correlation or coherence function 6-9 , however, we have the following.
1 It is assumed that the random field is homogeneous in space; it is obvious that reality does not fully conform to that assumption e.g., one bridge on rock and the other on alluvium .
2 The amplitude variability of seismic ground motion ought to be incorporated into the description of spatial variation of seismic motion.However, the correlation or coherence function only can reflect the phase variation 9 .
In the analysis of most spatial-temporal processes studies, there are few reasons to expect spatial covariance structures to be stationary over the spatial scales of interest.
Earthquake is a spatial-temporal field; the focus of the researchers are usually on the temporal nonstationary 10, 11 .The lessons learned during the devastating earthquakes of the last decade tell us the following: just as the character of temporal nonstationary has big influence on seismic response of the structure, the influence of spatial nonstationary on some spatially extended structure also cannot be ignored.In the design of structures, the reasonable input should take the feature of spatial nonstationary of earthquake into account.It is well known that nearly every significant discovery of feature of earthquake will bring huge breakthrough of seismic design of structure 12 ; the research of spatial nonstationary is of great significance.
In this study, a new method is proposed to describe the spatial variation of earthquake ground motion.This methodology was motivated by problems of how to investigate and visualize the variations of earthquake ground motion in space if the field is inhomogeneous.
The fruit of our research can be applied at least in two aspects: 1 to help design "optimal" networks of sampling stations for observing these distribution; 2 to investigate the feature of earthquake and as a fundamental prerequisite for the problems of producing properly correlated motions in inhomogeneous space to help the design of the long-span structure.
This paper explains our method, with descriptions of the application of the three principal tools: dynamic time warping distance, multidimensional scaling, and thin-plate spline interpolation, and we also present an application of the method to characterize the local spatial correlation of the ground motions which have been recorded and obtained from closely spaced arrays SMART-1 arrays Event-40 .

Outline of the Method
Our approach requires two tools: the Multidimensional Scaling MDS and the technique of spatial interpolation.The former generates a low-dimensional space two or three representation visualizing proximities of the sampling stations; the latter provides smooth mappings of the geographic representation of the sampling stations into their MDS representation.

The Dynamic Time Warping Distance
For the accelerogram record a ij in n station sample point , each has N independent realizations.The accelerogram record can be given in a n × N data matrix, The coherency/correlation function or coherency coefficient is often used as similarity measure of earthquake ground motion 13-16 ; the Pearson coherency function is defined as: where but some researchers 15, 16 have found that the simple and crude definitions of coherence imply high degree of variability of strong motion even for short separation.It may be from the limitations of coherency function or correlation coefficient itself.The correlation coefficient required that two series have the same length, the values of two series have pointto-point correspondence, and the weight of each pair of difference is equal.Due to such a correspondence, it may not be suitable to be applied to the similarity measurement of complex series with shift and stretching of amplitude.The dynamic time warping distance DTW proposed by Berndt and Clioffrd 17 is designated to depict the greatest similarity between series by calculating the minimum distance between them, which is defined as follows.
Let A a 1 , a 2 , . . ., a n and B b 1 , b 2 , . . ., b m be two series with the length of n and m, respectively, and an n × m matrix M can be defined to represent the point-to-point correspondence relationship between A and B, where the element M ij indicates the distance d a i , b j between a i and b j .Then the point-to-point alignment and matching relationship between X and Y can be represented by a time warping path: W w 1 , w 2 , . . ., w k , max m, n ≤ K < m n−1, where the element w k i, j indicates the alignment and matching relationship between a i and b j .If a path is the lowest cost path between two series, the corresponding dynamic time warping distance is required to meet DTW A, B min Mathematical Problems in Engineering d k d a i , b j is the Minkowski distance of the k element of the path W and is denoted as follows: when p 2 the distance between two series is called Euclidean distance.
Then the formal definition of dynamic time warping distance between two series is described as where indicates empty series, 2 : − indicates a subarray whose elements include the second element to the final element in an one-dimensional array, and d a i , b j indicates the distance between points a i and b j .The DTW distance of two time series can be calculated by the dynamic programming method based on accumulated distance matrix 17 , whose algorithm mainly is to construct an accumulated distance matrix: Any element r i, j in the accumulated matrix indicates the dynamic time warping distance between series A 1:i and B 1:j .Series with high similar complexity can be effectively identified because the best alignment and matching relationship between two series is defined by the dynamic time distance.
The different between DTW a and correlation coefficient can be observed in Figure 1.
To demonstrate the different performance of the DTW and the correlation coefficient r, an example is given here.In this example we calculate the DTW and the correlation coefficient r see Figure 2 and the result can be seen in Table 1 : Conditions A-B show that DTW can measure the variation of amplitude, However the correlation cannot.Conditions A-C show that correlation requires that the signal be synchronization or point-to-point correspondence, and the DTW has no such requirement.Conditions A-C show that correlation cannot calculate two signals with different length; however, the DTW can.Table 1: DTW and correlation coefficient under four condition.

Optimal Scaling
Multidimensional Scaling MDS is applied as a statistical technique to visualize dissimilarity data in this section.Let Δ δ ij and D d ij be two N × N matrices indexed by i and j, where the proximity or data value connecting object i with object j is represented by Δ, d ij refers to the Euclidean space composition map between objects i and j.The main objective of MDS is to represent these dissimilarities as distances between points in a low-dimensional space or called composition such that the distances d ij correspond as closely as possible to the dissimilarities δ ij 18 .

Mathematical Problems in Engineering
Classical Metric Multidimensional Scaling is a basic form of MDS.Classical MDS is called metric methods because the relationship between δ ij and d ij depend on the numerical or metric properties of the dissimilarities.It only works under the assumption that the geometrical model fits the data perfectly 18, 19 .However, it is often not possible to construct an explicit functional form D such that the mapped dissimilarities D of an empirical data set match sufficiently well metric distances.Therefore, we would like to assign numerical values to the optimal approximations of the transformed proximities to the distances in the geometrical representation.These numerical values are usually called disparities, and they are denoted by d and d ij ≈ d ij f δ ij ; it should be restricted by the monotonic constraint: The coordinates in the distance function in composition map and the function f which allows transforming the proximities into distances are estimated by minimizing the following badness of fit function usually called stress or S-function in the context of MDS or called optimal scaling.The stress function is given by

2.8
Optimal scaling aims to find a transformation of the data that fits as well as possible the distances in the MDS solution and find a matrix D by some iterative algorithms to make the S minimum.x 1i , x 2i , . . ., x pi , where p is the dimensional of the final composition map p 1, 2, . . ., n, etc. .The method is summarized as follows.
1 Select the initial matrix x 0 1i , x 0 2i , . 4 For the arbitrary i, j, k, l if where n, θ are respectively the number of the node the and iterative step.
6 Use the coordinate of node in step 5 to calculate the Euclidean d ij .
7 Calculate the S according to formula 2.8 .Finally, the result x 1i , x 2i , . . ., x pi i i 1, 2 is the solution of MDS.If we chose p 2, the coordinate can be expressed as x 1i , x 2i .We refer to the plane of geographical y 1i , y 2i of the sampling stations as G plane, and we refer to the plane of composition coordinates x 1i , x 2i as C plane see Figure 3 .

The Body-Fitted Grid-Based Depiction of Mapping
Classical research of investigating the variation of earthquake in space is based on the assumption that the ground motion is homogeneous in space or the function to describe the variation of earthquake is the functions of the separation distance between stations, but independent of absolute location .This simplifying assumptions may not always capture reality.How to capture the inhomogeneous feature of earthquake in space is a challenge but a foundation problem, to solve this problem; this paper introduces a body-fitted coordinatesbased method to represent the covariance structure of seismograph array.This method need not to set a reference point and can be taken as an alternative way to reflect the variation of seismic ground motion in inhomogeneous field.
Just as discussed above, the coordinate x 1i , x 2i in composition map C plane can reflect the difference or correlation of seismic ground motion; here we establish the bivariate mapping function of the C plane x 1 , x 2 into the G plane y 1 , y 2 by the sampling point of station and set a rectangle grid in the C plane.By the map function we mapping the grid in the C plane into the G plane the covariance structure of ground motion can be visualized through the density of grid in the G plane.
The main outline of this method can be described ad follows.
1 Calculate the DTW of the position y 1i , y 2i G plane .
2 Calculate the composition map by the MDS and we get the coordinates x 1i , x 2i in the C plane i 1, 2, . . ., N .
3 Establish the relationship of mapping by a bivariate functionf:  We define a roughness criterion and compute the bivariate function f f 1 , f 2 to minimize for specified smoothing parameter λ where R 2 is the domain of interest.
4 Define the point on the boundary line of the rectangle grid in the composition map and use the mapping function f to get the corresponding point in the G space: by the interpolating spline we can generate curve line and establish the grid.
Figure 4 summarizes the algorithm.

Results and Discussions
In this section we will do a preliminary analysis of seismograph array date by using the proposed nonparametric estimation method.This example manifests a somewhat extreme, but easily explained, form of nonstationary in the spatial covariance structure of the earthquake.The earthquake ground motion date is from the Event-40 recorded by SMART-1 array which is located in Lotung, Taiwan.The stations selected are located in a twodimensional surface array consisting of a center C00 and three concentric circles inner I,   5 .The station applied in this paper is I6, I9, I11, M10, M07, M03, O05, O06, O10, M06.The epicentral direction is 15 • to the north.The magnitude M L of the earthquake is 6.5, the location of earthquake source relative to the center station C00 is 22 km, and the depth of the source is 10 km.The recorded motion has a sampling rate of 100 Hz and a total of 1201 values.The computation uses an Ormsby filter to eliminate long periods from the acceleration recordings.The cut-off and rolloff frequencies are f C 0167 Hz and f T f T 0.03 0.197 Hz.By the aforementioned method in Section 2.1.1 we compute the observed proximity matrix δ ij DTW i, j of sampling in Table 2.By the aforementioned method in Section 2.1.2and the proximity matrix δ ij in Table 2, we can get the coordinates of the sampling station in the C plane in Figure 6.For example, stations records that are located close to each other in C plane are perceived as being similar such as I06 and I11.In contrast, stations positioned far away from each other indicate a large difference in perception, such as I06 and M10.
Figure 7 is the distance-disparities scares map.dispersions is nearly a straight line.It shows that this two-dimensional MDS solution accurately reflects the date in the matrix Δ δ ij , and the optimal metric has better linear relation than the classic metric method.
Estimated contours plots of the correlation coefficient are presented based on the strong motion recorded by SMART-1 array in Figure 8.What information conveys in Figure 8 is that, even in the alluvium alloy where seismic array located, the covariance structure of the earthquake is inhomogeneous in space choice or varies from location to location.The different choose of reference point may get different result, so the contours plot of the correlation coefficient is not a suitable way of representing the spatial variation of seismic ground motion.
Figure 8 is the figure of the method in this paper; it depicts the thin-plate spline mapping between the G plane and C plane representations of monitoring stations using the image of a rectangular of points located on C plane.Different from contour plot in Figure 9, it needs not to choose the reference point.The clearest feature in Figure 9 is that there is a relatively denser spacing of curves in the northeast; it is the epicentral direction.The area where the grid spacing of curves concentrated is the place where the correlative structure of earthquake changes significantly; the denser the grid the greater the variation and vice versa.It should be emphasized that this paper concentrates to illustrate how to apply our method to despite the corrective structure of earthquake ground motion, a more exact conclusion need more sufficiently dense spatial sampling.x (m)

Figure 9:
The grid for the spatial correlation structure of earthquake ground motion in geography space the x and y axes correspond to different geography direction .

Conclusions and Remarks
This study addresses the following topic: the description and estimation of spatial variation of seismic ground motions.A nonparametric method to the estimation and graphical depiction of the local spatial correlation for the earthquake ground motion is presented.The method paves a way how to use the multidimension scaling and body-fitted grid to represent the variation of ground motion.It can be a promising tool in detecting the correlation structure of spatial ground motion; further research includes the investigation of the influence, on the spatial correlation structure, of the seismologic parameters such as the epicenter the magnitude and location.Such covariance structure estimate by this method can be used to reveal the inhomogeneous feature of earthquake and help design large-span structure for example a suggested way to do so is to modify the traditional correlation model.By observing the density of the grid in the geography space we can introduce different weight coefficients at different locations.In this way we can produce more properly correlated motions in space for the design of engineering structure with character of spatial extension.The problem of how to establish the relationship between the grid by our method and weight coefficients of correlation model should be addressed in future work before the application of our methods to problems of artificial earthquake simulation.
The covariance estimates can also be used for some monitoring network design problems.For example, grid concentrated is the place more monitoring points stations should be placed and vice versa.

Figure 1 :Figure 2 :
Figure 1: The difference between DTW a and correlation coefficient b .

Figure 3 :
Figure 3: The mapping from the C plane a to G plane b .

Figure 4 :
Figure 4: Outline of the method.

Figure 7 aFigure 6 :
Figure 6: Coordinates of sampling station in the C plane.

Figure 7 :
Figure 7: Disparities versus MDS interstation distances for the seismograph array date.

Figure 8 :
Figure 8: Contour plot of correlation coefficient of sample record for Event 40: a take the I09 as the reference point and b take M05 as the reference point the unit of the coordinates is m . 10

Table 2 :
The DTW of sampling station.