1 INTRODUCTION

Earthquake location plays a fundamental role in seismic researches on seismicity, nuclear test monitoring and Earth’s interior structure (Zeng, 1991; Ding et al., 1999; Richards et al., 2006; Chen, 2009). To accurately determine sources’ positions and origin times, seismologists have developed various location methods (Tian et al., 2002; Cai et al., 2004; Yang et al., 2005), such as the Geiger method (Geiger, 1912), the double-difference algorithm (Waldhauser et al., 2000), the grid-search scheme (Shearer, 1997), the intersection method (Pujol, 2004; Lian et al., 2006; Zhou et al., 2012) and the time-reversal imaging technique (Xu et al., 2013a and 2013b). These location methods have their own characteristics. Compared with others, the intersection location method has advantages including excellent simplicity, visuality and robustness, and thereby is widely used in seismic networks (Meng et al., 2001; Song et al., 2001). Its basic idea is to determine the hypocenter using hypocentral loci. As a result, valuable information on the hypocenter can be obtained even when there is a few of seismic records available for use (Zhao et al., 2008). On the other hand, the intersection method is deficient in location accuracy, especially for focal depth (UdÍas, 1999). This is mainly attributed to the fact that the reference model is assumed to be homogeneous or laterally homogeneous (the focal loci in these cases have analytical solutions). In the real earth, however, there are large lateral as well as radial velocity variations. The discrepancy between the reference model and the real Earth inevitably results in different degrees of source bias (Billings et al., 1994; Chen et al., 2006). So, modifying the intersection method to be suitable to complex velocity models is critical for improving its location accuracy. Besides, highly accurate location based on 3-D models has become a reachable target because more and more regional three dimensional velocity models are constructed with considerable accuracy and the resolution of global structure image is continually improved, which greatly benefits from the development of seismic networks and inversion techniques.

Expansion of the conventional intersection method to a three dimensional velocity model firstly involves the calculation of hypocentral loci, which can rarely be expressed analytically when the velocity varies in three dimensions. As for the problem of calculating focal loci in complex media, there are presently two main schemes. One is representing a hypocentral locus with the points that satisfy the locus equation within a prescribed tolerance (Zhou, 1994; Font et al., 2004; Theunissen et al., 2012). This scheme is straightforward and adaptive, and hence has no restrictions on the geometry of hypocentral loci. However, resultant rough loci do not intersect at one point but to form a region even when the velocity model and picked arrival times are both accurate, especially for a large tolerance. The other one is based on a minimum travel time algorithm for tracing rays: A hypocentral locus is represented by the ray paths from the minimal point (namely the initial point) to low valued points (referred as hypocentral locus reference points, HLRPs) in its residual field (Zhao et al., 2008). The ray-tracing based scheme can produce such fine hypocentral loci that their directional distribution allows us to directly estimate the location quality and visually analyze the effects of the velocity model, the observed arrivaltime and the station distribution on the location result. Its major disadvantages are that (1) it can not correctly deal with multi-segmented loci: When a hypocentral locus consists of two or more separate segments, the scheme will produce false locus, namely, the rays connecting individual locus segments because only the minimal point in the residual filed is set as the initial point which is connected to all the HLRPs by rays; (2) it does not satisfy the requirement of processing large batches of data because the selection of HLRPs is rather laborious: HLRPs are chosen at a percentage of the number of all model nodes, and the percentage is required to vary with the length of hypocentral locus while the locus length can hardly be known in advance; (3) it is poorly adapted for producing both fine and complete loci: Even though the locus length is known, HLRPs selected purely according to the residual value are either so insufficient that they do not characterize the locus completely or so excessive that the resultant locus is rough because the residual differs at the nodes of the elements that a hypocentral locus goes through. Actually, earthquakes recorded by and observed data from modern seismic networks are both very large in quantity. For accurately locating large numbers of seismic events, hypocentral loci are required to be calculated with high efficiency and accuracy.

Thus it is clear that the solution for calculating hypocentral loci in complex media needs to be improved further. In this study, we attempt to modify the calculation method based on a ray tracing algorithm (Zhao et al., 2008) to overcome its disadvantages. For convenience, the methodology is illustrated using a 2D model.

2 EQUATIONS OF HYPOCENTRAL LOCI

Presume that the velocity model of a region is known and on the ground there are M seismic stations, R_i(i=1, 2, 3, …, M). To locate earthquakes within the region, we divide the region’s velocity model into equal square elements. Within an element, the seismic velocity is constant and the central point of the cell is a node of the model. We employ a minimum traveltime tree algorithm for tracing rays (Zhao et al., 2003, 2004) to calculate travel times of the first-break primary (P) and shear (S) waves from a station R_i(i=1, 2, 3, …, M) (the shot points are placed at stations) to each node with coordinates x=(x, z) of the region, T (x; R_i, P) and T (x; R_i, S)(i=1, 2, 3, …, M). As for the earthquake with the hypocenter at x₀=(x₀, z₀) and the origin time of t₀ in the region, the P-and S-wave first arrivals from the hypocenter to the seismic station R_i are denoted by τ_Ri^P and τ_Ri^S respectively. According to the hypocenter-station symmetry, we can construct two kinds of equations of hypocentral loci constrained with arrival time differences and arrival times, respectively.

2.1 Hypocentral Loci Constrained with Arrival Time Differences

Every two observed arrival times constitute one arrival time pair. For the jth arrival time pair, its two components are respectively denoted by τ_Rj1^Wj1 and τ_Rj2^Wj2, where R_j1 and R_j2 represent the seismic stations recording the two arrival times; W_j1 and W_j2 are the types of seismic phases, i.e., P or S. The hypocentral locus constrained with the difference between τ_Rj1^Wj1 and τ_Rj2^Wj2 meets the following equation (Zhao et al., 2008):

(1)

The equation above indicates that the focal locus is of quite small values near zero, taking calculation errors into account, in the following residual field with respect to travel time difference:

(2)

where I(x)=1.

2.2 Hypocentral Loci Constrained with Arrival Times

Using the observed arrival times, we can construct an origin time field of the event as follows (Zhao, 2014):

(3)

where G_j and W_j are respectively the station and the seismic phase type (P or S wave) of the jth arrival time, and K is the number of the arrival times observed. Hence, the hypocentral locus constrained with the jth arrival time meets the following equation:

(4)

The Eq.(4) indicates that the focal locus is of quite small values near zero, taking calculation errors into account, in the following residual field with respect to arrival time:

(5)

where I(x)=1.

Equations (1) and (4) have exact analytical solutions only when the velocity model is very simple. Different from hypocentral loci constrained with arrival time differences, the ones constrained by arrival times are rather complex even for homogeneous media, especially when the origin time field is constructed using many arrival times. Let us consider a homogeneous model with the dimensions of 150 km×90 km. The P-and S-wave velocities are 6.0 km·s⁻¹ and 3.6 km·s⁻¹, respectively. We divide the model into square cells of 3 km×3 km. Assume that an earthquake is triggered at x₀(37.5 km, –31.5 km) and t₀=0 s. The P-and S-waves from the event are recorded by five stations, R_i(13.5 + (i − 1) × 30.0 km, –1.5 km) (i=1, 2, 3, …, 5). Fig. 1 shows the origin time fields of the event and the arrival time residual fields corresponding to the P-wave arrival at the station R₃. The grey scale of Fig. 1c and 1d is proportional to the residual value. In order to clearly show the focal loci, the residual values have been normalized and converted into base-10 logarithms. From Fig. 1, it can be seen that 1) the origin time field constructed with only P-wave arrivals (Fig. 1a) is similar in structure to that with both P-and S-wave data (Fig. 1b), which results from the constant velocity ratio of the two modes of waves; 2) the hypocentral locus constrained with only P waves (Fig. 1c) and that with both P-and S-waves (Fig. 1d) are not circular or hyperbolic, but much more complex, especially the latter consists of two separate segments; 3) for different origin time fields, hypocentral loci constrained with the same arrival time are not same, even are entirely different.

3 CALCULATION OF HYPOCENTRAL LOCI

In a residual field described by Eq.(2) or (5), points on the hypocentral locus bear extremely small values, as Fig. 1 shows. The low residual (e.g., smaller than 0.2 s) area is divided in accordance with the geometry of the hypocentral locus into one (Fig. 1c) or more (Fig. 1d) connected regions. Assume that sources are exploded simultaneously at the minimum points in these separate connected regions. The seismic waves certainly will travel fast along the focal locus, but spread slowly in areas beyond it, especially outside the connected regions, where the field value is regarded as that of seismic slowness. We characterize the hypocentral locus using some nodes (referred as the reference points of the hypocentral locus, HLRPs), and assign the point with minimal residual of every connected region as an initial point. In this way, ray paths from the initial points to their respective HLRPs that lie in the same connected regions with themselves approximately represent the hypocentral locus, and there are no false hypocentral loci. For the convenience of treatment, nodes with residual below a given threshold are usually assigned as hypocentral locus reference points (Zhou, 1994; Zhao et al., 2008). As is illustrated by different gray scales in Figs. 1c and 1d, the residual varies considerably at the nodes of the cells that the hypocentral locus runs through. Therefore, it is rather difficult to obtain both fine and complete loci. When nodes of the cells crossed by the hypocentral locus are chosen as the HLRPs, the resultant hypocentral locus is quite fine and complete. Based on this analysis, we propose an improved method for calculating hypocentral loci in complex media.

3.1 Selection of Reference Points

The hypocentral locus lies in the middle of connected regions with low residuals, just like one or more cores of a fruit, as can be seen from Figs. 1c and 1d. As for nodes of the cells that the hypocentral locus goes through, their residuals are different considerably but are respectively minimal, compared with ones at their two sides beyond the locus. According to the characteristic of such nodes, we can select them in a way similar to the one that we get the core of a fruit by repeatedly cutting off its outside layer. For the sake of description, nodes with residuals below R_FLRP are referred to as the primary HLRPs and their set is denoted by Θ. The scheme for selecting HLRPs from the set Θ can be expressed as follows:

(1) Classification

A node i ∈ Θ is classified as one inner point of the set Θ if its adjacent points all belong to Θ; otherwise it is classified as one outside point of the set Θ. The sets of inner and outside nodes are denoted by I and O, respectively.

(2) Updating

A node i ∈ O is retained if among its adjacent points there are no or only such inner points that they bear bigger residuals than the node i; otherwise it is removed from the set O. After this operation is done for all members of the set O, the set Θ is updated as the union of the sets O and I.

(3) Iteration check

If the set I becomes empty go to (4); otherwise go to (1).

(4) Simplification

The set O sometimes includes such members that they are distributed in two lines side by side, along the hypocentral locus. In this case, we simplify the set O, for only one of the two lines is required for characterizing the locus. As for the node i ∈ O and its adjacent outside points, when they are more than three, we retain only the three minimum-residual ones.

(5) Purification

When the residual field is complex, the set O may include a few of outliers, which lie away from the hypocentral locus. In this study, a node with the residual larger than the average of all the points in O by two standard deviations is regarded as one outlier, and removed from O.

In the algorithm above, step (1)-(3) stands for the layer-by-layer removal of the outside members of the set of primary HLRPs in a similar way of peeling a fruit. Step (4) stands for the removal of superfluous reference points and step (5) for that of false ones. The points retained finally in the set O are used as HLRPs. In Figs. 1c and 1d, the yellow points are HLRPs selected from the primary ones with residuals of less than 0.2 s. It can be seen that the selected HLRPs characterize the hypocentral locus fairly well and hence are sufficient for calculating the locus although they are just a part of nodes that the hypocentral locus goes through. Actually, for the calculation of a hypocentral locus represented by ray paths, it is preferable to choose only the minimal residual node and the two end ones of each segment of the locus as the HLRPs but is difficult to realize.

3.2 Calculation of Ray Paths

Ray paths representing the hypocentral locus are calculated with an improved minimum traveltime tree algorithm (Zhao et al., 2004). The minimum traveltime tree algorithm for tracing rays (Nakanishi et al., 1986; Moser, 1991) is based on Huygens’ Principle and Fermat’s Principle in Seismology, and has an advantage of excellent robustness suitable for complex velocity models. Travel times and ray paths from the shot point to all nodes in the model are calculated simultaneously by constructing a minimum traveltime tree (Cao et al., 1993). As for a hypocentral locus consisting of multiple segments, in order to prevent from producing false hypocentral locus, the minimal residual point in each connected region that contains one segment of the hypocentral locus is required to be assigned as an initial point of seismic waves. Such initial points of all the connected regions are difficult to determine at the same time because the distribution of the hypocentral locus can hardly be known in advance. On the other hand, the minimum traveltime tree algorithm runs within a connected region, which is concluded by its ending condition. So, we restrict the calculation of ray paths to a low residual area and then employ the minimum traveltime tree algorithm to calculate the hypocentral locus by assigning in turn the initial points of the connected regions. As HLRPs include the minimal residual points of all the connected regions, initial points of seismic waves can be chosen from them.

The calculation domain of ray paths is defined as the area Ω with the residual of less than R_FL, which can be realized by setting the nodes outside Ω as the points at which the travel time has been minimum. In this subsection, the set of HLRPs selected from the nodes with the residual of less than R_FLRP is denoted by Z. The parameters R_FL and R_FLRP can be set same or different. The algorithm for tracing rays representing a hypocentral locus can be described as follows:

(1) Selection of the initial point

Find the node s ∈ Z with the minimum residual and remove it from Z.

(2) Calculation of ray paths

Take s as the initial point and run the minimum traveltime tree algorithm to calculate travel times and ray paths of the seismic waves in the residual field until it stops; if a node j ∈ Z has been traced then remove it from Z.

(3) Iteration check

If Z becomes an empty set then stop; otherwise go to (1).

3.3 Modification of Hypocentral Loci

Most of the reference points of a hypocentral locus do not exactly but approximately satisfy the locus equation. When their amount is excessive, the calculated focal locus is probably rather crude and needs some modifying.

Ray paths from the HLRPs to their perspective initial points of seismic waves form one or more tree structures, as Fig. 2 shows. Thus, we regard the rays as tree branches (directed line segments in Fig. 2) and the nodes that they go through as branch origins. A branch origin is also one branching point when the number of branches sprouting from it (figure in the circle) is not less than two. The trunk of the ray tree (namely the part constituted by the rays from the initial point to the HLRPs that there are not less than n_FLP rays going through, such as line segment BCD) certainly lies on the actual hypocentral locus. However, selecting the part only (Zhao et al., 2008) probably harms the intactness of a hypocentral locus. For example, hypocentral locus segment GHI in Fig. 2 will be lost for n_FLP=2. Therefore, for the sake of the completeness and the fineness of the calculated hypocentral locus, we remove the twigs without branching and with a length less than L_twig, such as the dashed lines in Fig. 2.

The blue lines in Figs. 1c and 1d are the corresponding calculated hypocentral loci for n_FLP=2 (in unit of the length of the element side). Clearly, the calculations are quite consistent with the real loci (the red dashed lines). Certainly, the consistence between them is better when the model is divided into smaller (e.g., 1 km×1 km) cells.

4 EXAMPLES OF CALCULATION 4.1 Virtual Event 4.1.1 Velocity model

The virtual event is designed against a background of Yunnan area, China. Yunnan area is characterized by active seismicity, and its deep structure has been investigated by many geophysical explorations (Kan et al., 1986; Wang et al., 2002). As a result, there are a number of high-resolution velocity profiles constructed from the deep seismic sounding data (Bai and Wang, 2003; Zhang et al., 2011; Teng et al., 2013). We slightly smoothed the Zhefang-Binchuan velocity profile (Bai and Wang, 2003) and obtained a model as shown by Fig. 3. The P-wave to S-wave velocity ratio vS : vP is taken as a constant 1:1.73. It is clear that the model is strongly heterogeneous in lateral and vertical directions. We divide the model into 1 km×1 km elements and set the velocity constant within an element. In the top surface, there are 9 stations R_i(i=1, 2, 3, …, 9), whose positions are listed in Table 1. Assume that an earthquake occurs at x₀(130.5 km, –27.5 km) with the origin time t₀=0 s. Traveltimes and ray paths of the P-and S-waves from the hypocenter to the stations are calculated with an improved minimum traveltime tree algorithm (Zhao et al., 2004). The ray paths are shown as the black lines in Fig. 3.

As for the earthquake in Fig. 3, we apply the modified method to calculate its two kinds of hypocentral loci, which are constrained with arrival times and arrival time differences, respectively. Arrival time difference constrained hypocentral loci each involve only two arrivals and hence are fairly independent. They have also been studied in detail (Bai et al., 2003; Lian et al., 2006; Zhao et al., 2007, 2008, 2009; Theunissen et al., 2012). Contrastingly, arrival time constrained ones are not only dependent upon their respective arrival times but also related to those constructing their origin time fields. Furthermore, they are poorly researched. Therefore, as for the former, we calculate them in only the case of all stations, a correct velocity model and precise arrival times; whereas, as for the latter, we investigate them in different observation cases (for example, all stations, sparse stations, near stations, far stations and one-side stations) when the velocity model and the arrival times are all correct and in the case that all the stations are employed, the velocity model is accurate but the arrival times are randomly perturbed. All the hypocentral loci are calculated with a set of parameters: R_FLRP=0.1 s, R_FL=0.2 s, L_twig=5.

4.1.2 Hypocentral loci constrained with arrival time differences

In theory, P-wave arrivals at the 9 stations can be pairwise combined into C₉²=36 arrival time differences. Among the hypocentral loci constrained with these differences, some ones are much alike in directional distribution. Herein we select ten representative hypocentral loci as shown in Fig. 4a. The character pair “i − j” beside a locus signifies that the locus is constrained with the arrival time difference τ_Rj^P − τ_Rj^P and the signal of red “+” represents the real hypocenter. One can see that the resultant hypocentral loci are all fine and quite complete without adjusting the calculation parameters; their directional distribution is fairly isotropic and complete near the hypocenter and results in a tight constraint upon the hypocenter; as for station twins symmetrical with respect to the epicenter, the arrival time difference constrained hypocentral loci (i.e. locus “1–9”) trend in nearly vertical directions near the hypocenter, and tightly constrain the epicenter while as for those away from the epicenter, the corresponding ones (i.e. locus “1–3”) have strong horizontal components near the hypocenter, and appear well capable of constraining the hypocentral depth. The perfect intersection of hypocentral loci at the hypocenter results from the correctness of the velocity model and the arrival times.

Hypocentral loci constrained with arrival time differences between P and S waves at the same stations are shown in Fig. 4b. The character “i” beside a locus expresses that the locus is constrained with arrival time difference τ_Ri^S − τ_RI^P. Obviously, the loci corresponding to stations R₂ and R₇ respectively consist of two separate segments, and both have been correctly calculated without producing false loci even though the velocity model is fairly complex. Similar to the hypocentral loci constrained with arrival time differences between P waves at different stations, those constrained by P and S wave differences at the same stations also have a complete directional distribution and sufficiently constrain the hypocentral location, but different from the former, their near-station loci tightly constrain the hypocentral depth and their far-station ones firmly constrain the epicenter.

When the observation is not complete, for example, only stations R₅~R₉ are employed, neither of the two types of hypocentral loci adequately constrains the hypocenter because of the loss of many distribution directions. On the other hand, they are strongly complementary in directional distribution and thus combining them can significantly reduce the effect of incomplete observation. When the arrival times and/or the velocity model are not accurate, it can be concluded that the resultant hypocentral loci will not cross exactly at the hypocenter but intersect into an area (See Zhao et al., 2008).

4.1.3 Hypocentral loci constrained with arrival times

Hypocentral loci constrained with arrival times of P waves at all the stations are shown in Fig. 5a. Their directional distribution near the hypocenter is not as complete and isotropic as that of ones constrained with arrival time differences (Fig. 4a). They mainly have strong vertical components near the hypocenter but there are two going through the hypocenter at large angles from the vertical direction. As a result, they also sufficiently constrain the hypocentral location. When the origin time field is constructed with both P-and S-wave arrivals, the arrival time constrained hypocentral loci are shown in Fig. 5b. The blue and the green lines correspond to the P-and S-waves, respectively. As can be seen, the additional use of S-wave data resulted in a slightly improved azimuthal distribution, which is much similar to the previous one. But then it is beneficial to decreasing the influence of random noise because more hypocentral loci are available to constrain the hypocenter.

In the case of all stations, the station interval is 30 or 20 km, whereas in practice it is two or more times longer than this distance. Hence we reserve only odd stations R_2i−1(i=1, 2, 3, …, 5). Fig. 5c and 5d show the hypocentral loci constrained with alone P-wave arrivals and with both P-and S-wave ones at these sparse stations, respectively. As one can see, the hypocentral loci resulting from the sparse observation have a somewhat narrower azimuthal distribution near the source, but still sufficiently constrain the hypocenter position. Similarly, there are hypocentral loci constrained with the arrival times at near stations (R₂~R₆), which are shown in Fig. 5e and 5f. When there are only far stations (R₁, R₂, R₇, R₈ and R₉), the arrival time constrained hypocentral loci for P-waves and for both P-and S-waves trend in nearly vertical directions near the hypocenter; thereby they can hardly constrain the depth of the focus.

In the practice of location, events often take place outside the seismic network. To simulate this situation, we keep only five stations R₅~R₉ at the right side of the epicenter. Hypocentral loci constrained with the P-wave arrival times are shown in Fig. 5i. Notably, they trend mainly in the horizontal direction, and appear capable of tightly constraining the hypocentral depth. When the S-wave data are additionally used, the directional distribution near the hypocenter becomes sharply different from those in the previous cases. It is both quite complete and isotropic, as Fig. 5j shows. This implies that the joint use of P-and S-wave arrival times is especially important for the earthquake location under the condition of partial observation. In Fig. 5j, it also can be seen that not only all the hypocentral loci intersect exactly at the hypocenter but also partial ones intersect in the area up and right of the hypocenter. From this we conclude that the function of total arrival time residual has local minimum points. Accordingly, when an earthquake is located using the Geiger method (its objective function is total arrival time residual), the input initial values of hypocentral parameters are required to be close to the real ones.

Arrival times always include more or less errors for they are affected by various factors such as environment, instruments and operators. To investigate the effect of arrival time picking errors on the distribution of hypocentral loci, the theoretical arrival times at all the nine stations are added random perturbations within the range of±0.5 s. The perturbations are generated by a random function of FORTRAN 90 compiler. The arrival time standard deviation is 0.356 s for P waves and 0.321 s for S waves. The resultant hypocentral loci constrained with alone P-arrivals and with both P-and S-arrivals are shown in Fig. 5(k, l). Clearly, they do not cross at the hypocenter but intersect into a small area. Compared with Fig. 5k, the point where the hypocentral loci intersect most densely in Fig. 5l is nearer to the real hypocenter.

4.2 Real Event 4.2.1 Event and velocity model

The real event comes from North China, a significant seismic monitoring area covered with a dense and nearly uniform seismic network. Since hypocentral loci are shown in a two-dimensional profile, an earthquake is desirable when its epicenter is located in a line through several stations. In the east of North China, there is an event of M_L4.2 basically meeting the requirement, which occurred on August 26, 2012 (Beijing time). Reported by Beijing Digital Telemetry Seismic Network (abbreviated to BJSN), the event is at 39.61◦N, 117.42◦E, and 13 km deep; its origin time is 7 h 13 m 34.3 s. Fig. 6 shows the epicentral location and its surrounding stations. It is worth to point out that there are more stations within a bigger area. Obviously, the event has a good azimuthal coverage. Furthermore, its epicenter closely approaches a station (the epicentral distance of station CAD is only 3 km). From the analysis of epicenter accuracy based on seismic network criteria (Bondár et al., 2004), the location result given by BJSN can be concluded to be reliable. The epicentral accuracy may be high up to 2 km or above, which is indicated by a map of BJSN horizontal location errors derived from repeating events (Jiang et al., 2008). In BJSN, earthquakes are located with a simplex method, based on a 3-D traveltime model. The determined traveltimes and arrival-time residuals are displayed in Table 2. The small location residuals in Table 2 also corroborate the conclusion that the source parameters from BJSN are reliable. In addition, the location result by BJSN is highly consistent with that by China Earthquake Networks Center: they are different by less than 1.5 km in epicenter, 3 km in depth and 0.1 s in origin time. Thus, it can be believed that the event’s epicenter stands approximately in the straight line AB through stations EWZ, CAD and XAZ.

The 3-D crustal velocity structure beneath North China has been investigated intensely by seismologists such as Jia et al. (2005), Qi et al. (2006) and Yu et al. (2010). Their investigations reveal that the lateral heterogeneity along the line AB is weak. Besides, one-dimensional models are widely used in modern seismic networks for locating earthquakes. Therefore, we also employ a 1-D model, which consists of three layers. The thicknesses (H) and velocities of P (v_P) and S (v_S) waves of these layers from top to bottom are: H₁=5 km, v_P1=5.0 km·s⁻¹, v_S1=2.6 km·s⁻¹; H₂=15 km, v_P2=6.1 km·s⁻¹, v_S2=3.5 km·s⁻¹; H₃=15 km, v_P3=6.8 km·s⁻¹, and v_S3=3.9 km·s⁻¹. The construction of the model is based on the observed traveltimes and by reference to a model of eastern North China (Teng et al, 1979). We divide the model into 0.1 km×0.1 km square elements. For the hypocenter determined by BJSN, the travel-time residuals with respect to our model are shown in Table 2. These residuals suggest our constructed model is acceptable.

4.2.2 Hypocentral loci

Stations EWZ, CAD and XAZ are combined into pairs. In this way, there are three P-wave arrival time difference constrained hypocentral loci as shown in Fig. 7a. Different loci are drawn in different colors. As can be seen, these loci are approximately similar to hyperbolas, and intersect at one point slightly to the right and below the hypocenter (the red “+”) located by BJSN. Fig. 7b shows the hypocentral loci constrained with arrival time differences between P and S waves at the same stations. The red, green and blue lines respectively correspond to the three stations in left-to-right order. These lines are basically circular in despite of the vertical heterogeneity. Their intersection points spread out. Noticeably, the locus to station EWZ (the red) runs away from the hypocentral location given by BJSN. That is probably related to the large travel time residual of Sg phase at this station. As for the phenomenon that the loci to stations EWZ and XAZ have gaps at the first interface, it results from the fact that seismic waves from the two stations are refracted and travel along the interface, and consequently their travel times change sharply at the refraction interface.

The hypocentral loci constrained with P-wave arrival times are shown in Fig. 7c. The red, green and blue lines respectively correspond to stations EWZ, CAD and XAZ. Among these loci, two are similar to hyperbolas and one near to be circular. The three loci not only intersect at one point but also have a good directional distribution. Especially they impose a tight constraint upon the hypocentral depth. Thus, we infer that the actual hypocenter is somewhat deeper than and to the northeast of that by BJSN. Fig. 7d shows the hypocentral loci constrained with both P-and S-wave arrivals. The red, green and blue lines respectively correspond to stations EWZ, CAD and XAZ; the solid and dashed lines correspond to P and S waves, respectively. Similar to Fig. 7b, the loci in Fig. 7d do not intersect at one point either. However, their intersection area contains the hypocenter located by BJSN. Why hypocentral loci constrained with both P-and S-waves intersect sparsely, different from those with alone P-waves, is perhaps that there are large errors in the picked arrivals or/and the velocities of S waves.

5 DISCUSSIONS 5.1 Stability of Hypocentral Loci

The calculation results of the complex model show that in the view of the directional distribution of hypocentral loci, far stations sometimes seem more adequate for constraining the hypocentral depth (see Fig. 4a and Fig. 5i). That directly contradicts with the common sense. Actually, as for hypocentral loci of an event, how tightly they constrain the hypocenter depends on their stability as well as their directional distribution near the source. When the residual field is regarded as a topographic map, the hypocentral locus in it accordingly lies in a valley. The steeper the valley is, the stabler the locus is; the gentler the valley is, the unstabler the locus is.

Figure 8a presents the arrival time residual field of the hypocentral locus marked with a red circle in Fig. 5a, which corresponds to station R₆. It is noticeable that the locus has two (left and right) segments. Compared with the left part, the low residual (for example, smaller than 0.1 s) zone containing the right segment is rather wide, thereby the latter is decidedly less stable. The blue lines express the calculated hypocentral locus for R_FLRP=0.01 s. By comparison with Fig. 5a, the new calculation result is slightly different and better represents the truth at the part marked with the red circle. When the observations are disturbed by random noise, the right segment with weak stability changes dramatically, as is shown by the part in Fig. 5k marked with the red circle; whereas the right one changes slightly because of its strong stability. This means that unstable hypocentral loci are easily influenced by disturbances.

For the sake of more clearly illustrating the actual role of far stations in constraining the hypocentral depth, we take into account the hypocentral loci constrained with arrival times at stations R₇, R₈ and R₉. Fig. 8b shows the locus for P-wave arrival at station R₉ and its arrival time residual field when the origin time field is constructed with P-wave data alone. Although the locus extends almost in the horizontal direction near the source, it can not tightly constrain the hypocentral depth because its stability is rather weak, especially in the vicinity of the hypocenter; accordingly, it is probably changed significantly in position and shape even by small disturbances. When the S-wave arrival times are additionally used in the construction of the origin time field, the resultant locus corresponding to the same arrival is much stabler, as Fig. 8c shows. This case is similar to those of the loci corresponding to stations R₇ and R₈. By this token, the effect of additionally using S-wave arrival times on the strengthening of the constraint on the hypocentral location is reflected in the improvement not only in the directional distribution of hypocentral loci, but also in their stability.

5.2 Selection of the Calculation Parameters for Hypocentral Loci

The parameter R_FLRP for selecting HLRPs controls the number of primary HLRPs. To ensure the intactness of the locus, the primary HLRPs are required to include the nodes of the model elements where individual endpoints of all the locus segments exist. In other words, the maximum of the residuals at the endpoints is the lower limit of R_FLRP. The lower limit usually increases with increasing element size and decreasing seismic velocity. On the other hand, for we repeatedly “peel off” the skin of the set of primary HLRPs to retain its core as the final HLRPs, R_FLRP must not be so large that the resultant primary HLRPs form connected regions with more than one core. Namely, it is desirable that each connected region contains only one locus segment. Actually, R_FLRP has a quite wide range of values meeting the above requirements. Fig. 9 shows the variation of the percentages of primary (dashed line) and final (solid line) HLRPs in the total model nodes with increasing R_FLRP for the P-wave arrival time constrained hypocentral locus in the homogeneous model. The step of R_FLRP is taken to be 0.005 s. As can be seen, the number of the final HLRPs changes very little when R_FLRP varies within a fairly wide range from about 0.03 s to 0.5 s. This stability suggests that the parameter R_FLRP does not need adjusting with the geometry of hypocentral loci, which is beneficial to the automated processing of large batches of data. The value of R_FLRP, under the condition that the calculated locus is complete, is suggested to be as small as possible. This not only saves computational time but also ensures good representativeness of the selected HLRPs for weakly stable hypocentral loci.

The parameter R_FL associated with the calculation of ray paths controls the ray tracing area. To ensure the intactness of the hypocentral locus and prevent false hypocentral loci, the connected regions with residuals less than R_FL are required that they each completely contain one and only one locus segment. Namely, R_FL must be no less than the maximum of residuals at the model nodes crossed by the hypocentral locus but must not be so large that there is one or more connected region with more than one locus segment. When the hypocentral locus is one-segmented, the upper limit of R_FL can be large enough that the ray tracing area extends to the whole model space. As for a multi-segment hypocentral locus, the value range of R_FL is fairly wide. Taking Fig. 8a as an example, it is proper to choose R_FL within the range of 0.1~0.5 s. For a small value of R_FL, the ray tracing area is one or more narrow zones of which the percentage in the model region is very low, and accordingly the calculation of the hypocentral locus is fast.

The modification parameter L_twig of hypocentral loci controls the length of ray paths to be removed. Its value is related to the selected HLRPs. When HLRPs are distributed along the hypocentral locus, resulting in short ray paths without branching, the value of L_twig can be chosen relatively small. When the hypocentral locus is weakly stable and the parameter R_FLRP is rather large, the final HLRPs probably contain poorly representative points. In this case, L_twig should be chosen a large value so as to remove the rays to the points away from the hypocentral locus.

5.3 Extension of the Algorithm to Three Dimensions

Earthquake location is a three-dimensional problem. In a three-dimensional residual field, nodes of the model elements crossed by the hypocentral locus are local minimum points. These nodes lie in the middle of the low residual regions, and are more similar to one or more cores of a fruit. Therefore, the peeling algorithm for selecting HLRPs is adapted to three dimensions except for slight difference in the elimination of superfluous HLRPs. As for ray paths, they are calculated with the minimum traveltime tree algorithm which has the same frame for two-and three-dimensional models. In the case of three dimensions, the calculated ray paths also form tree structures, and hence they can be modified to be finer in the similar way to that of trimming branches. Thus, the improved method presented in this study can be easily extended to three-dimensional models by means of 3-D minimum traveltime tree algorithms (Klimeš et al., 1994; Wang et al., 2000; Zhao et al., 2004). Compared with two-dimensional models, the calculation of hypocentral loci in 3-D models costs much more computational time, mainly because it is time consuming to model the three-dimensional travel time fields that shot points are located at stations. However, the time consuming step is needed only once. Therefore, the method is quite efficient in calculating hypocentral loci of numerous events within a region that has a stable velocity structure, recorded by the same seismic network.

6 CONCLUSIONS

The improved method put forward in this article for calculating hypocentral loci is as robust as the original one and hence suitable for complex models. Moreover, it works well in maintaining the completeness and the fineness of calculated hypocentral loci. There are no restrictions on the number of segments that constitute a hypocentral locus. The ray tracing is confined to low residual regions, and that saves a large amount of calculation time. The parameters for calculating hypocentral loci are easy to choose and do not need to be adjusted frequently with geometry of the loci. This is favorable for the automated processing of large batches of data. Additionally, the improved method can be easily extended to three dimensions. Its calculation efficiency is quite high for more than one event from a region with a stable velocity structure, recorded by the same seismic stations because the most time consuming travel time calculations have to be done only once.

The numerical tests indicate that hypocentral loci constrained with arrival time differences are different in form from those with arrival times, and the latter are quite complex even for homogeneous media. The capability of focal loci to constrain the hypocenter depends on their directional distribution near the source and their stability. Arrival times at near stations are strong in constraining the hypocentral depth. When both the velocity model and the picked arrival times are fairly accurate, the additional use of S-wave data can significantly improve the constraint upon the hypocentral location, especially for events outside the seismic network. Using more hypocentral loci can decrease the influence of random factors on earthquake location.

ACKNOWLEDGMENTS

This work is supported by the National Natural Science Foundation of China (41374098 and 40974050) and the Special Foundation for Basic Professional Scientific Research (DQJB11C03). We thank Professor Zhongjie Zhang, Institute of Geology and Geophysics, Chinese Academy of Sciences, who provided support and advice. Our thanks are also due to Dr. Baofeng Tian and Dr. Lihua Fang, Institute of Geophysics, China Earthquake Administration, for their help in obtaining the real event data from Beijing Digital Telemetry Seismic Network.