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

Authors: Yingmin Li; Zheqian Wu; Huiguo Chen
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

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 ground motion in space is a challenge for researcher. A body-fitted grid-coordinates-based method 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.


---

## Body

## 1. 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).


## 2. 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.


### 2.1. Multidimensional Scaling Technique

#### 2.1.1. The Dynamic Time Warping Distance

For the accelerogram recordaij in n station (sample point), each has N independent realizations. The accelerogram record can be given in a (n×N) data matrix,(2.1)A=[a11a12⋯a1Na21a22⋯a2N⋮⋮⋱⋮an1an2⋯anN].
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:(2.2)ri,j=Sij(SiiSjj)1/2,
where
(2.3)Sij=1N-1[∑k=1Naikajk-1N(∑k=1Naik)(∑k=1Najk)],
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 point-to-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.LetA(a1,a2,…,an) and B(b1,b2,…,bm)be two series with the length of nand m, respectively, and ann×m matrixM can be defined to represent the point-to-point correspondence relationship between A and B, where the element Mij indicates the distance d(ai,bj)between ai and bj. Then the point-to-point alignment and matching relationship between X and Y can be represented by a time warping path: W=〈w1,w2,…,wk〉, max(m,n)≤K<m+n-1, where the element wk=(i,j) indicates the alignment and matching relationship between aiand bj. If a path is the lowest cost path between two series, the corresponding dynamic time warping distance is required to meet(2.4)DTW(A,B)=min⁡W{∑k=1Kdk,W=〈w1,w2,…,wK〉},dk=d(ai,bj) is the Minkowski distance of the k element of the path W and is denoted as follows:(2.5)d(A,B)=(∑i=1n|ai-bj|p)1/p,
whenp=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(2.6)DTW(〈〉,〈〉)=0,DTW(A,〈〉)=DTW(〈〉,B)=∞,DTW(A,B)=d(ai,bj)+min⁡⁡{DTW(A,B[2:-],DTW(A[2:-],B),DTW(A[2:-],B[2:-]),
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(ai,bj) indicates the distance between points ai and bj. 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:(2.7)r(i,j)=min⁡{r((i-1,j),r(i,j-1),r(i-1,j-1))}+d(ai,bj).
Any element r(i,j) in the accumulated matrix indicates the dynamic time warping distance between series A1:i andB1: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 Figure1.The difference between DTW (a) and correlation coefficient (b).
(a)(b)To demonstrate the different performance of the DTW and the correlation coefficientr, 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):A:
y=Asin(t) (0≤t≤2π),B:
y=2Asin(t) (0≤t≤2π),C:
y=Asin(t-0.5) (0≤t-0.5≤2π),D:
y=Asin(t) (0≤t≤11/6π).Table 1
DTW and correlation coefficient under four condition.
ConditionDTWrDTWrA-A01A–C00.87A-B2A1A–D0/The curve under four conditions.
(a)(b)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.


#### 2.1.2. Optimal Scaling

Multidimensional Scaling (MDS) is applied as a statistical technique to visualize dissimilarity data in this section.LetΔ=⌊δij⌋ and D=⌊dij⌋ 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 Δ, dij refers to the Euclidean space (composition map) between objects iandj. 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 dijcorrespond as closely as possible to the dissimilaritiesδij [18].Classical Metric Multidimensional Scaling is a basic form of MDS. Classical MDS is called metric methods because the relationship betweenδij and dij 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̂ anddij≈d̂ij=f(δij); it should be restricted by the monotonic constraint: pij≤plk⇒D(pij)≤D(plk).The coordinates in the distance function (in composition map) and the functionf which allows transforming the proximities into distances are estimated by minimizing the following badness of fit function (usually called stressor S-function in the context of MDS) or called optimal scaling. The stress function is given by(2.8)S=(∑∑(d̂ij2-dij2)2∑∑dij4)1/2.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 matrixD̂ by some iterative algorithms to make the S minimum. (x1i,x2i,…,xpi), where pis thedimensional of the final composition map (p=1,2,…,n,etc.).The method is summarized as follows.(1)
Select the initial matrix(x1i0,x2i0,…,xpi0).(2)
Calculate the distance matrixDk=[dij0] and elements of matrix Dk; dij0=∥xi,xj∥ is the Euclidean space.(3)
FromΔ=⌊δij⌋ and Dk=[dij0] use a monotonous regression Pool-Adjacent Violators (PAVs) [20] to get D̂k=[d̂ij].(4)
For the arbitraryi, j, k, l if dijk<dlkk⇒d̂ijk<d̂lkk then dijk=d̂ijk, otherwise d̂ijk=d̂lkk=(dijk+dlkk)/2.(5)
Calculate the new coordinate of node:(x1ik+1,x2ik+1,…,xqik+1,…,xpik+1)(2.9)xqik+1=xqik+θn-1∑j≠i(1-d̂ijkdijk)(xjk-xik),
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 Euclideandij.(7)
Calculate theS according to formula (2.8).(8)
IfS<ε, the calculation ends; otherwise return to step 3.Finally, the result(x1i,x2i,…,xpi)i (i=1,2) is the solution of MDS. If we chose p=2, the coordinate can be expressed as (x1i,x2i). We refer to the plane of geographical (y1i,y2i) of the sampling stations as G plane, and we refer to the plane of composition coordinates (x1i,x2i) as C plane (see Figure 3).The mapping from the C plane (a) to G plane (b).
(a)(b)


### 2.2. 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 coordinates-based 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(x1i,x2i) 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 (x1,x2) into the G plane (y1,y2) 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(y1i,y2i) (G plane).(2)
Calculate the composition map by the MDS and we get the coordinates(x1i,x2i) in the C plane (i=1,2,…,N).(3)
Establish the relationship of mapping by a bivariate functionf:
(2.10)(y1iy2i)=f(x1ix2i).
We define a roughness criterion
(2.11)Ji(f)=∫R2[(∂2fi∂x12)2+2(∂2fi∂x1∂x2)2+(∂2fi∂x22)2]dx1dx2
and compute the bivariate function f=(f1,f2) to minimize
(2.12)L=∑j=12∑i=1N(yij-fj(xij))2+λ(J1(f1)+J2(f2))
for specified smoothing parameter λ where R2 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 functionf to get the corresponding point in the G space: by the interpolating spline we can generate curve line and establish the grid.Figure4 summarizes the algorithm.Figure 4
Outline of the method.


## 2.1. Multidimensional Scaling Technique

### 2.1.1. The Dynamic Time Warping Distance

For the accelerogram recordaij in n station (sample point), each has N independent realizations. The accelerogram record can be given in a (n×N) data matrix,(2.1)A=[a11a12⋯a1Na21a22⋯a2N⋮⋮⋱⋮an1an2⋯anN].
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:(2.2)ri,j=Sij(SiiSjj)1/2,
where
(2.3)Sij=1N-1[∑k=1Naikajk-1N(∑k=1Naik)(∑k=1Najk)],
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 point-to-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.LetA(a1,a2,…,an) and B(b1,b2,…,bm)be two series with the length of nand m, respectively, and ann×m matrixM can be defined to represent the point-to-point correspondence relationship between A and B, where the element Mij indicates the distance d(ai,bj)between ai and bj. Then the point-to-point alignment and matching relationship between X and Y can be represented by a time warping path: W=〈w1,w2,…,wk〉, max(m,n)≤K<m+n-1, where the element wk=(i,j) indicates the alignment and matching relationship between aiand bj. If a path is the lowest cost path between two series, the corresponding dynamic time warping distance is required to meet(2.4)DTW(A,B)=min⁡W{∑k=1Kdk,W=〈w1,w2,…,wK〉},dk=d(ai,bj) is the Minkowski distance of the k element of the path W and is denoted as follows:(2.5)d(A,B)=(∑i=1n|ai-bj|p)1/p,
whenp=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(2.6)DTW(〈〉,〈〉)=0,DTW(A,〈〉)=DTW(〈〉,B)=∞,DTW(A,B)=d(ai,bj)+min⁡⁡{DTW(A,B[2:-],DTW(A[2:-],B),DTW(A[2:-],B[2:-]),
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(ai,bj) indicates the distance between points ai and bj. 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:(2.7)r(i,j)=min⁡{r((i-1,j),r(i,j-1),r(i-1,j-1))}+d(ai,bj).
Any element r(i,j) in the accumulated matrix indicates the dynamic time warping distance between series A1:i andB1: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 Figure1.The difference between DTW (a) and correlation coefficient (b).
(a)(b)To demonstrate the different performance of the DTW and the correlation coefficientr, 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):A:
y=Asin(t) (0≤t≤2π),B:
y=2Asin(t) (0≤t≤2π),C:
y=Asin(t-0.5) (0≤t-0.5≤2π),D:
y=Asin(t) (0≤t≤11/6π).Table 1
DTW and correlation coefficient under four condition.
ConditionDTWrDTWrA-A01A–C00.87A-B2A1A–D0/The curve under four conditions.
(a)(b)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.


### 2.1.2. Optimal Scaling

Multidimensional Scaling (MDS) is applied as a statistical technique to visualize dissimilarity data in this section.LetΔ=⌊δij⌋ and D=⌊dij⌋ 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 Δ, dij refers to the Euclidean space (composition map) between objects iandj. 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 dijcorrespond as closely as possible to the dissimilaritiesδij [18].Classical Metric Multidimensional Scaling is a basic form of MDS. Classical MDS is called metric methods because the relationship betweenδij and dij 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̂ anddij≈d̂ij=f(δij); it should be restricted by the monotonic constraint: pij≤plk⇒D(pij)≤D(plk).The coordinates in the distance function (in composition map) and the functionf which allows transforming the proximities into distances are estimated by minimizing the following badness of fit function (usually called stressor S-function in the context of MDS) or called optimal scaling. The stress function is given by(2.8)S=(∑∑(d̂ij2-dij2)2∑∑dij4)1/2.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 matrixD̂ by some iterative algorithms to make the S minimum. (x1i,x2i,…,xpi), where pis thedimensional of the final composition map (p=1,2,…,n,etc.).The method is summarized as follows.(1)
Select the initial matrix(x1i0,x2i0,…,xpi0).(2)
Calculate the distance matrixDk=[dij0] and elements of matrix Dk; dij0=∥xi,xj∥ is the Euclidean space.(3)
FromΔ=⌊δij⌋ and Dk=[dij0] use a monotonous regression Pool-Adjacent Violators (PAVs) [20] to get D̂k=[d̂ij].(4)
For the arbitraryi, j, k, l if dijk<dlkk⇒d̂ijk<d̂lkk then dijk=d̂ijk, otherwise d̂ijk=d̂lkk=(dijk+dlkk)/2.(5)
Calculate the new coordinate of node:(x1ik+1,x2ik+1,…,xqik+1,…,xpik+1)(2.9)xqik+1=xqik+θn-1∑j≠i(1-d̂ijkdijk)(xjk-xik),
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 Euclideandij.(7)
Calculate theS according to formula (2.8).(8)
IfS<ε, the calculation ends; otherwise return to step 3.Finally, the result(x1i,x2i,…,xpi)i (i=1,2) is the solution of MDS. If we chose p=2, the coordinate can be expressed as (x1i,x2i). We refer to the plane of geographical (y1i,y2i) of the sampling stations as G plane, and we refer to the plane of composition coordinates (x1i,x2i) as C plane (see Figure 3).The mapping from the C plane (a) to G plane (b).
(a)(b)


## 2.1.1. The Dynamic Time Warping Distance

For the accelerogram recordaij in n station (sample point), each has N independent realizations. The accelerogram record can be given in a (n×N) data matrix,(2.1)A=[a11a12⋯a1Na21a22⋯a2N⋮⋮⋱⋮an1an2⋯anN].
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:(2.2)ri,j=Sij(SiiSjj)1/2,
where
(2.3)Sij=1N-1[∑k=1Naikajk-1N(∑k=1Naik)(∑k=1Najk)],
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 point-to-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.LetA(a1,a2,…,an) and B(b1,b2,…,bm)be two series with the length of nand m, respectively, and ann×m matrixM can be defined to represent the point-to-point correspondence relationship between A and B, where the element Mij indicates the distance d(ai,bj)between ai and bj. Then the point-to-point alignment and matching relationship between X and Y can be represented by a time warping path: W=〈w1,w2,…,wk〉, max(m,n)≤K<m+n-1, where the element wk=(i,j) indicates the alignment and matching relationship between aiand bj. If a path is the lowest cost path between two series, the corresponding dynamic time warping distance is required to meet(2.4)DTW(A,B)=min⁡W{∑k=1Kdk,W=〈w1,w2,…,wK〉},dk=d(ai,bj) is the Minkowski distance of the k element of the path W and is denoted as follows:(2.5)d(A,B)=(∑i=1n|ai-bj|p)1/p,
whenp=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(2.6)DTW(〈〉,〈〉)=0,DTW(A,〈〉)=DTW(〈〉,B)=∞,DTW(A,B)=d(ai,bj)+min⁡⁡{DTW(A,B[2:-],DTW(A[2:-],B),DTW(A[2:-],B[2:-]),
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(ai,bj) indicates the distance between points ai and bj. 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:(2.7)r(i,j)=min⁡{r((i-1,j),r(i,j-1),r(i-1,j-1))}+d(ai,bj).
Any element r(i,j) in the accumulated matrix indicates the dynamic time warping distance between series A1:i andB1: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 Figure1.The difference between DTW (a) and correlation coefficient (b).
(a)(b)To demonstrate the different performance of the DTW and the correlation coefficientr, 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):A:
y=Asin(t) (0≤t≤2π),B:
y=2Asin(t) (0≤t≤2π),C:
y=Asin(t-0.5) (0≤t-0.5≤2π),D:
y=Asin(t) (0≤t≤11/6π).Table 1
DTW and correlation coefficient under four condition.
ConditionDTWrDTWrA-A01A–C00.87A-B2A1A–D0/The curve under four conditions.
(a)(b)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.


## 2.1.2. Optimal Scaling

Multidimensional Scaling (MDS) is applied as a statistical technique to visualize dissimilarity data in this section.LetΔ=⌊δij⌋ and D=⌊dij⌋ 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 Δ, dij refers to the Euclidean space (composition map) between objects iandj. 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 dijcorrespond as closely as possible to the dissimilaritiesδij [18].Classical Metric Multidimensional Scaling is a basic form of MDS. Classical MDS is called metric methods because the relationship betweenδij and dij 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̂ anddij≈d̂ij=f(δij); it should be restricted by the monotonic constraint: pij≤plk⇒D(pij)≤D(plk).The coordinates in the distance function (in composition map) and the functionf which allows transforming the proximities into distances are estimated by minimizing the following badness of fit function (usually called stressor S-function in the context of MDS) or called optimal scaling. The stress function is given by(2.8)S=(∑∑(d̂ij2-dij2)2∑∑dij4)1/2.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 matrixD̂ by some iterative algorithms to make the S minimum. (x1i,x2i,…,xpi), where pis thedimensional of the final composition map (p=1,2,…,n,etc.).The method is summarized as follows.(1)
Select the initial matrix(x1i0,x2i0,…,xpi0).(2)
Calculate the distance matrixDk=[dij0] and elements of matrix Dk; dij0=∥xi,xj∥ is the Euclidean space.(3)
FromΔ=⌊δij⌋ and Dk=[dij0] use a monotonous regression Pool-Adjacent Violators (PAVs) [20] to get D̂k=[d̂ij].(4)
For the arbitraryi, j, k, l if dijk<dlkk⇒d̂ijk<d̂lkk then dijk=d̂ijk, otherwise d̂ijk=d̂lkk=(dijk+dlkk)/2.(5)
Calculate the new coordinate of node:(x1ik+1,x2ik+1,…,xqik+1,…,xpik+1)(2.9)xqik+1=xqik+θn-1∑j≠i(1-d̂ijkdijk)(xjk-xik),
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 Euclideandij.(7)
Calculate theS according to formula (2.8).(8)
IfS<ε, the calculation ends; otherwise return to step 3.Finally, the result(x1i,x2i,…,xpi)i (i=1,2) is the solution of MDS. If we chose p=2, the coordinate can be expressed as (x1i,x2i). We refer to the plane of geographical (y1i,y2i) of the sampling stations as G plane, and we refer to the plane of composition coordinates (x1i,x2i) as C plane (see Figure 3).The mapping from the C plane (a) to G plane (b).
(a)(b)


## 2.2. 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 coordinates-based 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(x1i,x2i) 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 (x1,x2) into the G plane (y1,y2) 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(y1i,y2i) (G plane).(2)
Calculate the composition map by the MDS and we get the coordinates(x1i,x2i) in the C plane (i=1,2,…,N).(3)
Establish the relationship of mapping by a bivariate functionf:
(2.10)(y1iy2i)=f(x1ix2i).
We define a roughness criterion
(2.11)Ji(f)=∫R2[(∂2fi∂x12)2+2(∂2fi∂x1∂x2)2+(∂2fi∂x22)2]dx1dx2
and compute the bivariate function f=(f1,f2) to minimize
(2.12)L=∑j=12∑i=1N(yij-fj(xij))2+λ(J1(f1)+J2(f2))
for specified smoothing parameter λ where R2 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 functionf to get the corresponding point in the G space: by the interpolating spline we can generate curve line and establish the grid.Figure4 summarizes the algorithm.Figure 4
Outline of the method.


## 3. 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 two-dimensional surface array consisting of a center C00 and three concentric circles (inner I, middle M, and outer O) with radii of 200, 1000, and 2000 m, respectively, each with 12 strong-motion seismographs having common time base, (see Figure5). 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 (ML) 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 roll-off frequencies are fC=0167 Hz and fT=fT+0.03=0.197 Hz.Figure 5
The array of SMART-1.By the aforementioned method in Section2.1.1 we compute the observed proximity matrix [δij]=DTW(i,j) of sampling in Table 2.Table 2
The DTW of sampling station.
I6I9I11M10M07M03O05O06O10O11M06I600.80360.86421.04371.07591.18911.12581.00091.29651.30181.1234I90.803600.57790.9391.1970.99751.07811.01361.19921.3330.9159I110.86420.577900.77551.28551.14391.07051.05461.37851.38171.0058M101.04370.9390.775501.20881.27640.99471.00471.4081.50891.1313M071.07591.1971.28551.208801.2260.91121.05941.10461.50461.0792M031.18910.99751.14391.27641.22601.26361.34580.87581.34131.1292O051.12581.07811.07050.99470.91121.263601.41431.41361.48471.414O061.00091.01361.05461.00471.05941.34581.414301.26221.42571.0249O101.29651.19921.37851.4081.10460.87581.41361.262201.1911.2415O111.30181.3331.38171.50891.50461.34131.48471.42571.19101.3824M061.12340.91591.00581.13131.07921.12921.4141.02491.24151.38240By the aforementioned method in Section2.1.2 and 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 6
Coordinates of sampling station in the C plane.Figure7 is the distance-disparities scares map. Figure 7(a) uses by the classic metric methods and Figure 7(b) is the result of optimal metric method; from Figure 7(b) it can be observed that the scatter plot of distances dij in the C plane versus the observed spatial 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.Disparities versus MDS interstation distances for the seismograph array date.
(a)(b)Estimated contours plots of the correlation coefficient are presented based on the strong motion recorded by SMART-1 array in Figure8. 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.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).
(a)(b)Figure8 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.Figure 9
The grid for the spatial correlation structure of earthquake ground motion in geography space (thex and y axes correspond to different geography direction).


## 4. 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.



---

*Source: 420309-2012-05-14.xml*