An Improved Vicinity Ray Tracing in Tuo-Dimensional Laterally Inhomogei:ieous Velocity Structure

We have developed a modified vicinity ray tracing (VRT) method that has improved the travel time approximations and strengthened the theoretical base of the original one. The vicinity ray tracing (Kim and Cormier, 1990) method is a high-frequency asymptotic procedure to compute the wavefield in laterally inhomogeneous velocity structures. However, in calculating travel times the original VRT method only considers the central rays on the far side of the station but misses those on the side near the source. The modified version extends the. travel time approximation to the rays on both sides of the station. This is facilitated by assuming that the ray paths near the surface in the vicinity of the station are parallel to each other. The modified approach has been examined with a variety of velocity models. For the model with a low velocity layer, the synthetic seismograms demonstrate that the dependence of the amplitudes on the source-receiver distance in the shadow zone is as clear as that illustrated using a classical Gaussian Beam Method. When strong lateral heterogeneity and thus multiple caustic regions are present, the synthetic wavefield is still shown to be stable. As another improvement, the calculation of wavefront parameters is furnished in kinematic ray tracing system without invoking Kim and Cormier's dynamic VRT system. This alternative procedure can save large computational time without sacrificing the accuracy in synthetic seismograms.


INTRODUCTION
There are primarily two major trends for solving the wave equation.The numerical approach includes typically finite difference, finite element, and pseudo-spectrum methods.They provide a full solution to the wave equation, but the longer calc ulation time and more required computer space are two major drawbacks of them.The Gaus sian beam method (GBM) and WKBJ/Maslov method belong to another category that usually employ high frequency asy mptotic approximation and can be applied in laterally inhomogeneous structure.
Both the GBM and WKBJ/Maslov are based on asymptotic ray tracing (ART) but avoid the major defect of the ART; the ART cannot deal with caustic regions, such as shadow zone, triplication zone and critical regions.There are intrinsically two parts in the algorithm of the GBM and WKBJ/Maslov methods: kinematic and dynamic.The kinematic ray tracing determines the paths of rays, and the dynamiC ray tracing, which Cari be derived from Eikonal equation (Cerveny and Hom, 1980; Cerveny, 198 5 ) or.parabolic wave equation (Cerveny et aL, 1982; Cerveny and Psencik, 1983; Cerveny, 1987), calculates the dynamic properties of the wavefront in ray center coordinate system (Figure 1 ).The amplitude of the seismograms at a station is obtained from the summation of the contributions from all rays around it Fig. 1.Ray center coordinate system.t , e1 and e2 are unit vectors of this coordinate.ti s tangent to the ray and ei and e2 are normal and binormal to t. (s,qi,q2) describes a point at the vicinity of ray.
Recently Kim and Cormier (1990) developed another high-frequency asymptotic method, namely the vicinity ray tracing (VRT).A vicinity ray is defined as a ray in the neighbourhood of a central ray and each vicinity ray has a slightly different initial take-off angle with respect to the central ray (Figure 2).The VRT is based on the GBM, but is more accurate and simple than the GBM since it does not involve paraxial approximation and there is no speculation of how wide the spatial bound of each central ray should be.The locus of vicinity ray is determined much more accurately by integrating a new set of differential equations, which is referred to as the vicinity ray tracing system.The amplitudes of the rays in a bundle follow a Gaussian distribution centered at the central ray with the beamwidth defined by beam Fresnel volume.The beam Fresnel volume is in tum defined by a bundle in which the rays at the center and the boundary differ in travel times by half the characteristic period of this bundle (Kim and Cormier, 1990).
In spite of the elegance of the VRT theory, there is, however, lack of thorough consid erations in their original fomlUlation for synthetics.Since the summation of all central rays is according to the estimated travel time of each ray, proper travel time calculation should

485
persist for the whole ray set.Unfortunately, the travel time calculation in the VRT method fails when the normal intersection from station to a central ray cannot be found (details in section 2).This occurs for those rays between the station and the source, or the near side rays.In order to furnish a complete theoretical base for the VRT approach, a modified travel time calculation technique is introduced for the near side rays using the concept of paraxial ray.In addition, since the normal to the ray path is always confined in the vertical plane in a 2-D velocity structure, the basic wavefront parameters can be obtained directly from kinematic ray tracing system, which do not need complicated numerical manipulations as done by Kim and Cormier (1990) in their original algorithm.In the present paper, we first review the VRT method and point out the neglect of near side rays in Kim and Cormier's version (1990).We then introduce our new approach.This is then followed by synthetic wavefield computations that demonstrate the validity of this new approach.

VICINITY RAY TRACING
As proposed by Kim and Cormier (1990), in the VRT system wave equation in frequency domain can be written as : ( 1) where i=l ,2,----,n, represents an n-dimensional space.The high frequency zero-order asymptotic solution of equation (1) can be expressed in the form ( 2 ) where A and Tare the amplitude and phase function of Qi• By Fermat principle and variational calculus, we obtain Hamiltonian-Jacoby partial differential equation from which the eikonal equation in ray centered coordinate can be derived : ( 3 ) where is the scale factor (Aki and Richard, 1980).
Define a new variable hi as the angular difference between the tangential vector of a central ray and a vicinity ray in the t-ei plane in ray centered coordinate, and qi as the noIJilal distance from a central ray to a vicinity ray along the ei plane (Figure 2) (Kim and Cormier, 1990).Then Pi can be written as :

Vi
where Vl = v{s,ql,0) and V2 = v(s,0,q2) Some manipulation of (3) and (4) leads to where Equation ( 5) is called the VRT system. (4) (5) (6) In the VRT system, the curvature of the wavefront is defined in terms of qi and W.
From Figure 3 the curvature of the wavefront at point s(s,qi) is written as : where Ri is the redius of curvature of wavefront. (7) Alternatively, the curvature of the wavefront can be written in the wavefront coordinate system in which qi represents the distance from the. central ray to vicinity ray along the wavefront.The relation between q� and qi can be written as: and The radius of the curvature of wavefront (Ri) is used to calculate both the kinematic and dynamic properties such as travel times of stations and beam Fresnel volume in the VRT system.
The amplitude contribution of a ray to a station depends on the beamwidth of each ray.
In the VRT system, the principle amplitude and phase are contained in the ray within the beam Fresnel volume (Kim and Cormier, 1990).The Fresnel volume is first defined in optics (Jenkins and White, 1937;Stone, 1963) to explain diffraction.The definition of the beam Fresnel volume is that, for a given frequency, the cross-section of the beam Fresnel volume perpendicular to the central ray at point S is to contain all the vicinity rays that are within a half period difference (ahead or behind) in travel time relative to the central ray (Kim and Cormier, 1990).Mathematically, the beam Fresnel volume is given as : ( 11 where I is the half-period. The amplitude described by equation (2) then can be written as : where .D is the nonualization factor for radiation pattern of source, (]} is the product of reflection transmission coefficients (Aki and Richard, 1980), and y(t) is .the source time function .. The parameter K denotes the KMAH index (Chapman and Drummond, 1982) that represents the 90° phase shift when the ray hits an x-caustic point.It signals changes of the sign of the radius of the wavefront curvature.
The whole wavefield now can be obtained by the superposition of all displacements and be given as : where 6 is the vertical take-off angle and </> the azimuth of a single ray.
The travel time at a station near a central ray in the VRT system is calculated using an approximation method from the arrival time of each central ray (Kim and Cormier, 1990).
From Figure 3, the travel time of the station at N(s,ni) could be specified by : where .6.n, the normal to the wavefront, is the distance measured from the wavefront to the station.e denotes a plus/minus sign representing convex/concave or planar wavefront to the station.V is the vclocily of vicinity ray at S(s,qi).Let ni be equal to the normal distance between S(s,0) and station, .6.r then can be written as : (16) Equation ( 16) is only valid for the velocity structure for which the locus of S(s,0) is regular enough to avoid the rapid change of wavefront.
Although this approximation method is accurate compared with the DRT system, it works only in two situations.The first situation is that the central ray is normal incident to the ground and the vicinity ray is not (Figure 3).In this case, the distance perpendicular to the central ray from station N(s,ni) to the central ray at S(s,0) can be found and 'rf i at S(s,0) is not equal to zero to ensure that the radius of the the wavefront Ri is finite.In the second situation when the central ray is not nonnal incident, this approximatf o n can only be applied for the station on the near side of the central ray (the side close to the source) and for the vicinity ray on the far side of the central ray (Figure 4).In the above two situations , it is ensured that the S(s,q;) on the vicinity ray and S(s,0) on the central ray can be found.For the station that is on the far side of the central ray, this method fails since the point S(s,0) corinecting to the station can never be found (Figure 5).This drawback will be eliminated in the present work.

MODIFIED TRAVEL TIME CALCULATION
Owing to the fact that the travel time at a station estimated from each central ray is used for sununing the contributions of all rays in order to calculate synthetic seismogram, it is necessary to have a general method of travel time calculation applied prolJerly to all central rays distributed in a broad region encompassing the station.We propose a modified method to calculate the travel time.We assume that ray.paths near the ground surfaee are nearly parallel to each other.Taking a line (NA) through station (N) which parallels the local tangent of the central ray at surface and normal to ray coodinate axis (ei) of central ray at point A (Figure 6), the travel time t(s, ni) can be written as : (17) where ()is the angle between the incident ray and the Z-axis, s1 is the distance between S(s,qi) and S(s,ni), and r(s,0) is the travel time of central ray to the ground surface.
In this modified method, the main source of error is due to the determination of point A.
Point A defined as mentioned before is acceptable based on the following two reasons.One is that the seismic velocity near the surface is always slower than those of deeper structure, and consequently the incidence angle of ray path is nearly vertical and ds1 is small.A short distance from N to A will reduce the approximation error.The other reason is that the radius of wave front is very large compared to the wavelength when the wavefield propagates a long distance from source.The large radius of wave front also reduces .theerror.Therefore, the error can still be tolerated when the incidence angle is not nearly horizontal or the radius of the wavefront Ri is substantially larger than Oi.A numerical test of this modified method will be shown in the next section.These equations can be solved by standard numerical tec hni que such as Runge-Kutta method.
For determining the two q uanti ties qi and T/i of each ray, Runge-kutta method is also used to solve equation (5).In this way, the q uantities from equation ( 5) have to be stored at every step of ray path in order to calcu la t e the wavefront curvature Ri of a specified point S on the central ray.This is a time-consuming p rocedure.Instead, we take advantage of the 2-dimensionality of the velocity structure, and estimate the vicinity ray from equation (20) with the initial angle d iffere nce T/o with respect to the correspondi ng central ray.The quantity "qi" is es tima ted from the distance between central ray and vicinity ray and "T/" is obtained by estimating the angle difference between the local tangent of the central ray and • vicinity ray.This is much more efficient than solving eq uatio n (5) directly.It is noticed that the vi cini ty ray should be on the �ar side of the central ray in order to guarantee that the qi and T/i could be found along the whole ray path.
Following the ab ove improvemen t, the travel time at a station and the beam Fresnel volume of rays can now be c alcula ted using equation ( 17).Then the amplitude of a ray can be ob tained from equation ( 12) and so is the wavefield of a ray to a station.The KMAH index K in equation (12) records how many times the ray touches the caustics along the whole ray which is manifested by a phase shift of 90°.
The synthetic seismograms then can be obtained by sum min g up the contributions of rays to a station (equation (13)).

TESTS
Several numerical tests have been carried out to assess the accuracy of modified travel time calc ulatio n and alternative procedure.One test is for error estimation of travel time.The others are for s tability of modified VRT system by calculating synthetic seismograms in various velocity structures.The prev ailin g freq ue ncy of the Gaussian source wavelet is assumed to be 5 Hz and the width of each wavelet is ch ara ct erized by the parameter of /=4 as defined in (Cerveny et al., 1977) for all tests.

Travel Time
In this test (A), the structure is a homogeneous half space with ve loci ty of 4 km/sec.An explosive line source is l ocated at a depth of 50 km and a receiver is at the horizontal distance of 50 km.Figure 7 is the ray diagram by two-point ray tracing method in which 9 rays are emanated from the source and arrive at the surface every 5 km from 25 km to 75 km.The initial values of q0 and 770 are 0 and 2°0 respectively.Equations ( 14) and ( 17) are used for calculating the travel time of the rays on the far side and the near side of th�_ receiver, respectively.The result of travel time calculated from central rays is shown in Figure 8.The travel time of ray 5 is taken as the standard value, since this ray arrives right at the receiver.
From Figure 8, we can see clearly that the travel times calculated from rays on the fa r side (rays 6 to 9) are more accurate than those on the near side (rays 1 to 4).Besides, due to the fact that the radius of the wavefront of the ray on the near sic!e is smaller than that on the far side, the amplitude contribution of each ray to receiver decreases more rapidly on the far side.The synthetic seismogram at 50 km with summations of all rays (1 to 9) is also shown in Figure 8.
It should be noted that even more approximations are entertained for rays 1 to 4, the fa rthest ray away from the receiver only suffers from an error of 0.5 percent.In addition, correct radius of wavefront for amplitude computation can be obtained by this modified method.From this test, the error of the travel time calculation is still satisfactory since the waveform remains clean and is free of high frequency wiggles or noise in comparison with standard Gaussian waveform proposed by Cerveny et al. (1977).Figure 12 is the corresponding synthetic seismograms obtained by Kim and Cormier (1990) for comparison.These two synthetic seismograms are quite similar except that .thedecrease of amplitude near the caustic is more clear in Figure 11 than in Figure 12.This is because all contributions of rays have been summed up using Equation (17). Figure 13 shows the synthetic seismograms with different 170 ( 0.5°).and they are almost the same as those in Figure 11.It indicates that the initial angle difference between central ray and vicinity ray (170 ) becomes an insensitive parameter in computing synthetic seismogram when it reaches a certain small value.
In Test C, the velocity structure model is designed to simulate the common geological situation in which a low velocity zone is present (Figure 14).The low velocity zone is set at a depth between 10 km and 15 km .The ray diagram is shown in Figure 15 in which 30 rays emanate from source and the angle step of central rays is 2° .The initial condition is 0 for fJo and 1° for 170• Owing to the positive and negative velocity gradient, a shadow zone and a triplication zone are located in the regions fr om 48.3 km to 67.4 km and from 67.4 km to 79.4 km respectively (Figure 15).The synthetic seismograms are shown in Figure 16.The amplitudes fr om 50 km to 65 km, do not vanish abruptly but decrease as the source receiver distances increase.At 70 km, there exists a great amplitude signal, resulting from the triplication of rays due to a caustics.This example demonstrates that the VRT system is able to solve the problem with both shadow zone and triplication zone present simultaneously.
In Test D, the model contains both vertical and horizontal gradient discontinuities to simulate the presence of faults.The velocity gradient is shown in Figure 17.Forty rays are emanated fr om the source with angle step of 1.5° between each central ray.The initial condition is Tf o of 0 and T/o of 1.5°.From the ray diagram (Figure 18), strong refraction is dottedlin e:vic.ray X distan ce [km] solid line : cen.ray 0.
5. l O. 15. 20. 2 . 5. 30,35. 40. 45. 50. 55. BO. 85. 70. 75.BO.TRAVEL TIME(SEC) 8.25 9.62 . 11.00 12.37 13.74 15.12 16.49 17.86 19.23       evident when rays propagate through the vertical velocity discontinuity at 20 km.Besides, the ray-focusing occurs at about 20 km.These two phenomena are due to the strong lateral heterogeneity as posed by the large gradient changes.Figure 19 shows the VRT synthetic seismograms.Multiple caustics in the ray diagram and the three main branches on the seismogram are noticeable due to the laterally inhomogeneous structure.The first branch is contributed from the rays that do not penetrate through the high velocity region.The second is contributed from the rays with KMAH of 0 but which do penetrate through the high velocity region.The amplitude of the first branch is greater than that of the second.This is because that the amplitude contribution of each ray to the first branch is larger with shorter propagating distance.The third one is contributed from the rays with KMAH of 1, of which the 90° phase shift of waveform is clear.

5.
The results of all the tests indicate that the modified VRT method in this study is capable of generating reliable synthetic seismograms in a complicated 2-D velocity structure.It should be emphasized that the above calculations employ the new procedure which utilizes kinemetic ray tracing system instead ofVRT system to determine the basic wavefront parameters, which saves large amounts of computing times without losing accuracy.

TRAVEL TIME (SEC)
� t----.1 " � 1----------------¥-.J 3l: Laterally heterogeneous structures are of great interest in geophysics; the earth crust, decending lithospheric slabes in the upper mantle, and the region of core-mantle boundary are of this catagory.The VRT system provides a new high-frequency asymptotic procedure in calculating the wavefield in this kind of medium.Calculating synthetic seismogram in two or three dimensional velocity structure by applying VRT system has many advantages such as it is more accurate for no paraxial approximation, faster calculation and giving clear definition of beam width.

•
The approximation method of travel time in VRT system developed by Kim and Cormier (1990), however, can not be used in certain cases.The modifi ed method of travel time calculation proposed in this study is much more useful and can be applied to the whole ray set for obtaining better result.In spite of the possible extra errors accompanying with the additional approximations made in the modified approach, we do not see any degradation in the quality of the computed synthetic waveforms.In this study, we advocate a direct calculation procedure operated in kinemetic ray tracing system instead of solving dynamic VRT system.This improved procedure has proven to be efficient in saving computation times and maintaining accuracy simultaneously.

Fig. 3 .
Fig. 3. Travel time of station can be estimated by central ray.(central ray is normal incidence)

Fig. 4 .Fig. 5 .
Fig. 4. Travel time of station can be estimated by central ray .(central ray is not normal incidence, see text)

5. 2
Synthetic Seismogram Three tests of different models B, C, and D are carried out for synthetic seismogram.

Figure 9
Figure 9 is the velocity model for Te st B. The P-wave velocity increases monotonically with two different gradients .A discontinuity is set at the depth of 10 km.64 rays are emanated from source with 1° angle step.The ray diagram with 170 of 1° and q0 of 0.0 is shown in Figure10.From this diagram, two caustics exist at 32.4 km and 46.3 km, between which • the triplication will take place.The synthetic seismograms of vertical component calculated in this study is shown in Figure11.Figure12is the corresponding synthetic seismograms obtained byKim and Cormier (1990) for comparison.These two synthetic seismograms are quite similar except that .thedecrease of amplitude near the caustic is more clear in Figure11than in Figure12.This is because all contributions of rays have been summed up using Equation (17).Figure13shows the synthetic seismograms with different 170 ( 0.5°).and they are almost the same as those in Figure11.It indicates that the initial angle difference between central ray and vicinity ray (170 ) becomes an insensitive parameter in computing synthetic seismogram when it reaches a certain small value.In Test C, the velocity structure model is designed to simulate the common geological situation in which a low velocity zone is present (Figure14).The low velocity zone is set at a depth between 10 km and 15 km .The ray diagram is shown in Figure15in which 30 rays emanate from source and the angle step of central rays is 2° .The initial condition is 0 for fJo and 1° for 170• Owing to the positive and negative velocity gradient, a shadow zone and a triplication zone are located in the regions fr om 48.3 km to 67.4 km and from 67.4 km to 79.4 km respectively (Figure15).The synthetic seismograms are shown in Figure16.The amplitudes fr om 50 km to 65 km, do not vanish abruptly but decrease as the source receiver distances increase.At 70 km, there exists a great amplitude signal, resulting from the triplication of rays due to a caustics.This example demonstrates that the VRT system is able to solve the problem with both shadow zone and triplication zone present simultaneously.In Test D, the model contains both vertical and horizontal gradient discontinuities to simulate the presence of faults.The velocity gradient is shown in Figure17.Forty rays are emanated fr om the source with angle step of 1.5° between each central ray.The initial condition is Tf o of 0 and T/o of 1.5°.From the ray diagram (Figure18), strong refraction is

Fig. 7 .Fig. 8 .
Fig. 7. Ray diagram of Test A square is position of source, solid lines are central rays and dotted lines are vicinity rays.

Fig. 9 .
Fig. 9. Velocity model of Test B. A velocity gradient discontinuity occurs at depth 10 km.
Ray digram of Test B. A triplication zone is located in the range between 32.4km and 46.3km.