Rapid mapping of seismic intensity assessment using ground motion data calculated from early aftershocks selected by GIS spatial analysis

Abstract

Following a major earthquake, disaster information services must deliver accurate damage assessment results during the emergency ‘black box’ phase when data is scarce. Seismic intensity maps contain crucial information for determining the damage in the affected area. For earthquakes with Mw between 5.5 and 7, this study proposes using GIS analysis to mine aftershock events in early aftershock sequences that are closely related to the mainshock fault, and then using these events to generate seismic intensity assessment maps. Regression curves were first obtained using a nonparametric method (rLowess) to analyse the geographical coordinates of early aftershocks. Then, a buffer of 1 or 1.5 km radius was made for the curve, and the aftershocks in the buffer were used to calculate the predicted peak ground velocity (PGV) values over a specific km-grid range. Finally, rapid mapping of seismic intensity was assessed based on the intensity scale. This straightforward and repeatable method employs seismic station data obtained shortly after the mainshock. The assessed seismic intensity accurately reflects the location and extent of the hardest hit areas and can be cross-referenced with geophysical results to accurately assess the damage in the affected areas.

Keywords:

• locally weighted regression
• GIS spatial analysis
• early aftershock sequence
• seismic intensity assessment map
• disaster information service

1. Introduction

Destructive earthquakes frequently cause serious human casualties and property damage. Earthquake emergency response is a paramilitary operation to mitigate earthquake damage that aims to make scientific decisions and implement them in the shortest amount of time, which is the most visible manifestation of the urgency and necessity of disaster information services (Nie et al. 2012). Rapid access to trustworthy damage assessment results is critical for driving and supporting earthquake emergency response operations (Poggi et al. 2021). Map products containing earthquake data provide information about the seismogenic region to a variety of audiences. Seismic intensity, as one of the important parameters for a comprehensive measure of earthquake damage, provides a quick description of the extent of potentially destructive shaking after an earthquake, which can assist decision-makers and rescuers in clarifying the extent of damage and conducting on-site work (Wald et al. 2006; Xu et al. 2020).

With the expansion of global air and land monitoring data and the advancement of computer hardware and software, earthquake disaster damage assessment tools are showing a more diverse development trend. Data from satellites and aircraft carrying optical sensors in affected areas are primarily used to assess building damage and secondary earthquake hazards like landslides, weirs, and mudslides (Dell’Acqua and Gamba 2012). Additionally, more radar data applications and studies on the rapid assessment of building damage have been conducted (Plank 2014; Polcari et al. 2017). Ground-based monitoring arrays and intensity meters provide valuable observational data for studying surface rupture processes and assessing disaster damage quickly, and instrumented intensity maps can be drawn quickly in areas with dense station populations (Okada et al. 2004; Wang and Mori 2011). Internet-based research for earthquake disaster damage assessment has grown in popularity in recent years. ‘Did you feel it?’ has collected ample new data on earthquake ground motion and impacts from online citizen responses as a relatively mature international Internet research platform for earthquake hazards, providing a new data resource for qualitative and quantitative earthquake research (Atkinson and Wald 2007). However, the representativity of the feedback collected by the system is still debated (Hough and Martin 2021; Hough and Martin 2022; Wald 2022). Some authoritative international agencies and developed countries with a high frequency of earthquakes have established earthquake damage assessment systems, which have assumed important responsibilities in regional and even global earthquake disaster information services. The results of ground shaking predictions are important outputs of these systems and are required input data in other damage assessments (Hosseinpour et al. 2021).

According to the amount of data available on the damage situation obtained from the seismic region after an earthquake, the earthquake response can be divided into three periods: the ‘black box’, ‘gray box’, and ‘white box’ periods. The ‘black box’ period refers to the time between the main earthquake and the first damage information being sent out from the disaster area, usually 2–3 h after the earthquake (Nie and An 2013; Xia et al. 2019). Decision-makers require reliable information about the disaster area during this time, but the data available for providing disaster information services are limited primarily to seismic stations, the Internet, and empirical model assessment results (Nie et al. 2012). Maps based on empirical model results are important during the ‘black box’ period (Xu et al. 2016), and seismic intensity assessment maps are necessary for reference. Therefore, various ground motion prediction equations and seismic intensity attenuation relationships have been proposed by many countries and regions. For example, the intensity attenuation relations in the earthquake damage assessment system in China frequently uses an elliptical model (Wang et al. 2013); the USGS ShakeMap system sets up appropriate attenuation relations for different regions and magnitudes (Wald et al. 2006; Worden et al. 2020). The advancement of artificial intelligence repurposes the massive amount of data and historical earthquake information accumulated by seismic monitoring stations. Neural networks trained by a large amount of data can capture patterns in complex geophysical relationships and achieve real-time or near-real-time assessment and prediction of seismic intensity, which has recently become a research hotspot (Mousavi et al. 2020; Otake et al. 2020; Fornasari et al. 2022). In particular, it has gained wide attention in the research of earthquake catalogue pickup and aftershock pinpointing, as they can quickly obtain a large number of aftershock sequences after an earthquake (Fang et al. 2017; Jiao and Alavi 2020).

The spatial distribution of aftershocks is intimately tied to the position, nature, and structure of causative faults, and it is widely accepted that stress changes across faults cause aftershocks (Kisslinger 1996; Kilb et al. 2000; Shebalin and Narteau 2017). The area of the aftershock zone might approximate the area of the surface rupture zone, and the very early aftershocks are a good indicator of the degree of rupture of the mainshock (Neo et al. 2021; Ozawa and Ando 2021). Chen et al. (2022a, 2022b) presented and verified a method for estimating fast seismic intensity based on surface rupture data from far-field array inversion combined with the ground-motion prediction equation (GMPE) based on fault distance. They provided reliable data support in the disaster information service during the emergency response phase of the Mw7.3 magnitude earthquake in Maduo, Qinghai Province, China, in 2021. Limited by the inversion accuracy of the surface rupture process, the location and extent of the hardest hit areas assessed by this method may have large deviations compared with the actual situation when the moment magnitude is relatively small (generally below Mw7.0) (Yang et al. 2022). Therefore, we propose replacing the GMPE-required fault data with aftershock sequences collected within 2 h of the earthquake to calculate peak ground acceleration (PGA) and PGV data. This method was tested and yielded reasonable intensity rapid assessment results in the 2022 Menyuan Mw6.6 earthquake in Qinghai Province, China (Zhao et al. 2022a). Subsequently, we introduced a locally weighted regression method for analysing the spatial distribution trend of aftershocks, used the regression results in the calculation of GMPE, and then assessed the seismic intensity, which was found to be more applicable to earthquakes with Mw ≥7.0 by testing over 50 earthquake cases (Zhao et al. 2022b).

Under certain conditions, earthquakes with Mw around 6 can cause severe casualties and property damage, such as the Mw 6.6 earthquake that occurred on 5 September 2022, in Luding County, Sichuan Province, China, killing over 90 people and causing severe damage to multiple infrastructures (Han et al. 2022). Therefore, timely transmission of high-quality intensity assessment results is required in emergency response activities of this magnitude. Here, we address the issue of using the local weighted regression results of early aftershock data to assess the shortcomings in the application of seismic intensity in this magnitude range. We propose an improved scheme in which a buffer of a certain distance is created for the regression curve after performing locally weighted regression on the geographic coordinates of early aftershocks to obtain a curve that can depict the spatial distribution information of the aftershock sequence. The information from the mainshock rupture is represented by aftershocks in the buffer. This portion of the aftershocks is then used to compute ground motion data and swiftly map the seismic intensity evaluation based on the ground motion prediction results. The main contribution of this article is the discovery of valuable aftershock events using simple spatial analysis methods, as well as the rapid mapping of seismic intensity assessment shortly after an earthquake, which enriches data sources for disaster information services during the emergency response ‘black box’ period.

2. Materials and methods

2.1. Materials

Because the GMPE used in this study was created by fitting seismic data from Japan and proved to be applicable in western China (Si and Midorikawa 1999; Chen et al. 2022a), all earthquake cases needing GMPE computations originated in China and Japan. We chose four earthquakes with substantial impact in China since 2021 because the country’s seismic network has steadily improved in recent years, and seismic data obtained from monitoring is more abundant and of higher quality than before (Peng et al. 2021). The macroscopic seismic intensity maps released by the China Earthquake Administration (CEA) following on-site investigation was digitized using ArcGIS software and used as the actual damage in the disaster area to test the intensity assessment results. Aftershock sequences within 2 h of the mainshock were obtained from the China Earthquake Network Center (https://news.ceic.ac.cn/). According to the seismic intensity scale used in the ShakeMap system (MMI), VIII degree intensity causes intense vibration and moderate to severe damage (García et al. 2012; Worden et al. 2020). We chose three earthquakes based on this. Earthquake magnitude, epicentre position, focal depth, and isoseismal line were downloaded from ShakeMap, and the International Earthquake Catalogue Center (http://www.isc.ac.uk/) provided the aftershock data. In addition, the 2016 Mw 7.8 earthquake in New Zealand was chosen for discussion. The VS30 data used for site correction were obtained from the USGS (https://earthquake.usgs.gov/data/vs30/), and all earthquake onset times were recorded in universal time coordinated. Both the Chinese Seismic Intensity Scale (CSI) and the MMI describe regions with intensity VIII as moderate or severe damage. The hardest-hit areas in this study are defined as having an intensity of VIII or higher. The information of all earthquakes is shown in Table 1.

Table 1. All earthquake information used in this study. The values in maximum intensity brackets are the field survey results released by the CEA.

Date	Country	Latitude	Longitude	Magnitude (Mw)	Aftershock number	Maximum Intensity
20001006	Japan	35.456°N	133.134°E	6.7	152	VIII
20080613	Japan	39.030°N	140.881°E	6.9	227	VIII
20161113	New Zealand	42.737°S	173.054°E	7.8	106	IX
20180905	Japan	42.686°N	141.929°E	6.6	162	IX
20210521	China	34.598°N	98.251°E	7.3	70	IX (X)
20220107	China	37.828°N	101.290°E	6.6	43	VIII (IX)
20220609	China	32.25°N	101.82°E	5.9	96	VII (VIII)
20220905	China	29.59°N	102.08°E	6.6	78	VIII (IX)

2.2. Methods

2.2.1. Aftershock data pre-processing

Secondary fault activation and remote triggering of stress result in a discrete distribution of some aftershocks at locations far from the epicentre (Felzer and Brodsky 2006). These aftershocks influence the regression results, moving the regression curves toward the sites of the aftershocks and interfering with the accuracy of determining aftershock spatial distribution trends. Even the robust regression program utilized in this study had difficulties identifying this fraction of aftershocks as outliers; thus, raw aftershock data must be pre-processed to exclude aftershocks with anomalous spatial distribution. We used 1.5 times the interquartile range (IQR) to look for outliers in the geographic coordinates of aftershocks. If the longitude or latitude of an aftershock event was outside that range, the aftershock was judged to have an anomalous spatial distribution and was eliminated. This approach detects the large majority of aftershocks with aberrant distribution and, if necessary, can be complemented with expert experience for manual exclusion (Zhao et al. 2022b). The range used to detect outliers is expressed as follows: (1) $$\left[\mathrm{Q}1-1.5\left(\mathrm{Q}3-\mathrm{Q}1\right),\mathrm{ }\mathrm{Q}3+1.5\left(\mathrm{Q}3-\mathrm{Q}1\right)\right]$$(1) where Q1 and Q3 are the first and third quartiles of aftershock longitude or latitude values, respectively.

2.2.2. Locally weighted regression

We anticipated obtaining a curve from the early aftershock sequence that reflected the mainshock rupture information. Cleveland (1979, 1981) proposed a locally weighted regression method for smoothing discrete point data, which is a nonparametric regression method whose application is flexible, does not require a specified function to fit a model to obtain global sample data, and is well suited for complex processes without theoretical modeling, and can also be used in the field of geophysics for analysing seismic data (Mariani and Basu 2014). In a previous study, we found that the curves fitted by this method showed some similarities to surface rupture in length and linear mean direction, and that for earthquakes with Mw ≥7.0, the fitted results could be used directly in the GMPE calculation (Zhao et al. 2022b). Lowess, rLowess, Loess, and rLoess are the four forms of the method, where Lowess and rLowess are locally weighted regression methods for univariate regressions, Loess is a development of Lowess for multivariate regressions (Cleveland and Devlin 1988), and rLowess and rLoess denote robust versions of the original model that can mitigate the interference. Figure 1 shows the regression effects of the four methods. Because using rLowess to regress the geographic coordinates of aftershocks is the most effective and simple method, we employed the rLowess method to analyse the spatial distribution information of early aftershocks in this study, and implemented this process using the Lowess function in R, which employs a robust version of Lowess by default. A detailed explanation of the difference between the Lowess and Loess functions in the R language can be found at (https://support.bioconductor.org/p/2323/).

Figure 1. On the same dispersed data, the regression outcomes of four locally weighted regression algorithms were compared. The blue scatter dots are aftershock events that occurred within 2 h of the 2022 Mw 6.6 earthquake in Menyuan County, Qinghai Province, China, while the red solid line represents the regression results.

The extent of surface rupture generated by earthquakes of low magnitude may be minimal (Wells and Coppersmith 1994). Most of the aftershocks distributed along both sides of the fault are con-fined within 1–1.5 km from the mainshock fault (Yukutake and Iio 2017), which is supported by Ozawa’s simulation study of mainshock and aftershock sequences in geometrically complex rupture zones (Ozawa and Ando 2021). Figure 2 depicts the spatial distribution of aftershocks from the 2022 Luding Mw6.6 earthquake, as well as the distribution of aftershock occurrence times along the fault strike. The aftershocks within a very early period after the main earthquake occurred mainly along the fault strike. This phenomenon is common, as evidenced by the Wenchuan Mw7.9 earthquake in 2008 (Yin et al. 2018). In this study, we created a buffer zone of 1 km or 1.5 km radius for the regression curves and considered that the spatial distribution of aftershocks within the buffer zone can accurately reflect the information of mainshock rupture. Their geographic coordinates were used in the GMPE calculation, and the seismic intensity assessment map was quickly drawn.

Figure 2. Aftershocks from the 2022 Mw6.6 earthquake in Luding, Sichuan Province, China: (a) spatial distribution of the aftershocks, with the red solid lines in the seismogenic region primarily being the Xianshui River and Anning River faults. (b) Temporal distribution of aftershocks, taking the logarithm of the time in seconds between the mainshock and the occurrence of the aftershock.

2.2.3. GMPE

We called the GMPE based on the shortest fault distance used in this study SM99, which is an empirical relationship obtained by fitting the data from 29 earthquakes in Japan. This equation has also been gradually applied in the seismic intensity assessment work in western China in recent years, and a preliminary methodological system for rapid seismic intensity assessment based on multiple data sources has been formed (Chen et al. 2022a; Chen et al. 2022b). In this study, we chose the regression equation with constraints, and the calculation process must generate a specific range of grid and grid centre points centred on the epicentre, with grid cell sizes of 1 km × 1 km, using the information needed by the grid centre points and the aftershock calculation equations. The GMPE is as follows: (2) $$\log \hspace{0.17em}\mathit{PGA}=0.50M_w+0.0043D+\sum d_i{S}_i-\log \hspace{0.17em}\left(X+0.0055\times {10}^{0.5M_w}\right)-0.003X+0.61+\epsilon$$(2) (3) $$\log \hspace{0.17em}\mathit{PGV}=0.58M_w+0.0038D+\sum d_i{S}_i-\log \hspace{0.17em}\left(X+0.0028\times {10}^{0.5M_w}\right)-0.002X-1.29+\epsilon$$(3) (4) $$X=\sqrt{{Rh}^2+{Zh}^2}$$(4) (5) $$lg\left(\mathit{Amp}\right)=1.83-0.66lg{V}_{V_{S30}}$$(5) (6) $${\mathit{PGV}}_{VS30}=\mathit{PGV}\times \mathit{Amp}$$(6) where D denotes the hypocentre depth and S denotes the type of faults, divided into three cases: crustal, interplate, and intraplate. In this study, the parameters corresponding to the case where S is crustal were used because the formula is more concise in this case and is also suitable for the prediction of ground shaking values for shallow source earthquakes. d denotes the regression coefficient, ε denotes the standard deviation, X is the shortest distance to the fault (km), and Rh denotes the distance between the grid point and the nearest aftershock point. Zh denotes the depth of this point and was set to 1 in the calculation. Amp is the PGV site correction coefficient, and V_S30 is the equivalent shear wave velocity at a depth of 30 m below the surface.

Figure 3 depicts the main steps for producing a quick seismic intensity assessment map, which is simple and repeatable. All the input data for GMPE can be obtained using ArcGIS spatial analysis functions (such as buffer, create fishnet, near and spatial join). Maps of earthquake intensity assessment results can also be created using ArcGIS software. R software can be used to perform pre-processing and local weighted regression operations on raw aftershock data.

Figure 3. Workflow of quickly drawing seismic intensity assessment map.

3. Results

3.1. Regression results of geographic coordinates of aftershocks

The IQR method was used to examine aftershock data within 2 h of the Luding earthquake, and seven aftershock events with abnormal spatial distribution were removed (Figure 4(a)). The processed aftershock sequences revealed an NW-SE strike distribution trend, which was consistent with the overall strike of the Xianshui River Fault Zone at this location (Han et al. 2022). This demonstrates that the spatial distribution of this part of the aftershocks reflected the spatial information of the earthquake’s rupture. Regression of the pre-processed aftershock data with rLowess produced the curve shown in Figure 4(b), which emphasizes the length and direction information among the aftershock data, indicating that the Luding earthquake showed a trend of bilateral rupture. The total length of the measured fitted curve was approximately 35.42 km, which was approximately 12.09 and 23.33 km long in the NW and SE directions of the epicentre, respectively. This indicates that the earthquake ruptured primarily in the SE direction of the epicentre, which is consistent with the source rupture inversion results released by the Institute of Geophysics, CEA (https://www.cea-igp.ac.cn/cxdt/279406.html). This earthquake’s surface rupture, computed using the Wells and Coppersmith (1994) surface rupture formula, was approximately 29.13 km, which is somewhat shorter than the fitted value. The earthquake rupture was further reflected by the aftershocks selected using buffers, meaning that the earthquake clearly ruptured in the SE direction, and the aftershocks in this direction were distributed within approximately 22 km from the epicentre. Additionally, a slight rupture may exist on the NW side of the epicentre, and the aftershocks in this direction were distributed within approximately 7.8 km from the epicentre. The range of aftershock spatial distribution was comparable to that determined by the near empirical formula, indicating that the aftershocks chosen by the buffer zone may accurately represent the earthquake rupture state. The number of chosen aftershocks increased with the buffer’s radius set to 1.5 km, primarily on both sides of the fitted curve. The range of spatial distribution did not vary substantially compared to the screening findings with the 1-km buffer’s radius.

Figure 4. Results of processing aftershock data within 2 h of the 2022 Luding Mw6.6 earthquake in Sichuan: (a) data pre-processing to remove aftershocks with abnormal spatial distribution; (b) rLowess regression results; (c) aftershocks selected using the regression curve’s 1-km radius buffer; and (d) aftershocks selected using the regression curve’s 1.5-km radius buffer.

3.2. Seismic intensity assessment results

We calculated the PGV and PGA values of the Luding earthquake using pre-processed aftershocks, rLowess regression results, and aftershocks within the 1-km and 1.5-km radius buffers of the regression curve, respectively, and compared the results with the macroscopic seismic intensity survey map released by the CEA (Figure 5). The seismic intensity results assessed directly using pre-processed aftershocks were consistent with the assessed values of the highest intensities of the isoseismic lines drawn after the field survey. Furthermore, the ranges of the two were relatively compatible, especially where the regional boundaries of degrees VI and VII overlapped with the isoseismic lines, but some aftershocks in the pre-processed data may interfere with the assessment result. The intensity ranges of the evaluations of degrees VIII and IX were substantially greater than the findings of the field survey, as shown in Figure 5(a). The seismic intensity assessed using pre-processed aftershock data can provide a preliminary understanding of the damage in the seismic region, but an accurate determination of the extent of the hardest hit areas can help mitigate casualties and property damage during the emergency ‘black box’ period, so further processing of the pre-processed data is required.

Figure 5. Map of seismic intensity assessment for the 2022 Luding Mw 6.6 earthquake in Sichuan Province, China. (a) The seismic intensity was calculated using pre-processed aftershocks; (b) seismic intensity was calculated using rLowess results; (c) seismic intensity was calculated using aftershocks within 1 km of the rLowess regression curve; and (d) seismic intensity was calculated using aftershocks within 1.5 km of the rLowess regression curve.

The range of seismic intensity VIII and IX that was assessed using the rLowess regression findings was substantially smaller, and the VIII area was more compatible with the isoseismic range in the SE direction of the epicentre. However, the rLowess regression results were impacted by aftershocks at both ends of the fault, the regression curves frequently exceeded the actual rupture, and the seismic intensities evaluated using the results frequently led to incorrect estimations of the size of the hardest-hit area when the magnitude of the earthquake was relatively small (Zhao et al. 2022b). The extent of the hardest-hit regions determined by aftershocks inside the buffer zone was more converging and matched the epicentre’s southeast direction. The intensity measured in the long-axis direction was more skewed toward the NW direction than the isoseismic line, and most of the assessment findings of the IX degree were found within the IX-degree range of the isoseismic line. The seismic intensity values assessed by the aftershocks within the 1-km and 1.5-km radius buffers did not substantially differ. When creating intensity maps using field surveys, several aspects were considered, including geological phenomena, population distribution, and damage in the disaster region. Therefore, it is reasonable that the isoseismic range was wider than the intensity assessment range of VI and VII degrees.

The PGA, PGV, and PGV_VS30 data derived from aftershocks inside the 1.5-km radius buffer zone and the station observation data were used to compute the log residuals. The average residuals of PGA and PGV, which are shown in Figure 6(a) as being extremely close to zero at 0.058 and 0.096, respectively, suggest that the ground shaking data derived using this approach were accurate. The distribution of residuals between the site-corrected PGV, PGV_VS30, and PGA data is shown in Figure 6(b). After site correction, the residuals of PGV values increased significantly, representing the large drop in PGV values after correction, suggesting that the site in the region mostly contributes to lowering PGV. As shown in Figure 7, the size of the disaster area (the region bigger than VI degrees) as determined by the PGA computed values was much greater than the findings of the field study. The PGV values were used in this research to assess the seismic intensity because the size of the impacted region determined by the PGV data without site correction was reasonable and more consistent with the findings of the field survey than the PGA results.

Figure 6. Residual distribution of predicted ground shaking values: (a) residual distribution of PGA and PGV; (b) residual distribution of PGA and PGV_VS30.

Figure 7. Seismic intensities were evaluated using ground shaking prediction values in two ways: (a) using PGA-predicted values and (b) using PGV-predicted values.

The Sichuan Basin, China, elliptical intensity attenuation relation was used to determine the intensity estimates for this earthquake (Lei et al. 2007). The rupture length was derived using the empirical Wells surface rupture relationship, the rupture site was assumed to be to the right of the epicentre, and the rupture direction was compatible with the rupture zone direction. This calculation considered the line source model. The results of the intensity assessment were regular ellipses with the long axis direction consistent with the spatial distribution of aftershocks, as shown in the isoseismic lines in Figure 8(a), and the highest intensity was VIII degrees, which is lower than the findings of the seismic damage field survey. All the ranges of IX degrees we assessed were located within the VIII-degree region of the elliptical relationship calculation results, further demonstrating the accuracy of our method. The locations and ranges of each intensity region obtained using this empirical relationship were generally similar to the intensity assessment results in this study.

Figure 8. Results of other methods used to assess the seismic intensity of the 2022 Luding Mw6.6 earthquake in Chinese Sichuan Province. (a) Results of an evaluation of the intensity decay relationship made using data from the Sichuan Basin in China; (b) the seismic intensity result of ShakeMap. The intensity obtained in this study is shown in the bottom panel.

ShakeMap, as an important international platform for earthquake disaster information services, has applicable ground motion attenuation relationships set for different earthquake magnitudes and different regions, and constantly updates the earthquake disaster damage assessment results with the accumulation of information on the affected areas (Wald et al. 2006; Worden et al. 2020). As shown in Figure 8b, we obtained the most recent version (version 5) of the ShakeMap intensity assessment findings and compared them to our own intensity evaluation results. ShakeMap’s greatest intensity rating is VIII degrees, which is compatible with our highest intensity value when measured using the MMI scale. Statistics of the hardest-hit areas for various intensity outcomes are shown in Table 2. The area of the hardest-hit regions predicted by ShakeMap was closest to the findings of the field survey, but the location and extent diverged significantly from the actual intensity distribution. The earthquake strength projected by ShakeMap, particularly the size of the hardest-hit regions, may include significant mistakes due to incomplete seismogenic region data. The location and size of the hardest-hit regions as assessed by aftershocks in the buffer were in excellent agreement with the real scenario, and the intensity zones as assessed by the 1.5-km radius buffer selection for aftershocks did not differ noticeably from the findings of the 1-km radius buffer screening. This demonstrates that the method proposed in this study may accurately pinpoint the size of the worst-hit regions.

Table 2. The maximum values of seismic intensity and the area of hardest hit areas were obtained by different methods.

Source of intensity result	Maximum intensity	Area of hardest hit region(km^2)
Site investigation	IX	1040.66
Pre-processed data + SM99	IX	2121
rLowess + SM99	IX	1281
Buffer (Radius = 1km) + SM99	IX	1127
Buffer (Radius = 1.5 km) + SM99	IX	1193
Elliptical attenuation relation	VIII	611.96
ShakeMap	IX (CSI)	1049.77

4. Discussion

4.1. Promotion of the method

Since the Wenchuan Mw7.9 earthquake, the Mw7.3 earthquake in Maduo, Qinghai, on 21 May 2021, has been the most powerful earthquake in China. Although there were no fatalities, the earthquake severely damaged the local infrastructure and generated a substantial surface rupture that was accompanied by different surface deformation and geological phenomena. This case was invaluable for earthquake study (Jin and Fialko 2021; Yuan et al. 2022). On 10 June 2022, an Mw5.9 earthquake struck Maerkang, Sichuan, China. Numerous places were definitely felt during this swarm-type quake (Yan et al. 2022). Figure 9 displays the results of evaluating the magnitudes of these two earthquakes using the methodology outlined in this research, with a 1.5-km radius buffer zone radius. Both the highest value of X projected for the Maduo earthquake and the maximum value of VIII predicted for the Maerkang earthquake intensity agreed with the findings of the field survey.

Figure 9. Results of the seismic intensity assessment: (a) 2021 Maduo, Qinghai Mw 7.3 earthquake intensity assessment result; (b) 2022 Maerkang Sichuan Mw5.9 earthquake intensity assessment result. The place where the profile maps were drawn is shown by the solid black line.

The projected range of the VII-degree zone of the Maduo earthquake fit the actual survey data well along the long axis of isoseismic lines. All isoseismic lines of VIII degree or higher were situated inside the locations that were projected to be the most impacted by the earthquake, according to the findings of the intensity prediction. According to Figure 9(b), the VIII-degree intensity zone of the Maerkang earthquake was essentially within the bounds of the survey findings, and the ranges of VI and VII degrees in the EN direction of the epicentre were compatible with the isoseismic lines. This forecast reflects that the site’s impact on the PGV mostly exhibits an amplification effect in the EN direction of the epicentre, which is also shown by the VI-degree isoseismic lines extending in that direction.

Plotted alongside the results of the field survey are the profiles of seismic intensities determined using the PGA, PGV, and PGV_VS30 data for the two earthquakes (Figure 10). The maximum intensities determined using the PGA predicted values were one degree different from the survey data, providing a rough description of the location and extent of the hardest hit areas. The PGV-based intensity assessment findings were more accurate than those utilizing PGA in identifying the location and size of the hardest-hit locations. However, the intensity maxima may be skewed when the earthquake’s magnitude is significant. The location and extent of intensity zones assessed using the PGV_VS30 data after site correction were close to the actual damage, and the intensity assessment results were more detailed than those obtained using the other two datasets, indicating potential intensity anomaly zones, demonstrating the accuracy of the method in quickly identifying the hardest hit areas of earthquakes. These two earthquakes serve as examples of how both relatively big and relatively minor magnitudes may be used to apply the approach. In light of the previous research, we believe that it is appropriate to use the aftershocks selected by the buffer to assess seismic intensity when Mw < 7.0; when 7.5 > Mw > 7.0, it is possible to use the locally weighted regression results directly or to gauge seismic intensity using the aftershock events determined by the buffer zone analysis, which can be based on the number of aftershocks available; and when Mw > 7.5, we suggest using the regression result.

Figure 10. Seismic intensity profiles are analysed based on predicted values of ground motion versus field survey results. (a)–(c) are profiles of expected vs. real earthquake intensities in Maduo in 2021; (d)–(f) are profiles of projected versus actual earthquake intensities in Maerkang in 2022. The actual value is shown by the solid black line, while the predicted value is shown by the red realization.

The seismic intensity evaluation findings of our technique, as shown in Figure 11, are realistic, particularly when considering the location and size of the hardest-hit regions, for the three earthquakes that occurred in Japan and the Menyuan Mw 6.6 earthquake in Qinghai, China, in 2022. The intensity assessment results of ShakeMap are reliable under the condition that the seismic data of these three Japanese earthquakes are more adequate and have been corrected for a long time. The intensity ranges we assessed were very similar to the results from ShakeMap, and the average residuals of the PGA and PGV_VS30 results for these four earthquakes and the station observations, as shown in Figure 12, all lied in the range of [−0.3,0.3]. This shows that the calculated values of ground shaking were close to the observed values, further validating the stability of our method. It is logical to use PGV_VS30 data to determine the earthquake intensity because the fluctuation of PGV_VS30 in Figure 12 was steadier than that of PGA. The location and size of the heaviest impacted regions may be identified using both PGA and PGV_VS30, which both exhibited a pattern of lower residuals close to the epicentre site in comparison with other places. Therefore, it can be used to determine the location and extent of the hardest-hit areas.

Figure 11. Seismic intensity assessment results and distribution of some seismic stations: (a) 2000 western Honshu, Japan Mw6.7 earthquake; (b) 2008 eastern Honshu, Japan Mw6.9 earthquake; (c) 2018 Tomakomai, Japan Mw6.6 earthquake; (d) 2022 Menyuan, China Mw6.6 earthquake.

Figure 12. Residuals of PGA predictions and PGVVS30 predictions: (a) 2000 western Honshu, Japan Mw6.7 earthquake; (b) 2008 eastern Honshu, Japan Mw6.9 earthquake; (c) 2018 Tomakomai, Japan Mw6.6 earthquake; (d) 2022 Menyuan, China Mw6.6 earthquake.

4.2. Applicable conditions

To use the approach for quick estimation of seismic intensity described in this study, a certain number of aftershock sequences must be collected shortly after the earthquake. The basis for the use of this approach has been laid by the installation of seismic monitoring stations and the advancement of seismic monitoring technology, which has consistently enhanced the amount and quality of the detected aftershock events. In our previous study, using 40–100 aftershock data yielded good weighted regression results, and the seismic intensity assessment time could be significantly shortened with enough aftershock data (Zhao et al. 2022b). However, whether the regression results match the information of the mainshock should be considered, such as whether the length of the fitted curve is consistent with the rupture scale at that magnitude, and some of the aftershock data can be removed when necessary with expert experience. The emphasis of this research is on earthquakes of Mw <7.0, particularly those with Mw close to 6.5. This size range of earthquakes may not result in surface rupture on a substantial scale, and geophysical methods used to quickly anticipate surface rupture processes may have major inaccuracies in their findings. When the magnitude is sufficiently high, the mainshock fault releases enough stress, and the seismic monitoring circumstances allow for the potential that few aftershocks are recorded in the first two hours following the earthquake, and the rupture scale at this magnitude may reach over 100 km (Wells and Coppersmith 1994). The 2016 Mw7.8 earthquake that struck New Zealand is shown in Figure 13 along with regression findings and buffers. The rupture length for this earthquake was determined to be 197.28 km using the Wells empirical formula for surface rupture at this magnitude. The buffer of 1.5 km radius had a very small extent at this rupture scale, so the number of selected aftershocks was also small. When the radius was set to 5 km, the number of selected aftershocks increased significantly, but the relationship between the multipart aftershock events and the mainshock rupture within this range may not be accurate (Zhao et al. 2022b).

Figure 13. Results of aftershock selection within 2 h after the 2016 Mw7.8 earthquake in New Zealand: (a) buffer radius set to 1.5 km; (b) buffer radius set to 5 km.

It is evident that using the proper GMPEs to compute ground shaking data may provide more accurate findings. Using merely the epicentre position, the SM99 formula utilized in this study can also compute ground shaking data, and the findings of the intensity evaluation that is produced are circular in shape. By integrating the approach suggested in this work, dynamic mapping of seismic intensity assessment data may be accomplished with the dynamic monitoring of aftershock occurrences.

5. Conclusions

Herein, we presented a technique to quickly map seismic intensity assessment based on geographic information included in the geographical distribution of early aftershock sequences during the ‘black box’ phase of earthquake emergency response. The geographic coordinates of the pre-processed early aftershocks were first regressed using the rLowess method to obtain a curve that can represent its spatial distribution trend, then a buffer of a certain radius was added to this curve to select the aftershock data for the GMPE calculation, and finally, the seismic intensity map was drawn based on the ground motion data calculation results. The contributions of the article are as follows:

1. RLowess is a nonparametric approach for analysing complicated physical pro-cesses without a theoretical model base. Early aftershock sequences are caused by a complicated geophysical process. We employed rLowess to perform regression analysis on the geographic coordinates of aftershocks, yielding a regression curve that roughly reflects the rupture information of the mainshock, simplifying the consideration of the elements influencing the spatial distribution of aftershocks, and quickly obtaining reliable and reasonable regression results, providing a new verification of the rupture results obtained using geophysical means.

2. The utility and effectiveness of data during the ‘black box’ period of earthquake emergency response was improved, and to some degree, information connected to the rupture of the mainshock contained in early aftershock geographic coordinates was mined. When compared to rLowess regression curves, early aftershock sequences selected by buffers could more accurately describe the length and direction information of the mainshock rupture, and seismic intensity assessment maps based on these aftershocks could accurately determine the location and extent of the hardest hit areas.

3. The application of geographic information technology in earthquake disaster information service was developed, and research results in seismology and solid geophysics were integrated with GIS spatial analysis tools, addressing the problem of large deviation of seismic intensity assessment results from the real damage situation in the case of small moment magnitude. This approach is easy to use and can rapidly generate a seismic intensity assessment map, which aids the earthquake disaster information service and provides more scientific data for early earthquake emergency response.

4. This work applied current research findings on the association between early aftershocks and seismogenic faults, and it validates the findings of previous studies to some degree. Simultaneously, the issues met using the method will give fresh ideas and entrance points for research into the interaction between aftershock sequences and mainshock rupture.

It should be noted that the goal of this research, which proposes a quick seismic intensity assessment map based on aftershock data, is to enrich the methods for evaluating the level of damage during the ‘black box’ phase of earthquake emergency response. However, the geophysical process of generating aftershocks and the relationship between aftershock sequences and causative faults are complex, and thus need to be cross-referenced with the disaster loss results assessed by geophysical means to assess loss more accurately. To acquire more precise intensity evaluation findings in a certain area, GMPE for that region must be established. In addition, other regression analysis methods can be tried for processing and analysing seismic geological data.

The authors would like to thank anonymous reviewers for their helpful comments on this article.

Author’s contributions

WC conceptualised the project, acquired funding, and supervised the project. HZ performed the investigation, deployed the software and code, and edited the manuscript. WC, YJ and HZ developed the methodology and revised the manuscript. CZ and DK provided input and assistance in the improvement of the methodology.

Data availability statement

The data that support the findings of this study are available from the corresponding author, Wenkai Chen, upon reasonable request.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This research was funded by the Fundamental Research Funds in the Institute of Earthquake Forecasting, China Earthquake Administration, grant number 2020IESLZ05, and the National Key Research and Development Program of China, grant number 2017YFB0504104.