CRS theory based on a data-driven homeomorphism

The main objective of this article is to present a new path to a solution for the common-reflection-surface (CRS) formula, which approximates, up to the second order in the midpoint-offset domain $({{\bf{x}}}_{m},h)$ of all source–receiver pairs, the time of flight of seismic impulses travelling along rays that, once reflected, carry to the surface the echo of the subsurface layering. The problem is first solved in an auxiliary homogeneous medium. Then, and this is the novelty of the present study, the resulting time of flight is synchronised, regardless of any earth velocity model, via a 3D homeomorphic transform, to match the signal traveltime in the prospected medium.

This result is made possible by first approximating the isochrones tangent to a small reflecting element embedded in a 3D isotropic homogeneous auxiliary medium with velocity v₀. The solution, parametrised by six geometric attributes, may be derived from the approximated primary response of the reflecting surface around P_NIP, with P_NIP the incidence point of the normal reference straight ray that emerges at the surface position ${{\bf{x}}}_{0}$ with elevation angle ${\alpha }_{0}$ and azimuth ${\beta }_{0}$ after a traveltime ${({v}_{0}{\kappa }_{{\rm{NIP}}})}^{-1}$. The derivation is solely based on Snell's law and ray geometry. The consequent paraxial traveltime requires the use of ${\kappa }_{{\rm{NIP}}}$, the curvature of the spherical diffraction wavefront leaving P_NIP, and of ${{\bf{K}}}_{N}$, the symmetric $(2\times 2)$-matrix related to the curvature of the reflector at the normal-incidence point P_NIP.

Then, in the case of a layered isotropic medium, the assignment to the diffraction front of one value of ${\kappa }_{{\rm{NIP}}}$ for each direction around ${{\bf{x}}}_{0}$ describes locally the front emerging with elevation angle ${\alpha }_{0}$ and azimuth ${\beta }_{0}$ after a traveltime t₀. By means of the symmetric $(2\times 2)$-curvature matrix ${{\bf{K}}}_{{\rm{NIP}}}$, the correspondence may be now established between reflecting elements of the two media, the real and the auxiliary, which respond, however, with different traveltimes. This 3D homemorphism provides the time correction to the traveltime in the auxiliary medium [4], giving rise to the well-known paraxial CRS traveltime formula $t({{\bf{x}}}_{m}-{{\bf{x}}}_{0},{\bf{h}})$ with eight parameters.

The approach here presented deeply differs from the one introduced by Bortfeld in his seminal work [8], where the traveltime near a central ray appears from the approximation of Hamilton's equation for rays, whose solution, given in terms of generic $(2\times 2)$-matrices, is based on reflection–transmission rules in a stratified medium. This result was later specialised by the German school of Karlsruhe [9] and recently summarised by Vanelle [29]. Afterward, in a joint venture with a major oil company [1, 10, 11], the theory became a mature numerical technology for the production of 3D high-resolution images from very large collections of seismic traces (up to ten Tbytes of input data). This application has turned out to be an fundamental tool [6, 23, 24, 33], especially for imaging seismic data with a weak signal-to-noise ratio, from areas characterised by highly fractured reflectors (crystalline or hardrock formations) and afflicted by a poor seismic response due to low impedance contrasts.

The inspiration of this study comes from a second seminal work written by Höcht, de Bazelaire, Majer and Hubral [16] on CRS theory, where the homeomorphic transform is restricted exclusively to the 2D case [3]. For about ten years, this promising methodology was stuck by the difficulty to construct in 3D the adequate correspondence between the real medium and the auxiliary one.

The correspondence illustrated in the present article between 3D quasi-similar media, besides providing the desired hyperbolic traveltime, opens a field of applications in time imaging that is no longer based on reflections but, more generally, on diffractions. This second case gives rise in 3D to a source-to-receiver paraxial traveltime in the form of a double-square root with only five geometric parameters.

It is worth mentioning a concurrent alternative to CRS, the multifocusing method developed early in the 1990s by Gelchisnky [14, 15], which applies a similar strategy in 2D as Höcht et al [16], using the same definition of auxiliary medium and the consequent traveltime synchronisation. However, in this case, there are even more difficulties for generalising the development of a quasi-similar earth model beyond the 2D problem.

Moreover, CRS theory in 3D yields a traveltime with eight unknown parameters or attributes, each one with a specific geometrical meaning, namely, two angles, ${\alpha }_{0}$ and ${\beta }_{0}$, and the six independent matrix elements of ${{\bf{K}}}_{N}$ and ${{\bf{K}}}_{{\rm{NIP}}}$, all associated to the pair $({{\bf{x}}}_{0},{t}_{0})$. For each 8-tuple, the adequacy of its numerical values to seismic data is measured by the semblance, the highly nonlinear multimodal objective function to be maximized. The optimisation phase, which is only sketched in this article, is obtained by shaping the CRS traveltime hypersurface so as to intercept the largest number of coherent signals in the volume of seismic traces gathered in the midpoint-offset domain.

With optimised attributes, a large number of correlated signals can be collected along the hypersurface profile and stacked, providing the effective reflectivity associated to a pixel of coordinate $({{\bf{x}}}_{0},{t}_{0})$ within the image volume (totally, a few Gbytes of output). Since the whole process is based on the mining of a massive number of coherent seismic signals of the midpoint-offset domain, the resolution method, which is correctly said data-driven, increases drastically the signal-to-noise ratio and the focusing of the medium layering.

Section 2 of this document introduces the concepts of seismic data stacking and seismic migration and the role of CRS processing. In section 3, before tackling the general case, the CRS paraxial traveltime approximation is derived from ab initio ray geometry in a 3D homogeneous medium, realising an important conceptual step. Section 4 presents the genuine 3D homeomorphism that yields, after the synchronisation of the traveltimes, the well-known 3D CRS expression for stratified media [21, 29]. Section 5 sketches the data-driven strategy for the numerical estimate of the eight CRS attributes. Three successful numerical applications of CRS processing are illustrated in section 6. Some conclusion are drowned in section 7, outlining the natural extension of this work.

The construction of a seismic image is a key issue in geophysical exploration. In principle, to obtain such an image, it is necessary to solve the wave equation in reverse time and reconstruct, from the gather of seismic traces, the sources of reflection and diffraction beneath the surface. This task is known as seismic data migration [2]. Obviously, to back-propagate the wave equation and focus seismic signals directly in depth, it is fundamental to know the medium seismic velocity, a (generally increasing) function of depth whose construction continue to be the main difficult and time-consuming problem of seismic imaging.

To obtain a seismic velocity model, the typical approach is to recover, from the collection of traces in the common-midpoint gathers, the traveltimes corresponding to the same subsurface point. From these traveltimes and their difference, it is possible to reconstruct an apparent velocity field, or stacking velocity, corresponding to the value of the seismic velocity obtained from the best fit of the traveltime curve by a hyperbola [32]. From this data redundancy, it is also possible to obtain another benefit: actually, stacking, or summing together the signals referring to the same event sensed at different surface positions, allows one to emulate a zero-offset source–receiver recording with a strongly improved signal-to-noise ratio. To form the pixel relative to this event, the summation across traces is performed along coherent patterns of signals with times of flight that must be properly modelled. Since these patterns are weakly dependent on seismic velocity, seismic data inversion is then by its nature a hard problem [30, 31]. To improve the reliability of the velocity model building techniques, service companies usually devise a processing chain to estimate a seismic velocity portrait at the cost of an intensive man–machine interaction [32]. This procedure is often time consuming and sometimes not accurate enough to give useful results, especially when the signal-to-noise ratio is too low and/or the target structure is complex. It is in this framework that the CRS traveltime formula plays its most important role.

Imaging strategies can be summarised in four categories: pre-stack time migration, post-stack time migration, pre-stack depth migration and post-stack depth migration [2, 32]. The word time or depth in the previous classification refers to the vertical axis of the output image. Of course, input data are recorded in time units, milliseconds, and at different positions with horizontal axes defined in length units, meters. However, the vertical dimension may be expressed in time units or, using a seismic velocity model, in length units. In the case of a good knowledge of the seismic velocity field, a situation that will not be discussed in this article, the optimal processing choice is to perform depth imaging having in mind that post-stack depth migration only has the advantage to require, with respect to pre-stack depth migration, a considerable reduced amount of computation but at cost of a less resolved image.

The situation is radically different with a noisy acquisition for which it is almost impossible to build a reliable seismic velocity model. With noisy data, in the case of layered media with mild lateral velocity variations and many faults with diffractors, CRS stacking becomes the optimal tool to improve the signal-to-noise ratio by reducing common midpoint data to zero-offset. As shown in the daily practice, section 6, a stacked image that is properly focused can be transformed in an interpretable migrated image with an acceptable distortion since the quality of the final time migrated volume may override the errors of the velocity field.

In 3D, CRS technique has three main advantages with respect to conventional stacking [32]. First, CRS processing does not require a velocity model but it can recognise coherent signal patterns via the optimisation method that will be sketched in sections 5. Second, the control of the Fresnel zone projected on the horizontal plane [19], which is fully exploited only by CRS processing, allows one to select, in quality and quantity around ${{\bf{x}}}_{0}$, gathers of seismic traces within a region of the midpoint-offset domain compatible with the CRS hyperbolic approximation. These traces enter in large numbers into the stack procedure, leading to a dramatic increase of the signal-to-noise ratio of the final result. The third advantage regards the level of details that this technique allows one to reach since every single pixel value of the 3D volume is independently optimised, while in the conventional stacking the normal move-out correction depends on a velocity field interpolated from a reduced number of control points [32].

Of course, CRS processing has one important restriction in media with a response too far from the traveltime hyperbolic approximation. This happens in presence of strong lateral velocity variations. In this case, the method collapses, thus requiring a further conceptual refinement that lead, for instance, to a double-square root formulation of the traveltime. A promising solution is a common-diffraction-point theory based on a data-driven homeomorphism [5, 7].

The trajectory of a ray in a homogenous medium, moving downward from the datum plane to a reflector and then upward, is illustrated in figure 1. According to Snell's law, in a homogeneous medium the normal ${\bf{n}}$ to the surface element at a reflection point P lies on the plane defined by P, the source S and the receiver G. As a consequence, the normal incidence ray, or NIP-ray [18], moving along ${\bf{n}}$ intersects the acquisition line spanned by vector ${\bf{u}}$. Let us denote by ${\bf{x}}$ the coordinate of the intersection point.

Zoom In Zoom Out Reset image size

Provided ${R}_{{\rm{NIP}}}({\bf{n}})$, the length of the normal incidence trajectory, for each half offset ${\bf{h}}=h{\bf{u}}$ of the source–receiver setup, one wants to determine the acquisition midpoint ${{\bf{x}}}_{m}$ and the traveltime t, two quantities defining together the isochrone ellipse tangent to the reflector at the normal incidence point P:

$$\begin{eqnarray}&&t=\displaystyle \frac{\overline{{PS}}+\overline{{GP}}}{{v}_{0}}.\end{eqnarray} \tag{ 1 }$$

The problem is intrinsically 2D, so, embedding its 2D solution [17] in 3D, one finds the following two fundamental exact relations:

$$\begin{eqnarray}{t}^{2}({\bf{n}},{\bf{h}}) & = & {\left(\displaystyle \frac{2h}{{v}_{0}}\right)}^{2}+2{\left[\displaystyle \frac{{R}_{{\rm{NIP}}}({\bf{n}})}{{v}_{0}}\right]}^{2}\\ & & \times \left\{\sqrt{{\left[\displaystyle \frac{h}{{\rm{\Gamma }}({\bf{n}})}\right]}^{2}+1}+1\right\},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{{\bf{x}}}_{m}({\bf{n}},{\bf{h}})={\bf{x}}({\bf{n}})-{\rm{\Gamma }}({\bf{n}})\left\{\sqrt{{\left[\displaystyle \frac{h}{{\rm{\Gamma }}({\bf{n}})}\right]}^{2}+1}-1\right\}{\bf{u}},\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&{\rm{\Gamma }}({\bf{n}})=\displaystyle \frac{{R}_{{\rm{NIP}}}({\bf{n}})}{2\ {\bf{n}}\cdot {\bf{u}}}=\displaystyle \frac{{R}_{{\rm{NIP}}}({\bf{n}})}{2\sin \alpha ({u}_{1}\cos \beta +{u}_{2}\sin \beta )},\\ &&\quad \parallel {\bf{n}}\parallel =\parallel {\bf{u}}\parallel =1,\end{eqnarray} \tag{ 4 }$$

v₀ is the ray velocity and α is the elevation angle of ${\bf{n}}$ with respect to the vertical $\overrightarrow{{oz}}$, while β is its azimuth with respect to $\overrightarrow{{ox}}$.

The objective is to compute the source-to-receiver traveltime t as a function of ${{\bf{x}}}_{m}$ and ${\bf{h}}$. This is possible by eliminating ${\bf{n}}$ between equations (2) and (3), although the resulting system is clearly highly nonlinear. The idea is to develop a paraxial approximation of t with respect to a normal-incident (NIP) ray of reference, which is characterised by an exit direction ${{\bf{n}}}_{0}$.

As shown in figure 2, the reference NIP-ray leaves the datum plane downward at ${{\bf{x}}}_{0}={\bf{x}}({{\bf{n}}}_{0})$ and, in a two-way traveltime t₀, reaches the buried reflecting surface at the reference normal-incidence point P₀ and then goes back along the same path. Since around P₀ the reflecting surface is locally quadratic, for each point P close enough to P₀, one can state there is a plane, formed by the intersection of ${{\bf{n}}}_{0}$ and the normal ray direction ${\bf{n}}$, cutting the reflector along a curve ${\gamma }_{P}$. This curve is approximately an arc of circle centred at C with radius of curvature ${R}_{R}^{{\gamma }_{P}}$ along the normal direction ${{\bf{n}}}_{0}$. Let us denote by ${R}_{N}({\bf{n}})$, see figure 2, the length of the segment joining along direction ${\bf{n}}$ the centre of curvature C of ${\gamma }_{P}$ to the datum plane at point ${\bf{x}}({\bf{n}})$:

$$\begin{eqnarray}&&{R}_{N}({\bf{n}})={R}_{R}^{{\gamma }_{P}}+{R}_{{\rm{NIP}}}({\bf{n}}).\end{eqnarray} \tag{ 5 }$$

Since ${R}_{R}^{{\gamma }_{{P}_{0}}}\approx {R}_{R}^{{\gamma }_{P}}$ holds for small displacements from P₀ along ${\gamma }_{P}$, within this approximation it easy to see that

$$\begin{eqnarray}&&{R}_{N}({\bf{n}})\approx {R}_{N}^{{\gamma }_{P}}\displaystyle \frac{\cos {\alpha }_{0}}{\cos \alpha },\hspace{0.5cm}{R}_{N}^{{\gamma }_{P}}={R}_{N}({{\bf{n}}}_{0}).\end{eqnarray} \tag{ 6 }$$

Under the same hypothesis, the merge of this last expression and equation (5) gives rise to the following relation

$$\begin{eqnarray}&&{R}_{{\rm{NIP}}}({\bf{n}})\approx {R}_{{\rm{NIP}}}+{R}_{N}^{{\gamma }_{P}}\left(\displaystyle \frac{\cos {\alpha }_{0}}{\cos \alpha }-1\right),\\ &&\quad {R}_{{\rm{NIP}}}={R}_{{\rm{NIP}}}({{\bf{n}}}_{0}).\end{eqnarray} \tag{ 7 }$$

Next, given ${{\bf{x}}}_{0}$, ${{\bf{n}}}_{0}$ and ${R}_{N}^{{\gamma }_{P}}={R}_{R}^{{\gamma }_{P}}+{R}_{{\rm{NIP}}}$, three quantities characterising the reference NIP-ray, the goal is to estimate ${\bf{x}}({\bf{n}})$ for P close to P₀ . To start, it is necessary to establish the geometric relation between the center of curvature C and ${\bf{x}}({\bf{n}})$. That is

$$\begin{eqnarray}&&{\bf{x}}({\bf{n}})-{{\bf{x}}}_{0}\approx {\widehat{{\bf{x}}}}_{C}({{\bf{n}}}_{0})+\displaystyle \frac{{z}_{C}({{\bf{n}}}_{0})}{{n}_{3}}\ \widehat{{\bf{n}}},\end{eqnarray} \tag{ 8 }$$

where ${\widehat{{\bf{x}}}}_{C}$ and $\widehat{{\bf{n}}}$ denote the projected coordinates onto the datum plane:

$$\begin{eqnarray}&&{\widehat{{\bf{x}}}}_{C}({{\bf{n}}}_{0})=-{R}_{N}^{{\gamma }_{P}}\sin {\alpha }_{0}\left[\begin{array}{l}\cos {\beta }_{0}\\ \sin {\beta }_{0}\end{array}\right],\\ &&\quad {z}_{C}({{\bf{n}}}_{0})={R}_{N}^{{\gamma }_{P}}\cos {\alpha }_{0},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\widehat{{\bf{n}}}=\sin \alpha \left[\begin{array}{l}\cos \beta \\ \sin \beta \end{array}\right],\hspace{.3cm}{n}_{3}=\cos \alpha .\end{eqnarray} \tag{ 10 }$$

Substituting equations (9) and (10) into (8), one obtains the desired estimate:

$$\begin{eqnarray}\displaystyle \frac{{\bf{x}}({\bf{n}})-{{\bf{x}}}_{0}}{{R}_{N}^{{\gamma }_{P}}} & \approx & -\sin {\alpha }_{0}\left[\begin{array}{l}\cos {\beta }_{0}\\ \sin {\beta }_{0}\end{array}\right]\\ & & +\tan \alpha \cos {\alpha }_{0}\left[\begin{array}{l}\cos \beta \\ \sin \beta \end{array}\right].\end{eqnarray} \tag{ 11 }$$

To summarise, the two auxiliary equations (7) and (11) provide the paraxial description of the NIP-ray associated to ${\bf{n}}$ near the reference NIP trajectory. With equations (2)–(4), all the ingredients are now established to construct the second-order paraxial approximation of the traveltime $t({{\bf{x}}}_{m}-{{\bf{x}}}_{0},{\bf{h}})$, for ${{\bf{x}}}_{m}$ close to ${{\bf{x}}}_{0}$ and small acquisition offsets h. Figure 3 illustrates the geometric setup.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Since the reciprocity principle, when ${\bf{h}}\to -{\bf{h}}$, preserves the source–receiver traveltime, the power-series expansion of t, up to the second order, takes the form

$$\begin{eqnarray}&&t({{\bf{x}}}_{m}-{{\bf{x}}}_{0},{\bf{h}})=\displaystyle \frac{2{R}_{{\rm{NIP}}}}{{v}_{0}}+{{\bf{a}}}^{{\mathsf{T}}}({{\bf{x}}}_{m}-{{\bf{x}}}_{0})\\ &&\quad {+({{\bf{x}}}_{m}-{{\bf{x}}}_{0})}^{{\mathsf{T}}}{{\bf{K}}}_{m}({{\bf{x}}}_{m}-{{\bf{x}}}_{0})+{{\bf{h}}}^{{\mathsf{T}}}{{\bf{K}}}_{h}{\bf{h}}+O(3).\end{eqnarray} \tag{ 12 }$$

The consequence is that ${{\bf{K}}}_{m}$ and ${{\bf{K}}}_{h}$ can be calculated by setting separately ${\bf{h}}=0$ and ${{\bf{x}}}_{m}={{\bf{x}}}_{0}$ into the traveltime and the midpoint equations (2) and (3). Figure 4 displays the acquisition setup corresponding to these two source–receiver configurations.

Zoom In Zoom Out Reset image size

With equations (3) and (11), the new objective is to calculate, up to the second order, the power series of $\alpha \,=\alpha ({{\bf{x}}}_{m}-{{\bf{x}}}_{0},{\bf{h}})$ and $\beta =\beta ({{\bf{x}}}_{m}-{{\bf{x}}}_{0},{\bf{h}})$, with ${\alpha }_{0}=\alpha (0,0)$ and ${\beta }_{0}=\beta (0,0)$.

By setting ${\bf{h}}={\bf{0}}$ in the auxiliary vector equation (3) and subsequently ${\bf{x}}({\bf{n}})={{\bf{x}}}_{m}({\bf{n}},0)$ in (11), the zero-offset problem gives rise to the nonlinear system

$$\begin{eqnarray}{{\bf{x}}}_{m}({\bf{n}},0)-{{\bf{x}}}_{0} & = & {R}_{N}^{{\gamma }_{P}}\left\{-\sin {\alpha }_{0}\left[\begin{array}{l}\cos {\beta }_{0}\\ \sin {\beta }_{0}\end{array}\right]\right.\\ & & +\left.\tan \alpha \cos {\alpha }_{0}\left[\begin{array}{l}\cos \beta \\ \sin \beta \end{array}\right]\right\}.\end{eqnarray} \tag{ 13 }$$

Then, under the condition ${{\bf{x}}}_{m}({\bf{n}},h)={{\bf{x}}}_{0}$, merging equations (3) and (11) results in the following nonlinear system

$$\begin{eqnarray}&&\left\{\sqrt{{\left[\displaystyle \frac{h}{{\rm{\Gamma }}({\bf{n}})}\right]}^{2}+1}-1\right\}{\bf{u}}\\ &&\quad =\displaystyle \frac{{R}_{N}^{{\gamma }_{P}}}{{\rm{\Gamma }}({\bf{n}})}\left\{-\sin {\alpha }_{0}\left[\begin{array}{l}\cos {\beta }_{0}\\ \sin {\beta }_{0}\end{array}\right]+\tan \alpha \cos {\alpha }_{0}\left[\begin{array}{l}\cos \beta \\ \sin \beta \end{array}\right]\right\}.\end{eqnarray} \tag{ 14 }$$

Both equations (13) and (14) must be expanded in α and β around ${\alpha }_{0}$ and ${\beta }_{0}$, respectively, and then the resulting series must be inverted up to the second order, the first in ${{\bf{x}}}_{m}-{{\bf{x}}}_{0}$ and the second in ${\bf{h}}$.

3.1. Zero-offset acquisition

Expanding in power series the right-end side of (13), to the second order one obtains:

$$\begin{eqnarray}\displaystyle \frac{{{\bf{x}}}_{m}({\bf{n}},0)-{{\bf{x}}}_{0}}{{R}_{N}^{{\gamma }_{P}}} & = & \left[\begin{array}{c}\cos {\beta }_{0}\\ \sin {\beta }_{0}\end{array}\right]\displaystyle \frac{\alpha -{\alpha }_{0}}{\cos {\alpha }_{0}}\\ & & +\,\left[\begin{array}{c}-\sin {\beta }_{0}\\ \cos {\beta }_{0}\end{array}\right]\sin {\alpha }_{0}\ (\beta -{\beta }_{0})\\ & & +\,\left[\begin{array}{c}\cos {\beta }_{0}\\ \sin {\beta }_{0}\end{array}\right]\displaystyle \frac{\sin {\alpha }_{0}}{{\cos }^{2}{\alpha }_{0}}\ {(\alpha -{\alpha }_{0})}^{2}\\ & & +\,\left[\begin{array}{c}-\sin {\beta }_{0}\\ \cos {\beta }_{0}\end{array}\right]\displaystyle \frac{(\alpha -{\alpha }_{0})(\beta -{\beta }_{0})}{\cos {\alpha }_{0}}\\ & & -\,\left[\begin{array}{c}\cos {\beta }_{0}\\ \sin {\beta }_{0}\end{array}\right]\displaystyle \frac{\sin {\alpha }_{0}}{2}\ {(\beta -{\beta }_{0})}^{2}\,+\,O(3),\end{eqnarray} \tag{ 15 }$$

a system of two equations in α and β whose inversion relies on the following result.

• The inverse function theorem: suppose the dependence between the variables $(\alpha ,\beta )$ and (x, y) is implicitly defined by two functions, analytic at $(0,0)$, of the form

  $$\begin{eqnarray*}&&x=f(\alpha ,\beta ),\hspace{0.3cm}y=g(\alpha ,\beta ),\end{eqnarray*}$$

  with $f(0,0)=g(0,0)=0$ and $(\partial f/\partial \alpha )(0,0)$ $(\partial g/\partial \beta )(0,0)$ $\ne (\partial f/\partial \beta )(0,0)$ $(\partial g/\partial \alpha )(0,0)$. Then it is possible to invert or solve the system of equations for $(\alpha ,\beta )$, where $\alpha =F(x,y)$ and $\beta =G(x,y)$ with F and G analytic at the point $(0,0)$.

• Inversion procedure: because f and g are analytic at $(0,0)$, one can expand x and y in Taylor series in α and β:

  $$\begin{eqnarray}&&\left[\begin{array}{c}x\\ y\end{array}\right] & = & \left[\begin{array}{c}{f}_{\mathrm{1,0}}\\ {g}_{\mathrm{1,0}}\end{array}\right]\alpha +\left[\begin{array}{c}{f}_{\mathrm{0,1}}\\ {g}_{\mathrm{0,1}}\end{array}\right]\beta +\left[\begin{array}{c}{f}_{\mathrm{2,0}}\\ {g}_{\mathrm{2,0}}\end{array}\right]{\alpha }^{2}\\ & & +\left[\begin{array}{c}{f}_{\mathrm{1,1}}\\ {g}_{\mathrm{1,1}}\end{array}\right]\alpha \beta +\left[\begin{array}{c}{f}_{\mathrm{0,2}}\\ {g}_{\mathrm{0,2}}\end{array}\right]{\beta }^{2}+O(3).\end{eqnarray} \tag{ 16 }$$

  In virtue of the theorem, the series expansions of α and β exist and take the general form:

  $$\begin{eqnarray}&&\left[\begin{array}{c}\alpha \\ \beta \end{array}\right] & = & \left[\begin{array}{c}{F}_{\mathrm{1,0}}\\ {G}_{\mathrm{1,0}}\end{array}\right]x+\left[\begin{array}{c}{F}_{\mathrm{0,1}}\\ {G}_{\mathrm{0,1}}\end{array}\right]y+\left[\begin{array}{c}{F}_{\mathrm{2,0}}\\ {G}_{\mathrm{2,0}}\end{array}\right]{x}^{2}\\ & & +\left[\begin{array}{c}{F}_{\mathrm{1,1}}\\ {G}_{\mathrm{1,1}}\end{array}\right]{xy}+\left[\begin{array}{c}{F}_{\mathrm{0,2}}\\ {G}_{\mathrm{0,2}}\end{array}\right]{y}^{2}+O(3).\end{eqnarray} \tag{ 17 }$$

  After the substitution of (17) into (16), the right-hand side of (16) becomes a polynomial of ${x}^{n}{y}^{m}$. Comparing coefficient by coefficient both sides, one first finds the following two linear systems

  $$\begin{eqnarray*}{f}_{\mathrm{1,0}}\ {F}_{\mathrm{1,0}}+{f}_{\mathrm{1,0}}\ {G}_{\mathrm{1,0}} & = & 1\\ {f}_{\mathrm{1,0}}\ {F}_{\mathrm{0,1}}+{f}_{\mathrm{1,0}}\ {G}_{\mathrm{0,1}} & = & 0,\\ {g}_{\mathrm{1,0}}\ {F}_{\mathrm{1,0}}+{g}_{\mathrm{0,1}}\ {G}_{\mathrm{1,0}} & = & 0\\ {g}_{\mathrm{1,0}}\ {F}_{\mathrm{0,1}}+{g}_{\mathrm{0,1}}\ {G}_{\mathrm{0,1}} & = & 1,\end{eqnarray*}$$

  whose solution in $\{{F}_{\mathrm{1,0}},{G}_{\mathrm{1,0}}\}$ and $\{{F}_{\mathrm{0,1}},{G}_{\mathrm{0,1}}\}$ is then used to solve the sequence of linear systems in $\{{F}_{\mathrm{2,0}},{G}_{\mathrm{2,0}}\}$, $\{{F}_{\mathrm{0,2}},{G}_{\mathrm{0,2}}\}$ and $\{{F}_{\mathrm{1,1}},{G}_{\mathrm{1,1}}\}$, and so on.

Using Maple to solve symbolically the cumbersome inversion of (15), the implementation of the above procedure to the second order provides the paraxial change of elevation and azimuth for small variations of ${{\bf{x}}}_{m}$ near ${{\bf{x}}}_{0}$. The resulting two expressions are new:

$$\begin{eqnarray}\alpha -{\alpha }_{0} & = & \displaystyle \frac{1}{\tan {\alpha }_{0}}\left[{\widehat{{\bf{n}}}}_{0}^{{\mathsf{T}}}\displaystyle \frac{({{\bf{x}}}_{m}-{{\bf{x}}}_{0})}{{R}_{N}^{{\gamma }_{P}}}\right.\\ & & \left.+\displaystyle \frac{1}{2}\displaystyle \frac{{({{\bf{x}}}_{m}-{{\bf{x}}}_{0})}^{{\mathsf{T}}}}{{R}_{N}^{{\gamma }_{P}}}\,{\bf{A}}\,\displaystyle \frac{({{\bf{x}}}_{m}-{{\bf{x}}}_{0})}{{R}_{N}^{{\gamma }_{P}}}\right]+O(3),\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}\beta -{\beta }_{0} & = & \displaystyle \frac{1}{\sin {\alpha }_{0}^{2}}\left[\displaystyle \frac{{\rm{d}}{\widehat{{\bf{n}}}}_{0}^{{\mathsf{T}}}}{{\rm{d}}{\beta }_{0}}\displaystyle \frac{({{\bf{x}}}_{m}-{{\bf{x}}}_{0})}{{R}_{N}^{{\gamma }_{P}}}\right.\\ & & +\left.\displaystyle \frac{1}{2}\displaystyle \frac{{({{\bf{x}}}_{m}-{{\bf{x}}}_{0})}^{{\mathsf{T}}}}{{R}_{N}^{{\gamma }_{P}}}\,{\bf{B}}\,\displaystyle \frac{({{\bf{x}}}_{m}-{{\bf{x}}}_{0})}{{R}_{N}^{{\gamma }_{P}}}\right]+O(3),\end{eqnarray} \tag{ 19 }$$

where ${\widehat{{\bf{n}}}}_{0}$, the projected reference normal vector, is defined by (10) while

$$\begin{eqnarray*}{\bf{A}}=\left[\begin{array}{ll}-3{\cos }^{2}{\beta }_{0}+2{\cos }^{2}{\beta }_{0}{\cos }^{2}{\alpha }_{0}+1 & -\cos {\beta }_{0}\sin {\beta }_{0}(1+2{\sin }^{2}{\alpha }_{0})\\ -\cos {\beta }_{0}\sin {\beta }_{0}(1+2{\sin }^{2}{\alpha }_{0}) & 3{\cos }^{2}{\beta }_{0}+2{\sin }^{2}{\beta }_{0}{\cos }^{2}{\alpha }_{0}-2\end{array}\right]\end{eqnarray*}$$

and

$$\begin{eqnarray*}{\bf{B}}=\left[\begin{array}{ll}2\cos {\beta }_{0}\sin {\beta }_{0} & \sin {\beta }_{0}-{\cos }^{2}{\beta }_{0}\\ \sin {\beta }_{0}-{\cos }^{2}{\beta }_{0} & -2\cos {\beta }_{0}\sin {\beta }_{0}\end{array}\right].\end{eqnarray*}$$

With the help of the auxiliary equation (7) and under the assumption ${\bf{h}}=0$, the traveltime expression (2) becomes

$$\begin{eqnarray}&&t({\bf{n}},0)=\displaystyle \frac{2{R}_{{\rm{NIP}}}({\bf{n}})}{{v}_{0}}\approx \displaystyle \frac{2}{{v}_{0}}\left[{R}_{{\rm{NIP}}}+{R}_{N}^{{\gamma }_{P}}\left(\displaystyle \frac{\cos {\alpha }_{0}}{\cos \alpha }-1\right)\right],\end{eqnarray} \tag{ 20 }$$

a function independent of β.

Expanding (20) in power series of $(\alpha -{\alpha }_{0})$, this equation takes the form

$$\begin{eqnarray}t({\bf{n}},0) & \approx & \displaystyle \frac{2{R}_{{\rm{NIP}}}}{{v}_{0}}+\displaystyle \frac{2{R}_{N}^{{\gamma }_{P}}\sin {\alpha }_{0}}{{v}_{0}\cos {\alpha }_{0}}(\alpha -{\alpha }_{0})\\ & & +\displaystyle \frac{(1+{\sin }^{2}{\alpha }_{0}){R}_{N}^{{\gamma }_{P}}}{{v}_{0}{\cos }^{2}{\alpha }_{0}}{(\alpha -{\alpha }_{0})}^{2}+O(3).\end{eqnarray} \tag{ 21 }$$

Finally, the preliminary form of the inversion result is obtained by substituting $(\alpha -{\alpha }_{0})$ using equation (18) in the traveltime expansion (21)

$$\begin{eqnarray}{t({{\bf{x}}}_{m}-{{\bf{x}}}_{0},{\bf{h}})| }_{{\bf{h}}=0} & \approx & \displaystyle \frac{2{R}_{{\rm{NIP}}}}{{v}_{0}}+\displaystyle \frac{2{\widehat{{\bf{n}}}}_{0}^{{\mathsf{T}}}}{{v}_{0}}({{\bf{x}}}_{m}-{{\bf{x}}}_{0})\\ & & +\displaystyle \frac{{({{\bf{x}}}_{m}-{{\bf{x}}}_{0})}^{{\mathsf{T}}}{\bf{M}}({{\bf{x}}}_{m}-{{\bf{x}}}_{0})}{{v}_{0}{R}_{N}^{{\gamma }_{P}}}+O(3),\end{eqnarray} \tag{ 22 }$$

where

$$\begin{eqnarray}{\bf{M}}=\left[\begin{array}{cc}1-{\cos }^{2}{\beta }_{0}{\sin }^{2}{\alpha }_{0} & -{\sin }^{2}{\alpha }_{0}\cos {\beta }_{0}\sin {\beta }_{0}\\ -{\sin }^{2}{\alpha }_{0}\cos {\beta }_{0}\sin {\beta }_{0} & 1-{\sin }^{2}{\beta }_{0}{\sin }^{2}{\alpha }_{0}\end{array}\right].\end{eqnarray} \tag{ 23 }$$

${\widehat{{\bf{n}}}}_{0}$ is the projected normal reference vector defined by (10).

It is worth observing that the preliminary traveltime expression (22) depends, via the curvature ${\kappa }_{N}^{{\gamma }_{P}}=1/{R}_{N}^{{\gamma }_{P}}$, on the relative position of P with respect to P₀

$$\begin{eqnarray}&&{\kappa }_{N}^{{\gamma }_{P}}=\displaystyle \frac{1}{{R}_{R}^{{\gamma }_{P}}+{R}_{{\rm{NIP}}}},\end{eqnarray} \tag{ 24 }$$

as illustrated in figure 4. Considering all points P close to P₀ on the quadratic region of the reflecting surface, expression (24) describes a fictitious wavefront, or N-wave, reaching the acquisition plane at ${{\bf{x}}}_{0}$ and emitted at the center of curvature C (exploding reflector).

To obtain the ultimate form that does not depend explicitly on P, matrix (23) may be written using its diagonal form ${\bf{D}}$, so that ${\bf{M}}=({{\bf{VD}}}^{1/2})\ ({{\bf{D}}}^{1/2}{{\bf{V}}}^{{\mathsf{T}}})$

$$\begin{eqnarray}{\bf{M}} & = & \left[\begin{array}{cc}\cos {\beta }_{0} & -\sin {\beta }_{0}\\ \sin {\beta }_{0} & \cos {\beta }_{0}\end{array}\right]\left[\begin{array}{cc}\cos {\alpha }_{0} & 0\\ 0 & 1\end{array}\right]\\ & & \times \,\left[\begin{array}{cc}\cos {\alpha }_{0} & 0\\ 0 & 1\end{array}\right]\left[\begin{array}{cc}\cos {\beta }_{0} & \sin {\beta }_{0}\\ -\sin {\beta }_{0} & \cos {\beta }_{0}\end{array}\right].\end{eqnarray} \tag{ 25 }$$

With the help of (25), the traveltime (22) can now be written as follows

$$\begin{eqnarray}t({{\bf{x}}}_{m}-{{\bf{x}}}_{0},0) & \approx & \displaystyle \frac{2{R}_{{\rm{NIP}}}}{{v}_{0}}+\displaystyle \frac{2{\widehat{{\bf{n}}}}_{0}^{{\mathsf{T}}}}{{v}_{0}}({{\bf{x}}}_{m}-{{\bf{x}}}_{0})\\ & & +\displaystyle \frac{{\kappa }_{N}^{{\gamma }_{P}}}{{v}_{0}}\parallel {{\bf{D}}}^{1/2}{{\bf{V}}}^{{\mathsf{T}}}({{\bf{x}}}_{m}-{{\bf{x}}}_{0}){\parallel }^{2}+O(3),\end{eqnarray} \tag{ 26 }$$

where the matrix product ${{\bf{D}}}^{1/2}{{\bf{V}}}^{{\mathsf{T}}}$ implements the orthogonal projection on the acquisition plane of the 3D rotation of angle $(-{\beta }_{0})$ leaving axis $\overrightarrow{{oz}}$ fixed, see figure 2, followed by the 3D rotation of angle ${\alpha }_{0}$ leaving axis $\overrightarrow{{oy}}$ fixed (left-hand convention).

Remark. Since the composition of the two 3D rotations aligns ${{\bf{n}}}_{0}$ to $\overrightarrow{{oz}}$, the rotated N-wavefront near ${{\bf{x}}}_{0}$ can be represented as $z={({\bf{x}}-{{\bf{x}}}_{0})}^{{\mathsf{T}}}{{\bf{K}}}_{N}({\bf{x}}-{{\bf{x}}}_{0})$, where ${{\bf{K}}}_{N}$ denotes the curvature matrix and ${\bf{x}}$ a point of the datum plane.

Recalling $\kappa ={\widehat{{\bf{e}}}}^{{\mathsf{T}}}{{\bf{K}}}_{N}\widehat{{\bf{e}}}$ is the curvature value that results from the intersection of the quadratic form and a vertical plane along a unit direction of the acquisition surface, $\parallel \widehat{{\bf{e}}}\parallel =1$, then, as a by-product of the projected rotations, one finds:

$$\begin{eqnarray}&&{\kappa }_{N}^{{\gamma }_{P}}=\displaystyle \frac{{[{{\bf{D}}}^{1/2}{{\bf{V}}}^{{\mathsf{T}}}({{\bf{x}}}_{m}-{{\bf{x}}}_{0})]}^{{\mathsf{T}}}{{\bf{K}}}_{N}[{{\bf{D}}}^{1/2}{{\bf{V}}}^{{\mathsf{T}}}({{\bf{x}}}_{m}-{{\bf{x}}}_{0})]}{\parallel {{\bf{D}}}^{1/2}{{\bf{V}}}^{{\mathsf{T}}}({{\bf{x}}}_{m}-{{\bf{x}}}_{0}){\parallel }^{2}}.\end{eqnarray} \tag{ 27 }$$

It follows, substituting (27) into (26), that the quadratic term of the traveltime expression is ultimately defined, via ${{\bf{K}}}_{N}$, by the two principal curvatures at ${{\bf{x}}}_{0}$ of the N-wave:

$$\begin{eqnarray}t({{\bf{x}}}_{m}-{{\bf{x}}}_{0},0) & \approx & \displaystyle \frac{2{R}_{{\rm{NIP}}}}{{v}_{0}}+\displaystyle \frac{2{\widehat{{\bf{n}}}}_{0}^{{\mathsf{T}}}}{{v}_{0}}({{\bf{x}}}_{m}-{{\bf{x}}}_{0})\\ & & +\displaystyle \frac{1}{{v}_{0}}{({{\bf{x}}}_{m}-{{\bf{x}}}_{0})}^{{\mathsf{T}}}({\bf{V}}\,{{\bf{D}}}^{1/2}){{\bf{K}}}_{N}\\ & & \times ({{\bf{D}}}^{1/2}{{\bf{V}}}^{{\mathsf{T}}})({{\bf{x}}}_{m}-{{\bf{x}}}_{0})+O(3).\end{eqnarray} \tag{ 28 }$$

Equation (28) is then the traveltime paraxial approximation for the zero-offset acquisition problem in a uniform medium, predicting the time of flight of a signal traveling close to the reference NIP-ray, emitted and received at ${{\bf{x}}}_{m}$.

3.2. Acquisition with the midpoint on the exit point of the reference ray

To analyse the situation with the acquisition setup where ${{\bf{x}}}_{m}({\bf{n}},h)={{\bf{x}}}_{0}$, one starts from (14):

$$\begin{eqnarray}&&\left\{\sqrt{{\left[\displaystyle \frac{h}{{\rm{\Gamma }}({\bf{n}})}\right]}^{2}+1}-1\right\}{\bf{u}}-\displaystyle \frac{{R}_{N}^{{\gamma }_{P}}}{{\rm{\Gamma }}({\bf{n}})}\left\{-\sin {\alpha }_{0}\left[\begin{array}{l}\cos {\beta }_{0}\\ \sin {\beta }_{0}\end{array}\right]\right.\\ &&\quad +\left.\tan \alpha \cos {\alpha }_{0}\left[\begin{array}{l}\cos \beta \\ \sin \beta \end{array}\right]\right\}=0,\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray*}{\rm{\Gamma }}({\bf{n}}) & = & \displaystyle \frac{{R}_{{\rm{NIP}}}({\bf{n}})}{2\sin \alpha ({u}_{1}\cos \beta +{u}_{2}\sin \beta )},\hspace{0.25cm}{R}_{{\rm{NIP}}}({\bf{n}})\\ & \approx & {R}_{{\rm{NIP}}}+{R}_{N}^{{\gamma }_{P}}\left(\displaystyle \frac{\cos {\alpha }_{0}}{\cos \alpha }-1\right),\hspace{0.25cm}\ {\bf{h}}=h{\bf{u}}.\end{eqnarray*}$$

• Objective: calculate the power series of $\alpha =\alpha (0,{\bf{h}})$ and $\beta =\beta (0,{\bf{h}})$ by linearizing (29) in the variables $\alpha -{\alpha }_{0}$ and $\beta -{\beta }_{0}$ and solve the resulting linear system:

  $$\begin{eqnarray*}&&\left[\begin{array}{c}{a}_{\mathrm{0,0}}\\ {b}_{\mathrm{0,0}}\end{array}\right]+\left[\begin{array}{cc}{a}_{\mathrm{1,0}} & {a}_{\mathrm{0,1}}\\ {b}_{\mathrm{1,0}} & {b}_{\mathrm{0,1}}\end{array}\right]\left[\begin{array}{c}\alpha -{\alpha }_{0}\\ \beta -{\beta }_{0}\end{array}\right]=0.\end{eqnarray*}$$

Then the two solutions must be expanded in Taylor series of the half offset ${\bf{h}}$ up to the second order. Using intensively Maple, the paraxial change of elevation and azimuth for small variations of ${\bf{h}}$ around ${{\bf{x}}}_{0}$ is a new result that takes the form

$$\begin{eqnarray}&&\alpha -{\alpha }_{0}=\displaystyle \frac{\cos {\alpha }_{0}\sin {\alpha }_{0}}{{R}_{N}^{{\gamma }_{P}}{R}_{{\rm{NIP}}}}{({h}_{1}\cos {\beta }_{0}+{h}_{2}\sin {\beta }_{0})}^{2}+O(3),\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}\beta -{\beta }_{0} & = & -\displaystyle \frac{1}{{R}_{N}^{{\gamma }_{P}}{R}_{{\rm{NIP}}}}({h}_{1}\cos {\beta }_{0}+{h}_{2}\sin {\beta }_{0})\\ & & \times (-{h}_{1}\sin {\beta }_{0}+{h}_{2}\cos {\beta }_{0})+O(3).\end{eqnarray} \tag{ 31 }$$

Expanding in $\alpha -{\alpha }_{0}$ and $\beta -{\beta }_{0}$ the squared traveltime expression (2), one obtains

$$\begin{eqnarray}&&{t}^{2}({\bf{n}},{\bf{h}})={\left(\displaystyle \frac{2h}{{v}_{0}}\right)}^{2}+\displaystyle \frac{2{R}_{{\rm{NIP}}}^{2}({{\bf{n}}}_{0})}{{v}_{0}^{2}}\left[\sqrt{\displaystyle \frac{{h}^{2}}{{{\rm{\Gamma }}}^{2}({{\bf{n}}}_{0})}+1}+1\right]\\ &&\quad +\left\{\displaystyle \frac{4{R}_{{\rm{NIP}}}{R}_{N}}{{v}_{0}}\tan {\alpha }_{0}\left[\sqrt{\displaystyle \frac{{h}^{2}}{{{\rm{\Gamma }}}^{2}({{\bf{n}}}_{0})}+1}+1\right]\right.\\ &&\quad +\displaystyle \frac{8\sin {\alpha }_{0}}{\sqrt{\tfrac{{h}^{2}}{{{\rm{\Gamma }}}^{2}({{\bf{n}}}_{0})}+1}+1}\left(\cos {\alpha }_{0}+\displaystyle \frac{{R}_{N}}{{R}_{{\rm{NIP}}}}\sin {\alpha }_{0}\tan {\alpha }_{0}\right)\\ &&\quad \times {({\widehat{{\bf{n}}}}_{0}^{{\mathsf{T}}}{\bf{h}})}^{2}\}\displaystyle \frac{(\alpha -{\alpha }_{0})}{{v}_{0}^{2}}+\displaystyle \frac{8{\sin }^{2}{\alpha }_{0}}{\sqrt{\tfrac{{h}^{2}}{{{\rm{\Gamma }}}^{2}({{\bf{n}}}_{0})}+1}+1}\\ &&\quad \times ({\widehat{{\bf{n}}}}_{0}^{{\mathsf{T}}}\,{\bf{h}})\left(\displaystyle \frac{{\rm{d}}{\widehat{{\bf{n}}}}_{0}^{{\mathsf{T}}}}{{\rm{d}}{\beta }_{0}}\,{\bf{h}}\right)\displaystyle \frac{(\beta -{\beta }_{0})}{{v}_{0}^{2}}+O(2).\end{eqnarray} \tag{ 32 }$$

After the substitution of (30) and (31) into the traveltime expression (32), taking the square root and expanding in power series of ${\bf{h}}$ up to the fourth order one obtains:

$$\begin{eqnarray}{t({{\bf{x}}}_{m}-{{\bf{x}}}_{0},{\bf{h}})| }_{{{\bf{x}}}_{m}={{\bf{x}}}_{0}} & \approx & \displaystyle \frac{2{R}_{{\rm{NIP}}}}{{v}_{0}}+\displaystyle \frac{{\kappa }_{{\rm{NIP}}}}{{v}_{0}}\parallel {{\bf{D}}}^{1/2}{{\bf{V}}}^{{\mathsf{T}}}{\bf{h}}{\parallel }^{2}\\ & & +O(4),\hspace{0.4cm}{\kappa }_{{\rm{NIP}}}=\displaystyle \frac{1}{{R}_{{\rm{NIP}}}}.\end{eqnarray} \tag{ 33 }$$

The matrix product ${{\bf{D}}}^{1/2}{{\bf{V}}}^{{\mathsf{T}}}$ is the same composite 3D rotation, followed by the orthogonal projection onto the acquisition plane, defined by equation (25).

Remark. Traveltime (33) takes a very simple form independent of ${R}_{N}^{{\gamma }_{P}}$.

In (33), ${\kappa }_{{\rm{NIP}}}$ may be interpreted as the curvature at ${{\bf{x}}}_{0}$ of the spherical diffraction wavefront, or NIP-wave, initiated at P₀ and traveling at a speed v₀. In the following, an analogue picture of the NIP-wave will be useful to construct the homeomorphic mapping.

3.3. Traveltime approximation for isotropic homogeneous media

Assembling equations (28) and (33), the desired paraxial approximation of the traveltime appears in the form of the power series (12), which is valid, up to the second order, for an arbitrary acquisition setup with a short offset h and a midpoint ${{\bf{x}}}_{m}$ close to ${{\bf{x}}}_{0}$, the exit point of the reference NIP-ray:

$$\begin{eqnarray}&&t({{\bf{x}}}_{m}-{{\bf{x}}}_{0},{\bf{h}})\approx \displaystyle \frac{2{R}_{{\rm{NIP}}}}{{v}_{0}}+\displaystyle \frac{2{\widehat{{\bf{n}}}}_{0}^{{\mathsf{T}}}}{{v}_{0}}({{\bf{x}}}_{m}-{{\bf{x}}}_{0})+\displaystyle \frac{1}{{v}_{0}}{({{\bf{x}}}_{m}-{{\bf{x}}}_{0})}^{{\mathsf{T}}}\\ &&\quad \quad \quad \quad \quad \quad \times \,({\bf{V}}\,{{\bf{D}}}^{1/2}){{\bf{K}}}_{N}({{\bf{D}}}^{1/2}{{\bf{V}}}^{{\mathsf{T}}})({{\bf{x}}}_{m}-{{\bf{x}}}_{0})\\ &&\quad \quad \quad \quad \quad \quad +\,\displaystyle \frac{{\kappa }_{{\rm{NIP}}}}{{v}_{0}}\parallel {{\bf{D}}}^{1/2}{{\bf{V}}}^{{\mathsf{T}}}{\bf{h}}{\parallel }^{2}+O(3).\end{eqnarray} \tag{ 34 }$$

Although this result is known, the way to obtain it, based only on Snell's law and 3D ray geometry, is new. Traveltime formula (34) must be completed with (7) and (11), the equations that provide, in the neighbourhood of the reference NIP trajectory, the paraxial description of the NIP-ray associated to the exit vector ${\bf{n}}$, figure 4.

Let $t({{\bf{x}}}_{m}-{{\bf{x}}}_{0},{\bf{h}})$ be the paraxial source–receiver traveltime of a signal traveling in an layered isotropic medium with a known constant near-surface velocity v₀. For the same acquisition geometry, let ${t}_{A}({{\bf{x}}}_{m}-{{\bf{x}}}_{0},{\bf{h}})$ be the paraxial time of flight in the auxiliary homogeneous medium with $v={v}_{0}$, where equation (34) is valid. By analogy with the 2D case [16], it is assumed in 3D that, near the surface, the two media are homeomorphic, or quasi-similar, in the sense that—in the auxiliary medium—the emerging NIP- and N-wavefronts are quadratic approximations of the real fronts.

The design of a synchronization mechanism between $t({{\bf{x}}}_{m}-{{\bf{x}}}_{0},{\bf{h}})$ and ${t}_{A}({{\bf{x}}}_{m}-{{\bf{x}}}_{0},{\bf{h}})$ is the important added value of this work regarding CRS theory in 3D. The time delay must be thought as caused by the different action of the two media on the propagation of the emerging wavefronts.

Remark. Since in the real medium the transform ${\bf{h}}\to -{\bf{h}}$ preserves the trajectory of the source–receiver ray, the power-series expansion of the traveltime takes, up to the second order, a form similar to (12):

$$\begin{eqnarray}t({{\bf{x}}}_{m}-{{\bf{x}}}_{0},{\bf{h}}) & = & {t}_{0}+{{\bf{a}}}^{{\mathsf{T}}}({{\bf{x}}}_{m}-{{\bf{x}}}_{0})+{({{\bf{x}}}_{m}-{{\bf{x}}}_{0})}^{{\mathsf{T}}}\\ & & \times {{\bf{K}}}_{m}({{\bf{x}}}_{m}-{{\bf{x}}}_{0})+{{\bf{h}}}^{{\mathsf{T}}}{{\bf{K}}}_{h}{\bf{h}}+O(3).\end{eqnarray} \tag{ 35 }$$

Here t₀ denotes the two-way traveltime of the reference ray along a curved normal-incidence trajectory, while ${{\bf{K}}}_{m}$ and ${{\bf{K}}}_{h}$ are symmetric $(2\times 2)$-matrices that can be calculated by alternatively setting ${\bf{h}}=0$ and ${{\bf{x}}}_{m}={{\bf{x}}}_{0}$.

4.1. Acquisition with the midpoint on the exit point of the reference ray

Equation (35) provides the inequality that is essential for the definition of the auxiliary medium. Actually, setting ${{\bf{x}}}_{m}={{\bf{x}}}_{0}$, figure 5, one finds the following upper bound for the traveltime:

$$\begin{eqnarray}&&| t(0,{\bf{h}})-{t}_{0}| \leqslant \parallel {\bf{h}}{\parallel }^{2}\,\rho ({{\bf{K}}}_{h})+O(4),\end{eqnarray} \tag{ 36 }$$

where the induced matrix norm $\rho ({{\bf{K}}}_{h})\geqslant 0$ denotes the spectral radius of matrix ${{\bf{K}}}_{h}$. In a homogeneous medium, comparing (35) and (34), it easy to verify that $\rho ({{\bf{K}}}_{h})=| {\kappa }_{{\rm{NIP}}}| /{v}_{0}$. The upper bound of (36) decreases either for short offsets $\parallel {\bf{h}}\parallel$, deep reflectors (via ${\kappa }_{{\rm{N}}{\rm{I}}{\rm{P}}}$) or large velocities v₀. Tightening these conditions, the reflection point P of each paraxial ray, see figure 5, moves within a small neighbourhood of the reference NIP-point P₀, covering a surface that becomes as large as the first Fresnel zone of the normal incident wave. When these conditions are met, each point P within the Fresnel zone can be viewed from the surface as if it were P₀, namely P ≈ P₀.

Zoom In Zoom Out Reset image size

This is precisely the hypothesis that is adopted to define the auxiliary medium. To proceed, one considers the diffraction wavefront initiated at P₀, or NIP-wave, where three rays are essential; they emerge at S, G and ${{\bf{x}}}_{0}$ along ${{\bf{n}}}_{S}$, ${{\bf{n}}}_{G}$ and ${{\bf{n}}}_{0}$, the three directions tangent to their trajectories. Beneath to the datum plane, where the near-surface velocity is assumed constant, $v={v}_{0}$, the three trajectories become straight lines.

Each diffraction ray hits the datum plane at three different instants: t_S, t_G and ${t}_{0}/2$. Then, under the hypothesis P ≈ P₀, the traveltime of the reflected paraxial ray, from S to G, can be approximated as follows:

$$\begin{eqnarray}&&t(0,{\bf{h}})\approx {t}_{S}+{t}_{G}.\end{eqnarray} \tag{ 37 }$$

Let ${\gamma }_{{\bf{h}}}$ denote the curve lying on the intersection of the NIP-wavefront and the plane containing ${\bf{n}}$ and ${\bf{h}}=h{\bf{u}}$ and let ${R}_{{\rm{NIP}}}^{{\gamma }_{{\bf{h}}}}$ be the radius of curvature of ${\gamma }_{{\bf{h}}}$ at ${{\bf{x}}}_{0}$ and ${P}_{A}^{{\gamma }_{{\bf{h}}}}$ the center of curvature. Figure 6 illustrates these definitions.

Zoom In Zoom Out Reset image size

Consequently, for a fixed acquisition setup, ${P}_{A}^{{\gamma }_{{\bf{h}}}}$ is the image point of P₀ in the constant-velocity auxiliary medium where $v={v}_{0}$. Moving P₀ and keeping constant the offset direction ${\bf{u}}$, the buried reflector is mapped onto the auxiliary medium.

Considering the quadratic approximation of the NIP-wavefront along the NIP-trajectory, ${\gamma }_{{\bf{h}}}$ becomes locally an arc of circle, centred at ${P}_{A}^{{\gamma }_{{\bf{h}}}}$, with a radius ${R}_{{\rm{NIP}}}^{{\gamma }_{{\bf{h}}}}$ that expands near the surface at a speed $v={v}_{0}$. In this approximation of ${\gamma }_{{\bf{h}}}$, the three rays may be treated as if their exit directions ${{\bf{n}}}_{S}$, ${{\bf{n}}}_{G}$ and ${{\bf{n}}}_{0}$ were coplanar and pointing ${P}_{A}^{{\gamma }_{{\bf{h}}}}$, a reasonable hypothesis in media with a low impedance contrast.

Accordingly, the two wavefronts, one starting in the auxiliary medium at ${P}_{A}^{{\gamma }_{{\bf{u}}}}$ while the other at P₀, are locally indistinguishable in the constant-velocity layer along the arc of circle ${\gamma }_{{\bf{h}}}$. The arc ${\gamma }_{{\bf{h}}}$ carried by the wavefront reaches point ${{\bf{x}}}_{0}$ and the source and the receiver positions, as illustrated in figure 7, at three different instants: ${t}_{0}^{\prime }/2$ , ${t}_{S}^{\prime }={t}_{0}^{\prime }/2+{{\rm{\Delta }}}_{S}$ and ${t}_{G}^{\prime }={t}_{0}^{\prime }/2+{{\rm{\Delta }}}_{G}$, where ${t}_{0}^{\prime }$ is the two-way traveltime of the NIP-ray. Although the traveltimes ${t}_{0}^{\prime }/2$, ${t}_{S}^{\prime }$ and ${t}_{G}^{\prime }$ are delayed in one medium with respect to the other, since the time intervals ${{\rm{\Delta }}}_{S}$ and ${{\rm{\Delta }}}_{G}$ are measured in the constant-velocity layer, the quantity ${{\rm{\Delta }}}_{G}+{{\rm{\Delta }}}_{S}={t}_{G}^{\prime }+{t}_{S}^{\prime }-{t}_{0}^{\prime }$ remains the same in both media.

Zoom In Zoom Out Reset image size

By assigning in each medium the proper value to ${t}_{{SG}}={t}_{G}^{\prime }+{t}_{S}^{\prime }$ and ${t}_{0}^{\prime }$, with the help of (37) one obtains the temporal synchronization of the two NIP-wavefronts:

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{G}+{{\rm{\Delta }}}_{S}={t}_{A}^{{\gamma }_{{\bf{h}}}}(0,{\bf{h}})-\displaystyle \frac{2{R}_{{\rm{NIP}}}^{{\gamma }_{{\bf{h}}}}}{{v}_{0}}\approx t(0,{\bf{h}})-{t}_{0}.\end{eqnarray} \tag{ 38 }$$

In the auxiliary medium, the traveltime ${t}_{A}^{{\gamma }_{{\bf{h}}}}(0,{\bf{h}})$ is given by (33) with ${\kappa }_{{\rm{NIP}}}^{{\gamma }_{{\bf{h}}}}$ equal to:

$$\begin{eqnarray}&&{\kappa }_{{\rm{NIP}}}^{{\gamma }_{{\bf{h}}}}=\displaystyle \frac{{[{{\bf{D}}}^{1/2}{{\bf{V}}}^{{\mathsf{T}}}{\bf{h}}]}^{{\mathsf{T}}}{{\bf{K}}}_{{\rm{NIP}}}[{{\bf{D}}}^{1/2}{{\bf{V}}}^{{\mathsf{T}}}{\bf{h}}]}{\parallel {{\bf{D}}}^{1/2}{{\bf{V}}}^{{\mathsf{T}}}{\bf{h}}{\parallel }^{2}}.\end{eqnarray} \tag{ 39 }$$

${{\bf{K}}}_{{\rm{NIP}}}$ is the symmetric $(2\times 2)$-matrix defined by the two principal curvatures at ${{\bf{x}}}_{0}$ of the emerging NIP-wavefront.

Finally, after the substitution of (33) and (39) into (38), the approximated traveltime of the reflected paraxial ray, for ${{\bf{x}}}_{m}={{\bf{x}}}_{0}$, becomes in the real medium:

$$\begin{eqnarray}&&t(0,{\bf{h}})\approx {t}_{0}+\displaystyle \frac{1}{{v}_{0}}{[{{\bf{D}}}^{1/2}{{\bf{V}}}^{{\mathsf{T}}}{\bf{h}}]}^{{\mathsf{T}}}{{\bf{K}}}_{{\rm{NIP}}}[{{\bf{D}}}^{1/2}{{\bf{V}}}^{{\mathsf{T}}}{\bf{h}}]+O(4).\end{eqnarray} \tag{ 40 }$$

The case of a zero-offset acquisition and the study of N-wave ray trajectories will provide the complete traveltime formula.

4.2. Zero-offset acquisition

When ${\bf{h}}=0$, source and receiver have the same coordinate ${{\bf{x}}}_{m}$, so that, moving ${{\bf{x}}}_{m}$ around ${{\bf{x}}}_{0}$, only seismic events that travel along normal rays in the vicinity of the normal reference ray are detected. Assuming that all these events are simultaneously initiated on the reflector surface close to the reference NIP-point (exploding reflector model), the resulting upgoing front is just an N-wave.

By considering the N-wave initiated in the neighborhood of P₀ where the quadratic approximation of the reflector surface is valid, then P₀ and a point P are connected by an arc of circle centred at C along ${{\bf{x}}}_{0}$. In figure 8, two normal rays are shown; one, the central, leaves the reflector from P₀ and reaches ${{\bf{x}}}_{0}$ along the emerging direction ${{\bf{n}}}_{0}$, while the second, leaving P, reaches ${{\bf{x}}}_{m}$ along direction ${\bf{n}}$. Each ray hits the datum plane at two different instants: ${t}_{0}/2$, the central, and ${t}_{m}/2$, the other. Within the quadratic approximation, ${{\bf{n}}}_{0}$ and ${\bf{n}}$ are coplanar. Let γ be the curve lying on the intersection of the N-wavefront and the plane containing ${{\bf{n}}}_{0}$ and ${\bf{n}}$, and let ${R}_{N}^{\gamma }$ denote the radius of curvature of γ at ${{\bf{x}}}_{0}$ and ${C}_{A}^{\gamma }$ the center of curvature.

As illustrated in figure 8, ${C}_{A}^{\gamma }$ is the image point of C in the constant-velocity auxiliary medium with $v={v}_{0}$. Then, by moving ${{\bf{x}}}_{m}$, the reflector is mapped onto the auxiliary medium.

Zoom In Zoom Out Reset image size

Considering the quadratic approximation of the N-wavefront along the reference NIP-trajectory, γ becomes locally an arc of circle centered at ${C}_{A}^{\gamma }$ with a radius ${R}_{N}^{\gamma }$ that expands near the surface at a speed $v={v}_{0}$. In this approximation of γ, the two rays may be treated as if their exit directions ${{\bf{n}}}_{0}$ and ${\bf{n}}$ were coplanar and pointing ${C}_{A}^{\gamma }$.

In the near surface constant-velocity layer, the two wavefronts, one starting in the auxiliary medium at ${C}_{A}^{\gamma }$ while the other at C, are locally indistinguishable, hence homeomorphic, along the arc of circle γ.

The arc γ reaches ${{\bf{x}}}_{0}$ and ${{\bf{x}}}_{m}$, as illustrated in figure 8, at two different instants: ${t}_{0}^{\prime }/2$ and ${t}_{m}^{\prime }/2$, where ${t}_{0}^{\prime }$ and ${t}_{m}^{\prime }$ are the two-way traveltimes of the NIP-rays under examination. These traveltimes are, however, delayed in one medium with respect to the other. Since the time interval ${{\rm{\Delta }}}_{m}={t}_{m}^{\prime }/2-{t}_{0}^{\prime }/2$ is measured in the constant-velocity layer, this quantity is preserved in both media.

By assigning in each medium the proper value to ${t}_{0}^{\prime }$ and ${t}_{m}^{\prime }$, one obtains the temporal synchronization of the two N-wavefronts:

$$\begin{eqnarray}&&2{{\rm{\Delta }}}_{m}={t}_{A}^{\gamma }({{\bf{x}}}_{m}-{{\bf{x}}}_{0},0)-\displaystyle \frac{2{R}_{{\rm{NIP}}}^{\gamma }}{{v}_{0}}=t({{\bf{x}}}_{m}-{{\bf{x}}}_{0},0)-{t}_{0}.\end{eqnarray} \tag{ 41 }$$

In the auxiliary medium, the traveltime ${t}_{A}^{\gamma }({{\bf{x}}}_{m}-{{\bf{x}}}_{0},0)$ is given by equation (28), in which ${{\bf{K}}}_{N}$ is the symmetric $(2\times 2)$-matrix defined by the two principal curvatures at x₀ of the homeomorphic N-wavefront.

After substitution of equation (28) into (41), the approximated time of flight of the two-way normal-incidence ray traveling in the stratified medium, emitted and received at ${{\bf{x}}}_{m}$, takes the form:

$$\begin{eqnarray}t({{\bf{x}}}_{m}-{{\bf{x}}}_{0},0) & \approx & {t}_{0}+\displaystyle \frac{2{\widehat{{\bf{n}}}}_{0}^{{\mathsf{T}}}}{{v}_{0}}\,({{\bf{x}}}_{m}-{{\bf{x}}}_{0})\\ & & +\displaystyle \frac{1}{{v}_{0}}{({{\bf{x}}}_{m}-{{\bf{x}}}_{0})}^{{\mathsf{T}}}({{\bf{VD}}}^{1/2}){{\bf{K}}}_{N}\\ & & \times ({{\bf{D}}}^{1/2}{{\bf{V}}}^{{\mathsf{T}}})({{\bf{x}}}_{m}-{{\bf{x}}}_{0})+O(3).\end{eqnarray} \tag{ 42 }$$

This last expression completes, for h = 0, the correction of the total traveltime (34) valid in the homeomorphic auxiliary medium.

4.3. Approximation of the complete traveltime in a layered medium

Assembling (40) and (42), one obtains the desired paraxial formula of the traveltime. It is a power series calculated up to the second order by assuming an acquisition setup near the exit point ${{\bf{x}}}_{0}$ of the reference normal-incidence ray travelling in an isotropic layered media with reasonable lateral velocity contrasts:

$$\begin{eqnarray}&&t({{\bf{x}}}_{m}-{{\bf{x}}}_{0},{\bf{h}})\approx {t}_{0}+\displaystyle \frac{2{\widehat{{\bf{n}}}}_{0}^{{\mathsf{T}}}}{{v}_{0}}\,({{\bf{x}}}_{m}-{{\bf{x}}}_{0})+\displaystyle \frac{1}{{v}_{0}}{({{\bf{x}}}_{m}-{{\bf{x}}}_{0})}^{{\mathsf{T}}}\\ &&\quad \times ({\bf{V}}\,{{\bf{D}}}^{1/2}){{\bf{K}}}_{N}({{\bf{D}}}^{1/2}{{\bf{V}}}^{{\mathsf{T}}})({{\bf{x}}}_{m}-{{\bf{x}}}_{0})\\ &&\quad +\displaystyle \frac{1}{{v}_{0}}\,{{\bf{h}}}^{{\mathsf{T}}}({\bf{V}}\,{{\bf{D}}}^{1/2}){{\bf{K}}}_{{\rm{NIP}}}({{\bf{D}}}^{1/2}{{\bf{V}}}^{{\mathsf{T}}}){\bf{h}}+O(3).\end{eqnarray} \tag{ 43 }$$

This expression is called the parabolic CRS-traveltime formula. Squaring (43) and keeping the same approximation order, one obtains the hyperbolic CRS traveltime, the well-known expression originally formulated starting from Bortfeld's work [8, 29]:

$$\begin{eqnarray}&&{t}^{2}({{\bf{x}}}_{m}-{{\bf{x}}}_{0},{\bf{h}})\approx {\left[{t}_{0}+\displaystyle \frac{2{\widehat{{\bf{n}}}}_{0}^{{\mathsf{T}}}}{{v}_{0}}({{\bf{x}}}_{m}-{{\bf{x}}}_{0})\right]}^{2}\\ &&\quad +\displaystyle \frac{2{t}_{0}}{{v}_{0}}{({{\bf{x}}}_{m}-{{\bf{x}}}_{0})}^{{\mathsf{T}}}({\bf{V}}\,{{\bf{D}}}^{1/2}){{\bf{K}}}_{N}({{\bf{D}}}^{1/2}{{\bf{V}}}^{{\mathsf{T}}})({{\bf{x}}}_{m}-{{\bf{x}}}_{0})\\ &&\quad +\,\displaystyle \frac{2{t}_{0}}{{v}_{0}}\,{{\bf{h}}}^{{\mathsf{T}}}({\bf{V}}\,{{\bf{D}}}^{1/2}){{\bf{K}}}_{{\rm{NIP}}}({{\bf{D}}}^{1/2}{{\bf{V}}}^{{\mathsf{T}}}){\bf{h}}+O(3).\end{eqnarray} \tag{ 44 }$$

These results are well known but, as already said, the path for reaching them, which is based on Snell's law, 3D ray geometry and the use of a quasi-similar earth model, is not only new, but it also opens a wide field of applications in seismic processing that are no longer based on reflections but, more generally, on diffractions [5, 7].

With the numerical value of the eight CRS attributes, both equations (43) and (44) provide—in the midpoint-offset domain surrounding ${{\bf{x}}}_{0}$—the moveout of all source–receiver traveltimes with respect to t₀. For each coordinate $({{\bf{x}}}_{0},{t}_{0})$ of the zero-offset volume, the moveout formulae allow to identify, collect and stack all primary events reflected and recorded within the Fresnel aperture in the midpoint-offset domain, so as to form globally a time-domain image that emulates seismic traces based on the specular reflection of normal incident rays. With these data, the subsequent step is to form a poststack-migrated reconstruction of the subsurface [12, 22].

CRS data stacking makes use of a very large number of traces with arbitrary midpoint-offset locations in the neighbourhood of ${{\bf{x}}}_{0}$, giving rise to an image with a superior signal-to-noise ratio and a higher resolution when compared to the standard common-midpoint stacking [32], which instead uses only traces around ${{\bf{x}}}_{m}={{\bf{x}}}_{0}$.

A computational strategy is however needed to efficiently and accurately estimate from seismic data the eight fields of CRS attributes. The solution to this problem represents a further obliged step for the industrial use of a CRS tool, designed to assist the interpretation of seismic data.

A strategy for the estimate of the eight CRS attributes $\xi =\{{\alpha }_{0},{\beta }_{0},{{\bf{K}}}_{{\rm{NIP}}},{{\bf{K}}}_{N}\}$ can be devised by solving concurrently a nonlinear optimisation problem for each coordinate $({{\bf{x}}}_{0},{t}_{0})$ of the zero-offset volume [1, 11]. Assuming that the near surface velocity v₀ is known, the attributes are determined by means of the coherence analysis [18], which essentially consists in maximizing the following semblance coefficient

$$\begin{eqnarray}&&S(\xi | {{\bf{x}}}_{0},{t}_{0})=\displaystyle \frac{1}{M}\ \displaystyle \frac{{\displaystyle \sum }_{k=-N/2}^{N/2}{\left|{\displaystyle \sum }_{i=1}^{M}{a}_{i,\tau (i)+k}\right|}^{2}}{{\displaystyle \sum }_{k=-N/2}^{N/2}\ {\displaystyle \sum }_{i=1}^{M}{| {a}_{i,\tau (i)+k}| }^{2}},\\ &&\quad \tau (i)=t[{{\bf{x}}}_{m}(i)-{{\bf{x}}}_{0},{\bf{h}}(i)],\end{eqnarray} \tag{ 45 }$$

where ${a}_{i,j}$ is the jth amplitude sample of the ith trace belonging to the midpoint-offset aperture around ${{\bf{x}}}_{0}$. Coefficient (45) is evaluated for a window of width N centred at $\tau (i)$, the traveltime prediction measured in sample units, either implementing (43) or (44). It is clear that ${t}_{0}=t(0,0)$. The value for the window width is evaluated running some initial test with $5\leqslant N\leqslant 15$, while the number of traces M, which depends on the extension of the projected Fresnel zone [19], can easily reach several tens of thousands.

It is worth noting that $0\leqslant S(\xi | {{\bf{x}}}_{0},{t}_{0})\leqslant 1$, the upper bound $S(\xi | {{\bf{x}}}_{0},{t}_{0})=1$ being achieved only if there is a perfect agreement between the CRS traveltime and the actual arrival time, and if all the amplitudes ${a}_{i,j}$ have the same value [18]. With seismic data, this upper bound is never attained because of the signal attenuation for long offsets.

It is an usual practice to solve a totally equivalent problem, minimising

$$\begin{eqnarray}&&f(\xi | {{\bf{x}}}_{0},{t}_{0})=1-S(\xi | {{\bf{x}}}_{0},{t}_{0}).\end{eqnarray} \tag{ 46 }$$

Since for each new ξ the estimate of the semblance has a high computational cost, the search process, if not correctly designed, might be very time consuming because the optimisation is performed in an eight-dimension volume. A quick but inaccurate greedy approach consists in splitting the problem in a nested sequence of approximations where for each optimisation step the search is limited to a reduced number of attributes on specific data collections [22, 26]. As illustrated with three parameters in 2D [13, 20], a local optimisation method for the simultaneous estimate of the eight parameters would certainly overcome all these deficiencies.

Given that the derivatives of the semblance are not explicitly available, a gradient-based minimization scheme cannot be implemented. An efficient way to keep good convergence properties without computing the gradients is to adopt a conjugate-direction method, together with a reliable line search algorithm [6]. Conjugate-direction methods, of which conjugate gradient is a particular case, are born to speed up the convergence rate of steepest descent. In case of convex cost functions, the iterations converge quadratically, starting from any real initial guess. Under these conditions, the fundamental property of these methods is that the result of the nth iteration is exactly the minimizer of (46) over the set ${{\bf{r}}}_{i}$ of conjugate directions in the attribute space. More precisely, if n nontrivial vectors ${{\bf{r}}}_{i}$, $i=0,1,\ldots ,n-1$, are mutually conjugate, the exact minimum of $f(\xi )$ can be obtained by a sequence of n one-dimensional searches [27]. Starting at point ${\xi }_{0}$, the final result ${\xi }_{\min }={\xi }_{n}$ is extracted from

$$\begin{eqnarray}&&f({\xi }_{i+1})=\mathop{\min }\limits_{{\alpha }_{i}\geqslant 0}f({\xi }_{i}+{\alpha }_{i}{{\bf{r}}}_{i}),\hspace{0.5cm}i=0,1,2,\ldots ,n-1.\end{eqnarray} \tag{ 47 }$$

However, it can be proved that these methods are locally supported by the same good convergence properties even if, as in (45) and (46), the cost function $f(\xi )$ is multimodal.

Recalling that a conjugate set ${{\bf{r}}}_{i}$, $i=0,1,\ldots ,n-1$ can be generated iteratively starting from a family of n linearly independent vectors ${{\bf{q}}}_{i}$, it comes out that the construction of the conjugate directions and the minimizing process can be carried ahead in the same iteration loop. For each coordinate $({{\bf{x}}}_{0},{t}_{0})$ of the zero-offset volume, the minimization algorithm, based on Powell's search method [28], is shown in figure 9.

Zoom In Zoom Out Reset image size

The goal of the line search is to minimize $h(\alpha )\ =f(\xi +\alpha {\bf{q}})$ along the direction ${\bf{q}}$ varying the steplength ${\alpha }_{{LB}}\leqslant \alpha \leqslant {\alpha }_{{UB}}$. To determine a new point ${\boldsymbol{\xi }}+{\boldsymbol{\alpha }}{\bf{q}}$, an effective adaptive selection rule for the steplength must be carefully implemented. A good choice is the adoption of the strong Wolfe–Powell conditions [25], a robust and reliable two-sided rule coupling a necessary and a sufficient condition. The only pre-condition to be checked is that ${h}^{\prime }(0)\leqslant 0$, so that ${\bf{q}}$ is locally a descent direction.

Unfortunately, the search for a good solution minimizing a multimodal cost function has to deal with the possibility of a premature convergence to a local minimum. Since semblance (45) is characterized by a large number of relative extrema, no optimization method based on a descent algorithm will prevent the trapping into local minima without allowing escape movements along the opposite direction to the descent lines.

The introduction of an uphill movement along each conjugate line towards regions of lower elevation greatly increases the capability of the line search algorithm to seek for solutions of lower cost. The idea [6] is to steer the search using a modified Wolfe–Powell rule to frame admissible solutions also along the counter-descent direction. In particular, to determine whether for a steplength $\alpha \lt 0$ a new point $\xi +\alpha {\bf{q}}$ is significantly better than the current one ξ, the necessary condition of Wolfe–Powell scheme has been properly modified, keeping, however, inalterate the sufficient one.

So far, this approach assumes that only one seismic event contributes to a specific pixel location of the zero-offset volume. More precisely, the optimisation does not provide the CRS attributes ξ for each conflicting event, such as diffractions that intersect reflections at $({{\bf{x}}}_{0},{t}_{0})$. However, initialising for each pixel a set of ${\xi }_{0}$ computed by a fast global optimisation heuristic, the modified conjugate-direction algorithm can concurrently reach a set of very good solutions ξ, each one of them resolving a different conflicting event. Sorting by importance these events, the weighted summation of the resulting pixel's values provides a practical solution to form an image [23, 26].

Compared to depth migration and seismic inversion, CRS processing is a powerful ersatz to form images in time domain when input traces are plagued by a low signal-to-noise ratio and the seismic velocity model, which is the fundamental ingredient for imaging in depth, may be very hard or even impossible to design. In similar situations, CRS imaging becomes the battle horse used to extract and enhance weak signals reflected by thin structures in fragmented and smooth geological background sedimentations. Unfortunately, small faults are not adequately focused: this deficiency is overcome by the formulation of a theory for imaging diffractions [7].

Since a CRS application implements a very large collection of concurrent optimisation problems, imaging works as an autofocusing method that collects and adaptively maps, independently of any velocity model, seismic signals from the midpoint-offset domain into the zero-offset stacked volume before time migration.

The quality of the image is fundamental for the interpretation of the subsurface layering. This activity makes use of the 3D images to trace geological strata horizons and faults and, finally, to build a conceptual model of the reservoir that will be used to locate the position of new wells. Horizons and faults are interpreted following the continuity or the discontinuity of the signal, having in mind a picture that is based on different geological sources of information that must be refined in order to sketch the potential reservoir location. Since the sketch is based on horizons and faults, it is only a geometrical model. From seismic attributes, it is hard to extract information regarding petrophysical quantities such as lithology, porosity, saturation, density and so on.

The figures shown refers to 3D land data with a low signal-to-noise ratio, which in part depends on the bad coupling between receivers and soil, and in part is caused by velocity inhomogeneities present in the near subsurface. These problem are normally handled by sources and receivers static corrections [32], whose estimation is one the most difficult task in the seismic imaging workflow. The three images illustrate the benefits of CRS technology and refer to data in time migrated domain. For every comparison with the conventional technique, the residual statics correction—applied before stacking—and the final seismic velocity model are the same, so that the remarkable difference in the final results is only due to the use of CRS stacking.

Figure 11 deals with a zone near the surface characterised by a high frequency content and a very low signal-to-noise ratio impeding the building of a reliable stacking velocity field. In the conventional approach, this is the cause of the lack of horizon continuity. The strong difference between the two images is indeed due to the intrinsic data-driven nature of CRS theory and to the consequent huge number of traces involved in the stacking procedure. It is evident that the interpretation on the CRS volume is incomparable easier than on the conventional result: the striking improvement is, apart from the global quality of the image, the better continuity of the layering and the identification of two faults in the lower part of the image (yellow arrows).

In figure 10, which refers to a deeper region with respect to figure 11, several faults are present on the CRS output that are totally invisible on the conventional image. The presence of these faults completely changed the macro geological model of the area, allowing a better understanding of the reservoir location. Before migration, these events appear in the stacked image as diffractions but not easily distinguishable from reflections. Although CRS is a theory based on reflections, diffraction traveltimes may be approximated to the second order by setting ${{\bf{K}}}_{{\rm{NIP}}}={{\bf{K}}}_{N}$ in (43) and (44). With this condition, the reflector element shrinks into a diffractor since its curvature goes to infinity. Normally, diffraction signals are ten or twenty times weaker than reflections, so that it is very hard to recover them using the conventional processing. The eight-parameter search driven by data allows one to intercept diffractions only when the optimisation algorithm is able to properly match the value of ${{\bf{K}}}_{{\rm{NIP}}}$ and ${{\bf{K}}}_{N}$.

Figure 12 shows a different use of the output volumes. The pictures display a depth slice in the 3D volume obtained by applying a coherence filter, via the semblance function (45), on the post-stack migrated image. The result allows one to understand the continuity of the signal in 3D emphasizing the shape of the discontinuities of the horizontal slices as well as the presence of faults or of incoherent conglomerates. The three red lines represent the positions of existing wells. The problem was to obtain a better comprehension of the subsurface characterised by a big fault (north–south direction) which delimited the reservoir. The target was to identify the reservoir limits in order to plan some new well and then enhance the oil extraction. From the CRS volume it was possible to obtain a better definition of the fault and a more detailed structural model of the whole area including the detection of a second big fault clearly visible on the right of the CRS slice. As a final comment on the cases illustrated by the three figures, where the construction of a reliable velocity model is impossible, the impressive achievement reached by the CRS application was fully validated by expert geologists.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The presented study exploits the full potential of the homeomorphic transform by constructing for the first time in 3D the adequate correspondence between a real layered medium and an auxiliary one. In this idealisation, a homogeneous earth model is designed to emulate the seismic response of the layered medium as the quadric approximation of the true wavefront emerging around a reference point ${{\bf{x}}}_{0}$. The two fronts reach the surface with different times of flight but with the same exit direction and the same radii of curvature. These two fronts are then quasi-similar around ${{\bf{x}}}_{0}$ and travel at the same speed v₀ in the near surface. As a consequence, this conceptual framework yields the estimate of the time correction to synchronise each seismic event with the exact traveltime expression calculated in the auxiliary medium. The eight-parameter CRS formula arises then as the final product.

CRS theory is a paradigm and in practice is still a cornerstone of time processing in seismic exploration, but it has some intrinsic theoretical and methodological limits. One is related to the difficulty, for a theory exclusively based on reflections, of properly modelling the response of a diffracting point. Another is due to the second-order truncation of the traveltime formula in the midpoint-offset domain, which limits the ability of a CRS application to resolve details. Finally, CRS theory assumes that only one seismic event contributes to form a specific pixel location of the stack volume, causing the so-called conflicting dips.

To overcome these restrictions, the natural extension of this work is the formulation of a theory capable to deal not only with reflectors but also with small structures such as faults and cracks that generate diffraction waves. Actually, the use of an auxiliary medium, homeomorphic to a fine-grained 3D earth model composed of diffracting points, leads, without series expansion, to a double-square-root traveltime with only five geometric parameters. In 2D, where the theory requires only two parameters, some remarkable results are already published and show the feasibility of a time-migration scheme solely driven by data and velocity-independent [5, 7]. A prototypal study, not yet published, demonstrates the huge potential of the time-migration application ported in 3D.

Moreover, the advantage of a data-driven theory based on diffractions is the degree of freedom that is obtained by the independent control of the two branches, descending and ascending, of the double-square-root expression, either in isotropic or anisotropic media. This feature allows as well to model the time of flight of converted rays or to obtain an effective traveltime formula for stacking data to a fixed source-to-receiver offset by assigning a reference ray to each branch of the double square-root.

This project was funded by Eni Upstream and Technical Services. The author thanks Eni's management for the images and the permission to present this work. I am also grateful to Paolo Marchetti and Daniela Theis for useful conversations and suggestions that helped to improve this article.