Study on Site Amplification Factor Correction for Earthquake Early Warning System

Study on Site Amplification Factor Correction for Earthquake Early Warning System Quancai Xie, Qiang Ma, Jingfa Zhang, Haiying Yu Key Laboratory of Earthquake Engineering and Engineering Vibration, Institute of Engineering Mechanics, China Earthquake Administration, Harbin, 150080, China 5 Institute of Crustal Dynamics, China Earthquake Administration, Beijing 100085, China Correspondence to: Quancai Xie (xiequancai@iem.ac.cn)


Introduction
In recent decades, real-time strong ground motion prediction has become an important part of earthquake early warning systems.In the world, there are many countries deployed operational earthquake early warning systems like Mexico (Espinosa-Aranda et al. 1995;Espinosa-Aranda et al. 2009), Japan (Kamigaichi et al. 2004;Hoshiba 2008;Nakamura et al. 2009), Taiwan (Wu et al. 2002;Hsiao et al. 2009), Turkey (Erdik et al. 2003;Alcik et al. 2009), and Romania (Wenzel et al. 1999;Ionescu et al. 2007) Also there are many some earthquake early warning system under developing and testing like the Unites States (Allen et al .2003;Allen et al. 2009;Bose et al. 2009;), Italy (Zollo et al. 2006;Zollo et al. 2009), and China (Peng et al. 2011).
Usually, the earthquake early warnings systems can be categorized as offsite (also known as regional EEW) and onsite warnings.The offsite warning utilizes a few seconds of seismogram observed at the first station, and then estimate the source parameter, such as magnitude, epicentre distance or others.Then according to the parameter estimated, a warning can be manually or automatically issued based on some rules made before an earthquake occurs.Hoshiba et al., 2008 mentioned that the Japan earthquake early warning system has been operational nationwide since October, 2007 by Japan Meteorological Agency.Approximately 1,100 stations from JMA network and the high sensitivity seismograph network (Hinet) were used to determine the hypocenter of Japan Meteorological Agency earthquake early warning system.In order to disseminate the warning quickly, hypocentre estimation should be done just after the first detection of the P phase at a single station.In order to ensure the reliability of the estimation, the B-Delta Method (Odaka et al., 2003) and Network Method (Horiuchi et al., 2005;Horiuchi et al., 2009) are used in combination.Usually, the current Japan Meteorological Agency earthquake early warning system works well.But after the main shock of the 2011 great Tohoku Earthquake, the earthquake early warning system did not work well due to high aftershocks activity and high background noise, as well as power failure and wiring disconnections (Hoshiba, 2011).Earthquakes that occurred nearby simultaneously in different locations also made the system provide false information.The site amplification factor was usually considered as scalar values, such as amplification of peak ground acceleration or peak ground velocity, increments of seismic intensity, in the conventional earthquake early warning system.There are some research papers on how to improve the site amplification factor for considering more accuracy calculating Japan Meteorological Agency seismic intensity of Japan earthquake early warning system (e.g., Iwakiri et al., 2011).Among them we chose the new idea proposed by Hoshiba (2013) and used it to design casual Filter for modelling the site amplification factor of Kik-net stations .We focus on full evaluate the performance of this method by selecting large amount of Kik-net data triggered in more than one thousand earthquakes, then we make simulation from borehole to surface and also simulation from front-detection station to far-field station.Then compare the statistical simulation result with other methods considering the accuracy of the seismic intensity prediction and clarify the advantages and some problem need to be considered when utilizing it to the earthquake early warning system.

Data
The hypocentre parameters including origin time, location of hypocentre, and magnitude were obtained from the JMA seismic catalogue.The strong motion data were downloaded from the website (http://www.kyoshin.bosai.go.jp/).The advantage of this network is that all stations have a borehole of 100 m or more in depth, with accelerographs installed both on the ground surface and at the bottom of boreholes.The site information measured in the boreholes includes soil type along with P and S wave velocity profiles.The sampling frequency is 200Hz for the records before November 2007 and is 100Hz thereafter.In this analysis, we use records observed at 2 stations.One of them with site code IBRH10 has been in operation since Septemper.1,2000 and the other with site code IBRH19 since May 15, 2004, respectively. Until Decmber.31, 2012, IBRH10 and IBRH19 recorded 1119 and 910 events respectively.I selected 673 strong ground motion recorded at the surface and borehole sensor when both of these two stations were triggered by an earthquake.The inner distance between IBRH10 and IBRH19 is 14.6 km.I selected data with hypocentre distance larger than at least three times of the inner distance.The number of earthquakes is up to 553 and the range of magnitude was between 3.3 and 9.The recording time spans from May.16, 2004to Dember.31, 2012.Excluding the earthquake which occurred in the sea area, the number of earthquakes ranging from Magnitude 3.3 to Magnitude 7 adds up to 208 (Figure 1).There exists 20 meter soft sediment at station IBRH10.The layer of appears at the depth of 518m.The IBRH19 is almost a completely rock site station.The site profiles can be downloaded from Kik-net website.

Theory and Methodology
Source parameters including hypocentre location and magnitude are determined within a few seconds after an earthquake occurrence.Then the ground motions are estimated based on these parameters.While the EEW using a few parameters, parameter uncertainty leads to another error in the ground motion prediction.A new method was proposed by Hoshiba (2013).The method predicts ground motion using ground motions observed at front stations in the direction of incoming waves.The idea of this method were shown in Figure 2. In this method, the observed information is sent directly forward to the target point.The core idea of this method is that the frequency-dependent site amplification factor can be reproduced by a casual recursive filter based on historical relative spectrum ratio between two stations.In the method the far-field simulated waveform can be obtained by real-time filtering of the observed waveform recorded in the front detection station.
Seismic ground motions are often modelled by convolution of source, propagation and site amplification factors.Site amplification factor plays important role in determine seismic wave amplitude except for propagation effect and source effect.Usually, the site amplification factor was evaluated in frequency domain.However, for earthquake early warning systems it is not suitable as this procedure needs some length windowed waveform for FFT in frequency domain.In many previous studies, site amplification factors are estimated using following equation.(e.g., Iwata and Irikura (1988)) Where f is frequency in HZ,   (),   (),   (), and  ()represent the observed seismic spectrum from event k at site l, the source spectrum characterizing the event k, the site amplification factor at site l, and the propagation factor between event k and site l respectively, and f is the frequency of the seismic waves.The frequency-dependent relative site amplification factors are assumed to be modelled by a following linear system of firstand second -order filters, ) .
+ 1 + 2  =1 . ∏ ( (2) where N and M stand for the numbers of the first and second-order filters, respectively, and s = iω.Here ω are angular frequencies and h are damping factors that characterize the frequency dependence, respectively. 2 + 2ℎ +   2 represents a damping oscillation (Hoshiba, 2013).0 ,ω 1 , ω 2 , ω 1 , ω 2 are estimated for given values of N and M by using the leastsquares method in logarithmic scales.We focus on the amplitude characteristics, ignoring phase characteristics.The bilinear transform (also known as Tustin's method) is introduced as Which is used in digital signal processing and discrete-time control theory to transform continuous-time system representations to discrete-time.Then the pre-warping equation is applied to ω 1 , ω 2 , ω 1 , ω 2 , Then the transfer functionF(z) is obtained, whereΔT is the sampling interval of the digital waveforms and z = exp (sΔT).Eq. (3) and Eq. ( 4) are the necessary procedures to obtain the coefficients of a causal recursive filter for real time processing.

Spectral Ratios
We use the strong motion data recorded by IBRH10 and IBRH19 during these 208 earthquakes.The spectral ratio results obtained are shown in Figure 3

Simulation from borehole to surface
Firstly, we make the simulation from borehole to surface, although it is not useful for earthquake early warning system.But it could be used to make full evaluation of this method.We use the strong motion data recorded by IBR10 borehole sensor to simulate the surface station acceleration waveforms and spectrum.For these two examples, compared the simulated acceleration and spectrum with the observed acceleration and spectrum of the surface records, the result simulated well.The different amplifications of maximum acceleration between Table 2 and Table 3 reflect the differences of the frequency contents of the incident waveforms that cannot be reproduced by a scalar site amplification factor (e.g., amplification of peak ground acceleration or peak ground velocity, or increment of seismic intensity).
The seismic intensity is calculated according to the method described in the paper (Yamazaki et al. 1998).Figure 8 shows the seismic intensity residuals.The average seismic intensity residuals of these 208 earthquakes is 0.139.The standard deviation of the difference is 0.254.98.6% of the seismic intensity residuals is less than 0.5.100% of the seismic intensity residuals is less than 1.As the minimum resolution for the seismic calculation is 0.1 degree, It is reasonable that we consider these simulations have good performance.

Simulation Between Front-detection Station and Far-flied Station
Then, Using the surface strong motion data of IBRH19,we get simulated waveforms for IBRH10.3. The simulation results for the M5.1 earthquake occurred on April 14,2011 is shown in Figure 10 summarized in Table 4. Compared the observed the surface record acceleration and spectrum with the simulated acceleration and spectrum, it shows that the simulation result is well.The acceleration amplification factor for the he M5.2 earthquake which occurred on Feb. 19, 2012 is 3.9,3.9,3.5 respectively, while The acceleration amplification factor for the he M5.1 earthquake which occurred on the M5.1 earthquake occurred on April 14,2011is 2.4,2.5,3.9 respectively.Comparing the simulation results for these two earthquakes, the different amplifications of acceleration between Table 3 and Table 4 shows the different frequent contents of the waveforms that cannot be reproduced by a scalar site amplification method.
Although most of the simulation result shows good performance, there exists the case that the simulation did not work well.
For example, Figure 11 The information about the site amplification factors and the increment of seismic intensity are summarized in Table 5.Compared the observed acceleration and spectra with simulated acceleration and spectra, it indicate that this simulation did not work well for the M5.3 earthquake occurred on March 16, 2011.The seismic intensity is calculated according to the method described in the paper (Yamazaki et al. 1998).Figure 12 shows the proposed the station correction method.The station correction method based on site amplifications obtained empirically from observed seismic intensity data.They compared the performance of the ARV method and station correction method.For the ARV method and station correction method, The seismic intensity residual within +-0.5 is 55% and 59% respectively.
The seismic intensity residual within +-1 is 84% and 93% respectively for ARV method and station correction method.The comparison result between different methods are shown in Table 6.The statistical 1 degree seismic intensity error of the current JMA EEW system (JMA,2018) was shown in figure 13.The average 1 degree seismic intensity error is 74.74% for all the eleven years data.The best result is 93.7% in 2017,the best result is 34.6% in 2010.From the analysis mentioned above, we can conclude that this method could improve the accuracy of the seismic intensity estimation.It will have good potential application in the future EEW system.

Figure 1 .
Figure 1.The station and epicentre distribution map used for this research.

Figure 2 .
Figure 2. Comparison of the method proposed by Hoshiba with the network method.(Modified after Hoshiba, 2013)

Figure 7
Figure 7(a) to Figure 7(d) show the simulation results for the M4.6 earthquake which occurred on December 7, 2012.The

Figture 6 .Figture 7 .Figture 8 .
Figture 6.An example of the results of simulation for IBRH100911211539: (a) the observed accerlation at the borehole, (b) the observed surface acceleration waveform (blue) compared with the simulated one (red), (c) the spectra of the observed surface record, (d) the spectra of the observed surface record (blue) compared with the simulated one (red).
Figure 9(a) through Figure 9(d) show the simulation results for the M5.2 earthquake which occurred on Feb. 19, 2012.The information about site amplification factors and the increment of seismic intensity are summarized in Table (a) to Figure 10(d) .The comparison information about the site amplification factors and the increment of seismic intensity are Nat.Hazards Earth Syst.Sci.Discuss., https://doi.org/10.5194/nhess-2018-400Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 4 February 2019 c Author(s) 2019.CC BY 4.0 License.
(a) to Figure 11(d) show the simulation results for the M5.3 earthquake occurred on March 16, 2011.
Nat. Hazards Earth Syst.Sci.Discuss., https://doi.org/10.5194/nhess-2018-400Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 4 February 2019 c Author(s) 2019.CC BY 4.0 License.seismic intensity residual.The average seismic intensity residual is 0.35 for all our data set.The standard deviation of seismic intensity residual is 0.36.69.7% of the seismic intensity residual is less than 0.5.98.1% of the seismic intensity residual is less than 1.Japan Meteorological Agency used ARV method (amplitude ratio of peak ground velocity at the ground surface relative to the engineering bedrock of averaged S-wave velocity 700 m/s) based on topographic data.Iwakiri et al. 2011

Figure 12 .
Figure 12.The seismic intensity residuals between the observed data and simulated data.

Figure 13 .
Figure 13.The percent ratio diagram for 1 degree seismic intensity error in the current Japan EEW System