Asymptotic Ray Method for the Double Square Root Equation

: The parabolic wave equation describes wave propagation in a preferable direction, which is usually horizontal in underwater acoustics and vertical in seismic applications. For dense receiver arrays (receiver spacing is less than signal wavelength), this equation can be used for propagating the recorded wavefield back into the medium for imaging sources and scattering objects. Similarly, for multiple source and receiver array acquisition, typical for seismic exploration and potentially beneficial for ocean acoustics, one can model data in one run using an extension of the parabolic equation—the pseudo-differential Double Square Root (DSR) equation. This extended equation allows for the modeling and imaging of multi-source data but operates in higher-dimensional space (source, receiver coordinates, and time), which makes its numerical computation time-consuming. In this paper, we apply a faster ray method for solving the DSR equation. We develop algorithms for both kinematic and dynamic ray tracing applicable to either data modeling or true-amplitude recovery. Our results can be used per se or as a basis for the future development of more elaborated asymptotic techniques that provide accurate and computationally feasible results.


Introduction
The parabolic wave equation, also known as the one-way or Single Square Root equation, describes wave propagation in a preferable direction, which is usually horizontal in underwater acoustics [1] and vertical in seismic applications [2].There is another use for this equation that relates to data recorded by dense receiver arrays (receiver spacing is less than the signal wavelength): hydrophone arrays in water column or ocean bottom cables in seismic studies.Then, the parabolic wave equation can be used for propagating the recorded wavefield back from the array towards the source (backwards in time and the preferable direction).Such back-propagation can be used for the localization of active sources and scattering of objects in acoustics [3] or for the imaging of geological boundaries in seismic applications [4].Studying 3D heterogeneous models requires more data from multiple source and receiver arrays.Such acquisition geometry is widely used in seismic exploration [5] and may be beneficial for ocean acoustics [6], where it would necessitate a vertical receiver array with an active source firing at several depths along the same array.The data recorded in this way depend on both the source and the receiver coordinates, and one can back-propagate these data using the extended parabolic equation-the pseudodifferential Double Square Root (DSR) equation [2].
The kinematics of the DSR equation was first described in [7], and its pseudo-differential version was later introduced by John Claerbout [2].He used the DSR equation for the downward continuation of reflection-seismic data implementing the seismic imaging of subsurface geological structures ("survey-sinking" migration).This approach to seismic imaging is advantageous with respect to standard migration due to the simultaneous processing of all the data coming from different physical experiments, thus providing artifact-free images of the reflectors [8].Alternatively, the DSR equation can be used for upward continuation and modeling seismic data following the "exploding reflector" concept [9].Later on, a true-amplitude version of the DSR equation was derived from one-way wave equations [10][11][12], and its numerical implementation was discussed in [13,14].The DSR migration was generalized to a VTI anisotropic case [15] and applied as a tool for migration velocity analysis [16]-a technique for the simultaneous recovery of the velocity model and the scatterer/reflector position.However, till now, DSR equation-based methods have found limited popularity due to prohibitively high computational costs, as they operate in higher-dimensional space (source, receiver coordinates, and time), especially in 3D settings [9].
For this paper, we applied the ray theory to find a high-frequency asymptotic solution to the DSR equation.We started from the true-amplitude DSR equation from [11] and the high-frequency asymptotics of the pseudo-differential operators from [10].In our derivations, we followed Hamiltonian formalism, as presented in [17].We developed formulae for high-frequency up/down data continuation that can be used for seismic modeling and imaging.We demonstrated the validity of our results with a series of numerical experiments.We restricted our considerations to a 2D case, postponing generalization to a 3D case for future work.

DSR Equation-Asymptotic Solution
In [11], the authors introduced a true-amplitude version of the DSR equation: where the Greek capital letters denote pseudo-differential operators, with the principal symbols being Here and in the following, we shall use subscripts "s" and "r" for "source" and "receiver" coordinates, respectively, and slashes will indicate that either index can be chosen (self-consistently for the whole equation), meaning that is velocity at the point (x s/r , z) T .Superscript T stands for transposition.We denote angular frequency as ω, the source and receiver wavenumbers as k s/r , and the imaginary unit as i.Function u(x s , x r , z, t) can be interpreted as the single-reflected acoustic pressure field initiated at point source (x s , z) T and recorded at the receiver (x r , z) T .
Our derivation closely follows that from the first chapters in [17], which we shall further refer to as the "standard ray method", and we recommend seeing [18] for more detailed consideration.We omit some proofs identical to those in the references above.
We seek the solution to (1) in the form of an asymptotic ray series: where A k represents amplitude coefficients; τ represents two-way travel time, i.e., the time required for a wave to travel downward from the source to the reflecting interface and upward to the receiver array (different from the standard time variable in the wave equation); and functions F k denote the signal shape.The dynamics of the leading term of this ray series is characterized by the amplitude A 0 .Let us substitute solution (4) into the DSR Equation ( 1) and rearrange the coefficients corresponding to different F k .This leads to an infinite system of equations for τ and A k .The first equation controls the two-way travel time and will be referred to as the DSR eikonal equation: The second equation controls the principal amplitude coefficient A 0 and will be referred to as the DSR transport equation: where we use the shortened notations ∇ DSR τ and R DSR : Angle brackets ⟨•, •⟩ denote the scalar product of two vectors, and symbol ∇→ x stands for the gradient operator acting on both the source and receiver coordinates: The rigorous derivation of these equations is complicated, and we refer the reader to the Supplementary Materials for details.Here, we just note that in all our derivations, we assume non-horizontal wave propagation, i.e., all expressions under square roots are strictly positive.We also emphasize that the eikonal Equation ( 5) has been known for a long time [2,7,19], while the transport Equation ( 6) is our original result.

DSR Equation-Rays and Hamiltonians
Consider a 2D space with the Z-axis pointed downward.Consider two pointsx 0 s , z 0 T , x 0 r , z 0 T -and two nowhere horizontal curves connecting them with the source and receiver points (x s , z) T and (x r , z) T , respectively, so that z < z 0 .We denote these paths as γ s and γ r and introduce a two-way Fermat functional: where dγ s/r are length elements of the curves, and (ξ s/r , ζ) T denote coordinates along these curves.Let σ be a monotonically increasing parameter parametrizing both curves.The length elements mentioned above have the following form: We consider the paths γ s/r as one curve γ in a formal 3D extended space consisting of horizontal source and receiver coordinates and depth: Original two-dimensional paths become plane projections of a curve γ starting from x 0 s , x 0 r , z 0 T in this extended space, and the Fermat functional can be written as a single integral: where L is the Lagrangian.One can show that T(x s , x r , z; . ..) satisfies the DSR eikonal Equation ( 5) with respect to x s , x r , and z when evaluated on the extremal paths ∼ γ (which minimize or maximize the integral value).The proof is similar to that from the standard ray method in [17].We shall refer to these extremals as DSR rays hereafter.Following standard Hamiltonian formalism, we introduce the DSR equation slowness vector: Note the DSR slowness is not tangent to the DSR rays: ( The inverse relations to (14) depend on the choice of parameter σ in (10).In this paper, we consider two options: two-way travel time τ itself and negative depth −z.The negative sign is required to make the parameter monotonically increase along the rays.

Expressing
. → x as a function of → p using (15) or ( 16), we write down two Hamiltonians as and the Ray Tracing System in Hamiltonian form as: where ∇→ p is the gradient with respect to the slowness vector's components: . Explicit formulae for the right-hand side of this equation can be found in Appendix A.
By solving (19), one can trace extremals of the two-way Fermat functional (9) and then obtain the solution to the DSR eikonal Equation ( 5) by evaluating this functional on these extremals.However, proper initial conditions are required to obtain meaningful solutions corresponding to physical rays.We shall return to this question in a separate section.

DSR Equation-Dynamic Ray Tracing
Consider a DSR equation ray and corresponding slowness vector → p = ∇→ x τ.Let us substitute its components from (14) into the ∇ DSR τ vector (7): Note that ∇ DSR τ is tangent to the DSR ray, while the slowness vector ( 13) is not.We denote the scaling coefficient before .→ x as C DSR .This coefficient depends on the parametrization: Using (20), one can recast the DSR transport Equation ( 6) as an ODE along the chosen ray: and write down its formal solution: The last integral in this formula can be evaluated immediately as the DSR ray is traced.The first one requires a more elaborated technique.
Let us consider a family of DSR rays covering some volume in the extended sourcereceiver space (x s , x r , z).This family can be parametrized by three parameters: (θ 1 , θ 2 , σ).The first two determine a particular ray, and the third one defines the position on it.We introduce two Jacobian matrices: Following the standard notion, the square root of the determinant |det(Q)| can be interpreted as the geometrical spreading in the extended space.Using the same rationale as in [17], one can show that as soon as |det(Q)| ̸ = 0, the divergence in ( 23) is given by and the solution can be recast into This is the continuation formula for the DSR equation for ray amplitude.It closely resembles its conventional ray method counterpart [18], yet it includes an additional multiplier depending on the vertical inhomogeneity of the velocity model through R DSR .
Following the Hamiltonian approach, we present the Dynamic Ray Tracing System describing the evolution of these matrices along a ray: where the operators acting on the Hamiltonians are given by We provide explicit formulae for these matrices in Appendix A.

Initial Conditions for Ray Tracing and Dynamic Ray Tracing
Both the Ray Tracing System (19) and the Dynamic Ray Tracing System (27) require appropriate initial conditions to describe singly reflected waves.As in the wave equation ray method, one imposes Point Source initial conditions to model the respective wavefield and Wavefront initial conditions to extrapolate an existing one [18], we propose Exploding Reflector initial conditions for the DSR equation modeling and Survey Sinking initial conditions for data extrapolation.

Exploding Reflector Initial Conditions
The DSR equation can be used in seismic data modeling using the "exploding reflector" concept [9].To apply this approach to the DSR equation asymptotic solution, we demand that the DSR equation rays originate in one physical point on a boundary and obey Snell's law.
Consider a reflecting interface in 2D space specified explicitly as a twice-differentiable function f (x): Let us place a source and a receiver at a small distance above the reflector so that we can treat both incident and reflected rays as straight lines.Let (x 0 , z 0 ) T ≡ (x 0 , f (x 0 )) T be the reflection point, as in Figure 1.The colored vectors represent source and receiver projections of a DSR equation ray starting at (x 0 , x 0 , f (x 0 )) T , and their directions correspond to the desired ray tangent vector .
x r , .z .We demand that Snell's law is satisfied, i.e., that both vectors make the same angle with the reflector's normal.Using geometric reasoning, we obtain: where α is the reflection angle, and γ is the reflector's dip angle at the reflection point.The reflection angle α is counted from the reflector's normal, clockwise for the receiver and anticlockwise for the source branch.The dip angle γ is counted clockwise from the horizontal level, as in Figure 1.Initial conditions (30) define a ray coinciding with a physical one at any point  ,   on the reflector.One can recast them in terms of the slowness vector's components rather than the ray's tangent vector  ⃗ using (14).Tending both the source and receiver to the reflection point, we obtain the Exploding Reflector conditions for  0: or, in terms of derivative ′ ≡ rather than dips: Initial conditions (30) define a ray coinciding with a physical one at any point (x 0 , f (x 0 )) T on the reflector.One can recast them in terms of the slowness vector's com- ponents rather than the ray's tangent vector .→ x using (14).Tending both the source and receiver to the reflection point, we obtain the Exploding Reflector conditions for τ = 0: or, in terms of derivative f ′ ≡ d f dx rather than dips: The latter form explicitly depends on the reflection point's horizontal coordinate x 0 and the reflection angle α, allowing one to use them as the ray parameters θ 1 and θ 2 from (24).Hence, we write down the Exploding Reflector initial conditions for the Dynamic Ray Tracing System (27): where we once again use angle notation to simplify the expressions and introduce new multipliers, namely: The last columns in Q| τ=0 and P| τ=0 can be found from the Ray Tracing System (19) using the initial coordinates and slowness vector (31) or (32).

Survey Sinking Initial Conditions
The DSR equation was originally suggested for wavefield extrapolation [2,7].In this approach, the recorded data represent a boundary condition at the Earth's surface, with straightforward interpretation in terms of the DSR equation rays and slowness vectors.
Let us consider an idealized survey carried out at depth z = z 1 and covering all possible source-receiver pairs.
Assume that the recorded two-way travel time ) is a twice-differentiable function of the horizontal source's and receiver's coordinates.Then, for a fixed source-receiver pair x 1 s , x 1 r , the initial conditions for the Ray Tracing System (19) take the following form: where we use the DSR eikonal Equation ( 5) to obtain the vertical slowness p z .Survey plane coordinates x 1 s and x 1 r can be interpreted as ray parameters θ 1 and θ 2 from ( 24), and the initial conditions for the Dynamic Ray Tracing System (27) read as follows: Once again, the last columns are given by the Ray Tracing System (19), while the horizontal partial derivatives of p z read as follows: Despite the label of "initial", conditions (35) and (36) are more likely to be used as "final" conditions for the Ray Tracing and Dynamic Ray Tracing Systems (19) and (27).They allow one to trace DSR equation rays backwards in time or depth down to the reflectors, recovering their shape and reflectivities, through using measured travel time and its derivatives.

Asymptotic Ray Amplitudes and True-Amplitude Recovery
The initial conditions derived in the previous section are not sufficient for evaluating the ray amplitude using (26) since the geometrical spreading vanishes at the reflecting interface due to (33), and thus, the amplitude A 0 (τ 0 ) tends to infinity as τ 0 approaches zero.In this section, we derive expressions for ray amplitude computation and reflectivity recovery from the recorded data.

Asymptotic Ray Amplitudes
Consider the limit of (26) as σ 0 = τ 0 tends to zero: Let us evaluate the expression under limit assuming small values of τ 0 .Velocity variations can be neglected under this assumption, and the incident and reflected rays can be treated as straight lines, as in the Figure 1.We set A 0 (τ 0 ) to be equal to the standard ray method solution [18], which reads as follows: where v 0 is given by (34), R(x 0 , α) is the angle-dependent reflection coefficient, M is the source's magnitude, and L(τ 0 ; x 0 , α) is the relative geometrical spreading [18]: The Wolfram Language [20] was used to derive this formula from the general expressions provided in [18].On the other side, the exact formula for C DSR (τ 0 )|det(Q(τ 0 ))| from ( 27) and (33) reads as follows: Evidently, and the final expression for the DSR equation ray amplitude reads as follows: This amplitude meets the standard ray method in homogeneous media by construction, and we shall demonstrate its validity in inhomogeneous models in a series of numerical examples later.

True-Amplitude Recovery
The Survey-Sinking initial conditions (35) and (36) are designed to extrapolate the observed data downward to the reflectors.Thus, the observed amplitude A obs x 1 s , x 1 r can be substituted in place of A 0 (σ 0 = −z 1 ) in ( 26): We note that the target reflection coefficient is implicitly included in this expression through A obs , and according to (36), det(Q(−z 1 )) is equal to one.
Based on the previous subsection's result, we expect this amplitude to tend to the wave equation asymptotic amplitude in homogeneous medium (39) as the extrapolation depth approaches the reflector: Using L(z; x 0 , α), we denote the same relative geometrical spreading as in (40) but expressed in terms of depth rather than two-way travel time.Their relation can be easily obtained from Figure 1.With no loss in accuracy, we neglect the interface's curvature κ(x 0 ) and develop the latter limit: where we collect all nonsingular multipliers in a scaling coefficient C reg : Both the numerator and the denominator in (46) vanish at the boundary, and hence, we propose using L'Hôpital's rule to evaluate the limit.By differentiating both with respect to depth and substituting the result in (45), we obtain: We underline that |det(Q(−z))| is analogous to the geometrical spreading and that its derivative does not become zero when the back-propagated wavefront focuses on the reflector.This derivative can be found using matrix Q(−z)'s elements and the right-handside of the Dynamic Ray Tracing System (27).
The source's magnitude M is not involved in the relation (48), and a proportionality is used instead of exact equality.The reason is that the source's magnitude is unavailable in practice, and the reflectivity is recovered with accuracy up to an unknown scaling factor.

Results
In this section, we provide several examples supporting the validity of our formalism.We consider both ray modeling and true-amplitude recovery and use the finite-difference wave equation modeling from the Madagascar 3 package [21] as a reference.We implemented the DSR ray tracing engine in the Python 3.8 language using NumPy 1.21 [22] and SciPy 1.8 [23] libraries.For visualization routines, we used Matplotlib 3.5 [24] and Plotly 2.9 [25] packages.
We propose the following algorithm: 1.
Trace a family of DSR rays from the reflector to the surface using the Ray Tracing System (19) with initial conditions (31).

2.
Compute ray amplitudes using formula (43) and form synthetic seismograms using the technique presented in [18].

3.
Compare the DSR ray-synthetic seismograms with finite-difference wave equation gathers computed in the same model.

4.
Back-propagate the DSR rays from the survey plane down to the reflector using the Ray Tracing System (19) and initial conditions (35).

5.
Extract the reflectivities from the finite-difference data using formula (48) and compare them to the true ones.
In our experiments, we passed "true" travel time derivatives computed at the ray modeling step to the initial conditions (35), (36) when we back-propagated the wavefield.In a more realistic case, one ought to somehow estimate them from the recorded data.We shall briefly touch on this problem in the Discussion section below.
We consider three laterally varying velocity models given by and differing in the reflector's shape f (x).We consider a horizontal reflector, a dipping reflector, and a curvilinear reflector: In all three models, the overburden velocity is greater than that of the bedrock, which precludes critical incidents where the standard ray method fails.There are no caustic regions in these models, which also favors ray method's validity.All the models are shown in Figure 2. In our experiments, we passed "true" travel time derivatives computed at the ray modeling step to the initial conditions (35), (36) when we back-propagated the wavefield.In a more realistic case, one ought to somehow estimate them from the recorded data.We shall briefly touch on this problem in the Discussion section below.
We consider three laterally varying velocity models given by  , and differing in the reflector's shape   .We consider a horizontal reflector, a dipping reflector, and a curvilinear reflector: In all three models, the overburden velocity is greater than that of the bedrock, which precludes critical incidents where the standard ray method fails.There are no caustic regions in these models, which also favors ray method's validity.All the models are shown in Figure 2. We set up a survey consisting of 51 sources and 51 receivers ranging from  700 meters up to  700 meters with step of 28 meters at zero depth.To trace the rays, we used the shooting method, iteratively searching for an appropriate reflection point and angle to hit the desired source-receiver pair.In the lower row of Figure 3, we show examples of common-shot ray diagrams.Evidently, our rays satisfy Snell's law, and they are curved in the direction of the velocity gradient, which is in accordance with the standard ray method.We set up a survey consisting of 51 sources and 51 receivers ranging from x = −700 m up to x = 700 m with step of 28 m at zero depth.To trace the rays, we used the shooting method, iteratively searching for an appropriate reflection point and angle to hit the desired source-receiver pair.In the lower row of Figure 3, we show examples of common-shot ray diagrams.Evidently, our rays satisfy Snell's law, and they are curved in the direction of the velocity gradient, which is in accordance with the standard ray method.
Moreover, the DSR ray method's travel time and dynamics also match the wave equation modeling, as illustrated in the upper row of Figure 3.In this figure, two sets of wiggle traces are plotted parallelly.The darker traces are the DSR equation's ray method solutions, and the lighter traces are finite-difference solutions of the wave equation.The direct wave was muted before plotting.
The focusing of the back-propagated wavefield can be better illustrated in the 3D extended volume.In Figure 4, we present ray diagrams of the back-propagated wavefield.It can clearly be seen that the rays are focusing on the reflector x s = x r = x 0 , z = f (x 0 ) in all three cases.We applied our true-amplitude recovery formula (48) to the finite-difference data and obtained partial reflectivity maps of the reflectors.These maps consist of scaled reflection coefficients covering the illuminated part of the interfaces.We normalized them, making the reflectivity in the survey's center a unit, and compared them to their theoretical values.We present the results of true-amplitude recovery in Figure 5.We plotted the recovered reflectivities as 3D scatterplots superimposed on the true reflectivity maps in the "reflection point-reflection angle-depth" space.We outlined the illuminated area with red tri- We applied our true-amplitude recovery formula (48) to the finite-difference data and obtained partial reflectivity maps of the reflectors.These maps consist of scaled reflection coefficients covering the illuminated part of the interfaces.We normalized them, making the reflectivity in the survey's center a unit, and compared them to their theoretical values.We present the results of true-amplitude recovery in Figure 5.We plotted the recovered reflectivities as 3D scatterplots superimposed on the true reflectivity maps in the "reflection point-reflection angle-depth" space.We outlined the illuminated area with red triangles.Evidently, the recovered reflectivity accurately matches the true one in all three models.We applied our true-amplitude recovery formula (48) to the finite-difference data and obtained partial reflectivity maps of the reflectors.These maps consist of scaled reflection coefficients covering the illuminated part of the interfaces.We normalized them, making the reflectivity in the survey's center a unit, and compared them to their theoretical values.We present the results of true-amplitude recovery in Figure 5.We plotted the recovered reflectivities as 3D scatterplots superimposed on the true reflectivity maps in the "reflection point-reflection angle-depth" space.We outlined the illuminated area with red triangles.Evidently, the recovered reflectivity accurately matches the true one in all three models.Though our numeric tests above are but "toy" examples, we designed them to highlight the validity of the developed formalism only, and we do not intend to address fullscale seismic processing challenges at this moment.Furthermore, we do not have access to any DSR equation solver and thus cannot compare the performance directly.However, scaling our ray tracing engine to more realistic 2D applications, we reached a performance as high as 15,000 DSR rays per second on a modern PC.We emphasize that we relied on Though our numeric tests above are but "toy" examples, we designed them to light the validity of the developed formalism only, and we do not intend to address full-scale seismic processing challenges at this moment.Furthermore, we do not have access to any DSR equation solver and thus cannot compare the performance directly.However, scaling our ray tracing engine to more realistic 2D applications, we reached a performance as high as 15,000 DSR rays per second on a modern PC.We emphasize that we relied on a Python implementation, which is not the best choice from the perspective of performance, and significant improvements should be expected from more professional programs.

Discussion and Conclusions
In our paper, we applied an asymptotic ray method approach to the true-amplitude Double Square Root equation.First, we reproduced the DSR eikonal equation known from other publications.Then, we derived an original transport equation describing ray amplitudes of the DSR equation solutions.We applied Hamiltonian formalism to solve these two equations and derived formulae for ray tracing, data modeling, and inversion.On the one hand, our results give further validation to the idea that the considered true-amplitude DSR equation preserves the ray amplitudes of acoustic waves.On the other hand, they form the basis for the further implementation of amplitude-preserving imaging methods, Gaussianbeam construction, and velocity estimation techniques.These asymptotic methods have proven to be reliable tools, significantly reducing computational costs in the context of wavefield modeling and imaging [17,18].We expect similar results for the DSR equation applied to seismic or acoustic data acquired by source-receiver arrays.For instance, we plan to address the slope tomography [26] application of the ray perturbation theory for scattering object localization together with the velocity model update.
All the limitations of the standard ray method also apply to our DSR equation asymptotics.The derived formulae are inapplicable in the vicinity of caustics and critical-incidence rays.However, we think that uniform asymptotics developed for these regions in the standard ray theory should have their analogs for the DSR equation.Another limitation of the parabolic equation approach is that we assume no wave propagation in the direction orthogonal to the chosen direction, which may be horizontal or vertical for ocean acoustics, sediment layer studies, or seismic studies.We believe that this assumption may be satisfied in moderately inhomogeneous media.In addition, this difficulty can be mitigated by introducing local curvilinear coordinates where the waves will not turn to the undesired direction [27].
Note that for our asymptotic formulas, we need travel times and their derivatives as input for data back-propagation and true-amplitude recovery.Wave arrival travel times and their slopes are routinely estimated in the slope tomography method, and several techniques have been developed to extract them from seismic data [26,28,29].Estimating the travel time second derivatives will require some extra work, but theoretically, they can be found by the slant stacking technique described in [30].(A3) We also consider that the Hamiltonians zero out along the DSR equation rays, and thus we omit them whenever they appear during differentiation by the chain rule.

19 Figure 1 .
Figure 1.Ray diagram of a reflection event.

Figure 1 .
Figure 1.Ray diagram of a reflection event.

Figure 2 .
Figure 2. Velocity models.The black lines denote the reflectors' position.

Figure 2 .
Figure 2. Velocity models.The black lines denote the reflectors' position.

Figure 3 .
Figure 3. Upper row: common-shot gathers.Darker lines correspond to the DSR equation ray modeling; lighter lines correspond to finite-difference wave equation modeling.Lower row: commonshot ray diagrams.Moreover, the DSR ray method's travel time and dynamics also match the wave equation modeling, as illustrated in the upper row of Figure3.In this figure, two sets of wiggle traces are plotted parallelly.The darker traces are the DSR equation's ray method solutions, and the lighter traces are finite-difference solutions of the wave equation.The direct wave was muted before plotting.The focusing of the back-propagated wavefield can be better illustrated in the 3D extended volume.In Figure4, we present ray diagrams of the back-propagated wavefield.It can clearly be seen that the rays are focusing on the reflector    ,    in all three cases.

Figure 3 .
Figure 3. Upper row: common-shot gathers.Darker lines correspond to the DSR equation ray modeling; lighter lines correspond to finite-difference wave equation modeling.Lower row: commonshot ray diagrams.

Figure 4 .
Figure 4. Ray diagrams in the extended volume.Colored rectangles represent survey points and recorded amplitudes.The red lines mark the true interface position, the blue lines mark the DSR equation rays, and the black circles represent the reflection points.

Figure 4 .
Figure 4. Ray diagrams in the extended volume.Colored rectangles represent survey points and recorded amplitudes.The red lines mark the true interface position, the blue lines mark the DSR equation rays, and the black circles represent the reflection points.

Figure 4 .
Figure 4. Ray diagrams in the extended volume.Colored rectangles represent survey points and recorded amplitudes.The red lines mark the true interface position, the blue lines mark the DSR equation rays, and the black circles represent the reflection points.

Figure 5 .
Figure 5. Reflectivity maps of the reflectors.The red triangles outline the scatterplots of the recovered reflectivity, while the background surfaces represent true reflectivity.

Figure 5 .
Figure 5. Reflectivity maps of the reflectors.The red triangles outline the scatterplots of the recovered reflectivity, while the background surfaces represent true reflectivity.