A real-time prediction model for macroseismic intensity in China

The macroseismic intensity spatial distribution is an important input for most rapid loss modeling and emergency work. Data from a total of 175 earthquakes (Ms ≥ 5.0) in China from 1966 to 2014 were collected, and the rapid assessment method of macroseismic intensity distribution was studied. First, simple relationships among the epicentral intensity, magnitude, and focal depth were established. A greater amount of database is used in this study than that in a previous work (Fu and Liu in Sci R 4(5): 350-354 (1960), Mei in Chin J Geophys 9(1): 1–18 (1960), and Yan et al. in Sci Chin 11: 1050-1058 (1984)), and the studied earthquakes all occurred in the last 50 years, providing more accurate and uniform parameter information. As the seismic intensity-attenuation relationship is traditionally used to estimate the intensity distribution, the macroseismic intensity-attenuation relationship for mainland China was fitted by the earthquake data collected in this region. The deviation of the intensity assessment by the macroseismic intensity-attenuation relationship was examined for 43 earthquakes (Ms ≥ 6.0). In addition, seismic damage emergency assessment work requires the isoseismal lines to be constantly modified according to the updated information. Therefore, an improved ellipse intensity-attenuation model was proposed in this study, completed by the establishment of a semimajor axis and semiminor axis length matrix. Based on the initial value of the length matrix obtained by the regression of historical data and survey data from the site, the least mean squares (LMS) algorithm is used to revise the length matrix. In the end, the practicability of this method is verified by a case study of the Lijiang 7.0 earthquake.


Introduction
Macroseismic intensity has traditionally been used worldwide as an index to describe the damage effects of earthquakes Vitor and Nick 2019), especially in regions where little or no measured strong motion data is available. As macroseismic intensity could be related directly to the damage degree, most rapid loss modeling and emergency work conducted around the globe use macroseismic intensity as an input Kamer et al. 2009;Trendafiloski et al. 2011). There are three main ways to achieve the rapid assessment of macroseismic intensity after an earthquake: by the relationship between macroseismic intensity and ground motion parameters at strong-motion stations (Bilal and Askan 2014;Boatwright et al. 2001;Wu et al. 2003), by the macroseismic intensityattenuation relationship (Allen et al. 2012;Saman 2018; Chen and Liu 1989), and by simulating the focal mechanism and propagation effect with a seismological method (Graves and Pitarka 2010;Rhie et al. 2009). The first two methods are the most commonly used so far. Numerous attempts have been made to present the relationships between macroseismic intensity and engineering ground motion parameters in the past. As the peak ground acceleration (PGA) is significant in seismic design, researchers have related the macroseismic intensity to PGA (Gutenberg and Richter 1956;Atkinson and Kaka 2007;Worden et al. 2012). However, the peak ground velocity (PGV) is related to the kinetic energy, which further influences the damage to a structure. Studies have demonstrated that high intensities (> VII) correlate best with the PGV (Kaka and Atkinson 2004;Saman et al. 2011;Wald et al. 1999). The US Geological Survey (USGS) adopts the relationship between macroseismic intensity and peak ground acceleration and velocity proposed by Wald (Wald et al. 1999), which is based on eight historical earthquakes in California. The intensity is calculated by the PGV when the intensity is higher than VII, while the intensity is calculated by combining PGA and PGV when the intensity is between V and VII. The average distance among ground motion stations in Japan is approximately 10 km. The Japan Meteorological Agency (JMA) and National Research Institute for Earth Science and Disaster Resilience (NIED) calculate the intensity using three-   1996). At present, the average distance among ground motion stations in China is approximately 90 km, mainly concentrated in the north-south seismic belt. The total number of ground motion stations is currently not sufficient for the rapid assessment of macroseismic intensity. However, historical earthquake data are very abundant in China. The rapid assessment of macroseismic intensity based on the macroseismic intensityattenuation relationship has a wide application in China (China Earthquake Administration 1992; Wang and Yu 2008;Cui et al. 2010). This method needs to know the macro-epicentral intensity and its location, the direction of the semimajor axis, the value of the semimajor axis, and the semiminor axis radius length to complete the isoseismal map. In addition, the empirical isoseismal map based on the attenuation relationship is expected to be modified in real time as the actual intensity is surveyed by field investigation, social software, and so on. For this purpose, an improved method based on the macroseismic intensity-attenuation relationship is proposed for the rapid assessment of seismic intensity, and a real-time correction method for isoseismal mapping using site investigation is studied.

The center and direction of isoseismal map
The location of instrumental epicenter could be determined in a few minutes based on the seismic monitoring network. However, many earthquakes show that there is a large deviation between the instrumental epicenter and macro-epicenter. The instrumental epicenter is the projection on the ground from the initial rupture, while the macro-epicenter is the most damaged area. Possible location of macro-epicenters could be acquired by the analysis of the relationship between the instrumental epicenter and nearby faults. Li (2000)  133 earthquake main shocks in China to analyze the deviation distances between the macro-epicenter and the instrumental epicenters. The study showed that the deviation distances of 88% of the earthquakes are below 35 km and that those of the others are from 35 to 75 km. A practical method to identify the most possible location of the macro-epicenter was given.
It is generally recognized that most strain energy is released by the main shock, and the residual energy is sequentially released by the aftershocks. The concentrated region of aftershocks may suggest the location of the macro-epicenter (Yan 2010). As the rupture develops along a fault, the aftershocks extend outward. In general, the distribution range of the aftershocks within 24 h is much larger than the meizoseismal area. A limited number of aftershocks occur within 2 h, but these aftershocks concentrate around the meizoseismal area. The distribution of aftershocks within 4 or 8 h could significantly reflect the meizoseismal area. Therefore, the distribution of aftershocks within 4 or 8 h could be used to estimate the location of the macro-epicenter (Yan 2010).
Intensity data are in agreement with major-minor axes. There are three main ways to define the direction of the semimajor axis: (1) The tectonic map of the major active faults could be used to speculate this direction. When there is only one active fault passing through the epicenter, the direction of the fault is considered to be the direction of the semimajor axis. (2) Using the focal mechanism solution to define the direction of the semimajor axis is another method (Ekstrom et al. 2012;Takaki 2008). (3) The distribution of aftershocks basically aligns with the direction of the semimajor axis. Therefore, in actual work, we can consider the direction of the active faults, the focal mechanism solution, and the distribution of the aftershocks to make a comprehensive judgment of the direction of the semimajor axis.

Database
The data used in this study consisted of 175 earthquakes that occurred in China between 1966 and 2014 with a surface-wave magnitude (M s ) varying from 5.0 to 8.0. A list of the location, epicentral intensity, and magnitude for each earthquake is summarized in Table 1. For earthquakes occurring in China, M s is largely used in size estimation because of the distribution of the seismic stations. Moment magnitude (M w ) is directly connected to earthquake source processes, does not saturate, and thus provides the most robust estimate of the magnitude of large earthquakes. It has been known that there is no significant difference between the two scales for earthquakes with magnitude ≤ 4.5 in China (Liu et al. 2007). Tang et al. (2016) derived an empirical M w -M s conversion in China for events with magnitude > 4.5. However, the M w results converted by Tang's relationship are Table 2 The coefficient of determination (R 2 ) and standard error results of Eqs. (4)- (8) Equation (4) E q u a t i o n( 5) E q u a t i o n( 6) Equation (7) Equation ( The isoseismal maps were obtained by detailed postearthquake survey and site evaluations after the earthquakes in the meizoseismal areas by reconnaissance teams. All the maps used in this study were prepared by the China Earthquake Administration (CEA). The dataset of macroseismic intensity employed in this study comprises the data reported in Zhang (1988aZhang ( , b, 1990Zhang ( , 1999Zhang ( , 2000, Chen (2002aChen ( , b, c, 2008 and Jiang (2014Jiang ( , 2018a. Moreover, in this study, the macroseismic intensities of earthquakes that occurred between 2012 and 2014 were obtained from the CEA website.

Simple relationships among epicentral intensity, magnitude, and focal depth
The rapid assessment of the epicentral intensity has great significance in earthquake relief work. Figure 1 shows the correspondence of the magnitude versus epicentral intensity with events' dots at different colors for different depths. The trendlines generally suggest that deeper events give lower epicentral intensity. But the trendline with depth > 20 km has intersection points with other lines, which does not satisfy the conventional relations. We found that Epicentral intensity Equation (5) Equation (6) Equation (7) Equation ( Equation (5) Equation (6) Equation (7) Equation (8) Gutenberg-Richter, 1956 Li, 1957Lu et al., 1981 Magnitude Magnitude Epicentral intensity Fig. 5 Comparison of the relationship between magnitude and epicentral intensity in this study and other studies. a Comparison between Eqs. (5), (6), (7), and (8). b Comparison between this study and previous studies there are only 6 events with epicentral intensity ≥ 8 in the total 32 events with depth > 20 km. And the slope of the trendline is influenced by the 6 individual events. Figure 2 shows the distribution of the focal depth. The focal depth of earthquakes occurred in mainland China is mainly distributed between 5 and 20 km. The expected relationships among magnitude, epicentral intensity, and focal depth are shown in the following equations.
where M s is the surface-wave magnitude; I e is the epicentral intensity; h is the focal depth; and a, b, and c are the fitting coefficients.
Under a wide range of conditions, earthquakes satisfy G-R frequency-magnitude scaling given by (Gutenberg and Richter 1954): where N(≥ m) is the cumulative number of earthquakes with magnitudes greater than m occurring in a specified area and time window, and a and b are constants. This relation is valid for earthquakes both regionally and globally with magnitudes above some lower cutoff m c , which defines the completeness of the catalog. Figure 3 shows the distribution of the magnitudes in the prepared database. This uneven distribution will cause the magnitude range with more data points to have a higher weight when the ordinary least squares method is used. Therefore, the weighted least squares method is used for regression fitting in this study. The weight values of 1, 1.72, 3.86, 5.06, 13.50, and 20.25 are assigned to magnitudes between 5.0-5.5, 5.6-6.0, 6.1-6.5, 6.6-7.0, 7.1-7.5, and 7.6-8.0, respectively. The results are shown in Eqs. (4)-(5) and Fig. 4. Furthermore, the database is divided into three classes by the focal depth (h): 0 ≤ h ≤ 10 km, 10 < h ≤ 20 km, and h > 20 km. The weighted least squares method is also used for regression fitting, and the results are shown in Eqs. (6)-(8). The coefficient of determination (R 2 ) and standard error are presented in Table 2. As mentioned in Fig. 1, there is only 6 events with epicentral intensity ≥ 8 in the total 32 events with depth > 20 km. And the slope of the trendline is influenced by the 6 individual events. So, these 6 events were removed to obtain Eq. (8). Thus, (Eq. 8) is only suitable for events with epicentral intensity ranging from 6 to 7.
The simple predictive relationships among the magnitude, epicentral intensity, and focal depth in this study were compared with similar relationships developed in the following studies: Fu and Liu (1960), Mei (1960), and Yan and Guo (1984). Similarly, the relationship between the magnitude and the epicentral intensity was compared with those from studies by Gutenberg and Richter (1956), Li (1957), and Lu et al. (1981). The corresponding equations derived by these researchers for comparison with this study are given in Eqs.  Figure 5a shows the comparison between Eqs. (5) and (8). The results show that deeper events give lower epicentral intensity. Figure 5b shows the comparison between the results of this study and those of Eqs. (12)-(14). Equation (5) is mainly located in the middle of the regression lines of Eqs. (12)-(14) and is most similar to the study by Gutenberg and Richter (1956).

Deviation of intensity distribution evaluated by seismic intensity-attenuation relationship
Howell and Schultz (1975) derived the functional relationship of seismic intensity attenuation as in Eq. (15), based on seismology.
where I is the intensity; R represents the distance from the epicenter; M denotes the magnitude; R 0 is the presupposed constant; ɛ denotes the uncertainty, the mean of which is zero; and A, B, C, and D are the regression constants. Musson (2005) obtained the seismic intensityattenuation relationship in the UK based on Eqs. (15) and (16), as shown in Eqs. (17) and (18).
Taking east longitude 105°as the bound, the earthquake intensity zoning map of China gives the seismic intensity-attenuation relationship of the eastern and western regions in China as Eq. (20)-(23). Equations (20) and (21) denote the attenuation relationships in the directions of the semimajor axis and semiminor axis in eastern China; Equations (22) and (23) denote the attenuation relationships in the directions of the semimajor axis and semiminor axis in western China. I a ¼ 6:046 þ 1:480M −2:081ln R þ 25 ð Þ ð20Þ I a ¼ 5:643 þ 1:538M −2:109ln R þ 25 ð Þ ð22Þ Based on the collected earthquake data and Eq. (19), the seismic intensity-attenuation relationships in the directions of the semimajor axis and semiminor axis in mainland China are established by the weighted least squares method and are given as Eqs. (24) and (25). A comparison of these results with the earthquake intensity zoning map of China is shown in Figs. 6 and 7. The coefficient of determination (R 2 ) and standard error are presented in Table 3. I a ¼ 6:1709 þ 1:3343M −1:9119ln R þ 30 ð Þ ð24Þ According to Eq. (5), the epicentral intensities for M s = 7.0, 6.0, and 5.0 are 9.1, 7.2, and 5.4. Compared with the results shown in the earthquake intensity zoning map of China, the results from this study present a more reasonable relationship between magnitude and epicentral intensity in the near-field. However, in the far-field, the results in this study decay more slowly. This may be because the data of more recent earthquakes were adopted in this study. In recent years, the evaluated disaster area is often larger than the actual disaster area. The area bounded by isoseismal line with intensity I = I i , could be regarded as the damage area with intensity ≥ I i . The area ratio, defined as the statistical model evaluation results based on Eqs. (24)-(25) divided by the actual data is chosen as the evaluation index to study the deviation of the intensity distribution evaluated by the seismic intensity-attenuation relationship. The compared results are shown in Table 4.
The statistical results with intensity values greater than 10 are smaller than the actual area. As the seismic data with intensity values greater than 10 are very limited, they are not listed in Table 4. The maximum, minimum, median, and average area ratios with intensity ≥ I i are shown in Table 5. The ratio between the maximum value and the minimum value could be at hundreds times, which denotes the large uncertainty in the earthquake data. However, the variation tendencies of the median and average values are basically consistent. The assessment area is always smaller than the actual area with high seismic intensity, while the assessment area is always larger than the actual area with low seismic intensity. This is possibly because all the isoseismal lines are derived by only the three regression coefficients in the attenuation Eq. (19). However, the only three coefficients are not enough to reveal an isoseismal map with several intensities.
6 The intensity attenuation in the matrix model and correcting algorithm 6.1 The intensity attenuation in the matrix model An improved seismic intensity-attenuation relationship based on the ellipse model is proposed in this study, which is completed by the establishment of a semimajor  axis and semiminor axis radius length matrix. Based on the initial value of the length matrix obtained by the regression of the historical data and survey data from the site, the least mean squares (LMS) algorithm is used to revise the length matrix and draw the intensity isoseismal lines, which is called the intensity attenuation in the matrix model. Based on the collected earthquake data shown in Table 1, the relationships between the radius length and the seismic magnitude are given in Tables 6, 7, and 8, which could be used to initialize the semimajor axis and semiminor axis radius lengths for the isoseismal lines.

Correcting algorithm
The least mean squares (LMS) algorithm (Javed and Ahmad (2020)) is an adaptive filtering method for solving the online estimation problem that is modeled by: where A n is n × N(n ≫ N) data matrix of rank r ≤ N, formed by a sequence of input signals θ f n ð Þg ∞ n¼1 . The vector b n ∈ R n consists of desired signals, and x n = [x n (0), …, x n (N − 1)] T ∈ R N is tap-weight vector of length N. The  objective is to predict desired signal s(n) with an online learning algorithm at time n, by estimated output signal.
where Θ n = [θ(n)θ(n − 1)…θ(n − N + 1)] T is the input vector at instant n. So that the error incurred is given by This process requires an adaptive algorithm to update filter tap-weight vector x n recursively as new signal comes in. The LMS algorithm does so by the update equation.
where η is the step-size parameter and controls the convergence speed and stability of the algorithm.
Then, based on the intensity survey data from the disaster area, the LMS algorithm could be used to revise the length matrix after the earthquake. The rules are as follows: (1) The solid lines in Fig. 8 are the initial evaluated isoseismal lines. For a survey point with intensity of I, if it is located at the position of point 1, the isoseismal lines do not need to be revised; if it is located at the position of point 2, the isoseismal line of intensity I + 1 needs to be revised; and if it is located at the position of point 3, the isoseismal line of intensity I needs to be revised.
(2) Based on the semimajor axial and semiminor axial projected length of the survey point and LMS algorithm, the radius length matrix could be revised. Here, the semimajor axis radius length of the isoseismal line of intensity I is taken as an example. The modification principle is shown in Eq. (30) based on the LMS algorithm. Taking 0.01 as the calculation step, a total of 101 operations were conducted for η values ranging from 0 to 1. When all the intensity survey points are substituted  into the modified procedure, the macroseismic intensity result with minimum total error is chosen as the modified isoseismal lines. And the total error is defined as the sum of the errors of all the intensity survey points.
where R Ia is the modified semimajor axis radius of the isoseismal line of intensity I, R ' Ia is the uncorrected semimajor axis radius of the isoseismal line of intensity I, R * Ia is the semimajor axis radius length of the survey point located in the ellipse, and η is the learning speed (the value of η ranges from 0 to 1).

Example
A case study of the Lijiang 7.0 earthquake (February 3, 1996) was conducted to verify the practicability of the modified method. The Yulong fault is the main causative structure of this earthquake. The direction of the semimajor axis is nearly aligned with the direction of the Yulong fault. As mentioned in section 2, in actual work, we can consider the direction of active faults, the focal mechanism solution, and the distribution of the aftershocks to make a comprehensive judgment of the direction of the semimajor axis. As the updated correction is the focus of this section, the orientation of the axes was specified based on the known facts. The 20 intensity survey points in Fig. 9 were selected at random to verify the practicability of this modified method.
Based on the initial radius length matrix in Table 7 and the modified method presented in this paper, the radius length with different amounts of intensity survey data could be calculated, as shown in Tables 9 and 10. Then, the isoseismal lines determined by considering different numbers of intensity survey points are drawn in Fig. 10. The area of isoseismal line with each intensity could be calculated by the radius length, as shown in Table 11. Table 12 shows the area errors, which is defined as Eq. (31).
where A i is the evaluated area, and the A is the actual one. With increasing number of survey points, the isoseismal lines determined by using this method become more similar to the actual ones.

Conclusions
The simple relationships among epicentral intensity, magnitude, and focal depth were established based on 175 earthquakes (M s ≥ 5.0) in China. Compared to the previous work, this study includes a greater amount of database, and the studied earthquakes all occurred in the last 50 years, resulting in a more accurate and uniform parameter information. The seismic intensity attenuation for mainland China was fitted. The deviation of intensity assessment by intensity-attenuation relationship was examined by 43 earthquakes (M s ≥ 6.0). The results show that the assessment area is always smaller than the actual area with high seismic intensity, while the assessment area is always larger than the actual area with low seismic intensity. This is possibly because all the isoseismal lines are derived by only the three regression coefficients in the attenuation equation. However, the only three coefficients are not enough to reveal an isoseismal map with several intensities.
An improved ellipse intensity-attenuation model, which is completed by the establishment of a semimajor axis and semiminor axis length matrix, was proposed in this paper. Based on the survey points and LMS algorithm, the isoseismal lines could be modified in real time. The result of the Lijiang 7.0 earthquake example shows that the proposed method is simple and practical in emergency assessment work.