STUDY ON FREQUENCY-DEPENDENT SITEAMPLIFICATION FACTORS FOR EARTHQUAKEEARLY WARNING SYSTEMXIE Quancai ∗MEE12601

Supervisor:

Toshiaki YOKOI∗∗

ABSTRACTIn the current Earthquake Early Warning (EEW) systems, the site amplification factor is usually given
as a scalar value empirically obtained. This study aims to study the method that estimates thefrequency-dependent site amplification factors for EEW: one is site effects for B-Δ method, another isto design the casual recursive infinite impulse response (IIR) filter to model these factors and tosimulate the spectra and accelerograms at target locations and also to evaluate the performance of IIR
method.
First B-Δ relation is tested using 208 earthquakes recorded at IBH10 and IBRH19 of Kik-Net.
The B-Δ relation obtained by curve fitting to the real data shows that the site condition does affect therelation. Second based on the idea of data assimilation, IIR filter method is tested using the samestrong motion data: the surface ground motion at both sites are simulated from the corresponding
borehole records using the IIR filter designed to fit the average spectral ratio of surface record toborehole’s one; and the surface motion at IBRH10 is simulated from the surface records of IBRH19
using the IIR filter designed to fit the average spectral ratio of surface record of IBRH10 to that of
IBRH19.The method using IIR filter shows a good performance in simulating acceleration waveforms
and spectra for most of the event that I used and highly improve the accuracy of the seismic intensity
estimation.
Keywords: B-Δ, Frequency dependent Site Amplification Factors, Casual Recursive IIR Filter.

1. INTRODUCTION
There are many papers which describe the science, engineering, and societal considerations of theactive warning systems in the world. In the current Earthquake Early Warning (EEW) systems, the site
amplification factor is usually given as a scalar value empirically obtained. I studied the frequency
dependent site amplification factors in this paper. I got the data from the Kik-net (NIED) recorded atsite IBRH10 and IBRH19 for 208 earthquakes. B-Δ relation is tested using 208 earthquakes recorded
at IBH10 and IBRH19 of Kik-Net. I get the simulated results by using the IIR method for the samedata.

2. RESULTS AND DISCUSSION
2.1. B-Δmethod∗

Institute of Engineering Mechanics (IEM), China Earthquake Administration (CEA), Harbin, China.

∗∗

Building Research Institute (BRI), Tsukuba, Japan.

1

By fitting Bt*exp (-At) function to the log-transformed acceleration waveform (UD component)
envelope of the first two seconds or three seconds from the P onset, we extracted the envelope feature
characteristics by 'A' and 'B' value. (Odaka et al., 2003; Kamigaichi, et al., 2004). I calculated the B
value for the vertical component in the borehole (UD1) and on the ground surface (UD2) at IBRH10
and IBRH19. All the calculated results are shown in Figure 1 to Figure4. I calculated the B value fromthe P-wave onset time 2 seconds (Figure 1 and Figure 2) and 3 seconds (Figure 3 and Figure 4). The
black curves in Figure 1 through Figure 4 show the currently used JMA formula:
log10B=-2.0008log10Δ+3.9458 where Δ is in km, B in gal/s (Noda et al., 2011). This formula is applied
to all stations of the JMA EEW system, uniformly. The red curves in Figure 1,Figure 2,Figure
3,Figure 4 are those fitted to the data obtained from the borehole records, whereas the blue ones fromthe surface records. Comparison Figure 1 with Figure 2 and also comparison Figure 3 with Figure 4
show that there are not any significant differences among the fitted curves of the B values: one
calculated using the time window 2 seconds and another 3 seconds. This is consistent with Odaka et al.
(2003). In Figure 3 and Figure 4, the JMA B-Δ formula (black curves) are almost parallel to thoseestimated from the borehole and surface data in this station (red and blue ones, respectively) and are
drawn above them. Especially to blue ones (surface records) the black curve shows a similarity.
10

IBRH10-2

2

Borehole

Surface

IBRH10-3

2

10
Current JMA

BoreholeFitting

Borehole

Surfacefitting

log10(B)=-2.008log10(Delta)+3.9458

10

Current JMA

BoreholeFitting

Surfacefitting

log10(B)=-2.008log10(Delta)+3.9458
1

1

10

0

B v a lu e

B v a lu e

10

Surface

-1

10

0

10

-1

10

10

log10(B)=-0.6362log10(Delta)+1.5611

log10(B)=-0.9086log10(Delta)+1.315

log10(B)=-0.5674log10(Delta)+1.4651
40

50

60

70

80

90

100

110

log10(B)=-0.9158log10(Delta)+1.2985

-2

10

-2

120

40

130 140

50

60

70

Figure 1. The B-Δ relationship at IBRH10
station: using 2 seconds after P wave arrival

Borehole

Surface

Current JMA

BoreholeFitting

10

Surfacefitting

log10(B)=-2.008log10(Delta)+3.9458

1

10

0

10

-1

10

10

log10(B)=-2.0638log10(Delta)+3.8407

110

120

130 140

IBRH19-3

2Borehole

Surface

Current JMA

BoreholeFitting

Surfacefitting

log10(B)=-2.008log10(Delta)+3.9458

1

50

60

0

-1

log10(B)=-1.9841log10(Delta)+3.6901

log10(B)=-1.7834log10(Delta)+2.8652

-2

40

100

B v a lu e

B v a lu e

10

10

10

90

Figure 2. The B-Δ relationship at IBRH10
station: using 3 seconds after P wave arrival.

IBRH19-2

2

10

80

Distance(Km)Distance(Km)

70

80

90

100

110

120

130

10

140

Distance(Km)log10(B)=-1.5951log10(Delta)+2.4807

-2

40

50

60

70

80

90

100

110

120

130

Distance(Km)Figure 3. The B-Δ relationship at IBRH19
station: using 2 seconds after P wave arrival

Figure 4. The B-Δ relationship at IBRH19
station: using 3 seconds after P wave arrival.

2

140

The formula fitted to the data of IBRH19 on the surface (log10B=-2.0638log10Δ+3.8407 or
log10B=-1.9841log10Δ+3.6901) shows much similarity with the current JMA formula(log10B=-2.0008log10Δ+3.9458). On other hands, in Figure 1 and Figure 2, this parallelism is lost. The
black curves tend to stick to the blue one at short distance and to the red one at the long distance. The
formula fitted by the data recorded by borehole sensor in IBRH10 (log10B=-0.9068log10Δ+1.3150 or
log10B=-0.9158log10Δ+1.2985) shows much difference from the current JMA formula. Even in thesame station, the blue curves calculated using the data recorded in the surface sensor are always drawn
above the red curves calculated using the data recorded in the borehole sensor. This shows the site
amplification factor did affect the B-Δ relationship as well as Nakamura et al. (2006) and JMA (2011)
reported. It means that the parameter estimated in EEW systems is different when using the dataobtained from different station or site condition even for the same earthquake. So for IBRH10, thisdifference cannot be ignored. We may use different B-Δ relationship according to different site
conditions for more accurate calculation. More accurately, we may use the empirical B-Δrelation ateach station for the correction of the B-Δ relationship for the EEW system. We may use the method
proposed by Hoshiba (2013a and b) to skip the source estimate procedure as well in the future.
2.2. Simulation Between Two Stations
Next, we simulated waveforms at station IBRH10 using strong motion data recorded by the surface
sensor at IBRH19. I show two examples of the 208 simulation results. The red curves are simulated
results while the blues ones are observed results in Figure 5 through Figure 10. Figure 5 and Figure 6
show the simulation results for the M5.2 earthquake which occurred on Feb. 19, 2012. The
amplification factors are about 4 (17.7→69.7 and 13.6→53.2) for horizontal component and about 4(11.2→38.9) for vertical component. The increment of seismic intensity is 1.5(2.2→3.7). The seismic
intensity difference between the observed surface acceleration and the simulated surface accelerationfrom the observed borehole acceleration is just 0.2 (3.5-3.7). When compared the observed
acceleration and spectrum of the surface record with the simulated acceleration and spectrum, we caneasily draw the conclusion that they are approximately simulated well. Figure 7 and Figure 8 show thesimulation results for the M5.1 earthquake which occurred on April 14, 2011. The amplification
factors are about 2 (7.0→16 and 7.0→17) for horizontal component and about 4 (4.2→16) for verticalcomponent. The increment of seismic intensity is 1.4 (1.3→2.7). The seismic intensity residuals
between the observed surface acceleration and the simulated surface acceleration from the observed
borehole acceleration are 0.0 (2.7→2.7). When compared the observed acceleration and spectrum of
the surface record with the simulated acceleration and spectrum, we can easily draw the conclusionthat they are approximately simulated well. Compared these results during two earthquakes, thedifferent amplifications of maximum acceleration reflect the differences of the frequent contents of theincident waveforms that cannot be reproduced by a scalar site amplification factor.
Although most of the simulation results show good performance of this method, some simulation
results did not provide reasonable estimations. For example, Figure 9 and Figure 10 show thesimulation results for the M5.3 earthquake which occurred on March 16, 2011. The amplification
factors are about 6(6→34 and 5→36) for horizontal component and about 6 (3.8→22 ) for verticalcomponent. The increment of seismic intensity is 2.3 (1.5→3.8). The seismic intensity difference
between the observed surface acceleration and the simulated surface acceleration from the observed
borehole acceleration is 1.1 (2.7→3.8). Comparison between observed and simulated accelerograms
and spectra indicate that this simulation did not provide reasonable estimation. The observed
amplification factors are about 2 (6→10 and 5→12) for horizontal component and about 1.5 (3.9→6.2) for vertical component.

3

20

40

100

60
Time(s)80

100

Amplitude(cm/s)
120

Amplitude(cm/s)-100
0

SMFromIBRH19EW2 PGA=53.2097 gal
IBRH10EW2 PGA=45.7466 gal

0
-100
0

20

40

50

60
Time(s)

80

100

120

Amplitude(cm/s)

Acceleration(gal)SMFromIBRH19NS2 PGA=69.7177 gal
IBRH10NS2 PGA=50.9617 gal

Ijma=3.7
Ijma=3.50Acceleration(gal)Acceleration(gal)

IBRH191202191454.kwin.surfsimulated10--Acceleration
100

SMFromIBRH19UD2 PGA=38.9243 gal
IBRH10UD2 PGA=23.7674 gal

0
-50
0

20

40

60
Time(s)

80

100

120

Ijma=2.7
Ijma=2.7

Amplitude(cm/s)

IBRH191104140735.kwin.surfsimulated10--Acceleration
50

SMFromIBRH19NS2 PGA=16.6861 gal
IBRH10NS2 PGA=25.4155 gal

0

20

40

20

60
Time(s)80

100

SMFromIBRH19EW2 PGA=17.6228 gal
IBRH10EW2 PGA=18.1371 gal

20

40

20

60
Time(s)

80

100

40

60
Time(s)

80

100

50

60
Time(s)

80

100

120

Amplitude(cm/s)

40

120

120

SMFromIBRH19EW2 PGA=36.1425 gal
IBRH10EW2 PGA=12.569 gal

0
-50
0

20

40

50

80

100

SMFromIBRH19UD2 PGA=22.1543 gal
IBRH10UD2 PGA=6.2254 gal

0
-50
0

60
Time(s)20

40

60
Time(s)80

100

SMFromIBRH19EW2

IBRH10EW201

10Frequency(Hz)10

SMFromIBRH19UD2

IBRH10UD2

2

0

-2

10 -1
10

0

10
Frequency(Hz)

10

1

2

201104140735--AmplitudeSpectrum

10

0

10

-2

10 -1
10

SMFromIBRH19NS2

IBRH10NS2

0

1

10

10

2

10

0

10

-2

10 -1
10

SMFromIBRH19EW2

IBRH10EW2

0

1

10

10

2

10

0

10

-2

10 -1
10

SMFromIBRH19UD2

IBRH10UD2

0

1

10

10

Figure 8. Comparison between the simulated
and observed spectra for the earthquake
1104140735

Amplitude(cm/s)

IBRH191103162239.kwin.surfsimulated10--Acceleration
Ijma=3.8
SMFromIBRH19NS2 PGA=34.8717 gal
Ijma=2.7
IBRH10NS2 PGA=10.7058 gal

20

1Frequency(Hz)Amplitude(cm/s)

Acceleration(gal)Acceleration(gal)Acceleration(gal)

0
-50
0

10

10

0

-2

10

IBRH10NS2010Frequency(Hz)

2

10 -1
10

120

Figure 7. Comparison between the simulated
and observed accelerograms for the earthquake
1104140735
50

10

SMFromIBRH19NS2

Frequency(Hz)SMFromIBRH19UD2 PGA=16.7134 gal
IBRH10UD2 PGA=14.6136 gal

20

10

120

0
-20
0

-2

Frequency(Hz)

0
-20
0

0

10 -1
10

120
Amplitude(cm/s)

-50
0

10

201202191454--AmplitudeSpectrum

2

Figure 6. Comparison between the simulated
and observed spectra for the earthquake
1202191454

Amplitude(cm/s)

Acceleration(gal)Acceleration(gal)Acceleration(gal)Figure 5. Comparison between the simulated
and observed accelerograms for the earthquake
1202191454

10

Figure 9. Comparison between the simulated
and observed accelerograms for the earthquake
1103162239

4

10
10

2

0

-2

10 -1
10
10
10

0

10
Frequency(Hz)IBRH10NS2

10

1

0

-2

10

SMFromIBRH19NS2

2

10 -1
10
10

201103162239--AmplitudeSpectrum

SMFromIBRH19EW2010Frequency(Hz)IBRH10EW2

10

1

2

0

-2

10 -1
10

SMFromIBRH19UD2010Frequency(Hz)IBRH10UD2

101

Figure 10. Comparison between thesimulated and observed spectra for theearthquake 1103162239

I selected 40 earthquakes randomly from 208 earthquakes to calculate the observed seismic
intensity and simulated one. The seismic intensity residuals are shown in Figure 11. The average
seismic intensity difference of these 40 earthquakes is 0.31. The standard deviation of the residuals is
0.23. 67.5% of the difference is less than 0.3. 82.5% of the residuals are less than 0.5. 97.5% of thedifference is less than 1. According to the paper (Iwakiri et al. 2011), JMA used ARV method based
on topographic data, they investigated the station correction method based on empirical site
amplifications obtained from recent observed seismic intensity data. Through the simple replacementof ARV by empirically estimated station corrections, they conclude that the number of residuals
within+-1.0 was increased from 84% to 93%, and the number of the residuals within +-0.5 wasincreased from 55% to 59%. The comparison results between different methods are shown Table 1.
Table 1. The comparison result
Methodmean residual standard deviation +-0.5 +-1.0
ARV(Iwakiri et al. 2011),
0.25
0.63 55% 84%
Station Correction (Iwakiri et al. 2011),
0.19
0.55 59% 93%
Casual IIR Filter
0.31
0.23 82.5% 97.5%
Comparing with these results, we can conclude that the method used in this paper could improve theaccuracy of the seismic intensity estimation and have good performance.

Figture 11. The residuals between the simulated seismic intensity and observed seismic intensity.
I evaluated the performance of this method using Kik-net data, the advantage of IIR method asshown below:
Comparison between the observed and simulated accelerograms and spectra indicate our method
shows good performance for simulating acceleration and spectrum for most of the data used in thisanalysis. Compared the simulated results of different earthquakes, the different amplifications ofmaximum acceleration reflect the differences of the frequency content of the incident waveforms thatcannot be reproduced by a scalar site amplification factor. This method highly improves the accuracyof the seismic intensity estimations in comparison with the current EEW systems. As the records ofthe front stations includes the real source factor and real wave propagation factor. The effects ofrupture directivity, source extent and simultaneous multiple events are included in this method. An
error of the source parameter estimate does not affect to the accuracy of the ground motion estimatesby this method. It would have more excellent performance than current EEW systems, especially when
considering some earthquakes which occurred simultaneously.
However, as shown in Figure 9 and Figure 10, the simulations did not provide reasonableestimations. The differences of site amplification factors are large. It may mean that the spectral ratio
for this earthquake is not suitable. There exists the possibility that the IIR filter parameter may be5sensitive to some effect not considered in this method. As the accuracy of this method depends on thefilter design, we can characterize the IIR filter parameter by the iteration procedure more accurately
mentioned in detail in Section 2.2. This problem may be overcome by more accurate design of thefilter for these data. This method focused on the amplitude characteristics of site amplification factor
and ignored the phase characteristics. When this method is applied to low frequency waves, phasecharacteristics may not be ignored. There had been very few researches on how to evaluate the phase
information of site amplification factors. Although there were some researches, operation was
performed in the frequency domain and then transformed to time domain, which is not suitable for
EEW purpose. This difficulty may be improved by seismic interferometry using earthquake records
triggered by small and intermediate size earthquake (e.g. Yamada et al., 2011) and more researches areneeded to improve the accuracy of EEW systems.

3. CONCLUSIONS
The fitted B-Δ formulas for site IBRH10 and IBRH19 (Kik-net) were obtained using strong motion
data recorded by the borehole and surface sensor. By the calculation, 12 spectral ratios were obtained.
Using the casual recursive filter designed to fit the averaged spectral ratio and borehole records at
IBRH10, I got 208 simulated surface accelerograms and spectra for IBRH10. Also, I got 208
simulated surface accelerograms and spectra for IBRH10 from surface records of IBRH19.
Our results suggest that the site condition does affect the B-Δ relationship. I suggested somesolutions that can improve the current B-Δ Method. We can use different B-Δ relationship according
to different site conditions for more accurate calculation. More accurately, we may use the empirical
B-Δ relation at each station for the correction of the B-Delta relationship for the EEW system. We may
also use the IIR method proposed by Hoshiba (2013a and b) to skip the source parameter estimate
procedure. IIR method has good performance in simulating acceleration and spectra for most of thedata used in my analysis and highly improve the accuracy of the seismic intensity estimation whencompared with the current EEW systems. The accuracy of IIR Filter method depends on the accuracyof the input spectral ratio and the design of the IIR filter. A possibility for accuracy improvement byusing larger number of first order and second order filter, however, the design procedure of this filter
will be time consuming. I have made an action plan to do some research and education about the EEWin China.
ACKNOWLEDGEMENTS
I would like to express my heartiest gratitude to Dr. Mitsuyuki Hoshiba (MRI/JMA) for supervisingme and giving me very useful advices and suggestions during my individual study. Thanks to NIED
for providing the strong motion data. The JMA Seismic Intensity is calculated by ViewWave
developed by Dr. Toshihide Kashima.

REFERENCES
Hoshiba, M., 2013a, J. Geophys. Res.: Solid Earth, 118, 1–13.
Hoshiba, M., 2013b, in press.
Iwakiri,K., et al.., 2011, Earth, Planets and Space.63, 57-69.
Japan Meteorological Agency, 2011, JMA Report, in Japanese.
Kamigaichi, O., et al., 2009, Seism. Res. Lett., 80, 717–726.
Nakamura, H., et al., 2009, Geophys. Res. Lett., 36.
Noda, S., et al., 2011, RTRI REPORT Vol. 25, No. 7, Jul. (in Japanese with English Abstract).
Odaka, T., et al., 2003, Bill. Seism. Soc. Am., 93, 526–532.
Yamada,M., et al. 2010, J. Geophys. Res., 115, B03302.

6