Physics of the Earth and Planetary Interiors 113 Ž1999. 355–365
A two-point, three-dimensional seismic ray tracing using genetic algorithms Hossein Sadeghi ) , Sadaomi Suzuki 1, Hiroshi Takenaka
2
Department of Earth and Planetary Sciences, Kyushu UniÕersity, Hakozaki 6-10-1, Fukuoka 812-8581, Japan Received 5 January 1998; received in revised form 4 May 1998; accepted 4 May 1998
Abstract Ray tracing between two fixed points is a boundary value problem. An initial path connecting two points, the source and the receiver, is algorithmically perturbed until it converges to a solution giving a minimum travel-time path. In multi-pathing cases, it is quite possible for an algorithm to converge to a ray path with local minimum travel time. Furthermore, when the velocity structure has discontinuities, there is another serious problem, i.e., how to determine the correct intersection between the ray and the surface of discontinuity. In this paper, a new approach for two-point ray tracing is presented, which uses genetic algorithms ŽGAs. to overcome these problems. Since GAs are efficient global search techniques, the proposed method guarantees to find a global minimum travel-time path and also the correct intersections. A micro-GA implementation is adopted to further enhance computational efficiency. This approach is suitable for tracing the first arriving seismic waves through a 3-D heterogeneous medium with discontinuities. The method can also find the minimum time paths of later arrivals, such as reflected and refracted waves, by constraining the ray to hit a prescribed interface. The accuracy of the method is demonstrated on the numerical examples. q 1999 Elsevier Science B.V. All rights reserved. Keywords: Discontinuity; Genetic algorithms; Minimum travel time; Ray tracing
1. Introduction Seismic ray tracing is a powerful tool for studying complicated heterogeneous structures. There are three methods for tracing rays through a velocity structure: shooting, exhaustive grid search and bending. In the shooting method, ray tracing is an initial value problem. An initial point is fixed and the ray is
) Corresponding author. Fax: q81-92-642-2685; e-mail:
[email protected] 1 E-mail:
[email protected]. 2 E-mail:
[email protected].
propagated by using some first guess of take-off angles. The take-off angles are iteratively improved until the ray hits the target location. The target point is frequently an ill-behaved function of the take-off angles. It is then impossible to trace some rays ŽLangan et al., 1985.. The exhaustive grid search tries to find the first arrival rays by finite-difference solution of the eikonal equation Že.g., Vidale, 1990., or by using network theory Že.g., Moser, 1991.. The implementation requires computation and storage of travel times to all spatial locations in a 3-D grid. It can be efficient if travel paths for a large set of sources and receivers are needed.
0031-9201r99r$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S 0 0 3 1 - 9 2 0 1 Ž 9 9 . 0 0 0 1 1 - 4
356
H. Sadeghi et al.r Physics of the Earth and Planetary Interiors 113 (1999) 355–365
In the bending method, ray tracing solves a boundary value problem by iteratively refining a proposed path between fixed source and receiver positions. The ray path is refined by searching a minimum time path based on Fermat’s principle. Exact bending ray-tracing methods Že.g., Wesson, 1971; Julian and Gubbins, 1977; Pereyra et al., 1980. are very time-consuming and formidable for routine use in seismological applications ŽThurber and Ellsworth, 1980.. Um and Thurber Ž1987. introduced an accurate approximate algorithm, the pseudo-bending technique, as a fast alternative to the bending method. The pseudo-bending is also more reliable than the exact bending with respect to numerical stability ŽKoketsu and Sekine, 1998.. However, this algorithm has still two inherent limitations. For complicated velocity structures, there can be many paths connecting two points of interest. When multi-arriving is possible, the solution depends on the starting path, and it is easy to miss the first arrival. Updating the ray path around discontinuities is another problem because in general, bending requires a gradient direction to update the path. We propose a new ray bending approach for tracing the first arrival rays through a heterogeneous structure with discontinuities. A kind of genetic algorithm ŽGA. is used as a global search technique for minimum time ray tracing. While large-population GAs have already been used in several other seismological field Že.g., Sen and Stoffa, 1992; Nolte and Frazer, 1994; Xie et al., 1996., in this paper a micro-population GA is employed to increase the computing speed. This GA, termed a micro-genetic algorithm Žmicro-GA., was first proposed by Krishnakumar Ž1989.. In our new approach for ray tracing, the micro-GA introduces a starting ray path that is set to be in the vertical plane containing two endpoints. The initial path is then refined three-dimensionally in continuous regions by the pseudobending, and in discontinuities by the micro-GA. We call this method the micro-GA bending. The method also allows the calculation of later arrivals, such as reflected, refracted and converted waves. In Section 2, we give a brief review of the micro-GA, and in Section 3, we describe our approach to ray tracing in detail. Then, in Section 4, we discuss its accuracy and efficiency. Conclusions are mentioned in Section 5.
2. Genetic algorithms GAs are stochastic search algorithms based on the mechanics of natural selection and natural genetics ŽGoldberg, 1989.. Unlike other stochastic techniques, GAs work simultaneously with a population of solutions in the search space. A possible solution to a problem is referred to as an indiÕidual. The population sizes are quite often in excess of 50 individuals ŽJohnson and Abushagur, 1995.. Each individual is coded as a string Ž chromosome., where variables are sub-strings sitting next to each other on the chromosome. In the most general form, a chromosome consists of a string of binary bits. Bits of the string are called genes and the two distinct values of a gene, 0 and 1, are called alleles. The number of bits in a chromosome depends on the required precision of the associated variables. If N is the number of bits assigned to a variable, then there are 2 N possible values which will be evenly spread out in the search internal of that particular variable. If our problem has M variables to each of which N bits are assigned, the length of chromosome will be M = N. In each iteration, called a generation, the randomly selected initial population is updated so that it approaches the optimum. A single iteration of a typical GA proceeds in the following three stages: reproduction, crossoÕer and mutation ŽGoldberg, 1989.. We briefly describe these stages below. For comprehensive descriptions, the reader can refer to Sambridge and Drijkoningen Ž1992.. Reproduction is a process in which individual strings are copied according to the values of their objective function Žcalled the fitness .. This means that strings with a higher fitness have a higher probability of being reproduced and passed down into the next generation. The probability of reproduction is calculated for each individual by: fi Pi s Q , Ž 1.
Ý
fk
ks1
where f i is the fitness of the individual i, and Q is the population size. This probability is used as the basis of the selection of a new population from the old. For selection we use a stochastic tournament algorithm Že.g., Goldberg, 1989.. In this algorithm,
H. Sadeghi et al.r Physics of the Earth and Planetary Interiors 113 (1999) 355–365
successive pairs of individuals are drawn at random from the entire old population, but the probability of choice is proportional to the individual’s probability Pi . A random number, r, is generated between 0 and 1, and compared to Pi . If Pi ) r, then the individual is drawn. After drawing a pair, the string with higher fitness is inserted in the new population and another pair drawn. This process is continued until the population has the same size as the old one. After reproduction, the individuals in the population are randomly paired. The crossover combines each pair with probability Pc . If the crossover is performed, two offsprings replace their two parents. Otherwise, both parents follow to the next generation unchanged. There are different types of crossover. The simplest crossover is single-point crossover: a position is chosen at random along the bit-strings, and the portions after this position are swapped. Other types are multiple-point Že.g., Goldberg, 1989. and uniform crossover ŽSyswerda, 1989.. Multiplepoint crossover selects several positions at random and swaps every other portion between these positions. Uniform crossover randomly exchanges bits of each pair with probability Pc s 1r2. Fig. 1 illustrates single-point, two-point, and uniform crossover. At the end of a single iteration, mutation is applied to each gene of each offspring with very small
357
probability Pm . The value of the mutated gene is changed from 0 to 1 or vice versa. As a contrast to the large population GA mentioned above, several authors have introduced a micro-GA ŽKrishnakumar, 1989; Johnson and Abushagur, 1995.. The micro-GA relies on the same principals as the conventional approach; however, it uses a population of only five estimates with no mutation. Although small populations converge very fast to local optima, restarting provides the diversity required to achieve the global optimum. With large populations, the diversity is provided from the beginning, but the processing of the genetic operations is much slower. The following pseudocode shows the micro-GA structure we use: Create an initial population of five chromosome at random Do for i s 1 to MAX_ITERATIONS until Žstopping criterion. while not Žconverged. evaluate the fitness of each chromosome save the best one select parent chromosomes for reproduction crossover chromosome4 save the best estimate replace remaining four members at random 4.
Fig. 1. Crossover comparison for single-point, two-point and uniform operation. Vertical bars denote crossover points.
358
H. Sadeghi et al.r Physics of the Earth and Planetary Interiors 113 (1999) 355–365
An initial set of five individuals is selected at random. Each individual is then evaluated in term of its fitness, and the best estimate is saved. The remaining four individuals compete for the reproduction by the tournament selection method. The selected parents are mated and crossed-over by means of the uniform crossover operation. This cycle is repeated until convergence is realized. The convergence condition is that the number of different bits from the best member in the micro-population is less than 5%. Once this condition is achieved, the best individual is passed onto the next iteration and an additional four individuals are selected at random. When the best estimate of a convergence does not differ significantly from the previous one, we stop the micro-GA search.
3. Micro-GA and two-point ray tracing Our algorithm Žthe micro-GA bending. in the first step uses the micro-GA for determining a good estimate of the starting ray path. The procedure may be illustrated easily with an example ŽFig. 2.. Suppose the micro-GA wants to introduce a starting path from the source at Xs to the receiver at X r . Let a small number of equidistant points Že.g., three points Xa , X b , and X c . discretize the straight path connecting the source and the receiver. The positions of these points Ždepth, h; latitude, u ; and longitude, f . are used to define the path. If we suppose these three points are in the vertical plane containing the source and receiver, then the depth h is the only variable. After defining the search areas for the depth of each point, they are coded as binary numbers into a chromosome. The fitness of the chromosome is eval-
Fig. 3. Illustration of the perturbation scheme which uses the micro-GA to find the intersection point between a ray path and a discontinuity.
uated by the travel time of the associated ray path. The travel time is calculated in three dimensions by numerical integration along the estimated ray path. The path is divided into n short segments and the travel time integral is approximated by: n
Ts
Ý Ui D X i ,
Ž 2.
is1
where Ui is the average slowness of the two endpoints of the ith segment and D X i is the length of the ith segment. The goal is to find a set of depth values or chromosomes, which minimize the travel time. This minimization is performed through the micro-GA; however, each cycle of convergence through the micro-GA may yield a local minimum. In order to ensure the global minimum, we continue the iterations until for two successive convergence, the travel time T no longer decreases, or the relative decrease DTrT is less than a preset tolerance limit. After completing the micro-GA search, the best estimate is
Fig. 2. Each cycle through the micro-GA converges a population of randomly selected paths to a minimum time path.
H. Sadeghi et al.r Physics of the Earth and Planetary Interiors 113 (1999) 355–365
359
Fig. 4. The initial path guessed by the micro-GA is iteratively perturbed by the joint use of micro-GA and pseudo-bending technique.
selected as the starting path for the pseudo-bending technique. With a reasonable starting path, the pseudo-bending technique works well, but still fails to work in cases where the velocity structure has discontinuities. When velocity discontinuities exist, internal boundary conditions have to be perturbed. We then again use the micro-GA for this perturbation. Consider a case, as shown in Fig. 3, in which two media with continuous velocity structures are separated by an irregular-shaped boundary Ždiscontinuity surface.. A and B are two points along a seismic ray, which are located on different sides of the discontinuity. C is the intersection of the seismic ray from A to B and the discontinuity. We use the micro-GA to find the location of the point C. We then consider a point CX on the straight line AB so that its vertical projection onto the discontinuity is C. Actually, the micro-GA searches to find the point CX instead of C. The search interval of A to B is coded into a chromosome. To evaluate the fitness of a chromosome, at first its real value is vertically projected onto the discontinuity, and then the travel time is calculated for the estimated ray path AB through the projected point as: TAB s TAC q TBC , TAC s
