 Research
 Open Access
 Published:
Investigating the spatial correlations in univariate random fields of peak ground velocity and peak ground displacement considering anisotropy
Geoenvironmental Disasters volume 8, Article number: 24 (2021)
Abstract
The results of seismic risk assessment of spatially distributed infrastructure systems are significantly influenced by spatial correlation of earthquake intensity measures (IM). The assumption of isotropy is a basis for most of the existing correlation models of earthquake IMs. In this study, the isotropy assumption of intraevent residuals of peak ground velocity (PGV) and peak ground displacement (PGD) is investigated by implementing a nonparametric statistical test. Using recorded IMs of 9 earthquakes, it is concluded that there is not sufficient evidence to support the assumption of isotropy in general, and the set of intraevent residuals of PGV and PGD should be considered as the realization of anisotropic random fields. Investigations show that the anisotropy properties of intraevent residuals of PGV and PGD are related to anisotropy properties of local soil characteristics indicated by average shear wave velocity of soil profile from the 30 m depth to the surface (V_{s30}). Finally, predictive models are proposed based on obtained results in order to simulate the correlated univariate random fields of PGV and PGD considering anisotropy.
Introduction
In order to evaluate the seismic risk of spatially distributed systems, we need spatially correlated ground motion intensity measures (IMs). Previous research works show that neglecting the spatial correlation of IMs causes significant miscalculations in seismic risk assessment of infrastructure systems or portfolio of buildings (Lee and Kiremidjian 2007; Park et al. 2007; Weatherill et al. 2015; Garakaninezhad and Bastami 2017). Park et al. (2007) by considering two types of large and small portfolios of buildings, concluded that underestimating or neglecting spatial correlations of earthquake IMs result in overestimating frequent losses and underestimating rare losses. In recent years, there has been a lot of interest in studying the spatial correlations of IMs, and different models have been proposed in this regard that among them the models of Boore et al., (2003), Du and Wang (2013), Goda and Hong (2008), Goda and Atkinson (2009), Jayaram and Baker (2009), Wang and Du (2013), Wang and Takada (2005), Markhvida et al. (2018), Abbasnejadfard et al. (2020), Garakaninezhad et al. (2017), Garakaninezhad and Bastami (2019), Zafarani et al.(2021) can be mentioned. Boore et al. (2003) proposed a spatial correlation function of peak ground acceleration (PGA) using the recorded PGA values of the 1994 Northridge earthquake. Using California and ChiChi ground motion records, K. Goda and Hong (2008) investigated the spatial correlation of PGA and spectral accelerations (SAs) at three periods of 0.3, 1, and 3 s. The researchers concluded in the mentioned study that the separation distance of interested locations and the ground motion parameters are influential factors in the spatial correlation of IMs' intraevent residuals. Goda and Atkinson (2009) used the ground motion database of Japan and studied the spatial correlation of univariate random fields of intraevent residuals of different SAs with different periods. They proposed an exponential decay function of lag distance as the correlations model. Wang and Takada (2005), showed that the correlation range of PGV values of different earthquakes recorded in Japan and Taiwan is between 60 to 120 km. In addition, the spatial correlation characteristics of IMs' intraevent residuals are found to be unaffected by various ground motion prediction equations (GMPE). Jayaram and Baker (2009) used geostatistical tools along with the records of seven earthquakes to assess the univariate spatial correlation of different SAs. The authors of the mentioned study also used the 1994 Northridge and the 1999 ChiChi earthquakes to explore the isotropy of spatial correlations of IMs and concluded that the isotropy is a reasonable assumption for both earthquakes. Du and Wang (2013), using the recorded data of earthquakes in California, Taiwan, and Japan, investigated the effect of regional site condition on spatial correlation of Arias intensity (I_{a}), cumulative absolute velocity (CAV), and SAs at different periods, and proposed predictive models of spatial correlation range of the mentioned IMs. Wang and Du (2013) used the linear model of coregionalization (LMC) method to investigate the spatial crosscorrelation of the vectors of IMs (PGA, PGV, and I_{a} as the first set and SAs at different periods as the second set of IMs). The results show that the spatial characteristics of the regional site condition significantly affect the spatial crosscorrelation of IMs. To simulate the multivariate random field of SAs at different periods, Markhvida et al. (2018) utilized principal component analysis (PCA) method and proposed a predictive model in this regard.
In all of the previously mentioned research works, the isotropy and stationarity assumptions are considered as the basis for spatial correlation investigations. Abbasnejadfard et al. (2021) by utilizing different spatial correlation models, including isotropic, anisotropic and uncorrelated models, depicted that “ignorance of anisotropy of spatial correlations of earthquake IMs causes unrealistic loss estimation and leads to inaccurate resilience assessment of spatially distributed assets and systems”. Garakaninezhad and Bastami (2017) used the statistical test introduced by Bowman and Crujeiras (2013) to investigate the isotropy hypothesis of intraevent residual of PGA and SAs and concluded that the hypothesis of isotropy is not valid in general. Besides, the properties of anisotropic spatial correlations of intraevent residues of PGA and SAs are evaluated. In this context, Abbasnejadfard et al. (2020) investigated the spatial correlations of the multivariate random field of intraevent residuals of multiple earthquake IMs. The authors of the mentioned study considered two sets of IMs (PGA, PGV, and PGD as the first set, and SAs at three periods of 0.5, 1, and 2 s as the second set) and utilized the latent dimensions (LD) method proposed by Apanasovich and Genton (2010) to determine predictive models of multivariate random fields of IMs considering anisotropy.
Visual inference of the results of Abbasnejadfard et al., (2020) shows that the random field of intraevent residuals of PGV and PGD are anisotropic. Although the visual inference of spatial data is an acceptable approach, this method generally does not include the intrinsic uncertainty of sample data and is open to interpretations and can be misleading (Weller et al. 2016). In this regard, the isotropy assumption of the PGV and PGD is examined in the current research work using a statistical test of isotropy. Moreover, the effect of regional site conditions and sourcetosite distance of record stations on the spatial correlation of the mentioned IMs are investigated. It should be noted that some of the infrastructure systems (specifically the burred networks) are influenced by only PGV or PGD. Although the results of Abbasnejadfard et al., (2020) can be utilized to model multivariate random fields of earthquake IMs, including PGV and PGD, it will be more helpful to use a simplified univariate anisotropic model to estimate correlated IMs where a single IM is of interest. In this regard, predictive models that consider the anisotropy of spatial correlations are proposed to simulate the univariate random field of intraevent residual of PGV and PGD for practical use.
Methods
Selected ground motions
In this study, the recorded data of nine earthquakes listed in Table 1 are used to investigate the isotropy assumption of spatial correlations of intraevent residuals of PGV and PGD. It should be mentioned in this section that selecting earthquake data for use in spatial statistics investigations is subject to various limitations, especially when anisotropy properties is of interest. Indeed, we require earthquake events recorded at a large number of record stations. Furthermore, the record stations should be appropriately distributed throughout the study region to comply with spatial statistics requirements. The reader is referred to Chiles and Delfiner (2012) for more information. Most of the existing recorded earthquake events are unable to be chosen due to the mentioned restrictions.
The geological and seismological characteristics of the selected records are obtained from the Pacific Earthquake Engineering Research (PEER) ground motion database (PEER NGA Database 2019). The selected earthquakes have Mw greater than 5. Figure 1 shows the distribution of collected data in terms of magnitude and rupture distance. The rupture distance herein is defined as the closest distance between the interested location and rupture plane. The utilized intensity measures are the geometric averages of the two orthogonal horizontal components recorded in the record stations.
It is worth mentioning that the doubleintegration of padstripped accelerations may cause distortions in some of the PGD values acquired from the PEERNGA database, making the estimation of PGD values inaccurate in some situations. However, because the results derived from PGD values are consistent with those of PGA (Garakaninezhad and Bastami 2017) and PGV, the investigations of PGD values are also included in the current research work.
The Campbell and Bozorgnia (2014) (CB14) and Campbell and Bozorgnia (2008) (CB08) GMPEs are utilized to estimate the intraevent residuals of PGV and PGD, respectively. The selected earthquake events have been used in the development of CB14 GMPE and recorded in the sufficient number of record stations so that they are proper to be used in the context of spatial statistics.
Modeling spatial correlations
Generally, ground motion IMs at distributed points constitute a twodimensional random field that can be modeled as:
where \(\overline{{Y_{ij} }}\) is the median value of IM that is estimated by GMPE, \(\delta_{ij}\) is the small scale variation or residual and \(Y_{ij}\) is the IM at site i, given jth earthquake event that is composed of interevent (\(\eta_{j}\)) and intraevent (\(\varepsilon_{ij}\)) residuals (see Eq. (2)). \(\eta_{j}\) represents the eventtoevent random variability of IM and consequently is constant for all sites in a single earthquake event j. The value of \(\eta_{j}\) is treated as a normal random variable with a mean of zero and a standard deviation of \(\tau\) (Jayaram and Baker 2009; Wang and Du 2013). \(\varepsilon_{ij}\) represents the sitetosite variability of IM. It is assumed that the vector of \({{\varvec{\upvarepsilon}}}_{{\mathbf{j}}} = (\varepsilon_{1j} ,\varepsilon_{2j} ,...,\varepsilon_{nj} )\) for n sites follows a multivariate normal distribution with a mean of zero and a standard deviation of \(\sigma\) (Jayaram and Baker 2008).
According to the previous studies, the standard deviation of intraevent residuals does not depend on the location of interested sites (Jayaram and Baker 2009; Du and Wang 2013; Garakaninezhad and Bastami 2017). Considering this assumption and substituting Eq. (2) into the Eq. (1), the normalized intraevent residual in an earthquake event j at site i can be calculated as:
In the current research work, geostatistical tools are used to investigate the spatial correlations of the set of normalized intraevent residuals of \(\varepsilon^{\prime}_{ij}\) at multiple sites. In this regard, the variogram function is defined as the variance of the difference of intraevent residuals of two points as presented in the Eq. (4) (Cressie 1993).
where \(\varepsilon^{\prime}_{i}\) and \(\varepsilon^{\prime}_{k}\) are the values of normalized intraevent residuals at two points i and k with coordinates of s_{i} and s_{k,} respectively. The random field of intraevent residuals is assumed to be secondordered stationary (Jayaram and Baker 2009), which implies that the random field's mean value is constant at all points, and the covariance of random field is only a function of the relative position vector (lag distance) of points (see Eqs. (5) and (6)) (Cressie 1993).
Considering the assumption of second orderedstationary, the Eq. (4) can be rewritten as:
where \(C({\mathbf{0}})\) is the variance of the secondordered stationary random field and \(\gamma ({\mathbf{h}})\) is known as semivariogram. In the case of normalized intraevent residuals, the m value of Eq. (5) equals zero, and \(C({\mathbf{0}})\) in Eq. (7) equals to 1. Consequently, the Eq. (7) can be modified in the following form:
In order to define the correlated values of \(\varepsilon^{\prime}_{ij}\) in n locations, a \(n \times n\) correlation matrix should be defined as:
where \(\rho ({\mathbf{h}}_{ik} )\) is the correlation coefficient (correlogram in terms of geostatistics) between values in points i and k which has the lag distance of \({\mathbf{h}}_{ik}\). The correlogram can be calculated as (Chilès and Delfiner 2012):
Substituting the Eq. (7) into the Eq. (10), the correlogram can be obtained as:
For a set of recorded data, the empirical semivariogram for lag vector (h) can be calculated as (Chilès and Delfiner 2012):
where \(N_{{\mathbf{h}}}\) represents all pairs of points that have the relative separation vector of h. A random field is called isotropic if its variogram is only a function of the norm of lag vector (\(h = {\mathbf{h}}\)) and is independent of the direction of the lag vector. In the case of an isotropic random field, the Eq. (12) can be represented as:
which \(N_{h}\) is the number of all pairs of points that have the relative separation distance of h. In order to minimize the unfavorable effects of outliers on estimating the variogram, Cressie and Hawkins (1980) introduced a robust estimator of variogram as provided in the Eq. (14). This estimator is utilized in the current study.
Different valid semivariogram models can be addressed in the literature to fit a set of empirical semivariogram values. A valid semivariogram model is a model that leads us to a positivedefinite covariance matrix. The exponential model ((15)) is shown to be the best fitting model for empirical semivariograms of intraevent residuals of ground motion IMs.
In the Eq. (15), a is the sill of the semivariogram model and shows the variance of data, b is the range of the semivariogram model. The \(\varepsilon^{\prime}(s_{\beta } )\) and \(\varepsilon^{\prime}(s_{\alpha } )\) are considered uncorrelated beyond the range value. For the random field of normalized intraevent residuals, the data variance and consequently the sill of semivariogram equal to 1 and correlogram can be obtained as:
Isotropy of the semivariograms of the intraevent residuals
If the semivariogram of a secondordered stationary random field does not vary with the direction, the random field is said isotropic (Chilès and Delfiner 2012). In contrast, when the semivariogram is the function of both direction and norm of lag distance h, the random field is anisotropic. Geometric and zonal anisotropy are two typical forms of anisotropy (Chilès and Delfiner 2012). In geometric anisotropy, the directional semivariograms in different directions have the same sill value while their ranges vary with direction. The ratio of maximum range to the minimum range in different directions is called anisotropy ratio (AR), and the orientation of the maximum range of the semivariogram is called anisotropy angle (AA). In the cases where both the sill and range of the semivariogram change by direction, the random field is called zonal anisotropic.
Generally, test methods of isotropy can be classified as graphical methods and statistical hypothesis test methods. Graphical methods are based on obtaining directional semivariograms in different directions and comparing their characteristic. This method is used by Jayaram and Baker (2009) to examine the assumption of isotropy of intraevent residuals of SAs using the 1994 Northridge and the 1999 ChiChi earthquakes, and it is concluded that the semivariograms of selected earthquakes are isotropic. On the other hand, (Abbasnejadfard et al. 2020) utilized recorded data of ten earthquakes and presented the twodimensional covariance and crosscovariance functions of different earthquake IMs (including PGV and PGD). Visually, it can be concluded from the findings of the mentioned study that the spatial correlations of earthquake IMs are anisotropic in some cases. Although the graphical investigation of isotropy is an acceptable approach (Goovaerts 1997), because of the visual constitute of this method, it generally does not include the intrinsic uncertainty of sample data and is open to interpretations and can be misleading (Weller et al. 2016). To avoid these problems, methods that are based on the statistical test of the hypothesis have been developed. These methods are classified as parametric and nonparametric tests that are reviewed and discussed in Weller et al. (2016). Based on the recommendations of Weller et al. (2016), the nonparametric test, which is introduced by Bowman and Crujeiras (2013), is utilized in the current study. This method compares the variograms which are estimated incorporating distance and angle (\(\hat{\gamma }_{1}^{*}\)) and the variograms which are estimated incorporating distance only (\(\hat{\gamma }_{0}^{*}\)) in order to probe the assumption of the isotropy. Considering a vector of mean values of distance bins (b) over a twodimensional space \((h,\theta )\), estimated values of variograms can be obtained using the Bsplines method introduced by Eilers and Marx, (1996) (see Bowman and Crujeiras (2013) for more details). The variograms values at bin means of b derived from this smoothing technique are demonstrated as \(\hat{\gamma }_{1}^{*} = {\mathbf{S}}_{1} {\mathbf{b}}\) and \(\hat{\gamma }_{0}^{*} = {\mathbf{S}}_{0} {\mathbf{b}}\) where S_{1} is smoothing matrix considering the variability with both distance and angle and S_{0} is smoothing matrix considering the variability with distance only. The difference between two smoothed variograms,\(\hat{\gamma }_{1}^{*}  \hat{\gamma }_{0}^{*} = ({\mathbf{S}}_{1}  {\mathbf{S}}_{0} ){\mathbf{b}}\), is a zeromean variable under the assumption of isotropy. The test statistic for assessing the isotropy is:
where \({\mathbf{V}}_{0}\) is covariance matrix for the isotropic case and \({\mathbf{S}} = ({\mathbf{S}}_{1}  {\mathbf{S}}_{0} )\). Smoothed matrices imply that the test statistic (17) does not follow a simple \(\chi^{2}\) distribution. However, singleshifted and scaled \(\chi^{2}\) distribution is a good approximation of test statistics (Bowman and Azzalini 1997), and the pvalue of interest can be approximated accurately. The pvalue is the measure that quantifies the strength of the evidence against the null hypothesis of isotropy. So the case of p < 0.05 provides sufficient evidence that the null hypothesis is not correct. By implementing the abovementioned nonparametric test of isotropy, Garakaninezhad and Bastami (2017) concluded that the hypothesis of isotropy is not valid in general for the intraevent residuals of PGA. and SAs.
The isotropy of intraevent residuals of PGV and PGD is investigated in the following sections using the isotropy test of Bowman and Crujeiras (2013).
Results and discussion
Investigating the isotropy of intraevent residuals of PGV
In the presented study, the package sm (Bowman and Azzalini 2014) for the R software (R Core Team 2014) is used to examine the isotropy assumption of intraevent residuals of the PGV that recorded in earthquakes presented in Table 1. The pvalues of the isotropy test are less than 0.001 except for Northridge and Parkfield earthquakes which, are 0.307 and 0.028, respectively. It can be concluded from these values that the null hypothesis (isotropy of the intraevent residuals of the PGV) is rejected at a significance level of 1% for all cases except for Northridge and Parkfield earthquakes.
The smoothed variogram of residuals of PGV in the fourthroot scale (concerning the distance and angle) for the Anza earthquake is shown in Fig. 2. In this figure, the distance can be measured as the length between each point and the center point in centimeter, and the angle can be measured as the azimuth of points’ radius from the center point. The relative distance between sampling points is a real positive number, and consequently, the smoothed variogram should be presented over (0, π], but because of simplicity, the smoothed variogram is showed over (0, 2π] by reflection (Bowman and Crujeiras 2013). The variable values of smoothed variogram relative to the angle near the origin represent the anisotropy. Contours in this figure represent the distance between the variogram as a function of distance only and the variogram as a function of both angle and distance in units of standard error (see Bowman and Crujeiras 2013 for more details). The smoothed variograms of residuals of PGV at record stations of other earthquakes are shown in Fig. 3.
A predictive model for intraevent residuals of PGV considering anisotropy
According to the previous section, it can be inferred that the set of intraevent residuals of PGV in an earthquake event should be considered as a realization of an anisotropic random field, which means that the correlation between the residuals values of points s and s + h is a function of both norm and direction of vector h. The main objective of this section is to introduce a predictive model of intraevent residuals of PGV considering anisotropy.
Anisotropy ratio (AR) (the ratio of widest to the shortest range of semivariograms in different directions), anisotropy angle (AA) (orientation of the widest range of semivariograms), and the value of the widest range of semivariograms (b_{max}) are the main parameters that are needed to model an anisotropic random field. In this study, the directional variogram method is used to determine the characteristics of variograms in different directions and obtain AR, AA, and b_{max} of selected earthquake events.
Directional variogram can be calculated based on Eqs. (13) or (14) except that the data which falls within the azimuth limits of interest will be considered in the calculation of variograms. The azimuth of the direction of interest, azimuth tolerance, and sampling bin separation are the main parameters to determine the directional variogram. In this study, azimuth directions of 0°, 45°, 90°, 135°, and azimuth tolerance of 45° are considered. Furthermore, a bin separation distance of 6 ~ 8 km is used. Figure 4 shows experimental and fitted directional semivariogram of residuals of PGV for the 2008 Iwate, 2004 Niigata, and 2004 Parkfield earthquakes. The exponential semivariogram model is fitted to the experimental data. This model is the most suitable model for data and is consistent with previous studies. Various methods are recommended to fit curves, such as the least square method, the weighted least square method, or the manual fitting method. The choice of each of these methods depends on the configuration of the available data. If there is enough experimental data at different separation lags (h), and semivariogram data has no outliers, the least square method is the best suitable method. Generally, the values of semivariograms at short separation lags are obtained based on the limited number of data and may have considerable deviations. In the case which there is not enough experimental data at different separation lags (h) or semivariogram data has outliers, assigning equal weights to the semivariogram data at all separation lags will result in an inaccurate fitted model, and consequently using the weighted least square method is recommended in this case (Jayaram and Baker 2009). In this study, the value of \(N(j)/\left( {h(j)} \right)^{2}\) is used as weighting factors for the semivariogram value of separation bin j where \(N(j)\) is the count of data pairs and \(h(j)\) is the separation distance corresponding to bin j.
As shown in Fig. 4, no semivariogram model can be fitted to experimental directional semivariogram values of residuals of PGV at Parkfield earthquake due to the large scatter. This mode is also observed in the Northridge earthquake and implies that the data are uncorrelated or the correlation range is very short.
A summary of ranges of directional semivariograms of intraevent residuals of PGV of selected earthquakes is presented in the rose diagrams of Fig. 5. Furthermore, this figure shows the rose diagrams of the ranges of directional semivariograms of the normalized V_{s30} values at the record stations of these earthquakes. The term “normalized V_{s30}” refers to scaling the set of V_{s30} values so that the scaled values have the zero mean and unit standard deviation. As shown in Fig. 5 no diagram is presented for the PGV and the normalized V_{s30} values of Parkfield and Northridge earthquakes because no valid directional semivariogram function can be fitted to data of these earthquakes. Figure 5 also presents the AR and b_{max} of the PGV residuals and normalized V_{s30} values of each earthquake event.
According to Fig. 5, almost for all cases, the AA of the residuals of PGV and normalized V_{s30} values are consistent with each other. Thus, it can be inferred that if a directional correlation exists between the V_{s30} values, there would be a directional correlation between the intraevent residuals of PGV. In addition, the AA of the random field of the PGV residuals would be compatible with the AA of the random field of the normalized V_{s30} values.
Figure 6 depicts the AR of residuals of PGV (AR_{PGV}) versus AR of the normalized V_{s30} values (AR_{Vs30}) of different earthquakes. As shown in this figure, there is a positive correlation between the AR_{PGV} and AR_{Vs30} with a correlation coefficient of 0.86. The linear relationship between AR_{PGV} and AR_{Vs30} is obtained by linear regression analysis as:
Investigations show that b_{max} has a good correlation with the average sourcetosite distance of record stations. Various sourcetosite distance definitions can be addressed in literature like epicentral distance (R_{EPI}), hypocentral distance (R_{HYP}), JoinerBoor distance (R_{.J.B.}) and minimum distance from recording site to rupture area (R_{RUP}). Moreover, because of the spatial essence of the problem, different definitions of average value can be used, like arithmetic mean ((19)), geometric mean ((20)), harmonic mean ((21)), and quadratic mean ((22)). The use of harmonic and quadratic mean values can represent the pattern of locations of points. Moreover, the harmonic mean can reduce the effect of outlier points on the total mean value.
Table 2 presents the coefficient of determination (R^{2}) and the pvalue of Ftest significance of linear regression model between b_{max} of intraevent residuals of PGV as the dependent variable and different definitions of the average sourcetosite distance of record stations as the independent variable.
As shown in this table, the linear regression based on the harmonic mean of R_{RUP} (\(\overline{{R_{RUP}^{hr} }}\)) fits the data better than other definitions of average distances. The relationship between b_{maxPGV} and \(\overline{{R_{RUP}^{hr} }}\) is shown in Fig. 7 and Eq. (23). The regression model between b_{maxPGV} and \(\overline{{R_{RUP}^{hr} }}\) has a coefficient of determination of 0.57 and is statistically significant in the significance level of 5%.
In Eq. (23), b_{maxPGV} is the maximum range of directional semivariograms of intraevent residuals of PGV and \(\overline{{R_{RUP}^{hr} }}\) is the harmonic mean of R_{RUP} of the set of interesting points.
Investigating the isotropy of intraevent residuals of PGD
The assumption of isotropy of intraevent residuals of PGD is tested in this section. The CB08 GMPE is used to calculate the residuals of PGD. Since some earthquakes considered in the present study are not used in developing the CB08 GMPE, the residuals of PGD of these earthquakes have trend against the predictive parameters like R_{rup} and V_{s30}. This issue is also mentioned in Du and Wang (2013) and Wang and Du (2013) for other intensity measures. Wang and Du (2013) used the following equation to remove the trend in residuals.
In Eq. (24), \(\varepsilon_{ij}^{corr}\) represents the detrended residual and \(\varphi_{1}\),\(\varphi_{2}\) and \(\varphi_{3}\) are coefficients that are obtained for the residuals of each event. This method is used in the present study for detrending the residuals of PGD. Figure 8 shows the residuals of PGD (\(\varepsilon_{PGD}\)), trend line, and detrended residuals (\(\varepsilon_{PGD}^{corr}\)) versus R_{rup}. Based on the reasons mentioned above, the detrended residuals are normalized to have zero mean and unit standard deviation. The Bowman and Crujeiras (2013) test of isotropy is implemented on the normalized detrended intraevent residuals of the PGD (\(\overline{{\varepsilon_{PGD}^{corr} }}\)) that resulted in pvalues less than 0.01 for all earthquake events except Anza and Northridge earthquakes, which have the pvalue of 0.195 and 0.064, respectively. Based on these values, there is statistically significant evidence to reject the null hypothesis of isotropy of intraevent residuals of PGD for all earthquake events except Anza and Northridge earthquakes. Figure 9 shows the smoothed estimate of semivariogram in the forthroot scale for \(\overline{{\varepsilon_{PGD}^{corr} }}\) of Chuetsu, Niigata, Tottori, and Northridge earthquakes concerning distance and angle. The direction of the slower rate of increase of semivariogram (AA) can be detected by the contours near the center point of graphs.
The directional semivariogram of \(\overline{{\varepsilon_{PGD}^{corr} }}\) for selected earthquakes is obtained based on the method mentioned in the previous sections, and the rose diagrams of the normalized ranges of the directional semivariograms are presented in Fig. 10. Scattering in the experimental directional semivariograms of \(\overline{{\varepsilon_{PGD}^{corr} }}\) for Parkfield and Northridge earthquakes makes it impossible to fit a reliable semivariograms model to these data. As shown in Fig. 10 the ranges of the directional semivariograms are not equal in different directions, and consequently, the assumption of isotropy can be rejected for the earthquakes presented in this figure. It should be noted that although the pvalue of the statistical test of isotropy for the Anza earthquake is greater than 0.05, which means the hypothesis of isotropy of \(\overline{{\varepsilon_{PGD}^{corr} }}\) is not rejected in significance level of 5%, but based on the results of directional semivariograms, the \(\overline{{\varepsilon_{PGD}^{corr} }}\) of Anza earthquakes cannot be considered as a realization of an isotropic random field.
Predictive model for intraevent residual of PGD considering anisotropy
The process used to develop a predictive model of intraevent residuals of PGV considering anisotropy is used in this section to develop a predictive model of intraevent residuals of PGD. The directional semivariogram method is used, and the ranges of directional semivariograms of \(\overline{{\varepsilon_{PGD}^{corr} }}\) of selected earthquakes are obtained. Figure 10 presents a summary of characteristics of directional semivariograms of \(\overline{{\varepsilon_{PGD}^{corr} }}\) and normalized V_{s30} values, including anisotropy angle, anisotropy ratio, and maximum range of directional semivariograms. As for intraevent residuals of PGV and normalized V_{s30}, no reliable directional semivariograms can be fitted on the experimental directional semivariograms of \(\overline{{\varepsilon_{PGD}^{corr} }}\) for Northridge and Parkfield earthquakes.
As shown in Fig. 10 the AA of \(\overline{{\varepsilon_{PGD}^{corr} }}\) is in good agreement with the AA of normalized V_{s30} values. Moreover, given a low correlation range of the normalized V_{s30} values of the Northridge and Parkfield record stations, no correlation range is reported for the \(\overline{{\varepsilon_{PGD}^{corr} }}\) of these two earthquakes.
Figure 11 depicts the plot of AR of \(\overline{{\varepsilon_{PGD}^{corr} }}\) (AR_{PGD}) versus AR of the normalized V_{s30} (AR_{Vs30}) of different earthquakes. According to this figure, a positive correlation can be identified between AR_{PGD} and AR_{Vs30} with a correlation coefficient of 0.96. The linear relationship between AR_{PGV} and AR_{Vs30} obtained by linear regression analysis as:
The maximum range of \(\overline{{\varepsilon_{PGD}^{corr} }}\)(b_{maxPGD}) of different earthquake events have a mean value of 118 km and standard deviation of 23 km and are not correlated with the average sourcetosite distance of recorded points (despite what is discussed in previous sections for \(\varepsilon_{PGV}\)). More investigations about the anisotropy properties of random fields of earthquake IMs need more ground motion data recorded in sufficient and regular spatial locations. At the current level of available data, the authors propose to consider the b_{maxPGD} as a random variable with a mean value of 118 km and a standard deviation of 23 km.
The procedure for using the proposed predictive models
In this section, a stepbystep procedure for using the proposed predictive models for PGV and PGD is presented. To estimate the spatially correlated values of PGV or PGD in an earthquake event i and n locations, the following steps can be utilized:

1.
Obtain the median value of intensity measure (PGV or PGD) at n points of interest using GMPE. According to the Wang and Du (2013), “very similar results can be derived by using different GMPEs, indicating that spatial correlations are not dependent on specific GMPEs”.

2.
Obtain the spatial correlation properties of V_{s30} values in the site of interest. If it is possible to fit a reliable directional semivariogram model of V_{s30} values, go to the next step. Otherwise, suppose that the correlation range of IM values is approximately 0 km and simulate the intraevent residuals (ε_{ij}) using a univariate normal distribution model with zero mean and standard deviation of σ and then go to step 8. The value of σ is determined in respective GMPE (Campbell and Bozorgnia (2014) for PGV and Campbell and Bozorgnia (2008) for PGD).

3.
In the case of reliable directional semivariogram models for V_{s30} values, find the azimuth of direction with the maximum range of correlation (AA_{Vs30}) and calculate the anisotropy ratio of directional semivariograms (AR_{Vs30}). The anisotropy angle of IM (AA) is aligned with the AA_{Vs30}.

4.
Calculate the anisotropy ratio (AR) of IM using the AR_{Vs30} and Eqs. (18) and (25) for PGV and PGD, respectively,

5.
Calculate the harmonic mean of R_{RUP} of interested points and determine the maximum range (b_{maxPGV}) of PGV using the Eq. (23). The maximum range of PGD (b_{maxPGD}) can be considered as a random variable with a mean value of 118 km and a standard deviation of 23 km.

6.
Construct an n × n correlation matrix \({\mathbf{R}} = [\rho ({\mathbf{h}}_{ik} )]\) based on previously determined anisotropy parameters where \(\rho ({\mathbf{h}}_{ik} )\) is the correlogram between values in points i and k.

7.
Simulate the intraevent residuals (ε_{ij}) using the correlation matrix R and standard deviation value of σ, determined by respective GMPE.

8.
Simulate the interevent residuals (η_{j}) using a univariate normal distribution with zero mean and standard deviation of τ, determined in respective GMPE.

9.
Combine the median values of intensity measure (determined in step 1), intraevent residuals (determined in step 2 or 7 depending on the case), and interevent residuals (determined in step 8) using Eqs. (1) and (2) to obtain spatially correlated values of intensity measures at n points.
Conclusion
A nonparametric test of isotropy has been implemented to examine the hypothesis of isotropy of spatial correlations of PGV and PGD intraevent residuals. Using more than 2800 recorded ground motions from 9 earthquakes, it has been concluded that there is no sufficient evidence to support the assumption of isotropy. Consequently, the set of intraevent residuals of PGV and PGD should be considered as realizations of anisotropic random fields.
The Campbell and Bozorgnia (2014) and Campbell and Bozorgnia (2008) GMPEs were used for computing the intraevent residuals of PGV and PGD, respectively. The obtained residuals were normalized to have unit standard deviation and zero mean. The directional semivariogram method was employed in four azimuths of 0°, 45°, 90°, and 135° to quantify the anisotropic correlation properties of normalized intraevent residuals of PGV, PGD, and normalized V_{s30} values as an indicator of regional site condition.
It is shown that the anisotropy properties of normalized intraevent residuals of intensity measures are consistent with the anisotropy properties of regional site conditions. In the case of the heterogeneous site where there is no spatial correlation of regional site conditions, intraevent residuals of the IMs can be supposed uncorrelated. In contrast, in the case of the homogenous site, which implies more correlation range of V_{s30} values, the anisotropic angle of IMs is aligned with the anisotropic angle of normalized V_{s30} values, and the anisotropic ratio of IMs have good linear correlations with the anisotropic ratio of normalized V_{s30} values. Based on the conducted investigations, the maximum range of directional semivariograms of normalized intraevent residuals of PGV has a good correlation with the average sourcetosite distance of interested points. A linear regression model is proposed to predict the maximum range of normalized intraevent residuals of PGV as a function of the harmonic mean of the sourcetosite distance of interested locations. However, it is identified that the effect of independent factors like sourcetosite distance, earthquake magnitude, and the focal mechanism is not significant on the maximum range of directional semivariograms of normalized intraevent residuals of PGD. Consequently, it is proposed that the maximum range of intraevent residuals of PGD is considered as a random variable with a mean value of 118 km and a standard deviation of 23 km. By determining the anisotropic angle, anisotropic ratio, and maximum correlation range, predictive models are proposed to simulate the correlated values of intraevent residuals of PGV and PGD considering anisotropy.
Finally, it should be noted that other predictive parameters, such as earthquake event magnitude or fault style, may have effects on the anisotropic spatial correlation properties of PGV and PGD; however, due to a lack of available earthquake data that can be used in the spatial statistics context, we are unable to establish a reliable relationship between the mentioned parameters. Using newly recorded data and earthquake simulation methods will aid in determining the relationship between predictive parameters and anisotropic properties of earthquake IMs.
Availability of data and materials
The data underlying this article will be shared on reasonable request to the corresponding author.
Abbreviations
 AA:

Anisotropy angle
 AR:

Anisotropy ratio
 CAV:

Cumulative absolute velocity
 GMPE:

Ground motion prediction equations
 Ia:

Arias intensity
 IIEES:

International Institute of Earthquake Engineering and Seismology
 IM:

Intensity measure
 LD:

Latent dimensions
 LMC:

Linear model of coregionalization
 PCA:

Principal component analysis
 PEER:

Pacific Earthquake Engineering Research
 PGA:

Peak ground acceleration
 PGD:

Peak ground displacement
 PGV:

Peak ground velocity
 SA:

Spectral acceleration
References
Abbasnejadfard M, Bastami M, Fallah A (2020) Investigation of anisotropic spatial correlations of intraevent residuals of multiple earthquake intensity measures using latent dimensions method. Geophys J Int 222:1449–1469. https://doi.org/10.1093/gji/ggaa255
Abbasnejadfard M, Bastami M, Fallah A, Garakaninezhad A (2021) Analyzing the effect of anisotropic spatial correlations of earthquake intensity measures on the result of seismic risk and resilience assessment of the portfolio of buildings and infrastructure systems. Bull Earthq Eng. https://doi.org/10.1007/s1051802101203z
Apanasovich TV, Genton MG (2010) Crosscovariance functions for multivariate random fields based on latent dimensions. Biometrika 97:15–30. https://doi.org/10.1093/biomet/asp078
Boore DM, Gibbs JF, Joyner WB, Tinsley JC, Ponti DJ (2003) Estimated ground motion from the 1994 Northridge, California, earthquake at the site of interstate 10 and La Cienega Boulevard bridge collapse, West Los Angeles, California. Bull Seismol Soc Am 93:2737–2751. https://doi.org/10.1785/0120020197
Bowman AW, Azzalini A (1997) Applied smoothing techniques for data analysis: the kernel approach with SPlus illustrations, vol 18. OUP, Oxford
Bowman AW, Azzalini A (2014) R package ‘sm’: nonparametric smoothing methods (version 2.25.4)
Bowman AW, Crujeiras RM (2013) Inference for variograms. Comput Stat Data Anal 66:19–31
Campbell KW, Bozorgnia Y (2008) NGA ground motion model for the geometric mean horizontal component of PGA, PGV, PGD and 5% damped linear elastic response spectra for periods ranging from 0.01 to 10 s. Earthq Spectra 24:139–171. https://doi.org/10.1193/1.2857546
Campbell KW, Bozorgnia Y (2014) NGAWest2 ground motion model for the average horizontal components of PGA, PGV, and 5% damped linear acceleration response spectra. Earthq Spectra 30:1087–1114. https://doi.org/10.1193/062913EQS175M
Chiles JP, Delfiner P (2012) Modeling spatial uncertainty. https://doi.org/10.1007/s1339801401737.2
Chilès JP, Delfiner P (2012) Geostatistics: modeling spatial uncertainty, 2nd edn. Wiley, Hoboken
Cressie N (1993) Statistics for spatial data: Wiley series in probability and mathematical statistics. Wiley, Hoboken
Cressie N, Hawkins DM (1980) Robust estimation of the variogram: I. J Int Assoc Math Geol 12:115–125
Du W, Wang G (2013) Intraevent spatial correlations for cumulative absolute velocity, arias intensity, and spectral accelerations based on regional site conditions. Bull Seismol Soc Am 103:1117–1129. https://doi.org/10.1785/0120120185
Eilers PHC, Marx BD (1996) Flexible smoothing with Bsplines and penalties. Stat Sci 11:89–102
Garakaninezhad A, Bastami M (2017) A novel spatial correlation model based on anisotropy of earthquake groundmotion intensity. Bull Seismol Soc Am 107:2809–2820. https://doi.org/10.1785/0120160367
Garakaninezhad A, Bastami M (2019) Intraevent spatial correlation model for the vertical component of response spectral accelerations. J Seismol 23:853–867
Garakaninezhad A, Bastami M, Soghrat MR (2017) Spatial correlation for horizontal and vertical components of acceleration from northern Iran seismic events. J Seismol. https://doi.org/10.1007/s1095001796798
Goda K, Atkinson GM (2009) Probabilistic characterization of spatially correlated response spectra for earthquakes in Japan. Bull Seismol Soc Am 99:3003–3020
Goda K, Hong HP (2008) Spatial correlation of peak ground motions and response spectra. Bull Seismol Soc Am 98:354–365. https://doi.org/10.1785/0120070078
Goovaerts P (1997) Geostatistics for natural resources evaluation. Oxford University Press on Demand, Oxford
Jayaram N, Baker JW (2008) Statistical tests of the joint distribution of spectral acceleration values. Bull Seismol Soc Am 98:2231–2243. https://doi.org/10.1785/0120070208
Jayaram N, Baker JW (2009) Correlation model for spatially distributed groundmotion intensities. Earthq Eng Struct Dyn 38:1687–1708
Lee R, Kiremidjian AS (2007) Uncertainty and correlation for loss assessment of spatially distributed systems. Earthq Spectra 23:753–770
Markhvida M, Ceferino L, Baker JW (2018) Modeling spatially correlated spectral accelerations at multiple periods using principal component analysis and geostatistics. Earthq Eng Struct Dyn 47:1107–1123. https://doi.org/10.1002/eqe.3007
Park J, Bazzurro P, Baker JW (2007) Modeling spatial correlation of ground motion intensity measures for regional seismic hazard and portfolio loss estimation. In: Applications of statistics and probability in civil engineering. Taylor & Francis Group, London, pp 1–8
PEER NGA Database (2019) Pacific Earthquake Engineering Research Center Next Generation Attenuation Database
R Core Team (2014) R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna
Wang G, Du W (2013) Spatial crosscorrelation models for vector intensity measures (PGA, Ia, PGV, and SAs) considering regional site conditions. Bull Seismol Soc Am 103:3189–3204. https://doi.org/10.1785/0120130061
Wang M, Takada T (2005) Macrospatial correlation model of seismic ground motions. Earthq Spectra. https://doi.org/10.1193/1.2083887
Weatherill GA, Silva V, Crowley H, Bazzurro P (2015) Exploring the impact of spatial correlations and uncertainties for portfolio analysis in probabilistic seismic loss estimation. Bull Earthq Eng 13:957–981. https://doi.org/10.1007/s1051801597305
Weller ZD, Hoeting JA et al (2016) A review of nonparametric hypothesis tests of isotropy properties in spatial data. Stat Sci 31:305–324
Zafarani H, Ghafoori SMM, Adlaparvar MR (2021) Spatial correlation of peak ground motions and pseudo spectral acceleration based on the Iranian multievent datasets. J Earthq Eng. https://doi.org/10.1080/13632469.2021.1911882
Acknowledgements
The authors acknowledge supports and the funding provided by the International Institute of Earthquake Engineering and Seismology (IIEES). The authors also acknowledge the strong ground motion data provided by Pacific Earthquake Engineering Research Center (PEER) Next Generation Attenuation (NGA) database (PEER NGA Database 2019).
Funding
This research was under the support and funding provided by the International Institute of Earthquake Engineering and Seismology (IIEES) under Grant Number 7552.
Author information
Affiliations
Contributions
MA: Methodology, Software, Formal analysis, Investigation, Writing—Original Draft, Visualization. MB: Conceptualization, Methodology, Validation, Resources, Writing—Review & Editing, Supervision, Project administration, Funding acquisition. AF: Methodology, Validation, Review & Editing. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Abbasnejadfard, M., Bastami, M. & Fallah, A. Investigating the spatial correlations in univariate random fields of peak ground velocity and peak ground displacement considering anisotropy. Geoenviron Disasters 8, 24 (2021). https://doi.org/10.1186/s4067702100196w
Received:
Accepted:
Published:
Keywords
 Seismic anisotropy
 Statistical seismology
 Earthquake
 Seismic hazard assessment