# Dynamically Focused Gaussian Beams for Seismic Imaging

Authors: Robert L. Nowack
Journal: International Journal of Geophysics
 (2011)

---

## 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.


---

## Body

## 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δm(x̲). Thus, (1)δm(x̲)~∫dx̲sIs(x̲),
where the adjoint image for each shot point can be written as
(2)Is(x̲)=∫dω2πB(ω)∫dx̲gg̅(x̲,x̲s,ω)g̅(x̲,x̲g,ω)us(x̲g,x̲s,ω)B(ω)=ω2S̅(ω),
where ω is the radial frequency, S(ω) is the source time function, x̲g is the geophone location, x̲s is the source location, and us(x̲g,x̲s,ω) is the scattered field recorded at the geophones. g(x̲,x̲g,ω) and g(x̲,x̲s,ω) are the Green’s functions, from the geophones and source to the subsurface scattering points x̲. From reciprocity of the Green’s function g(x̲,x̲′,ω)=g(x̲′,x̲,ω) and g̅ indicates the complex conjugate of g.In the 2D case, the Green’s function can be written in terms of a summation of Gaussian beams as(3)g(x̲,x̲′,ω)=-i4π(εv0)1/2∫dγugb(x̲,x̲′,ω)=-i4π(εv0)1/2∫dp1rp3rugb(x̲,x̲′,p̲r,ω),
where ε is the complex beam parameter [3] and(4)ugb(x̲,x̲′,ω)=[v(s)q(s)]1/2exp⁡{iωτ(s)+iω2M(s)n2}.The coordinates(s,n) correspond to the subsurface position x̲ in ray-centered coordinates, τ(s) is the travel time along the central ray, v(s) is the velocity along the central ray, and the horizontal component of the ray parameter vector at the initial point is p1r=sinγ/v0, where γ is the take-off angle from the vertical.The complex second derivative of the travel time field with respect to the transverse coordinaten can be written as(5)M(s)=MR(s)+iMI(s)=p(s)q(s),
where MR(s) is related to the wavefront curvature K(s) of the beam by MR(s)=K(s)/v(s).  To form a bounded Gaussian beam, then MI(s)>0. The variables p(s) and q(s) 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
(6)π(s,s0)=[q1(s)q2(s)p1(s)p2(s)],
where π(s0,s0)=I. q1(s), p1(s) are solutions for an initial plane wave and q2(s), and p2(s) are for an initial line source. The inverse of π(s,s0) is(7)π-1(s,s0)=π(s0,s)=[p2(s)-q2(s)-p1(s)q1(s)].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](8)q(s)=εq1(s)+q2(s),p(s)=εp1(s)+p2(s),
where ε is the beam parameter. The variable q(s) is related to the complex geometrical spreading along the beam. Since for the fundamental solution matrix Det(π(s,s0))=1 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 ton can be written as (9)M(s)=MR(s)+iMI(s)=εp1(s)+p2(s)εq1(s)+q2(s).Since at the source pointπ(s0,s0)=I and ε=εR-iεI, then (10)M(s0)=MR(s0)+iMI(s0)=1ε=εRε*ε+iεIε*ε,
where ε*ε is the magnitude squared of ε. Alternatively, ε can be written in terms of M(s0) as(11)ε=MR(s0)M*(s0)M(s0)-iMI(s0)M*(s0)M(s0),
where M*(s0)M(s0) is the magnitude squared of M(s0). Thus, the complex beam parameter can be specified either directly in terms of ε or in terms of M(s0).The complex beam parameter can also be written(12)ε=εr-iεi=v0S0-iv0L02,
where v0 is an initial velocity. In a homogeneous medium, S0 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(13)q(s)=q1(s)+M(s0)q2(s),p(s)=p1(s)+M(s0)p2(s),
where(14)M(s0)=MR(s0)+iMI(s0)=1ε.In either case, we can obtainM(s) at the scattering point, from the initial value of M(s0) along the source and receiver array by solving the dynamic ray equations for p(s) and q(s), and then forming M(s)=p(s)/q(s). Alternatively, we can specify M(s) at the scattering point and then use the solution matrix for the dynamic ray equation π(s0,s) from s to s0 to obtain p(s0) and q(s0) at the initial point of the ray. From these, one can obtain the corresponding value of beam curvature M(s0) at the beginning point of the ray. At a point along the beam, M(s) can be written as(15)M(s)=MR(s)+iMI(s)=K(s)v(s)+iMI(s),
where K(s) is the wavefront curvature in ray coordinates, v(s) is the velocity along the ray, and MI(s) 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(16)exp⁡⁡{-ω2MI(s)n2}=exp⁡⁡{-n22L2(s)},
where L(s)=(ωMI(s))-1/2 is the beam half-width transverse to the ray. At the initial point of the beam, L(s0)=(ωMI(s0))-1/2=(ε*ε/ωεI)1/2. For the case of the beam waist at the initial point of the beam MR(s0)=0, and L(s0)=(v0/ω)1/2L0 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 MR(s0)=0 where s≠s0, the beam front is planar at some point s along the ray and generally curved at the initial position s=s0. 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 x coordinate can be written as
(17)exp⁡⁡{iω2Kx(s0)v0(x-xL)2}exp⁡⁡{-(x-xL)22Lx-ref2(s)},
where Kx(s0) is the initial real-valued horizontal beam curvature and Lx-ref(s0) is the initial horizontal beam half-width at the reference frequency ωref. To match this with the initial parameters of the beam propagated into the medium, then the transverse coordinate of the beam n=cosγ(x-xL) where γ is the angle of the beam with respect to the vertical. Given the initial values Kx(s0) and Lx-ref(s0) along the source and receiver aperture, then the initial values for MR(s0) and MI(s0) for the beams are (18)MR(s0)=Kx(s0)v0cos⁡2⁡(γ),MI(s0)=(ωrefcos⁡2⁡(γ)Lx-ref2(s0))-1,
and then the initial beam parameter ε=εr-iεi 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(19)1~12πΔLσ∑m=-∞∞e-(x-mΔL)2/2σ2,
where (20)ΔL≪2σ.Assuming a regularly spaced set of beam centersxL=mΔL 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 p1g(x1g-mΔL)+Kx(s0)(x1g-mΔL)2/2v0. Then,
(21)g(x̲,x̲g,ω)~Cg∫dp1gp3gugb(x̲,x̲L,p̲g,ω)eiω(p1g(x1g-mΔL)+Kx(s0)(x1g-mΔL)2/2v0),
where Cg represents the terms in front of the last integral in (3) for the receiver Green’s functions and p1g is the horizontal component of the ray parameter vector along the receiver aperture. The imaging formula can then be written as(22)Is(x̲)=∑m=-∞∞∫dω2πA1(ω)∫dx̲g×∫dp1gp3gu̅gb(x̲,x̲L,p̲g,ω)×e-iω(p1g(x1g-mΔL)+Kx(s0)(x1g-mΔL)2/2v0)×g̅(x̲,x̲s,ω)us(x̲r,x̲s,ω)e-(x1g-mΔL)2/2σ2,
where (23)A1(ω)=(+i2π)(ε̅v0g)1/2ω2S̅(ω)12πΔLσ.The source Green’s function now needs to be decomposed into Gaussian beams, but for simplicity here will be just referred to asg̅(x̲,x̲s,ω), where the over-bar signifies the complex conjugateThe common-shot imaging formula for nonplanar, focused beams along the aperture plane can then be written as(24)Is(x̲)=∑m=-∞∞∫dω2π∫dp1gp3gA1(ω)g̅(x̲,x̲s,ω)×u̅gb(x̲,x̲L,p̲g,ω)Dp(x̲L,x̲s,p̲g,ω),
where (25)Dp(x̲L,x̲s,p̲g,ω)=∫dx1gus(x̲g,x̲s,ω)e-(x1g-mΔL)2/2σ2×e-iω(p1g(x1g-mΔL)+Kx(s0)(x1g-mΔL)2/2v0).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 (26)x1L=mΔL,ΔL≪2σ,
and the initial horizontal beam widths are
(27)σ=Lx-ref(s0)(ωrω)1/2,
where Lx-ref(s0) is the horizontal half-width of an initial Gaussian function at the reference frequency ωref. 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 ofv(z)=v0+k(z-zb), where v0=5000 ft/sec and k=.15. 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 .008 sec, and the peak frequency of the data is 5 Hz.Figure 1
Computed common-shot data is shown for a shot point on the surface at a position of 40,000 ft. The receiver array is from 25,000 ft to 55,000 ft. Diffractions from 5 compact scatterers are shown each with a horizontal position of 40,000 ft and depths of 8,000, 12,000, 16,000, 20,000, and 24,000 ft.Figure2 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 2
The partial image of the common shot data in Figure1 using a single vertical Gaussian beam at a receiver position at 40,000 ft using a beam with a planar beam waist at the surface.Figure 3
The partial image of the common-shot data in Figure1 using a single vertical Gaussian beam at a receiver position at 40,000 ft using a focused beam with a beam waist at a depth of about 12,000 ft.Figure4 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 4
The partial image of the common shot data in Figure1 using a vertical dynamically focused beam.Figure5 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.Figure 5
The Gaussian beam image is shown using the single-shot gather in Figure1 with beams from all beam positions launched at a range of angles, each with dynamically focused beams. This results in focused images of all 5 scatterers and indicates that the imaging is being properly applied even with dynamically focused individual beam components.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 athttp://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 v(z)=v0+k(z-zseafloor), where v0=5000 ft/sec and k=.30. 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 Figure6, 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 6
This shows the partial imaging result for a single Gaussian beam with the beam waist at the surface along the receiver array for one-shot gather from the Sigsbee2A dataset.Figure7 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 7
This shows the partial imaging result for a single Gaussian beam with the beam waist shifted to about 20,000 ft in depth for a single-shot gather from the Sigsbee2A dataset.Figure8 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.Figure 8
This shows the partial imaging result for a single dynamically focused beam for one-shot gather from the Sigsbee2A dataset.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.



---

*Source: 316581-2011-05-26.xml*