ESTIMATION OF HIGH FREQUENCY ENERGYRADIATION (HFER) OF THE 2003 TOKACHI-OKI
EARTHQUAKE USING EMPIRICAL GREEN’S
FUNCTIONS
Wang Dun
MEE08155

Supervisor: Tatsuhiko HARAABSTRACTWe determined the centroids and stopping point of the high frequency energy radiation (HFER) of the2003 Tokachi-Oki earthquake (Mw 8.1). We used tele-seismic data and basically followed theprocedure of Gusev et al. (2007), who used the empirical Green’s functions (EGFs) technique. We firstchose several small events that occurred in the nearby region of the source. The waveform data wereretrieved from the Data Management Center of the IRIS. We calculated correlation among theirsmoothed envelopes in the frequency bands, 0.4-1.2 Hz, 1.2-2.0 Hz, 2.0-3.0 Hz, and 3.0-4.0 Hz toobtain suitable EGFs. We calculated envelopes of the P-wave waveforms of the main shock and EGFs
in the above four frequency bands, and smoothed them. Then, we performed deconvolution of thesmoothed envelopes of the main shock using those of EGFs to obtain the HFER signals in each bandfor each station. A non negative constraint was applied to this calculation. Finally, we calculated thetimes and locations of the centroids and the stopping point of the HFER.
The obtained centroids of the HFER are consistent with those obtained by previous studies
using different data sets and methods. The centroid times and source time are well constrained, whileconstraint for the centoid locations and stopping point are relatively weak.
Our results show that this method is applicable to determination of centroids of HFER oflarge earthquakes for which strong motion data is not available. One of the advantages of this method
is that it is possible to obtain good azimuth coverage using tele-seismic waveform data.
Keywords: High frequency energy radiation, empirical Green’s functions, Tokachi-Oki earthquake.
1. INTRODUCTION
China is one of the earthquake-prone countries, and has many damaging earthquakes in its history.
One of the recent devastating earthquakes is the May 12, 2008 Wenchuan earthquake (Mw 8.1).
Nearly 70,000 people have died, thousands are missing, and more than 1.5 million people lost theirhomes by this earthquake. Therefore, studies on earthquakes are very important in China.
Among seismological data analysis on earthquake sources, earthquake rupture studies canprovide useful information to improve understanding of earthquake physics. Also, they provide useful
information for strong ground motion simulation, which is effective for earthquake disaster mitigation.
In the present study, in order to investigate the character of its high frequency energy radiation, weanalyzed high frequency signals (0.4-4 Hz) of tele-seismic broadband seismograms of the September
25, 2003 Tokachi-Oki earthquake (Mw 8.1); we determined the centroids and stopping point of thehigh frequency energy radiation.
Institute of Seismology, China Earthquake Administration, China.
International Institute of Seismology and Earthquake Engineering, Building Research Institute, Japan.

12. DATA
We selected several events that occurred in the source area as candidates of EGFs in this study (Table1). In this research, we need to choose EGFs with proper magnitude and good S/N ratio in order toobtain good accuracy in the following inversion. For this purpose, we checked the data quality visuallyand numerically after splitting them into the four bands and envelopes calculating.
Broadband vertical velocity
Table 1 The candidate events for the EGFs of the 2003
seismograms(P waves only) of our target
Tokachi-Oki earthquake (after NEIC).
earthquake (25 September, 2003,
Time
Latitude Longitude Depth (km) Mw
hypocenter=41.78 ° N, 143.86 ° E,
2003/10/09/06:53
41.56
143.8731.0
5.6
depth=27km, and magnitude=8.1 after
2007/02/17/00:03
41.65
143.97
38.2
6.0
USGS) and the four earthquakes, are2009/06/05/03:30
41.82
143.45296.4retrieved from the Data Management
2008/09/11/00:21
41.89
143.7525.0
6.8
Center of the IRIS; their epicentral
distances are constrained between 30°and 85°.
3. THEORY AND METHODOLOGY
3.1. Basic theory
The EGFs method is one of the most powerful methods for predicting strong ground motions of large
earthquakes (e.g., Irikura, 1986). It is also applied to determining HFER of huge earthquakes (Gusevet al., 2007). It uses observed waveforms of small events as Green’s functions to consider the effectsof seismic wave propagation. Suitable events that can be EGFs of certain earthquake are not always
available, which consequently put limitations on its application. In this study, there are some suitable
aftershocks and small events that occurred close to the source region.
After records for each station for the main shock and the small event chosen as EGFs aresmoothed properly, deconvolution is performed using the following convolution equation:

w1 (t )  wa (t )  wm (t ) ,(1)where wm (t ) and wa (t ) are squared envelopes of the main shock and the EGFs time series at eachstation, and w1 (t ) is the source signal of the main shock to be determined. By deconvolution for theequation (1), we can obtain w1 (t ) . Since source radiations are always positive, we performed nonnegative inversion following Okamoto and Takenaka (2009).
We solved the damped least squares to avoid the numerical instabilities using(GT G+ 2I)m＝GT d ,(2)where G is the kernel matrix consisting of the Green’s functions, I is the identity matrix,  2 is thedamping factor, m is the model vector, and d is the data vector of the main shock (Hara and Yokoi,2008).
Because of some noisy records, results obtained from deconvolution need to be checked
carefully, and some of them should be discarded manually by visual check. Once we obtain the source
radiations at each station, the following step is determining the centroid time and the signal duration.
For the latter, we need to set a threshold for all the signal time series, and obtain the longest duration
for all bands in each station. We set suitable threshold by trial and error.

2

3.2. Relative location method
Following Gusev et al. (2007), we use the relative location method to locate the centroids and thestopping point of the HFER. The following equations are used in this procedure:M t－(r1k M 1 +r2k M 2 )/c=e kfin

M t －(r1k M 1fin +r2k M 2fin )/c=Tfin ,

(3)
(4)

where M t is the average retardation of radiation with respect to the origin time, and M 1 and M 2 donate
relative location of centroid of HFER in the north south and eastwest directions, respectively. rik is theunit ray vector from the hypocenter to the k-th station. ek and Tfin are the centrod time of HFER, andthe signal duration used for determining stopping point, respectively. The terms marked by “fin” arecorresponding to those for the stopping point. C is P wave’s velocity here. As was mentioned by
Gusev and Pavlov (1991), it is difficult to determine both Mt and the shift of the centroid in thevertical direction simultaneously, and the term corresponding to the latter is neglected in the above
equations. In our analysis, we consider its effect to correct the observed centroid times by considering
arrival time differences from the hypocenter to the centroid.
In coordinates {x, y, z} following N, E and Z directions, using formulas:

r1k / c  dT / dx  cos( Az )dT / d 
r2 k / c  dT / dy  sin( Az )dT / d  ,(5)(6)

we calculate the terms rik / c neglecting the geometrical non-linearity of the earth (see Figure 1). Here and Az are station’s epicentral distance and azimuth ofthe main shock, respectively.
4. RESULTS AND DISCUSSION
4.1. Choosing EGFs
Choosing EGFs is a very important step in this study,
and should be carried out carefully. In this section, weFigure 1. Map view of the relative location
explain how we chose EGFs.
method. The red star denoates an epicenter andIn order to select EGFs, first, we compared
the purple circle represents a centroid of HFER
envelopes computed from seismograms observed at
or a stopping point. The blue arrow shows aeach station for some events that occurred in or near
unit vactor from epicenter to a station.
the source region.
We applied band-pass filtering (the fourfrequency bands), then we calculated theirP-waveenvelopes, and smoothed them (using5s’ smoothing window). The procedure we
conducted for the event pair between the 2007
and the 2009 earthquakes for station VLC isshown in Figure 2. We calculated all correlation
coefficients among all event pairs in eachfrequency band for each station.
Figure 3 shows that the correlation
coefficients between the smallest event andevents with large magnitudes are smaller thanFigure 2. (a) The observed waveforms of the 2007
and the 2009 events recorded at station VLC,
respectively. (b) Envelopes of the observed andfiltered waveform. (c) Smoothed envelopes.

3

those among other pairs. This suggests that the S/N ratio of the data from the smallest event is lower,
and that the HFER signals from relatively larger events are consistent with each other.
The data fromthe smallest event do
not have clear HF
signalsforsomestations due to highnoises, while large
magnitudes imply thatthey have long source
duration. Both of themcannot be good EGFs
inthisresearch.
Therefore, we chosethe 2009 earthquake as
the EGFs in this study.
4.2. HFER signals
determination
As we presented above,
raw HFER signals arealways distorted bywave scattering due toheterogeneities in theearth.
In order toFigure 3. The correlation coefficients of all four bands of the events pairs
obtaintheHF signals of
among the four eventsthe main shock, weapplied envelope inversion using EGFs chosen in the previous section, that is, the 2009 event.
We split the original waveform into four bands, which are 0.4-1.2, 1.2-2.0, 2.0-3.0, and3.0-4.0 Hz. The procedure we conducted is as follows:
1. Retrieve data from IRIS DMC;
2. Filter them into the four bands using band pass filters;
3. Perform envelope calculation;
4. Smooth them using 5s’ smoothing window;
5. Perform deconvolution using NNLS method;
6. Normalize the signals of the HFER.
First, we retrieved broadband vertical velocity data with epicentral distances between 30°and85°from IRIS DMC. Then, thefour bandpass filters were usedto obtain the signals in eachband. After that, envelope
calculation and smoothing
were done to the band passseismograms. Figure 4 shows
the centroid times and signal
durations for four bands
obtained for station VLC.
Note that there are
Figure 4. Schematic centroid times and signal durations for four bands
severalapparent
impulses inin station VLC
the obtained HF signals,

4which are likely to be strong energy radiation. This suggests that we may extract more informationfrom these results, not only the centroids of the HFER, but also the locations of sub-radiators.
4.3. Times and locations of the centroids and the stopping point of the HFER
In the preliminary analysis in which the effect of the vertical centroid shift is not considered, thecentroid location was determined
Table 2. The times and locations of the centroids and the stoppingnear the northwestern edge of thepoint of the HFER. Mt, M1, M2 are the average retardation of
rupture area suggested by previous
radiation respect to the origin time, relative locations of centroids ofHFER in the north south and eastwest directions, respectively. Thestudies. This indicates that theerrors are twice of their standard deviations.
vertical centroid shift is about 30km based on the assumed plateFrequency band
Mt
M1M2
boundary in this region. Therefore,
0.4-1.2
32.89±0.40
82.92±8.68-87.80±9.28
we corrected the observed centroid
1.2-2.0
29.69±0.44
78.71±9.60-48.89±9.90
and stopping times considering the2.0-3.0
30.00±0.50
83.44±9.60-65.30±9.90
vertical centroid shift of 30 km.
3.0-4.0
31.58±0.50
108.42±9.60-87.32±9.90
Figure 5 and Table 2 show theStopping point
59.67±0.42
75.35±8.98-148.87±9.90
centroids of the HFER signals, andthe stopping point of this event. It indicates that140 kmthe centroid times of the HFER are about 30 s
after the starting of the rupture, and the total
120
HFER duration is about 60 s. The locations of the100
0.4-1.2
centoids of the HFER are all huddled in the80
1.2-2.0
northwest of the epicenter, and the stopping point
2.0-3.060
is located in the even farther to the west. This
3.0-4.040
indicates that the rupture propagated northwest,
stopping point
and stopped at the northwest side of the epicenter.
20
Using the parameters we obtained from the0km
inversion, we calculated the theoretical centroid-180 -160 -140 -120 -100 -80-60-40-200times of the HFER for each frequency band in each
Figure 5. Spacial locations of the centroids thestation (Figure 6). We found overall agreement
starting point of the HFER. The red star shows
between the theoretically calculated and observed
the starting point of the rupture, and the error bars
centroid times, although there are considerable
show twice of the standard deviations.
scatteres between them.
Also, we found that the times, thatis, Mts, are very robust and keep consistenteven the damping factor and data checking
were not carefully treated. On the contrary,locations of the centroids of the HFER andthe stopping point are affected by data
selecting and the threshold we chose fordetermining time duration.
Several studies have already been
done on the source process of thisearthquake using different methods anddata sets. In terms of direction of therupture propagation, centroid times andsource duration, our results are consistentwith them. As for centroid locations, ourFigure 6. Comparison between the centroid times of thecentroids are located near the northwesternforward calculation using the parameters from theedge of the rupture area, which seems
inversion and those of the observed ones in the four
different from other results, most of which
frequency bands.

5suggests that the centroid of moment release is located around the center of the rupture area. Thisdifference may indicate the difference between slip distribution, which is determined by other studies,and high frequency energy radiation distribution, which is obtained in this study. However, since, as
was mentioned, the constraint for the centroid location is relatively weak, further analysis to improve
its accuracy is required to discuss this difference.
5. CONCLUSIONS
We determined the centroids and stopping point of the HFER of the 2003 Tokachi-Oki earthquake.
We used the empirical Green’s functions (EGFs) technique following the procedure of Gusev et al.
(2007). We found that the P-wave envelopes from small events in the nearby source region are similar
to each other for most of the stations, which demonstrated the applicability of this approach. Also, wefound that careful investigation and selection of EGFs are important for obtaining reliable results.
We performed deconvolution of the smoothed envelopes of the main shock using those ofEGFs to obtain the HFER signals in each band for each station. A non negative constraint is applied tothis calculation. Then, we calculated the times and locations of the centroids of the HFER and thestopping point. The centroid times and source duration are estimated to be about 30 and 60 s,respectively. The centroid locations suggest the northwestward propagation of the rupture. Theseresults are consistent with other results obtained using lower frequency data. Our results show that it ispossible to obtain the centroids of the HFER of large earthquakes using carefully chosen EGFs. Theresult shows that the centroid times can be determined reliably. The centroid locations are moredifficult to be constrained as in the same accuracy as that for centroid time.
Compared to methods using strong motion records, the advantages of this method are asfollows. First, due to the topographical features of some targets (such as that of in Tokachi-Oki), thestrong motion records cannot surround the epicenter. This will affect results of inversions. In themethod of the present study, we use tele-seismic waves, and the azimuth coverage of the data set isgood. Systematic bias due to poor azimuth coverage of the data set is expected to be small. Also, it iseasy to retrieve data and to implement the inversion. We plan to apply this method to a series of large
earthquakes and discuss the characteristics of the HFER.
6. RECOMMENDATIONA set of large earthquakes should be analyzed if we are going to explore the physical meaning ofHFERs. Furthermore, comparisons among earthquakes which are sorted by geodynamic environments
are appreciated and will be one in the list of our future plan.
AKNOWLEDGEMENT
I would like to express my sincere gratitude to Dr. N. Hurukawa and Dr. B. Shibazaki for their usefulcomments in improving this study and manuscript.
REFERENCES
Gusev, A. A., and Pavlov, A.A., 1991, Pure Appl. Geophys., 136, 235–244.
Gusev A.A., Guseva, E. M., and Panza, G. F., 2007, Geophys. J. Int., 170, 1119-1128.
Hara, T. and Yokoi, T., 2008-2009, Lecture notes on Data processing, IISEE/BRI.
Irikura, K., 1986, 7th Japan Earthq. Eng. Symp., Tokyo, 151-156.
Okamoto, T., and Takenaka, H., 2009, Earth Planets and Space, 61, e17-e20.
Web site: Global CMT web page, http://www.globalcmt.org/.
Web site: Incorporated Research Institutions for Seismology (IRIS), http://www.iris.washington.edu.

6