# Target detection using underwater acoustic communication links

Lu Shen1, Yuriy Zakharov1, Benjamin Henson1, Nils Morozs1,
Benoît Parrein2, Paul D. Mitchell1

1School of Physics, Engineering and Technology, University of York, U.K
2Nantes Université, Polytech Nantes, CNRS, LS2N, F-44000 Nantes, France

###### Abstract

Underwater monitoring and surveillance systems are essential for underwater
target detection, localization and classification. The aim of this work is to
investigate the possibility of target detection by using data transmission
between communication nodes in an underwater acoustic (UWA) network, i.e, re-
using acoustic communication signals for target detection. A new target
detection method based on estimation of the time-varying channel impulse
response between the communication transmitter(s) and receiver(s) is proposed
and investigated. This is based on a lake experiment and numerical experiments
using a simulator developed for modeling the time-varying UWA channel in the
presence of a moving target. The proposed detection method provides a clear
indication of a target crossing the communication link. A good similarity
between results obtained in the numerical and lake experiments is observed.

###### Index Terms:

distributed networks, lake experiment, target detection, underwater acoustics

##  I Introduction

Underwater monitoring and surveillance systems are essential for underwater
target detection, localization and classification. Dedicated passive and
active sonar systems are commonly used for this purpose. The detection
performance of passive systems is limited when the acoustic sound produced by
the target is of a low intensity, in which case active sonars, monostatic or
bistatic, are preferable. Monostatic sonars are based on backscattering of
acoustic signals from the target. It is however known that the forward
scattering from the target, exploited in bistatic sonars, may have a higher
intensity [1].

Active multiple-input multiple-output (MIMO) sonar systems can significantly
improve target detection performance [2, 3]. When the number of MIMO links is
large, it is possible to detect the target with a single MIMO snapshot [2].

A number of lake experiments were conducted between 1997 and 1999 to
investigate methods of detecting forward scattering from a moving target [4,
5]. In these experiments, tone signals are simultaneously transmitted at a set
of carrier frequencies within the frequency range 1–3 kHz. A signal source
fixed at the lake bottom or a vertical antenna array is used for transmission.
Both vertical and horizontal receiving arrays are deployed at a distance of
about 400 m. Targets of several meters in length and up to 1 m height are
towed underwater across the area between the transmitter and receiving arrays.
The aim is to detect the diffracted signals against the fluctuating direct
signal. In [4], fluctuations of the direct signal are removed by a high-pass
filter, the cut-off frequency of which is computed based on a priori
information, including the size and velocity of target. The signal then goes
through a filter matched to a theoretical model of the diffracted signal and
incoherently accumulated over the array elements and the carrier frequencies.
The target detection is performed based on the variation of the output signal.
The performance is evaluated by the output signal-to-noise ratio
($\mathrm{SNR}_{\mathrm{out}}$) which is defined as the ratio between an
increment in the output signal and the standard deviation of the output signal
without a target. It is found in [5] that $\mathrm{SNR}_{\mathrm{out}}$ can be
improved by jointly processing the output signals of two vertical arrays.
Another finding is that it is possible to indicate the direction of the target
movement based on the shape of the peaks in the processed signal. A
comparative analysis of detection methods is conducted in [6] using the
experimental data described in [4, 5]. It is found that the use of a vertical
array shows more robust target detection performance compared to a horizontal
array. Further analysis is conducted in [7] to investigate the detection
performance with coherent, partially coherent and incoherent spatial
processing. It is indicated that the choice of spatial processing method
depends on the vertical distribution of the acoustic field. Aforementioned
works focus on investigating matched filtering based methods for detecting the
forward scattering from a target using transmit and/or receive antenna array
and low-frequency pulse transmission.

In [8], an adaptive cancellation algorithm is applied to the received signals
to enhance the detection of forward scattering from the strong direct path in
a lake experiment. The source is located at 4 m depth. A vertical array with
10 hydrophones is deployed at a distance of 900 m. The dimension of the target
is 6 m$\times$1 m. Successive pulses at a 10 kHz carrier frequency are
transmitted. An adaptive filter is used for signal processing. The input of
the adaptive filter is a segment of a signal received when the target is
absent. Other segments of the received signal are used as the desired signal,
with which the filter output is compared to produce the adaptive filter error.
As the forward scattered signal is weakly correlated with the direct signal,
it is assumed that the algorithm will only cancel the direct signal, and the
adaptive filter error will be higher when the target is present. This method
also requires a vertical antenna array.

In this work, we consider scenarios with a distributed network of multiple
unsynchronized acoustic sensor nodes capable of communicating with each other,
instead of using source/receiver arrays. Such a network provides flexibility
in the deployment of the acoustic sensor nodes based on the environment.
Wideband acoustic signals at a relatively high frequency (32 kHz carrier
frequency in our experiments) are used for data transmission. The aim of this
work is to investigate the possibility of re-using the communication signals
in the network for target detection. The target detection is performed based
on the measurement of the time-varying channel impulse responses (CIRs)
between communication nodes in the network. Note that the channel estimation,
i.e. CIR measurement, is one of the main steps of signal processing in the
receiver, hence an CIR estimate is normally available at the receiver. The
idea is to detect changes in the CIR caused by the target movement. The use of
wideband signal provides high resolution CIRs, which increases the chance of
target detection. In this paper, we focus on the target detection performed
for each individual network link. Designs of multiple access control and
network layer protocols for such a target detection network are presented in
[9, 10].

Figure 1: Acoustic propagation between two nodes (A and B) with presence of a
target (T) in a distributed network (Top view).

Fig. 1 illustrates a scene of a distributed network with six Tx/Rx nodes and a
target. We consider the scenario when there is a ‘silent’ target moving
underwater. When the target is absent, the acoustic multipath propagation from
node A to node B only includes the node-to-node link. When there is a target,
apart from the link between two nodes, scattering from the target will also be
present. The acoustic propagation between two nodes in the presence of the
target will correspond to a different CIR compared to the ‘no target’ case.
Therefore, it is possible to perform target detection based on changes in the
CIR in time. A challenge is to distinguish between changes in the CIR due to
the environmental factor, such as node drift or surface waves, and change
caused by a moving target. To this end, a detection algorithm is proposed in
this paper which reduces the impact of the environmental variation on the CIR
and subsequent detection performance.

For performance evaluation of the target detection scheme, an acoustic
simulator for conducting virtual signal transmission experiments is developed;
it is capable of modelling the multipath propagation and Doppler effects for
different UWA environments taking into account the motion of the
transmitter/receiver/target. One of the key requirements of the simulator is
to generate highly accurate acoustic field information for the communication
links based on the trajectories of the nodes and the target. Unlike the
typical static/slowly moving transmitter and receiver nodes, the target might
be moving at a high speed. To track the fast variation of the acoustic
channel, a high sampling rate of channel is required. Therefore, the simulator
computes the acoustic field at the baseband sampling interval of the
communication signals. To reduce the computational time, the simulator
exploits the principles of the Waymark and VirTEX simulators [11, 12], where
the acoustic field is pre-computed within an area of interest on a grid of
space points and then interpolated based on the pre-computed information.

The field interpolation method proposed in [12] assumes that the receiver is
in the far field with respect to the signal source. However, the target can
reflect/re-radiate a transmission and could be positioned closely to the
communication nodes, thus the interpolation method used in [11, 12] might not
be accurate enough. To improve the field interpolation in the near field, an
improved interpolation approach is proposed. The proposed simulator is used to
generate the time-varying CIR and the signal at the receiver.

To verify the target detection performance and validate the simulator, a lake
experiment is conducted with several nodes and a target whose trajectory is
accurately estimated using acoustic navigation. In the lake experiments, the
CIR estimates are obtained from the channel estimator in the receiver and then
used for the target detection. After the lake experiment, we run the simulator
with the same node configuration, and use the target trajectory estimate as
input to generate the CIRs and received signals. A good match between the
experimental and simulation results is observed.

The contributions of this paper are as follows.

  * •

A target detection algorithm is proposed based on measuring changes in the
time-varying CIRs between communication nodes, which are estimated using the
existing communication signals in the network.

  * •

An acoustic simulator is developed, which takes into account the motion of the
target and the communication nodes. The output of the simulator includes the
signals received at nodes and time-varying CIRs between each transmitter-
receiver pair. An approach is proposed and implemented to improve the acoustic
field interpolation in the near field. The simulator is validated using the
lake experiment.

  * •

The performance of the proposed detection method is investigated using
numerical and lake experiments.

The structure of the paper is as follows. Section II describes the target
detection algorithm. Section III introduces the acoustic simulator. Section IV
describes the acoustic field interpolation approaches. Numerical simulation
and lake experiments are presented in Section V and VI, respectively. The
paper is concluded in Section VII.

##  II Target detection based on CIR variation

The key idea of target detection is to identify changes in the CIR caused by a
moving target. The changes in the CIR can be measured by the mean squared
deviation (MSD) [13]:

$\mathrm{MSD}=\sum_{n=0}^{L-1}|h_{0}(n)-h(n)|^{2},$ (1)

where $L$ is the length of the CIR, $h_{0}(n)$ is the reference CIR when the
target is absent, $h(n)$ is the CIR of interest, obtained at a time interval
where a target might be present. The MSD in (1) can be computed using the
frequency-domain representation of the channel:

$\mathrm{MSD}=\sum_{k=0}^{K-1}|H_{0}(k)-H(k)|^{2},$ (2)

where $H_{0},H$ are the frequency responses of filters $h_{0}$ and $h$,
respectively, and they are computed using the fast Fourier Transform (FFT) of
a size $K$. When the MSD is a small value, we can assume that the target is
absent. For high MSD values, we can assume that the target is present.

However, even when the target is absent, the CIR will be time-varying due to
the time-varying nature of the UWA channel, thus resulting in a high MSD.
Therefore, we need to distinguish between the change of the impulse response
due to the time-varying environment and the change caused by a moving target.
Instead of using $h_{0}$ in (1), we can use $h_{0}*g$, where $*$ denotes the
convolution operation and $g$ is an impulse response representing a filter
with a small delay spread. In such a case, the change in delay and amplitude
caused by the slow motion of the source and the receiver, the variation of the
environment and the clock mismatch between the transmitter and the receiver
can be eliminated from influencing the MSD. When a target is crossing the
communication link, the multipath structure of the propagation channel will
change significantly, and these changes cannot be compensated by the filter
$g$.

Small delays introduced by the filter $g$ in the time domain correspond to
complex exponentials of long periods in the frequency domain. Therefore, the
frequency response of $g$ can be expressed as a combination of complex
exponentials with expansion coefficients $c(p)$ [14]:

$G(k)=\sum_{p=0}^{2P}c(p)e^{j2\pi(-P+p)k/K},\,\,k=0,\ldots,K-1,$ (3)

where the number of complex exponentials $2P+1$ corresponds to the delay
spread of the filter $g$.

Therefore, to find the filter $g$, we solve the following optimization
problem:

$\hat{G}=\arg\\!\min_{G}\sum_{k=0}^{K-1}|H_{0}(k)G(k)-H(k)|^{2}.$ (4)

We denote the combined frequency response as:

$\begin{split}\tilde{H}(k)&=H_{0}(k)G(k)\\\
&=\sum_{p=0}^{2P}c(p)\tilde{G}(k,p),\end{split}$ (5)

where $\tilde{G}(k,p)=H_{0}(k)e^{j2\pi(-P+p)k/K}$, and find the expansion
coefficients $\mathbf{c}=[c(0),\ldots,c(2P)]^{T}$ by minimizing the cost
function:

$J(\mathbf{c})=\sum_{k=0}^{K-1}|\tilde{H}(k)-H(k)|^{2}.$ (6)

This is achieved by solving the normal equation:

$(\mathbf{R}+\varepsilon\mathbf{I})\mathbf{c}=\mathbf{b},$ (7)

where $\mathbf{I}$ is an identity matrix,
$\mathbf{R}=\mathbf{\tilde{G}}^{H}\mathbf{\tilde{G}}$,
$\mathbf{b}=\mathbf{\tilde{G}}^{H}\mathbf{H}$, $\mathbf{\tilde{G}}$ is a
$K\times(2P+1)$ matrix with elements $\tilde{G}(k,p)$,
$\mathbf{H}=[H(0),\ldots,H(K-1)]^{T}$ and $\varepsilon>0$ is a regularization
parameter. The solution of (7) is given by:

$\hat{\mathbf{c}}=(\mathbf{R}+\varepsilon\mathbf{I})^{-1}\mathbf{b}.$ (8)

After obtaining $\hat{\mathbf{c}}$, we have
$\tilde{H}(k)=\sum_{p=0}^{2P}\hat{c}(p)\tilde{G}(k,p)$. The deviation in the
frequency response is computed as: $\Delta H(k)=\tilde{H}(k)-H(k)$. Finally,
the normalized MSD is computed by:

$\mathrm{MSD}_{\mathrm{norm}}=\dfrac{\sum_{k=0}^{K}|\Delta
H(k)|^{2}}{\sum_{k=0}^{K}|H_{0}(k)|^{2}}.$ (9)

The normalized MSD is used as the statistic for the target detection. When it
is small, it means that the communication link is not distorted, i.e, the
target is absent. When it is high, it can be inferred that the link distortion
is caused by a target.

##  III Simulator for acoustic links

Underwater acoustic channel simulators, such as Waymark [11, 15, 16], can
model the signal transmission with a moving transmitter and receiver. To
reduce the computational complexity, the acoustic field is sampled at an
interval, which is typically significantly longer than the baseband sampling
interval [11, 16]. In this work, we introduce a moving target crossing the
communication links, and to guarantee accurate modelling of the signal
transmission in the presence of a fast-moving target, the acoustic field is
computed at the baseband sampling rate. Even at baseband frequencies,
generating the acoustic field for prolonged scenarios can be impractical for
meaningful numerical experiments. Therefore, the simulator exploits the
principles used in [11, 12] and performs the task in two stages: pre-
processing and received signal computation.

(a)

(b)

Figure 2: General structure of the simulator: (a) Received signal computation;
(b) Filtering for the $l$th link, $l=1,2,3$, $i$ is the index of the baseband
sample, where $s_{l}$ is the baseband input signal and $r_{l}$ is the baseband
output.

###  III-A Pre-processing

The acoustic field between a pre-defined source position and possible receiver
positions within an area of interest on a predefined set of grid points is
pre-computed. The BELLHOP3D ray-tracing program [17] is used to generate the
acoustic field information, including the delay $\tau$, amplitude $A$,
departure angle $\theta_{d}$ and incident angle $\theta_{i}$ of each arrival
for each grid point. We consider a simplified environment with a flat surface,
flat bottom and range-invariant sound speed profile. The pre-computation is
performed for a set of source depths required to cover the range of possible
source depths. For each source depth, acoustic field computation is performed
on a two-dimentional (depth, range) grid map of receiver positions. For the
source parameters, the anticipated source positions and target trajectories
define the minimum depth, the maximum depth and the resolution (step size) in
depth $\Delta D_{s}$. The parameters of the grid map include the minimum and
maximum range, the minimum and maximum depth, and the resolution of the grid
map in range and in depth.

###  III-B Received signal computation

To obtain the received signal, the acoustic propagation of three links need to
be considered as illustrated in Fig. 2(a). Link 1 is the direct link between
the source and the receiver when the target is absent; Link 2 is the link
between the source and target; and Link 3 is the link between the target and
receiver. The target plays the role of a receiver when modelling link 2, and
it acts as a source when modelling link 3. A target filter is used to
incorporate a target strength model [18].

The acoustic field of all three links are computed by acoustic field
interpolation. The trajectories of the source, receiver and target are
interpolated to the baseband sampling rate. The field interpolation is
performed according to the source, target and receiver positions with respect
to the closest grid points within the grid maps using the field information at
these grid points. More details on the acoustic field interpolation are given
in Section IV.

As shown in Fig. 2(b), the output of the $l$th link is the signal:

$r_{l}(i)=\sum_{m=1}^{M}\tilde{A}_{m}(i)s_{l}(i-\tau_{m}(i)),$ (10)

where $i$ is the index of a baseband sample, $s_{l}(i)$ is the baseband input
of the $l$th link, $M$ is the number of arrivals, $\tau_{m}(i)$ is the delay
of the $m$th arrival obtained from the field interpolation, and
$\tilde{A}_{m}(i)$ is the baseband equivalent amplitude of the $m$th arrival
computed according to beampatterns of the source and receiver:

$\tilde{A}_{m}(i)=G_{s}(\theta_{d}^{m})G_{r}(\theta_{i}^{m})A_{m}(i)e^{-j\omega_{c}\tau_{m}(i)},$
(11)

where $\omega_{c}=2\pi f_{c}$, $f_{c}$ is the carrier frequency of the
transmitted signal, $A_{m}(i)$ is the amplitude of the $m$th arrival obtained
from the field interpolation, $G_{s}$ and $G_{r}$ are the source and receiver
beampatterns, $\theta_{d}^{m}$ is the departure angle of the $m$th arrival,
$\theta_{i}^{m}$ is the incident angle of the $m$th arrival, obtained from the
field interpolation. The baseband input sample $s_{l}(i-\tau_{m}(i))$ is
obtained using the local spline interpolation as described in [15]. The signal
at the receiver is then given by:

$y(i)=r_{1}(i)+r_{3}(i).$ (12)

In this paper, the target is modelled as a rigid sphere of radius $a$. The
bistatic target strength can be written as a function of the angle $\alpha(i)$
between the incident and scattered rays. Here we compute $\alpha(i)$ as:

$\alpha(i)=|\alpha_{\mathrm{ST}}(i)-\alpha_{\mathrm{TR}}(i)|,$ (13)

where $\alpha_{\mathrm{ST}}(i)$ is the azimuth angle between the source and
the target and $\alpha_{\mathrm{TR}}(i)$ is the azimuth angle between the
target and receiver and $\alpha_{\mathrm{TR}}(i)\in(-\pi,\pi]$. Therefore, the
target filter is modelled as a single-tap filter with a target strength
coefficient [19]:

$\sigma_{TS}(i)=\begin{cases}\dfrac{a^{2}}{4}\left[1+\tan^{2}\left(\dfrac{\alpha(i)}{2}\right)J_{1}^{2}\left(wa\sin\alpha(i)\right)\right]&\alpha(i)\neq\pi\\\
\dfrac{a^{2}}{4}(1+w^{2}a^{2})&\alpha(i)=\pi\end{cases}$ (14)

where $w=2\pi/\lambda$ is the wavenumber, $\lambda$ is the acoustic
wavelength, $J_{1}(\cdot)$ is the Bessel function of the first order.

###  III-C Channel impulse response computation

Apart from the signal at the receiver, the simulator also provides the time-
varying CIR between the source and receiver. For every position of the signal
source, receiver and target, the channel frequency response is computed as:

$H(i,k)=H_{1}(i,k)+\sigma_{TS}(i)H_{2}(i,k)H_{3}(i,k),$ (15)

where $H_{1}$, $H_{2}$, $H_{3}$ represent the baseband frequency responses of
the three links in Fig. 2(a). For instance, the frequency response
$H_{1}(i,k)$ is computed as:

$H_{1}(i,k)=\sum_{m=1}^{M}\tilde{A}_{m}(i)e^{-j2\pi k\Delta
f\tau_{m}(i)},\,\,\,k\in[-B/2,B/2],$ (16)

where $B\Delta f$ is the frequency bandwidth, $\Delta f$ is the frequency
resolution. The channel impulse response $\mathbf{h}(i)$ is then computed
using the inverse FFT of $H$.

##  IV Acoustic field interpolation

For pre-computation, a set of grid maps are generated, each corresponding to a
particular source depth. When the source/receiver is not on a grid point,
acoustic field interpolation is performed. Plane wave acoustic field
approximation works well when the receiver is in the far field of the signal
source. In simulation scenarios, the target (as a receiver and signal source)
can be arbitrarliy close to a node; in such cases, the spherical wave
approximation improves the interpolation [20]. In this section, we present
acoustic field interpolation based on plane wave and spherical wave
approximations.

###  IV-A Source at a depth with precomputed grid map and receiver between
grid points

We first consider a scenario when the source is at a depth for which a grid
map has been precomputed and the receiver is located at a point $(x,y)$
between grid points on the map.

####  IV-A1 Plane wave interpolation

Figure 3: An illustration of plane wave interpolation for a receiver position
$(x,y)$ based on the ray information at the grid point $(x_{1},y_{1})$, where
$\theta_{i}$ is the incident angle of a ray arriving at $(x_{1},y_{1})$.
Figure 4: An illustration of spherical wave interpolation for a receiver
position $(x,y)$ based on the ray information at the grid point
$(x_{1},y_{1})$, where $\theta_{i}$ is the incident angle of a ray arriving at
$(x_{1},y_{1})$, $\gamma_{1,1}$ is the angle between $d_{1}$ and $d_{2}$, and
$\beta$ is the angle bewteen $d_{2}$ and $x$ positive axis.

In [12], a plane wave assumption is used for acoustic field interpolation as
shown in Fig. 3. The arrivals at the receiver position $(x,y)$ is the
combination of the arrivals at the four neighbouring grid points
$(x_{a},y_{b}),a,b\in[1,2]$, with weighted amplitudes and delay adjustments.
The delays are adjusted by the difference in the wave travel time between the
neighbouring points and the position of the receiver. Assuming that a ray is
arriving at a grid point $(x_{a},y_{b})$ with an incident angle $\theta_{i}$,
the difference in the travel time is given by [12]:

$\Delta\tau=(\Delta x\cos{(\theta_{i})}+\Delta y\sin{(\theta_{i})})/c,$ (17)

where $c$ is the speed of sound, $\Delta x=x-x_{a}$ and $\Delta y=y-y_{b}$.
The delay of arrival for the $m$th ray from $(x_{a},y_{b})$ is updated as:
$\tau_{m}^{{}^{\prime}}=\tau_{m}+\Delta\tau$.

The amplitude of the arrivals for the four neighbouring grid points are
weighted by:

$\displaystyle a_{1}:\,(1-w_{1})(1-w_{2}),$ (18) $\displaystyle
a_{2}:\,(1-w_{1})w_{2},$ $\displaystyle a_{3}:\,w_{1}w_{2},$ $\displaystyle
a_{4}:\,w_{1}(1-w_{2}),$

where $a_{1}$, $a_{2}$, $a_{3}$, $a_{4}$ represent amplitudes of the arrivals
at the four grid points. The weights are computed based on the position of the
receiver $(x,y)$ with respect to the grid points:

$\displaystyle w_{1}$ $\displaystyle=(x-x_{a})/(x_{2}-x_{1}),$ (19)
$\displaystyle w_{2}$ $\displaystyle=(y-y_{b})/(y_{2}-y_{1}).$

The main steps of the plain wave interpolation can be summarized as:

  1. (i)

Adjust the delay of the arrivals at four grid points based on (17).

  2. (ii)

Apply weighting to the amplitudes of the arrivals at four grid points
according to (18) and (19).

(a)

(b)

Figure 5: Illustrations of interpolation schemes for new source depth based on
the ray information at the known source depth; (a) plane wave interpolation,
$\theta_{d}$ is the departure angle of a ray travelling from the original
source; (b) spherical wave interpolation, $\tilde{\theta}_{i}$ is an incident
angle due to the change in source depth.

####  IV-A2 Spherical wave interpolation

An illustration of delay adjustment when using spherical wave approximation is
shown in Fig. 4. As the delay adjustment for a direct path is straightforward,
we use a surface reflection as an example. The same procedure applies for the
bottom reflections.

We first find the image source position ($\mathrm{Source}^{{}^{\prime}}$)
using the image source method (ISM) [21]. As shown in Fig. 4, $d_{1}$ is the
distance a ray travels to arrive at $(x_{1},y_{1})$ and $d_{3}$ is the
distance it needs to travel to arrive at $(x,y)$. The delay adjustment is made
based on the difference between $d_{1}$ and $d_{3}$.

For the $m$th arrival at a grid point $(x_{a},y_{b})$, the distance between
the image source position and the grid point can be computed as:
$d_{1}=c\tau_{m}$, where $\tau_{m}$ is the delay of the $m$th arrival. The
distance between the receiver position $(x,y)$ and $(x_{a},y_{b})$ can be
found as: $d_{2}=\sqrt{(\Delta x)^{2}+(\Delta y)^{2}}$. By denoting the angle
between $d_{1}$ and $d_{2}$ as $\gamma_{a,b}$, and assuming the angle to be
known, distance $d_{3}$ can be computed based on the cosine rule as:

$d_{3}=\sqrt{d_{1}^{2}+d_{2}^{2}-2d_{1}d_{2}\cos{(\gamma_{a,b})}}.$ (20)

Thus, we can obtain the change in the delay $\tau_{m}$ as:

$\Delta\tau=(d_{3}-d_{1})/c,$ (21)

and the delay of the $m$th arrival at $(x,y)$ can be updated as:
$\tau_{m}^{{}^{\prime}}=\tau_{m}+\Delta\tau$.

The angle $\gamma_{a,b}$ can be computed using the incident angle of the
arrival $\theta_{i}$ and the angle between $d_{2}$ and $x$ positive axis
(denoted as $\beta$). Note that the incident angle $\theta_{i}$ is within
$[-\pi/2,\pi/2]$. For a grid point $(x_{a},y_{b})$:

$\beta=\arctan{\left(\dfrac{|y-y_{b}|}{|x-x_{a}|}\right)}.$ (22)

For the four grid points, we have:

$\displaystyle\gamma_{1,1}$
$\displaystyle=\begin{cases}\pi-\mathrm{sign}(\theta_{i}+\beta)(\theta_{i}+\beta)&\theta_{i}<0\\\
\pi-(\theta_{i}+\beta)&\theta_{i}\geq 0\end{cases}$ (23)
$\displaystyle\gamma_{1,2}$
$\displaystyle=\begin{cases}\pi-\mathrm{sign}(\theta_{i}-\beta)(\theta_{i}-\beta)&\theta_{i}>0\\\
\pi+\theta_{i}-\beta&\theta_{i}\leq 0\end{cases}$ $\displaystyle\gamma_{2,1}$
$\displaystyle=\begin{cases}\mathrm{sign}(\theta_{i}-\beta)(\theta_{i}-\beta)&\qquad\theta_{i}>0\\\
\theta_{i}-\beta&\qquad\theta_{i}\leq 0\end{cases}$
$\displaystyle\gamma_{2,2}$
$\displaystyle=\begin{cases}\mathrm{sign}(\theta_{i}+\beta)(\theta_{i}+\beta)&\qquad\theta_{i}<0\\\
\theta_{i}+\beta&\qquad\theta_{i}\geq 0\end{cases}$

The amplitudes of the arrivals are weighted according to (18) and (19).

###  IV-B Source at a depth without precomputed grid map and receiver at a
grid point

Now we consider the scenario when a source depth deviates from the source
depth, for which a grid map exists, and the receiver position is on a grid
point of the map. The idea is to adjust the delay of each arrival based on the
deviation of the source depth from the nearest source depth with a precomputed
grid map.

####  IV-B1 Plane wave interpolation

As illustrated in Fig. 5(a), the difference in the travel time can be computed
as:

$\Delta\tau=-\Delta D\sin{(\theta_{d})}/c,$ (24)

where $\Delta D=D_{\mathrm{off}}-D_{\mathrm{on}}$, $D_{\mathrm{off}}$ is a
source depth not with a pre-computed grid map and $D_{\mathrm{on}}$ is the
nearest source depth with a grid map.

####  IV-B2 Spherical wave interpolation

We first find the image receiver position based on ISM as shown in Fig. 5(b).
Assuming that there is a ray from the depth $D_{\mathrm{on}}$ with a departure
angle of $\theta_{d}$ arriving at the receiver with a delay $\tau_{m}$, the
distance $d_{1}$ between the source and the image receiver can be computed as:
$d_{1}=\tau_{m}c$. Based on the cosine rule, we have:

$d_{2}=\begin{cases}\sqrt{d_{1}^{2}+(\Delta D)^{2}-2d_{1}|\Delta
D|\cos{(\pi/2-\theta_{i})}}&\theta_{d}>0,\\\ \sqrt{d_{1}^{2}+(\Delta
D)^{2}-2d_{1}|\Delta D|\cos{(\pi/2+\theta_{i})}}&\theta_{d}\leq 0.\end{cases}$
(25)

The delay can be computed as:

$\tau_{m}^{{}^{\prime}}=\tau_{m}+\Delta\tau=\tau_{m}+(d_{2}-d_{1})/c.$ (26)

For the plane wave approximation, the incident angle remains the same when
there is a deviation in source position. However, for spherical wave
approximation, there will be a change in the incident angle as illustrated in
Fig. 5(b). The new incident angle $\tilde{\theta}_{i}$ is computed based on
the angle $\gamma_{i}$ between $d_{1}$ and $d_{2}$. The angle $\gamma_{i}$ is
given by:

$\gamma_{i}=\arccos{\left(\dfrac{d_{1}^{2}+d_{2}^{2}-|\Delta
D|^{2}}{2d_{1}d_{2}}\right)}.$ (27)

Finally, $\tilde{\theta}_{i}$ is computed as:

$\tilde{\theta}_{i}=\theta_{i}+\mathrm{sign}(\theta_{i})\mathrm{sign}(\Delta\tau)\gamma_{i}.$
(28)

###  IV-C Arbitary position of source and receiver

Now we consider the general scenario when both the source and the receiver
positions deviate from a known depth/grid point.

The main steps of interpolation are:

  1. 1.

Adjust delays of the arrivals at the four grid points based on the deviation
of the source depth. For spherical interpolation, use (25) and (26). For plane
wave interpolation, $\Delta\tau$ is computed using (24).

  2. 2.

For the spherical interpolation only, adjust the incident angles of the
arrivals at four grid points based on the deviation of source depth using (27)
and (28).

  3. 3.

Adjust delays of the arrivals based on the receiver position between grid
points. For spherical interpolation, use (20) - (23). For plane wave
interpolation, use (17).

In subsection IV-A and IV-B, we consider specific circumstances when the
source/receiver is at a known depth/grid point to breakdown the description of
the calculations for arbitary source/receiver positions. In practice, for any
source and receiver positions, we follow the interpolation steps summarized in
this subsection.

##  V Numerical simulation

In this section, we investigate the performance of the acoustic field
interpolation and show an example of the time-varying CIR when a target
crosses a link between a signal source and receiver.

###  V-A Acoustic field interpolation

In this subsection, the performance of the plane wave and spherical wave
acoustic field interpolation methods are compared.

We consider an environment with a flat surface and a flat bottom. The sea
depth is 200 m, and the sound speed profile is the same as in the SWellEx-96
experiment [22]. Six grid maps are pre-computed for source depths from 15 m to
20 m with a depth step size of $\Delta D_{s}=1$ m. The signal source depth in
the simulation trials is random and independently uniformly distributed within
$[15,20]$ m. Each grid map covers the receiver range from 1 m to 50 m and
receiver depths from 5 m to 40 m. The resolution of the grid maps in both
range and depth is 0.5 m. In every simulation trial, the receiver range and
depth are independent and uniformly distributed within $[15,20]$ m and
$[5,40]$ m, respectively.

The interpolation accuracy is evaluated by computing the mean squared
deviation:

$\mathrm{MSD}_{\mathrm{interp}}=10\log_{10}\left(\dfrac{\sum_{n=0}^{K-1}\left(|\hat{h}(n)|-|h(n)|\right)^{2}}{\sum_{n=0}^{K-1}|h(n)|^{2}}\right),$

where $\mathbf{\hat{h}}=[\hat{h}(0),\ldots,\hat{h}(K-1)]^{T}$ is a CIR
estimate and $\mathbf{h}=[h(0),\ldots,h(K-1)]^{T}$ is the true CIR.

The CIR estimate and the true CIR are computed as the inverse FFT of the
channel frequency response with the interpolated and true field information,
respectively. The interpolated channel frequency response is computed using
(15) and (16). For computation of the true CIR, the amplitude and delay of
arrivals are generated at the precise positions of the source and receiver by
running the BELLHOP3D program. For the CIR estimates, the acoustic field is
computed using the plane or spherical wave interpolation as described in
Section IV. In the simulations, the carrier frequency is $f_{c}=32$ kHz and
the frequency bandwidth is $F_{d}=6$ kHz. In total, 10000 simulation trials
are run.

(a)

(b)

Figure 6: CDF of $\mathrm{MSD}_{\mathrm{interp}}$ for plane wave and spherical
wave interpolation: (a) the receiver is at a short range (within 50 m from the
signal source); (b) the receiver is at a long range (further than 1000 m away
from the signal source).

Fig. 6(a) shows the cumulative distribution function (CDF) of the
$\mathrm{MSD}_{\mathrm{interp}}$ when using the plane wave and spherical wave
interpolation. It can be seen that the performance of the spherical wave
interpolation is significantly better than that of the plane wave
interpolation. An $\mathrm{MSD}_{\mathrm{interp}}$ better than -15 dB is
achieved in $99\%$ of the simulation trials for the spherical interpolation,
while only in $86\%$ of the simulation trials for the plane wave
interpolation.

In Fig. 6(a), the spherical interpolation shows a significantly better
performance than the plane wave interpolation when the receiver is close to
the signal source (range smaller than 50 m). To investigate the performance at
a longer distance, we generate a single grid map at 20 m source depth for a
longer range scenario. The signal source depth is randomly chosen within
$[19.5,20.5]$ m. The range of the receiver varies from 1000 m to 1030 m, and
its depth varies from 15 m to 25 m. The resolution of the grid map in range
and depth are both 1 m. Again, the receiver position is randomly chosen to be
within the grid map. In total, 5000 simulation trials are run. As expected and
demonstrated in Fig. 6(b), almost the same performance is achieved by both
plane wave and spherical wave interpolation.

To conclude, the spherical interpolation is accurate regardless of the range
between the source and receiver, which is not the case for the plane wave
interpolation. As the trajectory of the target can be close to the source or
receiver, the spherical interpolation is a better candidate for the acoustic
field interpolation.

Figure 7: Time-varying CIR generated by the simulator based on the acoustic
field information.  Figure 8: Time-varying CIR estimate computed using the
received signal.

###  V-B Time-varying channel impulse response

In this subsection, we show an example of the time-varying CIR generated by
the simulator using the spherical wave interpolation when the target is moving
between the signal source and receiver.

This simulation scenario is setup based on one of the communication links in a
lake experiment described in Section VI. For the simulation, we use the first
355 s of the estimated target trajectory shown in Fig. 11. The target is
moving at a speed of 1 m/s at 1.5 m depth. The distance between two nodes is
14.2 m. The depth of the lake is set as 7 m. The transmitter is at 4.6 m depth
and the receiver is at 4.4 m depth. In this simulaton, two grid maps are
generated for a source depth of 1.5 m and 4.6 m depth. For each grid map, the
range of the receiver varies from 0 to 35 m, and its depth varies from 0 to 7
m. The resolution of the grid map in range and depth is 1 m and 0.5 m,
respectively.

As described in Section III-C, the CIR is computed as the inverse FFT of the
channel frequency response, given by (15). The CIR variation is due to the
second term in (15), where both the links, link 2 and link 3, have time-
varying channels. The magnitude of the CIR variation depends on the target
strength $\sigma_{TS}$, which is computed using (14). For the purpose of
observing the additional reflections introduced by the moving target, we set
$\sigma_{\mathrm{TS}}(i)=1$, i.e., at a high value, to highlight these
reflections in the CIR. The time-varying CIR generated by the simulator is
shown in Fig. 7. A group of static multipath components can be observed
between the $250$th and $500$th taps, which represents the acoustic
propagation between the source and receiver without the target. A large group
of weaker multipath components can be observed with time-varying delays, these
are the additional reflections introduced by the moving target. In this
simulation scenario, the target crossed the link between the source and the
receiver twice. This can be observed from the time-varying delays of the group
of reflections, which reduce when the target is approaching the link and
increase as the target moving away from the link.

The time-varying CIR allows us to understand the changes of the CIR structure
when a target is crossing the communication link. In practice, channel
estimates are obtained by a receiver. A detailed description of the receiver
used in this paper can be found in [23]. One of the features of the receiver,
useful for the target detection, is the use of turbo iterations for channel
estimation and demodulation, resulting in channel estimates obtained using
both the pilot and demodulated symbols, thus significantly improving the
estimation accuracy compared to the use of the pilot symbols alone. The
baseband received signal generated by the simulator is used for channel
estimation. The time-varying CIR estimate is shown in Fig. 8. Similar delay
variations of the group of refelctions is observed when the target is crossing
the communication link. This demonstrate that the channel estimator in the
receiver is capable of tracking the change of the CIR due to a moving target.
In Section VI, this receiver is used for generating the time-varying channel
estimate.

The simulator provides a mean to evaluate different algorithms for target
detection and has a flexibility of investigating different node configurations
within an area of interest. Note that organizing and conducting field
experiments with multiple nodes is complicated and expensive; therefore, a
simulator for investigation of target detection by a network of nodes is a
very important research tool.

With the two-stage processing, the computational time of simulating the
acoustic field for one communication link is significantly reduced. For this
simulation, two grid maps are pre-computed. The computation of each grid map
takes 25 minutes. For a simulation of 1 s duration, it takes 1.3 minutes to
generate the received signal at 12 kHz baseband sampling rate. If the acoustic
field is generated directly using BELLHOP3D, it takes about 9 hours to
generate 1 s of the received signal with the same sampling rate.

##  VI Lake experiment

In this section, we describe the lake experiment. Subsection VI-A describes
the configuration of the lake experiment. Subsection VI-B introduces the
transmitted signal. Subsection VI-C describes the technique used for target
trajectory estimation. In subsection VI-D, we demonstrate the target detection
capability. Subsection VI-E demonstrates the target detection performance with
the simulated data.

###  VI-A Experiment setup

Figure 9: Configuration of the lake experiment.

(a)

(b)

Figure 10: Target deployment: (a) floating platform with the acoustic
transmitter operating at the carrier frequency 24 kHz for the target
navigation; (b) underwater system attached to the floating platform, including
the yellow buoy with a weight and projector T2 connected to the transmitter on
the top of the floating platform by a cable.

The purpose of the lake experiment is to demonstrate the detection of a target
crossing a communication link between two nodes. Fig. 9 shows the
configuration of the lake experiment. Two projectors (T1 and T2) and five
hydrophones (H1 to H5) are deployed. T1 is the source node transmitting data
packets of 0.1 s duration and 1 s period at $32$ kHz carrier frequency, while
T2 is fixed on the target (the yellow buoy in Fig. 10(b)) transmitting data
packets at $24$ kHz carrier frequency for the purpose of the acoustic
navigation. T1 is fixed at 4.6 m depth. Hydrophones H2 to H5 are fixed at a
depth of 1.5 m above the lake bottom. Hydrophone H1 is fixed at a depth of 3 m
above the lake bottom. The depth of the lake site varies from 5 m to 8 m
depth.

The target is fixed underneath a $1\mathrm{m}\times 1\mathrm{m}$ floating
platform shown in Fig. 10(a). A weight is attached to the target to ensure it
is at approximately 1.5 m depth. During the experiment, a rope is used to pull
the floating platform with the target attached as shown in Fig. 10(b)
diagonally from the point $(22.7,0)$, near the hydrophone H4, to the point
$(1,7.2)$, or in the opposite direction. The speed of the target movement is
about 1 m/s. This target movement was repeated several times during the
experiment.

###  VI-B Data packet transmission

We consider data packet transmission in the network. A superimposed pilot and
data packet structure is used, with pilot symbols $p(i)$ in the real part and
data symbols $d(i)$ in the imaginary part [23]. Both the pilot and data are
binary phase-shift keying (BPSK) symbols. The data symbols are encoded by a
rate 1/3 convolutional code and interleaved. For the transmitted signal,
single-carrier modulation is considered. The superimposed pilot and data
symbols are pulse-shape filtered and up-sampled by a root raised-cosine (RRC)
filter with a roll-off factor of 0.2. After pulse-shaping, the signal is up-
converted to the carrier frequency.

###  VI-C Target trajectory estimation

To evaluate the target detection performance, the ground truth of the target
trajectory is required. To accurately estimate the trajectory of the target
during the experiment, the projector T2 attached to the target transmits the
target navigation signal in a frequency band outside of that of the
communication signals for the target detection (see Fig. 10(b)). A simple
navigation algorithm is then used based on estimating delays of the first
arrivals at the multiple receiving nodes H1 to H5. With the signal bandwidth
$F_{d}=6$ kHz and oversampling at the receiver to 12 kHz, the delay resolution
is 1/12 ms, thus the distance resolution is about 0.1 m, which is good enough
for the purpose of the experiment.

(a)

(b)

(c)

Figure 11: Coordinate estimates of the target; (a) trajectory of the target
during the experiment; (b) $x$ coordinate estimates; (c) $y$ coordinate
estimates.

Fig. 11(a) shows the estimated trajectory of the target during the experiment
as a black curve. The green circle shows the starting point and the green
diamond shows the final position of the target. Fig. 11(b) and 11(c) show how
the target coordinates $x$ and $y$, respectively, vary in time. Each slope in
Fig. 11(b) and 11(c) corresponds to a crossing of the target between the links
T1 and H1/H5. The other links are not useful, since the target moves very
closely to the other receivers (H2 to H4). Therefore, we use only these two
communication links for target detection.

###  VI-D Target detection

Figure 12: Channel impulse response between T1 and H5.

The first step of target detection is to obtain the CIR estimates between T1
and H1 and between T1 and H5, these are available at the receivers. As an
example, we show the estimates of the CIR between T1 and H5 in Fig. 12. It can
be observed that there is a shift in the delay of the multipath structure over
time. This is caused by the unsynchronized clocks on recorders used for data
transmission and reception. During the whole transmission interval of about 20
minutes, the delay of the direct path shifted by 0.8 ms. From Fig. 12, it is
not possible to visually identify the CIR changes corresponding to the target
crossing the communication link.

(a)

(b)

(c)

Figure 13: Normalized MSD in the experiment; (a) $x$-coordinate estimates; (b)
link between T1 and H1; (c) link between T1 and H5.

(a)

(b)

Figure 14: Normalized MSD with the simulated data: (a) Data transmission
starts from zero time instant; (b) Data transmission starts from 0.3 s.

The normalized MSD during the experiment for all 1200 transmitted data packets
according to (9) is plotted in Fig. 13(b) and 13(c) for the receivers H1 and
H5, respectively. We use the first CIR in the experiment, when the target was
in the position $(22.7,0)$, far from the communication links of interest, as
the reference impulse response. The parameters used in the filter $g$ in (3)
are: $P=15$, $K=2048$. The regularization parameter is $\varepsilon=10^{-3}$.
A large number of taps in filter $g$ is used to compensate for the delay shift
due to the unsynchronized clocks in the recorders and environment variation.

The potential time window of the target crossing is marked with red dotted
lines and arrows in Fig. 13. These time windows are chosen based on the
changes in the $x$ coordinate estimates. For the channel between T1 and H1,
two crossings are detected at the 263th and 617th data packets. For the
channel between T1 and H5, two crossings are detected at the 100th and 440th
data packets. Due to the one second gap in the data packet transmission, it is
possible to miss the detection of the target crossing. In this experiment,
four out of five target crossings are clearly detected based on the CIRs from
H1 and H5. This could be further improved by deploying more hydrophones at
different depths and/or reducing the gap between data packet transmissions.

###  VI-E Experiment-based simulation

This subsection shows the target detection performance obtained with the
simulated data generated based on the setup in lake experiment.

We simulate the channel between T1 and H1 for demonstration. During the
experiment, T1 and H1 are located at 4.6 m and 5.9 m depth, respectively. We
assume that the target is moving at a depth of 1.5 m. The coordinate estimates
shown in Fig. 11 are used for interpolating the trajectory of the target
during the experiment. The depth of the lake is set to 7 m. Two grid maps are
generated at source depths of 1.5 m and 4.6 m. Each grid map covers an area of
$35~{}\mathrm{m}$ in range and $7~{}\mathrm{m}$ in depth. The resolution of
the grid map is 1 m in range and 0.5 m in depth. The sound speed in the lake
is set based on the temperature of the water ($10^{\circ}$) as $c=1443$ m/s.

The target strength $\sigma_{TS}$ is computed using (14) with the radius of
the target set to 0.1 m. Random Gaussian noise is added to the time-varying
CIRs to limit the CIR estimation accuracy at a level of -20 dB in terms of
MSD. As there is no issue with the unsynchronized clock in the simulation, we
use $P=1$.

Fig. 14 shows the normalized MSD in two experiments. In the first experiment,
the data packets are transmitted from zero time instant, while for the second
experiment, the transmission starts from 0.3 s. It can be seen that four of
five target crossings are detected in the first experiment. In the second
experiment three crossings are detected. This is similar to what is observed
in the lake experiment. The simulation results demonstrate that the target
detection performance is sensitive to the timing of data packet transmissions
with respect to the target crossing. To increase the possibility of target
detection, more frequent data packet transmission is required.

##  VII Conclusions

In this paper, we have investigated the possibility of target detection by
using data transmission between communication nodes in an UWA network. A
target detection method has been proposed based on the changes in the time-
varying CIRs with respect to a reference CIR obtained without a target. An
acoustic simulator has been developed, which is capable of virtual
transmission of acoustic signals through a communication link in the presence
of a moving target and an approach has been proposed to improve the acoustic
field interpolation in the near field. The performance of the target detection
method has been investigated using numerical simulation and lake experiments.
The proposed detection metric allows clear indication of a moving target
crossing a communication link.

## References

  * [1] R. Urick, _Principles of Underwater Sound_. Peninsula Publishing, 1983.
  * [2] Y. Pailhas, Y. Petillot, K. Brown, and B. Mulgrew, “Spatially distributed MIMO sonar systems: principles and capabilities,” _IEEE Journal of Oceanic Engineering_ , vol. 42, no. 3, pp. 738–751, 2016.
  * [3] Y. Pailhas, J. Houssineau, Y. R. Petillot, and D. E. Clark, “Tracking with MIMO sonar systems: applications to harbour surveillance,” _IET Radar, Sonar & Navigation_, vol. 11, no. 4, pp. 629–639, 2017.
  * [4] V. Zverev, P. Korotin, A. Matveev, V. Mityugov, D. Orlov, B. Salin, and V. Turchin, “Experimental studies of sound diffraction by moving inhomogeneities under shallow-water conditions,” _Acoustical Physics_ , vol. 47, pp. 184–193, 2001.
  * [5] A. Matveev and V. Mityugov, “Determination of the parameters of motion for an underwater object,” _Acoustical Physics_ , vol. 48, no. 5, pp. 576–583, 2002\.
  * [6] A. Matveev, D. Orlov, A. Rodionov, B. Salin, and V. Turchin, “Comparative analysis of tomographic methods for the observation of inhomogeneities in a shallow sea,” _Acoustical Physics_ , vol. 51, pp. 218–229, 2005.
  * [7] A. L. Matveev, R. C. Spindel, and D. Rouseff, “Forward scattering observation with partially coherent spatial processing of vertical array signals in shallow water,” _IEEE Journal of Oceanic Engineering_ , vol. 32, no. 3, pp. 626–639, 2007.
  * [8] B. Lei, Y. Yang, K. Yang, and Y. Ma, “Detection of forward scattering from an intruder in a dynamic littoral environment,” _The Journal of the Acoustical Society of America_ , vol. 141, no. 3, pp. 1704–1710, 2017.
  * [9] N. Morozs, P. D. Mitchell, and Y. Zakharov, “Target detection using underwater acoustic networking,” in _IEEE OCEANS 2023-Limerick_ , 2023, pp. 1–5.
  * [10] N. Morozs, B. Parrein, L. Shen, B. T. Henson, and P. D. Mitchell, “Network protocols for simultaneous underwater acoustic communication and target detection,” in _Proceedings of the Underwater Acoustics Conference and Exhibition (UACE’23), Kalamata, Greece_ , 2023.
  * [11] L. Liao, B. Henson, and Y. Zakharov, “Grid waymark baseband underwater acoustic transmission model,” in _Underwater Acoustics Conference and Exhibition_ , 2017, pp. 343–350.
  * [12] M. Siderius and M. B. Porter, “Modeling broadband ocean acoustic transmissions with time-varying sea surfaces,” _The Journal of the Acoustical Society of America_ , vol. 124, no. 1, pp. 137–150, 2008.
  * [13] A. H. Sayed, _Fundamentals of adaptive filtering_. John Wiley & Sons, 2003.
  * [14] H. Yu, A. Song, M. Badiey, F. Chen, and F. Ji, “Iterative estimation of doubly selective underwater acoustic channel using basis expansion models,” _Ad Hoc Networks_ , vol. 34, pp. 52–61, 2015.
  * [15] C. Liu, Y. V. Zakharov, and T. Chen, “Doubly selective underwater acoustic channel model for a moving transmitter/receiver,” _IEEE Transactions on Vehicular Technology_ , vol. 61, no. 3, pp. 938–950, 2012.
  * [16] B. Henson, J. Li, Y. V. Zakharov, and C. Liu, “Waymark baseband underwater acoustic propagation model,” in _2014 Underwater Communications and Networking (UComms)_. IEEE, 2014, pp. 1–5.
  * [17] M. B. Porter, “BELLHOP3D User Guide,” Heat, Light, and Sound Research, Inc., Tech. Rep., 2016.
  * [18] R. P. Goddard, “The sonar simulation toolset, release 4.6: Science, mathematics, and algorithms,” Washington University Seattle Applied Physics Lab, Tech. Rep., 2008.
  * [19] “Sonar Acoustics Handbook,” Science and Technology Organization Centre for Maritime Research and Experimentation., Tech. Rep., 2016.
  * [20] K. Li and M. Chitre, “Data-aided underwater acoustic ray propagation modeling,” _IEEE Journal of Oceanic Engineering_ , vol. 48, no. 4, pp. 1127–1148, 2023.
  * [21] J. B. Allen and D. A. Berkley, “Image method for efficiently simulating small-room acoustics,” _The Journal of the Acoustical Society of America_ , vol. 65, no. 4, pp. 943–950, 1979.
  * [22] N. O. Booth, A. T. Abawi, P. W. Schey, and W. S. Hodgkiss, “Detectability of low-level broad-band signals using adaptive matched-field processing with vertical aperture arrays,” _IEEE Journal of Oceanic Engineering_ , vol. 25, no. 3, pp. 296–313, 2000.
  * [23] L. Shen, B. Henson, N. Morozs, Y. Zakharov, and P. D. Mitchell, “Rake receiver with interference cancellation for an underwater acoustic modem,” in _Proceedings of the Underwater Acoustics Conference and Exhibition (UACE’23), Kalamata, Greece_ , 2023.