\documentclass[11pt,reqno]{amsart}
\usepackage[T2A]{fontenc}
\usepackage[utf8]{inputenc}
% \usepackage[cp1251]{inputenc}
\usepackage[english,russian]{babel}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{graphicx}
\usepackage{xcolor}
\usepackage[colorlinks]{hyperref}

\newcounter{aequation}
\newcounter{bequation}
\newcounter{cequation}
\newenvironment{aequation}{
  \[
  \refstepcounter{aequation}
  \tag{\arabic{aequation}A}
}{\]}
\newenvironment{bequation}{
  \[
  \refstepcounter{bequation}
  \tag{\arabic{bequation}B}
}{\]}
\newenvironment{cequation}{
  \[
  \refstepcounter{cequation}
  \tag{\arabic{cequation}C}
}{\]}

\usepackage{titlesec}
\titlespacing{\section}{0pt}{2.5ex}{1.5ex}
\titlespacing{\subsection}{0pt}{1.5ex}{1ex}
\titlespacing{\subsubsection}{0pt}{1.5ex}{1ex}
\titleformat{\section}{\large\bfseries\centering}{\thesection}{1em}{}
\titleformat{\subsection}[runin]{\bfseries}{\thesubsection.}{0.5em}{}[.\mbox{\ }]
\titleformat{\subsubsection}[runin]{\bfseries}{\thesubsubsection.}{0.4em}{}[.\mbox{\ }]





%%%%%---------- ORCiD
\newcommand\orcidicon[1]{\href{https://orcid.org/#1}{\includegraphics[scale=0.02]{orcid.pdf}}}


\def\udcs{550.34.013.4} %Here the author places classificators of the paper according to Russian classification system
\def\mscs{86-08} %Here the author places  classificators according to the AMS classification list.
\setcounter{page}{144}

\newtheorem{lemma}{Lemma}
\newtheorem{theorem}{Theorem}
\newtheorem{definition}{Definition}
\newtheorem{corollary}{Corollary}
\newtheorem{proposition}{Proposition}



\hypersetup{linkcolor=blue,filecolor=black,urlcolor=blue, citecolor=blue}

\def\logo{\textcolor{blue}{{\bf\Huge S\raisebox{0.2ex}{\hspace{0.55ex}\raisebox{0.05ex}e\hspace{-1.65ex}$\bigcirc$}MR}}}

\def\semrtop
     {
  \vbox{
     \noindent\logo \hspace{-3mm}  \parbox{12cm}{     \begin{center}
     {\LARGE \textsf{СИБИРСКИЕ ЭЛЕКТРОННЫЕ}} \\[2mm]
     {\LARGE \textsf{МАТЕМАТИЧЕСКИЕ ИЗВЕСТИЯ}} \\[2mm]
     {\Large \textsf{Siberian Electronic Mathematical Reports}} \\[1mm]
     {\Large\tt{\textcolor{blue}{http://semr.math.nsc.ru}}}\\[1mm]
     \raisebox{1ex}{ISSN 1813-3304 }
     \vspace*{-3mm}
%     {\small Институт математики им. С.Л.Соболева СО РАН}\\[-1mm]
%     {\small Sobolev Institute of Mathematics SB RAS}
     \end{center}
}
     \vspace{-1mm}
     \noindent
      \begin{tabular}{c}
     \hphantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa!} \\
     \hline\hline
     \end{tabular}

     \vspace{3mm}
     {\flushleft\it Vol. ??, No. ?, pp ???--??? (2025) \hspace{\fill}{\rm\small УДК \udcs \hspace{-4mm}}}
     \newline  %Volume, pages in Russian, filled by Editorial board
     	{\rm\small \textcolor{blue}{https://doi.org/10.33048/semi.2025.20.???}}\hphantom{aaaaaaaaaaaaaaaaaaaaaaaa\;}{\rm\small MSC\ \ \mscs }
  }%\vbox
}




\begin{document}
\renewcommand{\refname}{References}
\renewcommand{\proofname}{Proof.}
\renewcommand{\figurename}{Fig.}




\thispagestyle{empty}

\title[CDR tomography by the DSR equation asymptotics]{\large Controlled Directional Reception tomography based on the ray method asymptotics of the Double Square Root equation}


\author[N.N. Shilov]{{\bf N.N. Shilov}\protect\orcidicon{0000-0002-7618-1229}}



\address{Nikolay Nikolaevich Shilov  % Контактные данные всех авторов и место работы указываются только на английском языке
\newline\hphantom{iii} Novosibirsk State University,
\newline\hphantom{iii} ul. Pirogova, 2,
\newline\hphantom{iii} 630090, Novosibirsk, Russia}%
\email{\textcolor{blue}{n.shilov@g.nsu.ru}}%

\thanks{\sc Shilov, N.N.,
CDR tomography by the DSR equation asymptotics}
\thanks{\copyright \ 2025 Shilov N.N.}
\thanks{\rm The work is supported by the Ministry of Science and Higher Education of the Russian Federation (grant FSUS-2025-0015)}
\thanks{\it Received ???, ?, 2025, Published ???, ??, 2025}%


\semrtop \vspace{1cm}
\maketitle {\small

\vspace{-12pt}
\centerline{{\it Communicated by} {\sc P.P. Petrov}
%({\it filled by Editorial board})
}

\bigskip
\bigskip

\begin{quote}
\noindent{\bf Abstract:} Controlled Directional Reception (CDR) is a reflection tomography technique that accepts seismic traveltimes and slopes (traveltime derivatives w.r.t source and receiver coordinates) and returns a velocity model fitting this data. In contrast to other slope-based methods, it uses parsimonious model parametrization and relies on ray tracing thus being computationally efficient and fairly general. However, it is unstable w.r.t. data errors. In this paper we revisit the CDR method and develop a formalism to mitigate its instability. Our approach is based on linearized estimates of ray tracing errors allowing for suboptimal regularization of the inverse problem. We apply our original ray method asymptotics of the pseudodifferential Double Square Root equation to parametrize the wavefield and test our formulation of the CDR method on two benchmark synthetic datasets. We demonstrate that it provides competitive results suitable for depth migration. We restrict ourselves to 2D settings although the approach can be generalized to 3D problems as well.\medskip

\noindent{\bf Keywords:} slope tomography, reflection tomography, Controlled Directional Reception, Double Square Root equation.
 \end{quote}
}

\bigskip

\section{Introduction}

Slope tomography is a distinct group of reflection tomography techniques. They accept a set of so-called events (vectors containing source and receiver coordinates, reflection traveltimes and slopes i.e. traveltime derivatives w.r.t. source and receiver coordinates) and adjust the velocity model to fit this data. The resulting model is suitable for depth migration or full waveform inversion \cite{1}. The earliest slope tomography methods date back to 1950s and were based on straight-rays considerations \cite{2}. In 1987 Charles Sword reformulated the original approach for heterogeneous media with curved rays, dubbing his formulation the Controlled Directional Reception (CDR) method \cite{3}. However, it lacked stability w.r.t. data errors, and in 1998 Frederic Billette et al. introduced more stable and general method of stereotomography \cite{4}, which became a classic reference for coming algorithms \cite{5, 6}. Nowadays slope tomography is widely adopted and developed for both 2D and 3D data, isotropic and anisotropic models \cite{7, 8, 9}. Notable industrial examples are Slb’s CIP tomography \cite{10} and Viridien’s slope tomography \cite{11}, and a number of academic studies have been carried out (see \cite{1, 12} for overview).

Each event in seismic data corresponds to a reflecting facet (a fragment of an interface specified by its coordinates, orientation, and the reflection angle) in the subsurface. In classic stereotomography one iteratively traces incident and reflected rays upward from trial facets and then updates both the velocity and the facets to minimize data mismatch \cite{1}. In contrast, the CDR method deals with the velocity model only. Instead of tracing rays upward from unknown reflectors, one traces them downward from fixed sources and receivers, passing the observed slopes as the initial conditions for ray tracing. The tracing is stopped as soon as the total traveltime along the incident and the reflected rays equals the observed one. If the rays do not meet each other as the tracing stops, the velocity model is updated until they converge. The reflecting facets are determined then automatically \cite{3}. On the one hand, this concept benefits from lower dimensionality of the inverse problem. On the other hand, it suffers from instability of one-point ray tracing \cite{1}. However, this drawback formally can be attributed to some more recent slope tomography algorithms \cite{5, 6} as well. Hence, the CDR method still can be a good starting point for new slope tomography algorithms.

Double Square Root (DSR) equation is a one-way wave equation describing singly scattered wavefield as a function of source and receiver coordinates. It was introduced in \cite{13} although its kinematics had been described earlier \cite{14}. It acts on the recorded wavefield and allows one to extrapolate the observed data to another depth level, performing so-called “survey-sinking” migration \cite{13, 15}. In this migration, one evaluates the extrapolated wavefield at all depths at zero time and extracts zero offset traces (i.e. with coinciding source and receiver) only to obtain the subsurface image. Significant energy at nonzero offsets provides a basis for velocity-building techniques known as migration velocity analysis (MVA) \cite{16, 17}. The “survey-sinking” migration and corresponding MVA much resemble the CDR method where one “sinks” the source and the receiver and updates the velocity model until the rays converge (i.e. the energy focuses at zero offsets). This similarity has been employed in \cite{18} to reformulate the Sword’s approach using a ray method asymptotics of the DSR equation \cite{19}. However, the authors limited themselves to noise-free data inversion. In this paper we consider more realistic cases and stress on the method’s regularization.

The paper is organized as follows. First, we briefly outline theoretical aspects of the proposed formalism and build up a set of tools for iterative velocity update following the CDR approach. Then we linearize the inverse problem in simplified settings to derive suboptimal estimates of regularization parameters which we use later in the original nonlinear problem. Finally, we provide some technical information on our implementation and proceed with numerical tests. In the Discussion and Conclusion section, we weigh our results versus alternative approaches and summarize the paper’s contents. We bring more technical features in three Appendices in the end of the paper.


\section{Method's theory}

\subsection{DSR equation rays}

Let us consider a smooth 2D velocity distribution $v(x,z)$ with z axis being the depth, i.e. pointed downward. Let $\tau(x_s,x_r,z)$ be the traveltime of a reflected wave emitted from a source $(x_s,z)^T$ and recorded at a receiver $(x_r,z)^T$. Superscript $T$ denotes transposition. Traveltime $\tau$ is described by an equation \cite{13, 14}
\begin{equation}\label{f1}
     \dfrac{\partial\tau}{\partial z}=-\sqrt{\dfrac{1}{v_s^2}-\left(\dfrac{\partial\tau}{\partial x_s}\right)^2}-\sqrt{\dfrac{1}{v_r^2}-\left(\dfrac{\partial\tau}{\partial x_r}\right)^2},
\end{equation}
where $v_s=v\left(x_s,z\right)$ and $v_r=v\left(x_r,z\right)$ are the velocities at the source and receiver points respectively. We assume that the expressions under square roots keep positive, which implies nowhere-horizontal wave propagation. We shall briefly touch this restriction in the Discussion.

By solving \eqref{f1} one can reconstruct the surface $\tau\left(x_s,x_r,z\right)=0$ which defines all possible source-receiver pairs where the reflected wave is recorded just as the source fires. We shall search for such pairs using ray method \cite{19, 20}. We define a phase space consisting of 6 coordinates:
\begin{equation}\label{f2}
    \vec{x}=\left(x_s,x_r,z\right)^T,\quad\vec{p}=\left(p_s,p_r,p_z\right)^T=\left(\dfrac{\partial\tau}{\partial x_s},\dfrac{\partial\tau}{\partial x_r},\dfrac{\partial\tau}{\partial z}\right)^T.
\end{equation}
Vector $\vec{x}\left(\tau\right)$ defines a source $\left(x_s,z\right)^T$ and a receiver $\left(x_r,z\right)^T$ in physical space. Thus, each point of the phase space describes a data configuration: a reflected wave will be recorded at the time $\tau$ if one places a source-receiver pair at $\vec{x}\left(\tau\right)$ while the traveltime derivatives will be described by the vector $\vec{p}\left(\tau\right)$.

To trace the rays we introduce a Hamiltonian \cite{19}:
\begin{equation}\label{f3}
    \begin{gathered}
        H\left(\vec{x},\vec{p}\right)=-\dfrac{p_z+\sqrt{\dfrac{1}{v_s^2}-p_s^2}+\sqrt{\dfrac{1}{v_r^2}-p_r^2}}{\dfrac{1}{v_s^2\sqrt{\dfrac{1}{v_s^2}-p_s^2}}+\dfrac{1}{v_r^2\sqrt{\dfrac{1}{v_r^2}-p_r^2}}}= \\ = -\left(p_z+\sqrt{\dfrac{1}{v_s^2}-p_s^2}+\sqrt{\dfrac{1}{v_r^2}-p_r^2}\right)C_\tau\left(\vec{x},\vec{p}\right).
    \end{gathered}
\end{equation}

The DSR equation ray tracing system reads
\begin{equation}\label{f4}
    \begin{cases}
        \dfrac{d\vec{x}}{d\tau}=\nabla_{\vec{p}}H, \\[6pt]
        \dfrac{d\vec{p}}{d\tau}=-\nabla_{\vec{x}}H,
    \end{cases}
\end{equation}
where $\nabla$ is the gradient operator and its subscript indicates the variables of differentiation. We impose the initial conditions for this system at the Earth’ upper surface $\hat{z}=0$:
\begin{equation}\label{f5}
    \begin{cases}
        \left.\vec{x}\right|_{\tau=\hat{\tau}}=\left({\hat{x}}_s,{\hat{x}}_r,\hat{z}=0\right)^T=\widehat{\vec{x}}, \\
        \left.\vec{p}\right|_{\tau=\hat{\tau}}=\left({\hat{p}}_s,{\hat{p}}_r,{\hat{p}}_z=-\sqrt{\dfrac{1}{{\hat{v}}_s^2}-{\hat{p}}_s^2}-\sqrt{\dfrac{1}{{\hat{v}}_r^2}-{\hat{p}}_r^2}\right)^T=\widehat{\vec{p}},
    \end{cases}
\end{equation}
with ${\hat{v}}_s=v\left({\hat{x}}_s,\hat{z}=0\right)$ and ${\hat{v}}_r=v\left({\hat{x}}_r,\hat{z}=0\right)$. The “capped” variables represent conventional slope tomography data including actual source and receiver coordinates ${\hat{x}}_s$ and ${\hat{x}}_r$, observed traveltime $\hat{\tau}$ and its derivatives ${\hat{p}}_s=\dfrac{\partial\hat{\tau}}{\partial{\hat{x}}_s}$ and ${\hat{p}}_r=\dfrac{\partial\hat{\tau}}{\partial{\hat{x}}_r}$.

System \eqref{f4} with conditions \eqref{f5} define a DSR ray, i.e. a trajectory $\left(\vec{x}\left(\tau\right),\vec{p}\left(\tau\right)\right)$. By tracing it backwards in time one can reconstruct a source-receiver pair $\vec{x}\left(\tau=0\right)$ where the reflected wave is recorded just as the source fires. We introduce the subsurface midpoint and half-offset to describe this pair:
\begin{equation}\label{f6}
    m=\left.\dfrac{x_r+x_s}{2}\right|_{\tau=0},\quad h=\left.\dfrac{x_r-x_s}{2}\right|_{\tau=0}.
\end{equation}
A physical wave cannot instantly get across a finite distance between the source and the receiver. Thus, assuming error-free data and accurate velocity model one can expect $h$ to be zero and $m$ to be the reflection point’s horizontal coordinate.

\subsection{DSR equation ray perturbations}

Let us assume that the velocity model is inaccurate so that $h\neq0$. We are going to build a framework for iterative velocity update aiming to reduce $h$ to zero. To do so, we parametrize the model by a parameters’ matrix $\mathbf{V}$:
\begin{equation}\label{f7}
    v=v\left(x,z;\mathbf{V}\right),\quad \mathbf{V}=\left(\begin{matrix}v_{11}&\cdots&v_{1N}\\\vdots&\ddots&\vdots\\v_{M1}&\cdots&v_{MN}\\\end{matrix}\right).
\end{equation}
We assume that $v\left(x,z;\mathbf{V}\right)$ is twice smooth w.r.t. spatial coordinates and linear w.r.t. $\mathbf{V}$. By differentiating \eqref{f4} and \eqref{f5} w.r.t. a parameter $v_{mn}$ we obtain:
\begin{equation}\label{f8}
    \begin{cases}
        \dfrac{d}{d\tau}\dfrac{\partial\vec{x}}{\partial v_{mn}}=\nabla_{\vec{p}}\nabla_{\vec{x}}H\dfrac{\partial\vec{x}}{\partial v_{mn}}+\nabla_{\vec{p}}\nabla_{\vec{p}}H\dfrac{\partial\vec{p}}{\partial v_{mn}}+\nabla_{\vec{p}}\dfrac{\partial H}{\partial v_{mn}}, \\[6pt]
        \dfrac{d}{d\tau}\dfrac{\partial\vec{p}}{\partial v_{mn}}=-\nabla_{\vec{x}}\nabla_{\vec{x}}H\dfrac{\partial\vec{x}}{\partial v_{mn}}-\nabla_{\vec{x}}\nabla_{\vec{p}}H\dfrac{\partial\vec{p}}{\partial v_{mn}}-\nabla_{\vec{x}}\dfrac{\partial H}{\partial v_{mn}},
    \end{cases}
\end{equation}
and
\begin{equation}\label{f9}
    \begin{cases}
        \left.\dfrac{\partial\vec{x}}{\partial v_{mn}}\right|_{\tau=\hat{\tau}}=\vec{0}, \\
        \left.\dfrac{\partial\vec{p}}{\partial v_{mn}}\right|_{\tau=\hat{\tau}}=\left(\begin{matrix}0\\0\\\dfrac{\dfrac{\partial{\hat{v}}_s}{\partial v_{mn}}}{{\hat{v}}_s^3\sqrt{\dfrac{1}{{\hat{v}}_s^2}-{\hat{p}}_s^2}}+\dfrac{\dfrac{\partial{\hat{v}}_r}{\partial v_{mn}}}{{\hat{v}}_r^3\sqrt{\dfrac{1}{{\hat{v}}_r^2}-{\hat{p}}_r^2}}\\\end{matrix}\right),
    \end{cases}
\end{equation}
where the symbol $\vec{0}$ stands for a vector of zeros. We refer the reader to \cite{18, 19} for the Hamiltonian derivatives from \eqref{f4} and \eqref{f8} and discuss the model parametrization in the Practical implementation section.

The solution of \eqref{f8} determines the ray’s partial derivatives w.r.t. velocity perturbations caused by the parameter $v_{mn}$. By solving it for all the elements of $\mathbf{V}$, one can obtain the gradient of the subsurface half-offset \eqref{f6} w.r.t. velocity parameters:
\begin{equation}\label{f10}
    \nabla_{\mathbf{V}} h=\left.\dfrac{\nabla_{\mathbf{V}}x_r-\nabla_{\mathbf{V}}x_s}{2}\right|_{\tau=0}=\left(\dfrac{\partial h}{\partial v_{11}},\dfrac{\partial h}{\partial v_{12}},\ldots,\dfrac{\partial h}{\partial v_{MN}}\right)^T
\end{equation}
and perform an iterative velocity model update aiming to minimize the absolute value of $h$.

\section{Inverse problem statement}

\subsection{Linearized CDR method for a single ray in homogeneous medium}

Let us restrict ourselves to the class of homogeneous velocity models parametrized by a single parameter $\hat{v}$:
\begin{equation}\label{f11}
    v\left(x,z;\mathbf{V}\right)\equiv\hat{v}=const.,\quad\mathbf{V}=\left(\hat{v}\right).
\end{equation}
In such media the subsurface half-offset \eqref{f6} reads:
\begin{equation}\label{f12}
    \hat{h}=\dfrac{{\hat{x}}_r-{\hat{x}}_s}{2}-\left(\dfrac{{\hat{p}}_r}{\sqrt{\dfrac{1}{{\hat{v}}^2}-{\hat{p}}_r^2}}-\dfrac{{\hat{p}}_s}{\sqrt{\dfrac{1}{{\hat{v}}^2}-{\hat{p}}_s^2}}\right)\dfrac{C_\tau\left(\widehat{\vec{x}},\widehat{\vec{p}}\right)\hat{\tau}}{2},
\end{equation}
with $C_\tau\left(\widehat{\vec{x}},\widehat{\vec{p}}\right)$ defined in \eqref{f3} and “capped” variables in \eqref{f5}. We introduced $\hat{h}$ to distinguish it from the general case and emphasize its dependence on $\hat{v}$.

The CDR method consists of zeroing out $\hat{h}$ by appropriate choice of $\hat{v}$. However, errors in the observed data will cause a perturbation of the subsurface half-offset $\delta\hat{h}$, making it nonzero in the precise velocity model. Thus, we propose regularizing the inverse problem to avoid overfitting. Let us search for a small perturbation $\delta\hat{v}$ minimizing a linearized and regularized loss functional:
\begin{equation}\label{f13}
    \hat{L}\left(\delta\hat{v}\right)=\left(\hat{h}+\delta\hat{h}+\dfrac{\partial\hat{h}}{\partial\hat{v}}\delta\hat{v}\right)^2+\hat{\alpha}\dfrac{{\delta\hat{v}}^2}{{\hat{v}}^2}\rightarrow\min_{\delta\hat{v}},
\end{equation}
where $\hat{\alpha}$ is the Tikhonov’s regularization parameter \cite{21}. Given $\delta\hat{h}$, its optimal value can be found analytically:
\begin{equation}\label{f14}
    \hat{\alpha}=\sup{\left|\delta\hat{h}\right|}\hat{v}\dfrac{\partial\hat{h}}{\partial\hat{v}}.
\end{equation}

To estimate $\sup{\left|\delta\hat{h}\right|}$, we suggest considering the underlying data errors as independent normally distributed random variables with zero expectations and known variances:
\begin{equation}\label{f15}
    \delta\hat{\tau}\sim\mathcal{N}\left(0,{\hat{\sigma}}_\tau^2\right),\quad\delta{\hat{p}}_s\sim\mathcal{N}\left(0,{\hat{\sigma}}_s^2\right),\quad\delta{\hat{p}}_r\sim\mathcal{N}\left(0,{\hat{\sigma}}_r^2\right).
\end{equation}
Assuming small data errors, $\delta\hat{h}$ also can be estimated as a normal variable:
\begin{equation}\label{f16}
    \delta\hat{h}\sim\mathcal{N}\left(0,{\hat{\sigma}}_h^2\right),\quad{\hat{\sigma}}_h^2=\left(\dfrac{\partial\hat{h}}{\partial\hat{\tau}}{\hat{\sigma}}_\tau\right)^2+\left(\dfrac{\partial\hat{h}}{\partial{\hat{p}}_s}{\hat{\sigma}}_s\right)^2+\left(\dfrac{\partial\hat{h}}{\partial{\hat{p}}_r}{\hat{\sigma}}_r\right)^2.
\end{equation}
Finally, we apply the “three sigma rule” to estimate the optimal $\hat{\alpha}$ value:
\begin{equation}\label{f17}
    \hat{\alpha}\approx3{\hat{\sigma}}_h\hat{v}\dfrac{\partial\hat{h}}{\partial\hat{v}}.
\end{equation}

We conclude that given a single event in homogeneous medium one can estimate the subsurface offset’s variance \eqref{f16} and the optimal regularization parameter \eqref{f17}. We suggest considering the homogeneity restriction \eqref{f11} formally and interpret $\hat{v}$ just as a realistic velocity estimate. We expect \eqref{f16} and \eqref{f17} to be accurate at least in order of magnitude in heterogeneous models. We provide more details on the errors' distribution in the Numerical test section and list the explicit formulae for the $\hat{h}$'s derivatives in the Appendix A.

\subsection{Nonlinear CDR method for multiple rays in heterogeneous media}

Let us return to the general case of heterogeneous velocity models parametrized by a matrix $\mathbf{V}$ \eqref{f7}. We shall build an iterative algorithm of velocity updates and denote the parameters’ matrix at i-th iteration by $\mathbf{V}^{\left(i\right)}$. In particular, $\mathbf{V}^{\left(0\right)}$ will denote the initial guess matrix. Unlike the previous section, we shall deal with $K$ events, each defining a DSR ray. We shall mark all the quantities related to $k$-th event with a corresponding sub- or superscript. For instance, $h_k$ will denote the $k$-th subsurface half-offset.

We suggest minimizing weighted and regularized loss functional:
\begin{equation}\label{f18}
    L\left(\mathbf{V}\right)=\sum_{k=1}^{K}{w_k^2h_k^2}+\alpha\dfrac{\Vert v\left(x,z;\mathbf{V}-\mathbf{V}^{\left(0\right)}\right) \Vert^2}{\Vert v\left(x,z;\mathbf{V}^{\left(0\right)}\right) \Vert^2}\rightarrow\min_{\mathbf{V}}
\end{equation}
with the weights $w_k$ and the regularization parameter $\alpha$ chosen so that the $k$-sum is dominated by high-accuracy events while the regularizing term prohibits excessively complex models.

We are going to use estimates \eqref{f16} and \eqref{f17} in their “formal” interpretation. For each ray we prescribe data errors \eqref{f15} and define ${\hat{v}}_k$ as the average of the source and receiver initial velocities:
\begin{equation}\label{f19}
    {\hat{v}}_k=\dfrac{v\left({\hat{x}}_s^k,\hat{z}=0;\mathbf{V}^{\left(0\right)}\right)+v\left({\hat{x}}_r^k,\hat{z}=0;\mathbf{V}^{\left(0\right)}\right)}{2},
\end{equation}
Then we proceed with the following multipliers for \eqref{f18}:
\begin{equation}\label{f20}
    w_k=\dfrac{\dfrac{1}{{\hat{\sigma}}_h^k}}{\sqrt{\sum_{k'=1}^{K}\left(\dfrac{1}{{\hat{\sigma}}_h^{k'}}\right)^2}},\quad\alpha=\sum_{k'=1}^{K}{w_{k'}^2{\hat{\alpha}}_{k'}}.
\end{equation}
Essentially, we divided each $h_k$ by its expected variance ${\hat{\sigma}}_h^k$ \eqref{f16} and formed up $\alpha$ by individual rays’ contributions ${\hat{\alpha}}_k$ \eqref{f17}. Sub- and superscripts $k$ and $k'$ refer to the ray’s number, not power. Such weights $w_k$ favor high-accuracy rays while $\alpha$ decreases with data errors. We expect the loss functional \eqref{f18} to be well-behaved, allowing one to minimize it by gradient-based methods. We emphasize that $w_k$ and $\alpha$ are estimated prior to inversion.

\section{Practical implementation}

\subsection{Model parametrization}

We required the model $v\left(x,z;\mathbf{V}\right)$ \eqref{f7} to be twice smooth w.r.t. $x$ and $z$ and linear w.r.t. $\mathbf{V}$. We shall use gridded velocity models supplied with bicubic interpolation \cite{22}. We introduce a fixed inversion grid and interpret $\mathbf{V}$ as a matrix of grid velocities
\begin{equation}\label{f21}
    \vec{X}=\left(x_1,\ldots,x_M\right)^T,\quad\vec{Z}=\left(z_1,\ldots,z_N\right)^T,\quad v_{mn}=v\left(x_m,z_n;\mathbf{V}\right).
\end{equation}
In this paper we assume that the grid step is constant:
\begin{equation}\label{f22}
    \Delta x = x_{m + 1} - x_m = const., \quad \Delta z = z_{n + 1} - z_n = const.
\end{equation}

A bicubic spline at $x\in\left[x_m,x_{m+1}\right]$, $z\in\left[z_n,z_{n+1}\right]$ is given by
\begin{equation}\label{f23}
    \left.v\left(x,z;\mathbf{V}\right)\right|_{\substack{x\in\left[x_m,x_{m+1}\right] \\ z\in\left[z_n,z_{n+1}\right]}}=\Delta\vec{X}_m^T \mathbf{C}_{mn} \Delta\vec{Z}_n
\end{equation}
with
\begin{equation}\label{f24}
    \begin{gathered}
        \Delta\vec{X}_m=\left(\begin{matrix}x-x_m\\x_{m+1}-x\\\left(x-x_m\right)^3\\\left(x_{m+1}-x\right)^3\\\end{matrix}\right),\quad\Delta\vec{Z}_n=\left(\begin{matrix}z-z_n\\z_{n+1}-z\\\left(z-z_n\right)^3\\\left(z_{n+1}-z\right)^3\\\end{matrix}\right), \\ \mathbf{C}_{mn}=\left(\begin{matrix}c_{mn}^{11}&\cdots&c_{mn}^{14}\\\vdots&\ddots&\vdots\\c_{mn}^{41}&\cdots&c_{mn}^{44}\\\end{matrix}\right).
    \end{gathered}
\end{equation}
Coefficient matrices $\mathbf{C}_{mn}$ linearly depend on the grid velocities $\mathbf{V}$, see the Appendix B for details.

Bicubic splines belong to the Sobolev space of twice differentiable functions with a norm defined as an integral of their squared values and derivatives. We propose a finite-sum variant of this norm:
\begin{equation}\label{f25}
    \begin{gathered}
        \Vert v(x, z; \mathbf{V}) \Vert^2 = \sum_{m=1}^{M}\sum_{n=1}^{N}{v_{mn}^2 + \left(\left.\dfrac{\partial v}{\partial x}\right|_{\substack{x = x_m \\ z = z_n}} \Delta x\right) ^2 +} \\
        + \left(\left.\dfrac{\partial v}{\partial z}\right|_{\substack{x = x_m \\ z = z_n}} \Delta z\right) ^2 + \left(\left.\dfrac{\partial^2v}{{\partial x}^2}\right|_{\substack{x = x_m \\ z = z_n}} \Delta x^2\right) ^2 + \\
        + \left(\left.\dfrac{\partial^2v}{{\partial z}^2}\right|_{\substack{x = x_m \\ z = z_n}} \Delta z^2\right) ^2 + 2 \left(\left.\dfrac{\partial^2v}{{\partial x \partial z}}\right|_{\substack{x = x_m \\ z = z_n}} \Delta x \Delta z\right) ^2.
    \end{gathered}
\end{equation}
When inserted in \eqref{f18}, this norm will penalize large deviations from the initial model, fast growth and fast oscillations in both directions, like the “damping” and “smoothing” terms in other slope tomography algorithms \cite{3, 4, 23}. In homogeneous models this norm simplifies to the squared velocity used in \eqref{f13} which partially justifies the use of \eqref{f16} and \eqref{f17} in heterogeneous media.

\subsection{Gauss-Newton minimization}

Given a set of events and an initial model, we shall use the Gauss-Newton method to minimize the loss functional \eqref{f18}, following \cite{3, 4}. In this method, the problem is linearized at each iteration, and a perturbation $\delta\vec{V}=\left(\delta v_{11},\delta v_{12},\ldots,\delta v_{MN}\right)^T$ is searched as the best-fit solution for a system of linear equations:
\begin{equation}\label{f26}
    \mathbf{G}\delta\vec{V}=-\vec{l},
\end{equation}
where $\vec{l}$ consists of the loss terms and $\mathbf{G}$ consists of their gradients w.r.t. velocity parameters, see the Appendix C for explicit expressions. System \eqref{f26} is solved in the least-squares sense:
\begin{equation}\label{f27}
    \delta\vec{V}=-\left(\mathbf{G}^T\mathbf{G}\right)^{-1}\mathbf{G}^T\vec{l}
\end{equation}
and the model is updated in the estimated direction. However, due to nonlinearity of the original problem, a line search step is required. Assuming $i$-th iteration, matrix $\mathbf{V}$ is updated as follows:
\begin{equation}\label{f28}
    \begin{gathered}
        \mathbf{V}^{\left(i+1\right)}=\mathbf{V}^{\left(i\right)}+\varepsilon\cdot\delta\mathbf{V}, \\
        \varepsilon=\underset{0\le\varepsilon\le1}{\mathrm{argmin}}\;{L\left(\mathbf{V}^{\left(i\right)}+\varepsilon\cdot\delta\mathbf{V}\right)},\quad\delta\mathbf{V}=\left(\begin{matrix}\delta v_{11}&\cdots&\delta v_{1N}\\\vdots&\ddots&\vdots\\\delta v_{M1}&\cdots&\delta v_{MN}\\\end{matrix}\right).
    \end{gathered}
\end{equation}
The iterations are halted as soon as the loss stops decreasing.

\subsection{Data picking and selection}

We shall use the slant stacks to pick the events in seismic data, following \cite{7} with some minor modifications. This method scans the data and returns events ${\vec{d}}_k$ specified by source and receiver coordinates, traveltimes, their derivatives, so-called semblances (events' coherency \cite{24}) and the numbers of traces falling in the analysis window (the folds):
\begin{equation}\label{f29}
    {\vec{d}}_k=\left({\hat{x}}_s^k,{\hat{x}}_r^k,{\hat{\tau}}_k,{\hat{p}}_s^k,{\hat{p}}_r^k,{\hat{s}}_k,{\hat{Q}}_k\right),\quad{\hat{s}}_k\in\left[0,1\right],\quad{\hat{Q}}_k>1.
\end{equation}
While the first five components form up the initial conditions for the ray tracing system \eqref{f5}, the last two bring information on the picking error. We expect the events with higher semblances and folds to be more accurate, see the next section for details.

As highlighted in \cite{1}, it is important to sort out some events during the inversion due to non-uniform ray coverage of the model. In our parametrization the model’s variability is controlled by the grid spacing \eqref{f22}. We suggest an ad-hoc selection algorithm:
\begin{enumerate}
    \item Trace all the DSR equation rays downward and estimate the subsurface midpoints’ positions \eqref{f6};
    \item Select at most $10$ events whose subsurface midpoints fall into each model cell, prioritizing higher weight \eqref{f20} events;
    \item Proceed by the Gauss-Newton scheme using the selected events only;
    \item Repeat the selection at the beginning of each inversion iteration.
\end{enumerate}
The factor $10$ at the second step is heuristic and implies that the number of processed events should be greater than that of velocity parameters. We emphasize that the proposed selection is performed at each iteration, leading to possible instability of the loss-based stopping criterion.

\subsection{Programming issues}

We implemented our algorithm using NumPy \cite{25}, SciPy \cite{26} and Numba \cite{27} libraries of the Python language. We employed FteikPy \cite{28} library to compute traveltimes for the Kirchhoff migration and Matplotlib \cite{29} library to plot the results, see the next section for details.

\section{Numerical Tests}

\subsection{Tests' layout}

To test our formulation of the CDR method we designed a hierarchy of velocity models revolving around the Marmousi dataset \cite{30} which is a classic benchmark for reflection seismology and slope tomography in particular \cite{5, 6, 7}. The models’ hierarchy is displayed on the Fig. \ref{fig:models}. It contains four models: the original Marmousi model, its smoothed version (Marmousoft), its best linear fit (Marmoulin) and a step-like approximation of the latter one (Marmouflat). The first model is publicly available \cite{31}, the second one was proposed in \cite{7} and kindly provided to us by Gilles Lambaré, and the last two models were created by ourselves. We shall test our method on the Marmousi and Marmousoft models, passing the Marmoulin model as the initial guess, and using the Marmouflat model to estimate the picking errors.

\begin{figure}[h]
\includegraphics[width=12.5cm]{models_hierarchy_eng_article.png}
\caption{Models' hierarchy}
\label{fig:models}
\end{figure}

We borrowed the survey geometry and time sampling from \cite{7}. The survey line consisted of $261$ sources with $96$ receivers per source. The sources and receivers spacing was $25$ m, with sources ranging from $2600$ to $9100$ m and source-receiver offsets ranging from $-2475$ to $-100$ m. Time sampling interval was $4$ ms. We  simulated the Marmousi data in the Madagascar package \cite{32}, received the Marmousoft data from Gilles Lambaré \cite{7} and simulated the Marmouflat data using kinematic ray tracing followed by a convolution with the Ricker wavelet \cite{33}. Three common-shot gathers are displayed on the Fig. \ref{fig:data}. Despite being oversimplified, the Marmouflat model exhibits certain similarities with its more complex counterparts and thus can be used to estimate the picking accuracy.

\begin{figure}[h]
\includegraphics[width=12.5cm]{models_gathers_eng_article.png}
\caption{Common-shot gathers in different models}
\label{fig:data}
\end{figure}

\subsection{Picking accuracy}

We picked events in the data and filtered out low-semblance and low-fold events \eqref{f29}. We display the picking results on the Fig. \ref{fig:events}. We can describe the picking quality as poor in the Marmousi model, satisfying in the Marmousoft model and good in the Marmouflat model.

\begin{figure}[h]
\includegraphics[width=12.5cm]{models_gathers_events_eng_article.png}
\caption{Common-shot gathers in different models with picking results}
\label{fig:events}
\end{figure}

We compared the picked Marmouflat events to the ray-modelled ones. We estimate the traveltime error’s variance \eqref{f15} as a single time sampling interval for all the events at once. In contrast, we link the slopes’ errors to events’ semblances and folds \eqref{f29}, individual for each event. We propose the following estimates:
\begin{equation}\label{f30}
    {\hat{\sigma}}_\tau^k\approx4\cdot{10}^{-3}\ \left[s\right],\quad{\hat{\sigma}}_s^k={\hat{\sigma}}_r^k\approx\dfrac{1}{{\hat{v}}_k}\dfrac{34.2{\hat{s}}_k}{\left(-4.9+1.6{\hat{Q}}_k\right)}\ \left[\dfrac{s}{m}\right],
\end{equation}
where the sub- and superscripts $k$ refer to the event’s number, not power. We obtained numeric multipliers in ${\hat{\sigma}}_s^k={\hat{\sigma}}_r^k$ using maximum likelihood method. To check them, we divided the actual slope errors by their estimated variances and built a histogram fitting theoretical probability density function of the standard normal distribution, see the Fig. \ref{fig:errors}. We multiplied the obtained estimates ${\hat{\sigma}}_s^k={\hat{\sigma}}_r^k$ by a factor of $3$ while working with the Marmousi and Marmousoft models to consider their complexity. Average estimated slope error variance in these models amounted for $3.6\cdot{10}^{-5}$ s/m and $1.5\cdot{10}^{-5}$ s/m respectively which accords with the a-priori estimate $1\cdot{10}^{-5}$ s/m used in \cite{7}.

\begin{figure}[h]
\includegraphics[width=6.5cm]{histogram_eng_article.png}
\caption{Normalized slope error distribution}
\label{fig:errors}
\end{figure}

\subsection{Inversion results}

We passed the Marmousi and Marmousoft events to nonlinear inversion, using inversion grids \eqref{f21} with $100$ m spacing in both directions. In both cases we started from the Marmoulin model and performed $6$ Gauss-Newton iterations to reach the stopping criterion. We present the results on the Fig. \ref{fig:marmousi} and Fig. \ref{fig:marmousoft}.

\begin{figure}[h]
\includegraphics[width=12.5cm]{marmousi_results_eng_article.png}
\caption{Marmousi model inversion}
\label{fig:marmousi}
\end{figure}

\begin{figure}[h]
\includegraphics[width=12.5cm]{marmousoft_results_eng_article.png}
\caption{Marmousoft model inversion}
\label{fig:marmousoft}
\end{figure}

Evidently, our tomography failed to reconstruct the Marmousi model, although it did image some of its principal structures. We attribute this failure to the model’s blocky structure, which is not representable in the class of bicubic splines \eqref{f23}, and poor picking quality. The bottom part of the model was almost unilluminated by rays, obstructing its successful resolution. In contrast, smooth Marmousoft model is fairly well-resolved down to its bottom part, thanks to its dense ray coverage. In this model our results are competitive with \cite{5, 6, 7}.

To further assess the tomography results, we performed first-arrival Kirchhoff depth migration \cite{15} of the Marmousi and Marmousoft datasets. We present the resulting images on the Fig. \ref{fig:migration}, along with the true model migration results as a reference. Although being inexact, the inverted models allow for interpretative imaging even in the Marmousi case. We admit that the bottom part of the Marmousi image remains unfocused after inversion, but we remind that we used first-arrival migration, which does not consider all the raypaths and cannot image reflectors in such complex media \cite{34}, as the true model migration suggests.

\begin{figure}[h]
\includegraphics[width=12.5cm]{migration_eng_article_3.png}
\caption{Migration results in different models}
\label{fig:migration}
\end{figure}

Finally, we checked our regularization parameter’s estimate \eqref{f20} using the L-curve method \cite{35}. To do so, we performed a sequence of Marmousi data inversions, setting $\alpha$ to different fractions of its estimated value. Then we evaluated the models’ norms \eqref{f25} and residuals $\sum_{k=1}^{K}{w_k^2h_k^2}$ at final iteration of each inversion run and plotted the squared norms as a function of the residuals, see Fig. \ref{fig:alpha}. We attribute the non-monotonous behaivour of the curve to the events selection algorithm. The residuals were computed for different sets of rays, producing such a non-regular output. However, our $\alpha$ estimate appears to be close to the curve’s bending point, thus being at least suboptimal.

\begin{figure}[h]
\includegraphics[width=8cm]{marmousi_alpha_assessment_eng_article.png}
\caption{L-curve for the Marmousi model}
\label{fig:alpha}
\end{figure}

\section{Discussion and conclusion}

In our paper we have presented a reformulation of the CDR method based on our asymptotics of the DSR equation. We provided theoretical aspects of our approach and developed several heuristic or semi-heuristic tricks to obtain a stable slope tomography algorithm competitive with other techniques. Our most important contribution is the suboptimal regularization of the CDR inverse problem which mitigates the notorious instability of the Sword’s method. The key steps towards this result were the DSR equation-based wavefield parametrization, homogeneous medium assumption and linearization w.r.t. both velocity model and data errors. Our approach is parsimonious (compared to stereotomography \cite{4}), does not rely on the single-arrival assumption (compared to adjoint-state method \cite{5, 6}) and does not require multiple data migrations (compared to residual moveout-based methods \cite{36, 37}). On the other hand, it does not process horizontal rays which limits its applicability in complex media with steep reflectors and contrastive layers \cite{1}. However, it manages to construct acceptable velocity models for depth migration and allows for improvements like more elaborated picking engine \cite{38, 39, 30}, directional smoothing \cite{23} and multiscale inversion \cite{1, 5, 6}. In addition, the DSR parametrization theoretically allows for lifting the nowhere-horizontal rays restriction by means of a special coordinate system \cite{41}. As a conclusion, we add that it is scalable to 3D settings and anisotropic problems \cite{42} as well.

\section{Acknowledgements}

This research is supported by grant FSUS-2025-0015 of the Ministry of Science and Higher Education of the Russian Federation. I thank my scientific supervisor, Anton Duchkov, for the problem statement.

\bigskip
\begin{thebibliography}{1} % Ссылки на литературу приводятся только на английском языке (даже если у источника нет англоязычной версии).


\bibitem{1} G. Lambaré,
\href{https://doi.org/10.1190/1.2952039}{\it Stereotomography}, GEOPHYSICS, {\bf 73} (2008), VE25–VE34.

\bibitem{2} L.A. Ryabinkin,
{\it Voprosy Reguliruemogo Napravlennogo Priyoma (RNP) Sejsmicheskih Voln}, Gostoptekhizdat: Moskva, 1957

\bibitem{3} C.H. Sword,
{\it Tomographic Determination of Interval Velocities from Reflection Seismic Data: The Method of Controlled Directional Reception}, Stanford University, 1987.

\bibitem{4} F. Billette, G. Lambaré,
\href{https://doi.org/10.1046/j.1365-246X.1998.00632.x}{\it Velocity Macro-Model Estimation from Seismic Reflection Data by Stereotomography}, Geophys. J. Int., {\bf 135} (1998), 671–690.

\bibitem{5} F.B. Tavakoli, S. Operto, A. Ribodetti, J. Virieux,
\href{https://doi.org/10.1093/gji/ggx111}{\it Slope Tomography Based on Eikonal Solvers and the Adjoint-State Method}, Geophys. J. Int., {\bf 209} (2017), 1629–1647.

\bibitem{6} S. Sambolian, S. Operto, A. Ribodetti, F.B. Tavakoli, J. Virieux,
\href{https://doi.org/10.1093/gji/ggz150}{\it Parsimonious Slope Tomography Based on Eikonal Solvers and the Adjoint-State Method}, Geophys. J. Int., {\bf 218} (2019), 456–478.

\bibitem{7} F. Billette, S.le Bégat, P. Podvin, G. Lambaré,
\href{https://doi.org/10.1190/1.1581072}{\it Practical Aspects and Applications of 2D Stereotomography}, Geophysics, {\bf 68} (2003), 1008–1021.

\bibitem{8} E. Chalard, P. Podvin, S.le Bégat, P. Berthet, B. David,
\href{https://doi.org/10.1190/1.1817421}{\it 3D Stereotomographic Inversion On Real Data Set}, SEG Technical Program Expanded Abstracts, (2002), 946-948.

\bibitem{9} B.S.S. Barbosa, J.C. Costa, E.N.S. Gomes, J. Schleicher,
\href{https://doi.org/10.1190/1.2370243}{\it Sensitivity Analysis for Stereotomography in Elliptic and Anelliptic Media}, SEG Technical Program Expanded Abstracts, (2006), 3413-3417.

\bibitem{10} SLB,
\href{https://www.slb.com/products-and-services/innovating-in-oil-and-gas/reservoir-characterization/seismic/seismic-imaging/earth-model-building/cip-tomo}{\it Common Image Point (CIP) Tomography} (accessed on 18 July 2025).

\bibitem{11} Viridien,
\href{https://www.viridiengroup.com/expertise/subsurface-imaging/geovation#1702591434-2582096304}{\it Geovation} (accessed on 18 July 2025).

\bibitem{12} S. Sambolian,
{\it Kinematically Consistent Slope Tomography Using Eikonal Solvers and the Adjoint-State Method: Theory and Applications to Velocity Model Building and Event Location}, Université Côte d’Azur, 2021.

\bibitem{13} J.F. Claerbout,
{\it Imaging the Earth’s Interior}, Blackwell Scientific Publications, 1985.

\bibitem{14} A.V. Belonosova, A.S. Alekseev,
{\it Ob Odnoj Postanovke Obratnoj Kinematicheskoj Zadachi Sejsmiki Dlya Dvumernoj Neodnorodnoj Sredy}, Nekotorye metody i algoritmy interpretacii geofizicheskih dannyh, (1967), 137–154.

\bibitem{15} B.L. Biondi,
{\it 3D Seismic Imaging}, Society of Exploration Geophysicists, 2006.

\bibitem{16} P.C. Sava, S. Fomel,
\href{https://doi.org/10.1190/1.1581078}{\it Angle‐domain Common‐image Gathers by Wavefield Continuation Methods}, GEOPHYSICS, {\bf 68} (2003), 1065–1074.

\bibitem{17} A. Khoury, W. Symes, P. Williamson, P. Shen,
\href{https://doi.org/10.1190/1.2370028}{\it DSR Migration Velocity Analysis by Differential Semblance Optimization}, SEG Technical Program Expanded Abstracts, (2006), 2450–2454.

\bibitem{18} N.N. Shilov, A.A. Duchkov,
\href{https://doi.org/10.1134/S1990478924010137}{\it Migration Velocity Analysis Using a Ray Method Asymptotics of the Double Square Root Equation}, J. Appl. Ind. Math., {\bf 18} (2024), 150–166.

\bibitem{19} N.N. Shilov, A.A. Duchkov,
\href{https://doi.org/10.3390/jmse12040636}{\it Asymptotic Ray Method for the Double Square Root Equation}, J. Mar. Sci. Eng., {\bf 12} (2024), 636.

\bibitem{20} A.A. Duchkov, M.V. de Hoop,
\href{https://doi.org/10.1190/1.3459955}{\it Extended Isochron Rays in Prestack Depth (Map) Migration}, Geophysics, {\bf 75} (2010), S139–S150.

\bibitem{21} A.N. Tikhonov, B.Ya. Arsenin,
{\it Metody Resheniya Nekorrektnyh Zadach, 3-e ispr. izd.}, Nauka. Glavnaya redakciya fiziko-matematicheskoj literatury: Moskva, 1986.

\bibitem{22} G.I. Marchuk,
{\it Metody Vychislitel'noj Matematiki}, Nauka. Glavnaya redakciya fiziko-matematicheskoj literatury: Moskva, 1980.

\bibitem{23} J.C. Costa, F.J. da Silva, E.N. Gomes, J. Schleicher, L.A. Melo, D. Amazonas, 
\href{https://doi.org/10.1190/1.2967499}{\it Regularization in Slope Tomography}, GEOPHYSICS, {\bf 73} (2008), VE39–VE47.

\bibitem{24} M.T. Taner, F. Koehler, R.E. Sheriff, 
\href{https://doi.org/10.1190/1.1440994}{\it Complex Seismic Trace Analysis}, GEOPHYSICS, {\bf 44} (1979), 1041–1063.

\bibitem{25} C.R. Harris, K.J. Millman, S.J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N.J. Smith et al., 
\href{https://doi.org/10.1038/s41586-020-2649-2}{\it Array Programming with NumPy}, Nature, {\bf 585} (2020), 357–362.

\bibitem{26} P. Virtanen, R. Gommers, T.E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright et al., 
\href{https://doi.org/10.1038/s41592-019-0686-2}{\it SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python}, Nat. Methods, {\bf 17} (2020), 261–272.

\bibitem{27} S.K. Lam, A. Pitrou, S. Seibert, 
\href{https://doi.org/10.1145/2833157.2833162}{\it Numba: a LLVM-based Python JIT compiler}, LLVM '15: Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC, (2015), 1–6.

\bibitem{28} GitHub,
\href{https://github.com/keurfonluu/fteikpy}{\it Fteikpy: Accurate Eikonal Solver for Python} (accessed on 18 July 2025).

\bibitem{29} J.D. Hunter, 
\href{https://doi.org/10.1109/MCSE.2007.55}{\it Matplotlib: A 2D Graphics Environment}, Comput. Sci. \& Eng., {\bf 9} (2007), 90–95.

\bibitem{30} A. Brougois, M. Bourget, P. Lailly, M. Poulet, P. Ricarte, R. Versteeg, 
\href{https://doi.org/10.3997/2214-4609.201411190}{\it Marmousi, Model and Data}, Proceedings of the EAEG Workshop - Practical Aspects of Seismic Data Inversion; European Association of Geoscientists \& Engineers, 1990.

\bibitem{31} SW3D,
\href{https://seis.karlov.mff.cuni.cz/software/marmousi/marmousi.htm}{\it Marmousi model and data set} (accessed on 18 July 2025).

\bibitem{32} S. Fomel, P. Sava, I. Vlad, Y. Liu, V. Bashkardin, 
\href{https://doi.org/10.5334/jors.ag}{\it Madagascar: Open-Source Software Project for Multidimensional Data Analysis and Reproducible Computational Experiments}, J. Open Res. Softw., {\bf 1} (2013), e8.

\bibitem{33} V. \u{C}erven\'{y},
{\it Seismic Ray Theory}, Cambridge University Press: Cambridge, 2001.

\bibitem{34} M.S. Operto, S. Xu, G. Lambaré, 
\href{https://doi.org/10.1190/1.1444814}{\it Can We Quantitatively Image Complex Structures with Rays?}, GEOPHYSICS, {\bf 65} (2000), 1223–1238.

\bibitem{35} P.C. Hansen, D.P. O’Leary, 
\href{https://doi.org/10.1137/0914086}{\it The Use of the L-Curve in the Regularization of Discrete Ill-Posed Problems}, SIAM J. Sci. Comput., {\bf 14} (1993), 1487–1503.

\bibitem{36} H. Chauris, M.S. Noble, G. Lambaré, P. Podvin, 
\href{https://doi.org/10.1190/1.1500382}{\it Migration Velocity Analysis from Locally Coherent Events in 2‐D Laterally Heterogeneous Media, Part I: Theoretical Aspects}, GEOPHYSICS, {\bf 67} (2002), 1202–1212.

\bibitem{37} H. Chauris, M.S. Noble, G. Lambaré, P. Podvin, 
\href{https://doi.org/10.1190/1.1500383}{\it Migration Velocity Analysis from Locally Coherent Events in 2‐D Laterally Heterogeneous Media, Part II: Applications on Synthetic and Real Data}, GEOPHYSICS, {\bf 67} (2002), 1213–1224.

\bibitem{38} S. Nguyen, R. Baina, M. Alerini, G. Lambaré, V. Devaux, M. Noble, 
\href{https://doi.org/10.1111/j.1365-2478.2008.00699.x}{\it Stereotomography Assisted by Migration of Attributes}, Geophys. Prospect., {\bf 56} (2008), 613–625.

\bibitem{39} T. Barros, R. Lopes, H. Chauris, 
\href{https://doi.org/10.1111/1365-2478.13181}{\it Towards a More Robust Input for Stereotomography}, Geophys. Prospect., {\bf 70} (2022), 502–524.

\bibitem{40} F. Brandolin, M. Ravasi, T. Alkhalifah, 
\href{https://doi.org/10.1190/geo2023-0323.1}{\it PINNslope: Seismic Data Interpolation and Local Slope Estimation with Physics Informed Neural Networks}, GEOPHYSICS, {\bf 89} (2024), V331–V345.

\bibitem{41} C.C. Stolk, M.V. de Hoop, W.W. Symes, 
\href{https://doi.org/10.1190/1.3256285}{\it Kinematics of Shot-Geophone Migration}, GEOPHYSICS, {\bf 74} (2009), WCA19–WCA34.

\bibitem{42} C. Wu, B. Feng, H. Wang, T. Wang, 
\href{https://doi.org/10.1007/s11600-020-00450-z}{\it Three-Dimensional Angle-Domain Double-Square-Root Migration in VTI Media for the Large-Scale Wide-Azimuth Seismic Data}, Acta Geophys., {\bf 68} (2020), 1021–1037.

\end{thebibliography}

\appendix{\textbf{Appendix A}}

In this Appendix we provide explicit formulae for the subsurface offset’s $\hat{h}$ \eqref{f12} derivatives appearing in \eqref{f13}–\eqref{f17}. We define
\begin{aequation}\label{a1}
    {\hat{C}}_\tau=C_\tau\left(\widehat{\vec{x}},\widehat{\vec{p}}\right).
\end{aequation}
With this shorthand:
\begin{aequation}\label{a2}
    \dfrac{\partial\hat{h}}{\partial\hat{\tau}}=-\left(\dfrac{{\hat{p}}_r}{\sqrt{\dfrac{1}{{\hat{v}}^2}-{\hat{p}}_r^2}}-\dfrac{{\hat{p}}_s}{\sqrt{\dfrac{1}{{\hat{v}}^2}-{\hat{p}}_s^2}}\right)\dfrac{{\hat{C}}_\tau}{2},
\end{aequation}
\begin{aequation}\label{a3}
    \dfrac{\partial\hat{h}}{\partial{\hat{p}}_s}=\left(\dfrac{1}{\sqrt{\dfrac{1}{{\hat{v}}^2}-{\hat{p}}_s^2}}+\dfrac{{\hat{p}}_s^2}{\left(\dfrac{1}{{\hat{v}}^2}-{\hat{p}}_s^2\right)^\frac{3}{2}}\right)\dfrac{{\hat{C}}_\tau\hat{\tau}}{2}+\dfrac{{\hat{p}}_s}{\sqrt{\dfrac{1}{{\hat{v}}^2}-{\hat{p}}_s^2}}\dfrac{\partial{\hat{C}}_\tau}{\partial{\hat{p}}_s}\dfrac{\hat{\tau}}{2},
\end{aequation}
\begin{aequation}\label{a4}
    \dfrac{\partial\hat{h}}{\partial{\hat{p}}_r}=-\left(\dfrac{1}{\sqrt{\dfrac{1}{{\hat{v}}^2}-{\hat{p}}_r^2}}+\dfrac{{\hat{p}}_r^2}{\left(\dfrac{1}{{\hat{v}}^2}-{\hat{p}}_r^2\right)^\frac{3}{2}}\right)\dfrac{{\hat{C}}_\tau\hat{\tau}}{2}-\dfrac{{\hat{p}}_r}{\sqrt{\dfrac{1}{{\hat{v}}^2}-{\hat{p}}_r^2}}\dfrac{\partial{\hat{C}}_\tau}{\partial{\hat{p}}_r}\dfrac{\hat{\tau}}{2},
\end{aequation}
\begin{aequation}\label{a5}
    \begin{gathered}
        \dfrac{\partial\hat{h}}{\partial\hat{v}}=-\dfrac{1}{{\hat{v}}^3}\left(\dfrac{{\hat{p}}_r}{\left(\dfrac{1}{{\hat{v}}^2}-{\hat{p}}_r^2\right)^\frac{3}{2}}-\dfrac{{\hat{p}}_s}{\left(\dfrac{1}{{\hat{v}}^2}-{\hat{p}}_s^2\right)^\frac{3}{2}}\right)\dfrac{{\hat{C}}_\tau\hat{\tau}}{2}- \\
        - \left(\dfrac{{\hat{p}}_r}{\sqrt{\dfrac{1}{{\hat{v}}^2}-{\hat{p}}_r^2}}-\dfrac{{\hat{p}}_s}{\sqrt{\dfrac{1}{{\hat{v}}^2}-{\hat{p}}_s^2}}\right)\dfrac{\partial{\hat{C}}_\tau}{\partial\hat{v}}\dfrac{\hat{\tau}}{2},
    \end{gathered}
\end{aequation}
where
\begin{aequation}\label{a6}
    \begin{gathered}
        \dfrac{\partial{\hat{C}}_\tau}{\partial{\hat{p}}_s}=-\dfrac{{\hat{C}}_\tau^2}{{\hat{v}}^2}\dfrac{{\hat{p}}_s}{\left(\dfrac{1}{{\hat{v}}^2}-{\hat{p}}_s^2\right)^\frac{3}{2}},\quad\dfrac{\partial{\hat{C}}_\tau}{\partial{\hat{p}}_r}=-\dfrac{{\hat{C}}_\tau^2}{{\hat{v}}^2}\dfrac{{\hat{p}}_r}{\left(\dfrac{1}{{\hat{v}}^2}-{\hat{p}}_r^2\right)^\frac{3}{2}}, \\
        \dfrac{\partial{\hat{C}}_\tau}{\partial\hat{v}}=-\dfrac{2{\hat{C}}_\tau^2}{{\hat{v}}^3}\left(\dfrac{1}{\sqrt{\dfrac{1}{{\hat{v}}^2}-{\hat{p}}_s^2}}+\dfrac{1}{\sqrt{\dfrac{1}{{\hat{v}}^2}-{\hat{p}}_r^2}}\right)- \\
        - \dfrac{{\hat{C}}_\tau^2}{{\hat{v}}^5}\left(\dfrac{1}{\left(\dfrac{1}{{\hat{v}}^2}-{\hat{p}}_s^2\right)^\frac{3}{2}}+\dfrac{1}{\left(\dfrac{1}{{\hat{v}}^2}-{\hat{p}}_r^2\right)^\frac{3}{2}}\right).
    \end{gathered}
\end{aequation}

We emphasize that \eqref{a3}–\eqref{a5} linearly grow with the reflection traveltime $\hat{\tau}$. Thus, when inserted in \eqref{f16}, \eqref{f17} and \eqref{f20}, these multipliers will reduce the weights of deeper events and magnify corresponding regularization contributions. 

\appendix{\textbf{Appendix B}}

In this Appendix we present the remaining formulae for bicubic splines’ construction. For the theory see \cite{22}. We shall follow the main text’s notation \eqref{f21}–\eqref{f24}.

The coefficient matrices $\mathbf{C}_{mn}$ from \eqref{f24} linearly depend on the grid velocities:
\begin{bequation}
    \mathbf{C}_{mn}=\mathbf{D}_{x,m}^T\ \mathbf{V}\mathbf{D}_{z,n},
\end{bequation}
with
\begin{bequation}
    \begin{gathered}
        \mathbf{D}_{x,m}=\begin{pmatrix}
        \dfrac{1}{\Delta x} \begin{pmatrix} {\vec{e}}_{m+1}^T \\[6pt]
        {\vec{e}}_{m}^T \\[6pt]
        0 \\[6pt]
        0 \end{pmatrix} +
        \dfrac{\Delta x}{6} \begin{pmatrix} {- \vec{e}}_{m+1}^T \\[6pt]
        {- \vec{e}}_{m}^T \\[6pt]
        {\vec{e}}_{m+1}^T \\[6pt]
        {\vec{e}}_{m}^T \end{pmatrix}
    \end{pmatrix} ^ T, \\[6pt]
    {\vec{e}}_m=\left(\begin{matrix}0\\\vdots\\1\\\vdots\\0\\\end{matrix}\right)\gets m^{th}\,row.
    \end{gathered}
\end{bequation}
Matrices $\mathbf{A}_x$ and $\mathbf{B}_x$ are sparse matrices with nonzero elements defined as
\begin{bequation}
    \begin{gathered}
        a_{ij}^x = \begin{cases}
        1,\quad\begin{split} i = j = 1, \\ i = j = M,\end{split} \\[6pt]
        \dfrac{2 \Delta x}{3},\quad1 < i = j < M, \\[6pt]
        \dfrac{\Delta x}{3},\quad\begin{split} 2 < i = j + 1 < M, \\ 1 < i = j - 1 < M - 1,\end{split}
        \end{cases}
        \\
        b_{ij}^x = \begin{cases}
            - \dfrac{2}{\Delta x},\quad1 < i = j < M, \\[6pt]
            \dfrac{1}{\Delta x},\quad\begin{split} 2 < i = j + 1 < M, \\ 1 < i = j - 1 < M - 1.\end{split}
        \end{cases}
    \end{gathered}
\end{bequation}
Analogous expressions for $\mathbf{D}_{z,n}$, $\mathbf{A}_z$ and $\mathbf{B}_z$ are obtained by a change of indices and subscripts.

The coefficient matrices $\mathbf{C}_{mn}$ are required to evaluate the spline or its derivatives w.r.t spatial coordinates. Matrices $\mathbf{D}_{x,m}$ and $\mathbf{D}_{z,n}$ are required to evaluate its derivatives w.r.t. grid velocities. The coefficients $\mathbf{C}_{mn}$ are updated during the inversion along with the grid velocities while $\mathbf{D}_{x,m}$ and $\mathbf{D}_{z,n}$ remain fixed. 

\appendix{\textbf{Appendix C}}

In this Appendix we provide explicit expressions for the matrix $\mathbf{G}$ and vector $\vec{l}$ from \eqref{f26}. They read
\begin{cequation}
    \mathbf{G}=\nabla_{\vec{\mathbf{V}}}\vec{l}=\left(\frac{\partial\vec{l}}{\partial v_{11}},\frac{\partial\vec{l}}{\partial v_{12}},\ldots,\frac{\partial\vec{l}}{\partial v_{MN}}\right),\quad\vec{l}=\left(\begin{matrix}{\vec{l}}_h\\{\vec{l}}_v\\{\vec{l}}_{\partial_xv}\\{\vec{l}}_{\partial_zv}\\{\vec{l}}_{\partial_x^2v}\\{\vec{l}}_{\partial_z^2v}\\{\vec{l}}_{\partial_{xz}^2v}\\\end{matrix}\right).
\end{cequation}
We note that matrix $\mathbf{G}$ is formed by vector-columns, and no transposition is required. We shall list all the $\vec{l}$’s subvectors below. To avoid large repetitive multipliers, we denote
\begin{cequation}
    w_\alpha = \dfrac{\sqrt\alpha}{\Vert v\left(x,z;\mathbf{V}^{\left(0\right)}\right) \Vert}.
\end{cequation}
With this shorthand we have:
\begin{cequation}
    \begin{gathered}
        {\vec{l}}_h=\left(\begin{matrix}w_1h_1\\\vdots\\w_Kh_K\\\end{matrix}\right),\quad{\vec{l}}_v=w_\alpha\left(\begin{matrix}v_{11}\\\vdots\\v_{MN}\\\end{matrix}\right), \\
        {\vec{l}}_{\partial_xv}=w_\alpha \Delta x \begin{pmatrix}
        \left.\dfrac{\partial v}{{\partial x}}\right|_{\substack{x = x_1 \\ z = z_1}} \\
        \vdots \\
        \left.\dfrac{\partial v}{{\partial x}}\right|_{\substack{x = x_M \\ z = z_N}}
        \end{pmatrix}, \quad
        {\vec{l}}_{\partial_zv}=w_\alpha \Delta z \begin{pmatrix}
        \left.\dfrac{\partial v}{{\partial z}}\right|_{\substack{x = x_1 \\ z = z_1}} \\
        \vdots \\
        \left.\dfrac{\partial v}{{\partial z}}\right|_{\substack{x = x_M \\ z = z_N}}
        \end{pmatrix}, \\
        {\vec{l}}_{\partial_x^2v}=w_\alpha {\Delta x} ^ 2 \begin{pmatrix}
        \left.\dfrac{\partial^2 v}{{\partial x}^2}\right|_{\substack{x = x_1 \\ z = z_1}} \\
        \vdots \\
        \left.\dfrac{\partial^2 v}{{\partial x}^2}\right|_{\substack{x = x_M \\ z = z_N}}
        \end{pmatrix}, \quad
        {\vec{l}}_{\partial_z^2v}=w_\alpha {\Delta z} ^ 2 \begin{pmatrix}
        \left.\dfrac{\partial^2 v}{{\partial z}^2}\right|_{\substack{x = x_1 \\ z = z_1}} \\
        \vdots \\
        \left.\dfrac{\partial^2 v}{{\partial z}^2}\right|_{\substack{x = x_M \\ z = z_N}}
        \end{pmatrix}, \\
        {\vec{l}}_{\partial_{xz}^2v}=\sqrt2 w_\alpha \Delta x \Delta z \begin{pmatrix}
        \left.\dfrac{\partial^2 v}{{\partial x \partial z}}\right|_{\substack{x = x_1 \\ z = z_1}} \\
        \vdots \\
        \left.\dfrac{\partial^2 v}{{\partial x \partial z}}\right|_{\substack{x = x_M \\ z = z_N}}
        \end{pmatrix}.
    \end{gathered}    
\end{cequation}
Evidently, the squared Euclidean norm of $\vec{l}$ gives the loss functional \eqref{f18}.

\end{document}
