TY - JOUR
AU - Symes, William W
AB - SUMMARY Seismic migration in the angle-domain generates multiple images of the earth's interior in which reflection takes place at different scattering-angles. Mechanically, the angle-dependent reflection is restricted to happen instantaneously and at a fixed point in space: Incident wave hits a discontinuity in the subsurface media and instantly generates a scattered wave at the same common point of interaction. Alternatively, the angle-domain image may be associated with space-shift (regarded as subsurface offset) extended migration that artificially splits the reflection geometry. Meaning that, incident and scattered waves interact at some offset distance. The geometric differences between the two approaches amount to a contradictory angle-domain behaviour, and unlike kinematic description. We present a phase space depiction of migration methods extended by the peculiar subsurface offset split and stress its profound dissimilarity. In spite of being in radical contradiction with the general physics, the subsurface offset reveals a link to some valuable angle-domain quantities, via post-migration transformations. The angle quantities are indicated by the direction normal to the subsurface offset extended image. They specifically define the local dip and scattering angles if the velocity at the split reflection coordinates is the same for incident and scattered wave pairs. Otherwise, the reflector normal is not a bisector of the opening angle, but of the corresponding slowness vectors. This evidence, together with the distinguished geometry configuration, fundamentally differentiates the angle-domain decomposition based on the subsurface offset split from the conventional decomposition at a common reflection point. An asymptotic simulation of angle-domain moveout curves in layered media exposes the notion of split versus common reflection point geometry. Traveltime inversion methods that involve the subsurface offset extended migration must accommodate the split geometry in the inversion scheme for a robust and successful convergence at the optimal velocity model. Image processing, Numerical approximations and analysis, Numerical solutions, Tomography, Wave propagation, Wave scattering and diffraction INTRODUCTION An image of the earth's structure is traditionally produced by pre-stack depth migration applied on surface seismic data. The lack of an accurate velocity model makes the imaging process kinematically wrong. Meaning that, the migration of data space reflections misplaces model space imaged reflectors, and results with poor quality. Hence, pre-stack migration itself is being considered as a tool for velocity model optimization via techniques regarded as migration velocity analysis (MVA; Liu & Bleistein 1995; Sava & Biondi 2004; Shen & Symes 2008). The velocity optimization process is cast iteratively by generating a loop between migration, MVA, and velocity updates. To estimate the velocity update at each iteration, MVA requires to extend the model space by one or more redundant parameters, and to compute an extended image in the form of common-image gathers (CIGs). The kinematic compatibility between data and velocity model is measured in the CIGs by some semblance conditions that maximize when velocity is chosen correctly (Symes 2008). The most common condition is flatness, applied with respect to various types of CIGs, such as acquisition offset, scattering-angle or shot index (Al-Yahya 1989; Xu et al. 2001; Huang et al. 2016). The more flat the CIGs, the more accurate the velocity model. Another useful semblance condition is focusing. Space-shift or time-shift image extensions (Sava & Fomel 2003, Sava & Fomel 2006) are usually associated with such condition that indicates the true velocity when the image is focused at the zero-lag trace. The concept of MVA has been proposed as a picking-free application, by Symes & Carazzone (1991) in the framework of differential semblance optimization, and was found useful as a regularization term in full waveform inversion schemes (Symes 2008; Biondi & Almomin 2014). MVA is also popular as a picking-based method. Hand- or automatic-made picking of wrongly migrated reflections in the CIGs is followed by a traveltime tomography inversion to back project the imaging errors into the velocity field and obtain a compatible velocity update (Kosloff et al. 1996). Either way, parametrizing and adequately predicting the moveout (or defocusing) of reflections in CIGs is a common objective of many studies. Overall, it is usually accomplished by relating the geometry of reflection with the physical model properties (i.e. velocity). The objective of this work is to study CIGs, computed by space-shift extended migration, and their angle-domain transformation (Rickett & Sava 2002; Sava & Fomel 2003). We regard the space-shift as subsurface offset to emphasize its geometrical role. It turns out this role is what differentiates the subsurface offset CIGs from others. The subsurface offset extended migration has a peculiar kinematic description, where incident and scattered waves interact at a distance. It is governed by a non-physical split reflection geometry, and contradicts some basic continuum mechanics principles relating strain to stress only at a common point in space. Therefore, the kinematics of subsurface offset domain imaging is obviously not intuitive and waiting to be explored. What makes it even more confusing is the angle-domain transformation of the subsurface offset CIGs (Dafni & Symes 2016a). Describing the direction of scattering, based on a split, rather than common, reflection geometry, is somewhat contradictory. As noted by Bartana et al. (2006) and Montel & Lambare (2013), angle-domain CIGs that are derived from the subsurface offset extended image shows a distinctive moveout shape, different from other angle-domain imaging methods. Our work is inspired by these two preliminary studies. We provide the reasoning behind that moveout dissimilarity and further claim that considering the subsurface offset to angle domain transform as a true representation of the reflection's half-opening angle is misleading under some specific, but general, conditions. However, such transformation is still usable for velocity optimization. In the following, we describe the kinematics of two variations of angle-domain imaging techniques. The first is regarded as conventional since it decomposes the angles at common reflection points across the imaged reflector. The second is associated with the subsurface offset extended migration, where the split reflection geometry is underlay. We highlight the fundamental differences between the two descriptions and how it amounts to a dissimilar angle-domain moveout behaviour. A semi-analytic simulation is derived accordingly, to support the observations in the CIGs, by predicting asymptotic angle-domain moveout curves in layered media. The simulation produces as a byproduct the defocusing pattern of the subsurface offset extended image as well. This work strongly asserts that, although being non-physical by definition, the split reflection geometry of subsurface offset domain imaging is as valuable for MVA purposes, as the conventional common reflection geometry. EXTENDED MIGRATION IN THE ANGLE-DOMAIN Migration operators are generally regarded as the adjoint of Born-type modelling operators. The output of the adjoint operator gives an image of the subsurface I(x), and has the integral representation   \begin{eqnarray} I({{\bf x}}) &=& \int{{{\rm{d}}{{{\bf x}}_{{\bf r}}}\int{{{\rm{d}}{{{\bf x}}_{{\bf s}}}\int{{{\rm{d}}t\frac{{{\partial ^2}}}{{\partial {t^2}}}D({{{\bf x}}_{{\bf r}}},t;{{{\bf x}}_{{\bf s}}})}}}}}}\nonumber\\ &&\times\,\int{{{\rm{d}}\tau G ({{\bf x}},t - \tau ;{{{\bf x}}_{{\bf r}}})G({{\bf x}},\tau ;{{{\bf x}}_{{\bf s}}})}}, \end{eqnarray} (1)where G(x, t) is Green's function, D(xr, t; xs) stands for the seismic data subjected to shot-receiver pairs (xs, xr), and τ is the migration time. Extended migration operators have similar adjoint formulation, after extending the definition of reflectivity to depend on more degrees of freedom in the modelling formula. The output of the extended migration operators is the pre-stack image. It depends on the same extra degrees of freedom as the extended reflectivity in the Born modelling construction (Symes 2008). The pre-stack image is usually sorted into CIGs, allowing migration velocity and amplitude analysis to exploit the redundancy in the seismic data. In the following we introduce two versions of extended migration operators that compute angle-domain pre-stack images. For simplicity, we restrict our study to 2-D only. Explicit angle-domain extension An angle-domain extended migration is explicitly derived via an angular decomposition of the wavefield. The pre-stack image I(x, γ), extended by the scattering-angle γ, is formulated as   \begin{eqnarray} I({{\bf x}},\gamma ) &=& \int{{{\rm{d}}{{{\bf x}}_{{\bf r}}}\int{{{\rm{d}}{{{\bf x}}_{{\bf s}}}\int{{{\rm{d}}t\frac{{{\partial ^2}}}{{\partial {t^2}}}D({{{\bf x}}_{{\bf r}}},t;{{{\bf x}}_{{\bf s}}})}}}}}}\nonumber\\ &&\times\,\!\int\!\!{{{\rm{d}}\tau \tilde{G} }}({{\bf x}},{\xi _r} + \tan \gamma ,t - \tau ;{{{\bf x}}_{{\bf r}}})\tilde{G}({{\bf x}},{\xi _s} - \tan \gamma ,\tau ;{{{\bf x}}_{{\bf s}}}).\nonumber\\ \end{eqnarray} (2) In the inner integrand, the incident and scattered wavefields (represented by the Green's functions) are decomposed into the angle-domain at each migration time step with a local Radon transform (indicated by the tilde symbol)   \begin{equation}\tilde{G}({{\bf x}},\xi ,t) = \frac{1}{{\Delta x}}\int\limits_{{ - \Delta x}}^{{\Delta x}}{{G(z + \xi x',x + x',t){\rm{d}}x'}}.\end{equation} (3) The transform's slope is defined as: ξ = – ∂z/∂x and indicates the wavefield's local angle of propagation. The additional variable x΄ represents the locality of the Radon transform around the x-coordinate, and Δx stands for the effective range where x΄ is sampled. The scattering-angle is defined as the half-opening angle between the incident and scattered waves. Therefore, the angle-dependent image extension is obtained by shifting the Radon transformed Green's functions accordingly   \begin{equation}\tan \gamma = \frac{{{\xi _r} - {\xi _s}}}{2}.\end{equation} (4) Several variations of this angle-domain extended migration have been proposed over the years based on wave-equation methods (de Bruin et al. 1990; Prucha et al. 1999; Mosher & Foster 2000; Wu & Chen 2006; Yoon & Marfurt 2006) or ray-theory implementations (Koren & Ravve 2011; Ravve & Koren 2011). They all decompose the angles explicitly from the propagating wavefields before the imaging condition in applied. Furthermore, their imaging condition is purely physical: Interaction between incident and scattered wavefields at the same migration time and at a common point in space (i.e. free of temporal or spacial extensions). Hence, the kinematic description of the explicit angle-domain extended migration is governed by a common-reflection point (CRP) geometry. Implicit angle-domain decomposition via subsurface offset extension One common choice of image extension is the subsurface offset, defined as the space-shift between a sunken shot and receiver pair in the subsurface. We restrict here the subsurface offset to the horizontal direction only, as proposed by Claerbout (1985) in the framework of survey-sinking migration. The pre-stack image I(x, h), extended by the horizontal subsurface half-offset h, is formulated as   \begin{eqnarray} I({{\bf x}},h) &=& \int{{{\rm{d}}{{{\bf x}}_{{\bf r}}}\int{{{\rm{d}}{{{\bf x}}_{{\bf s}}}\int{{{\rm{d}}t\frac{{{\partial ^2}}}{{\partial {t^2}}}D({{{\bf x}}_{{\bf r}}},t;{{{\bf x}}_{{\bf s}}})}}}}}}\nonumber\\ &&\times\,\int{{{\rm{d}}\tau G }}({{\bf x}} + h,t - \tau ;{{{\bf x}}_{{\bf r}}})G({{\bf x}} - h,\tau ;{{{\bf x}}_{{\bf s}}}). \end{eqnarray} (5) It describes an action with an offset between the incident and scattered wavefields (Stolk et al. 2009). Any offset distance (or space-shift) other than zero is considered as non-physical since applying stress at one point cannot cause strain at a distance. Nevertheless, the non-zero subsurface offsets are essential to capsulize the traveltime and amplitude variants of the seismic data. The non-physical nature of the subsurface offset extension implies that a unique kinematic description is involved. Reflection takes place with an offset, rather than at a common-point in space. We categorize this non-physical reflection geometry as a split-reflection point (SRP) setting. The subsurface offset extended image is transformed into the angle-domain with a conventional Radon transform, defined by the transform's slope p = –∂z/∂h  \begin{equation}I({{\bf x}},\gamma ) = \int{{I(z + ph,x,h){\rm{d}}h}},\end{equation} (6)where the equality p = tan γ is assumed. Hence, the angle-dependent pre-stack image results from a post-migration transformation of the subsurface offset extended image. The angles are evaluated post the imaging condition and do not relate directly to the incident and scattered wavefields. The relations between the subsurface offset and the angle-domain have been studied extensively by many authors (Rickett & Sava 2002; Sava & Fomel 2003; Biondi & Symes 2004; Sava & Vasconcelos 2011). One fundamental aspect, which is still controversial, is whether the transform's slope in eq. (6) truly defines the scattering-angle γ as the half-opening angle between incident and scattered waves. We address and clarify this confusing issue latter in this paper. KINEMATIC DESCRIPTION OF ANGLE-DOMAIN IMAGING Both of the methods provided in the preceding section yield the pre-stack image with respect to the same angle-domain quantity, namely the scattering-angle of reflection. However, the latter of the two methods do so implicitly, as it involves an extended imaging condition in space coordinates (i.e. subsurface offset). Moreover, it is still arguable if this implicit formulation truly derives the scattering-angle as the half-opening angle of reflection. Due to the spatial extension, the kinematic description of the two methods is fundamentally different and amounts to two different geometry settings we consider as CRP and SRP settings. We use high-frequency asymptotics of wave propagation to depict these kinematic differences and describe the geometry of reflections and imaged reflectors. For the purpose of this paper, we assume that all traveltimes between shot or receiver and reflecting point are singled-valued, that is, that no caustics occur in the scattering domain. Kinematics of CRP reflections For a particular shot-receiver pair (xs, xr), the two-way traveltime isochron T of a CRP reflection is defined in terms of the one-way traveltimes tr and ts as   \begin{equation}T(z,x;{{{\bf x}}_{{\bf r}}},{{{\bf x}}_{{\bf s}}}) = {t_r}(z,x;{{{\bf x}}_{{\bf r}}}) + {t_s}(z,x;{{{\bf x}}_{{\bf s}}}).\end{equation} (7) Migration relates the phase space coordinates of reflection events in the seismic data with those of the imaged reflectors (Stolk et al., 2009). Asymptotically, a reflection event in the data volume is migrated to the image volume along a plane-wave component (or ray) in the acoustic field. This plane-wave in phase space amounts to the two-way traveltime gradient (i.e. slowness vector)   \begin{equation}\nabla {{\bf T}}(z,x;{{{\bf x}}_{{\bf r}}},{{{\bf x}}_{{\bf s}}}) = ({\nabla _z}T,{\nabla _x}T).\end{equation} (8) At the reflection's two-way traveltime, an image of a reflector would be constructed if the propagating plane-wave and the reflector match in their phase space co-ordinates. This establishes Snell's law of reflection: The reflecting wave front is locally parallel to the reflector interface, and the bisector of the angle subtended by incident and scattered rays is the reflector normal. We characterize the normal to the reflector in phase space with the wavenumbers   \begin{equation}{{\bf n}} = ({k_z},{k_x}).\end{equation} (9) The traveltime gradient ∇T obeys Snell's law with the reflector normal n when their cross product vanishes   \begin{equation}\nabla {{\bf T}} \times {{\bf n}} = 0.\end{equation} (10) From which it follows that   \begin{equation}{\nabla _z}T{k_x} = {\nabla _x}T{k_z}.\end{equation} (11) Accordingly, the angle subtended from the vertical axis defines the reflector's local slope q and indicates its normal direction   \begin{equation}- \frac{{{\nabla _x}T}}{{{\nabla _z}T}} = - \frac{{{k_x}}}{{{k_z}}} = \frac{{\partial z(x)}}{{\partial x}} = q.\end{equation} (12) Hence, the reflector's image is a consequence of CRP reflections oriented in the direction of the reflector normal q. In Fig. 1(a) we plot an image of a –5° dipping up reflector. The image was computed by using 10 per cent too-slow migration velocity. Therefore, the reflector is misplaced in the image section as indicated by the blue cross that marks the reflector's true depth. The red arrow indicates the normal to the imaged reflector, whereas the angle subtended from the vertical is the angle q. The green curve represents the two-way traveltime isochron. It is tangent to the imaged reflector at the reflection point. The image was constructed by CRP reflection events as illustrated on the left side of the figure. The green arrows mark an incident and scattered ray-pair, interacting at a common-point in space. Figure 1. View largeDownload slide (a) A conventional and (b) a subsurface offset extended image of a –5° dipping up reflector, migrated by using 10 per cent too-slow velocity. The traveltime isochron (green) is perpendicular to the reflector normal at the reflection point, as indicated by the angles q and p. The isochron represents incident and scattered ray pairs that interact at a (a) common or (b) split reflection point, as illustrated on the left. Figure 1. View largeDownload slide (a) A conventional and (b) a subsurface offset extended image of a –5° dipping up reflector, migrated by using 10 per cent too-slow velocity. The traveltime isochron (green) is perpendicular to the reflector normal at the reflection point, as indicated by the angles q and p. The isochron represents incident and scattered ray pairs that interact at a (a) common or (b) split reflection point, as illustrated on the left. Eq. (12) reveals angular information about the constructive part of the two-way traveltime isochron. However, it does not give rise to the hitting direction of the one-way incident and scattered rays at the CRP location. To capture this information and compute the angle-dependent image, one must decompose the incident and scattered wave fronts into the angle-domain before they interact to construct the reflector's image, as suggested by eqs (2) and (3). Kinematics of SRP reflections Similarly to the CRP description, the mutual relations in phase space between reflections and imaged reflectors establish the basis for the kinematic depiction of the SRP reflection. However, the subsurface offset that splits the reflection point between the incident and scattered rays must be taken into account. Therefore, we extend the definition of the two-way traveltime isochron in eq. (7) with respect to the subsurface offset split   \begin{equation}T(z,x,h;{{{\bf x}}_{{\bf r}}},{{{\bf x}}_{{\bf s}}}) = {t_r}(z,x + h;{{{\bf x}}_{{\bf r}}}) + {t_s}(z,x - h;{{{\bf x}}_{{\bf s}}}),\end{equation} (13)and derive an extension for the traveltime gradient   \begin{equation}\nabla {{\bf T}}(z,x,h;{{{\bf x}}_{{\bf r}}},{{{\bf x}}_{{\bf s}}}) = ({\nabla _z}T,{\nabla _x}T,{\nabla _h}T).\end{equation} (14) We also extend eq. (9) reflector normal wavevector in accordance   \begin{equation}{{\bf n}} = ({k_z},{k_x},{k_h}).\end{equation} (15) Snell's law is extended respectively by substituting eqs (14) and (15) into eq. (10). It assures that the extended reflecting wave front and the extended reflector interface are locally parallel. From which the following relations hold:   \begin{eqnarray} {\nabla _z}T{k_x} = {\nabla _x}T{k_z}\nonumber\\ {\nabla _h}T{k_z} = {\nabla _z}T{k_h}{\rm{ }}{\rm{.}} \end{eqnarray} (16) Two angles, subtended from the vertical axis, are defined by eq. (16) to describe the normal direction of the SRP extended reflector   \begin{eqnarray} - \frac{{{\nabla _x}T}}{{{\nabla _z}T}} &=& - \frac{{{k_x}}}{{{k_z}}} = \frac{{\partial z(x)}}{{\partial x}} = q \nonumber\\ - \frac{{{\nabla _h}T}}{{{\nabla _z}T}} &=& - \frac{{{k_h}}}{{{k_z}}} = \frac{{\partial z(h)}}{{\partial h}} = p{\rm{ }}{\rm{.}} \end{eqnarray} (17) The first angle, indicated by the ratio q, has already been introduced in eq. (12). It is defined in the z, x subspace in the same way as in the CRP description. The second angle is defined in the z, h subspace and is denoted by the ratio p. Hence, the extended reflector's image is a consequence of extended SRP reflections oriented in the normal direction, given by the angles q and p. In Fig. 1(b) we plot the subsurface offset extended image of the same dipping up reflector as in Fig. 1(a), where too-slow velocity was used. On the right-hand side we extract two slices from the extended image. We show the z, x subspace at the offset h = 250 m to the left of the z, h subspace at the location x = 4000 m. The reflector is clearly misplaced in the image volume, due to the incorrect migration velocity, and defocuses away from the zero-offset position. The normal to the extended image is indicated by the red arrows. It denotes the angles q and p in the two subspaces, respectively. The green curve represents the two-way traveltime isochron in the extended image space. It is tangent to the extended reflector's image. On the left-hand side of the figure we illustrate the SRP geometry, where incident and scattered rays (shown in green) interact at a distance. In their study of the subsurface offset extended Born modelling approximate inverse operator, Hou & Symes (2015) provide an asymptotic analysis of the normal operator. They derive stationary phase conditions with respect to the SRP kinematics that coincide with the phase space relations in eq. (17). Moreover, they invoke these conditions to express the horizontal slowness of the one-way incident and scattered rays in terms of the phase space properties of the two-way pair   \begin{eqnarray} {\nabla _x}{t_r} = - \frac{1}{2}(q + p){\nabla _z}T \nonumber\\ {\nabla _x}{t_s} = - \frac{1}{2}(q - p){\nabla _z}T{\rm{ }}{\rm{.}} \end{eqnarray} (18) The obvious implication of the equation is that q is a bisector of the horizontal slowness. They apply the eikonal equation several times and rearrange the terms to obtain an expression for ∇zT, under the assumption of no turning waves (i.e. ∇zT > 0) . The net result is   \begin{equation}{\left| {{\nabla _z}T} \right|^2} = \frac{{b + \sqrt {{b^2} - ac} }}{a}\end{equation} (19)in which   \begin{eqnarray} a = ({q^2} + {\rm{1)}}({p^2} + 1){\rm{, }} b &=& (S_ + ^2 - S_ - ^2)qp + (S_ + ^2 + S_ - ^2){\rm{, }} \nonumber\\ c &=& {(S_ + ^2 - S_ - ^2)^2}. \end{eqnarray} (20) The slowness S± = V− 1(z, x ± h) is the reciprocal of the velocity V at the evaluation points of the SRP incident and scattered rays. Hence, the subsurface offset extension adds another degree of freedom to the image space. It enables to extract the hitting direction of the one-way incident and scattered rays at the SRP location (i.e. eq. (18)) after the construction of the extended reflector's image. Therefore, the angle-dependent image may be computed from the subsurface offset extended image itself, as suggested by eqs (5) and (6), without having to explicitly measure the direction of the one-way wavefields at each migration step. From subsurface offset to local angle domain The essential result evolving from this section of our study is that angle-domain quantities may be derived with respect to the normal direction of the SRP extended reflector. It is depicted in phase space by the ratios q and p in eq. (17). Conventionally, the reflector's image is represented in the angle-domain via a local angle domain (LAD) description (Ravve & Koren 2011). In 2-D, it consists of two angles, namely the dip-angle ν, and the scattering-angle γ. In the following, we will investigate under what conditions q and p give rise to this classical LAD description and truly amounts to the local dip and scattering angles (i.e. ν and γ). The dip-angle ν is defined as the angle subtended from the vertical axis to the normal of the reflecting interface in the z, x space. This coincides with our definition of the ratio q, therefore we state that: q = tan ν. When migration velocity is perfectly known, this definition is obviously true. The extended image is focused correctly at the zero subsurface offset, and the normal to the reflector is indeed in the direction of the structural dip. On the other hand, when velocity errors are present, the image is misplaced in the extended space and defocuses away from the zero offset trace. However, q and ν still describe the same direction normal to the reflector's fake image, although it may give a misleading indication about the structural dip of the real reflector and may become a subsurface offset dependent. The scattering-angle γ is one-half of the angle subtended by incident and scattered rays. It is commonly argued that the ratio p is by definition a measure of the same angle quantity: p = tan γ (Sava & Fomel 2003). Therefore, the angle-domain image we compute implicitly by eq. (6) is indeed a representation of the scattering-angle, and an indication of its normal in the z, h subspace. We assert that this conclusion holds only under the additional restriction that   \begin{equation}{S_ + } = {S_ - }.\end{equation} (21) In other words we claim that the ratio p defines the scattering-angle γ only when the slowness (or velocity) at the evaluation points of the SRP incident and scattered rays are identical. We demonstrate this by substituting the proposed restriction into eq. (19), and use S = S±:   \begin{equation}{\left| {{\nabla _z}T} \right|^2} = 4{S^2}\frac{1}{{({q^2} + 1)({p^2} + 1)}}.\end{equation} (22) We now sum eq. (18) to obtain the magnitude square of the SRP two-way traveltime gradient in the physical z, x subspace   \begin{equation}{\left| {{\nabla _{z,x}}T} \right|^2} = {\left| {{\nabla _z}T} \right|^2} + {\left| {{\nabla _x}T} \right|^2} = ({q^2} + 1){\left| {{\nabla _z}T} \right|^2},\end{equation} (23)and substitute eq. (22)   \begin{equation}{\left| {{\nabla _{z,x}}T} \right|^2} = 4{S^2}\frac{1}{{({p^2} + 1)}}.\end{equation} (24) On the other hand, we use the law of cosines to geometrically evaluate γ as the half-opening angle between incident and scattered rays   \begin{equation}{\left| {{\nabla _{z,x}}T} \right|^2} = {\left| {\nabla {t_r} + \nabla {t_s}} \right|^2} = S_ + ^2 + S_ - ^2 + 2S_ + ^2S_ - ^2\cos 2\gamma ,\end{equation} (25)and introduce once again our restriction by using S = S±,   \begin{equation}{\left| {{\nabla _{z,x}}T} \right|^2} = 4{S^2}{\cos ^2}\gamma .\end{equation} (26) Comparing eqs (24) and (26) clearly indicates that p = tan γ only when S+ = S− is imposed. This restriction holds for laterally invariant velocity models (i.e. assuming horizontally layered earth). But, in the general case of arbitrary velocity function, it holds only at the zero subsurface offset. Hence, it can be taken for granted only when migration velocity is correct and the image focuses at the zero offset. As the velocity estimation deviates from true, the image defocuses, and the angle denoted by the ratio p may no longer be a valid description of the scattering-angle γ. We further emphasize this issue with the illustrations in Fig. 2. An interaction between incident and scattered rays takes place in a hypothetic model of horizontal reflector and two vertical velocity blocks. Velocity in the right block is larger than in the left block. Fig. 2(a) presents a SRP interaction that violates eq. (21) restriction: Velocity changes at the split reflection coordinates, and the slowness vectors do not have the same length. Therefore, the dip direction q (indicated in red) is not a bisector of the opening angle. The two angles subtended from q to the incident and scattered rays are not the same, what leads to an ambiguous definition of the scattering-angle γ. However, the dip q does act as a bisector of the ray's horizontal slowness, as implied by eq. (18). Fig. 2(b) shows another SRP interaction where velocity is the same at the split reflection coordinates, and the slowness vectors have equal length. In such case, the dip q is a bisector of the opening angle as well as a bisector of the horizontal slowness, and the scattering-angle γ is by definition the direction p. In Fig. 2(c) a CRP interaction is illustrated. The dip q in this conventional common-point geometry is always a bisector of the opening angle between incident and scattered rays. Figure 2. View largeDownload slide An interaction between incident and scattered rays in a model with single horizontal reflector and two vertical velocity blocks. When velocity is the same at the split reflection point coordinates, the reflector normal is a bisector of the opening angle (panels b and c). Otherwise, it acts as a bisector of the slowness vectors (panel a). Figure 2. View largeDownload slide An interaction between incident and scattered rays in a model with single horizontal reflector and two vertical velocity blocks. When velocity is the same at the split reflection point coordinates, the reflector normal is a bisector of the opening angle (panels b and c). Otherwise, it acts as a bisector of the slowness vectors (panel a). To conclude this discussion we stress that the angle-dependent image, derived implicitly from the subsurface offset extension in eq. (6), still carries valuable traveltime information for MVA. Despite of being non-physical, the SRP kinematic description gives rise to a unique geometry of split incident and scattered rays (see eq. 18), and should be adopted accordingly by conventional traveltime inversion methods for velocity optimization. We devote the rest of this paper to emphasize the uniqueness of the SRP description when compared against the conventional CRP description. ANGLE-DOMAIN MOVEOUT ANALYSIS We have introduced two angle-domain imaging techniques in this work to compute scattering-angle CIGs. An explicit method for angle decomposition is described in eqs (2) and (3), whereas an implicit transformation is suggested by eqs (5) and (6) via a subsurface offset image extension. Each one of the methods exploits different reflection geometry we regard as CRP or SRP respectively, and has a unique kinematic description. Thus, although describing the same angle quantity (i.e. the scattering-angle), the CRP and the SRP related angle-domain CIGs are not expected to show the same moveout shape when incorrect migration velocity is used (Bartana et al. 2006; Montel & Lambare 2013). We first demonstrate the moveout dissimilarity with the Marmousi synthetic model (Bourgeois et al. 1990). Figs 3(a) and (b) show the CRP (left side) and SRP (right side) angle-domain images, computed with true and erroneous migration velocity, respectively. Vm(z) = 1500 + 1.33 · z m s−1. The scattering-angle CIGs are extracted in the figures around the x = 2500 m mark. When migration velocity is true, the CRP and SRP angle-domain images are essentially the same. Asymptotically, the same structure of the subsurface stratigraphy can be extracted from the two, and they both produce flat CIGs. On the other hand, the CRP and SRP images are unlike when false migration velocity is used. Although the two images seem to similarly misplace the reflectors in the image section (i.e. at zero angle of incidence), they are strongly distinguished by their angle-dependent moveout. The moveout bends in the same direction but not with the same curvature or extent. Some evident moveout differences are shown at shallow depth and in the deep section as well. We highlight in yellow one particular reflection event for comparison. Its moveout is extracted by automatic picking, and displayed in Fig. 4. The figure shows a map view of the moveout picks. The CRP (left side) and SRP (right side) moveout maps are clearly distinguishable, especially in the mid and far angle range. The following analysis studies the essence of this moveout dissimilarity in a more quantitative manner. The net result implies that in spite of the peculiar discrepancy, both types of CIGs provide essential information about the velocity error as long as it is extracted through the appropriate CRP or SRP geometry. Figure 3. View largeDownload slide Marmousi model example. The scattering-angle extended image, computed by CRP (left) and SRP (right) angle-domain imaging operators, by using (a) true or (b) erroneous migration velocity. The zero-angle image section is shown together with several CIGs, extracted at the locations: x = 1300, 1700, 2100, 2500, 2900, 3300 and 3700 m. Figure 3. View largeDownload slide Marmousi model example. The scattering-angle extended image, computed by CRP (left) and SRP (right) angle-domain imaging operators, by using (a) true or (b) erroneous migration velocity. The zero-angle image section is shown together with several CIGs, extracted at the locations: x = 1300, 1700, 2100, 2500, 2900, 3300 and 3700 m. Figure 4. View largeDownload slide Marmousi model example. A map view of the moveout picks demonstrates the dissimilarity between CRP (left) and SRP (right) angle-domain imaging methods. Figure 4. View largeDownload slide Marmousi model example. A map view of the moveout picks demonstrates the dissimilarity between CRP (left) and SRP (right) angle-domain imaging methods. Consider a single reflector in a homogenous media of velocity V. Seismic data, acquired above this reflector will asymptotically show a reflection's traveltime curve with respect to the midpoint location xm and the half acquisition offset H. The migration of the acquired data maps this traveltime information to the image space along isochrons. We have briefly introduced these isochrons in eqs (7) and (13). The migration map may be incorrect when the migration velocity Vm is other than the true velocity V. In homogenous media, and with respect to the CRP geometry, the isochron curves are a family of ellipses in the z, x image space (Dafni & Symes 2016a) and take the form   \begin{equation}T(z,x;H) = \frac{{{{(x - {x_m})}^2}}}{{{\varepsilon ^2}(z_0^2 + {H^2})}} + \frac{{{z^2}}}{{{\varepsilon ^2}(z_0^2 + {H^2}) - {H^2}}} - 1 = 0,\end{equation} (27)where ε = Vm/V is the relative velocity error, and z0 is the zero-dip imaging depth. Each member in eq. (27) is subjected to a given half-offset value H, and the resulting image is obtained by the envelope of all family members. The CRP related angle decomposition technique stores the scattering-angle information of each one of the isochron curves right before they interact to form the image envelope. Under the assumption of constant velocity, the angles have a simple map relation with the acquisition offset. For simplicity, we assume that a horizontal reflector is embedded. This makes the problem laterally invariant, hence the image point and the midpoint coincide: x − xm = 0, and the offset to angle map reads: H = zp. Substituting these two realizations into eq. (27), and rearranging the terms, yields an expression for the CRP related angle-domain moveout curve   \begin{equation}z(p) = \frac{{\varepsilon {z_0}}}{{\sqrt {1 - {p^2}({\varepsilon ^2} - 1)} }}.\end{equation} (28) Dafni & Symes (2016a) extended the family of isochrons in eq. (27) with the subsurface offset extra degree of freedom to represent the SRP kinematic description. Their proposal for the extended isochrons takes the elliptic form   \begin{eqnarray} T(z,x,h;H) &=& \frac{{{{(x - {x_m})}^2}}}{{{\varepsilon ^2}(z_0^2 + {H^2})}}\,{ +}\, \frac{{{z^2}}}{{{\varepsilon ^2}(z_0^2 + {H^2}) - {{(h - H)}^2}}}\,{ -}\, 1 = 0.\!\!\!\!\!\!\nonumber\\ \end{eqnarray} (29) Note that each family member becomes a function of the additional subsurface offset coordinate. Therefore, enveloping eq. (29) members constructs the subsurface offset extended image. We derive an expression for such envelope in Appendix  A, by considering a horizontal reflector is under study. The envelope reads   \begin{equation}z(h) = \varepsilon \sqrt {z_0^2 - \frac{{{h^2}}}{{({\varepsilon ^2} - 1)}}} ,\end{equation} (30)and describes the defocusing of the reflector's image in subsurface offset CIGs. Dafni & Symes (2016b) transformed eq. (29) into the angle-domain with a parametric Radon transformation, to derive the scattering-angle variant of the extended traveltime isochron. They obtained the family of hyperbolas   \begin{eqnarray} T(z,x,p;H) &=& z - \sqrt {\left({\varepsilon ^2}\left(z_0^2 + {H^2} \right) - {{(x - {x_m})}^2}\right)(1 + {p^2}{\eta ^2})}\nonumber\\ && +\, pH = 0, \end{eqnarray} (31)where the stretch factor η is given by   \begin{equation} \eta = \sqrt {\frac{{{\varepsilon ^2}\left(z_0^2 + {H^2} \right)}}{{{\varepsilon ^2}\left(z_0^2 + {H^2}\right) - {{(x - {x_m})}^2}}}} . \end{equation} (32) The envelope of eq. (31) members represents the angle-domain image, derived in relation with the subsurface offset extension and the SRP kinematic description. We derive an expression for such envelope in Appendix  A, under the assumption of horizontal reflector that gives raise to the SRP related angle-domain moveout curve   \begin{equation}z(p) = {z_0}\sqrt {{\varepsilon ^2} + {p^2}({\varepsilon ^2} - 1)} .\end{equation} (33) Eqs. (33) and (30) maintain the following parametric Radon transform relation:   \begin{equation}z(p) = z(h) - ph = z(h) - \frac{{\partial z(h)}}{{\partial h}}h.\end{equation} (34) Since eq. (33) moveout curve represents homogenous media, the velocity at the evaluation points of the incident and scattered rays are the same. Therefore, we reassure by eq. (21) restriction that the SRP moveout description is indeed of the scattering-angle, and that p = tan γ holds. Comparing eqs (28) and (33) immediately shows the angle-domain moveout dissimilarity. Despite describing depth with respect to the same scattering-angle of reflection by definition, two different expressions evolve. In Figs 5(a) and (b), we draw the CRP and SRP ray trajectories of a single shot-receiver ray-pair, respectively. The dashed green rays represent propagation with true velocity. The SRP coincides with the CRP in this special case, and reflections occur at the true reflector's depth z0. Physical reflectors that are being imaged with the true velocity cannot be explained by any subsurface offset other than zero. The red ray trajectories correspond to a too-fast (left side) and too-slow (right side) migration velocity. When velocity is incorrect, the ray geometry is wrong. However, it is wrong differently with regards to the CRP and SRP description. The red CRPs are misplaced in depth only, whereas the red SRPs are misplaced in depth and in the lateral position. Moreover, the geometry differences amount to a different evaluation of the scattering-angle as the half-opening angle of reflection. With respect to the CRP geometry it is calculated by the ratio H/z, whereas with respect to the SRP geometry it is the ratio (H − h)/z. We reemphasize again that both evaluations are correct by definition, as long as the correct geometry is kept adequate. The SRP misleading reflection depth is transformed into its equivalent angle-domain depth by the parametric Radon transform in eq. (34). It is illustrated in Fig. 5(b) with the purple dashed segments. Figure 5. View largeDownload slide The geometry of (a) CRP and (b) SRP reflections. The two geometries coincide when migration velocity is true (green rays). They differ considerably when migration velocity is too-fast (left) or too-slow (right) (red rays). The purple trajectory demonstrates the angle transformation of the SRP reflection. Figure 5. View largeDownload slide The geometry of (a) CRP and (b) SRP reflections. The two geometries coincide when migration velocity is true (green rays). They differ considerably when migration velocity is too-fast (left) or too-slow (right) (red rays). The purple trajectory demonstrates the angle transformation of the SRP reflection. In Fig. 6 we plot the envelope curves that have been derived in this section for the flat reflector example. Eq. (30) subsurface offset domain defocusing curve is shown in Fig. 6(a), whereas eqs (28) and (33) angle-domain moveout curves are in Fig. 6(b). We evaluate the equations with 10 per cent too-fast, true, and 10 per cent too-slow migration velocity (left to right). The dashed lines demonstrate the transformation of the point ‘A’ on the subsurface offset domain envelope to its equivalent point ‘B’ on the SRP angle-domain envelope (i.e. eq. (34) transformation), and vice versa. Note that when true migration velocity is used, there is no defocusing in subsurface offset, and the scattering-angle domain becomes redundant. Yet, when migration velocity is wrong defocusing develops, and the angle-domain moveout curves coincide only at γ = 0, that corresponds to h = 0, and to the migration of the H = 0 data trace. Figure 6. View largeDownload slide (a) Subsurface offset and (b) angle-domain envelope curves of horizontal reflector in homogenous media. 10 per cent too-fast, true, and 10 per cent too-slow velocities are used from left to right, respectively. The dashed line maps between points on the subsurface offset and angle-domain envelopes. Figure 6. View largeDownload slide (a) Subsurface offset and (b) angle-domain envelope curves of horizontal reflector in homogenous media. 10 per cent too-fast, true, and 10 per cent too-slow velocities are used from left to right, respectively. The dashed line maps between points on the subsurface offset and angle-domain envelopes. ANGLE-DOMAIN MOVEOUT CURVES IN LAYERED MEDIA We provide a semi-analytic method to simulate the CRP and SRP angle-domain moveout curves in a piecewise constant velocity model, consisting of L horizontal layers. We assume the true velocity Vi ∈ V is known, and evaluate the moveout with respect to a given migration velocity error εi ∈ ε and a fixed set of acquisition offsets Hj ∈ H. In our notation, the index i specifies a layer within the model and the index j indicates one of the offset traces. For simplicity, we restrict the acquisition surface to z = 0. The simulation scans over the take-off angle of incident and scattered rays in the search for the direction that would give rise to a reflection traveltime in the seismic data. It can be regarded as a traveltime map between the data and model spaces. A two-point ray tracing in the z, x subspace connects a data point td(0, H) at the recording surface with a subsurface model point tm(z, h) in the extended image volume, and vice versa. Since the model is laterally invariant, the choice for the origin x = 0 is arbitrary, and set to the location where the midpoint and image point coincide. The forward map models the data traveltimes with a given velocity V, and by restricting the model points to h = 0 at the true depth of the reflecting layers   \begin{equation}{t_m}(z = {z_0},h = 0)\stackrel{V}{\longrightarrow}{t_d}(0,H).\end{equation} (35) Fig. 7(a) illustrates the forward map. A fan of rays is shot from a model point at the bottom interface of the ith layer. The rays propagate upward until they reach the acquisition surface. The traveltime recorded at the jth acquisition offset data trace is calculated along the ray hitting the surface exactly at that offset. Note that the incident and scattered rays are symmetric with respect to the midpoint location at the origin x = 0. They propagate with the same traveltimes and opposite angles. Hence, one can search for the incident ray solely and duplicate the scattered ray in accordance. Figure 7. View largeDownload slide Traveltime map between data and model spaces by two-point ray-tracing. (a) The forward modelling map, and the (b) CRP and (c) SRP backward migration maps. Figure 7. View largeDownload slide Traveltime map between data and model spaces by two-point ray-tracing. (a) The forward modelling map, and the (b) CRP and (c) SRP backward migration maps. The backward map migrates the data traveltimes to the extended image volume with the migration velocity Vm, and may be wrong when Vm ≠ V  \begin{equation}{t_m}(z,h)\stackrel{V_m}{\longleftarrow}{t_d}(0,H).\end{equation} (36) Eq. (36) has been derived as an SRP migration map. Imposing h = 0 in the equation yields the equivalent CRP migration map. Figs 7(b) and (c) illustrate the CRP and SRP migration maps, respectively. A fan of rays is shot at the acquisition surface with respect to the jth offset. The rays are traced downward towards the ith reflecting interface, until the traveltime along the rays matches its equivalent in the data: tm = td. The migration map is complete by those rays that satisfy a certain geometric condition with respect to the subsurface offset split. This condition is simply h = 0 for the CRP migration map (see Fig. 7b). As for the SRP migration map, the split h at the reflection point has to do with the envelope construction of the extended image. It is determined by the contact relation of the envelope theorem, as introduced in eq. (A6) for a single layer model. We provide Appendix  B to formulate a generalized contact relation for multilayer models. The net result is presented in eq. (B6). Once the migration map rays are found, the CRP and SRP angle-domain moveout curves are constructed by superposing the end points of these rays. The corresponding scattering-angle at the end point of the rays is evaluated by means of Snell's law with respect to the take-off angle at the starting point. We summarize the scheme of our method in Table 1. Note that the method predicts the asymptotic shape of the angle-domain moveout curves without explicitly executing a time consuming and computationally expensive migration job. However, to validate the results in the following examples we do compute angle-domain CIGs via migration operators extended according to the CRP and SRP description. We employ a ray-based variation of eqs (2) and (3) as our CRP related migration operator (Koren & Ravve 2011; Ravve & Koren 2011), and a subsurface offset extended Born-type RTM as our SRP related operator (see eqs 5 and 6). Table 1. Simulation scheme. Asymptotic calculation of CRP and SRP angle-domain moveout curves.     View Large The subsequent subsections exemplify the implementation of the proposed moveout simulation technique. Although assuming a horizontally layered earth, the significance of the proceeding examples is indicated by the asymptotic and semi-analytic analysis of reflections in CIGs. The examples gradually expose the narrative of this paper and demonstrate the fundamental difference between the kinematics of the CRP and SRP reflections. The quantitative approach enables to easily distinguish between the two, and understand how differences in the reflection geometry amount to dissimilar ray paths and unlike moveout in angle-domain CIGs. We start with a simple single-layer over a half-space example to introduce the concept, and then proceed with multilayer examples where the moveout behaviour is non-trivial and more complicated, but still illuminative in terms of the reflection configuration contrast. Single-layer example The proposed simulation technique is demonstrated first with a single-layer model, consisting of a horizontal reflecting interface at z0 = 2000 m, and constant velocity V = 2000 m s−1. The forward map is followed to calculate reflection traveltimes with respect to the acquisition half-offsets H = {–2500, …, 2500} m. Fig. 8(a) illustrates in green the ray trajectories of this map, connecting shot-receiver pairs with a reflection point at the reflector. For the backward migration map an erroneous migration velocity is employed. Figs 8(b) and (c) shows the migration rays of the CRP (in blue) and SRP (in magenta) maps. 10 per cent too-fast and 10 per cent too-slow velocity is used in these figures, respectively. The CRP rays are restricted to meet at the midpoint location where incident and scattered rays share a common point of reflection (i.e. h = 0). Each of the corresponding reflection points is wrongly placed in depth, implying that the migration map is wrong and offset (or angle) dependent. The SRP rays, on the other hand, are allowed to defocus horizontally from the midpoint location, and split the incident from the scattered ray at the reflection point. The wrong velocity SRP migration map misplaces the reflection points in depth, as well as in the lateral position. Figure 8. View largeDownload slide Simulation of ray trajectories—single-layer example. (a) The forward modelling map, and the CRP (blue) and SRP (magenta) backward migration maps by using (b) 10 per cent too-fast and (c) 10 per cent too-slow velocities. Figure 8. View largeDownload slide Simulation of ray trajectories—single-layer example. (a) The forward modelling map, and the CRP (blue) and SRP (magenta) backward migration maps by using (b) 10 per cent too-fast and (c) 10 per cent too-slow velocities. Asymptotically, the collection of the rays end-points at the reflection point constructs the subsurface offset defocusing pattern and angle-domain moveout curve, one should expect to observe in the corresponding CIGs. On the right-hand side of Figs 9 and 10 we extract the end-points of Fig. 8 ray trajectories and plot them with respect to the subsurface offset and the scattering-angle, respectively. Note that these end points match eqs (28), (30) and (33) analytic curves, in the same way as in Fig. 6 illustrations, even though the extent of the curves is bounded by the acquisition geometry. On the left-hand side of Figs 9 and 10 we plot this example's subsurface offset and scattering-angle CIGs. The images were computed with angle-domain CRP and SRP migration operators as mentioned previously. Imaging with a subsurface offset extension often leads to the appearance of some artefacts in the CIGs. It is evident in the most left gather in Figs 9(a) and (b) as indicated by the red arrows. Dafni & Symes (2016b) regard these artefacts as kinematic artefacts that emerge due to the truncation of the seismic data at a maximum acquisition offset. In their paper, they also propose a dip-domain specularity filter to suppress the artefacts. We apply their filter and show the results in the middle panel of Figs 9(a) and (b). The filtered CIGs are clearly free of the kinematic artefacts. Thus, the SRP angle-domain CIGs in Fig. 10 are the angle-domain transformation of those filtered subsurface offset CIGs. Figure 9. View largeDownload slide Simulation of subsurface offset defocusing—single-layer example. The simulated curves (right) and the computed subsurface offset CIGs (left), obtained by using (a) 10 per cent too-fast and (b) 10 per cent too-slow migration velocities. A specularity filter is applied in the middle panels to suppress truncation artefacts. Figure 9. View largeDownload slide Simulation of subsurface offset defocusing—single-layer example. The simulated curves (right) and the computed subsurface offset CIGs (left), obtained by using (a) 10 per cent too-fast and (b) 10 per cent too-slow migration velocities. A specularity filter is applied in the middle panels to suppress truncation artefacts. Figure 10. View largeDownload slide Simulation of angle-domain moveout—single-layer example. The simulated curves (right) and the computed subsurface offset CIGs (left), obtained by using (a) 10 per cent too-fast and (b) 10 per cent too-slow migration velocity. The CRP gathers are presented to the left of the SRP gathers. Figure 10. View largeDownload slide Simulation of angle-domain moveout—single-layer example. The simulated curves (right) and the computed subsurface offset CIGs (left), obtained by using (a) 10 per cent too-fast and (b) 10 per cent too-slow migration velocity. The CRP gathers are presented to the left of the SRP gathers. We overlay the images in Figs 9 and 10 with the simulated curves to confirm that the right kinematic behaviour is predicted. The curves certainly follow the wavelet's zero-crossing across the CIGs traces, and stretch to the right extent. It is also evident that the CRP and SRP angle-domain curves are different in shape and size. The same set of acquisition offsets in the data illuminates a different range of scattering-angles in the images of Fig. 10. Once again, this is a clear implication of a split, rather than common, reflection point geometry. Multilayer examples The next two examples demonstrate our method with blocky multilayer models (i.e. constant velocity in each layer). Here, we expand the acquisition coverage and calculate the forward map reflection traveltimes with respect to H = {–3750, …, 3750} m. The first out of the two examples consists of three layers above a half space. The layers true velocity and depth are presented in the left columns of Table 2. The migration velocity and the relative velocity error are shown to the right. Figs 11(a) and (b) show the ray trajectories of reflections from the second and third layers, respectively. The rays of the forward modelling map and the backward migration maps are coloured the same as in the preceding example. Here we would like to stress that erroneous migration velocity not only misplaces the reflection point in space, it also wrongly propagates the rays and amounts to a misleading ray path. The modelling ray trajectories on the left part of Fig. 11 clearly do not match those of the migration maps. Moreover, since the CRP and SRP migration maps do not share the same reflection geometry, their rays do not share the same path as well. In Fig. 12 we extract the CRP and SRP reflection points from the trajectories in Fig. 11, and plot them as a function of the subsurface offset and scattering-angle. The set of all points results as the defocusing/moveout curves of this example's simulation. To the left, we compute the corresponding CIGs, and overlay with the simulated curves. Fig. 12(a) presents the original (left panel) and the specularity filtered (middle panel) subsurface offset CIGs. Some of the detracting kinematic artefacts are successfully suppressed, and the defocused events are kept intact. Fig. 12(b) shows the CRP (left panel) and SRP (middle panel) scattering-angle CIGs. The simulated curves obviously follow the three reflection events in the CIGs to their full extent. Figure 11. View largeDownload slide Simulation of ray trajectories—three-layer example. The forward modelling map (green), and the CRP (blue) and SRP (magenta) backward migration maps. Reflections from the bottom of (a) the second and (b) the third layers are shown. Figure 11. View largeDownload slide Simulation of ray trajectories—three-layer example. The forward modelling map (green), and the CRP (blue) and SRP (magenta) backward migration maps. Reflections from the bottom of (a) the second and (b) the third layers are shown. Figure 12. View largeDownload slide Simulation of (a) subsurface offset defocusing and (b) angle-domain moveout—three-layer example. The simulated curves (right) and the computed CIGs (left). A specularity filter is applied in the upper middle panel to suppress truncation artefacts. Figure 12. View largeDownload slide Simulation of (a) subsurface offset defocusing and (b) angle-domain moveout—three-layer example. The simulated curves (right) and the computed CIGs (left). A specularity filter is applied in the upper middle panel to suppress truncation artefacts. Table 2. Three layer example. Depth–velocity model and migration velocity error. Layer No.  z0 (m)  V (m s−1)  Vm (m s−1)  ε  1  1000  2500  2250  0.9  2  2000  2300  2645  1.15  3  3000  3200  2720  0.85  Layer No.  z0 (m)  V (m s−1)  Vm (m s−1)  ε  1  1000  2500  2250  0.9  2  2000  2300  2645  1.15  3  3000  3200  2720  0.85  View Large The second multilayer example demonstrates a more peculiar case. It consists of four layers above a half-space. We design the 1-D velocity function here to oscillate between the model layers, and set the migration velocity as constant. See Table 3 for the velocity values and errors. Note that this example's errors are quite large, so we anticipate the migration maps to be completely off with respect to the modelling map. Figs 13(a) and (b) present the ray trajectories of reflections from the second and fourth layers, respectively, with the same colour scheme as in the previous examples. The CRP and SRP migration maps obviously project wrong trajectories when compared against the modelling map, and misplace the reflection point in space. However, the reflection from the second layer (Fig. 13a) is somewhat strange. The CRP and SRP migrations clearly do not map the reflection point to its true depth at z = 2000 m, although it seems like these maps are still angle (or offset) independent. In other words, the second layer reflection event in all data traces is mapped by migration to about the same, but wrong, depth position. There is no sign of significant defocusing or moveout. This gives rise to the familiar ‘velocity–depth ambiguity’ of depth imaging: Reflections in the data can match several depth–velocity models, other than the true model. The traveltime error along the rays is accumulated as the rays travel through the layers. This example's depth–velocity model is designed in such a way that the traveltime error is averaged and effectively vanishes at the second layer interface. Moreover, since the third layer's migration velocity is correct, we expect the same ambiguity to occur at the third layer interface as well. Figure 13. View largeDownload slide Simulation of ray trajectories—four-layer example. The forward modelling map (green), and the CRP (blue) and SRP (magenta) backward migration maps. Reflections from the bottom of (a) the second and (b) the fourth layers are shown. Figure 13. View largeDownload slide Simulation of ray trajectories—four-layer example. The forward modelling map (green), and the CRP (blue) and SRP (magenta) backward migration maps. Reflections from the bottom of (a) the second and (b) the fourth layers are shown. Table 3. Four layer example. Depth–velocity model and migration velocity error. Layer No.  z0 (m)  V (m s−1)  Vm (m s−1)  ε  1  1000  1500  2000  1.33  2  2000  2500  2000  0.8  3  3000  2000  2000  1.0  4  4000  3000  2000  0.66  Layer No.  z0 (m)  V (m s−1)  Vm (m s−1)  ε  1  1000  1500  2000  1.33  2  2000  2500  2000  0.8  3  3000  2000  2000  1.0  4  4000  3000  2000  0.66  View Large We extract the simulated defocusing and moveout curves as a function of the subsurface offset and scattering-angle, and plot them on the right-hand side of Fig. 14. The ambiguous focusing of the second and third reflections is clearly predicted in the subsurface offset domain. It translates into an ambiguous flattening in the angle-domain. This example's CIGs are computed to the left in Fig. 14. Fig. 14(a) presents the original (left panel) and the specularity filtered (middle panel) subsurface offset CIGs, whereas Fig. 14(b) presents the CRP (left panel) and SRP (middle panel) scattering-angle CIGs. The same ambiguity is shown in the computed gathers, although resolution is limited to the wavelength size. We should also mention that our simulation successfully match the other two events in the CIGs (first and fourth reflections) to their full extent. While looking more closely at the second reflection event, we note that the ambiguous focusing (or flatness) is not perfectly complete. It is demonstrated in Fig. 15 where the second event in extracted. The asymptotic simulation predicts a small-scale bowtie shaped defocusing pattern in the subsurface offset domain (right side of Fig. 15a). Its scale is way below the wavelength size, and therefore it cannot be noted in the subsurface offset CIG (left side of Fig. 15a). One can only qualitatively indicate a rather focused event. The bowtie shape implies that the subsurface offset representation is multivalued. Meaning that a given subsurface offset may be associated with more than one depth. Angle-domain imaging on the other hand, would be a preferable choice in such case. The scattering-angle representation uniquely describes reflections in depth, as shown in Fig. 15(b). The angle quantity unfolds the traveltime information, and shows some mild uplift (about 20 m) at the near angle range. There are no multivalued angles along the moveout curve of Fig. 15(b). Image-domain traveltime inversion methods would surely prefer the unpacked angle representation over the collapsed and multivalued subsurface offset representation. Either way, to avoid any kind of velocity-depth ambiguity, one should use a layer stripping approach and search for the optimal model in a top-down manner. Figure 14. View largeDownload slide Simulation of (a) subsurface offset defocusing and (b) angle-domain moveout—four-layer example. The simulated curves (right) and the computed CIGs (left). A specularity filter is applied in the upper middle panel to suppress truncation artefacts. Figure 14. View largeDownload slide Simulation of (a) subsurface offset defocusing and (b) angle-domain moveout—four-layer example. The simulated curves (right) and the computed CIGs (left). A specularity filter is applied in the upper middle panel to suppress truncation artefacts. Figure 15. View largeDownload slide Simulation of (a) subsurface offset defocusing and (b) angle-domain moveout—four-layer example. A closer look at the second event in Fig. 12. Figure 15. View largeDownload slide Simulation of (a) subsurface offset defocusing and (b) angle-domain moveout—four-layer example. A closer look at the second event in Fig. 12. DISCUSSION The obvious application of our study is angle-domain kinematic analysis for velocity model optimization. In particularly, we infer to adapt the SRP kinematic description for that matter. Modifying conventional traveltime inversion methods with respect to a split, rather than common, reflection point geometry would link between the defocusing/moveout patterns in the SRP image gathers and the migration velocity error. In spite of being described in 2-D, the SRP description extends similarly to 3-D. Practically speaking, it amounts to a computation of a multidirectional subsurface offset extension in both horizontal coordinates, and an azimuthal-dependent angle transformation (Dafni & Symes 2016a). Subsurface offset extended imaging is mostly associated with wave-based migration methods, such as RTM. Hence, deriving the appropriate SRP inversion methods would allow to optimize the velocity model by processing image gathers, produced by RTM. This may become beneficial in areas with complex geologic structure, where wave-based imaging is superior to ray-based imaging. Another benefit of the SRP description has to do with its computational cost. Computing the SRP angle-domain CIGs via an implementation of subsurface offset extended RTM (eqs 5 and 6) is computationally intensive, but not as much as computing the CRP angle-domain CIGs with the explicit angular extension of RTM (eqs 2 and 3). The latter computation decomposes the wavefields into the angle-domain at each migration time step, prior to the imaging condition, and therefore tends to be more demanding in computational time and storage (Vyas et al. 2011). One peculiar characteristic of the SRP kinematic description is the legitimacy of the angle-domain transform as a measure of the half-opening angle at the reflection point (i.e. the scattering-angle). We extensively discuss the validity of p = tan γ in eq. (6) transform and assert that it can only be true under eq. (21) restriction. Meaning that the angle transform's slope p is by definition a representation of the scattering-angle γ only when the velocity at the split reflection point is the same for the incident and scattered waves. In any other case, the reflector normal is a bisector of the wave-pair horizontal slowness (see eq. 18) rather than a bisector of the opening angle. This result must be addressed appropriately while adapting the SRP kinematic description for traveltime inversion. Nevertheless, there is a way to loose eq. (21) restriction and make p a definition of γ more generally. It follows the ambitious attempt of Biondi & Symes (2004) to compute the subsurface offset extension in the direction of the imaged dip, instead of taking the horizontal offset distance. Splitting the incident from the scattered waves along the local reflecting interface is more likely to result with the same velocity value at the split reflection point. However, the validity of the SRP angle-domain representation is not the only reason way the CRP and the SRP differ in their moveouts. In fact, we ensure that in all of this work numerical examples p is truly a definition of the scattering-angle γ. The other reason for the moveout dissimilarity is the geometry of reflection. The CRP and SRP geometries give rise to different travel paths of incident and scattered wave pairs, what results in two unlike kinematic descriptions. The moveout dissimilarity generally increases with scattering-angle, and is likely to be more prominent as the migration velocity error grows, or when large lateral velocity variations are present. Moreover, as a rule of thumb we find that the more extended the subsurface offset defocusing, the bigger the angle-domain moveout dissimilarity. This work's kinematic analysis is based on aspects from envelope theorem. We find the theorem a practical approach to unfold the seismic image into its fundamental building blocks. Migration is asymptotically regarded as a superposition of two-way traveltime isochrons, whereas the image is considered as the resulting envelope. It is demonstrated with horizontally layered models, although treating the image as an envelope of isochron curves still holds in general with more realistic subsurface models. Therefore, our analysis may be upgraded by numerical means to address more structurally complex models and some level of anisotropy. Subsequently, the kinematic analysis evolved into a simulation of angle-domain moveout curves. The obvious impact of the proposed simulation is the ability to extensively study the CRP and SRP moveout behaviour without a costly execution of a migration algorithm. It is a fast and computationally cheap technique to predict the shape and extent of reflections in CIGs. This a priori prediction may also be useful to better estimate the optimal image extension range (i.e. maximum subsurface offset or angle) in a consequent migration job, and reduce some of the associated cost. However, the simulation obviously fails to predict artefacts or coherent noise in the CIGs, like the truncation artefacts associated with the subsurface offset extension (Dafni & Symes 2016b). It also lacks the essential ability to handle caustics and multiarrivals. To validate the accuracy of the moveout simulation we apply migration and compute angle-domain CIGs. Then, we match the asymptotic simulation results with the band-limited events in the CIGs. This match has a certain uncertainty in the scale of the wavelength size. Its quality often drops at the far angle range, where stretch effects are likely to occur. More than that, truncation artefacts in the imaging process are another factor threatening to weaken the match quality, especially at the far angle range. Those artefacts may also impair the convergence of traveltime inversion methods. We apply a dip-domain specularity filter (Dafni & Symes 2016b) to suppress those artefacts and enhance the true kinematic behaviour. However, adaptive image domain filters, like the specularity filter, do not always result with complete removal of the artefacts and must be applied cautiously. Their application should be done with adequate level of quality control. CONCLUSIONS Angle-domain imaging is alternatively derived by associating reflections with split, rather than common, reflection point geometry. It is based on a subsurface offset extended migration followed by a conventional Radon transform. A semi-analytic simulation demonstrates that the resulting angle-domain moveout differs from the conventional moveout curve, what implies about a unique kinematic description. The split reflection geometry is clearly a non-physical description, but surely delineates a valid traveltime narrative for the benefit of traveltime inversion and velocity model optimization. Analysis of subsurface offset extended reflections and imaged reflectors suggests that they match in their phase space coordinates as an indication of the underlying direction of scattering. It turns out that the phase space quantities relate to the local dip and scattering angles. They define those angles explicitly in the special case when the velocity does not change at the split reflection point. In the general case, this definition obviously does not hold, but still valuable. Meaning that, conventional angle-domain traveltime inversion would probably fail to converge without being modified accordingly with the subsurface offset split. Adapting the split geometry in the inversion scheme is a crucial step towards its general and robust success. The only requirement is a basic understanding of how the geometry difference translates into a compatible traveltime description. ACKNOWLEDGEMENTS We are grateful to the sponsors of The Rice University Inversion Project (TRIP) for their long-term support, in particular to Shell International Exploration and Production Inc. for its financial support. We thank Paradigm for using their full-azimuth angle domain decomposition and analysis system, and for partially supporting this work. We have also benefited greatly from the IWAVE package developed by TRIP to carry out our examples, and from the high performance computing resources provided by the Rice University Research Computing Support Group. REFERENCES Al-Yahya K., 1989. Velocity analysis by iterative profile migration, Geophysics , 54( 6), 718– 729. https://doi.org/10.1190/1.1442699 Google Scholar CrossRef Search ADS   Bartana A., Kosloff D., Ravve I., 2006. Discussion and Reply, Geophysics , 71( 1), X1– X3. https://doi.org/10.1190/1.2168013 Google Scholar CrossRef Search ADS   Biondi B., Almomin A., 2014. Simultaneous inversion of full data bandwidth by tomographic full-waveform inversion, Geophysics , 79( 3), WA129– WA140. https://doi.org/10.1190/geo2013-0340.1 Google Scholar CrossRef Search ADS   Biondi B., Symes W.W., 2004. Angle-domain common-image gathers for migration velocity analysis by wavefield-continuation imaging, Geophysics , 69( 5), 1283– 1298. https://doi.org/10.1190/1.1801945 Google Scholar CrossRef Search ADS   Bourgeois A., Bourget M., Lailly P., Poulet M., Ricarte P., Versteeg R., 1990. Marmousi, model and data, in EAEG Workshop—Practical Aspects of Seismic Data Inversion , doi:10.3997/2214-4609.201411190. Claerbout J.F., 1985. Imaging the Earth's Interior , Blackwell Scientific Publ. Dafni R., Symes W.W., 2016a. Scattering and dip angle decomposition based on subsurface offset extended wave-equation migration, Geophysics , 81( 3), S119– S138. https://doi.org/10.1190/geo2015-0528.1 Google Scholar CrossRef Search ADS   Dafni R., Symes W.W., 2016b. Kinematic artifacts in the subsurface-offset extended image and their elimination by a dip-domain specularity filter, Geophysics , 81( 6), S477– S495. https://doi.org/10.1190/geo2016-0115.1 Google Scholar CrossRef Search ADS   de Bruin C., Wapenaar C., Berkhout A., 1990. Angle-dependent reflectivity by means of prestack migration, Geophysics , 55( 9), 1223– 1234. https://doi.org/10.1190/1.1442938 Google Scholar CrossRef Search ADS   Hou J., Symes W.W., 2015. An approximate inverse to the extended Born modeling operator, Geophysics , 80( 6), R331– R349. https://doi.org/10.1190/geo2014-0592.1 Google Scholar CrossRef Search ADS   Huang Y., Nammour R., Symes W.W., 2016. Flexibly preconditioned extended least-squares migration in shot-record domain, Geophysics , 81( 5), S299– S315. https://doi.org/10.1190/geo2016-0023.1 Google Scholar CrossRef Search ADS   Koren Z., Ravve I., 2011. Full-azimuth subsurface angle domain wavefield decomposition and imaging Part I: directional and reflection image gathers, Geophysics , 76( 1), S1– S13. https://doi.org/10.1190/1.3511352 Google Scholar CrossRef Search ADS   Kosloff D., Sherwood J., Koren Z., Machet E., Falkovitz Y., 1996. Velocity and interface depth determination by tomography of depth migrated gathers, Geophysics , 61( 5), 1511– 1523. https://doi.org/10.1190/1.1444076 Google Scholar CrossRef Search ADS   Liu Z., Bleistein N., 1995. Migration velocity analysis: theory and an iterative algorithm, Geophysics , 60( 1), 142– 153. https://doi.org/10.1190/1.1443741 Google Scholar CrossRef Search ADS   Montel J.P., Lambare G., 2013. Wave equation angle domain common image gather asymptotic analysis, in 83rd Annual International Meeting , SEG, Expanded Abstracts, pp. 3757– 3761. Mosher C., Foster D., 2000. Common angle imaging condition for prestack depth migration, in 70th Annual International Meeting , SEG, Expanded Abstracts, pp. 830– 833. Prucha M., Biondi B., Symes W.W., 1999. Angle-domain common-image gathers by wave-equation migration, in 69th Annual International Meeting , SEG, Expanded Abstracts, pp. 824– 827. Ravve I., Koren Z., 2011. Full-azimuth subsurface angle domain wavefield decomposition and imaging: Part 2: local angle domain, Geophysics , 76( 2), S51– S64. https://doi.org/10.1190/1.3549742 Google Scholar CrossRef Search ADS   Rickett J.E., Sava P.C., 2002. Offset and angle-domain common image-point gathers for shot-profile migration, Geophysics , 67( 3), 883– 889. https://doi.org/10.1190/1.1484531 Google Scholar CrossRef Search ADS   Sava P., Biondi B., 2004. Wave-equation migration velocity analysis. I. Theory, Geophys. Prospect. , 52 6, 593– 606. https://doi.org/10.1111/j.1365-2478.2004.00447.x Google Scholar CrossRef Search ADS   Sava P., Fomel S., 2006. Time-shift imaging condition in seismic migration, Geophysics , 71( 6), S209– S217. https://doi.org/10.1190/1.2338824 Google Scholar CrossRef Search ADS   Sava P., Vasconcelos I., 2011. Extended imaging conditions for wave-equation migration, Geophys. Prospect. , 59( 1), 35– 55. https://doi.org/10.1111/j.1365-2478.2010.00888.x Google Scholar CrossRef Search ADS   Sava P.C., Fomel S., 2003. Angle-domain common-image gathers by wavefield continuation methods, Geophysics , 68( 3), 1065– 1074. https://doi.org/10.1190/1.1581078 Google Scholar CrossRef Search ADS   Shen P., Symes W.W., 2008. Automatic velocity analysis via shot profile migration, Geophysics , 73( 5), VE49– VE59. https://doi.org/10.1190/1.2972021 Google Scholar CrossRef Search ADS   Stolk C.C., de Hoop M.V., Symes W.W., 2009. Kinematics of shot-geophone migration, Geophysics , 74( 6), WCA19– WCA34. https://doi.org/10.1190/1.3256285 Google Scholar CrossRef Search ADS   Symes W.W., 2008. Migration velocity analysis and waveform inversion, Geophys. Prospect. , 56( 6), 765– 790. https://doi.org/10.1111/j.1365-2478.2008.00698.x Google Scholar CrossRef Search ADS   Symes W.W., Carazzone J.J., 1991. Velocity inversion by differential semblance optimization, Geophysics , 56( 5), 654– 663. https://doi.org/10.1190/1.1443082 Google Scholar CrossRef Search ADS   Vyas M., Du X., Mobley E., Fletcher R., 2011. Methods for computing angle gathers using RTM, in 73rd Annual International Conference and Exhibition , EAGE, Extended Abstracts, F020. Wu R.S., Chen L., 2006. Directional illumination analysis using beamlet decomposition and propagation, Geophysics , 71( 4), S147– S159. https://doi.org/10.1190/1.2204963 Google Scholar CrossRef Search ADS   Xu S., Chauris H., Lambare G., Noble M., 2001. Common-angle migration: a strategy for imaging complex media, Geophysics , 66( 6), 1877– 1894. https://doi.org/10.1190/1.1487131 Google Scholar CrossRef Search ADS   Yoon K., Marfurt K.J., 2006. Reverse-time migration using the Poynting vector, Explor. Geophys. , 37( 1), 102– 107. https://doi.org/10.1071/EG06102 Google Scholar CrossRef Search ADS   APPENDIX A: SUBSURFACE OFFSET AND ANGLE-DOMAIN IMAGE ENVELOPES An envelope of one-parameter family of surfaces T(z, x, y; σ), in a given plane, is a curve that is tangent to each member of the family at some contact point. Each member of the family T is subjected to the characteristic parameter σ. The envelope of such family is formed at the geometric place where its members are stationary with respect to a change in σ. This defines the envelope's contact condition:   \begin{equation}\frac{{\partial T(z,x,y;\sigma )}}{{\partial \sigma }} = 0.\end{equation} (A1) We regard seismic migration as a superposition of two-way traveltime isochrons, and the resulting image as the envelope of such isochrons. Eq. (29) in the body of this paper describes a family of isochrons in the extended subsurface offset domain   \begin{eqnarray} T(z,x,h;H)\,{ =}\, \frac{{{{(x - {x_m})}^2}}}{{{\varepsilon ^2}(z_0^2 + {H^2})}}\,{ +}\, \frac{{{z^2}}}{{{\varepsilon ^2}(z_0^2 + {H^2}) - {{(h - H)}^2}}}\,{ -}\, 1\,{ =}\, 0.\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\nonumber\\ \end{eqnarray} (A2) Eq. (31) is another family of isochrons, defined in the angle-domain as the scattering-angle variant of eq. (A2):   \begin{eqnarray}T(z,x,p;H) &=& z - \sqrt {({\varepsilon ^2}(z_0^2 + {H^2}) - {{(x - {x_m})}^2})(1 + {p^2}{\eta ^2})} \nonumber \\ && +\, pH = 0.\end{eqnarray} (A3) Each member in these two families is subjected to a characteristic acquisition offset value H. For simplicity, we assume that a horizontal reflector in under study. Therefore, the image point to midpoint distance vanishes: x − xm = 0, and the stretch factor minimizes to η = 1. Hence, the problem becomes laterally invariant, such that the described families are independent of the horizontal coordinate x. They take the form   \begin{equation}T(z,h;H) = {z^2} + {(h - H)^2} - {\varepsilon ^2}(z_0^2 + {H^2}) = 0,\end{equation} (A4)  \begin{equation}T(z,p;H) = z - \sqrt {{\varepsilon ^2}(z_0^2 + {H^2})(1 + {p^2})} + pH = 0.\end{equation} (A5) Applying the contact condition of eq. (A1) results in the following contact relations   \begin{equation}h = - H({\varepsilon ^2} - 1),\end{equation} (A6)  \begin{equation}p = \frac{{\varepsilon H}}{{\sqrt {z_0^2 - {H^2}({\varepsilon ^2} - 1)} }}.\end{equation} (A7) These relations give the contact point of each family member, characterized by the offset H, on the subsurface offset and the scattering-angle image envelopes, respectively. Substituting eqs (A6) and (A7) back into eqs (A4) and (A5) to eliminate H, yields the corresponding envelope expressions   \begin{equation}z(h) = \varepsilon \sqrt {z_0^2 - \frac{{{h^2}}}{{({\varepsilon ^2} - 1)}}} ,\end{equation} (A8)  \begin{equation}z(p) = {z_0}\sqrt {{\varepsilon ^2} + {p^2}({\varepsilon ^2} - 1)} .\end{equation} (A9) APPENDIX B: CONTACT CONDITION FOR SUBSURFACE OFFSET EXTENDED IMAGING The envelope theorem in Appendix  A explains the construction of the subsurface offset extended image as an envelope of isochron curves. It is established on a contact condition that amounts to a contact relation between the offsets h and H in eq. (A6). This contact condition is linked to the SRP kinematic description, and in some sense equivalent to the phase space behaviour of the waves. We discovered that under the assumption of horizontally layered media the contact condition asymptotically imposes a consistency in the horizontal slowness of the incident and scattered waves. The horizontal slowness is kept constant regardless of the migration velocity accuracy   \begin{equation}\begin{array}{@{}l@{}} {\nabla _x}{t_r}(V) = {\nabla _x}{t_r}({V_m})\\ {\nabla _x}{t_s}(V) = {\nabla _x}{t_s}({V_m}). \end{array}\end{equation} (B1) It is obtained for a single-layer model by substituting eq. (A6) contact relation into eq. (18) horizontal slownesses. Moreover, using Snell's law, it is straightforward to verify its validity for a general multilayer model. In the following, eq. (A6) contact relation is generalized with respect to a multilayer subsurface model. Considering a horizontally layered model and a reflection from the ith layer interface, we use the concept of datuming to extrapolate the data traveltimes from the acquisition surface downward to the layer interface i – 1 (see Fig. 7c). At this datum level, eq. (A4) family of isochron curves is given by   \begin{eqnarray}{\left. {T(z,h;H)} \right|_{{z_{i - 1}}}} = {(z - {z_{i - 1}})^2} + {(h - {H_{i - 1}})^2} - \varepsilon _i^2{({V_i}\Delta t)^2} = 0, \nonumber \\ \end{eqnarray} (B2)where Hi–1 is the datumed offset at the datum level zi–1. The third term in the equation asymptotically determines the residual length of the incident/scattered ray (both share a symmetric ray trajectory) from the datum level to the reflection point. The residual time Δt is the sum over the traveltime error down to the datum level plus the remaining traveltime in the ith layer   \begin{equation}\Delta t = \sum\limits_{k = 1}^{i - 1} {\left[ {{t_k}(V) - {t_k}({V_m})} \right]} + {t_i}({V_m}).\end{equation} (B3) We employ the fundamental envelope condition of eq. (A1) with respect to the datumed family of isochrons in eq. (B2) to obtain   \begin{equation}\frac{{\partial {{\left. {T(z,h;H)} \right|}_{{z_{i - 1}}}}}}{{\partial H}} = 2(h - {H_{i - 1}}) + 2\varepsilon _i^2V_i^2\Delta t\left( {\frac{{\partial \Delta t}}{{\partial H}}} \right) = 0,\end{equation} (B4)whereas from eq. (B3) we get   \begin{equation}\frac{{\partial \Delta t}}{{\partial H}} = \sum\limits_{k = 1}^{i - 1} {\left[ {{\nabla _x}{t_k}(V) - {\nabla _x}{t_k}({V_m})} \right] + {\nabla _x}{t_i}({V_m})} .\end{equation} (B5) Incorporating eq. (B1) horizontal slowness consistency eliminates the sum term in eq. (B5). Substituting eq. (B5) into eq. (B4) and rearranging the terms yields the generalized contact relation in a multilayer media   \begin{equation}h = {H_{i - 1}} - {\varepsilon _i}{V_i}\Delta t\sin \gamma .\end{equation} (B6) The equation coincides with eq. (A6) when i = 1 (i.e. a single layer model). The generalized contact relation gives the contact point of each member in the family of isochrons on the subsurface offset domain image envelope. As demonstrated in the body of this paper, it may be exploited to predict the defocusing pattern of the image in subsurface offset CIGs. © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society.
TI - Kinematics of reflections in subsurface offset and angle-domain image gathers
JF - Geophysical Journal International
DO - 10.1093/gji/ggy051
DA - 2018-05-01
UR - https://www.deepdyve.com/lp/oxford-university-press/kinematics-of-reflections-in-subsurface-offset-and-angle-domain-image-6W9Ln9Jf8O
SP - 1212
EP - 1230
VL - 213
IS - 2
DP - DeepDyve
ER -