Comparison of ground truth location of earthquake from InSAR and from ambient seismic noise: A case study of the 1998 Zhangbei earthquake

1.   Introduction

Earthquake locating is one of the first order issues in study on earthquake and seismic hazard. The improvement of location precision of natural and artificial events is helpful to studies on Earth structure. So seismology community makes large effort to provide a standard database of reference events whose location qualities are guaranteed. IASPEI collects and validates many events with high accuracy ( <10 km) provided by global institutes named after ground truth event (Bondár and McLaughlin, 2009). Ground truth events (referred to as GT, hereafter) are employed as input data of geophysical inversion and benchmark in model validation (Yang et al., 2004; Murphy et al., 2005). The calibrations obtained from GT events are also used to help improve the event locations (Flanagan et al., 2007). Most of the GT0/GT2 events are nuclear and chemical explosions whose quantity is limited. The tra- ditional methods used to obtain natural GT events depend on dense seismic network close to the hypocenter. Bondár et al. (2008) proposed a hybrid multiple events location method to obtain GT events which is widely adopted and provides most GT events. This kind of methods based on seismic traveltimes normally require close stations to reduce the error caused by 3D heterogeneity or cross-over of Pn and Pg. Recently, numerous groups used InSAR method to reveal the near source displacement field of large and moderate earthquakes (Massonnet and Feigl, 1998; Dawson et al., 2008.; Zha et al., 2008; Fialko et al., 2005). This geodesy method is not only used to study the focal mechanism, but also provide high precision centroid location. Dawson et al. (2008) identified a very shallow moderate earthquake in west Australia. The synergy of seismic and InSAR method is also used to provide GT events in Asia and Middle East (Saikia et al., 2002). But this method does not work well in rich vegetation area, or small (M_W < 5.0) and deep events. The perturbations in atmosphere and post seismic deformation also affect the precision of the InSAR method.

Besides the travel times of seismic body wave phases, the waveforms of body waves and surface waves also includes information about the source parameters including location (Ekström, 2006). Tan et al. (2006) proposed a hybrid method which used P wave traveltimes and calibration from surface waves to obtain GT5 events. Since the noise cross-correlation provides high quality EGF of surface wave between two sites (Weaver and Lobkis, 2001; Shapiro and Campillo, 2004; Snieder, 2004), two groups (Zhan et al., 2010; Barmin et al., 2010) used the noise cross-correlation functions (referred to as NCF, hereafter) between remote stations and station which is very close to the earthquake as calibration, then provided good centroid locations. Both approaches are based on waveform which is affected by focal mechanism uncertainties. We proposed a new traveltime-based method to obtain GT event and the 1998 Zhangbei earthquake was taken as example.

The M_W5.7 Zhangbei earthquake occurred at UTC 03:50:46.4 on January 10, 1998, with centroid (41.34°N, 114.34°E) (www.globalcmt.org, referred to as GCMT, hereafter) where is only 260 km away from Beijing. It is the strongest earthquake in northern China region after the 1976 Tangshan earthquake and caused 49 deaths, about 0.836 billion RMB economic loss (Cai et al., 1998). There is no historical damaging earthquake in this region. Because of the importance of this event, the source parameters were determined with different datasets (Gao et al., 2002; Yang et al., 1999; Zheng et al., 1998; Wang et al., 2000; Zhang et al., 2002; Shan et al., 2002). But the results are not consistent with each other in both focal mechanism and epicenter (Figure 1a). The accuracy of seismological results with traveltime method is limited by the available data recorded by close stations. The location of maximum intensity which is close to the "Gao aftershock" in Figure 1 is about 25 km far away from the location published by global CMT. Two groups using InSAR data also provided different locations with larger separation than the usual error (Zhang InSAR in Figure 1a, Zhang et al., 2002; Li InSAR in Figure 1a, Li et al., 2008). Two blind faults are discovered by satellite image (Xu et al., 1998) which crossover each other in the north of epicenter. It is still under controversy on which fault the earthquake occurred (Gao et al., 2002; Lai et al., 2007; Li et al., 2008), though Li et al. (2008)'s result prefers to north-east striking fault plane. Indeed it is difficult to resolve which plane is the ruptured plane for such moderate earthquakes with seismology only, except with small aftershocks as empirical Green's function (Luo et al., 2010). In this paper, we do not concern with which fault plane is the ruptured plane, and only employed data record at the closest permanent station ZHB and five IC/IRIS stations are used to compare the location with InSAR result.

Figure  1.  (a) The locations and mechanism obtained from different datasets and methods. (b) The stations location and the 20 s Rayleigh wave group velocity in East Asia. GCMT: global CMT solution, Li InSAR: result of Li et al. (2008), Zhang InSAR: location provided by Zhang et al. (2002), Gao aftershock: location provided by Gao et al. (2002).

DownLoad: Full-Size Img PowerPoint

2.   Data and methods

We collecte the continuous data of BHZ (a permanent station of Hebei Province Seismic Network) channel of ZHB and LHZ channel of five stations of IC/IRIS (www.iris.edu) during 2008 and 2009 (Figure 1a). The records of the earthquake are also ready after removing the linear trend, mean value and instrument response (Figure 2a). For comparison with the EGF obtained from ambient seismic noise, the records are within the band between 0.03 and 0.1 Hz in which band both the ambient noise and crustal earthquake signals are strong.

Figure  2.  (a) The record section of 1998 Zhangbei earthquake. (b) The NCFs between IC/IRIS stations and ZHB station. The solid line denotes the segments used in differential time measurements in following.

DownLoad: Full-Size Img PowerPoint

The procedure of calculation NCF is quite similar to that of Bensen et al. (2007). After removing the linear trend, mean value, instrument response, the records are divided into two hours segments. To reduce the earthquake signals, a running average absolute method is applied. The normalization weight function is defined as the envelope of record which is band-passed between 0.02 and 0.066 Hz. In spectral whitening, we set the pass band between 0.03 Hz and 0.1 Hz. After single station data preparation mentioned above, the cross correlation between two stations are obtained from a frequency domain method. The final NCFs are obtained from stacking cross correlations of every two hours segments. We use symmetry energy of NCF to improve the signal to noise ratio (referred to as SNR, hereafter) (Lin et al., 2008) (Figure 2b).

As well known, the dominant challenge of location with remote stations is the velocity heterogeneity of the Earth interior. There are two approaches to reduce such structural error. Richards et al. (2006) used the differential times between multiple events by cross-correlation of Lg wave to improve the relative location. The other ones aimed to obtaining high accuracy absolute location with empirical calibration (Zhu et al., 2006) or calibration derived from 3D model (Yang et al., 2004). Bermin et al. (2010) constructed the synthetic NCF between hypothetical location and remote stations from the NCFs of many neighboring stations of source and remote stations, and then the arrival time was obtained by cross-correlation of synthetic NCF and earthquake records which is normalized by NCF amplitude spectra. This method requires the dense temporal/permanent network in source region but it could reach a high accuracy ( <1 km). Zhan et al.(2010) proposed another approach which invokes the synthetic earthquake and NCF waveforms and it has a double difference form. It only requires a temporal station nearby the source and the accuracy reaches about 2 km with well established focal mechanism. We test this method in the following.

2.1   Waveform-based calibration

The idea of NCF calibration is quite simple. If there is a station close to the source, the propagating paths linking between source and remote station and between close station and remote station are very close too. The structural errors of 1D model are almost identical. This error could be considered as the misfit between observations (earthquake records and NCFs) and synthetic waveforms with 1D model. For NCF, the difference between observed and theoretical ones reflects the 3D structure effect and it is employed as calibration. Besides the structural error, the synthetic waveforms of earthquake are also affected by two error terms about earthquake: location and focal mechanism. After the reliable focal mechanism is ready, the remaining error is due to location. It is very convenient to use a grid search to obtain an optimal location with the minimum misfit. The differential time is obtained by the timeshift cross-correlation scheme in the same frequency band. In summary, this method could be explained as equation (1). In equation (1), the subscript "eq" denotes the earthquake terms while the NCF terms are marked by "ncf". The differential time "dT " is the time shift which the cross-correlation between two waveforms reaches the maximum.

The 1D model (Table 1) used in calculation of synthetic waveforms is extracted from CrustModel 2.0 (http//igppweb.ucsd.edu/∼gabi/rem.html). A vertical point force source with duration of 15 s at ZHB station is used to calculate the synthetic NCFs waveforms. The differences between focal mechanisms published by Global CMT and University of Tokyo, Ma et al. (1998) are very small (10). We adopted the Global CMT solution (strike: 207°, dip: 54°, rake: 135°). But the teleseismic surface wave is difficult to constrain the source depth, so the source depth is set at 5 km which is confirmed by InSAR and aftershock distribution.

Table  1.  The 1D model

 | Show Table

DownLoad: CSV

Since the synthetic waveforms of both NCF and earthquake are ready, we chose the Rayleigh wave windows which are identified by the PI/2 phase difference between R and Z components in synthetic waveforms and earthquake records (Figure 3). Figure 4 shows the waveform comparison between observations and synthetic ones with the differential times. The crosscorrelation of NCFs ranges from 0.46 to 0.88. The low cross correlation coefficient of MDJ and SSE due to low SNR in NCFs. The separation of these two station-pairs is larger than others and both two stations are close to the sea which generates a lot of seismic noise, so the low SNR is understandable. A grid search method is applied to the area around the GCMT location. The differential time caused by the location movement is computed with a numerical gradient. We calculate the synthetic earthquake waveform at a new location and the differential time are obtained by the same cross-correlation method mentioned above. Then the gradient could be represented as the differential time divided by location movement. First we use the differential time between observed and synthetic earthquake waveform to relocate. The optimal location is quite close to the GCMT location (∼10 km). The minimum time misfit is about 12 s (Figure 5). Then the NCF calibration is employed. The new location move to the southeast and is about 13 km far away from the location provided by InSAR study. The obvious benefit of NCF calibration is the time misfit which is reduced to 6.2 s (Figure 5).

Figure  3.  (a) The synthetic waveforms of earthquakes. (b) The synthetic waveforms of NCF. The solid line denotes the segments used in differential time measurements in following.

DownLoad: Full-Size Img PowerPoint

Figure  4.  The comparison between observed (solid lines) and synthetic (dashed lines) waveforms of earthquake (a) and NCF (b). The waveforms have been aligned with the shifted time of maximum cross-correlation.

DownLoad: Full-Size Img PowerPoint

Figure  5.  (a) The location result without NCF calibration. The optimal location with minimum misfit is marked with "1D" square. The contour lines denote the distance far from the InSAR location with interval of 5 km. (b) The location result with NCF calibration. Both (a) and (b) use the same color bar of misfit function.

DownLoad: Full-Size Img PowerPoint

2.2   Traveltime-based method

The arrival time of surface wave group is almost independent of the focal mechanism, as Levshin et al. (1999) found that source term for such moderate earthquake is not large. Using the group arrival time, the uncertainty of focal mechanism will not take error into the relocation. We define the new NCF calibration as the group velocity of inter-station Rayleigh wave in NCF. This velocity is equal to the average velocity in path linking ZHB and stations which is almost the same as the path linking between source and stations. The principle of this calibration is quite similar to the waveformbased one, but the uncertainty of focal mechanism has been eliminated. The group travel time is calculated by the distance divided by group velocity. We used the multiple filter method to measure the Rayleigh wave group velocity (Hermann and Ammon, 2002). Considering the different spectra of earthquake and ambient seismic noise, we use the frequency bands in which both earthquake and noise signals are strong. The bands are shown with squares in the Figure 6. The low SNR makes it unstable to measure the group velocity of inter-station Rayleigh wave with NCFs of ZHB-MDJ and SSE. So we only used the other three stations. At BJT, the most time differences at different periods are similar, which suggests the dispersion of group velocity is not strong in this narrow band and the difference between the location of ZHB and the source is dominant factor (Figure 6). The similar phenomenon is also observed at XAN. The dispersion effect is obvious at HIA, which causes the difference time ranges from positive to negative at different periods. The new misfit is defined as equation (2), where X and R denotes the location of hypothetical source location and station. T means the group traveltime while U denotes the group velocity.

Figure  6.  The group arrival times (see detail of definition). The solid lines show the ones of earthquake while the ones of NCF are shown in dashed lines. The squares show the period range used in this study.

DownLoad: Full-Size Img PowerPoint

The group arrival time of earthquake record is controlled by the location and origin time. Therefore, new grid search scheme will search three parameters. Finally we obtain a new location at only ∼3 km west to the InSAR location (Figure 7, Li et al., 2008) that is also consistent with aftershock distribution by dense temporal seismic network (Gao et al., 2002). The average misfit time is less than 4 s.

Figure  7.  The location result with NCF traveltime calibration. The optimal location is at the dot. The result obtained from waveform based method is denoted with ZWN calibration.

DownLoad: Full-Size Img PowerPoint

3.   Discussion and conclusions

In this paper, the waveform-based method provided a result with about 13 km error while the traveltimebased method works much better (error ∼3 km). There are two potential error sources in the waveform-based method. First this method requires a well known focal mechanism. The Zhangbei earthquake is an unusual case of which the global CMT solution is obviously different from the InSAR ones in focal mechanism. The large difference (∼40°) in rake angle caused more than 3 s error in most stations except MDJ, which could not be ignored comparing with the average misfit of 6 s in our result. The second one is NCF itself. In this case, the separation of station-pairs is much larger than the previous application of the method (Zhan et al., 2010). Although the NCFs during longer time duration is stacked to improve the SNR, it is more difficult to extract reasonable inter-station Rayleigh wave in short period. If the azimuth gap is large, the error of earthquake origin time will take error into relocation. For example, the optimal origin time obtained by traveltime-based method is about 6.6 s. The gap of our dataset is about 180° that may cause a shift of about 20 km to southeast. In contrast, the traveltime-based method is independent of focal mechanism and takes the origin time into inversion. But the different spectra of earthquake and ambient seismic noise still affect the available measured data. For example, the Rayleigh wave at BJT is quite weak at periods larger than 10 s which is affected by the radiation pattern of Rayleigh wave. In NCF of stationpair with large separation, the short period inter-station Rayleigh is difficult to obtain. So the precise location of small earthquakes may be difficult to obtain with this method.

In summary, two NCF calibration methods are employed to study the centroid of the 1998 Zhangbei earthquake. Our results confirm the location revealed by InSAR study (Li et al., 2008). The waveform-based method with five stations provides comparable location precision by global seismic network. We proposed a new NCF calibration based on traveltime to reduce the error caused by uncertainty of focal mechanism. The location is only 3 km far away from InSAR location which is consistent with maximum intensity area of this earthquake. This method could provide GT5 events with sparse regional seismic network. Such GT events could be used as reference event in real time seismology (Wan et al., 2009).