Real-time numerical shake prediction and updating for earthquake early warning

Ground motion prediction is important for earthquake early warning systems, because the region’s peak ground motion indicates the potential disaster. In order to predict the peak ground motion quickly and precisely with limited station wave records, we propose a real-time numerical shake prediction and updating method. Our method first predicts the ground motion based on the ground motion prediction equation after P waves detection of several stations, denoted as the initial prediction. In order to correct the prediction error of the initial prediction, an updating scheme based on real-time simulation of wave propagation is designed. Data assimilation technique is incorporated to predict the distribution of seismic wave energy precisely. Radiative transfer theory and Monte Carlo simulation are used for modeling wave propagation in 2-D space, and the peak ground motion is calculated as quickly as possible. Our method has potential to predict shakemap, making the potential disaster be predicted before the real disaster happens. 2008 MS8.0 Wenchuan earthquake is studied as an example to show the validity of the proposed method.


Introduction
Earthquake early warning (EEW) is used to prevent disaster before the strong ground shakes arrive. Ground motion prediction is an important task for EEW systems, and the peak ground motion [e.g., peak ground velocity (PGV)] is often predicted by EEW systems. If PGV is predicted to exceed a threshold, alarms can be triggered to prevent the potential disasters.
Based on the method of estimating source parameters, traditional methods can be classified as front detection, onsite, and network methods (Allen et al. 2009). Front detection methods are often used to protect the city off the coast by installing one or more stations along the coast, such as Istanbul EEW system proposed by Alcik et al. (2009) along the northern shore of the Marmara Sea. The front detection methods trigger a warning if the ground shaking exceeds a certain level, but cannot predict the intensity information for a specific region, thus cannot provide the severity order of the disaster.
On-site methods and network methods use P wave information of one or more stations to determine the epicenter and magnitude. They calculate peak ground motion by applying ground motion prediction equation (GMPE) of magnitude and distance between epicenter and the specific location. The amplitude of P wave (Wu and Kanamori 2005;Kanda et al. 2009), the dominant ground motion period, and average ground motion period (Nakamura 1988;Allen and Kanamori 2003;Wu et al. 2007) are proposed to predict the upcoming peak ground motions. The systems have operated in many countries such as Japan for the general public (Hoshiba et al. 2008;Hoshiba and Ozaki 2014) and have been investigated in the USA, Taiwan of China, the European Union, Turkey, and other countries and regions (e.g., Alcik et al. 2009;Hsiao et al. 2009;Gasparini and Manfredi 2014;Kuyuk et al. 2013).
Traditional methods fail to predict the peak ground motion when source parameters are predicted error. The EEW system provided by Japan Meteorological Agency (JMA) (Hoshiba et al. 2008) failed to predict the peak ground motion of the 2011 M W 9.0 Tohoku earthquake (Hoshiba et al. 2011). The earthquake with large fault rupture leads to underprediction of ground motions, and multiple simultaneous aftershocks may lead to overprediction of peak ground motion (Hoshiba and Aoki 2015).
In order to tackle the failure of traditional methods above, some approaches have been proposed for estimating the real-time ground motion without using the source parameters. In 2013, Hoshiba (2013) applied boundary integral equation to predict ground motion in real time. The method only requires the incoming seismic wave directions of front stations. The amplitude and the phase information of seismic waves should be assumed, leading to the limited performance. Hoshiba and Aoki (2015) investigated numerical shake prediction method to predict real-time ground motion for EEW system on dense seismic networks. Energy distribution is used to represent the amplitude of seismic wave, which removes the phase information of seismic waves. Energy propagation is represented as radiated transfer theory (RTT) (Sato et al. 2012). Data assimilation technique (Kalnay 2003) is used to estimate the current seismic wave energy distribution. The computation time of the method is fast enough for realtime estimation of ground motion. It is shown that the performance of numerical shake prediction method on M W 9.0 earthquake is good with the large fault rupture and complex source mechanism.
The real-time ground motion estimation method and traditional method are complementary to each other. The traditional method is faster, and real-time ground motion estimation method is more precise. In this article, we propose a ground motion prediction and updating framework, for determining the peak ground motion for a specific region as quickly and precisely as possible. We design a prediction and updating scheme for fast and accurate ground motion prediction. Firstly, when P wave is detected by one or more stations, we use traditional methods to calculate an initial prediction for a preliminary estimation. For the earthquake in the seismic network and the magnitude of that is not large, the initial estimation can predict the peak ground motions precisely. Secondly, we update the peak ground motion by simulating the real-time wave propagation. The purpose of the updating process is to correct the initial prediction as quickly as possible. Our updating scheme differs from the real-time ground motion prediction methods by Hoshiba's method (Hoshiba and Aoki 2015): our concern is to predict the peak ground motion of a specific region, instead of the real-time ground motion. RTT is adopted in this article for modeling the wave propagation. For efficient computation of RTT, direct simulation Monte Carlo is used to model the wave energy propagation as particles carrying specific energies in 2-D moving directions. Data assimilation technique is incorporated to calculate the ground motion in each region by using the real-time records of seismic stations. We model the wave propagation for 60 s in advance and calculate the peak ground motion. The peak ground motion predicted is shown to be stable and precise.
Our method has similar output with shakemap generation. Shakemap is used for estimating the severity of the disaster for each region, and it is often produced after the earthquakes. Our method can predict the shakemap before the earthquake ends, meaning that the areas potentially ruptured can be determined before the real shake occurs. If the shakemap can be produced before the peak ground motion comes, the potential disaster could be avoided.
We apply our method to 2008 M S 8.0 Wenchuan earthquake and analyze the results of initial prediction and realtime shake prediction based on intensity updating. The results show that our method is robust to the small hypocenter location errors determined by EEW systems and can predict the peak ground motion precisely and fast in a real-time updating manner.
2 Numerical shake prediction and updating method 2.1 The framework of numerical shake prediction and updating In this section, we describe the proposed framework of numerical shake prediction and updating. The ground motion predicted in this article is the maximum intensity for a specific region during the earthquake, which indicates the potential damage of that region. The framework is divided into two steps, as Fig. 1 shows. Firstly, the peak ground motions in different areas are calculated after the determination of earthquake parameters of seismic stations. Epicenter location can be calculated using the relationship of the travel time difference, and the magnitude of the earthquake could be estimated using the first several seconds of the P wave information. Given the hypocenter location and the magnitude of the earthquake, the peak ground motion can be calculated using the GMPE. The initial shakemap could be estimated after P wave detection.
Second, we use a real-time ground motion updating method to predict the upcoming ground motion. The updating process is designed to alleviate the failure of initial prediction caused by wrong hypocenter parameter, underestimation of earthquake with large magnitude, and so on. Wrong estimation or underestimation for large earthquake will result in potential disasters. The real-time updating method predicts the ground motion by wave propagation modeling. The ground motion is derived from data assimilation and wave propagation. The energy of wave propagation is modeled by radiative transfer theory given the real-time records of seismic stations. In order to compute efficiently, the Monte Carlo method is used for modeling the wave propagation.

The 2-D space model in the framework
For efficient computation, 2-D space model is used for modeling the wave propagation. The objective of our method is to calculate the peak ground motion for the region of interest. Let x, y be the space model with horizontal parameter I 1 , vertical parameter I 2 , and with its corresponding increments Dx, Dy, respectively. Figure 2 shows the space model, where the region of interest is divided into grids with identical size. The ith grid is calculated as i = I 1 9 (i 2 -1) ? i 1 , where i 1 = 1, 2,…,I 1 , i 2 = 1, 2,…,I 2. The total number of grids is I = I 1 9 I 2 . The center of each grid is used to represent the location of the grid, and the location at the ith component of U n is x i = ((i 1 -0.5) Dx, (i 2 -0.5)Dy). The observations of the stations at time n are represented as a vector y n = [y 1 n , y 2 n ,…,y Y n ], where Y is the number of stations in the region. Our task is to predict the peak ground motion of U n as quickly and precisely as possible.

Initial ground motion prediction based on GMPE
The initial prediction of ground motion is estimated after the first several seismic stations' detection of P waves. The epicenter location is determined by the maximum intersection (MAXI) method (Font et al. 2004). The MAXI method is the extension of the master station method (Zhou 1994), and it determines the absolute hypocenter location based on the concept of equal differential time (EDT) surface in 3-D velocity model. In this section, we briefly describe the method. For more information, please refer to (Font et al. 2004).
The EDT surface is the surface on which the difference between the travel times of stations j and k is equal to the difference between the corresponding observed arrival times of the two stations. It is defined as the collection of all spatial points within the velocity model that satisfies the time difference between the arrivals observed at two stations, as Eq. (1) shows:  where O rj means the observation time from epicenter R to the observation station j, and C rj means the travel time from the epicenter R to the observation station j. The observation time can be calculated using P wave detection method, and the P wave travel time from any location to the station can be computed using ray path tracking algorithm within the velocity model. Thus, the EDT surface could be established with efficient algorithms. In the ideal case, the intersection of EDT surface from a different pair of stations is the location of the earthquake. For practical consideration, the MAXI method constrains the hypocenter by searching the node within the velocity model in two steps. The first step is the predetermination step for preliminary hypocenter solution and outlier(s) removal. The intersection of EDT surfaces defines a small volume where the potential hypocenter solution is. The second step fixes the final solution by minimizing the RMS of the residual variance.
After hypocenter is detected using EDT, we derive the relationship of intensity with magnitude and the distance between the region and the epicenter. The magnitude can be estimated using the relationship of magnitude with P wave statistics and the focal distance. In this article, we use P d to calculate the magnitude using the following equation (Zhang 2013): where P d is the peak displacement in 3 s after P wave detection, and M is the magnitude, D is the distance of epicenter with the station. Using more stations can calculate more stable magnitude. We convert the P d into the normalized P d where D is 10 km and take the average of the P d as the parameter P d,m to calculate the magnitude M, using the following equation (Zhang 2013): M ¼ 0:46 log 10 ðP d;m Þ À 3:74: The magnitude M can be determined by Eq. (3). After the determination of magnitude M, we calculate intensity for a region by deriving its relationship with M and the distance D between the region and epicenter. We derive the equation from their relationship with PGV. The PGV has the relationship with P d with the following equation (Zhang 2013): The PGV in Eq. (4) is composed from three-component synthesis (up-down, east-west, north-south) of seismic record. Using Eqs. (2) and (4), we can derive Eq. (5) as follows: M ¼ 1:4 log 10 ðPGVÞ þ 0:48 log 10 ðDÞ þ 4:54: The intensity can be calculated using PGV. In this paper, we use the Chinese Seismic Intensity Scale (GB/T 17742-2008), which can be calculated by using the relationship of PGV and intensity from the official document of China Earthquake Administration (2015) No. 18 (in Chinese): I PGV ¼ 3:00 log 10 ðPGVÞ þ 9:77; ð6Þ where I PGV is the intensity calculated from threecomponent synthesis of PGV. For simplicity, the peak ground acceleration (PGA) is not used to calculate the intensity in Eq. (6). Using Eqs. (5) and (6), we can derive the relationship of intensity I PGV with the magnitude and the distance D as follows: I PGV ¼ 2:14M À 1:03 log 10 ðDÞ þ 0:04: Thus, the intensity of any location can be calculated using Eq. (7) when the epicenter is determined by the MAXI method and the magnitude is estimated using P d .

Real-time wave propagation modeling and computation
Ground motion updating procedure solves the problem of failure calculation of the intensity distribution by GMPE. The procedure models the energy of wave propagation using radiative transfer theory and use direct simulation Monte Carlo method for efficient computation. Given a distribution of the energy on the map, we can calculate the future energy distribution of the map in real time.

Radiative transfer theory (RTT) for wave propagation modeling
Wave propagation simulation solves the problem of calculating the energy distribution of regions U n at time n based on U n-1 (the energy distribution at time n -1). In this article, RTT is used to model wave propagation. RTT interpret seismogram envelopes of high-frequency seismic waves based on the ray theory and Boltzmann equation (Sato et al. 2012). Radiative transfer equation represents the scattering and attenuation processes, and the mean energy density is used to model wave propagation. The radiated transfer equation in the 2-D space model is formulated as follows: In Eq. (8), q represents the direction of wave propagation and f (x, t; q) represents the distribution of mean energy density in direction q. dX q 0 is a solid angle element in direction q 0 and g(x) is the scattering coefficient at x with the solid angle. v 0 is the velocity of P wave (v P ) or S wave (v S ) at x. In this article, we make several assumptions: the velocity structure and the attenuation structure are homogeneous; scattering does not cause wave conversion; the isotropic scattering strength and intrinsic absorption are independent of location. Therefore, g 0 (x) and h 0 (x) are set to constants and denoted as: g 0 (x) = g 0 , h 0 (x) = h 0 . Thus, Eq. (8) can be rewritten as follows: The radiative transfer equation models the dynamic process of wave energy propagation across the 2-D space heterogeneous medium. Attenuation and scattering process from q 0 to q is modeled in the radiative transfer equation. It is possible to use Eq. (9) to simulate the process of wave propagation and to predict the distribution of wave energy with the elapse of time.

Monte Carlo method for RTT simulation
The radiative transfer Eq. (9) can be solved by efficient approximation algorithms, and we use direct simulation Monte Carlo (DSMC) method proposed by Yoshimoto (2000) to simulate the propagation of wave energy. In DSMC method, the energy can be represented as particles carrying some energy. Scattering and absorption process is represented as the movement of a large number of particles (e.g., Gusev and Abubakirov 1987;Hoshiba 1995Hoshiba , 1997Yoshimoto 2000;Hoshiba and Aoki 2015). The integral of energy in a specific region denoted by Eq. (9) is calculated as the sum of the particle energies in that region.
For representing the particle movement in 2-D space, the particles are represented as M and the mth particle at time n locates at X n,m .
In the scattering process, the DSMC method models the probability that energy does not scatter during the time interval Dt is exp[g 0 v 0 (x)Dt], which can be approximated as g 0 v 0 (x)Dt. The scattering probability density from X 1 to X 2 is 1/2p. Direction X 1 is a unit vector, and it can be represented as the travel angle h which denotes the angle between the particle direction and x-axis. In order to implement the process efficiently, two independent uniform random variables (between 0 and 1) k 1 and k 2 are used to calculate the probability of scattering and the particle's direction if it scatters. If k 1 \ g 0 v 0 (x)Dt, which means scattering occurs, the new traveling direction is X 2 = 2pk 2 , and the particle location is X n,m = X n-1,m ? v 0 Dtq X2 . Otherwise, the traveling direction of the particle remains X 1 , and the location of the particle is X n,m = X n-1,m ? v 0 Dtq X1 .
The velocity v 0 (x) in 2-D space is calculated by simulating the travel time in 3-D space. As Fig. 3 shows, the particle propagates in 3-D space, and we derive the equivalent travel velocity v 0 (x) in 2-D spaces. The particle is assumed to propagate to the target place directly from the epicenter, and the S wave velocity in 3-D space is assumed as constant v 0 . Therefore, for particle which moves from location X to X 0 , the distance of them in Because the travel time is the same, we can get D(X, X 0 )/ v 0 (x) = |D(X, R)-D(X 0 , R)|/v 0 . By taking the derivative, the S wave velocity in 2-D space can be calculated as: where h is the angle between the vertical axis and the direction from the epicenter to the location of particle X. For modeling the absorption process, the process is shown in Fig. 4. Let the mth particle has energy q n,m at time n which attenuates due to absorption as follows: ( 1 is assumed, and the energy decreases slowly in each time step. Because h 0 is assumed to be homogeneous in this article, the amount of attenuation for each particle is a constant. The energy of the particle is attenuated when time elapses. Based on the procedure above, we can estimate the future situation from the present situation of particles with location X n,m , direction X 1 and energy q n,m . The procedure repeats for m = 1, 2, 3……N, as Fig. 4 shows. The energy propagation of the wave can be modeled using the process described in Sect. 2.4. For seismic stations, the difference between the prediction of energy and the real observation of amplitude can be integrated into the modeling process. Data assimilation is used in this article to improve the prediction accuracy of wave energies, by modeling the difference of the prediction of seismic stations and the real observation of them. Data assimilation has been widely used in many areas such as predicting weather, oceanology, and atmospheric forecast (Kalnay 2003;Hoshiba and Aoki 2015). The task of the data assimilation solves the calculation of the distribution U n of wave field at t = t n given the current ground motion of stations y n at time t n , and the distribution U n-1 of wave field at time t n-1 = t n -Dt. U n a is denoted as an optimum wave field for U n ,U n b is the wave field propagating from time t n-1 . Thus, the equation can be written as follows: where A is the wave propagation process described in Eq. (9). Data assimilation in this section is expressed as In which H is the operation of getting the stations' energy by interpolating the energy of grids U b n ; H can be represented as q Â p matrix. Thus, y n À HðU n b Þ is the difference of real energy observed and energy calculated by one-step wave propagation described in Sect. 2.4. W is p Â q weight matrix, and it is used to correct the total wave energy based on the difference of wave energy distribution. It is calculated as follows: B and R represent background error covariance and observational error covariance, respectively. BH T represents the background error covariance of grids with observation points. BH T is a p Â q matrix, which can be represented as follows.
l o ij is the correlation between the ith grid and jth station, and it can be calculated as l o ij ¼ exp Àa 2 l 2 , where a is distance between the ith grid and the jth station and l is correlation distance. r denotes the standard deviation of the background errors. HBH T represents the correlation of observation points with observation points. It is a q Â q matrix, which can be represented as follows.
l jk ¼ expðÀr 2 =l 2 Þ represents the correlation of background error between the kth and jth stations. r is the distance between the kth station and the jth station, and l is the correlation distance. The observational error R for this case is: In Eq. (17), r 0 denotes the standard deviation of observational error, and I is an identity matrix, which means that observational error at different locations are assumed to have no correlations. The error ratio r 0 /r is used to determine W instead of the individual values r 0 and r. The present distribution U n a can be calculated by using all past observations y nÀ3 ; y nÀ2 ; y nÀ1 . . .; y n by using Eqs. (12) and (13) repeatedly.
In the DSMC setting, the energy is represented as particles, and the integral of wave energy with respect to direction is used instead of arrays of observations. Denote f a and f b as the energy density after and before assimilation, respectively. U n b can be calculated by the integral of the particle energies in the corresponding regions: in which U ni b is the ith element of vector U n b , F is the energy distribution in the 2-D space, x i is the location of the ith grid. U n a is combined from the present observation y n 0 using data assimilation as shown in Eq. (13). The distribution of energy density is revised based on the comparison of U a and U b . If overprediction occurs at the ith grid, then U ni a \U ni b . It is necessary to reduce the energy, and the reducing method is shown in the following equation: If underprediction occurs, which means U ni a [ U ni b . In this article, we add an extra amount of energy U ni a À U ni b À Á radiated from the location of the ith grid. The radiation is assumed to be isotropic, and f a is expressed as: The one-step calculation of data assimilation process is as follows: Given the energy density distribution f a x i ; t nÀ1 ; q ð Þf b x i ; t n ; q ð Þcan be calculated from the DSMC process described in Sect. 2.4. U n b is calculated based on Eq. (19), and U n a can be obtained in the data assimilation process described in Eq. (13). Then, f a x i ; t n ; q ð Þis revised based on Eqs. (19) or (20). By repeating this process, the energy distribution of ground motion for any region can be calculated precisely.
F(x, t) represents the distribution of energy intensity and is assumed to be proportional to the envelop of the squared seismic wave amplitude A 2 NS x; t ð Þþ A 2 EW x; t ð ÞþA 2 UD x; t ð Þ.

Ground motion prediction
After determining the present U n a , the future situation can also be estimated using wave propagation modeling described in Sect. 2.4. That is U nþ1 means the prediction at time step n ? 1. This repeating process enables us to predict the situation after k time step: Given U n a , the further energy distribution can be modeled using Eq. (12). The past U t a $ U n a have been calculated, and W n is denoted as the past maximum energy at time n. The W n can be calculated as: Denote C n as the future maximum energy to be predicted at time n, and it can be calculated as: Denote H as the initial prediction shakemap, and we get the future ground motion by combining the initial shakemap H and the current maximum energy that has happened (W) and the predicted maximum energy (C) in the near k time interval. The equation for the peak intensity U n can be formulated as follows: Using Eqs. (22-24), for any time of the earthquake, we can compute and update the real-time ground motion in real time.

Wenchuan earthquake
We apply the method described above to actually observed data in this section and use the main shock of Wenchuan earthquake as an example in this article. At 14:28 on May 12, 2008 (Beijing time), a great earthquake occurred in Sichuan Province, China, with a surface wave magnitude M S 8.0. The earthquake epicenter was located at latitude 31.021°N, longitude 103.367°E, which is close to the town of Yingxiu (31.0°N, 103.4°E), with a focal depth of 14 km. During the Wenchuan earthquake, the National Strong Motion Observation Network System (NSMONS) of China has obtained a large number of strong motion records from the mainshock, including records from 460 free-field stations in 17 provinces, municipalities, and autonomous regions (Wen et al. 2012). We use the data recorded by NSMONS, with the region 600 km 9 600 km as the study object, where the epicenter is centered on the map, as shown in Fig. 5. The seismic stations of NSMONS are sparse compared to the Kiknet in Japan, which makes the data assimilation challenging in the study. The energy radiated by the Wenchuan earthquake was so tremendous that shaking could be felt in most parts of China, and with the most violent shaking being felt in Sichuan, Gansu, and Shanxi provinces. The rupture of Wenchuan earthquake occurred over a length of 270 km along a north-northeast-striking parallel to the northeaststriking Longmenshan thrust belt, as reported by the U.S. Geological Survey, National Earthquake Information Center, in 2008. However, the source mechanism is quite complex (Zhang et al. 2008;Chen et al. 2009). As Zhang (Zhang et al. 2009) analyzed through the P wave statistics of global seismic network, the great Wenchuan earthquake occurred on a fault more than 300 km long and had a complicated rupture process of about 90 s duration time. The slip distribution was highly inhomogeneous. The whole rupture process is divided into five stages. The first stage occurred in the first 14 s after the rupture initiation, with 14% of the total scalar seismic moment releasing. The second event is the main event, and it occurred from 14 to 34 s, releasing 60% of the total scalar seismic moment. The third, fourth, last event released 8%, 17%, 6% of the total scalar seismic moment, respectively, and they occurred in the time period of 34-43 s, 43-58 s, and 58-90 s, respectively (Zhang et al. 2009). The severe ground shaking was up to 1 g and caused heavy damage to many reinforced concrete buildings, bridges, and lifelines in the regions from Yingxiu to Beichuan.
In this article, we use only the epicenter information and magnitude. For epicenter localization using the MAXI method, we use the velocity model shown in Table 1 (Hartzell et al. 2013).
For modeling the real-time updating process using data assimilation, we set the following parameter listed in Table 2. The grid is partitioned by 3 km 9 3 km grids, and the S wave velocity is assumed to be 3.464 km/s as described in Wen et al. (2012). The attenuating scattering strength and attenuating absorption strength are set as the same as in that (Hoshiba 1993). The correlation parameter l is set to 14, and the time interval is set to 1 s, which means that our method calculates the updating process every 1 s. Larger time interval may cause inaccurate prediction, but smaller time interval will cause more computation time. 1 s interval is a good trade-off for the computation.

Site amplification correction
Site amplification is an important factor to determine seismic wave amplitude in addition to source and propagation factors, and local site condition has great effects on characteristics of ground motion and earthquake damage (Wen et al. 2012). Waveforms recorded at two adjacent stations may be very different if their site amplification factors are very different. The different site classifications can induce the amplifications of response spectra in different period ranges and have been estimated in the frequency domain in many studies (Nakamura 1989;Zhao et al. 2006;Wen 2010).
Due to the lack of sufficient seismic records and the soil profile data in Sichuan province, we divide the site into four classes according to the site natural periods and compute the mean response spectral ratio for four site classes. We use horizontal-to-vertical spectral ratio (HVSR) method to derive the site amplification factor for each site. The spectral ratio method was proposed firstly by Nakamura (1989), and it is used for estimating site characteristic based on horizontal-to-vertical Fourier amplitude spectral ratio of the microtremor. By studying site characteristics of HVSR method, Yamazaki and Ansary (1997) showed that HVSR is independent of magnitude, distance, and focal depth.
We use the ratios of 5% damped velocity response spectra instead of traditional Fourier spectral ratios for the site classification, because the Fourier spectra have many spikes; for extracting the spectral peaks corresponding to the site characteristics, smoothing is inevitable. In order to remove any such spikes causing unrelated site responses, more records are needed to average the spectral ratios of those records (Zhao et al. 2006).
We use the definition of NEHRP site classes (Wen et al. 2010) based on the site natural period (s). The site are categorized into four classes, named class SC I: (rock/stiff soil) T p \ 0.2 s, SC II: (hard soil) 0.2 s B T p \ 0.4 s, SC III: (medium soil) 0.4 s B T p \ 0.6 s, SC IV: (soft soil) T p C 0.6 s. The average response spectra H/V ratio for the four classes above is calculated by 204 records of Wenchuan earthquake data, and it is shown in Fig. 6. The site-corrected waveforms can be estimated by applying fast Fourier transform (FFT), multiplying the site factor for each frequency, and then applying the inverse FFT with small time interval (e.g., 1 s or less). It is possible to convert the seismic records in real time using modern processors. Figure 7 shows the original seismic record waves of station 051LXS and the converted seismic waves after site amplification correction. The wave record in UD direction is not changed, and the spectral of wave records in EW and NS directions are divided by the response spectra H/V ratio. It is observed that the peak values of wave records in EW, NS directions are reduced.

The results and analysis of the method
We use only the velocity model of Sichuan province for prediction. First, we show the prediction results of the initial prediction calculated by GMPE, which is shown in Fig. 8. Figure 8 shows that the initial estimated shakemap prediction using the MAXI method. The MAXI method determines the epicenter located at latitude 31.354°N, longitude 102.462°E, with a focal depth 20.8 km. The epicenter estimated by MAXI is the westnorth of the real epicenter, and the focal depth is larger than the real focal depth. The shakemap is a set of concentric circles where the circle center is the epicenter computed. It is shown that the initial shakemap is underestimated by a large gap. The first reason is that when using P d method to calculate the magnitude, large earthquake is often underestimated by the P d method. The second reason is that the rupture process is too complex: in the first stage of the rupture process, only 16% of the total energy is released.
After the initial ground motion prediction, our method will update the real-time shakemap using data assimilation. Figure 9 shows the updating predicted shakemap at the time 25, 35, 45, and 55 s after the earthquake using the epicenter computed by the MAXI method. Figure 9a shows the real-time observation intensity of the seismic stations, and Fig. 9b shows the real-time ground motion computed using data assimilation. Figure 9c shows the predicted shakemap using 60 s DSMC simulation, and it can be compared with the real shakemap accumulated from the beginning of the real-time updating process to 60 s after the current time. Figure 9d shows the intensity distribution calculated by the real-time data assimilation and wave propagation modeled by RTT at the time (t ? 60) s.   Fig. 9c is expected to be closed to the ground motion shown in Fig. 9d, which calculates the peak intensity from the start to time (t ? 60) s. As Fig. 9b shows, from 35 to 55 s, the distribution of maximum intensity occurred along the Longmenshan fault, which is consistent with the analysis of (Zhang et al. 2009). The rupture propagates from the epicenter to the northeast gradually, but by RTT method we Fig. 7 Comparison of strong motion data before and after the site amplification correction. a The seismic record of the 051LXS station from Wenchuan strong ground motion data (NS, EW, UD, respectively); b the seismic record after the site amplification correction. The spectral of UD is unchanged. The spectral of NS and EW is divided by the spectral ratio of Site Class II. The inverse FFT is used to reconstruct the signal treat it as wave propagation. At the 25 s after the origin time of the earthquake, the peak ground motion does not arrive at any seismic stations, making the shakemap the same as the initial prediction. At the 35 s, the peak ground motion occurred in several stations (e.g., station code 051LXM, 051LXT), and the shakemap is updated with a great change. Considering site amplification factor, most regions could be triggered for alarm. By comparing Fig. 9c, d, we can observe that for regions far away from the epicenter, the initial prediction can produce a shakemap in advance before the real seismic wave occurs. For the region where the stations are sparse, the initial prediction can give estimation in advance for a long time. The updating process corrects the shakemap in real time. In the time from 35 to 55 s, some regions in the Longmenshan fault are underestimated and updated in the coming time. This is because the rupture process is about 90 s (Zhang et al. 2009), and new energy often occurs before the rupture process finishes. When new energy is radiated, the updating process will receive the difference of the real station energy and the predicted energy.
In order to study the shakemap difference by the error location of MAXI method, we use the real epicenter to model the wave propagation, and the ground motion prediction results are shown in Fig. 10, which also shows the updating predicted shakemap at the same time as Fig. 9, only with the difference of hypocenter location. By comparing Figs. 9 and 10, we can see that small error of hypocenter location does not affect the prediction results too much. The predicted shakemap is nearly the same as that in Fig. 9. Because of the Longmenshan fault, the energy of seismic wave propagation is not identical for different directions. The energy propagating along the Longmenshan fault can be observed from the real-time evolution of ground motion. The intensity of region in Longmenshan fault is predicted to above 5 even at the 25 s after the origin time of Wenchuan earthquake, which is enough to trigger an alarm for that region.
In order to compare the prediction time in advance, we take 4 representative stations for comparison, as shown in Fig. 11. The station code is 051LXS, 051BXZ, 051JYH, 051JZZ, respectively. The stations 051LXS, 051JYH are located in the region of Longmenshan fault, with the distance to epicenter larger gradually, 051BXZ is near the epicenter but does not located in the region of Longmenshan fault, and 051JZZ is far from the epicenter. Figure 11a shows the comparison of maximum intensity of station 051LXS predicted by our method using the MAXI epicenter and the real epicenter. Because this station is one of the earliest stations that reach peak ground motion, our updating method cannot predict its maximum ground motion in advance for the data assimilation method is based on the real-time wave propagation. Figure 11b shows the comparison of the intensity in station 051BXZ. Because the station is near the station 051BXY, when the peak ground motion occurs in 051BXY, the peak ground motion of 051BXZ can be predicted by RTT theory of wave propagation. Overprediction occurs in 051BXZ station, and the possible reason is that the site amplification factor is not estimated accurately in the stations. The station 051JYH is located in the region of Longmenshan fault, and the result is shown in Fig. 11c. The predicted intensity above 7 is nearly 10 s in advance, but it fails to predict the maximum ground motion in advance. The reason is that when the peak ground motion of 051LXS occurs, the seismic wave propagated to 051JYH where the intensity is above 7. The failure of the maximum intensity prediction is because we set the direction of energy propagation uniform. Station 051JZZ is shown in Fig. 11d. Because the station is far away from the epicenter, we can observe that our system can predict the peak intensity about 30 s before the real peak ground value occurs. When the station is not near from the epicenter, the updating scheme can predict the ground motion much earlier than the real peak value. As Fig. 11 shows, our system can predict a stable shakemap at 30 s after the origin time of Wenchuan earthquake.
In order to compare the accuracy of shakemap generation, we draw the contour with the same intensity. We compare the shakemap generated by 45, 55, 65 s prediction with the seismic intensity computed by the interpolation of For Wenchuan earthquake, our method works as follows. At 21 s after the origin time of Wenchuan earthquake, the initial prediction step detects the epicenter and the magnitude and generates the initial shakemap. After the initial prediction (starts from 22 s after the origin time of Wenchuan earthquake), the updating process updates the shakemap every 1 s. The updating process updates the shakemap until the earthquake ends.
Computation time of our method can fit the needs of real-time prediction and updating. In this example, we use p = 400 3 400 = 160,000 grids, q = 98 stations and M = 10 6 particles to simulate the energy propagation process. At each time step (4t = 1), it takes 0.6 s CPU time for data assimilation and it takes 0.38 s for prediction of the 60 s future ground motion by using a powerful desktop PC. The total time in this method is 0.98 s at each time step, and the computation time could be shortened by using more powerful computers.

Discussion
The proposed method in this article combines the merits of traditional EEW systems and real-time updating method. If the initial ground motion can be estimated precisely by traditional EEW systems, the ground motion can be determined before the destructive S-wave of the earthquake arrives. If the initial ground motion cannot be predicted by Fig. 11 Comparison of real-time intensity from seismic station record, prediction intensity using MAXI epicenter localization method, and prediction intensity using real epicenter EEW systems, the updating method models the energy of wave propagation and predicts the upcoming maximum motion that occurs in tens of seconds. This means that the proposed method can predict the upcoming peak ground motion as quickly as possible.
Compared with the traditional EEW systems, our method can predict every region and determine if an alarm is triggered for some specific regions. The updating process provides the revision of ground motion in real time, to alleviate the wrong estimation of the earthquake parameters such as the epicenter, magnitude, the GMPE equation, etc. For a new region that little earthquake occurred, our method is applicable for the region. Our method is a general method for EEW, so more information can be integrated to improve the accuracy of the ground motion calculation.
Compared with Hoshiba's numerical shake prediction method (Hoshiba and Aoki 2015), our method has several differences. Firstly, in this article our objective is to predict the peak ground motion for a specific region, while Hoshiba's method aims to predict the upcoming real-time ground motion. Therefore, our method starts predicting the peak ground motion after P wave is detected by several stations in the network. The intensity for each region is derived from the relationship of magnitude and distance between the region and epicenter. Thus, epicenter information is estimated using P wave information in our method. Secondly, our method operates on seismic network where the seismic stations are not dense enough, while Hoshiba's method needs a dense network (e.g., Kiknet) to predict the future ground motion. In the initial prediction process of our method, the epicenter location is determined before the destructive S-wave come. Therefore, incorporating epicenter information is a natural way to compensate the sparseness of seismic network. Epicenter location is incorporated into the velocity of S-wave in 2-D spaces.
When the region to be calculated is near the hypocenter, we should change the velocity of S-wave on 2-D spaces to make the calculation of wave propagation speed more precisely. Our formulation is in accordance with the real observation. The small error of hypocenter location does not affect the predicting results of ground motion very much, as illustrated in Wenchuan earthquake. The result reveals that the accuracy of the peak ground motion can be predicted even using sparse stations. Our method is practical for predicting the ground motion intensity of the district that contains less seismic stations. The strong ground motion stations in China are relatively sparse compared with the stations in Japan, USA, etc. Our method has shown potential performance for predicting the peak intensity of the region even in sparse networks. For Wenchuan earthquake whose rupture process has several substages, it is difficult to precisely estimate the magnitude using only 3 s P wave information, which leads to the inaccurate shakemap prediction for some regions in the initial prediction process. The real-time updating process is a useful tool for earthquakes whose source mechanism is complex. The real-time updating process updates the shakemap when new energy is radiated, and it makes EEW system correct the prediction of shakemap in real time.
For wave propagation modeling, the RTT theory is modeled in 2-D space. One reason is that data assimilation in 3-D space is hard, because the borehole seismic recorder cannot be distributed under the earth surface in several to tens of kilometers. Another reason is that modeling in 2-D space can fit the requirement of real-time computation. In real situation, seismic waves propagate in 3-D space. The traveling path in 3-D space from one region to another is longer than that in 2-D space. Due to the attenuation of seismic waves, the longer traveling path will lead to less energy preservation. Therefore, the energy of seismic wave travelled in 2-D space is larger than that in 3-D space, and overprediction could be observed using RTT in 2-D space.
The computation time of our method can be determined in real time when using powerful computers, including the initial prediction and updating process. Compared with the real-time delay of the transmission of network signals, the computation time fits the needs of EEW.

Conclusions
We have proposed a numerical shake prediction and updating method for earthquake early warning, which combines the merits of traditional EEW systems and realtime prediction systems. Our method predicts and updates the peak ground motions of regions with limited seismic stations in real time, after several stations detect P wave. The method is divided into two processes: initial prediction and real-time updating. The initial prediction process can provide an initial intensity estimation using P wave information, by deriving the relationship of intensity with epicenter location and magnitude. The real-time updating process updates the seismic intensity using S-wave records of seismic stations, based on data assimilation and RTT by integrating epicenter location. This method has the potential of predicting shakemap during the process of the earthquake. We apply our method to 2008 M S 8.0 Wenchuan earthquake, whose rupture process is complicated. Prediction and updating results show that the proposed method can tolerate small epicenter location error and show promising predictions in seismic networks with limited seismic stations. The peak ground motion of some regions can be predicted from several seconds to tens of seconds earlier than the real peak ground motion occurs.