# Redatuming via an iterative primal-dual TRAC approach

Franck Assous, Frédéric Nataf

<figure><img src="image_1.png" /></figure>

## To cite this version:

Franck Assous, Frédéric Nataf. Redatuming via an iterative primal-dual TRAC approach. 2026. [hal-05440174v2](https://hal.science/hal-05440174v2)

HAL Id: hal-05440174

https://hal.science/hal-05440174v2

Preprint submitted on 10 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Redatuming via an iterative primal-dual TRAC approach

Franck Assous*,¹ and Frédéric Nataf†,²,³

¹Ariel University, Department of Mathematics, Israël

²Sorbonne Université, Université Paris Cité, CNRS, INRIA, Laboratoire Jacques-Louis Lions, LJLL, EPC ALPINES, F-75005 Paris, France

³Blavatnik School of Computer Science, Tel Aviv University, Israël

January 10, 2026

Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Problem setting</td><td>3</td></tr><tr><td>3</td><td>TRAC methods</td><td>4</td></tr><tr><td>3.1</td><td>TRAC in the time domain</td><td>4</td></tr><tr><td>3.2</td><td>TRAC method for the Helmholtz equation with one recorded field</td><td>6</td></tr><tr><td>3.3</td><td>TRAC Frequency domain variants</td><td>7</td></tr><tr><td>3.3.1</td><td>TRAC Frequency domain with two recorded fields</td><td>7</td></tr><tr><td>3.3.2</td><td>Dual TRAC frequency domain method</td><td>9</td></tr><tr><td>4</td><td>Primal and dual TRAC-least-square methods for the phase conjugate field</td><td>9</td></tr><tr><td>4.1</td><td>Primal TRAC-least-square</td><td>9</td></tr><tr><td>4.2</td><td>Dual TRAC-least-square</td><td>10</td></tr><tr><td>5</td><td>Primal-dual formulation of the redatuming problem</td><td>11</td></tr><tr><td>6</td><td>Analysis of the Fixed-Point algorithm</td><td>13</td></tr><tr><td>7</td><td>Numerical Results</td><td>15</td></tr></table>

<footnote>*assous@ariel.ac.il</footnote>

<footnote>†frederic.nataf@sorbonne-universite.fr</footnote>
Abstract

In inverse problems, redatuming data consists in virtually moving the
sensors from the original acquisition location to an arbitrary position. This
is an essential tool for target oriented inversion. This paper presents an
exact redatuming method based on the Time Reversal Absorbing Condi-
tions (TRAC) approach. Unlike conventional methods, our approach is
built on a convergent iterative process instead of a least-squares formula-
tion. Numerical results and comparisons with other approaches illustrate
the efficiency and convergence of our method applied to the stationary
Helmholtz equation.

1 Introduction

Inverse problems play a crucial role in imaging whether in the medical field or
in seismic exploration, and remains an active subject of research, with many
applications such as tumor detection or stroke prevention in the case of medical
imaging [24, 23]. In geophysics, seismic imaging is used as a tool for exploring
subsoil for oil, gas or other deposits [30, 14]. From the mathematical point
of view, parameter or signal estimation can be written as a partial differential
equation constrained optimization problem [29, 19, 22], which tries to minimize
the misfit between the recorded data and the reconstruction obtained by an
estimated parameter or signal.

At each step of the optimization process, this estimation is updated to get
closer to the original parameter to recover. This requires to solve at each it-
eration the so-called forward problem on the whole domain. The large size
of the domain and the number of solves make this operation time-consuming
and resource-intensive. Hence, redatuming, a target oriented approach which
consists in reducing the inverse problem to a smaller region of interest, is very
appealing. Redatuming data consists in virtually moving the sources and re-
ceivers from the original acquisition level to a new location. This method can be
employed when the original acquisition location is inadequate for imaging pur-
poses e.g., remoteness from target, irregular spatial sampling, etc. Redatuming
relies on some prior knowledge of the medium parameters between the original
and the new location which will be called overburden. The standard approach
consists in modeling the propagation of the direct waves between the surface
and the positions of the virtual sources/receivers that requires a velocity model
of the overburden as prior information, see [7, 28]. These approaches usually
are based on simplifying assumptions such as one-way wave propagation. The
standard model-based redatuming methods can successfully focus the direct
waves at the new location but do not account for surface and internal multiple
reflections in the overburden. These multiples can generate ghost arrivals in
the redatumed data set and thus ghost reflectors in the subsequent images. In
more recent papers (e.g., [25] and [27] for a review), redatuming is formulated
for the two-way wave equation. This means that events such as multiples and
refractions are included. This is made possible by formulating redatuming as an inverse problem for the full acoustic wave equation. In [26] the Marchenko equations, see also [31] and references therein, are combined with a one-way version of the Rayleigh integral representation to obtain a new redatuming scheme that handles internal as well as free-surface multiples. Our approach is different and based on the full wave equation but can be conveniently related to previously developed methods for target oriented waveform redatuming, see [17, 18] and references therein.

Inverse problems are in general ill-posed and need some regularization. In [3], we proposed an exact redatuming method which is robust with respect to noise in the data and avoids thus the not so obvious choice of a regularization function. This method is based on the TRAC (Time Reversal Absorbing Conditions, see [2, 16]) method which is a blend of time reversal techniques [15] and absorbing boundary conditions [13, 5]. In its original version, it addressed the redatuming issue for homogeneous infinite domains with inclusions and was used in [11] to solve target oriented inverse problems. In [3], we introduced its generalization to inhomogeneous media e.g., a stratified domain made of layers which creates surface and internal multiple reflections. It was based on a least-square approach.

Unlike conventional methods, the approach presented here is built on a convergent iterative process instead of a least-squares formulation. It is derived for the Helmholtz equation and utilizes the phase conjugation technique, which is the frequency-domain counterpart of the time reversal method. The principles are not limited to the frequency and can be derived in the time domain. The frequency-domain formulation is adopted here to ensure a straightforward implementation and to minimize computational costs.

More precisely, in § 2 we introduce the model redatuming problem, then in § 3 we describe some existing approaches. The novel iterative method is derived in § 4 and § 5. The convergence of the algorithm is proved in § 6. In § 7, numerical results and comparisons with other approaches illustrate the efficiency and convergence of our method.

## 2 Problem setting

Let $\Omega$ be an inhomogeneous domain (e.g. in Figure 1 with four scatterers). The domain is partitioned in a subdomain $\Omega_{top}$ with known acoustic properties, and a subdomain $B$ with unknown acoustic properties. We assume that the boundary $\partial\Omega$ of $\Omega$ is a Source-Receivers Array (SRA), denoted in the sequel by $\Gamma_{SRA}$ which gives access to the value of the solution as well as to its normal derivative. The acoustic phenomenon is modeled by the Helmholtz equation with a given frequency $\omega$. Let $f(x)$ be a given source term, we define the total field $u^T$ which satisfies in $\Omega$:

$$ -\frac{\omega^2}{c^2}u^T - \Delta u^T = f \quad \text{in } \Omega, \qquad (1) $$
<figure><img src="image_3.png" /><figcaption>Figure 1: Domain $\Omega = \Omega_{top} \cup B$, with the Source-Receivers Array (SRA) boundary $\Gamma_{SRA}$ and the redatuming boundary $\partial B = \Gamma_{RDT}$</figcaption></figure>

where $c = c(x)$ is a positive function that models the wave speed.

We consider the problem of recovering $u^T$ in the entire subdomain $\Omega_{top}$, given the source term $f$, the wave speed $c(x)$ in $\Omega_{top}$, and the Cauchy data $(u^T, \frac{\partial u^T}{\partial n})$ on the boundary $\Gamma_{SRA}$. The boundary $\partial B$, denoted hereafter as $\Gamma_{RDT}$, plays a central role in the redatuming procedure.

In the sequel, we denote $c_I$ a wave speed function in $\Omega$ which is equal to $c$ in $\Omega_{top}$ but may be different from $c$ in $B$.

# 3 TRAC methods

We recall here the TRAC method introduced in [2] for solving the redatuming problem which consists in recovering the field $u^T$ in the subdomain $\Omega_{top}$ included in $\Omega$. We start with the more intuitive time dependent wave equation in section § 3.1 and then in § 3.2 its counterpart in the frequency domain where time reversal is replaced by phase conjugation. In § 3.3, we introduce variants of these methods.

## 3.1 TRAC in the time domain

We recall here the rationale for the TRAC method as derived in [2]. In this section, the domain $\Omega$ is then the whole space $\mathbf{R}^3$. We consider the case of the acoustic wave equation with a propagation speed $c$ which is constant outside an inclusion $D$ which is a subset of a ball $B$. We also denote the time dependent total field as $u^T$ which is decomposed into an incident $u^I$ and a scattered field $u^S$ which satisfies outside the inclusion $D$:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial^2 u^S}{\partial t^2} - c^2 \Delta u^S = 0 \quad \text{in } (0, T_f) \times \mathbf{R}^3 \setminus B \\[1em] u^S \text{ satisfies a Sommerfeld condition at } \infty \\[1em] \text{homogeneous initial conditions.} \end{array} \right. \qquad (2) $$
where $T_f$ is a time large enough so that the solution vanishes.

This problem is under-determined since we do not specify any boundary condition on $\partial B$. Thus $u^S$ is not necessarily zero. We define the time-reversed total field $u_R^T$ by $u_R^T(t, \cdot) := u^T(T_f - t, \cdot)$.

We now search a boundary value problem satisfied by $u_R^T$ in $\mathbb{R}^3 \setminus B$. First, since the wave equation involves only second order time derivatives, this definition ensures that the reverse field $u_R^T$ is a solution to the wave equation as well. In order to derive a boundary condition on $\Gamma_{RDT}$, we consider, for the sake of simplicity, that the subdomain $B$ is a ball centered at the origin. Let $r$ be the radial coordinate, we consider the first order Bayliss-Turkel absorbing boundary condition (ABC) $BT^1$ [5, 4] assumed to be exact:

$$ \frac{\partial u^S}{\partial t} + c \frac{\partial u^S}{\partial r} + c \frac{u^S}{r} = 0 \text{ on } \Gamma_{RDT}. \quad (3) $$

For ease of the derivation in three dimensions, we use an equivalent form of equation (3):

$$ \frac{\partial}{\partial t}(ru^S) + c \frac{\partial}{\partial r}(ru^S) = 0 \text{ on } \Gamma_{RDT}. \quad (4) $$

We next wish to infer from (4) a boundary condition for $u_R^T$ on $\Gamma_{RDT}$. The total field $u^T = u^I + u^S$ satisfies

$$ \frac{\partial}{\partial t}(r u^T(t, \cdot)) + c \frac{\partial}{\partial r}(r u^T(t, \cdot)) = \frac{\partial}{\partial t}(r u^I(t, \cdot)) + c \frac{\partial}{\partial r}(r u^I(t, \cdot)). $$

Using $u_R^T(t, \cdot) = u^T(T_f - t, \cdot)$, we get

$$ \left(-\frac{\partial}{\partial t}(r u_R^T) + c \frac{\partial}{\partial r}(r u_R^T)\right)|_{T_f - t} = \frac{\partial}{\partial t}(r u^I(t, \cdot)) + c \frac{\partial}{\partial r}(r u^I(t, \cdot)), $$

or equivalently,

$$ -\frac{\partial}{\partial t}(r u_R^T(t, \cdot)) + c \frac{\partial}{\partial r}(r u_R^T(t, \cdot)) = \left(\frac{\partial}{\partial t}(r u^I) + c \frac{\partial}{\partial r}(r u^I)\right)|_{T_f - t}. $$

Note, that on $\Gamma_{RDT}$, $\partial/\partial r = -\partial/\partial n$ where $n$ is the outward normal to the restricted domain $\Omega \setminus B$. Multiplying by $-1$, we get

$$ \frac{\partial}{\partial t}(r u_R^T(t, \cdot)) + c \frac{\partial}{\partial n}(r u_R^T(t, \cdot)) = -\left(\frac{\partial}{\partial t}(r u^I) - c \frac{\partial}{\partial n}(r u^I)\right)|_{T_f - t}. \quad (5) $$

Another way to write this boundary condition is to introduce the time-reversed incident field

$$ u_R^I(t, \cdot) := u^I(T_f - t, \cdot), $$

so that (5) can be expressed as

$$ \frac{\partial}{\partial t}(r u_R^T(t, \cdot)) + c \frac{\partial}{\partial n}(r u_R^T(t, \cdot)) = \frac{\partial}{\partial t}(r u_R^I(t, \cdot)) + c \frac{\partial}{\partial n}(r u_R^I(t, \cdot)). \quad (6) $$

Since $\partial r/\partial n = -1$, relation (6) can be rewritten as:

$$ \frac{\partial}{\partial t}(u_R^T(t, \cdot)) + c \frac{\partial}{\partial n}(u_R^T(t, \cdot)) - c \frac{u_R^T(t, \cdot)}{r} = \frac{\partial}{\partial t}(u_R^I(t, \cdot)) + c \frac{\partial}{\partial n}(u_R^I(t, \cdot)) - c \frac{u_R^I(t, \cdot)}{r}. \quad (7) $$
So we define the boundary condition TRAC by:

$$ \mathrm{TRAC}(u_R^T) := \frac{\partial}{\partial t}(u_R^T(t, \cdot)) + c \frac{\partial}{\partial n}(u_R^T(t, \cdot)) - c \frac{u_R^T(t, \cdot)}{r}. \quad (8) $$

Note, that due to the minus sign before the term $u_R^T/r$, the TRAC (7) is not the $BT^1$ absorbing boundary condition. Assuming that the $BT^1$ absorbing condition (ABC) is exact, the time-reversed problem on the total field $u_R^T$ reads:

$$ \left\{ \begin{array}{l} \frac{\partial^2 u_R^T}{\partial t^2} - c^2 \Delta u_R^T = f(T_f - t, \cdot) \text{ in } (0, T_f) \times \Omega_{top} \\ \mathrm{TRAC}(u_R^T) = \mathrm{TRAC}(u_R^I) \text{ on } \Gamma_{RDT} \\ u_R^T(t, \cdot) = u^T(T_f - t, \cdot) \text{ on } \Gamma_{SRA} \\ \text{zero initial conditions,} \end{array} \right. \quad (9) $$

where $f$ is the given right-hand side. In practice, the $BT^1$ ABC is not exact and solving eq. (9) yields an approximation to the time reversed solution.

**Remark 1** The TRAC is not only not the standard $BT^1$ ABC but also has an “anti absorbing” term $(-cu_R^T/r)$. A natural concern arises about the well-posedness of problem (9). Although we have not developed a general theory, we prove an energy estimate for this problem in a special geometry, see [1].

The generalization of (8) to two space dimensions is straightforward, see [5, 4]. In the above computation, it is sufficient to replace $r$ by $\sqrt{r}$ and $1/r$ by $1/(2r)$ and (8) reads:

$$ \frac{\partial u^S}{\partial t} + c \frac{\partial u^S}{\partial r} + c \frac{u^S}{2r} = 0. \quad (10) $$

In the sequel, we introduce a parameter $\beta$ that depends on the dimension: it is equal to 2 in 2D and to 1 in 3D.

In the above derivation we have assumed that the boundary $\Gamma_{RDT}$ is a sphere or a circle. In the general case, it suffices to replace the term $r$ by the curvature of the boundary.

## 3.2 TRAC method for the Helmholtz equation with one recorded field

In the time-harmonic dependency, the analog to the time-reversal method is the phase conjugation technique, see e.g. [10]. Indeed, let $v(t, \boldsymbol{x})$ denote a time-dependent function and let $v_R := v(-t, \boldsymbol{x})$ denote its time reversed counterpart. Its Fourier transform denoted $\hat{v}_R(\omega, \boldsymbol{x})$ reads:

$$ \hat{v}_R(\omega, \boldsymbol{x}) = \int v(-t, \boldsymbol{x}) e^{-i\omega t} dt = \int v(t, \boldsymbol{x}) e^{i\omega t} dt = \overline{\int v(t, \boldsymbol{x}) e^{-i\omega t} dt} = \hat{v}(\omega, \boldsymbol{x}). $$

This identity shows that the phase conjugation of $\hat{v}(\omega, \boldsymbol{x})$ is the Fourier transform of the time-reversal process.
The TRAC method in the frequency domain is simply derived by taking the Fourier transform of the equations. In what follows and in order to simplify the notations, we use the same symbol for a function and its Fourier transform. Following section 3.1, we take the Fourier transform of equation (9) to get:

$$ egin{cases} -\omega^2 \overline{u^T} - c^2 \Delta \overline{u^T} = \overline{f} \quad \text{in } \Omega_{top} \\ i\omega \overline{u^T} + c \frac{\partial \overline{u^T}}{\partial n} - c \frac{\overline{u^T}}{r} = i\omega \overline{u^I} + c \frac{\partial \overline{u^I}}{\partial n} - c \frac{\overline{u^I}}{r} \quad \text{on } \Gamma_{RDT} \\ u_R^T = \overline{u^T} \quad \text{on } \Gamma_{SRA}, \end{cases} \quad (11) $$

where $\overline{u^I}$ denotes the phase conjugation of the incident field. As in the time dependent case, the solution to eq. (11) is only an approximation to $\overline{u^T}$. Here again, we emphasize that the phase conjugated absorbing condition TRAC

$$ \mathrm{TRAC}(\overline{u^T}) := i\omega \overline{u^T} + c \frac{\partial \overline{u^T}}{\partial n} - c \frac{\overline{u^T}}{r} \quad (12) $$

is not the $BT^1$ absorbing boundary condition

$$ i\omega \cdot +c \frac{\partial \cdot}{\partial n} + c \frac{\cdot}{r}. $$

As before, generalizations to two dimensions with

$$ \mathrm{TRAC} := i\omega \cdot +c \frac{\partial \cdot}{\partial n} - c \frac{\cdot}{2r} \quad (13) $$

and other shapes for $B$ are straightforward.

## 3.3 TRAC Frequency domain variants

### 3.3.1 TRAC Frequency domain with two recorded fields

In § 3.2, the recorded data is only the field $u^T$ itself. In this section, we introduce the additional use of the normal derivative of $u^T$.

More precisely, let $g^T$ be a linear combination of the recorded data:

$$ g^T := \frac{i\omega}{c} u^T + \frac{\partial u^T}{\partial n} \quad \text{on } \Gamma_{SRA}. \quad (14) $$

We define the incident field $u^I$ that satisfies:

$$ \begin{cases} -\frac{\omega^2}{c_I^2} u^I - \Delta u^I = f \quad \text{in } \Omega \\ \frac{i\omega}{c_I} u^I + \frac{\partial u^I}{\partial n} = g^T \quad \text{on } \Gamma_{SRA}. \end{cases} \quad (15) $$

We introduce the scattered field:

$$ u^S := u^T - u^I, \quad (16) $$
which satisfies the problem:

$$ 
\begin{aligned}
& -\frac{\omega^2}{c^2}u^S - \Delta u^S = 0 \text{ in } \Omega_{top} \\
& \frac{i\omega}{c}u^S + \frac{\partial u^S}{\partial n} = 0 \text{ on } \Gamma_{SRA}.
\end{aligned} 
\qquad (17) $$

Since we lack a boundary condition on the interior boundary $\Gamma_{RDT}$, this problem is under-determined. The first step of the TRAC method consists in taking the phase conjugation of equation (17):

$$ 
\begin{aligned}
& -\frac{\omega^2}{c^2}\overline{u^S} - \Delta\overline{u^S} = 0 \text{ in } \Omega_{top} \\
& -\frac{i\omega}{c}\overline{u^S} + \frac{\partial\overline{u^S}}{\partial n} = 0 \text{ on } \Gamma_{SRA}.
\end{aligned} 
\qquad (18) $$

In the following, we will make use of the fact that $\overline{u^S}$ satisfies a Helmholtz equation in $\Omega_{top}$. However, the boundary condition on $\Gamma_{SRA}$ can't be used since it corresponds to an ill-posed problem due to the minus sign before the term $\frac{i\omega}{c}$.

For this reason, the second step consists in approximating $\overline{u^S}$ as the solution $v_{PC}^S$ of a problem set in $\Omega_{top}$. For this, let us introduce $G^T$ a new linear combination of the recorded data:

$$ G^T := -\frac{i\omega}{c}u^T + \frac{\partial u^T}{\partial n} = \frac{i\omega}{c}\overline{u^T} + \frac{\partial \overline{u^T}}{\partial n} \quad \text{on } \Gamma_{SRA} \qquad (19) $$

and use the TRAC boundary operator on $\Gamma_{RDT}$ defined in eq. (13). Then, the field $v_{PC}^S$ is the solution to the following problem:

$$ 
\begin{aligned}
& -\frac{\omega^2}{c^2}v_{PC}^S - \Delta v_{PC}^S = 0 \text{ in } \Omega_{top} \\
& \frac{i\omega}{c}v_{PC}^S + \frac{\partial v_{PC}^S}{\partial n} = \left(\frac{i\omega}{c} + \frac{\partial \cdot}{\partial n}\right)(\overline{u^T} - \overline{u^I}) = G^T - \overline{g^T} - 2\frac{i\omega}{c}\overline{u^I} \quad \text{on } \Gamma_{SRA} \\
& \text{TRAC}(v_{PC}^S) = 0 \text{ on } \Gamma_{RDT}
\end{aligned} 
\qquad (20) $$

Finally, $\overline{v_{PC}^S}$ is the TRAC approximation to $u^S$ and has to be close to $u^S$.

In this paper, instead of considering scattered fields, we rather approximate the conjugate of the total field $u^T$ by solving the following phase conjugate problem which is the counterpart of equation (11) using again $u^I$ from eq. (15):

$$ 
\begin{aligned}
& -\frac{\omega^2}{c^2}v_{PC} - \Delta v_{PC} = \overline{f} \quad \text{in } \Omega_{top} \\
& \frac{i\omega}{c}v_{PC} + \frac{\partial v_{PC}}{\partial n} = G^T \quad \text{on } \Gamma_{SRA} \\
& \text{TRAC}(v_{PC}) = \text{TRAC}(\overline{u^I}) \quad \text{on } \Gamma_{RDT}
\end{aligned} 
\qquad (21) $$

so that $v_{PC}$ is an approximation to $\overline{u^T}$ in $\Omega_{top}$ or equivalently $\overline{v_{PC}}$ is an approximation to $u^T$ in $\Omega_{top}$.
### 3.3.2 Dual TRAC frequency domain method

Similarly, the conjugation of $u^T$ solution to (1) satisfies a Helmholtz equation with a conjugated right-hand side $\bar{f}$. Consequently, we can apply the same methodology for an approximate reconstruction. It is thus possible to apply the TRAC method to approximate $\bar{u}^T$. The first step also consists in defining $u_{PC}^I$ an incident field associated to the phase conjugate problem (26) solution to:

$$ 
\begin{aligned}
& -\frac{\omega^2}{c_I^2} u_{PC}^I - \Delta u_{PC}^I = \bar{f} \quad \text{in } \Omega \\
& \frac{i\omega}{c_I} u_{PC}^I + \frac{\partial u_{PC}^I}{\partial n} = G^T \quad \text{on } \Gamma_{SRA}.
\end{aligned}
\qquad (22) $$

Note that $u_{PC}^I$ is not the conjugate of $u^I$ introduced in eq. (15). Then, the counterpart of problem (21) consists in solving for $v$

$$ 
\begin{aligned}
& -\frac{\omega^2}{c^2}v - \Delta v = f \quad \text{in } \Omega_{top} \\
& \frac{i\omega}{c}v + \frac{\partial v}{\partial n} = g^T \quad \text{on } \Gamma_{SRA} \\
& \mathrm{TRAC}(v) = \mathrm{TRAC}(u_{PC}^I) \quad \text{on } \Gamma_{RDT}.
\end{aligned}
\qquad (23) $$

Finally, $v$ is an approximation to $u^T$ or equivalently, $\bar{v}$ is the TRAC approximation to $\bar{u}^T$ and has to be close to $\bar{u}^T$.

# 4 Primal and dual TRAC-least-square methods for the phase conjugate field

In this section, we introduce the primal and dual TRAC-least-square methods which allow for a possible exact redatuming at the expense of a minimization procedure.

## 4.1 Primal TRAC-least-square

Up to here, the procedures are one-shot methods that necessarily yield an error between the reconstructed and the exact fields. In order to improve this procedure, we introduced in [3] for the time dependent problem a least square formulation whose unknown is the right-hand side on $\Gamma_{RDT}$. The minimization was performed by the conjugate gradient algorithm, see [3].

In the frequency domain, it would amount to introducing the unknown function $H_{PC}$ on $\Gamma_{RDT}$ defined by

$$ H_{PC} := \mathrm{TRAC}(\bar{u}^T). \qquad (24) $$
And similarly to eq. (21), we look for the optimal value of $H_{PC}$ that minimizes the difference on $\Gamma_{SRA}$ between $u^T$ and $v_{PC}$ solution to

$$ 
\begin{aligned}
& -\frac{\omega^2}{c^2} v_{PC} - \Delta v_{PC} = \bar{f} \quad \text{in } \Omega_{top} \\
& \frac{i\omega}{c} v_{PC} + \frac{\partial v_{PC}}{\partial n} = G^T \quad \text{on } \Gamma_{SRA} \\
& \mathrm{TRAC}(v_{PC}) = H_{PC} \quad \text{on } \Gamma_{RDT}.
\end{aligned}
$$

The TRAC-ls method consists in minimizing the cost function obtained by applying the operator $\frac{i\omega}{c} \cdot + \frac{\partial}{\partial n}$ to $\bar{v}_{PC}$:

$$ J(H_{PC}) := \int_{\Gamma_{SRA}} |g_T - \left(\frac{i\omega}{c}\bar{v}_{PC} + \frac{\partial \bar{v}_{PC}}{\partial n}\right)|^2. \quad (25) $$

If the minimization is performed by an iterative solver, we see from (21) that a good initial guess is $H_{PC} = \mathrm{TRAC}(u^T)$.

## 4.2 Dual TRAC-least-square

As in § 3.3.2, we note that the conjugate of the total field also satisfies a Helmholtz equation similar to eq. (1) but with a conjugate right-hand side:

$$ -\frac{\omega^2}{c^2} \overline{u^T} - \Delta \overline{u^T} = \bar{f} \quad \text{in } \Omega. \quad (26) $$

It is thus possible to apply the same methodology as in § 4.1 in order to determine $\overline{u^T}$. Similarly to the introduction of the TRAC-ls method with the unknown $H_{PC}$ eq. (24), it is natural to introduce the unknown function $H$ on $\Gamma_{RDT}$ defined by

$$ H := \mathrm{TRAC}(u^T), \quad (27) $$

that minimizes the difference on $\Gamma_{SRA}$ between $u^T$ and $v$ solution to:

$$ 
\begin{aligned}
& -\frac{\omega^2}{c^2} v - \Delta v = f \quad \text{in } \Omega_{top} \\
& \frac{i\omega}{c} v + \frac{\partial v}{\partial n} = g^T \quad \text{on } \Gamma_{SRA} \\
& \mathrm{TRAC}(v) = H \quad \text{on } \Gamma_{RDT}.
\end{aligned}
\quad (28) $$

The dual TRAC-ls method consists in minimizing the following cost function:

$$ J(H) := \int_{\Gamma_{SRA}} |G_T - \left(\frac{i\omega}{c}\bar{v} + \frac{\partial \bar{v}}{\partial n}\right)|^2. \quad (29) $$

In the next section, we take a different path by combining primal and dual TRAC-ls methods to develop an iterative algorithm. The resulting error amplification operator has a norm less than or equal to 1. Furthermore, unlike traditional approaches, this derivation avoids the use of a cost function, reformulating the redatuming problem instead as a linear system of two unknowns on the boundary $\Gamma_{RDT}$.
5 Primal-dual formulation of the redatuming problem

The first step is to find a linear system solved by the two unknown functions $H_{PC}$ and $H$ on $\Gamma_{RDT}$ whose definitions are recalled here

$$H_{PC} := TRAC(\overline{u^T}) \text{ and } H := TRAC(u^T).$$

Recall that, $H_{PC}$ is such that the solution $v_{PC}$ of

$$\begin{gathered} -\frac{\omega^2}{c^2}v_{PC} - \Delta v_{PC} = \bar{f}(\omega, x) \text{ in } \Omega_{top} \\ \frac{i\omega}{c}v_{PC} + \frac{\partial v_{PC}}{\partial n} = G^T \text{ on } \Gamma_{SRA} \\ TRAC(v_{PC}) = H_{PC} \text{ on } \Gamma_{RDT} \end{gathered} \quad (30)$$

is equal to $\overline{u^T}$ or equivalently that $\overline{v_{PC}}$ is equal to $u^T$. Similarly, $H$ is such that the solution $v$ of

$$\begin{gathered} -\frac{\omega^2}{c^2(x)}v - \Delta v = f \text{ in } \Omega_{top} \\ \frac{i\omega}{c}v + \frac{\partial v}{\partial n} = g^T \text{ on } \Gamma_{SRA} \\ TRAC(v) = H \text{ on } \Gamma_{RDT} \end{gathered} \quad (31)$$

is equal to $u^T$ or equivalently that $\bar{v}$ is equal to $\overline{u^T}$.

In order to derive a linear (for real coefficient) system satisfied by $(H, H_{PC})$ for which we have a converging fixed point method, see § 6, we mix equations (30) and (31). We note that we would like that $v$ is equal to $\overline{v_{PC}}$ so that it makes sense to require that:

$$H = TRAC(\overline{v_{PC}}) \text{ on } \Gamma_{RDT}. \quad (32)$$

This yields a new equation linking $H$ and $H_{PC}$ since $\overline{v_{PC}}$ from (30) depends on $H_{PC}$.

In order to explicit the above equation, we introduce the operator $A_{TRAC}^{-1}$ which maps the right-hand sides $\overline{f}, G^T, H_{PC}$ to the solution of the boundary value problem (30). Using this notation, we can write $v_{PC}$ as:

$$v_{PC} = A_{TRAC}^{-1}(\overline{f}, G^T, H_{PC}).$$

By linearity of $A_{TRAC}^{-1}$, we have:

$$H - TRAC(\overline{A_{TRAC}^{-1}(0, 0, H_{PC})}) = TRAC(\overline{A_{TRAC}^{-1}(\overline{f}, G^T, 0)}) \text{ on } \Gamma_{RDT}.$$

Introducing the operator $\mathcal{S}$ defined by:

$$\mathcal{S}(\cdot) := TRAC(\overline{A_{TRAC}^{-1}(0, 0, \cdot)}), \quad (33)$$
the expression (32) reads:

$$H - S(H_{PC}) = TRAC(\overline{A_{TRAC}^{-1}(\overline{f}, G^T, 0)}) .$$

Similarly in order to find an equation for $H_{PC}$, we write:

$$H_{PC} = TRAC(\bar{v}).$$

Above, $\overline{u^T}$ is approximated by $\bar{v}$ from (31) which depends on $H$:

$$v = A_{TRAC}^{-1}(f, g^T, H)$$

We have a second equation that involves $H$ and $H_{PC}$:

$$H_{PC} - TRAC(\overline{A_{TRAC}^{-1}(0, 0, H)}) = TRAC(\overline{A_{TRAC}^{-1}(f, g^T, 0)}),$$

which reads:

$$H_{PC} - S(H) = TRAC(\overline{A_{TRAC}^{-1}(f, g^T, 0)}).$$

Overall we have found a square system for $H$ and $H_{PC}$:

$$\begin{align}
H - S(H_{PC}) &= TRAC(\overline{A_{TRAC}^{-1}(\overline{f}, G^T, 0)}) \tag{34} \\
H_{PC} - S(H) &= TRAC(\overline{A_{TRAC}^{-1}(f, g^T, 0)}). \nonumber
\end{align}$$

Note that due to the conjugation in the application of the operator $S$, the above system is not linear.

In order to simplify the implementation and bypass normal derivative computations on $\Gamma_{RDT}$, we make use of several formulas. First, on $\Gamma_{RDT}$ we have:

$$TRAC(\overline{A_{TRAC}^{-1}(\cdot, \cdot, H)}) = \overline{H} + 2 \frac{i\omega}{c} \overline{A_{TRAC}^{-1}(\cdot, \cdot, H)}|_{\Gamma_{RDT}} = \overline{H} + 2 \frac{i\omega}{c} \overline{B(\cdot, \cdot, H)}. \tag{35}$$

This simplifies the expression of $S(H)$:

$$S(H) = \overline{H} + 2 \frac{i\omega}{c} \overline{A_{TRAC}^{-1}(0, 0, H)}|_{\Gamma_{RDT}} = \overline{H} + 2 \frac{i\omega}{c} \overline{B(0, 0, H)}. \tag{36}$$

Next we simplify the right-hand sides of (34). Using the definition of the TRAC boundary condition, we get the following relation for any function $w$:

$$TRAC(\overline{w}) = \overline{TRAC(w)} + 2 \frac{i\omega}{c} \overline{w}|_{\Gamma_{RDT}}. \tag{37}$$

When the last argument of $A_{TRAC}^{-1}$ is zero, we have thus:

$$TRAC(\overline{A_{TRAC}^{-1}(\cdot, \cdot, 0)}) = 2 \frac{i\omega}{c} \overline{B(\cdot, \cdot, 0)}. \tag{38}$$

Using formulas (36) and (38), system (34) reads:

$$\begin{align}
H - \overline{H_{PC}} - 2 \frac{i\omega}{c} \overline{B(0, 0, H_{PC})} &= 2 \frac{i\omega}{c} \overline{B(\overline{f}, G^T, 0)} \tag{39} \\
H_{PC} - \overline{H} - 2 \frac{i\omega}{c} \overline{B(0, 0, H)} &= 2 \frac{i\omega}{c} \overline{B(f, g^T, 0)} \nonumber
\end{align}$$
Note that this system is not linear (w.r.t. complex numbers) due to the conjugations on the unknowns $H$ and $H_{PC}$. This prevents a direct use of Krylov methods. However, by substitution, we write a single equation on $H_{PC}$:

$$ 
\begin{aligned}
& -2 \frac{i\omega}{c} \mathcal{B}(0, 0, H_{PC}) + 2 \frac{i\omega}{c} \overline{\mathcal{B}(0, 0, \overline{H_{PC}})} + 2 \frac{i\omega}{c} \overline{\mathcal{B}(0, 0, H_{PC})} \\
& = 2 \frac{i\omega}{c} \mathcal{B}(\bar{f}, G^T, 0) - 2 \frac{i\omega}{c} \overline{\mathcal{B}(f, g^T, 2 \frac{i\omega}{c} \overline{\mathcal{B}(\bar{f}, G^T, 0)})},
\end{aligned} 
\quad (40) $$

Note that the left-hand side is linear (w.r.t. multiplication by complex numbers) in $H_{PC}$ since we take twice a conjugation. Thus equation 40 can be solved by a Krylov method, see § 7.

An alternative to Krylov methods is the fixed point method. It turns out that it requires to come back to the system on $H$ and $H_{PC}$, namely eq. (34) or equivalently eq. (39). We consider the following iterations:

$$ 
\begin{aligned}
H^{k+1} &= \overline{H_{PC}^k} + 2 \frac{i\omega}{c} \overline{\mathcal{B}(0, 0, H_{PC}^k)} + 2 \frac{i\omega}{c} \overline{\mathcal{B}(\bar{f}, G^T, 0)} = \overline{H_{PC}^k} + 2 \frac{i\omega}{c} \overline{\mathcal{B}(\bar{f}, G^T, H_{PC}^k)} \\
H_{PC}^{k+1} &= \overline{H_{PC}^k} + 2 \frac{i\omega}{c} \overline{\mathcal{B}(0, 0, H_{PC}^k)} + 2 \frac{i\omega}{c} \overline{\mathcal{B}(f, g^T, 0)} = \overline{H_{PC}^k} + 2 \frac{i\omega}{c} \overline{\mathcal{B}(f, g^T, H_{PC}^k)}
\end{aligned} 
\quad (41) $$

By substitution, we write a fixed point method only on the $H_{PC}$ variable:

$$ H_{PC}^{k+1} = H_{PC}^{k-1} - 2 \frac{i\omega}{c} \mathcal{B}(\bar{f}, G^T, H_{PC}^{k-1}) + 2 \frac{i\omega}{c} \overline{\mathcal{B}(f, g^T, \overline{H_{PC}^{k-1}})} + 2 \frac{i\omega}{c} \overline{\mathcal{B}(\bar{f}, G^T, H_{PC}^{k-1})} \quad (42) $$

At convergence, we recover equation (40).

**Remark 2** Although the operator $S$ (33) is non linear, its square $S^2$ is a linear operator. Indeed, we have the following factorization: $S(h) = L\overline{Rh} = L\overline{Rh}$ where $L$ corresponds to apply the TRAC condition and $R$ is the solution operator. Note that both $L$ and $R$ are linear operators that can be expressed as complex valued matrices. Then,

$$ S^2(h) = S \circ S(h) = S(L\overline{Rh}) = L\overline{RLR}h. $$

Thus, the operator $S \circ S$ is a linear operator whose matrix form is $L\overline{RLR}$. Consequently, $H$ is solution to an equation of the form

$$ (I - S \circ S)(H) = (I - L\overline{R}\overline{L}R)(H) = \mathbf{b} $$

where $\mathbf{b}$ is a right-hand side which depends on $f, g^T$ and $G^T$.

Apparently, this expression does not follow a least-squares structure of the type $A^*A$ because we do not define a minimization problem with respect to a given norm.

# 6 Analysis of the Fixed-Point algorithm

From eq. (34), it is easy to check that the convergence of the algorithm amounts to the convergence to zero of the following iterations on the difference $E_{PC}^{k+1}$
between the solution $H_{PC}$ and its approximation:

$$E_{PC}^{k+1} = S(S(E_{PC}^{k-1})). \qquad (43)$$

The convergence of the algorithm relies on the contraction property of the operator $S$ that is proved in Proposition 1.

Recall that the operator $S$ maps any function $g: \Gamma_{RDT} \to \mathbf{C}$ to $\mathcal{T}(g): \Gamma_{RDT} \to \mathbf{C}$ as follows:

Let $u$ be the solution to:

$$ \begin{aligned} & -\frac{\omega^2}{c^2(x)}u - \Delta u = 0 \quad \text{in } \Omega_{top} \\ & \frac{i\omega}{c(x)}u + \frac{\partial u}{\partial n} = 0 \quad \text{on } \Gamma_{SRA} \\ & \mathrm{TRAC}(u) = g \quad \text{on } \Gamma_{RDT} \end{aligned} \qquad (44) $$

and $S(g) := \mathrm{TRAC}(\bar{u})$.

**Proposition 1** For a non vanishing function $g$, we have:

$$ \int_{\Gamma_{RDT}} \frac{c}{\omega} |S(g)|^2 < \int_{\Gamma_{RDT}} \frac{c}{\omega} |g|^2. $$

**Proof.** We proceed as in [12] but in a formal way. We multiply the Helmholtz equation by $\bar{u}$ and integrate by parts:

$$ \int_{\Omega_{top}} -\frac{\omega^2}{c^2}|u|^2 + |\nabla u|^2 + \int_{\Gamma_{SRA}} \frac{i\omega}{c}|u|^2 = \int_{\Gamma_{RDT}} \frac{\partial u}{\partial n}\bar{u} $$

We use the identity:

$$ \frac{1}{4}(|a+b|^2 - |-a+b|^2) = \Re(\bar{a}b) $$

for $a, b \in \mathbf{C}$. We use it on the boundary $\Gamma_{RDT}$ with $a = \frac{i\omega}{c}u$ and $b = \frac{\partial u}{\partial n} - \frac{u}{\beta r}$ and we get:

$$ \frac{1}{4}(|\frac{i\omega}{c}u + \frac{\partial u}{\partial n} - \frac{u}{\beta r}|^2 - |\frac{-i\omega}{c(x)}u + \frac{\partial u}{\partial n} - \frac{u}{\beta r}|^2) = \Re(\frac{-i\omega}{c}\bar{u}(\frac{\partial u}{\partial n} - \frac{u}{\beta r})) = \frac{\omega}{c}\Im(\bar{u}\frac{\partial u}{\partial n}) $$

so that

$$ \frac{c}{4\omega}(|\frac{i\omega}{c}u + \frac{\partial u}{\partial n} - \frac{u}{\beta r}|^2 - |\frac{i\omega}{c(x)}\bar{u} + \frac{\partial \bar{u}}{\partial n} - \frac{\bar{u}}{\beta r}|^2) = \Im(\bar{u}\frac{\partial u}{\partial n}) $$

Taking the imaginary part of the energy estimate and substituting with the identity above:

$$ \int_{\Gamma_{SRA}} \frac{\omega}{c} |u|^2 = \frac{1}{4} \int_{\Gamma_{RDT}} \frac{c}{\omega} (|g|^2 - |\mathcal{T}(g)|^2). \qquad (45) $$

This proves the inequality in the non-strict sense. If there is the equality, then $u$ is zero on $\Gamma_{SRA}$ and thus due to the boundary condition we also have that $\frac{\partial u}{\partial n} = 0$ on $\Gamma_{SRA}$. By the unique continuation property, $u$ is zero everywhere in $\Omega_{top}$ and thus $g$ is zero as well.

We then deduce the uniqueness of the solution to eq. (40). ■
<figure><img src="image_4.png" /><figcaption>Figure 2: Relative errors for a gap $\delta = 1.5\lambda$ (left) and $\delta = 2.\lambda$ (right)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: Computational domain $\Omega$</figcaption></figure>

**Remark 3** At the discrete level with a finite element method, energy estimate (45) will hold. Nevertheless the concluding argument of the proof, the unique continuation property is no more valid.

# 7 Numerical Results

We illustrate our method by comparing it with the Time Reversal method, the TRAC method and our iterative primal-dual TRAC approach. We consider the following configuration, see Figure 3. The domain $\Omega$ is a disk of radius $10\lambda$, where $\lambda$ denotes the wavelength in the background medium characterized by a wave number equal to 1. It contains two scatterers: one of size $\lambda$ located at the center, and another of size $3\lambda/8$ separated by a variable gap $\delta$. The sub-domain $\Omega_{\text{top}}$ is defined as $\Omega \setminus B$, where $B$ is the disk of radius $\lambda$ centered at the origin. The external boundary $\Gamma_{\text{SRA}}$ corresponds to $\partial\Omega$. In the scatterers, the wavenumber is 0.2. The source is an impinging plane wave with changing incidence.

We compare the following methods:

* Time Reversal method (RTM), see [8, 9, 6]
<figure><img src="image_6.png" /><figcaption>Figure 4: True and reconstructed scattered fields without noise for a gap $\delta = 1.5 \lambda$ depicted at the same scale</figcaption></figure>

* original TRAC method see [2]

* the primal-dual TRAC method introduced in § 3.3.2.

The first two methods are one shot methods whereas the new method is iterative for which we present results for: the fixed point method (42) with various iteration counts and for the GMRES solver applied to system (40) with a stopping criterion equal to $\max\{10^{-3}, \sigma \times 0.1\}$ where $\sigma$ is the level of noise, see e.g. [20]. More precisely, we shall add multiplicative Gaussian noise by replacing a recorded data $g$ on $\Gamma_{SRA}$ by

$$g \leftarrow (1 + \sigma * \text{randn}) * g, \quad (46)$$

where $\text{randn}$ satisfies a centered reduced normal law and $\sigma$ is the level of noise taken from 0% to 32%.

For numerical simulations, the Helmholtz equations are approximated by the FreeFEM package, see [21].
<figure><img src="image_7.png" /><figcaption>Figure 5: True and reconstructed scattered fields with a noise level $\sigma = 0.08$ for a gap $\delta = 1.5 \lambda$ depicted at the same scale</figcaption></figure>

To quantitatively compare the accuracy of the different methods, we plot the relative $L^2$-norm of the error between the reconstructed field $Ru^T$ and the exact solution $u^T$ of eq. (1):

$$E_{\text{rel}} = \frac{\|Ru^T - u^T\|_{L^2(\Omega_{\text{top}})}}{\|u^T\|_{L^2(\Omega_{\text{top}})}}.$$

In figure 2, we plot the relative error $E_{\text{rel}}$ as a function of the noise level for different redatuming methods for two gaps $\delta = 1.5\lambda$ (left) and $\delta = 2\lambda$ (right). All TRAC methods outperform the RTM method applied to the redatuming problem. As expected, the iterative methods are more accurate than the one-shot TRAC method which is used as an initialization of the iterative primal-dual TRAC methods. In these tests, the convergence of the fixed point method is quite fast, 6 iterations are sufficient to nearly reach convergence. GMRES algorithm convergence in 4-5 iterations when there is no noise and in 1 iteration when the noise level $\sigma$ is high. As $\sigma$ increases, all methods deteriorate but
remain robust. Here the illumination angle is $\pi$ (see Figure 3) and varying the incidence does not change the hierarchy of the methods nor order of magnitudes of the errors.

In order to better qualify the different methods of redatuming for the same numerical experiments, we present in figures 4 (without noise) and 5 (noise level $\sigma = 0.08$), the real part of the scattered fields $u^S$: top-left the true one, top-right the primal-dual TRAC reconstruction with GMRES as an iterative solver, bottom-left the TRAC method [1] and bottom-right the RTM redatuming based method. The efficiency of the methods is in agreement with the error curves of fig. 2: the diffracted field reconstructed using the primal-dual TRAC method is almost indistinguishable from the exact diffracted field, whereas the TRAC and RTM redatuming based methods yield more significant errors. As expected for phase conjugation methods, they are all robust with respect to noise.

## 8 Conclusion

In this preprint, we have presented the primal-dual TRAC (Time Reversal Absorbing Conditions) method, which is an exact redatuming approach. Our method is built on a convergent iterative process rather than a least-squares formulation. Numerical results and comparisons with other approaches demonstrate the efficiency and convergence of the method when applied to the stationary Helmholtz equation.

More extensive numerical experiments will be added soon.

## References

[1] Franck Assous, Marie Kray, Frederic Nataf, and Eli Turkel. Time reversed absorbing conditions. *Comptes Rendus Mathematiques*, 348(19-20):1063-1067, 2010.

[2] Franck Assous, Marie Kray, Frédéric Nataf, and Eli Turkel. Time-reversed absorbing condition: application to inverse problems. *Inverse Problems*, 27(6):065003, 2011. http://hal.archives-ouvertes.fr/hal-00491912/fr/.

[3] Franck Assous and Frédéric Nataf. Full-waveform redatuming via a trac approach: a first step towards target oriented inverse problem. *Journal of Computational Physics*, 440:110377, 2021.

[4] Alvin Bayliss, Max Gunzburger, and Eli Turkel. Boundary conditions for the numerical solution of elliptic equations in exterior regions. *SIAM J. Appl. Math.*, 42(2):430-451, 1982.

[5] Alvin Bayliss and Eli Turkel. Radiation boundary conditions for wave-like equations. *Comm. Pure Appl. Math.*, 33(6):707-725, 1980.

[6] Edip Baysal, Dan D Kosloff, and John WC Sherwood. Reverse time migration. *Geophysics*, 48(11):1514-1524, 1983.
[7] J.R. Berryhill. Wave-equation datuming. *Geophysics*, 44(206):132944, 1979.

[8] Jon F Claerbout. Toward a unified theory of reflector mapping. *Geophysics*, 36(3):467–481, 1971.

[9] Jon F. Claerbout. *Imaging the Earth's interior*. Blackwell, 1985.

[10] J.C. Curlander and R.C. McDonough. *Synthetic Aperture Radar*. Wiley, 1991.

[11] Maya de Buhan and Marie Kray. A new approach to solve the inverse scattering problem for waves: combining the TRAC and the Adaptive Inversion methods. *Inverse Problems*, 29(8):085009, 2013.

[12] Bruno Després. Domain decomposition method and the helmholtz problem. *Mathematical and Numerical aspects of wave propagation phenomena*, pages 44–52, 1991.

[13] Bjorn Engquist and Andrew Majda. Absorbing boundary conditions for the numerical simulation of waves. *Math. Comp.*, 31(139):629–651, 1977.

[14] Ioannis Epanomeritakis, Volkan Akçelik, Omar Ghattas, and Jacobo Bielak. A Newton-CG method for large-scale three-dimensional elastic full-waveform seismic inversion. *Inverse Problems*, 24(3):034015, 2008.

[15] Mathias Fink. Time-Reversal Acoustics. *Journal of Physics: Conference Series*, 118:012001, 2008.

[16] Marcus J. Grote, Marie Kray, Frédéric Nataf, and Franck Assous. Time-dependent wave splitting and source separation. *Journal of Computational Physics*, 330:981–996, 2017.

[17] Qiang Guo and Tariq Alkhalifah. Datum-based waveform inversion using a subsurface-scattering imaging condition. *Geophysics*, 84(4):S251–S266, 2019.

[18] Qiang Guo and Tariq Alkhalifah. Target-oriented waveform redatum-ing and high-resolution inversion: Role of the overburden. *Geophysics*, 85(6):R525–R536, 2020.

[19] Eldad Haber, Uri M. Ascher, and Doug Oldenburg. On optimization tech-niques for solving nonlinear inverse problems. *Inverse Problems*, 16:1263, 2000.

[20] Martin Hanke. *Conjugate gradient type methods for ill-posed problems*. Chapman and Hall/CRC, 2017.

[21] Frédéric Hecht. New development in FreeFEM++. *Journal of Numerical Mathematics*, 20(3-4):1–14, 2012.
[22] Michael Hinze, Rene Pinnau, Michael Ulbrich, and Stefan Ulbrich. *Optimization with PDE Constraints*, volume 23 of *Mathematical Modelling: Theory and Applications*. Springer Netherlands, first edition, 2009.

[23] Felix Lucka, Nam Huynh, Marta Betcke, Edward Zhang, Paul Beard, Ben Cox, and Simon Arridge. Enhancing compressed sensing 4D photoacoustic tomography by simultaneous motion estimation. *SIAM J. Imaging Sciences*, 11(4):2224-2253, 2018.

[24] Puyan Mojabi and Joe LoVetri. Microwave Biomedical Imaging Using the Multiplicative Regularized Gauss-Newton Inversion. *IEEE Antennas Wireless Propag. Lett.*, 8:645-648, 2009.

[25] Wim A. Mulder. Rigorous redatuming. *Geophysical Journal International*, 161(2):401-415, 2005.

[26] M Ravasi. An overview of Marchenko-based redatuming: past, present, (and future). In *80th EAGE Conference and Exhibition 2018 Workshop Programme*, pages cp-556. European Association of Geoscientists and Engineers, 2018.

[27] Gerard T Schuster and Min Zhou. A theoretical overview of model-based and correlation-based redatuming methods. *Geophysics*, 71(4):SI103-SI110, 2006.

[28] V Shtivelman and A Canning. Datum correction by wave-equation extrapolation. *Geophysics*, 53(10):1311-1322, 1988.

[29] Albert Tarantola. *Inverse Problem Theory And Methods For Model Parameter Estimation*. Society for Industrial and Applied Mathematics, 2005.

[30] Jorn F. M. van Doren, Paul M. J. Van den Hof, Jan D. Jansen, and Okko H. Bosgra. Parameter identification in large-scale models for oil and gas production. In *Proc. 18th IFAC World Congress*, pages 10857-10862, Milano, Italy, 2011.

[31] Kees Wapenaar, Jan Thorbecke, Joost Van Der Neut, Filippo Broggini, Evert Slob, and Roel Snieder. Marchenko imaging. *Geophysics*, 79(3):WA39-WA57, 2014.