First arrival traveltime methods have proven to fail to image deeper structure, when most of the reflected energy is transported by later arrivals. The family of ray tracing methods is one of the few available alternatives to generate Green's functions with multiple arrivals. Dynamic ray methods also compute amplitudes and phase-shifts. Even in rough structures, ray methods can be used to compute useful Green's functions, if the model is properly smoothed. We discuss ways to obtain a good image with Green's functions computed with the paraxial ray tracing method in the 2-D Marmousi model. Our paraxial ray tracing implementation uses the paraxial and regularity conditions to limit the region where the paraxial solution for a specific ray is useful. Smoothing the model, does not only make the model more ``ray valid'', but also makes the high frequency wavefield more like the band-limited solution. We apply a merging of phases technique to reduce two nearby arrivals into one. This technique reduces the high frequency noise in the image significantly at a similar cost as using just the most energetic arrival. Smoothing over a length proportional to the local wavelength of the center-frequency, results in a better image than just smoothing over a constant length. Using the two most energetic arrivals, instead of just the most energetic, improves the image of the ``target'' zone.
INTRODUCTION During the past few years, there has been a significant increase in applying 3-D data recordings and 3-D imaging. 3-D imaging is computationally intensive, and much effort has been expended to find cheap and efficient ways to do the full 3-D prestack Kirchhoff depth migration. First arrival traveltime methods are not always sufficient, and the family of asymptotic ray methods is, especially in 3-D, one of the few alternatives. Dynamic ray tracing methods are not as fast as many first arrival traveltime methods. However, recent developments demonstrate that there are ways to make the ray computations more efficiently, e.g. Vinje et al . A more serious problem is that for rough models, the ray validity conditions are no longer fulfilled, e.g. Ben-Menahem and Beydoun . The ray methods are high frequency approximations, which are only valid if the model varies smoothly over a wavelength. To make rough models more ``ray-valid'', they have to be smoothed. Versteeg experimented with smoothing in the Marmousi model, and found that the velocity could be smoothed with a Gaussian bell shaped operator with radius up to 200 m, while still obtaining a reasonable image with shot profile finite-difference modeling. Another way to compute Green's function in a rough model, is to use band-limited methods. A recent method () computes band-limited Green's functions by solving the Helmholtz equation in 2-D polar coordinates for a few frequencies.
In this paper, we review the ray validity conditions, which are important to understand the advantages and shortcomings of the ray tracing methods. We have computed Green's functions with the paraxial ray tracing method in the smoothed Marmousi model. The ray validity conditions were not fulfilled, even in the smoothed model. However, the paraxial ray tracing method was still able to compute reasonable and accurate traveltimes for multiple arrivals. We have experimented with both a constant length smoothing operator, and an operator with length equal to the wavelength of a center frequency of 30 Hz. Smoothing is preferred, not only to make the model more ray valid, but to make the ray tracing solutions more like band-limited solutions, e.g, Biondi . We demonstrate, with finite difference modeling, what effect smoothing of the model has on a high frequency wavefield. We use the images, obtained from feeding the Green's functions into a Kirchhoff prestack depth migration algorithm, to evaluate the performance of the paraxial ray tracing solutions. Because high frequency finite difference modeling demonstrates that different arrivals often interfere, we use a simple technique to merge interfering ``ray-arrivals'' into one. This technique reduces some of the high frequency noise in the images.
RAY METHODS AND THEIR VALIDITY CONDITIONS
Asymptotic ray methods are based upon the high frequency assumption, i.e. both the structure and the amplitude of the wavefront are assumed to be slowly varying over a wavelength. In seismic applications, this condition is often violated. Ben-Menahem and Beydoun and Beydoun and Ben-Menahem have given a rigorous study of the ray methods validity conditions in an elastic isotropic medium, and they defined several conditions that should be fulfilled: 1) The mode decoupling condition expresses the condition for decoupled P- and S-waves. 2) The high frequency condition expresses when the traveltime can be computed from the eikonal equation. 3) The cumulative error condition expresses how far in traveltime a ray is accurate, which is also called the Fresnel condition. For paraxial rays, Ben-Menahem and Beydoun also introduced the paraxial condition and the regularity condition. In the paraxial ray tracing the solution is computed for points on, and also in the vicinity of, a central ray. We use the two latter conditions to decide how far from a central ray the paraxial solution is valid. Here we give a brief review of the ray validity conditions and discuss their importance.
In an isotropic medium, the elastodynamic equation can be written as:
In ray tracing the solution is approximated with the zeroth order term in the ray series Rayseri. Inserting this solution into equation Helmhol gives:
Solving equations Eikon and Transport gives the following solution, e.g. Ben-Menahem and Beydoun
To compute the traveltime at points in the vicinity of the central ray, we must compute the wavefront curvature. For this purpose it is convenient to introduce two new coordinate systems: the ray-centered coordinate system, which is a curvilinear coordinate system that follows the rays, and the ray coordinate system where and can be the rays take-off angles at the source. A system of first order linear differential equation for the solution of the eikonal equation in the ray-centered coordinate system, can be obtained by applying a paraxial approximation to this equation ():
The ray methods are far-field approximations, and both the amplitude and the unit slowness vector must be slowly varying. This condition is expressed in the far-field condition:
Sub-index denotes value taken at the source. C is a factor dependent upon the source, and is constant along the ray.
IMPLEMENTATION OF THE PARAXIAL RAY TRACING METHOD
In 2-D, the 14 kinematic and dynamic ray equations defined in equations Kineq1-Dynameq, reduce to 6. We chose a Runge Kutta method to solve the 2-D kinematic and dynamic equations. A bicubic spline representation is used for the P-velocity model. A more efficient choice would have been to represent the velocity model with constant gradient slowness- or velocity-field within triangles. (). Then the ray equations can be solved analytically. Also other velocity representations that give analytical solution, exist.
The paraxial and the regularity condition were used to find the maximum distance from the central ray where the solution was valid. An upper limit of this distance was set to 250 m. This was done to avoid problems in situations where the ray is straight and the wavefront curvature, computed from the second derivatives of traveltime in equation Hessian, is zero. Then the regularity and paraxial validity conditions are not sufficient to limit the region of validity. Figure one-ray shows an example of where the solution is valid according to conditions Paraxial and Regularity: the left-hand sides of the these inequalities were set to be less than 10 percent. The figure shows the amplitudes variations, computed along the central ray, within this region. The solution cannot be valid in the outer parts of this region, since a neighboring central ray will give another result here. The paraxial and regularity conditions are also approximations that are only valid in the vicinity of the central ray.
Since there are shadow zones, our implementation of the paraxial ray tracing always computes one arrival-time, even if that solution don't fulfill the paraxial or the regularity condition. If the distance between two neighboring rays exceeds a limit (75m in our experiments), a new ray was computed from the source with take-off angles in between. One or two of the most energetic arrivals, computed from rays closer than a distance of 1.5 times the distance from the actual point to the nearest ray, were picked. Traveltime, amplitude, and phase, were stored for each pick. Figure Green-tab shows the Green's functions computed with the paraxial ray tracing in the Marmousi model. The most energetic arrival has been picked. The source was located at the surface, with a distance of 5200 m from the left edge of the model. We clearly see that the traveltime functions are discontinuous. And for rays that pass through caustic regions, the amplitudes become far too high. In one implementation of the program, in which two arrivals near in time were merged together, we assumed a center frequency of 30 Hz to determine if two arrivals with different phases interfere. We assumed that the pulses had the simple form of a cosine with a sine as its Hilbert transform. If the pulses interfere, a new amplitude, phase, and traveltime were computed from the merged arrival. This allows us to operate with phase shifts that are not only in integer multiples of . This is demonstrated in Figure Green-tabc where we see more variations of the phase.
GREEN'S FUNCTION COMPUTATIONS
In this section, we have compared results of modeling with a 2-way finite difference method and the paraxial ray tracing method. Figure fd-seismo shows synthetic seismograms recorded at the surface for a shot location at depth 2500 m and 6000 m from the left edge of the model. The seismograms were computed with a two-way acoustic finite difference program. In Figure fd-seismoa, the seismic bandwidth was used. In Figure fd-seismob a high frequency bandwidth was chosen to simulate a high frequency method behavior in the model. We see some shadow zones in the wavefield, and the energy is distributed differently. We see more arrivals, and some of the high-energy late arrivals are not as clear as in the low-frequency solution. Figure fd-seismoc was obtained in a smoothed version of the Marmousi model. The model was smoothed using a triangle smoothing operator with a fixed length of 200 m. The smoothing reduced the number of arrivals, and made the dominating arrivals more clear. The solution became more like the seismic frequency band solution. Based upon a heuristic idea of smoothing over the local Fresnel-zone, we also applied a smoothing operator with a variable length. In Figure fd-seismod, a Gaussian operator with a half-width equal to the wavelength for a center frequency of 30 Hz was used to smooth the slowness field. This resulted into a smoother wavefield than when we smoothed over a constant length operator: the energy distribution is different. Note however that the traveltimes are more or less the same.
In Figure sm-comp, the traveltimes, amplitudes and phases are plotted for the two most energetic arrivals. The arrivals were obtained by picking the traveltime values obtained by shooting from different locations at the surface. They are computed from the paraxial ray tracing method in the smoothed versions of the Marmousi model. Note, that only in some areas has the program detected more than one arrival. We see that the paraxial approximation is quite successful in computing the most energetic arrivals, even in a rough model such as the Marmousi. The most significant difference between these two kinds of smoothing is the amplitude changes. And the amplitudes determine which traveltimes are the most energetic. This result implies that a better way to smooth the model, might improve the amplitude calculation.
IMAGING OF THE MARMOUSI MODEL The Marmousi dataset was generated at the Institut Français du Pétrole, and used for the workshop on practical aspects of seismic data inversion at the 1990 EAEG meeting in Copenhagen, where contractors, universities, and oil companies tried to determine a correct image from this data set. The geometry of the Marmousi is based on a profile through the North Quenguela in the Cuanza basin. The target zone is a reservoir located at a depth of about 2500 m (). The data set consists of 240 shots with 96 traces per shot. The shot and receiver spacing are 25m, and the near offset is 200 m. The first shot position was located at 3000 m. Recorded time was 2900 ms, with a sampling interval of 4 ms. The Green's functions computed from each shot and/or receiver location were pre-calculated with the paraxial ray tracing program described above. The traveltime, amplitude, and phase values were stored at grid-points, with 12.5m spacing in depth and 25m spacing in the horizontal direction. Values for 2000 m in each horizontal direction from each surface location and 3000 m in depth were stored. The imaging was done with a Kirchhoff prestack depth migration algorithm. Figure close-up shows close up images of the target obtained from different implementations of the paraxial ray tracing code, and different smoothing of the model. First, the slowness field was smoothed over a constant length triangle with a baseline of 200 m. Amplitude, phase, and traveltime information are used. Figure close-upa is computed using the most energetic arrival. We can clearly see the reservoir. However, the image is a bit blurred, and the target zone is not flat. In Figure close-upb, the two most energetic arrivals were used. We see that the target zone becomes more flat and visible.
In Figure close-upc, the merging of phases technique is used when we compute the most energetic arrival. This technique remove much of the high frequency noise in the image. The variable length smoothing operator was used in Figure close-upd for the most energetic arrival, and this also resulted in a more flat target-zone. Also here was the merging of phases technique applied.
In Figure images-1, we show the difference in the images obtained just by using the most energetic arrival and the image obtained using the merging of phases technique, variable length smoothing and the two instead of the one most energetic arrivals. The merging of phases technique removes a lot of the high-frequency noise in the figure. Using two instead of one arrivals flattens the image of the reservoir, and the base of the reservoir becomes more visible. The variable-length smoothing operator also improves the image.
In Figure images-2, we demonstrate that the amplitude and phase information computed from paraxial ray tracing is important. The amplitudes can not only be used to determine which arrival is the most energetic, but using Green's functions that have amplitude and phase, in addition to traveltime information, improve the image significantly. The merging of phases technique was applied to compute the most energetic arrivals with amplitude, phase, and traveltime values. The images demonstrate the improvements made in also exploiting phase and amplitude information. Also using both phase and traveltime improved the quality of the image compared to using traveltime only.
The differences are more visible, when we see what will be summed in the migration algorithm. In Figure crp-amp, the common reflector points image for the surface position 6450 m is shown for the two images in Figure images-1. In the best result obtained, the reflected energy at the target zone is much more coherent over offset than in the case where we just used the most energetic arrival in a model with a constant smoothing operator. This is most obvious for the near offsets. The values are obtained by stacking over the common offset diffraction hyperbolas, and dividing by the sum of the multiplied source and receiver Green's function amplitude over the same hyperbola. In our amplitude calculations, we did very little to correct for far too-high amplitudes due to caustics. Hence, a very high amplitude with a wrong phase could distort the image.
Figure crp-tpa clearly demonstrate the improvements possible by exploiting amplitude and phase information in the Green's function. In the case of traveltimes only, we can hardly see the target. When we exploit traveltime and phase information, the coherency at the target zone improves significantly. Even if the target zone becomes much clearer when we exploit traveltime, phase and amplitude, we loose some of the coherence at the target zone ( Figure crp-tpaa). The only difference between Figures crp-tpaa and crp-tpab is the amplitude factor in the summation. This suggests that a better control of the amplitude near caustics is important.
CONCLUSIONS In this paper, we have studied the potential of using the paraxial ray tracing method to compute Green's functions in rough models, where the ray validity conditions are not fulfilled. Finite difference modeling in a smoothed Marmousi velocity-model shows that the high frequency wavefield (80-200 Hz) is closer to the band-limited wavefield around 20-60 Hz in the un-smoothed model. Hence, proper smoothing can lead to a ray solution closer to a band-limited solution. Experiments with the paraxial ray tracing in smoothed Marmousi models demonstrate that this method is able to compute the traveltimes fairly accurate, and that the amplitude is accurate enough in most cases to successfully pick the most energetic arrival. We found that the amplitudes were more sensitive to the smoothing than the traveltimes. The images demonstrate that exploiting amplitude and phase information together with the traveltimes, improves the quality of the images. Comparisons of common reflector point gathers obtained from exploiting amplitude, traveltime, and phase and exploiting traveltime and phase, suggest there is a potential for an improved result by better control of the amplitudes near caustics. We also introduced a technique to merge two arrivals near in time into one, if they interfered. Results demonstrate that this merging of phases technique reduces much of the high frequency noise in the images. Using the two most energetic arrivals instead of just the most energetic arrivals improved the image of the target zone. Smoothing of the slowness field over one wavelength of the center frequency, instead of applying a constant-length smoothing operator, also improved the image. The most significant improvement was from the merging of phases technique combined with a variable length smoothing operator. The experiments with the simple smoothing operators used here, demonstrate the need for further research on that topic.
We wish to thank Dave Nichols for many helpful discussions, and to provide us with his migration code.