# Ocean surface radial velocity imaging in the AT-INSAR Velocity Bunching
Model. A functional approach

Fabricio Perez, Miguel Angel Moreles, Hector Morales
Centro de Investigación en Matemáticas
Jalisco s/n, Valenciana
Guanajuato, GTO 36240,Mexico
email: fabricio.perez<EMAIL_ADDRESS>
Mathematics Deparment
Universidad Autonoma Metropolitana
Ciudad de Mexico,Mexico
email<EMAIL_ADDRESS>
###### Abstract.

This work is concerned with the estimation of radial velocities of sea surface
elevations. The data, is a noisy along-track interferometric synthetic
aperture radar (AT-INSAR) image. We assume the Velocity Bunching Model. This
model relates the complex AT-INSAR image at a point in the image plane, with
the radial velocity of a scatterer point in the sea surface. The relation is
by means of a nonlinear integral operator mapping radial velocities into AT-
INSAR images. Consequently, the estimation of radial velocities, amounts to
the solution of nonlinear integral equations. Our proposal is to solve the
latter by Newton’s methods on function spaces, the optimize then discretize
approach. We show that this continuous version is accurate, and faster than
the classical discretize then optimize version. Also a physical comparison is
carried out with the interferometric velocities.

###### Contents

  1. 1 Introduction
  2. 2 Materials and Methods
    1. 2.1 Problem statement
    2. 2.2 Newton’s method for the nonlinear integral equation
    3. 2.3 A minimization approach
  3. 3 Synthetic data
    1. 3.1 Sea surface
    2. 3.2 The AT-INSAR data
    3. 3.3 Interferometric velocities
  4. 4 Numerical Results
    1. 4.1 Radial Velocity Imaging
    2. 4.2 Computational Efficiency
    3. 4.3 A Physical Comparison
    4. 4.4 Technical information
  5. 5 Conclusions

##  1\. Introduction

In recent decades, imaging of the surface of our planet Earth has increased
with the appearance and improvement of tele-detection systems, such as
Synthetic Aperture Radar (SAR). A SAR system is capable of constructing an
image from the information of electromagnetic waves, which are firstly emitted
by the radar and then backscattered by the observed region. See Moreira et al
(2013) for a review on the subject.

In the case of the ocean surface, Goldstein & Zebker (1987) developed the so
called, airborne along-track interferometric synthetic aperture radar (AT-
INSAR). It has been applied to measure ocean surface currents, ocean surface
waves, etc.

Consequently, the understanding of ocean-radar interaction, is of great
relevance both in theory and in practice. On the theoretical side, of interest
is to derive models of such an interaction. One of such models is the focus of
this work, namely, the AT-INSAR Velocity Bunching Model presented in Bao,
Brüning and Alpers (1997).

This model relates the complex AT-INSAR image at a point in the image plane,
with the radial velocity (line of sight velocity component) of a scatterer
point in the sea surface. In the mathematical jargon, this relation is by
means of an integral operator mapping radial velocities into AT-INSAR images.
If the radial velocity is known in a sea surface region, the AT-INSAR image is
readily obtain by quadrature. This is called the direct problem. The purpose
of this work, is to consider the following inverse (imaging) problem:

Given a noisy AT-INSAR image of an unknown scalar field of sea surface
elevations, estimate the scalar field of radial velocities of the sea surface
elevations.

A solution to this inverse problem is already contained in Goldstein & Zebker
(1987). The radar data is acquired by two antennas, the fore and aft, carried
by a flying platform in the along-track direction at a given velocity. While
the aft antenna transmits the signals, both antennas receive the backscattered
signals. These are processed separately, then combined interferometrically. It
follows that the phase difference caused by the motion of the surface, is
proportional to the interferometric velocity. The latter is an approximation
to the radial velocity, hence, yields a solution to the inverse problem.

The interferometric velocity is used in Hwang et al (2013), for observations
of wave breaking in swell-dominant conditions. For further motivation on the
imaging problem, see references therein.

A purpose of this work on solving the inverse problem, is to compare the
estimated radial velocity field with that of the interferometric velocities.
It serves as a query on the At-INSAR Velocity Bunching model.

As a first study, we consider an AT-INSAR image generated by a a swell sea.
Then, we estimate radial velocities by solving the AT-INSAR Velocity Bunching
integral equation.

The core of the paper is on developing numerical methods for solving the
underlying integral equation. It will become apparent that the latter is
nonlinear and oscillatory, which makes the solution challenging.

A classical approach is to discretize the integral equation and solve the
resulting nonlinear system by Newton’s methods. The so called discretize then
optimize approach.

Alternatively we opt to postpone discretization _until the last minute_ , that
is, the optimize then discretize approach. It is proven to be more efficient,
Stuart (2010), and sometimes necessary, Zuazua (2005).

We develop two modified Newton’s methods on function spaces for solution.
First a nonlinear system, second as a nonlinear least squares problem.
Derivatives are computed in the sense of Fréchet. See for instance Cheney
(2001), for the required Functional Analysis.

Both solutions are mutually cross-validating. For comparison we implement also
a discretize then optimize approach, in the case of minimization. As expected,
the former perform better and on execution time are considerably faster.

A physical comparison is also carried out between the bulk kinetic energy on
the ocean surface area under study, associated to the estimated fiel of radial
velocities and that of the interferometric velocities. The comparison is in
terms of relative errors, again, the latter is outperformed.

We work with synthetic data, we generate a field of surface elevations
following the classical variance spectra to surfaces approach. We obtain a
random 2-D realization of a sea surface. Following Mobley (2016), we develop
our own implementation.

##  2\. Materials and Methods

In this section we pose the imaging problem of interest and the modified
Newton’s methods for solution. We follow a functional approach, that is, we
optimize on normed vector spaces of functions.

More precisely, we shall consider all function spaces as subspaces of
$L^{2}\big{(}(a,b)\big{)}$, the space of square summable complex functions.
For two such functions, $\phi$, $\psi$ the inner product is

(1)
$\big{\langle}\phi\,,\psi\big{\rangle}\quad=\quad\int_{a}^{b}\\!\phi(x)\,\overline{\psi(x)}\,dx$

We shall use freely all well known hilbertian properties of $L^{2}$, see for
instance Cheney (2001).

###  2.1. Problem statement

The point of departure is the AT-INSAR configuration as in Goldstein & Zebker
(1987).

The AT-INSAR image is acquired by two antennas, the fore and aft, carried by a
flying platform in the along-track direction at velocity $V$. The antennas are
separated by a a $2B$ distance. We assume that the system operates in mode 1,
the aft antenna transmits radar signals, and both antennas receive the
backscattered signals. The wavenumber of the incoming electromagnetic wave is
denoted by $k_{r}$.

Let $\mathbf{x}=(x,y)$ be in the reference frame for the sea surface
$z(\mathbf{x})$. $x$ is the coordinate in ground range (cross-track), and $y$
the coordinate in azimuth (along-track).

Let $I(\mathbf{x}_{R})$ be the AT-NSAR image at the position
$\mathbf{x}_{R}=(x_{R},y_{R})$ that is associated with the scatterer
$P_{\mathbf{x}}$ at the point $(\mathbf{x},z(\mathbf{x}))$. The distance from
the median of the two antennas to the point $(\mathbf{x},0)$, is denoted by
$R$. Also, denote by $\tau_{s}$, the scene coherence time.

We assume the AT-INSAR Velocity Bunching Model for a complex AT-INSAR image
$I_{vb}$ as introduced in Bao, Bruning and Alpers (1997).

In this model, the AT-INSAR (single-look) integration time, $T_{0}$ is
regarded small compared to the period of the dominant ocean wave. Hence, the
normalized radar cross section (NRCS) $\sigma(\mathbf{x},t)$ and the radial
velocity $u_{r}(\mathbf{x},t)$ vary little and are approximated by quantities
independent of time, denoted by $\sigma_{0}(\mathbf{x})$, $u_{r}(\mathbf{x})$
respectively. In particular, for the radial velocity a first order
approximation in time iabout $t_{0}=\mathbf{x}/V$ is used,

$u_{r}(\mathbf{x},t)\approx u_{r}(\mathbf{x})+a_{r}(\mathbf{x})(t-t_{0}).$

Here $a_{r}(\mathbf{x})$ is the radial acceleration.

With these considerations, the expression for $I_{vb}(y_{R})\equiv
I_{vb}(\mathbf{x}_{R})$ is,

$\displaystyle I_{vb}(y_{R})$ $\displaystyle=$ $\displaystyle\frac{\pi
T_{0}^{2}\rho_{a}}{2}exp\bigg{[}-\frac{4B^{2}}{V^{2}T_{0}^{2}}\bigg{]}\
\int_{-\infty}^{+\infty}\frac{\sigma_{0}(\mathbf{x})}{\rho_{a}^{\prime}(\mathbf{x})}$
$\displaystyle\times\
exp\bigg{[}-2jk_{r}\frac{B}{V}u_{r}(\mathbf{x})\bigg{]}exp\bigg{[}\frac{4B^{2}\rho_{a}^{2}}{V^{2}\,T_{0}^{2}\,\rho_{a}^{\prime
2}(\mathbf{x})}\bigg{]}$ $\displaystyle\times\
exp\bigg{[}\frac{2jBk_{r}}{R}\bigg{(}\frac{2\rho_{a}^{2}}{\rho_{a}^{\prime
2}(\mathbf{x})}-1\bigg{)}\bigg{(}y_{R}-y-\frac{R}{V}u_{r}(\mathbf{x})\bigg{)}\bigg{]}$
$\displaystyle\times\ exp\bigg{[}-\frac{\pi^{2}}{\rho_{a}^{\prime
2}(\mathbf{x})}\bigg{(}y_{R}-y-\frac{R}{V}u_{r}(\mathbf{x})\bigg{)}^{2}\bigg{]}dy$
$\displaystyle\rho_{a}^{\prime}(\mathbf{x})$ $\displaystyle=$
$\displaystyle\Bigg{\\{}\rho_{a}^{2}+\bigg{[}\frac{\pi}{2}\,\frac{T_{0}R}{V}\,a_{r}(\mathbf{x})\bigg{]}^{2}+\frac{\rho_{a}^{2}T_{0}^{2}}{\tau_{s}^{2}}\Bigg{\\}}^{1/2}$

Here $\rho_{a}^{\prime}(\mathbf{x})$ denotes the degraded single-look
azimuthal resolution, and $\rho_{a}=\lambda_{r}R/(2VT_{0})$ is the full-
bandwidth, single-look azimuthal resolution for stationary targets, where
$\lambda_{r}$ denotes the radar wavelength.

The inverse (imaging) problem of interest is: Given AT-INSAR noisy data D of
an unknown sea surface $z$, and given the radar parameters, estimate the
radial velocities $u_{\text{r}}$ of $z$.

It is assumed that $D(\mathbf{x}_{R})$ is a complex AT-INSAR image
$I_{vb}(\mathbf{x}_{R})$ corrupted by additive noise $\eta$, namely
$D=I_{vb}+\eta$.

Notice that the problem amounts to solving an oscillatory nonlinear integral
equation for $u_{r}(\cdot)\equiv u_{r}(x,\cdot)$, for each fixed $x$ in the
cross-track coordinate in the observation area.

For later reference, let us define

(2) $A=\frac{\pi
T_{0}^{2}\rho_{a}}{2}\text{exp}\bigg{[}-\frac{4B^{2}}{V^{2}T_{0}^{2}}\bigg{]}$

Denoting the integrand by $f_{\text{vb}}$, a scalar complex-valued function,
we have.

(3)
$I_{vb}(y_{R})=A\int_{-\infty}^{+\infty}\\!\\!\\!f_{\text{vb}}(u_{r}(\mathbf{x}),a_{r}(\mathbf{x}),\mathbf{x},y_{R})\,dy$

#### Remark

In our case study, it is found a fortiori, that variations of this integral
operator with respect to the radial acceleration are negligible. Consequently,
the terms involving the latter in the Newton’s methods that follow, are
discarded. It is possible to show this mathematically, but we focus on the
numerical results.

###  2.2. Newton’s method for the nonlinear integral equation

The nonlinear integral equation above, implicitly defines a map between some
function spaces $\mathcal{V}$ and $\mathcal{W}$. Namely,

$\mathcal{V}\to\mathcal{W},\quad u_{r}\mapsto I_{vb}(u_{r}).$

To pose the inversion problem as the solution of a nonlinear integral
equation, consider the residual map

$\mathcal{F}:\mathcal{V}\rightarrow\mathcal{W},\quad\mathcal{F}(u_{r})=D-I_{vb}(u_{r}).$

The problem is to find $u_{r}$ such that

$\mathcal{F}(u_{r})=0.$

Assuming Frechet differentiability, we apply the Newton’s method.

Given the initial guess $u_{r}^{0}\in\mathcal{V}$, solve at each iteration $k$
for the function $h$

(4) $\mathcal{F}^{\prime}(u_{r}^{k})\,h\ =\ -\mathcal{F}(u_{r}^{k})$

and update

(5) $u_{r}^{k+1}=u_{r}^{k}\,+\,h.$

Here $\mathcal{F}^{\prime}$ is the Frechet derivative of $\mathcal{F}$. It
follows that $\mathcal{F}^{\prime}=-I_{vb}^{\prime}$, and

$I_{vb}^{\prime}(u_{r})h=A\int_{-\infty}^{\infty}\bigg{[}\frac{\partial
f_{v}b}{\partial u_{r}}\bigg{]}h(y)\,dy$

where $(\partial f_{vb}/\partial u_{r})$ is the vector calculus derivative of
$f_{vb}$ with respect to $u_{r}$,

(6) $\frac{\partial f_{vb}}{\partial
u_{r}}=\bigg{[}\frac{2\,\pi^{2}R\,C}{V\rho_{a}^{\prime}\,\\!\\!{}^{2}}\ -\
j\,\frac{4\,B\,k_{r}\,\rho_{a}^{2}}{V\rho_{a}^{\prime}\,\\!\\!{}^{2}}\bigg{]}\,f_{vb}$

We remark that the scheme (4),(5) is the continuous (infinite dimensional)
version of the Newton’s method.

We discretize with an appropriate quadrature to obtain a finite dimensional
residual $F:\mathbb{R}^{N_{y}}\rightarrow\mathbb{R}^{2N_{y}}$.

The Newton’s method reads: Given an initial guess
$\vec{u}_{r}\,\\!\\!^{0}\in\mathbb{R}^{N_{y}}$, solve at each iteration $k$:

(7) $\displaystyle F^{\prime}(\vec{u}_{r}\,\\!\\!^{k})\,\vec{h}^{k}\ $
$\displaystyle=$ $\displaystyle\ -F(\vec{u}_{r}\,\\!\\!^{k})$ (8)
$\displaystyle\vec{u}_{r}\,\\!\\!^{k+1}\ $ $\displaystyle=$ $\displaystyle\
\vec{u}_{r}\,\\!\\!^{k}\,+\,\vec{h}^{k}$

The linear system to solve in each iteration is overdetermined and ill
conditioned. Thus the numerical solution is somewhat challenging.

We apply Tikhonov regularization using the SVD decomposition of
$F^{\prime}(\vec{u}_{\text{r}}\,\\!\\!^{k})$. Namely

(9) $\displaystyle\vec{h}^{k}\ $ $\displaystyle\approx$ $\displaystyle\
\vec{h}_{\alpha_{T}}^{k}\quad=\quad-\sum_{i=1}^{N_{y}}\frac{\sigma_{i}}{\sigma_{i}^{2}+\alpha_{T}}\Big{[}\vec{u}_{i}\,\\!\\!\cdot
F(\vec{u}_{r}\,\\!\\!^{k})\Big{]}\vec{v}_{i}$

where $\sigma_{i}\in\mathbb{R}_{\geq 0}$ is the $i$-th largest singular value
of $F^{\prime}(\vec{u}_{r}\,\\!\\!^{k})$, $\vec{u}_{i}\in\mathbb{R}^{2N_{y}}$
is the corresponding $i$-th left singular vector, and
$\vec{v}_{i}\in\mathbb{R}^{N_{y}}$ is the corresponding $i$-th right singular
vector.

The regularization parameter is set to the square of the first (and largest)
singular value of $F^{\prime}(\vec{u}_{r}\,\\!\\!^{k})$:

(10) $\displaystyle\alpha_{T}$ $\displaystyle\leftarrow$
$\displaystyle\sigma_{1}^{2}$

###  2.3. A minimization approach

For comparison, we develop a descent method for the corresponding nonlinear
squares problem. The underlying functional,

$\mathcal{G}:\mathcal{V}\rightarrow[0,\infty),$

is given by

$\mathcal{G}(u_{\text{r}})=\frac{1}{2}\|\mathcal{F}(u_{\text{r}})\|^{2}=\frac{1}{2}\langle\mathcal{F}(u_{\text{r}})\,,\,\mathcal{F}(u_{\text{r}})\rangle$

It is readily seen that the Fréchet derivative of $\mathcal{G}$ at
$u_{\text{r}}$ applied to $h$ is:

$\mathcal{G}^{\prime}(u_{\text{r}})h=\frac{1}{2}\left[\langle
h,(\mathcal{F}^{\prime}(u_{\text{r}}))^{*}\mathcal{F}(u_{\text{r}})\rangle+\overline{\langle
h,(\mathcal{F}^{\prime}(u_{\text{r}}))^{*}\mathcal{F}(u_{\text{r}})\rangle}\right]$

where $(\mathcal{F}^{\prime}(u_{\text{r}}))^{*}$ is the adjoint operator of
$\mathcal{F}^{\prime}(u_{\text{r}})$.

Hence

$\mathcal{G}^{\prime}(u_{\text{r}})h=\langle
h,Re\left\\{(\mathcal{F}^{\prime}(u_{\text{r}}))^{*}\mathcal{F}(u_{\text{r}})\right\\}\rangle$

Substituting $\mathcal{F}^{\prime}(u_{\text{r}})$, it follows that

(11) $\displaystyle\mathcal{G}^{\prime}(u_{r})h\ $ $\displaystyle=$
$\displaystyle\
\Bigg{\langle}h,Re\,\bigg{\\{}\\!-A\int_{-\infty}^{\infty}\bigg{[}\frac{\partial
f_{vb}}{\partial
u_{r}}\bigg{]}\,\overline{\bigg{[}\big{[}\mathcal{F}(u_{r})\big{]}(y_{R}})\bigg{]}\,dy_{R}\,\bigg{\\}}\Bigg{\rangle}$

By the Riesz representation theorem, the second argument of the inner product
is the continuous gradient of $\mathcal{G}$ at $u_{r}$. That is, the function
$\nabla\mathcal{G}(u_{r})$, defined as:

(12) $\displaystyle\nabla\mathcal{G}(u_{r})\ $ $\displaystyle=$ $\displaystyle
Re\,\bigg{\\{}\\!-A\int_{-\infty}^{\infty}\bigg{[}\frac{\partial
f_{vb}}{\partial
u_{r}}\bigg{]}\,\overline{\bigg{[}\big{[}\mathcal{F}(u_{r})\big{]}(y_{R}})\bigg{]}\,dy_{R}\,\bigg{\\}}$

The minimization problem is solved by the BFGS method with known gradient. In
our case, a discrete version of (12).

##  3\. Synthetic data

The area under study is a square $Q$, with $1280m$ side length centered at the
origin. A $128\times 128$ uniform square mesh is considered, with computation
points at the vertices.

The purpose of this section is to specify the parameters to simulate the ocean
surface $z:Q\rightarrow\mathbb{R}$ and its associated AT-INSAR data
$D:Q_{R}\rightarrow\mathbb{C}$.

###  3.1. Sea surface

We follow the classical variance spectra to surfaces approach, to generate a
random 2-D realization of a sea surface. See Mobley (2016).

From knowledge from real ocean surfaces one starts with an omnidirectional
spectrum. For a swell sea we use at frequency $k$ the spectrum,Bao, Bruening
and Alpers (1997).

$\mathcal{S}_{S}(k)=\frac{\alpha_{S}}{2\,k^{3}}\,exp\Bigg{[}-\frac{5}{4}\bigg{(}\frac{k}{k_{S}}\bigg{)}^{-2}\Bigg{]}\,\,\gamma_{S}\,\\!^{G_{S}}$

where

$G_{S}=exp\Bigg{[}\\!\\!-\frac{1}{2}\,\frac{\big{(}k^{1/2}\ -\
k_{S}^{1/2}\big{)}^{2}}{\sigma_{S}^{2}\,k_{S}}\Bigg{]}$

Here,

  * •

$\alpha_{S}=0.212\times 10^{-3}$, the energy scale of $\mathcal{S}_{S}$.

  * •

$k_{S}\ =\ 2\pi\,(\lambda_{S})^{-1}$, the spatial peak frequency of $S_{S}$.

  * •

$\gamma_{S}=10$, the peak enhancement factor of  $\mathcal{S}_{S}$.

  * •

$\sigma_{\text{S}}$, the spectral width centered at $k_{\text{S}}$,

(13) $\sigma_{\text{S}}\quad=\quad\begin{cases}0.07&\quad\text{for}\ \
k\,\leq\,k_{\text{S}}\\\ 0.09&\quad\text{for}\ \
k\,>\,k_{\text{S}}\end{cases}$

Next a spreading function is used. In this case, the two-sided cosine-power
model,

$\Phi_{\text{cp2}}(k,\phi)\ =\
\frac{1}{2}\,N_{p}\,|\text{cos}\,(\phi-\phi_{w})|^{2p},\qquad|\phi-\phi_{w}|\,\leq\,\pi,$

Leading to the two-sided directional swell spectrum,

$\mathcal{\hat{S}}_{2C}(k_{x},k_{y})\ =\
\frac{1}{k}\,\mathcal{S}_{1}(k)\,\Phi_{\text{cp2}}(k,\phi).$

Where $(k,\phi)$  and  $(k_{x},k_{y})$  are the equivalent polar and
cartesian coordinates, respectively.

We obtain a particular instance of the ocean variance spectrum $\hat{z}$. The
ocean surface z is obtained by computing the discrete inverse Fourier
Transform of $\hat{z}$.

We have developed our own software for ocean surface simulation. Let us show
some graphics from omnidirectional swell spectrum to ocean and radial velocity
surfaces.

In the frequency domain, figures 1(a), 1(b), 1(c) depict the omnidirectional
spectrum $\mathcal{S}_{\text{S}}$, the sampled directional spectrum
$\hat{\text{S}}_{\text{2C}}$ and the sampled-variance spectrum
$|\hat{\text{z}}|^{2}$, respectively.

(a) Omnidirectional spectrum $\mathcal{S}_{\text{S}}$: continuous (line) and
sampled (dots).

(b) Sampled directional spectrum $\hat{\text{S}}_{\text{2C}}$.

(c) Sampled-variance spectrum $|\hat{\text{z}}|^{2}$.

Figure 1. Ocean spectra to construct and validate the simulated ocean surface.

The ocean surface z is shown in figure 2, where three properties of such
surface can be visualised: there is a regular pattern of waves whose
directions are very close to the wind direction $\phi_{w}=0$ [rad], the
majority of wavelengths are around $\lambda_{\text{S}}=100$ [m], and a big
amount of measured wave heights are well characterised by $H_{m_{0}}=0.586215$
[m].

Figure 2. The simulated ocean surface z.

Finally, the scalar field of radial velocities in Figure 3

Figure 3. Scalar fields of radial velocities.

###  3.2. The AT-INSAR data

The configuration of the AT-INSAR system is that of Bao, Bruening and Alpers
(1997). For each pointt $\mathbf{x}$ in the grid, $I_{vb}(\mathbf{x})$ is
approximated by quadrature. Then noise $\eta(\mathbf{x})$ is added. More
precisely,

$\eta(\mathbf{x})=2^{-1/2}\big{[}a_{\eta}(\mathbf{x})\,+\,jb_{\eta}(\mathbf{x})\big{]},$

where elements $a_{\eta}(\mathbf{x})$ and $b_{\eta}(\mathbf{x})$ are
independent real Gaussian random variables with mean $0$ and variance
$\sigma_{\eta}^{2}$.

The standard deviation $\sigma_{\eta}$, is taken as:

(14) $\sigma_{\eta}=\left[10^{SNR/20}\right]^{-1},$

where $SNR$ is the signal to noise ratio in [dB]. Below, we report results for
the value $SNR=174$ [dB].

###  3.3. Interferometric velocities

Given the AT-INSAR noisy image $D$, the interefometric phase $\Phi_{ATI}$ is
given by $\Phi_{ATI}=\angle D$. It is proportional to the inetrferometric
velocity $u_{ATI}$, Goldstein & Zebker (1987). Namely

$u_{ATI}(\mathbf{x})=-\frac{\lambda_{r}}{4\pi}\,\frac{V}{B}\,\Phi_{ATI}(\mathbf{x}_{R}).$

In applications, $u_{ATI}$ is used as an approximation of the radial velocity
$u_{r}$. We gauge this approximation in the results that follow.

##  4\. Numerical Results

###  4.1. Radial Velocity Imaging

For each $x$ in the cross-track coordinate we solve both, the Nonlinear System
(NL), and the Functional Minimization (FM) associated to the nonlinear
integral equation. In total there are $128$ problems ordered from left to
right, from $0$ to $127$.

We stress that we follow the optimize then discretize approach. For
comparison, we use the discretize then optimize in the minimization problem.
Results with the latter shall be referred as DFM.

The main difference is the computation of the gradient, as usual, it is
approximated with appropriate finite differences of the discretized (finite
dimensional) functional.

To avoid bias, we start all iterative methods with $u_{r}^{0}\equiv 0$ as
initial guess. Results are remarkable homogeneous for the 128 problems.

First we show the fit for problems $0$ and $64$ in figures 4 and 5,
respectively. The interferometric velocity is included.

Figure 4. TOP: radial velocity $\text{u}_{\text{r}}$ (black), initial point
$\text{u}_{\text{r}}\,\\!\\!^{0}$ (blue), interferometric velocity
$\text{u}_{\text{ATI}}$ (violet), estimated radial velocity
$\text{u}_{\text{r}}\,\\!\\!^{*}$ (red); MIDDLE:
$\text{Re}\big{\\{}\text{D}\big{\\}}$ (black),
$\text{Re}\big{\\{}\,\text{I}_{0}\big{\\}}$ (blue),
$\text{Re}\big{\\{}\,\text{I}_{\text{ATI}}\big{\\}}$ (violet),
$\text{Re}\big{\\{}\,\text{I}_{*}\big{\\}}$ (red); BOTTOM:
$\text{Im}\big{\\{}\text{D}\big{\\}}$ (black),
$\text{Im}\big{\\{}\,\text{I}_{0}\big{\\}}$ (blue),
$\text{Im}\big{\\{}\,\text{I}_{\text{ATI}}\big{\\}}$ (violet),
$\text{Im}\big{\\{}\,\text{I}_{*}\big{\\}}$ (red). Figure 5. TOP: radial
velocity $\text{u}_{\text{r}}$ (black), initial point
$\text{u}_{\text{r}}\,\\!\\!^{0}$ (blue), interferometric velocity
$\text{u}_{\text{ATI}}$ (violet), estimated radial velocity
$\text{u}_{\text{r}}\,\\!\\!^{*}$ (red); MIDDLE:
$\text{Re}\big{\\{}\text{D}\big{\\}}$ (black),
$\text{Re}\big{\\{}\,\text{I}_{0}\big{\\}}$ (blue),
$\text{Re}\big{\\{}\,\text{I}_{\text{ATI}}\big{\\}}$ (violet),
$\text{Re}\big{\\{}\,\text{I}_{*}\big{\\}}$ (red); BOTTOM:
$\text{Im}\big{\\{}\text{D}\big{\\}}$ (black),
$\text{Im}\big{\\{}\,\text{I}_{0}\big{\\}}$ (blue),
$\text{Im}\big{\\{}\,\text{I}_{\text{ATI}}\big{\\}}$ (violet),
$\text{Im}\big{\\{}\,\text{I}_{*}\big{\\}}$ (red).

The corresponding fitting results for problem 64 are shown in figures 6 and 7

Figure 6.  TOP: radial velocity $\text{u}_{\text{r}}$ (black), initial point
$\text{u}_{\text{r}}\,\\!\\!^{0}$ (blue), interferometric velocity
$\text{u}_{\text{ATI}}$ (violet), estimated radial velocity
$\text{u}_{\text{r}}\,\\!\\!^{*}$ (red); MIDDLE:
$\text{Re}\big{\\{}\text{D}\big{\\}}$ (black),
$\text{Re}\big{\\{}\,\text{I}_{0}\big{\\}}$ (blue),
$\text{Re}\big{\\{}\,\text{I}_{\text{ATI}}\big{\\}}$ (violet),
$\text{Re}\big{\\{}\,\text{I}_{*}\big{\\}}$ (red); BOTTOM:
$\text{Im}\big{\\{}\text{D}\big{\\}}$ (black),
$\text{Im}\big{\\{}\,\text{I}_{0}\big{\\}}$ (blue),
$\text{Im}\big{\\{}\,\text{I}_{\text{ATI}}\big{\\}}$ (violet),
$\text{Im}\big{\\{}\,\text{I}_{*}\big{\\}}$ (red). Figure 7.  TOP: radial
velocity $\text{u}_{\text{r}}$ (black), initial point
$\text{u}_{\text{r}}\,\\!\\!^{0}$ (blue), interferometric velocity
$\text{u}_{\text{ATI}}$ (violet), estimated radial velocity
$\text{u}_{\text{r}}\,\\!\\!^{*}$ (red); MIDDLE:
$\text{Re}\big{\\{}\text{D}\big{\\}}$ (black),
$\text{Re}\big{\\{}\,\text{I}_{0}\big{\\}}$ (blue),
$\text{Re}\big{\\{}\,\text{I}_{\text{ATI}}\big{\\}}$ (violet),
$\text{Re}\big{\\{}\,\text{I}_{*}\big{\\}}$ (red); BOTTOM:
$\text{Im}\big{\\{}\text{D}\big{\\}}$ (black),
$\text{Im}\big{\\{}\,\text{I}_{0}\big{\\}}$ (blue),
$\text{Im}\big{\\{}\,\text{I}_{\text{ATI}}\big{\\}}$ (violet),
$\text{Im}\big{\\{}\,\text{I}_{*}\big{\\}}$ (red).

In figure 8 and 9, we show the RMSE of the nonlinear system solution
$u_{r}^{*}$–NL, the Discrete Functional Minimization solution $u_{r}^{*}$–DFM,
and the interferometric velocity solution $u_{ATI}$. Each point in the
horizontal axis corresponds to a fitting problem, 128 in total.

It is apparent that in both cases the functional (infinite dimensional)
approach performs better that the DFM (finite dimensional) solution.

Figure 8. Blue: $\widetilde{\text{u}}_{\text{r}}^{*}$–NL Red:
$\widetilde{\text{u}}_{\text{r}}^{*}$–FM Black:
$\widetilde{\text{u}}_{\text{r}}^{*}$–DFM Violet: $\text{u}_{\text{ATI}}$
Figure 9. Blue: $\widetilde{\text{u}}_{\text{r}}^{*}$–NL Red:
$\widetilde{\text{u}}_{\text{r}}^{*}$–FM Black:
$\widetilde{\text{u}}_{\text{r}}^{*}$–DFM Violet: $\text{u}_{\text{ATI}}$

###  4.2. Computational Efficiency

In practice, imaging problems are computationally expensive. In the discretize
then optimize approach, the approximation of derivatives by finite differences
is costly. Having the exact derivative, and postponing discretization until
_the last minute_ is in general more efficient.

In figures 10 and 11, we show execution times for the 128 inversion problems.

It is noticed that the functional versions are at least three orders of
magnitud faster.

Figure 10. Blue: $\widetilde{\text{u}}_{\text{r}}^{*}$–NL Red:
$\widetilde{\text{u}}_{\text{r}}^{*}$–FM Black:
$\widetilde{\text{u}}_{\text{r}}^{*}$–DFM Violet: $\text{u}_{\text{ATI}}$
Figure 11. Blue: $\widetilde{\text{u}}_{\text{r}}^{*}$–NL Red:
$\widetilde{\text{u}}_{\text{r}}^{*}$–FM Black:
$\widetilde{\text{u}}_{\text{r}}^{*}$–DFM

###  4.3. A Physical Comparison

A first inspection of radial velocity fitting and AT-INSAR inages, figures 4 -
7, the gain with respect to the interferometric velocity $u_{ATI}$ may seem
marginal.

To compare in terms of a physical quantity, we compute the relative error of
associated kinetic energies. Results are shown in Table 1.

Estimated solution |  $\boldsymbol{\text{u}_{\text{r}}^{*}}$-NL |  $u_{r}^{*}$-FM | $u_{ATI}$
---|---|---|---
RE of KE of the estimated solution | 0.0582341 | 0.0219289 | 0.1180122
Table 1. Relative errors of kinetic energies of the estimated solutions.

The relative error of the interferometric velocity is much greater that the
fitted radial velocities. The error of using $u_{ATI}$ is about $11\%$. In
some applications this might be critical.

###  4.4. Technical information

For the functional versions of the Newton’s methods and the surface
simulations, we developed our own in house implementations. The comparative
performance of the methods above was carried out in a computer with the
following specifications:

#### Hardware

  * •

Processor: AMD © A10-5800B with Radeon(tm) HD Graphics, 3.80 GHz

  * •

Physical memory: 8.00 Gb (7.20 Gb usable).

  * •

Round-off unit (machine epsilon):
$\epsilon_{\mathcal{M}}=2.220446049250313\times 10^{-16}$.

#### Software

  * •

System type: 64-bit operating system.

  * •

Operating system: Windows 7 Professional © 2009 Microsoft Corporation with
Service Pack 1.

  * •

Programming language: Anaconda3 5.2.0 with Python 3.6.5 for 64 bits, Qt 5.9.4,
PyQT5 5.9.2.

  * •

IDE: The Scientific PYthon Development EnviRonment (Spyder) ©, version 3.2.8

For classical numerical methods, e.g., SVD decomposition and BFGS, we used the
Python’s routines.

##  5\. Conclusions

Assuming the AT-INSAR-VB model, we have posed the radial velocities imaging
problem, as the solution to a nonlinear integral equation. We have developed
functional (infinite dimensional) versions of Modified Newton’s methods, to
solve this integral equation. Namely, a nonlinear system method coupled with
Tikhonov regularization, and the BFGS method with known gradient for
functional minimization.

For each technique, we have formulated the solution on function spaces, where
the application of the Newton’s method requires the Fréchet derivative of the
objective functions Cheney 2001.

Discrete models and numerical algorithms have been implemented. The numerical
results are satisfactory. The functional approach leads to faster solutions in
comparison with the classical discretize-then-optimize strategy. The fitting
of the estimated radial velocity improves upon that of the $u_{ATI}$. More
over, the comparison of predicted Kinetic Energies, shows that in some
applications, a better fit other than interferometric velocities is required.

This research is manifold, ocean waves modelling, sea surface imaging,
computational methods, etc. On the modelling side, we have consider only the
swell spectrum. It is of interest to consider for instance the JONSWAP and
Pierson-Moskowitz spectra.

Research on methods for sea surface imaging is ongoing. A straightforward
computational continuation of this work, is the use of High Performance
Computing. In our application, an integral equation is solved for each point
in the cross track coordinate. By the ATI-SAR-VB model, each solution is
independent. Consequently, a parallel implementation in a low level computer
language shall lead to even faster solutions.

## References

  * [1] M. Bao, C. Bruening, and W. Alpers. Simulation of Ocean Waves Imaging by an Along-Track Interferometric Synthetic Aperture Radar. IEEE Transactions on Geoscience and Remote Sensing, 35(3):618,631, May 1997.
  * [2] E. W. Cheney. Analysis for Applied Mathematics. Springer-Verlag, New York, 2001.
  * [3] Goldstein, R. M., & Zebker, H. A. (1987). Interferometric radar measurement of ocean surface currents. Nature, 328(6132), 707-709.
  * [4] Hwang, P. A., Toporkov, J. V., Sletten, M. A., & Menk, S. P. (2013). Mapping surface currents and waves with interferometric synthetic aperture radar in coastal waters: Observations of wave breaking in swell-dominant conditions. Journal of physical oceanography, 43(3), 563-582.
  * [5] C. D. Mobley. Modeling Sea Surfaces. A Tutorial on Fourier Transform Techniques. Version 2.0. Sequoia Scientic, Inc., 2016.
  * [6] A. Moreira, P. Prats-Iraola, M. Younis, G. Krieger, I. Hajnsek, and K. P. Papathanassiou. A Tutorial on Synthetic Aperture Radar. IEEE Geoscience and Remote Sensing Magazine, 1(1):6, 43, March 2013.
  * [7] Stuart, A. M. (2010). Inverse problems: a Bayesian perspective. Acta numerica, 19, 451-559.
  * [8] Zuazua, E. (2005). Propagation, observation, control and numerical approximation of waves. SIAM Review, 47(2), 197-243.