Introduction
Typical problems associated with locating regional earthquakes are inaccurate arrival times, an inaccurate velocity model and therefore a less stable G matrix. Because epicentral distances are usually much larger than the distances among stations, unknown parameters in the G matrix (hypocenter and origin time of the event) can only be solved for epicenters and origin times with fixed focal depths (e.g., Palvis, 1992, Lienert, 1997). An additional problem arises when using a small aperture array to locate regional earthquakes because parameters of the epicenters become highly dependent on the origin time parameter in the G matrix. Common location programs based on the generalized inverse method like HYPO71 (Lee and Lahr, 1972) or SEISMOS (Hartse, 1991), which were originally designed for locating local earthquakes, tend to fail because no clear minima exists.
In this study, we approached the problem by using a modified G matrix containing only travel time intervals. Thus, no origin time parameter was included in the modified G matrix. Furthermore, we solved the modified G matrix using a forward method instead of an inverse method. During the process, the G matrix was regrouped into two matrices based on the characteristics of travel time differences: the matrix confining the epicentral distance (S-P time interval) and the matrix confining the distribution of azimuth (P-P and S-S time intervals). Individual matrices were mapped into the fuzzy logic space with measured uncertainties in terms of deviations in wave velocities. Locations of epicenters of earthquakes were then evaluated in the logic space.
A real earthquake was selected as a case study in this paper. The location of the epicenter of the event was evaluated based on both inverse and forward methods. In addition, a series of tests on the sensitivity of the location methods to faulty arrival times were conducted to evaluate the fitness of both methods.

Common Procedure for Earthquake Location
The generalized linear inversion method is the most commonly used procedure to achieve a high accuracy location of hypocenters for earthquakes. There are four hypocentral parameters (x, y, z, t), the spatial coordinates and the origin time to be solved by this method, and this is usually done by using Geiger's method or its variants.
Typically Geiger's method starts with the time residual between the arrival time at the ith station and the trial solution

(1)

By assuming the time residual is small, a Taylor expansion of it will give

(2)

where m is a known model vector composed of the source location and origin time

. (3)

Let G be equal to the partial derivative matrix

(4)

and the equation becomes

. (5)

Since this is often an overdetermined problem it can be solved with the generalized inverse of G or
G^-g so that

. (6)

As a result, this method is iterative and the new model

(7)

is repeatedly adjusted until the minimum misfit is reached.
This procedure works well for local earthquakes occurring within networks. However, for regional earthquakes occurring outside a network and hypocenter-station distances significantly larger than distances between stations, it is more difficult to constrain focal depths. People usually circumvent this problem by using fixed focal depths. Nonetheless, problems still arise when closely spaced stations have little azimuthal constraint on the epicenters of earthquakes. As a result, the epicenters are highly dependent on observed arrival times and initial epicenter estimates. Small errors in observed arrival times often result in huge shifts to epicenters and sometimes the epicenter is placed at a local minima far from the real location.

Modified G Matrix
Another approach to solving this problem is to remove the timing parameter from the equation of time residual by using only S and P travel time intervals. With little modification of the travel time-distance equation the time residual becomes

. (8)

The major disadvantage of this approach is that we will be unable to incorporate direct P and S arrivals into the G matrix, which are crucial losses to an already limited number of recorded seismic parameters.
It is possible to increase input data to equation (8) by expanding the partial derivative matrix G to accommodate both P-P and S-S travel time intervals. During the process, we assume a one dimensional crustal model and small azimuthal distributions of seismic stations. In other words, all P-P and S-S travel time intervals are required to have similar ray paths in order to be included into the G matrix. Near field variations of crustal structure can be taken into account by using station corrections. As a result, we retain the same forward model as shown in equation (8). However, for P-P and S-S travel time intervals become travel time differences between two stations through a specified layer and S-P travel time intervals are travel time differences due to differences in P and S wave velocities along the same raypath to a single station.
Figure 1 illustrates a simple example of travel time interval between two stations with respect to a trial epicenter. As shown in the figure, the distances from a trial epicenter to stations and are and , respectively. Both and are significantly larger than the spacing between the two stations, , and have similar raypaths. After station corrections for near field variations in crustal structure, the difference between arrival times can be attributed to the travel time interval through a specifed layer with wave velocity V. For a half-space velocity model the travel time interval between two stations can be simply expressed as

(9)

where

and . (10)

Note that the origin time term does not appear the equation.
Unlike relationships between travel distances and origin time derived from direct P and direct S arrival times, P-P and S-S travel time intervals result in relationships between D and as shown in Figure 2. The curve radiates outward approximately from the mid-point between two coupling stations with asymptotes of

. (11)

As distance D increases, changes in azimuthal angle with change in location decreases to minimum. Thus, is less sensitive to D when D is significantly larger than d. Also, the unconfined nature of curves makes the G matrix become less stable during the inversion process. Therefore it would not be practical to incorporate P-P and S-S travel time intervals into the G matrix and solve for earthquake locations using an inverse method.

Forward Modelling Plus Fuzzy Logic Approach
Distinctive Features
In this study we solve equation (8) using a forward modelling technique. Instead of solving the modified G matrix directly as in an inverse method, we modify the procedure for optimum location of regional earthquakes. Distinctive features of the modified procedure are listed below:
Quantifying uncertainties. Recorded arrival times of regional earthquakes often contain large uncertainties in both arrival times and in the velocity model. Traditional methods usually assume a perfect velocity model and attribute uncertainties of recorded phases to uncertainties in reading arrival times on seismic records. By excluding variations in velocity structure along their raypaths and assigning fixed numbers of uncertainties to phases with respect to the arrival time, earthquake locations are often incorrect. In this study we convert uncertainties in arrival times into uncertainties in the velocity model. With preset upper and lower bound velocities for individual phases, the distance weighing factor is easily taken into account in terms of uncertainties in wave velocities.
Regrouping G matrix. Derived curves from S-P, P-P, and S-S travel time intervals are divided into two groups based on their distinct nature. In other words, circular curves of S-P travel time intervals constrain epicentral distance and parabola like curves of P-P and S-S travel time intervals constrain azimuth. Joint evaluation of earthquake locations by combining results from a S-P matrix and a P-P and S-S matrix provide more stable results.
Logic operations using fuzzy logic. Fuzzy logic operation is primarily a forward modeling technique. This method is especially useful for handling data with large uncertainties. Deviations between trial and observed travel times in equation (8) for different time intervals can be mapped into the logic space. In the logic space, logic operations can be easily carried out according to relationships between individual G matrices and yield optimal results for earthquake location.
Fuzzy Logic Approach
The theory of fuzzy logic deals with two problems: fuzzy set theory, which deals with the ambiguity found in semantics, and fuzzy measurement theory, which deals with the ambiguous nature of judgements and evaluations.
A fuzzy set may be represented by a mathematical formulation often known as the membership function. This function gives a degree or grade of membership within the set. The membership function of a fuzzy set A, denoted by , maps the elements of the universe X into a numerical value within the range [0, 1], i.e.,

(12)

(13)

A simple graphical comparison of classic (or crisp) set theory and fuzzy set theory is shown in Figure 3 (Jamshidi et al., 1993).
By adapting fuzzy logic theory, we are able to incorporate uncertainties from measurements of arrival times into variations of velocity structure. Figure 4 shows an example of such assessment. Let equal the difference in epicentral distances with respect to station 1 and station 2 (Figure 1). Thus the uncertainty of travel time difference at these two stations can be expressed as the uncertainty in the velocity model

(14)

and

. (15)

In Figure 4a, the region bounded by theoretical wave velocity V plus or minus one standard deviation yields a fuzzy output of 1. Surrounding this region is the area bounded by the highest and the lowest possible wave velocity. Note that in these two areas, the fuzzy logic outputs are less than one, which represents possible but less likely results.
Joint determination of earthquake location from equation (8) for individual time intervals in logic space is a straight forward logic operation. Fuzzy logic operations are identical to classic logic operations such as union, intersection, complement, and difference. In this study we use S-P, P-P, and S-S time intervals which can be divided into two groups, an S-P group and a P-P and S-S group. Evaluation of a trial location in logic space can therefore be expressed as

. (16)

Earthquake Location Procedure for New Mexico
Problem Definition
In New Mexico earthquakes occur throughout the state, but recording is frequently by only one small aperture network (see Figure 6). Thus, many events are at large distances relative to the dimensions of the networks. SEISMOS, the currently used location program at New Mexico Tech (Hartse, 1991), uses the location of the closest seismic station as the initial estimate of epicenter. The assumption is valid for events occurring within or near a seismic network. However, at event distance several times the network aperture, the assumption seldom applies, and the corresponding initial estimates of hypocenters provide little to no help in the location process.
Another problem with locating regional earthquakes comes from variability in crustal structure (Stewart and Pakiser, 1962, Toppozada 1974, Sinno et al., 1986). Unlike local events, where we are able to solve for hypocenters using a detailed crustal model, the New Mexico area has large differences in crustal structure, for example Moho depths ranging from about 30 km to more than 50 km. Based on the available data, a crustal model applicable throughout the state is simply not possible. Sanford et al. (1991) circumvented the problem by ignoring Pn arrivals and by using only Pg and Sg arrivals with a half-space crustal model with P wave velocity of 6.15 km/s and a Poisson's ratio of 0.25. This procedure gives the best locations for earthquakes occurring within New Mexico and bordering areas.
Defining Fuzzy Sets
In this demonstration study, we fix focal depth for regional earthquakes and solve for the epicenter (x, y). We use a homogeneous half-space model with P wave velocity of 6.15 ± 0.05 (1 s.d.) km/s and Poisson's ratio of 0.25 ± 0.01 (1 s.d.). Four fuzzy sets are defined:
1. P-P travel time interval (PP). Theoretical and observed P wave travel time intervals at two arbitrary stations for a trial epicenter are compared using equation (14). Deviation of P wave velocity with respect to the velocity model is then mapped into logic space as shown in Figure 5a. According to the model, deviations within 6.15 ± 0.05 km/s yield a fuzzy output of 1. Fuzzy outputs for larger deviations decrease linearly until the maximum possible velocity of 6.4 km/s or minimum possible velocity of 5.9 km/s are reached, which both have fuzzy outputs of 0.
2. S-S travel time interval (SS). The S-S time interval fuzzy set is similar to that for P-P time intervals. Figure 5b shows the fuzzy model for S wave velocity. The theoretical S wave velocity of 3.551 km/s is based on a P wave velocity of 6.15 km/s and a Poisson's ratio of 0.25. We assume variations of Poisson's ratio ranging from 0.24 to 0.26. Thus, the maximum possible S wave velocity is 3.743 km/s (P wave velocity 6.4 km/s and Poisson's ratio 0.24) and minimum possible S wave velocity is 3.360 km/s (P wave velocity 5.9 km/s and Poisson's ratio 0.26). Deviations of S wave velocity 3.551 ± 0.03 km/s corresponding to a deviation between trial and observed time interval in equation (14) will yield a fuzzy output of 1.
3. S-P travel time interval (SP). S-P time intervals only apply to stations with both P and S readings. For a half-space model, S-P travel time interval can be expressed as

(17)

, (18)

where

. (19)

Maximum and minimum factor C were obtained by mixing P and S velocity models. Finding distance is the objective of this set yet there is no direct way to calculate deviations of the velocity model with respect to travel time differences. To simplify the process, we assigned a fixed width of uncertainty as shown in Figure 5c. As usual, differences between theoretical and observed hypocentral distances for a trial hypocenter will be examined and mapped into the logic space.
In summary, assume that a seismic network records n P phases, m S phases, and k P and S phases for an earthquake. The fuzzy output of a trial epicenter can be expressed as

. (20)

Example Earthquake
We chose one of the most recent felt earthquakes in New Mexico as an example to illustrate the procedure for locating earthquakes. This duration magnitude 3.0 earthquake occurred on 14 July, 1998 at 05:38UT near Logan, New Mexico. Figure 6 shows the geographical locations of the earthquake and the seismic network that recorded it. The nearest station of the seismic network to the epicenter (CPRX) is ~260 km and the farthest (GDL2) is ~360 km. The aperture of the network with respect to the earthquake is ~15^o. The epicenter is reasonably well constrained on the basis of felt reports.
Figure 7 shows the P-P and the S-S travel distance curves for stations CPRX and CL7. The two curves in the figure are theoretical relationships between epicentral D and azimuth based on comparisons of P and S travel time differences, respectively. A narrow band of fuzzy output 1 (shown in black) indicates maximum tolerance of uncertainty. Surrounding the black band are two gray bands in descending color, indicating decrease in fuzzy outputs. Logic operations applied on these outputs depend on the purpose of applications. Figure 8 shows two of the most commonly used operations: Union and intersection. It is clear that the outputs of union operation retains more information but converges slower. On the other hand, an intersection operation results in faster convergence but tends to lose more information. A comparison of the intersection operation using classic and fuzzy logic procedures is shown in Figure 9.
Figure 10 shows results of P-P, S-S, and S-P travel distance curves, respectively. The study area was divided into a matrix of 100 by 100 cells of 10 km on each side, and the center of each cell was used as a trial epicenter. For each travel time interval, the fuzzy outputs of trial epicenters are the results of union operations for all combinations of travel time intervals. As we expected, results of P-P (Figure 10a), and S-S (Figure 10b) travel distance curves reveal a prominent NNE trending azimuth. Circular travel distance curves of S-P time intervals (Figure 10c) constrain epicentral distance.
The matrix for deriving the final location of epicenter was reached by combining all three travel distance matrices based on equation (20). In Figure 11, unwanted or less likely trial epicenters were removed from the output matrix after logic operations. Therefore, the most likelihood location for the earthquake can then be easily resolved using a center of gravity method. For this example, the location derived by this method is almost identical to the location obtained from the traditional method with the addition of arrival times from three stations in another network.
In our example earthquake, we replaced some of the observed arrival times with erroneous arrive times to simulate typical errors caused by human operations or automated computer operations. Because uncertainties in the travel time differences in the G matrix are not redistributed throughout the matrix during the forward modeling process, the erroneous data are simply dismissed during the evaluation procedure. Figure 12 shows the location of the estimated epicenter for the Logan event after we included large and numerous erroneous data entries with respect to results from the actual observed data. It is clear that the forward process tolerated several erroneous arrival times in its G matrix and yielded approximately the same location as with the actual observed data. In contrast, a typical inverse method fails to solve the location of the epicenter as soon as a single faulty entry is included in the G matrix.

Discussions and Conclusions
The most unique feature of the modified G matrix we used in this study is that it contains only arrival time differences between phases. In solving the data matrix, we divided the matrix into two matrices for deriving hypocentral distance and azimuthal angles with respect to seismic stations. Uncertainties in arrival time differences were mapped into the fuzzy logic space in terms of uncertainties in velocity model. Locations of earthquakes were determined based on results of logic outputs in the logic space by searching a gridded area.
The forward modeling method used in this study has advantages over inverse method in handling uncertainties of arrival times and even man made errors. As shown in equation (6), in the process of deriving a generalized inverse matrix , errors and uncertainties are redistributed throughout the matrix. A wrong entry of arrival time is likely to destroy the whole matrix. On the other hand, errors of arrival times are isolated in the forward method and have less effect on other data entries. In our example earthquake, we were able to introduce several P and/or S arrival time errors out of 13 arrival time set without destroying the result. This technique can be adapted to applications such as automated triggering system which sometimes must contend with misidentified phases.
We have incorporated the fuzzy logic algorithm into the location program SEISMOS to increase stability in locating regional earthquakes. This technique was converted into a computer subroutine as an initial hypocenter estimator and incorporated into SEISMOS in early 1994 at New Mexico Tech (Lin and Sanford, 1994). The Fuzzy/SEISMOS combination has proved to be very effective in locating regional earthquakes.

References
Hartse, H.E., Ph.D. Thesis (1991). Simultaneous Hypocenter and Velocity Model Estimation Using Direct and Reflected Phases from Microearthquakes Recorded Within the Central Rio Grande Rift, New Mexico, New Mexico Institute of Mining and Technology, New Mexico.
Jamshidi, M., N. Vadiee , and T.J. Ross, Editors (1993). Fuzzy Logic and Control Software and Hardware Applications, PTR Prentice-Hall Inc., New Jersey.
Lee, W.H.K. And J.C. Lahr (1972). HYPO71: a computer program for determining hypocenter, magnitude, and first motion pattern of local earthquakes, U.S. Geol. Surv. Open-File Rept. 75-311.
Lienert, B.R. (1997). Assessment of earthquake location accuracy and confidence region estimates using known nuclear tests, Bull. Seism. Soc. Am., 87, 1150-1157.
Lin, K.W., A.R. Sanford (1994). Regional earthquake hypocenter location using a fuzzy logic algorithm enhanced seismos program, New Mexico Institute of Mining and Technology Geophysics Open File Report 74, 82p.
Palvis, G. (1992). Apprasing relative earthquake location errors, Bull. Seism. Soc. Am., 82, 836 859.
Sanford, A.R., L.H. Jaksha and D.J. Cash, Editors (1991). Neotectonics of North America: Boulder, Colorado, Geological Society of America, 229-244.
Sanford, A.R., K.W. Lin, and I.C. Tsai (1997). Preliminary seismicity map for New Mexico and bordering areas, New Mexico Institute of Mining and Technology Geophysics Open File Report 83, 2p.
Sinno, Y.A., P.H. Daggett, G.R. Keller, P. Morgan, and S.H. Harder (1986). Curstal structure of the southern Rio Grande rift determined from seismic refraction profiling,, J. Geophys. Res., 91, 6143-6156.
Toppozada, T.R. Ph.D Thesis (1974). Seismic Investigation of Crustal Structure and Upper Mantle Velocity in the State of New Mexico and Vicinity, New Mexico Institute of Mining and Technology, New Mexico.

Figure 1. Schematic diagram for seismic stations S₁ and S₂ and the trial epicenter O. Epicentral distance D represents the distance from O to the midpoint between two seismic stations.

Figure 2. Travel distance curves for S-P, P-P, and S-S travel time intervals. Note that S-P travel time interval resulted in circular curve and P-P and S-S travel time intervals in parabola like curves.

Figure 3. Comparisons of classic logic and fuzzy logic. For a given logic set A, complement, intersection and union operations based on classic logic and fuzzy logic are shown.

Figure 4. Results of P-P travel time interval between station CPRX and CL7. (a) P wave velocity model for mapping uncertainties in calculated wave velocities into logic space. (b) Derived travel distance curve. Theoretical travel distance curves for the velocity model are shown in red and pink lines, respectively.

Figure 5. Velocity models for mapping uncertainties in calculated wave velocities into logic space. (a) P wave velocity. (b) S wave velocity. (c) S-P travel time coefficient.

Figure 6. Geographical location of the Logan earthquake with respect to seismic networks in New Mexico.

Figure 7. Results of (a) P-P and (b) S-S travel time intervals between station CPRX and CL7 after mapping uncertainties in calculated velocities into logic space.

Figure 8. Results of fuzzy logic operation (a) Union and (b) Intersection for P-P and S-S travel time intervals between station CPRX and CL7.

Figure 9. Results of (a) classic and (b) fuzzy logic intersection operations for P-P and S-S travel time intervals between station CPRX and CL7.

Figure 10a. Results of P-P travel time intervals for all recorded stations (6 stations). A total number of 120 travel distance curves are plotted. The fuzzy logic operation union was applied to individual travel distance curves.

Figure 10b. Results of S-S travel time intervals for all recorded stations (7 stations). A total number of 720 travel distance curves are plotted. The fuzzy logic operation union was applied to individual travel distance curves.

Figure 10c. Results of S-P travel time intervals for all recorded stations (6 stations). A total number of 6 travel distance curves are plotted. The fuzzy logic operation union was applied to individual travel distance curves.

Figure 11. Resolution matrix for determining epicenter. This matrix was derived by combining the three individual matrices in Figure 10 using equation (19). The best location of epicenter was obtained using a center of gravity method.

Figure 12. Resolution matrices with various levels of arrival time errors: (a) One P phase reading offset by 10 seconds; (b) One P and one S phase reading offset by 10 seconds; (c) Two P and one S phase readings offset by10 seconds; and (d) Two P and two S phase readings offset by 10 seconds.