Abstract
An initial study is performed in which dynamically focused Gaussian beams are investigated for seismic imaging. Focused Gaussian beams away from the source and receiver plane allow the narrowest and planar portions of the beams to occur at the depth of a specific target structure. To match the seismic data, quadratic phase corrections are required for the local slant stacks of the surface data. To provide additional control of the imaging process, dynamic focusing is investigated where all subsurface points are specified to have the same planar beam fronts. This gives the effect of using nondiffracting beams, but actually results from the use of multiple focusing depths for each Gaussian beam. However, now different local slant stacks must be performed depending on the position of the subsurface scattering point. To speed up the process, slant stacking of the local data windows is varied to match the focusing depths along individual beams when tracked back into the medium. The approach is tested with a simple model of 5-point scatterers which are then imaged with the data, and then to the imaging of a single dynamically focused beam for one shot gather computed from the Sigsbee2A model.
1. Introduction
Dynamically focused Gaussian beams are investigated for the seismic imaging of common-shot reflection data. This extends the work of Nowack [1] in which focusing of Gaussian beams propagated from the source, and receiver surface is used to allow the narrowest portions of the beams to occur at the depth of a specific target structure. The beam fronts at the beam waists are also planar leading to more stable beam summations for imaging. To match with the surface data, quadratic phase corrections are required for the local slant-stacks of the surface data. However, in the earlier approach only a single focusing depth can be specified. To provide additional control on the imaging process using focused Gaussian beams, a dynamic process is investigated where all subsurface points are specified to have planar beams of the same width. This gives the effect of using nondiffracting beams at the scatterers, but actually results from the use of multiple imaging depths for each Gaussian beam. As a result, there is a tradeoff in speed where now different slant stacks of the data are required for each subsurface point. In order to speed up the imaging process, the slant stacking of the data window is varied to match the focusing depths along individual beams when tracked back into the medium. These focused slant stacks are then matched to Gaussian beams with beam waists located at different subsurface locations. The approach is first tested using 5-point scatterers which are then imaged with the data, and then using a single dynamically focused beam for one shot gather computed from the Sigsbee2A model.
Summations of Gaussian beams have been applied for the computation of high-frequency seismic wavefields in smoothly varying inhomogeneous media (see e.g., Popov [2]; Cerveny et al. [3]; Nowack and Aki [4]). Reviews of Gaussian beam summation have been given by Cerveny [5, 6], Babich and Popov [7], and more recently by Popov [8], Nowack [9], Cerveny et al. [10], and Bleistein [11]. An advantage of summations using Gaussian beams to construct more general wavefields is that the individual Gaussian beams have no singularities along their paths, no two-point ray tracing is required, and triplicated arrivals are naturally incorporated into either forward or inverse modeling. More recently overcomplete frame-based Gaussian beam summations have been developed based on window and wavelet transforms to address some of the issues related to completeness of beam summations [12]. In an overcomplete frame-based approach, the wavefield is decomposed into beam fields that are localized both in position and direction. Although an orthonormal basis cannot be formed using a Gabor frame, an overcomplete frame expansion can be constructed which has the added benefit of providing redundancy in the expansion [13–17]. Here curved initial beams are used to decompose the data and these are then propagated into the subsurface using dynamically focused Gaussian beams. An alternative is to shoot beams from the scattering locations upwards to the source and receiver plane as was done by Popov et al. [18], but this will be computationally slower than shooting beams from the surface.
In order to test the dynamically focused Gaussian beam approach, a model is constructed with 5 small scatterers with depth in a vertically varying medium resulting in five diffracted arrivals each with different move-outs with distance. The partial image is then given for a single beam with different focusing depths such that the planar beam waists at the subsurface points are the same. However, when summed over all angles and data windows, the scatterers are properly imaged. This is then tested for a single dynamically focused beam for one-shot gather computed from the Sigsbee2a model. The advantage of using dynamically focused beams is the consistency of having the same focused beams for all subsurface points. A disadvantage is that now multiple slant stacks are required for different points of the subsurface which increases the computational burden compared to a standard Gaussian beam algorithm or a beam algorithm with a single focus depth. This can in part be offset by precomputing the beam stacks at a set number of focusing depths and then interpolating the results required at different subsurface points.
2. Gaussian Beam Imaging with Dynamically Focused Beams
In common-shot migration, each shot gather is migrated separately and the results are summed to give the final image . Thus, where the adjoint image for each shot point can be written as where is the radial frequency, is the source time function, is the geophone location, is the source location, and is the scattered field recorded at the geophones. and are the Green’s functions, from the geophones and source to the subsurface scattering points . From reciprocity of the Green’s function and indicates the complex conjugate of .
In the 2D case, the Green’s function can be written in terms of a summation of Gaussian beams as where is the complex beam parameter [3] and
The coordinates correspond to the subsurface position in ray-centered coordinates, is the travel time along the central ray, is the velocity along the central ray, and the horizontal component of the ray parameter vector at the initial point is , where is the take-off angle from the vertical.
The complex second derivative of the travel time field with respect to the transverse coordinate can be written as where is related to the wavefront curvature of the beam by . To form a bounded Gaussian beam, then . The variables and are solutions to the dynamic ray equations and for a beam solution are also complex [19]. The dynamic ray equations in 2D have two real fundamental solutions which can be written as where . , are solutions for an initial plane wave and , and are for an initial line source. The inverse of is
The two real solutions of the dynamic ray equations must then be combined to form a Gaussian beam. There are a number of ways to combine the solutions, but one way is [3] where is the beam parameter. The variable is related to the complex geometrical spreading along the beam. Since for the fundamental solution matrix for all points along the ray, the complex geometric spreading can never be zero at any point along the beam if it is nonzero at any one point. Since the Gaussian beam amplitude is related to the inverse square root of the geometric spreading, the beam amplitudes are always finite, even at caustics for the ray solution. This is one of the useful features of Gaussian beam solutions, in contrast to ray solutions.
The second derivative of the time field with respect to can be written as
Since at the source point and , then where is the magnitude squared of . Alternatively, can be written in terms of as where is the magnitude squared of . Thus, the complex beam parameter can be specified either directly in terms of or in terms of .
The complex beam parameter can also be written where is an initial velocity. In a homogeneous medium, is the distance of the beam waist from the initial point of beam.
An alternate, and in some sense simpler way, to blend the solutions to the dynamic ray equations is where
In either case, we can obtain at the scattering point, from the initial value of along the source and receiver array by solving the dynamic ray equations for and , and then forming . Alternatively, we can specify at the scattering point and then use the solution matrix for the dynamic ray equation from to to obtain and at the initial point of the ray. From these, one can obtain the corresponding value of beam curvature at the beginning point of the ray. At a point along the beam, can be written as where is the wavefront curvature in ray coordinates, is the velocity along the ray, and is related to the transverse beam width, and these can be specified at the scattering point along the ray and used to determine the value at the beginning point of the ray.
The exponential term away from the central ray can be written as where is the beam half-width transverse to the ray. At the initial point of the beam, . For the case of the beam waist at the initial point of the beam , and is the beam half-width at the beam waist. In a homogeneous medium, this is the narrowest point along the beam and is also the only point where the beam front is planar. For the case when where , the beam front is planar at some point along the ray and generally curved at the initial position . Also, the beam waist is shifted along the beam away from the initial point of the beam.
Although the planar beam waist is often placed at the initial source point for forward modeling, it is also common to put the beam waist at the receiver location [5, 6]. This reduces the number of beams required for the summation at the receiver, and also planar beam fronts at the receiver provide more stable beam summations. Recent true amplitude migration formulations using Gaussian beams have used beams launched directly from the scattering points up to the surface with the beam waists specified at the scattering points [18, 20, 21]). However, it is more economical to launch beams from the source and receiver positions down into the subsurface since there are fewer source and receiver locations than subsurface scattering points, and this minimizes the amount of beam tracing required. In order to locate the beam waists in the subsurface when the beams are launched from the source and receiver aperture, then generally curved beam fronts are required along the source and receiver aperture.
For dynamic beam focusing in Gaussian beam imaging, the beams are launched from the initial source and receiver plane, but then are dynamically focused for each of the individual scattering points. Although the beams are generally nonplanar along the receiver plane, the initial beam parameters can be determined from the specified beam parameters at depth. However, now these curved beams on the initial surface must be matched to the local slant stacks of the data for beam parameters specified at each of the subsurface scattering points.
Assuming that the initial beam parameters have been determined from dynamic ray tracing for beam parameters specified at the scattering points, then the algorithm for imaging-focused Gaussian beams from Nowack [1] can be extended to dynamic focusing. For generally nonplanar initial beams at the source or receiver locations launched at some angle to the aperture plane, the quadratic part of the initial beam with respect to the horizontal coordinate can be written as where is the initial real-valued horizontal beam curvature and is the initial horizontal beam half-width at the reference frequency . To match this with the initial parameters of the beam propagated into the medium, then the transverse coordinate of the beam where is the angle of the beam with respect to the vertical. Given the initial values and along the source and receiver aperture, then the initial values for and for the beams are and then the initial beam parameter can be obtained and used to construct the beam solution propagated into the medium.
The 2D resolution of unity by Gaussian functions in the aperture plane is where
Assuming a regularly spaced set of beam centers along the receiver array, the initial locations of the Green’s functions at the receivers can be phase shifted to these beam center locations with a phase adjustment of . Then, where represents the terms in front of the last integral in (3) for the receiver Green’s functions and is the horizontal component of the ray parameter vector along the receiver aperture. The imaging formula can then be written as where
The source Green’s function now needs to be decomposed into Gaussian beams, but for simplicity here will be just referred to as , where the over-bar signifies the complex conjugate
The common-shot imaging formula for nonplanar, focused beams along the aperture plane can then be written as where
This is a local slant-stack of the data with a quadratic phase correction term to match the data with the beams launched into the medium. The standard Gaussian beam migration formulas with the planar beam waists along the aperture plane do not include this quadratic phase correction term [14–17]. The beam centers are spaced along the receiver aperture plane at and the initial horizontal beam widths are where is the horizontal half-width of an initial Gaussian function at the reference frequency . The spacing of the beams in horizontal position and launch angle can then be determined based either on physical reasoning [14–17] or by arguments based on frames [13]. However, for curved beams along the initial plane, the data windows need to be broader to account for the wider beams at the surface.
3. Applications of Dynamically Focused Gaussian Beam Migration
In order to test the dynamically focused beam migration formulation, two examples are given. The first application has 5 compact sources located at depths of 8,000, 12,000, 16,000, 20,000, and 24,0000 ft at a distance of 40,000 ft from the left side of the model. The background velocity model has two layers. The first layer has a thickness of 6000 ft with a constant velocity of 5000 ft/sec. The second layer goes from 6000 ft to 30,000 ft in depth with a vertical velocity gradient of , where ft/sec and . The shot position is located along the surface at a horizontal position of 40,000 ft from the left side of the model. The receiver array is from 25,000 ft to 55,000 ft on the surface. Figure 1 shows the computed wavefield from the 5 compact scatterers. The sampling rate is sec, and the peak frequency of the data is 5 Hz.
Figure 2 shows the partial image of the data from a single vertically propagated Gaussian beam with the planar wavefront at the surface. For simplicity, the source side Green’s function is constructed separately using Gaussian beams that are planar at the source location for all the examples given. Note on the figure that the images of the diffractors are curved and increase in width with depth. Figure 3 shows the partial image of the common shot data in Figure 1 using a single vertical Gaussian beam at a receiver position at 40,000 ft using a focused beam with the beam waist at a depth of 12,000 ft. Now the diffractor at a depth of 12,000 ft is the most focused using a single beam with the images of the other diffractors being broader and generally curved.
Figure 4 shows the partial image of the common shot data in Figure 1 using a single vertical Gaussian beam at a receiver position at 40,000 ft using dynamic focusing with the beam waist now specified at all depths. Now the partial image of the scatterers all have similar images with plane beam waists of the same beam widths for all scatterers. It would look as if we had been able to perform Gaussian beam migration with nondiffracting Gaussian beams. However, actually the image is formed with a number of different Gaussian beams for all scattering points in the medium such that in a dynamic fashion the beams have the same planar beam waists and beam widths at each scatterer.
Figure 5 shows the complete dynamically focused Gaussian beam image for the single-shot gather from Figure 1 using beams from all beam center locations launched at a range of angles. This results in focused images of all 5 scatterers and indicates that the imaging is being properly applied even with the shifted beam waists of the individual beam components for all imaging points of the medium. The small tails result from the limited angle range for the summation of the beams for the one-shot gather.
The dynamically-focusing beam approach is now applied to a single-shot gather from the Sigsbee2A data set distributed by SMAART (Subsalt Multiples Attenuation and Reduction Team) and available at http://www.delphi.tudelft.nl/SMAART/. In order to test the focused beam approach, a single-shot gather with a shot location at 6,325 ft from the left edge of the model is used. The receiver array starts at the shot location and has a maximum offset of 26,025 ft with a spacing of 75 ft. The background velocity model has the first layer from the surface down to the seafloor with a velocity of 5000 ft/sec. The second layer goes from seafloor to 30,000 ft in depth with a background velocity of , where ft/sec and . A salt dome exists with a velocity of 14,800 ft/sec in the middle and right parts of the model, but for the simple application here only one-shot gather away from the salt dome is used.
In Figure 6, the partial imaging for a single Gaussian beam is shown with the planar beam waist at the surface receiver depth. As in the earlier examples, the source side Green’s function is constructed separately using Gaussian beams that are planar at the source location. The receiver Gaussian beam has an initial location near the source and the partial image for a beam with a slight angle from the vertical is shown. As in the earlier example, when the beam waist is at the surface receiver depth, then curved beam fronts result which broaden with depth over the depth range of the model shown between 15,000 and 30,000 ft. This is typical of standard implementations of Gaussian beam migration. However, in regions of a complicated background medium, the medium itself can cause additional focusing of the beams.
Figure 7 shows the partial imaging results for a single Gaussian beam with the beam waist shifted to about 20,000 ft in depth. At this depth, the narrowest part of the beam image occurs. If the target structure were located at this depth, then fewer beams would be required to form a complete image. Also, the beam images would have planar beam fronts at this depth leading to more stable images. However, as shown in Figure 7 at other depths, the partial image results in curved and broader beam fronts.
Figure 8 shows the imaging results using dynamic focusing where all subsurface points have planar beams with the same beam width. This has the appearance of an image from a single nondiffracting beam, but since Gaussian beams generally diffract this is really a composite image for Gaussian beams with multiple imaging depths for a single central beam. Nonetheless, this type of imaging can be important with certain true amplitude formulations where planar and localized beams are required at the imaging points in the subsurface [20, 21]. These formulations originally involved launching beams directly from the scattering points in the subsurface up to the surface [18]. In the faster alternative proposed here, dynamically focused Gaussian beams are launched from the surface down to the imaging points at depth instead. However, the use of dynamically focused Gaussian beams is still slower than traditional Gaussian beam imaging or imaging using focused Gaussian beams with a single focusing depth. Precomputing dynamically focused slant stacks of the data can be implemented to improve the speed of the algorithm, although this would shift some of the computational burden into precomputing and storage of the slant stacks.
An additional area of concern when specifying the beam-widths at the imaging points is that if narrow beam-widths are specified then, a finer ray sampling would be required for adequate subsurface coverage of the beams. This strategy was used by Albertin et al. [22] when using windowed, planar Maslov propagators. However, in contrast to windowed Maslov propagators, Gaussian beams asymptotically satisfy the wave equation without further windowing. For dynamic focusing with narrow beam-widths at the imaging points, a finer ray sampling would be feasible when imaging more localized subsurface targets. Planar beams with larger beam-widths could also be specified in the dynamic focusing approach, and this would still provide for stable imaging results.
4. Conclusions
The application of dynamically focused Gaussian beams has been initially investigated for seismic imaging and is an extension of the focused beam approach of Nowack [1]. The shifting of the beam waists away from the source and receiver aperture adds flexibility to Gaussian beam algorithms allowing for the narrowest portions of the beams to occur at the depth of a specific target structure. This minimizes the number of beams required to form an image at the target depth. Also, at the beam waists the beam fronts are planar leading to more stable beam summations for imaging. To match with the surface data, quadratic phase corrections are required for the local slant-stacks of the data. Using dynamically focused Gaussian beams allows for beams to have planar wavefronts and the same beam widths at all scattering points in the subsurface. It has the appearance of using non-diffracting Gaussian beams at the scatterers, but since all Gaussian beams diffract, this is really the result of using beams with multiple focusing depths. The use of planar and localized Gaussian beams at all the subsurface scattering points also has advantages with regard to certain true amplitude imaging formulations. However, using dynamically focused Gaussian beams from the surface down to the scattering points avoids having to launch beams from every scattering point in the medium up to the surface. Nonetheless, dynamically focused Gaussian beam imaging is still slower than using either traditional Gaussian beam imaging or using a single focusing depth for imaging. Dynamically focused imaging was initially tested using a single-shot gather for a model with 5 scatterers at different depths, and then for a dynamically focused Gaussian beam for one-shot gather computed from the Sigsbee2A model.
Acknowledgments
This work was supported in part by the National Science Foundation EAR06-35611 and the Air Force Geophysics Laboratory contract FA8718-08-C-002 and partly by the sponsors of the Geo-Mathematical Imaging Group (GMIG) at Purdue University.