Temporal and Spatial Evolution Analysis of Earthquake Events in California and Nevada Based on Spatial Statistics

: Studying the temporal and spatial evolution trends in earthquakes in an area is beneﬁcial for determining the earthquake risk of the area so that local governments can make the correct decisions for disaster prevention and reduction. In this paper, we propose a new method for analyzing the temporal and spatial evolution trends in earthquakes based on earthquakes of magnitude 3.0 or above from 1980 to 2019 in California and Nevada. The experiment’s results show that (1) the frequency of earthquake events of magnitude 4.5 or above present a relatively regular change trend of decreasing–rising in this area; (2) by using the weighted average center method to analyze the spatial concentration of earthquake events of magnitude 3.0 or above in this region, we ﬁnd that the weighted average center of the earthquake events in this area shows a conch-type movement law, where it moves closer to the center from all sides; (3) the direction of the spatial distribution of earthquake events in this area shows a NW–SE pattern when the standard deviational ellipse (SDE) method is used, which is basically consistent with the direction of the San Andreas Fault Zone across the north and south of California; and (4) the spatial distribution pattern of the earthquake events in this region is found to be clustered using the global spatial autocorrelation analysis method. This study provides a new perspective for the exploration of the temporal and spatial evolution trends in earthquakes and understanding the earthquake risk in an area.


Introduction
Earthquakes are one of the main natural disasters on Earth [1,2]. An earthquake of magnitude 7.8 occurred in San Francisco, California, USA, on 18 April 1906, affecting most areas of California, from Oregon City to Los Angeles, causing more than USD 4 million in economic damage [3,4]. The prediction of earthquakes is an important task worldwide [5]. If the location, time, and magnitude of future large earthquakes can be accurately predicted, millions of lives may be saved. Although it is difficult to predict the occurrence of earthquakes, recent studies have shown that we can analyze the trends in the change in earthquake risk in a certain region by studying the temporal and spatial evolution processes of earthquake events [6,7], which is beneficial to disaster prevention and reduction.
Earthquake events can be regarded as data with geospatial attributes. Due to the spatial constraints of different directions and distances between the spatial data, the traditional statistical analysis methods cannot be used to describe constraint relationships between these spatial data. Spatial statistics are a more recent application of spatial data [8,9]; the discipline combines statistics with modern graphical computing technology and shows the implied spatial distribution, spatial pattern and spatial interactions in spatial data using intuitive methods [10]. After continuous development, the spatial statistics method has been widely applied in fields such as sociology, biology, demography, criminology and geology [11].
From the perspective of spatial statistics, the temporal and spatial evolution processes of things in geospatial data have certain rules. It is helpful for us to study and master these rules in all aspects of our lives. In the response to disasters, for example, Gan et al. [12] studied the temporal and spatial evolution of the vegetation in the Mianyuan River Basin using remote sensing images of vegetation from 1994 to 2017. The results showed that slopes of 30 to 40 degrees in this area are more susceptible to earthquake disturbances but are easier to repair, and the recovery period of vegetation in this area is eight to nine years after being damaged. Wang et al. [13] studied the temporal and spatial evolution trends in droughts in Northern Shanxi using the multi-model ensemble mean value; the results showed that the precipitation and temperature in this region will increase from south to north in the future, with the largest increase in the west. This study provided information that can be used to alleviate the threat of drought in this area. Based on the measured data of 56 rainfall stations and 11 hydrological stations in Central Guizhou, China, He et al. [14] developed the lagging index (LI) to analyze the temporal and spatial evolution characteristics of rainfall runoff in this area. The results showed that the lagged effect of runoff to rainfall in the area mainly demonstrated three lagged periods, and the order of temporal and spatial distribution differences of the lagged intensities was two lagged periods (Cv = 0.62). Zhang et al. [15] used the method of spatial statistics to conduct a systematic statistical analysis of the Iburi landslides, and determined several of its specific characteristics, such as high concentration and high mobility in the severedisaster area. Qi et al. [16] studied the remote sensing data of landslides caused by the Wenchuan earthquake, and found that landslides' incidence varies in different slopes with different structures.
Many researchers have explored the temporal and spatial evolution characteristics of earthquake events. For example, Dias et al. [17] proposed a method for studying the temporal and spatial probability distributions of earthquakes, where the influence of the earthquakes' data set on the temporal and spatial probability distribution of earthquakes was discussed based on different depths of hypocenters and thresholds. The experimental results showed that the increase in the jump forming the sequence between earthquakes affects the non-extensive characteristics in the temporary probability distribution. This analysis demonstrated the role of the robustness of the non-extensive statistical mechanics treatments in earthquake research, showing that it is an effective method with which to test memory and interactions. Yang et al. [18] used the temporal and spatial scan method to identify earthquake clusters based on the data of earthquakes with a moment magnitude (Mw) of not less than 5.6 from 1960 to 2014. The results indicated that earthquake clusters can be classified into two types based on duration: persistent clusters and burst clusters. Huang et al. [19] used the matched-filter detection method to obtain a relatively complete (magnitude of completeness ≈ 0.9) and precisely relocated earthquake catalog during the 2019 M 6.4 and M 7.1 Ridgecrest, California earthquake sequence, and discussed the influence of the spatial-temporal evolution process of the foreshock sequence on the preparation and nucleation process of the mainshock, which is important to improve our understanding of fault interactions, earthquake triggering, and evolving earthquake hazards during an ongoing earthquake sequence.
In recent years, with the rapid development of artificial intelligence, researchers have begun to use artificial intelligence to study earthquake risk. Jena et al. [20] proposed a new earthquake risk assessment method that combines artificial neural network cross-validation (four-fold ANN-CV) and the hybrid Analytic Hierarchy Process-Technique for Order of Preference by Similarity to Ideal Solution (AHP-TOPSIS), and applied it to Aceh, Indonesia. The experimental results showed that the accuracy of the model was 85.4% and the consistency ratio was 0.06. The model provides the advantages of a simple structure, convenient parameter adjustment, easy extension to other areas, and conduciveness to earthquake risk prevention and mitigation. Jia et al. [21] propose a machine learning method for seismic data interpolation. DeVries et al. [22] proposed a method of deep learning to predict the spatial distribution of aftershocks. Karimzadeh et al. [23] used different machine learning algorithms to predict the aftershock patterns. Mignan et al. [24] compared the performance of DNN and logistic regression in forecasting the spatial distribution of aftershocks in the aftermath of large seismic events.
If the temporal and spatial evolution rules of earthquake events in a certain region with frequent earthquakes are mastered, local governments can make correct decisions regarding disaster prevention and reduction. There are many fault belts in California and Nevada, which typically experience frequent earthquake events [25]. Unlike the above research, we used the weighted average center, standard deviational ellipse (SDE), and global spatial autocorrelation in spatial statistics to study the temporal and spatial evolution characteristics in the earthquake events in California and Nevada in the past 40 years, and analyzed the change in the trend of earthquake risk in this area based on experimental results from the above methods.

Study Area
California and Nevada are in the West United States and have many fault belts that are amongst the most active in the world. We attempted to analyze the trend in the change in the earthquake risk in this region by studying the temporal and spatial evolution characteristics of earthquake events in California and Nevada. Limited by the accuracy of seismic instruments and the technology of seismic location, earthquakes below a magnitude of 3.0 in different earthquake catalogs may be different, and there are few seismic data-from before 1980, so the earthquake events of magnitude 3.0 or above in California and Nevada from 1980 to 2019 were selected as the study objects. The research area (115.24 •~1 22.54 • W, 32.97 •~3 8.37 • N) was selected, as shown in Figure 1. A total of 15,762 earthquake events of magnitude 3.0 or above occurred in this region from 1980 to 2019, including 22 earthquake events of magnitude 6.0 or above, and they were obtained from the USGS earthquake catalog (https://earthquake.usgs.gov/). As the purpose of this study was to analyze the change in the trend of earthquake risk in a certain region by studying the temporal and spatial evolution characteristics of earthquake events in this region, the earthquake events were divided into multiple segments in the time dimension.

Method
The epicenter of an earthquake can be regarded as a point on the map, and its latitude and longitude as its coordinates of point (x, y). In this study, we used four main methods to examine the temporal and spatial evolution characteristics of earthquake events in different time slices and to analyze the risk of earthquake events in California and Nevada. First, the frequency rule of seismicity in the research region was analyzed using the number of earthquakes in different time slices. Secondly, we used the weighted average center method to calculate the center of the earthquake events in the different time slices and analyzed the temporal and spatial evolution of the spatial concentration trends in earthquake events in the research area. Then, the SDE model of spatial statistics was used to study the temporal and spatial evolution of the directional distribution of earthquake events in different time slices. Finally, global spatial autocorrelation was used to study the temporal and spatial evolution process of the spatial distribution patterns of earthquake events in different time slices.
The standard deviational ellipse and weighted average center methods used in this study need to be supported by enough earthquake events; otherwise, large calculation errors occur. The longer the single time segment, the more reliable the result. To study whether the earthquake events in this area have regularity in the time dimension, we needed more time slices. The larger the number of time segments, the more obvious the result. This created a contradiction. To meet the above research needs, after comprehensive consideration, we divided the 40 years (1980-2019) of earthquake events in the California-Nevada region into 10 time segments, each of which had 4 years' of earthquake events. Earthquake events were classified according to the international standard of Richter magnitude. The statistics of the earthquakes of each time slice are shown in Table 1.

Weighted Average Center
The weighted average center algorithm reflects the influence of the importance of each datum in a group of data on the overall concentration trend. The algorithm assumes that we have a finite data set f i ; different data in the data set are assigned different weights, and the weight represents the importance of the datum in the data set. Then, the weighted average value of the data set is the sum of the product of each number and its weight divided by the sum of all weights [26]. The size of the weighted average depends not only on the size of each value, but also on the number of times each value appears. The epicenter coordinates of earthquake events in different time slices were regarded as a data set; considering that many aftershocks often occurred near the major earthquake, we used earthquake magnitude as the weight of the earthquake events in the research region to calculate the weighted average center coordinate of earthquake events in the different time slices, which was useful for analyzing the temporal and spatial evolution process of spatial concentration trends of earthquake distribution. The weighted average center was calculated as follows: where x i and y i are the latitude and longitude coordinates of the earthquake event i, respectively; n is the total number of earthquake events in different time slices; and w i is the magnitude of earthquake event i. The weighted average center analysis tool in the ArcGIS 10.5 toolbox was used to analyze the earthquake events of the different time slices. The earthquake event subsets of different time slices were selected as the input elements, and the earthquake magnitude was used as the weight field to obtain the weighted average centers of the earthquake events of magnitude 3.0 or above in the different time slices.

Standard Deviational Ellipse
The SDE in spatial statistics was used to describe the spatial distribution characteristics of geographical elements in a certain region. It was first proposed by Lefever [27][28][29] in 1926 and has been widely applied in fields such as sociology, demography, criminology, geology and ecology. The spatial distribution characteristics of the studied object can be described quantitatively using the spatial distribution ellipse of the studied object with the major axis, minor axis and azimuth as basic parameters. The spatial distribution ellipse takes the average center of the spatial distribution of geographical elements as the center, the azimuth reflects the main trend direction of the distribution of the study object, the major axis direction of the ellipse represents the direction with more spatial distribution of geographical elements, and the minor axis direction represents the direction with less spatial distribution of the geographical elements [30,31]. The larger the difference between the major axis and the minor axis, the more significant the directionality of the geographical elements. If the length of the major axis is equal to the length of the minor axis, the distribution of geographical elements has no directional feature. These main parameters are calculated as follows: where x i and y i are the latitude and longitude coordinates of the earthquake event i, respectively; n represents total number of earthquake events in different time slices; and X and Y are the SDE centers in different time slices.
where θ is the azimuth of the ellipse, indicating the angle formed by the clockwise rotation of the north direction to the major axis of the ellipse; x i and y i denote the coordinate deviation from the coordinate of earthquake event i and j to the average center; σ x is the length of the ellipse major axis; and σ y is the length of the ellipse minor axis.
To create the SDE model for earthquake events of different time slices in the research region, the directional distribution tool in the ArcGIS 10.5 toolbox was used first. Then, the earthquake events of different time slices were selected as the input elements, the Ellipse Size was set to 1_STANDARD_DEVIATION, and the earthquake magnitude was used as the weight field to represent the importance of the earthquake in the event set. Finally, the standard deviation ellipses of different time slices were generated.

Global Spatial Autocorrelation Analysis
Global spatial autocorrelation is an evaluation criterion of the correlation between variables and their spatial positions, which measures the relationship between the attributes and positions of spatial features [32][33][34]. According to Tobler's first law of geography, everything is related to everything else, and near things are more related than distant things [35][36][37], which coincides with definition of global spatial autocorrelation of a point element. The global Moran's index is a popular statistical index of global spatial autocorrelation, which is used to analyze the spatial distribution pattern of research objects in the whole research area, whether it shows a cluster pattern, and the overall similarity.
We used the spatial analysis tools in ArcGIS 10.5 to evaluate the global spatial autocorrelation of earthquake events of the different time slices in the research area.
Generally speaking, the spatial distribution of point patterns is divided into three basic types: dispersed, random and clustered. The quadrat analysis method is the most commonly used intuitive method for studying the distribution types of spatial point patterns. The spatial distribution types of earthquake events in California and Nevada were analyzed by using the quadrat analysis method.
Firstly, we divided the research area into small grids by creating fishing nets, and then counted the number of earthquake events in the different grids. According to the Smith's experiment [38][39][40], the optimal calculation formula for creating fishing nets is: where S is the area of the grid, A represents the area of the research area, and r represents the total number of earthquake events in the different time slices in the study area. Although any closed figure can theoretically be used as the shape of the grid, the square is the best choice as it is easy to construct and merge into a large figure. Thus, we chose squares to construct the grid. We used Equation (9) to deduce a, which is the side length of the grid.  Secondly, we executed the Spatial Join Processing command, which is used to create statistics of the number of earthquakes of each grid in the fishing net, and the numbers of earthquakes of each grid were selected as the analysis field in the global spatial autocorrelation analysis method.
Finally, ArcGIS produced the analysis report, including the global Moran's index, p-value and z-score. The global Moran's index is expressed as follows [41]: where x i and x j represent the number of earthquake events in grids i and j, respectively; x represents the average number of earthquake events in each grid Z i and Z j represent the number of deviations between the number of earthquake events in grids i and j and the average number of earthquake points in each grid, respectively; w i,j represents the spatial weight between grid i and grid j; n represents the total number of earthquake points in the different time slices; and s 0 is the cluster of spatial weights for all grids. The value of the global Moran's index varies between 1 and −1. When the global Moran's index is greater than 0, the spatial distribution of earthquake events tends toward clustering; the larger the global Moran's index, the more obvious the clustering. When Moran's index equals 0, there is no obvious spatial relationship, which means that the spatial distribution of the earthquake events is random; otherwise, it is dispersed. [42].
In Equation (10), it is necessary to know the value of the spatial weight w when calculating the global Moran's index. ArcGIS 10.5 provides six methods to generate spatial weight. Among the six methods, the INVERSE_DISTANCE method supposes that the spatial relationship between the element and other elements is a relationship that represents attenuation with an increase in distance, which is coincident with the first law of geography and so is also applicable to geographical phenomena with obvious spatial relationships such as earthquakes. Thus, we chose the INVERSE_DISTANCE method to generate the spatial weight. In the CONTIGUITY_EDGES_ONLY method, the relationship between elements is only affected by the elements in a common boundary or overlapping, but the research goal of this paper was to analyze the spatial autocorrelation of earthquake events in the whole study area. Thus, this method was not suitable for our study, and none of the other methods were more suitable for the study than the INVERSE_DISTANCE method. In addition to the spatial weight, the Distance Method and Standardization should be set when the global spatial autocorrelation analysis is conducted. In this study, the Distance Method parameter was set to EUCLIDEAN_DISTANCE, and the Standardization parameter was set to ROW. After setting the above parameters, the global spatial autocorrelation of earthquake events in this area was analyzed, and the analysis report was generated.

Frequency Statistical Results
An earthquake is a phenomenon of the slow accumulation and then quick release of energy in the crust. When the energy in the Earth's crust in a certain area is in the accumulation period, the seismicity in the region is weak, the frequency of earthquake events is low, and the earthquake risk is low; by contrast, when the energy is in the release period, the seismicity in this area is strong, the frequency of earthquake events is relatively high, and the earthquake risk is high. According to the international standard of Richter magnitude, earthquake events of magnitude 4.5 or above are middle-strong earthquakes, which may cause harm to buildings and people's lives in the surrounding areas. Therefore, it is meaningful to study the frequency change of earthquake events with magnitude 4.5 or above. According to the statistical results (Table 1), we find that the frequency of earthquake events of magnitude 4.5 or above in California and Nevada presents a regular decreasing-rising trend as seen in Figure 3. Using four time slices as a group, the frequency of earthquake events of the first three time slices decreases gradually and increases in the last time slice; the next group presents the same change rule. As shown in Figure 3

Temporal and Spatial Evolution Process of the Weighted Average Center of Earthquake Events
We used the weighted average center analysis tool of ArcGIS 10.5 to obtain the weighted average centers of the earthquake events of magnitude 3.0 or above in different time slices in the research region and then connect them in chronological order. Figure 4 shows the temporal and spatial evolution trace of the weighted average center in earthquake events in different time slices and the epicenters of the earthquakes of magnitude 6.0 or above in the research region As shown in Figure 4, the weighted average center of the earthquake events in this area shows a conch-type movement law, where it moves close to the center from all sides. The initial weighted average center of earthquake events in Fresno is close to the Owens Fault Zone (time slice : 1980-1983). Then, the moving path of the weighted average center of the earthquake is like a conch, moving closer to the center from all sides. Figure 4 also shows that the epicenters of most historical great earthquakes move with the direction of the weighted average center of earthquake events in that time slice. For example, two earthquakes of magnitude 6.0 or above occurred near the Owens Fault Zone and one earthquake of magnitude 6.0 or above near the San Andreas fault zone in the 1980-1983 time slice, and two earthquakes with magnitudes of 7.1 and 6.4 occurred near the Garlock Fault Zone in 2016-2019 time slice.

Temporal and Spatial Evolution Process of the SDE
The SDEs of the earthquakes of magnitude 3.0 or above in the different time slices in the research area were generated using the directional distribution tool in the ArcGIS 10.5 toolbox. Figure 5 illustrates the temporal and spatial evolution processes of the SDE of the earthquake events in different time slices. The statistical parameter analysis results of the SDE in different time slices are listed Table 2. In the active periods of frequent earthquakes (e.g., [1980][1981][1982][1983][1992][1993][1994][1995], and 2016-2019 time slices), the SDE areas were relatively small because of the concentrated distribution of earthquakes, and the directions of earthquake distribution were obvious because of the large differences between the long axis and the short axis. The azimuth of the SDE reflects the main trend direction in earthquake distribution; from Table 2, the standard deviation ellipse rotation angles of the different time slices were all between 115 • and 137 • , indicating that the NW-SE spatial distribution of earthquakes was dominant in the research area, which is basically consistent with the direction of the San Andreas Fault Zone across the north and south of California.

Temporal and Spatial Evolution Process of Global Spatial Autocorrelation Analysis
The global spatial autocorrelation analysis report of the earthquake events of magnitude 3.0 or above in the different time slices in the research area were obtained using the spatial autocorrelation analysis tool in the ArcGIS 10.5 toolbox. In addition to the global Moran's index, the analysis results in the report also include z-scores and p-values. In general, the z-score and p-values in the report were analyzed first. The z-score being greater than 1.96 and the p-value being less than 0.05 indicate that the analysis result of the data set is statistically reliable [43]. Secondly, we analyzed the global Moran's index in the report, which represented the spatial distribution pattern of earthquake events in this area.
The analysis results of the time slice from 1980 to 1983, as an example (Figure 6), show that the p-value is less than 0.05, the z-score is greater than 1.96, and the global Moran's index is greater than 0, which shows that the spatial distribution pattern of earthquake events in this period of the research region was clustered. All the spatial autocorrelation analysis reports of the earthquake events in the different time slices are shown in Table 3. The p-values in the global spatial autocorrelation analysis results of the earthquake events in different time slices are less than 0.05, and z-scores are all greater than 1.96, and the global Moran's Indexes are greater than 0, which indicate that the spatial distribution pattern of the earthquake events of magnitude 3.0 or above in California and Nevada is clustered and the results are statistically reliable.

Discussion and Conclusions
In this study, we used statistical spatial analysis methods to study the temporal and spatial evolution trends in earthquake events of magnitude 3.0 or above in California and Nevada, and to explore the earthquake risk in this region based on the experimental results.
According to the statistical results of earthquake events, we found that the frequency of earthquakes of magnitude 4.5 or above in California and Nevada presents a regular change trend of decreasing-rising every 16 years, where the number of earthquakes gradually increases in the first 12 years, and then gradually decreases in the next four years. If this rule is true, the frequency of earthquakes of magnitude 4.5 or above in the next time slice (2020-2023) will be less than that of 2016-2019 as this time slice is expected to be in the decline phase of this cycle. However, it does not mean that no large earthquake will occur in the next four years in the research area.
The spatial and temporal evolution results of the weighted average centers method showed that the spatial concentration trend in earthquake events in California and Nevada shows a conch movement pattern, which means the epicenter of a large earthquake is approaching the center of the study area. If the law is true, the weighted average center of earthquake events in the next time slice is likely to move toward the south-west.
The result of the temporal and spatial evolution of SDE showed that the dominant direction of spatial distribution of earthquakes in California-Nevada is NW-SE, which is consistent with the direction of the San Andreas Fault Zone. This result agrees with the trend in the foreshock sequence during the 2019 M 6.4 and M 7.1 Ridgecrest earthquakes [19].
In summary, this study provides a new perspective for the exploration of the temporal and spatial evolution trends of earthquakes and for understanding the earthquake risk in an area. The experimental results of the temporal and spatial evolution trends in this paper are consistent with other research, which indicates the validity of using the spatial statistical method to study the spatial and temporal evolution characteristics of earthquake events. The approach could be extended to other regions in the future. The earthquake catalog is important for identifying the temporal and spatial evolution characteristics of earthquake events; the more complete the earthquake catalog, the more accurate the research results. However, considering the error of earthquake data recorded in the USGS catalog and the uneven distribution of stations, the experimental results obtained in this paper may show deviation.