††thanks: Published in the Journal of the Acoustical Society of America. The
published version of this preprint can be found at
https://doi.org/10.1121/10.0002974

# Resolution dependence of rough surface scattering using a power law
roughness spectrum

Derek R. Olson<EMAIL_ADDRESS>Oceanography Department, Naval Postgraduate
School, Monterey, CA 93943  Anthony P. Lyons  University of New Hampshire,
Durham, NH 03824

###### Abstract

Contemporary high-resolution sonar systems use broadband pulses and long
arrays to achieve high resolution. It is important to understand effects that
high-resolution sonar systems might have on quantitative measures of the
scattered field due to the seafloor. A quantity called the broadband
scattering cross section is defined, appropriate for high-resolution
measurements. The dependence of the broadband scattering cross section,
$\sigma_{bb}$ and the scintillation index, $SI$ on resolution was investigated
for one-dimensional rough surfaces with power-law spectra and backscattering
geometries. Using integral equations and Fourier synthesis, no resolution
dependence of $\sigma_{bb}$ was found. The incoherently-averaged frequency-
domain scattering cross section has negligible bandwidth dependence. $SI$
increases as resolution increases, grazing angle decreases, and spectral
strength increases. This trend is confirmed for center frequencies of 100 kHz
and 10 kHz, as well as for power-law spectral exponents of 1.5, 2, and 2.5.
The hypothesis that local tilting at the scale of the acoustic resolution is
responsible for intensity fluctuations was examined using a representative
model for the effect of slopes (inspired by the composite roughness
approximation). It was found that slopes are responsible in part for the
fluctuations, but other effects, such as multiple scattering and shadowing may
also play a role.

##  I Introduction

Theoretical treatment of wave scattering from rough interfaces is generally
performed using an incident monochromiatic plane wave, which has a single
direction and exists over infinite spatial extent. However, experimental
measurements of the scattered field often employ broadband pulses to achieve
high spatial resolution - desirable for seafloor mapping or target detection.
Performance of such systems typically depends on the mean intensity of the
scattered field from the seafloor, and more generally its probability density
function. For scattering from one-dimensional (1D) roughness, the mean
intensity is usually characterized in terms of the scattering cross-section
per unit length per unit angle111or per unit area per unit solid angle for
two-dimensional (2D) rough interfaces, $\sigma$ – hereafter referred to as the
“cross section”, “scattering cross section,” or “scattering strength” for the
decibel version. Variability of the scattered intensity is often characterized
using the scintillation index, $SI$ Tatarski (1961); Ishimaru (1978).

The scattering cross section for a frequency-domain incident field for one-
dimensional roughness is defined as Thorsos (1988)

$\displaystyle\sigma_{f}=\frac{\langle I_{s}\rangle R}{I_{i}L_{eff}}$ (1)

where $\langle I_{s}\rangle=\langle|p_{s}(f)|^{2}\rangle/(2\rho_{0}c)$ is the
mean scattered intensity in the far-field, $I_{i}=|p_{0}|^{2}/(2\rho_{0}c)$ is
the incident intensity in the direction of the incident wave vector,
$\rho_{0}$ is the ambient density, $p_{0}$ is the complex amplitude of the
incident plane wave, $L_{eff}$ is the effective ensonified length of the
incident field, and $R$ is the distance between a patch of rough interface and
the receiver location. Note that this definition is valid only for geometries
with well-defined incident and scattered field directions and for a single
frequency. Strictly, this definition of the scattering cross section is only
true in the limit as the ensonified length becomes large compared to all
length scales of interest (i.e. outer scale of the rough surface, or the
acoustic wavelength), since monochromatic plane waves interact with the entire
rough surface. This definition is often used for narrowband222In this work,
when the bandwidth of the signal is a tenth of the center frequency or
smaller, it is considered to be narrowband. incident fields that provide a
good approximation to a single-frequency tone.

In this work, a quantity termed the “broadband scattering cross section” is
investigated, that is appropriate for cases with short pulses333To be
consistent with the criterion for narrowband defined above, a “short” pulse is
defined to contain ten or fewer cycles.. For a plane-wave rectangular pulse of
length $\tau$, this is defined as

$\displaystyle\sigma_{bb}=\frac{R}{c\tau/(\cos\theta_{i}+\cos\theta_{s})}\frac{\langle
I_{s}\rangle}{I_{i}}$ (2)

where $c$ is the sound speed, $\theta_{i}$ is the incident grazing angle, and
$\theta_{s}$ is the scattered grazing angle. This quantity is discussed more
fully and derived in Sec. V. In this work, the frequency-domain version of the
scattering cross section is referred to as the “cross section”, and the
broadband version is $\sigma_{bb}$. The broadband scattering cross section may
exhibit pulse-length dependence if the properties of the ensemble of rough
surfaces vary with resolution, especially for high-resolution systems. This
dependence is not possible for the frequency-domain cross section.

The interface scattering cross-section (or its broadband version)
characterizes the mean scattered power from an interface, but a more general
property of the scattered field is the probability density function (pdf) of
the modulus of the complex pressure, termed the envelope pdf. The envelope pdf
is connected to performance of target detection systems, and has potential
utility for use in remote sensing of the environment using high resolution
systems Lyons _et al._ (2009, 2016); Olson _et al._ (2019). The statistical
distribution of pressures resulting from scattering from a rough, homogeneous
interface with Gaussian height statistics has commonly been assumed to be
Gaussian for the real and imaginary components, and a Rayleigh distribution
for the envelope Jakeman (1980)). In this situation, the scintillation index,
or normalized intensity variance is unity. For heavy-tailed statistics (with
more frequent large amplitude events), the scintillation index is greater than
unity.

Arguments for Rayleigh-distributed scattered pressure magnitude follow from
the assumption of a large number of independent surface elements contributing
to the scattered field Jakeman (1980); Abraham and Lyons (2002). So long as
the ensonfied area of a Gaussian, homogeneous rough interface is large (so
there are many independent scatterers contributing to the field), this
assumption holds true. Another argument for Rayleigh magnitudes follows from
perturbation theory and the interpretation in terms of Bragg scattering. In
this framework, the scattered pressure is proportional to the amplitude
spectrum of the roughness evaluated at the Bragg wavenumber,
$2k_{w}\cos(\theta_{i})$, where $k_{w}$ is the acoustic wavenumber in the
water column, and $\theta_{i}$ is the incident grazing angle. If the surface
has Gaussian statistics in the spatial domain, then the wavenumber components
will have Rayleigh-distributed magnitude via the central limit theorem.
Therefore, the acoustic envelope pdf will be Rayleigh-distributed as well.

Contemporary high-resolution seafloor imaging systems, such as synthetic
aperture sonar (SAS) have spatial resolutions on the order of the center
wavelength Fossum _et al._ (2008); Bellettini and Pinto (2009); Pinto (2011);
Dillon (2018); Sternlicht _et al._ (2016). Small resolution cell sizes may
result in ensembles that vary with the resolved area of the seafloor, thereby
causing a departure from Rayleigh statistics. The resolution dependence of the
scintillation index has implications for target detection performance,
synthetic aperture autofocus algorithms (e.g. Marston and Plotnick (2015)),
and preprocessing algorithms for automated target recognition Kwon and
Nasrabadi (2005); Williams (2015); Galusha _et al._ (2018).

It was observed in Lyons _et al._ (2016) that measurements of the
scintillation index from SAS images of homogeneous random rough interfaces had
a strong dependence on range, which was interpreted as a result of modulation
of the local slope by roughness components at the scale of the acoustic
resolution or larger, which is called the local-tilting hypothesis. This
interpretation uses a model for local slope modulation inspired by the
composite roughness approximation McDaniel and Gorman (1983). Combined with
interpretations in Lyons _et al._ (2016), this effect results in a dependence
of the scintillation index on the acoustic resolution, the underlying pixel
statistics, range (through grazing angle), as well as roughness spectrum
parameters. These interpretations, while plausible, suffer from a lack of
experimental confirmation. In the electromagnetics literature, slope
modulation was postulated to cause non-Rayleigh scattering by Valenzuela and
Liang (1971), and Li and Johnson (2017). In the specular direction, the
Fresnel zone (a form of resolved area) has been shown to affect the
scintillation index Yang and McDaniel (1991), although this situation contains
a coherent component and is not germane to the current problem.

In this work, the question is examined of whether there is a dependence on
resolution of the broadband scattering cross section and scintillation index.
The acoustic resolution, $\Delta X$ is defined as the full width half maximum
spatial extent of the square of the incident pulse envelope. For broadband
signals, the temporal resolution is set by $1/(2aB_{3dB})$, where $B_{3dB}$ is
the 3 dB full width bandwidth of the transmitted pulse444If phase coding and
pulse compression are employed, this analysis deals with the resolution of the
compressed, or matched-filtered pulse., and $a$ is a constant which depends on
the shape of the pulse used. The spatial resolution is $c/(2aB_{3dB})$ for
small grazing angles.

These questions were investigated through numerical solution of the Helmholtz-
Kirchhoff integral equation for the scattered pressure using the boundary
element method (BEM) Sauter and Schwab (2011); Wu (2000) using pressure-
release boundary conditions. This method is similar to that used by Thorsos
(1988). Fourier synthesis was used to construct the broadband scattered
pressure at various spatial resolutions, and metrics were computed based on
the scattered time-domain pressure. The numerical method detailed here can, in
general, treat bistatic geometries, but only the monostatic case was examined.
Comparisons were made to the ensemble averaged cross section performed in the
frequency domain (i.e. computed at a single frequency, which is a good
approximation of the monochromatic plane wave case). These simulations were
performed for center frequencies of both 100 kHz and 10 kHz, and for one-
dimensional rough surfaces with power law spectra, whose parameters are the
spectral strength and spectral exponent.

Comparisons between numerical simulations and field experiments would require
two-dimensional (2D) roughness with a three-dimensional (3D) scattering
geometry. However, the numerical method has very high computational and memory
requirements that made 3D simulations impossible at this time, and therefore a
2D geometry was used instead. Although the values of the scattering cross
section and scintillation index will be different in 2 and 3D geometries, the
fundamental scattering phenomena, such as Bragg scattering, local tilting,
multiple scattering, and shadowing will be present in both 2D and 3D.

Through these numerical experiments, it was found that the broadband
scattering strength does not vary as a function of bandwidth for the
parameters investigated in this study. The error of this comparison is within
Monte-Carlo error. For scintillation index it was found that it becomes
greater than one as resolution increases, grazing angle decreases, and
spectral strength increases. For larger spectral exponents, the scintillation
index is more sensitive to changes in spectral strength and grazing angle.

An overview of the geometry and roughness statistics is presented in Sec. II.
The integral equations and discretization methods are given in Sec. III, and
the incident field in Sec. IV. Methods to estimate the broadband scattering
cross section and scintillation index are given in Sec. V. A discussion on how
the parameters of the numerical simulations were selected is given in Sec. VI.
Results are presented in Sec. VII, with a discussion and some preliminary
hypotheses given in Sec. VIII. Conclusions are given in Sec. IX

##  II Geometry and Environment

The geometry of the scattering problem is presented in Fig. 1. The problem
takes place in two dimensions with position vector $\textbf{r}=(x,z)$. The
acoustic medium is above the rough surface, defined as $z_{s}=f(x_{s})$ and
shown as the thick black line in this figure. The coordinates $(x_{s},z_{s})$
are points on the rough surface. In this figure, the nominal incident and
scattered wave directions are shown with their grazing angles and nominal wave
vectors. The slope angle is defined as
$\epsilon=\tan^{-1}\left(df(x)/dx\right)$ and is shown in the figure, along
with the normal vector pointing out of the acoustic medium (into the lower
half-space). The sound speed in the upper medium is $c$, which is taken to be
1500 m/s, but these results can be applied to other sound velocities by
performing the appropriate dimensional scaling. The acoustic frequency is $f$,
and is related to the wavenumber by $k=2\pi f/c$. Simulations are performed at
a center frequency $f_{0}$, and 3 dB bandwidth $B_{3dB}$. The center
wavelength and wavenumber are $\lambda_{0}$ and $k_{0}$ respectively

Figure 1: Rough surface scattering geometry. The nominal incident and
scattered wave vector are shown, along with the slope angle, $\epsilon$, at
the origin, and the unit normal vector at that point, which is defined to
point out of the acoustic domain, which is defined as the upper medium.

The rough interface is assumed to have wide-sense homogeneity (spatial
stationarity) and a Gaussian pdf for both height and slope. By assuming
Gaussian statistics, the rough-interface second-order properties can be
completely described by its autocovariance function,

$\displaystyle B(x)=\langle f(y)f(y+x)\rangle$ (3)

and power density spectrum

$\displaystyle W(K)=\frac{1}{2\pi}\int B(x)e^{iKx}\,\textrm{d}x.$ (4)

Several second-order properties of this spectrum are useful for the analysis
performed here. In particular, the mean square height, $h^{2}$ is given by

$\displaystyle h^{2}=\int\limits_{-\infty}^{\infty}W(K)\,\textrm{d}K=B(0).$
(5)

The mean square slope $s^{2}$ is

$\displaystyle
s^{2}=\int\limits_{-\infty}^{\infty}K^{2}W(K)\,\textrm{d}K=-\left.\frac{\partial^{2}B(x)}{\partial
x^{2}}\right|_{x=0}.$ (6)

The power density spectrum used in this work is the truncated power law,

$\displaystyle W(K)=\frac{w}{|K|^{\gamma}}$ (7)

for $k_{l}\leq|K|\leq k_{u}$, and zero otherwise. The spectral strength is $w$
with units of m3-γ, and $\gamma$ is the dimensionless spectral exponent. The
lower wavenumber cutoff is $k_{l}=2\pi/L_{0}$, where $L_{0}$ is the outer
scale. The upper wavenumber cutoff is $k_{u}=\pi/\ell_{0}$, where $\ell_{0}$
is the inner scale. The extra factor of $1/2$ in defining $k_{u}$ is chosen
such that the interval $[-k_{u},k_{u}]$ has a total length of $2\pi/\ell_{0}$.
Random realizations are produced from this power spectrum using the Fourier
synthesis technique from Thorsos (1988), and is given for completeness in
Appendix A. The outer scale is specified independently of the surface length,
$L$, and is required to satisfy $L_{0}<L$. Similarly, the inner scale
satisfies $\ell_{0}>\delta x$, where $\delta x$ is the sampling interval of
the rough interface realization. Although not required in general, the inner
scale is smaller than the smallest wavelength with significant energy in these
broadband simulations.

For the power-law form used here, the non-dimensional mean square slope and
mean square height are

$\displaystyle s^{2}$
$\displaystyle=\frac{2k_{0}^{3-\gamma}w}{3-\gamma}\left[\left(\frac{k_{u}}{k_{0}}\right)^{3-\gamma}-\left(\frac{k_{l}}{k_{0}}\right)^{3-\gamma}\right]$
(8) $\displaystyle k_{0}^{2}h^{2}$
$\displaystyle=\frac{2k_{0}^{3-\gamma}w}{\gamma-1}\left[\left(\frac{k_{0}}{k_{l}}\right)^{\gamma-1}-\left(\frac{k_{0}}{k_{u}}\right)^{\gamma-1}\right].\,,$
(9)

where $k_{0}$ is the center wavenumber defined in Eq. (16). These parameters
have been expressed in a form where the terms outside and inside the brackets
parentheses are dimensionless. L’Hôpital’s rule can be used to show that the
mean square slope is finite for $\gamma=3$, and the mean square height is
finite for $\gamma=1$.

The true upper wavenumber is set by the inner scale, $k_{u}=\pi/\ell_{0}$.
However, the way in which the rough surfaces enter into the acoustical
simulations may be subject to an effective upper limit,
$k_{u}^{\prime}=\pi/\ell^{\prime}$, where $\ell^{\prime}$ is an effective
inner scale. Roughness wavelengths much less than $\lambda_{0}$ likely have an
insignificant effect on the scattered field. Thus, the effective upper
wavenumber is likely much less than that defined by the surface sampling. The
effective upper limit likely does not affect the scales causing scattering
near the Bragg wavelength, but rather sets the upper wavenumber limit for
computing the large-scale slope in the slope modulation model examined below.
For $\gamma>1$, root-mean-square height is insensitive to the upper cutoff,
and more sensitive to the low-wavenumber cutoff. For $\gamma<3$, the root-
mean-square slope is sensitive to the upper cutoff, and insensitive to the
lower cutoff (so long as it is sufficiently small). As $k_{u}$ becomes large,
$s$ grows without bound, and an effective upper limit can alleviate this
problem. To make the effective upper wavenumber limit explicit, the notation
$s_{\ell^{\prime}}$ is used to denote the rms slope computed using
$k_{u}^{\prime}=\pi/\ell^{\prime}$.

##  III Integral equations and discretization

This study was performed numerically using a discretized form of the 2D
Helmholtz-Kirchhoff integral equation for Dirichlet boundary conditions
Thorsos (1988). Although the motivation for this work is seafloor scattering,
the assumption of a Dirichlet boundary allows us to focus solely on the role
of the rough interface. For a single frequency, this integral equation,
defined on the rough interface is Wu (2000),

$\displaystyle p_{i}\left(\textbf{r}_{p}\right)=-\int\limits_{S}\frac{\partial
p\left(\textbf{r}_{s}\right)}{\partial
n_{s}}G_{k}\left(|\textbf{r}_{s}-\textbf{r}_{p}|\right)\textrm{d}S,$ (10)

where $p_{i}$ is the incident pressure, $\textbf{r}_{p}=(x_{p},z_{p})$ and
$\textbf{r}_{s}=(x_{s},z_{s})$ are points on the rough surface (with a
subscript $s$ denoting the integration variable), $\partial p/\partial n_{s}$
is the total pressure normal derivative at $\mathbf{r}_{s}$, with the normal
direction pointing out of the acoustic domain (into the lower space), and
$G_{k}(|\textbf{r}_{s}-\textbf{r}_{p}|)=(i/4)H_{0}^{(2)}(k|\textbf{r}_{s}-\textbf{r}_{p}|)$
is the 2D free-space Green function Wu (2000), where $H_{0}^{(2)}(z)$ is the
zeroth-order Hankel function of the second kind. Note that the Green function
used here differs from that used in Thorsos (1988) due to the differing time
convention. The normal vector points downward here, as opposed to upward in
Thorsos (1988), although in both cases it points out of the acoustic domain.
This integral equation can also describe electromagnetic scattering from 1D
corrugated surfaces with perfectly conducting boundary conditions subject to
an incident wave with transverse magnetic (also known as p) polarization
Toporkov _et al._ (1998).

The scattering problem is solved in two steps. First, Eq. (10) is numerically
solved for $\partial p/\partial n$ on the surface, through discretization of
the integral equation using the boundary element method Sauter and Schwab
(2011); Wu (2000). In particular, piecewise linear basis functions are used to
approximate $\partial p/\partial n$, and collocation to compare the true and
approximate solution at discrete points. These two steps convert the integral
equation into a linear system,

$\displaystyle Vu=b,$ (11)

where $u$ is the solution vector consisting of the basis function coefficients
used in the approximation for $\partial p/\partial n$, and $b=p_{i}$ evaluated
at the discrete collocation points $\textbf{r}_{m}=(x_{m},z_{m})$. The matrix
$V$ has elements

$\displaystyle V_{mn}=-\int
G_{k}\left(|\textbf{r}_{m}-\textbf{r}_{s}|\right)\phi(\xi_{n}(\textbf{r}_{s}))\textrm{d}S.$
(12)

Here, $\phi(\eta)$ is a linear basis function defined on the interval
$\eta\in[-1,\,1]$. Outside of the interval, $\phi$ is zero. The function
$\xi_{n}$ maps the basis function centered at the $n$-th point from physical
space, $\textbf{r}_{s}$ to the $\eta$ domain. In this case, the basis
functions are centered at the same collocation points $\textbf{r}_{m}$,
resulting in a square matrix. Integration is carried out using a 4 point
Gauss-Legendre quadrature rule Abramowitz and Stegun (1972) for nonsingular
elements. Due to the weak singularity in the Green function, the diagonal
elements of the matrix are computed using a 16 point quadrature rule combined
with a variable transformation whose Jacobian exactly cancels the singularity
Wu (2000). LAPACK routines were used to solve the linear system using LU
decomposition and back substitution Anderson _et al._ (1999).

Collocation points are defined on the rough surface, $(x_{m},z_{m})$ with
equal spacing, $\delta x$ on the horizontal axis. The method to generate these
points is given in Appendix A. From these points, a cubic approximation is
used to construct a continuous and smooth surface. This interpolation process
forces the surface normal, and thus $\partial p/\partial n$ to be continuous,
which improves the convergence rate of the discretization of the integral
operator Atkinson (1997). The interpolation scheme may extend the region of
wavenumber support beyond $\pm\pi/\ell_{0}$, affecting the estimate of the rms
slope and height. However, monotonic piecewise Hermite interpolation Fritsch
and Carlson (1980) was used, which does not suffer from overshoot and has a
negligible effect on the rms slope and height estimates.

Once the surface pressure normal derivative is found, the scattered pressure
at a field point in the domain, $p_{s}(\mathbf{r}_{f})$, is found using,

$\displaystyle p_{s}(\textbf{r}_{f})=\int\limits_{S}\frac{\partial
p\left(\textbf{r}_{s}\right)}{\partial
n_{s}}G_{k}\left(|\textbf{r}_{s}-\textbf{r}_{f}|\right)\textrm{d}S,$ (13)

where $\textbf{r}_{f}=(x_{f},z_{f})$. In this work, the field points are
equally spaced intervals of one degree at a distance $R$ from the rough
surface, which is approximately 25 times the Rayleigh distance from the
surface, $L^{2}/\lambda_{0}$, where $L$ is the total surface length and
$\lambda_{0}$ is the center wavelength. This criterion for the far-field is
very conservative (Jackson and Richardson, 2007, Appendix J), Winebrenner and
Ishimaru (1986); Lysanov (1973), although it enables the use of asymptotic
expansions for the Hankel function. The scattered angle, $\theta_{s}$ is
related to the field point locations using $x_{f}=R\cos\theta_{s}$, and
$z_{f}=R\sin\theta_{s}$. The variable $p_{s}(f,\theta_{s},\theta_{i},R)$ is
the complex pressure $p_{s}$ measured at a location $\mathbf{r}_{f}$, produced
by an incident field with frequency $f$ and nominal incident grazing angle
$\theta_{i}$.

##  IV Incident Field

The incident fields used in this work are broadband pulses whose spatial
dependence is an approximation of a plane wave. The nominal directions of the
incident and scattered wave vectors are indicated as arrows in Fig. 1. The
incident and scattered wave vector lengths vary due to the broadband nature of
the field, although the center wave vectors can be defined at the center
frequency by the expressions $\textbf{k}_{0i}=(k_{0ix},k_{0iz})$ and
$\textbf{k}_{0s}=(k_{0sx},k_{0sz})$. The components are defined in terms of
the grazing angles $\theta_{i}$ and $\theta_{s}$ (with respect to the
horizontal axis) by

$\displaystyle k_{0ix}$ $\displaystyle=-k_{0}\cos\theta_{i}$ $\displaystyle
k_{0sx}$ $\displaystyle=k_{0}\cos\theta_{s}$ (14) $\displaystyle k_{0iz}$
$\displaystyle=-k_{0}\sin\theta_{i}$ $\displaystyle k_{0sz}$
$\displaystyle=k_{0}\sin\theta_{s}.$ (15)

The incident unit wave vector is
$\hat{\mathbf{k}}_{0i}=\mathbf{k}_{0i}/k_{0}$, and the scattered unit vector
is $\hat{\mathbf{k}}_{0s}=\mathbf{k}_{0s}/k_{0}$.

The center wavenumber, $k_{0}$, is defined by an average of the wavenumber
weighted by power spectrum of the transmitted source

$\displaystyle k_{0}=\frac{\int\limits_{-\infty}^{\infty}\frac{2\pi
f}{c}S^{2}(f)\,\textrm{d}f}{\int\limits_{-\infty}^{\infty}S^{2}(f)\,\textrm{d}f},$
(16)

where $S(f)=\int s(t)\exp(-i2\pi ft)\,\textrm{d}t$ is the linear (amplitude)
spectrum of the transmitted pulse, $s(t)$. The transmitted pulse used here was
a complex exponential multiplied by a Gaussian envelope, with the time- and
frequency-domain forms,

$\displaystyle s(t)=\exp\left(-t^{2}/\tau^{2}+i\omega_{0}t\right)$ (17)
$\displaystyle
S(f)=\tau\sqrt{\pi}\exp\left(-(f-f_{0})^{2}\pi^{2}\tau^{2}\right)\,,$ (18)

where $\tau$ is a parameter of the pulse length, and $\omega_{0}=2\pi f_{0}$
is the center angular frequency. Using the Gaussian form for $S(f)$, the
numerator in Eq. (16) is $\pi\sqrt{2\pi}f_{0}\tau/c$, and the denominator is
$\sqrt{\pi/2}\tau$. The ratio is $2\pi f_{0}/c=k_{0}$.

The temporal resolution of the pulse, $\Delta\tau$ is defined by the duration
of the pulse envelope between its half power points. For the Gaussian pulse
used, this quantity can be obtained by solving the equation
$\exp\left(-(\Delta\tau/2)^{2}/\tau^{2}\right)=1/\sqrt{2}$, resulting in
$\Delta\tau=\tau\sqrt{2\ln 2}$. The 3 dB bandwidth, $B_{3dB}$, of the pulse is
defined as the full-width, half maximum of $|S(f)|^{2}$, namely
$|S(f_{0})|^{2}/2=|S(f_{0}\pm B_{3dB}/2)|^{2}$. These definitions of the
temporal resolution and bandwidth result $B_{3dB}\Delta\tau=2\ln(2)/\pi\approx
0.44$. For reference, if a rectangular function with full width of $B_{3dB}$
is used for $S(f)$, then $B_{3dB}\Delta\tau\approx 0.88$. The same
relationship is obtained if constant envelope pulse of length $\Delta\tau$ is
used. Although the rectangular pulse has a larger time-bandwidth product, the
Gaussian pulse has no sidelobes in the time-domain, but requires a
computational bandwidth much larger than $B_{3dB}$ to approximate a true
Gaussian function. The equivalent noise bandwidth (EQNB)Harris (1978) is also
needed, which for the Gaussian pulse is $B_{EQNB}=(\tau\sqrt{2\pi})^{-1}$. The
constant $A$ is defined in terms of the product

$\displaystyle\Delta\tau B_{EQNB}=A^{-1}\,,$ (19)

so that $A=\sqrt{\pi/\ln 2}$ for the Gaussian pulse. This constant is used
later to define an effective inner scale to the rough interface based on its
acoustic resolution.

Broadband fields are synthesized from single frequency approximations of a
plane wave. This narrowband incident field is the extended Gaussian beam
developed by Thorsos (1988) that provides tapering to guard against edge
effects entering into the scattering calculation. The form of this field
(adapted to this time convention) is given by

$\displaystyle\begin{split}p_{i}(\textbf{r}_{s},f)=p_{0}&\exp\left(-i\textbf{k}_{i}\cdot\textbf{r}_{s}\left(1+w_{t}(\textbf{r}_{s})\right)\right)\\\
\times&\exp\left(-\left(x_{s}-z_{s}\cot\theta_{i}\right)^{2}/g^{2}\right),\end{split}$
(20)

where

$\displaystyle
w_{t}(\textbf{r}_{s})=(kg\sin\theta_{i})^{-2}\left[2\left(x_{s}-z_{s}\cot\theta_{i}\right)^{2}/g^{2}-1\right],$
(21)

$g$ is a width parameter of the incident field, $\mathbf{r}_{s}$ is a point on
the rough surface, and $\theta_{i}$ is the nominal incident grazing angle. The
factor $p_{0}$ has units of Pa, and is helpful for keeping track of units when
estimating the scattering cross section, or broadband cross section. For
broadband simulations, Eq. (20) is used for each frequency. Since the Gaussian
function has an infinite domain of support, it must be truncated to use in
numerical simulations.

The function $w_{t}\left(\textbf{r}\right)$ improves the agreement between the
numerical solution of the Helmholtz-Kirchhoff integral equation of the first
kind, Eq. (10), and the integral equation of the second kind. Discrepancies
between these two solutions can result because the incident field satisfies
the Helmholtz equation approximately to order $(kg\sin\theta_{i})^{-2}$
Thorsos (1988). Good agreement between the two solutions was observed when
$kg\sin\theta_{i}$ is large. Therefore, to analyze low-grazing angles, (which
are important to contemporary synthetic aperture sonar systems, e.g. Fossum
_et al._ (2008); Bellettini and Pinto (2009); Pinto (2011); Dillon (2018);
Sternlicht _et al._ (2016)), the parameter $g$ must grow as $\theta_{i}$
approaches zero. This requirement can be thought of as enforcing the
constraint that the angular width of the incident beam (full width half max),

$\displaystyle\Delta\theta=\frac{2\sqrt{2\ln(2)}}{kg\sin\theta_{i}},$ (22)

should be small compared to $\theta_{i}$. When the relative angular width,
defined as $\Delta\theta/\theta_{i}$ is not small, the direction of the
incident field is spread over a large range of angles compared to the incident
grazing angle. This situation should be avoided at low grazing angles, since
the scattering cross section is commonly rapidly varying in that region.

##  V Estimating time-domain quantities of the scattered field

In this section, expressions are derived for the broadband scattering cross
sections and scintillation index of the scattered field due to a broadband
incident pulse. The time-domain pressure is computed by

$\displaystyle
p_{s}(t^{\prime},\theta_{s},\theta_{i},R)=\int\limits_{-\infty}^{\infty}S(f)p_{s}(f,\theta_{s},\theta_{i},R)e^{i2\pi
f(t^{\prime}+R/c)}\,\textrm{d}f$ (23)

where $p_{s}(f,\theta_{s},\theta_{i},R)$ is the scattered pressure measured at
grazing angle $\theta_{s}$, distance from the origin, $R$, incident field
frequency $f$, and nominal incident grazing angle $\theta_{i}$. The factor
$R/c$ in the exponent removes the time delay associated with propagation to
the far field when the incident pulse is centered on the origin, so that the
scattered pressure time series can be mapped to the mean plane of the rough
interface. This delayed time is denoted $t^{\prime}$, which is used below in
the calculation of the effective energy flux. The scattered grazing angle,
$\theta_{s}$ is computed using the location at which the pressure is
calculated in the far field, $\theta_{s}=\tan^{-1}(z_{f}/x_{f})$. In practice
the integral is computed using the fast Fourier transform.

The broadband scattering cross section, $\sigma_{bb}$, is computed directly
from the scattered pressure in the time domain. Although there is no general
framework for this kind of quantity, motivation is provided here by comparison
to the definition of the scattering cross section in the frequency domain. In
practice, finite resolution in scattering measurements is sometimes obtained
using short pulses Urick (1954, 1983), and this definition accords with the
sonar equation used in those situations.

The frequency domain version of the cross section, $\sigma_{f}$ for a
monochromatic plane wave from a rough interface of length $L$, was given in
Sec. I as Eq. (1), due to Thorsos (1988). This expression can be cast in terms
of the incident energy-flux (Landau and Lifshitz, 1987, p. 255),Thorsos (1988)
passing through the plane $z=0$

$\displaystyle\sigma_{f}=\frac{\langle I_{s}\rangle R\sin\theta_{i}}{E_{f}},$
(24)

where

$\displaystyle
E_{f}=\int\limits_{-\infty}^{\infty}\left[\left.\frac{1}{2}\mathrm{Re}\\{p_{i}\mathbf{v}_{i}^{\ast}\\}\cdot\hat{n}\right.\right|_{z=0}\,\mathrm{d}x,$
(25)

is the total energy flux passing through the $z=0$ plane, $p_{i}$ is the
incident complex pressure, $\mathbf{v}_{i}$ is the incident acoustic particle
velocity, and $\hat{n}$ is the unit normal vector of the line $z=0$. Since the
energy is directed downards towards the rough interface, $\hat{n}$ was chosen
to be $-\hat{z}$. In Thorsos (1988), the incident acoustic energy was directed
upwards towards the rough interface, and correspondingly, their normal vector
pointed up. Making this substitution, the energy flux is

$\displaystyle
E_{f}=-\int\limits_{-\infty}^{\infty}\left[\left.\frac{1}{2}\mathrm{Re}\\{p_{i}{v}_{iz}^{\ast}\\}\right.\right|_{z=0}\,\mathrm{d}x,$
(26)

where $v_{iz}$ is the vertical component of the incident acoustic particle
velocity. For an untapered plane wave incident on a rough surface of length
$L$, $E_{f}=I_{i}L\sin\theta_{i}$, and the definition in Eq. (1) is recovered.

This definition can be extended to an incident plane wave with a time
dependence set by $e^{i\omega_{0}t}$ times a rectangular pulse of length
$\tau$. Here, a new quantity is defined called the broadband scattering cross
section, $\sigma_{bb}$. It might seem reasonable to start with Eq. (24) with
$E_{f}$ evaluated in the time-domain. In this case,
$E_{f}=I_{i}c\tau\sin\theta_{i}/\cos\theta_{i}$, and thus
$\sigma_{bb}=(\langle I_{s}\rangle R)/(I_{i}c\tau/\cos\theta_{i})$. However,
this definition is insufficient because in the case of an incident pulse, the
(effective) ensonified length that is relevant for the scattering cross
section is set by both the incoming and outgoing angles. Physically, the
dependence on the scattered field point (or angle) is due to the fact that at
a given instant in time, only part of the ensonified surface contributes to
the total field measured at a point in the far field.

To account for this effect, the vertical component of the energy flux density
can be written as a function of space and time

$\displaystyle e_{f}(t,\textbf{r})=\frac{1}{2}Re\\{p_{i}v_{iz}^{\ast}\\}\,,$
(27)

which is the integrand of (26). Instead of integrating this quantity over $x$
to obtain the energy flux, a delayed version of $e_{f}$ is integrated over $x$
with $z=0$ to obtain the effective energy flux. This quantity is defined as

$\displaystyle
E_{f}^{\prime}=-\int\limits_{-\infty}^{\infty}e_{f}(t-t_{s},\mathbf{r})|_{z=0}\,\mathrm{d}x\,,$
(28)

where

$\displaystyle
t_{s}=|\mathbf{r}_{f}-\mathbf{r}|/c\approx(R-x\cos\theta_{s}-z\sin\theta_{s})/c\,,$
(29)

and $R=|\mathbf{r}_{f}|$. This time-delay takes into account the alteration in
the region that contributes to the scattered field at a given instant in time
due to the position of the field point. For a line source, or other compact
source configuration, a similar analysis can be performed by taking into
account the constant-time (isochronous) ellipse for the transmitted pulse and
geometry.

The broadband scattering cross section is defined in terms of the effective
incident energy-flux, analogous to the frequency-domain version

$\displaystyle\sigma_{bb}=\frac{\langle I_{s}\rangle
R\sin\theta_{i}}{E_{f}^{\prime}}\,.$ (30)

Note that $E_{f}^{\prime}$ may in general be a function of time, since the
pulse may be attenuated by a spatial taper (as it is here), or by the beam
pattern of a transducer in experiments. To use the entire time series
scattered by the rough interface due to the broadband pulse, the energy flux
must be calculated as a function of time. For a plane-wave (untapered) pulse
the energy flux is independent of time, for all pulse shapes.

First, this defintion will be demonstrated for a plane wave with a rectangular
pulse shape, to build intuition. The incident pressure in this case is given
by

$\displaystyle
p_{i}(\mathbf{r},t)=p_{0}e^{i\omega_{0}(t-\hat{\mathbf{k}}_{i}\cdot\mathbf{r}/c)}\Pi\left(\frac{t-\hat{\mathbf{k}}_{i}\cdot\mathbf{r}/c}{\tau/2}\right)$
(31)

where $\mathbf{r}$ is a general point in space, $\hat{\mathbf{k}}_{i}$ is the
incident unit wave vector, $\Pi(x)=1$ if $|x|\leq 1$ and zero otherwise, and
sets the incident pulse duration. The effective energy flux is then
$E_{f}^{\prime}=I_{i}c\tau\sin\theta_{i}/\nu_{x}$, where
$\nu_{x}=\cos\theta_{i}+\cos\theta_{s}$. For a broadband rectangular pulse,
the sine factors cancel, and

$\displaystyle\sigma_{bb}=\frac{\langle I_{s}\rangle R}{I_{i}c\tau/\nu_{x}}.$
(32)

The denominator now includes dependence on the scattered direction, and
contains the effective ensonified length, $L_{eff}$. The form
$L_{eff}=c\tau/\nu_{x}$ is consistent with the down-range resolution of
bistatic synthetic aperture radar, given in Moccia and Renga (2011) after
converting between different angle conventions.

The effective energy flux for a broadband Gaussian pulse and incident Gaussian
beam is derived in Appendix B. The result is

$\displaystyle E_{f}^{\prime}$
$\displaystyle=\frac{|p_{0}|^{2}}{2\rho_{0}c}\sin\theta_{i}L_{eff}D(t^{\prime})$
(33) $\displaystyle L_{eff}$
$\displaystyle=\frac{\sqrt{\pi/2}}{\sqrt{g^{-2}+(c\tau/\nu_{x})^{-2}}}$ (34)
$\displaystyle D(t^{\prime})$
$\displaystyle=e^{-\frac{2(ct^{\prime}/\nu_{x})^{2}}{g^{2}+(c\tau/\nu_{x})^{2}}}(1+J_{1}(1+J_{2}J_{3}))\,.$
(35)

where

$\displaystyle J_{1}$
$\displaystyle=\frac{\cot\theta_{i}\nu_{x}}{\sin\theta_{i}g^{2}(2g^{-2}+(c\tau/\nu_{x})^{-2})}$
(36) $\displaystyle J_{2}$
$\displaystyle=\frac{2\left(g^{-2}+(c\tau/\nu_{x})^{-2}\right)}{\left(2g^{-2}+(c\tau/\nu_{x})^{-2}\right)}$
(37) $\displaystyle J_{3}$
$\displaystyle=\frac{\left(\frac{t^{\prime}}{\tau}\right)^{2}\left(1-\chi^{2}\right)-\left(\frac{\omega_{0}\tau}{2}\right)^{2}\left(1+\chi^{2}\right)}{\left(\frac{t^{\prime}}{\tau}\right)^{2}\left(1-\chi^{2}\right)^{2}-\left(\frac{\omega_{0}\tau}{2}\right)^{2}\left(1+\chi^{2}\right)^{2}}$
(38) $\displaystyle\chi$
$\displaystyle=\frac{\nu_{x}/(c\tau)}{\sqrt{2g^{-2}+(c\tau/\nu_{x})^{-2}}}\,.$
(39)

$L_{eff}$ is the effective ensonified length, and $D(t^{\prime})$ can be
thought of as the directivity function that captures the effect of the
incident pulse traveling throughout space and changing amplitude due to the
incident beam. Away from specular, if $g\gg c\tau/2$, then the
$L_{eff}\approx(\pi/2)^{1/2}c\tau/\nu_{x}$. Note that when
$\sin\theta_{i}g\nu_{x}/(c\tau)\gg 1$, $J_{1}$ is small, $J_{2}$ is of order
unity, and $J_{3}$ tends to 1/2.

The broadband scattering cross section is defined here in terms of an
intermediate variable, the unaveraged broadband scattering cross section as a
function of time, incident angle and scattered angle,
$q(t^{\prime},\theta_{s},\theta_{i})$,

$\displaystyle
q(t^{\prime},\theta_{s},\theta_{i})=\frac{|p_{s}(t^{\prime},\theta_{s},\theta_{i},R)|^{2}R}{|p_{0}|^{2}L_{eff}D(t^{\prime})}.$
(40)

To reduce uncertainty, $q(t^{\prime},\theta_{s},\theta_{i})$ is computed for
different roughness realizations, and results concatenated. The broadband
scattering cross section can be computed by averaging $q$ over time,
$t^{\prime}$ and random realizations, $N_{E}$.

$\displaystyle\sigma_{bb}(\theta_{s},\theta_{i})=\langle
q(t^{\prime},\theta_{s},\theta_{i})\rangle_{t^{\prime},N_{E}}.$ (41)

In Sec. VII, the broadband cross section is compared to the frequency domain
scattering cross section, $\sigma_{f}$, calculated using Eq. (14) from Thorsos
(1988). Although the scattering cross section is exclusively a frequency-
domain quantity, the subscript $f$ is included for clarity.

Figure 2: (color online) Steps to estimate broadband scattering cross section.
The frequency domain scattered pressure, and the pressure weighted by the
source spectrum are given in (a) (angle arguments to $p_{s}$ have been
omitted). In $(b)$ the raw time-domain scattered-pressure magnitude-squared is
plotted, along with a version divided by $D(t^{\prime})$ (defined in Eq.
(35)). In (c), the time-domain unaveraged broadband scattering cross section
(defined in Eq. (40)) is plotted for five different pulse lengths. The
incident and scattered grazing angles were 20∘, $g$ was 3.75 m, and $R$ was
placed conservatively in the far-field of the rough interface.

An example realization of the scattered pressure in the frequency domain is
plotted in Fig. 2(a). The frequency domain pressure is plotted as the raw
scattered pressure, and also weighted by the amplitude spectrum, $S(f)$. The
time-domain pressure-squared after weighting by $S(f)$ and an inverse Fourier
transform, is plotted in Fig. 2(b). The time-domain squared pressure magnitude
contains fluctuations with a characteristic time-scale close to the pulse
length, as well as deterministic changes due to the incident beam used in the
Helmholtz integral calculations. One example of
$q(t^{\prime},20^{\circ},20^{\circ})$ is plotted in Fig. 2(c) for five
different pulse resolutions.

The mapping $x=-ct^{\prime}/\nu_{x}$, can be used to convert the time series
of the scattered pressure to the mean plane. This mapping gives the location
of the center of the incident pulse on the mean plane, as is sometimes
performed for imaging sonars. If single-scattering is assumed, then the time-
domain scattered field can be mapped to the horizontal position of the rough
interface. If the surfaces are very rough and multiple scattering is present,
then the scattering will occur from locations other than
$x=-ct^{\prime}/\nu_{x}$.

The scintillation index, $SI$ is also examined. It is the variance of the
scattered intensity divided by the square of the mean scattered intensity
(Ishimaru, 1978, p. 437). Since $SI$ is invariant under a multiplication of
the intensity by a constant, the unaveraged broadband cross section, $q$, may
be used instead of the scattered intensity, so that

$\displaystyle SI=\frac{\langle q^{2}\rangle-\langle q\rangle^{2}}{\langle
q\rangle^{2}}\,.$ (42)

The scintillation index characterizes the fluctuations in the scattered field.
If $SI=1$, then the magnitude of the complex pressure (known as the envelope)
has a Rayleigh distribution and its real and imaginary components are
Gaussian. If $SI>1$, then the pdf of the scattered field is heavy-tailed,
which means that there is a higher probability of occurrence of high-amplitude
events compared to the Rayleigh distribution.

##  VI Parameters for Numerical Experiments

###  VI.1 Signal Parameters

The objective of this work is to study the resolution (or bandwidth)
dependence of the scattered field. These experiments covered the resolutions
typically used in narrowband scattering experiments Jackson _et al._ (1986a);
Williams _et al._ (2002), with the resolution cell on the order of 10 or more
wavelengths, down to a value of one wavelength, which is on the order of what
is achievable by modern SAS systems. Specific values of $\Delta
X/\lambda_{0}=(1,2,4,8,16)$ were used. The proportional spatial resolutions
correspond to temporal resolutions $\Delta\tau f_{0}=(2,4,8,16,32)$ at small
grazing angles, since $\Delta X=c\Delta\tau/(2\cos\theta_{i})$ for
backscattering. In all cases, the resolution is defined for
$\theta_{i}=\theta_{s}=0$. At larger incident and scattered grazing angles,
the resolution will be somewhat larger than the values listed in the figures
and tables presented below.

High-frequency acoustic imaging systems provided the motivation for this work,
and thus the simulations used a center frequency of 100 kHz. However, the
parameters of the simulation were specified in a non-dimensional fashion. As
long as every dimensional quantity is scaled properly, results of these
simulations should be valid for lower frequencies with larger roughness
parameters, and longer surface lengths. To check whether the non-dimensional
scaling was valid, one of the simulations was performed at 10 kHz as well, and
the roughness parameters were scaled accordingly. A sound speed of 1500 m/s
was used for all simulations.

###  VI.2 Roughness parameters

Roughness parameters were specified using two dimensional constants, the
center wavenumber, $k_{0}$ (through a combination of $f_{0}$ and $c$) and
spectral strength, $w$. Three spectral exponents were used,
$\gamma=(1.5,2,2.5)$, since this parameter has been observed to vary for
measured seafloor roughness (Jackson and Richardson, 2007, Ch 6). A spectral
exponent of 2 was used for both 100 kHz and 10 kHz. The $\gamma=1.5$ and
$\gamma=2.5$ simulations used a center frequency of 100 kHz only.

It now remains to specify the spectral strength. For $\gamma=2$, spectral
strengths of $w_{2}=(1\times 10^{-6},1\times 10^{-5},2\times 10^{-5},3\times
10^{-5},4\times 10^{-5})$ m were used, where a subscript on the spectral
strength denotes that it is used for a specific value of the spectral
exponent. These values resulted in $SI\approx 1$ for the smallest $w_{2}$, and
$SI>1$ for larger values. These values are much smaller than the measurements
with exponent close to 2 summarized in Table 6.1 of Jackson and Richardson
(2007). The same is true for spectral strength for other $\gamma$ in Table 1.

Given our interest in investigating the local tilting hypothesis, the spectral
strength for other values of $\gamma$ were specified such that they resulted
in equal rms slope when the effective inner scale was held constant at
$\ell^{\prime}=A\lambda_{0}\approx 2\lambda_{0}$, where $A$ is defined in Eq.
(19). The factor of $A$ is included since $(A\Delta X)$ is the length scale
relevant for computing rms slope and height with an effective upper wavenumber
limit. With this requirement, values of $w$ for other $\gamma$ were computed
using

$\displaystyle
w_{\gamma}=w_{2}(3-\gamma)\frac{k_{u}^{\prime}-k_{l}^{\prime}}{k_{u}^{\prime
3-\gamma}-k_{l}^{\prime 3-\gamma}}.$ (43)

The effective upper cutoff was set to $k_{u}^{\prime}=k_{0}/(2A)$, since
$\ell^{\prime}=A\lambda_{0}$.

The true inner scale was $\ell_{0}=\lambda_{0}/6$ to ensure that the Bragg
condition was satisfied for all wavelengths in the broadband simulation with
significant support. This criterion ensured that $1.5f_{0}$ was the largest
frequency that satisfied the Bragg condition. The pulse with the largest
bandwidth has a spectrum that is 67 dB down at $1.5f_{0}$ compared to the the
peak value at $f_{0}$. The outer scale was $L_{0}=112.5\lambda_{0}=L/10$ to
ensure that a broad range of scales were included in the power law spectrum,
but the surface length was sufficiently larger than the outer scale.

Roughness parameters for the 100 kHz simulations are summarized in Table 1.
Root-mean-square slope is given in degrees for different values of
$\ell^{\prime}$, the effective inner scale. With
$\ell^{\prime}=\ell_{0}=\lambda_{0}/6$, the true rms slope of the rough
interface is given. Values of $s_{\ell^{\prime}}$ for $\ell^{\prime}$ between
$A\lambda_{0}$ and $16A\lambda_{0}$ are also given, each of which corresponds
to the acoustic resolutions used in the numerical simulations. For the
$\gamma=2$ case, the true rms slope varied between $2.9^{\circ}$ and
$18^{\circ}$, depending on the spectral strength. For the $\gamma=1.5$ case
the true rms slope varied between $5.3^{\circ}$ and $31^{\circ}$, and for the
$\gamma=2.5$ case it varied between $1.6^{\circ}$ and $10^{\circ}$. The
effective rms slopes were much smaller than the true values. The large values
of true rms slope result because this quantity increases without bound as the
inner scale tends to zero. For the 10 kHz simulations, the same $s_{\delta x}$
and $kh$ were chosen. This condition can be satisfied if the spectral
strength, surface length, inner scale, and outer scale, for $100$ kHz and
$\gamma=2$ are all multiplied by 10. For other values of $\gamma$, a different
scaling must be used.

$\gamma$ $w$ $s_{\lambda_{0}/6}$ $s_{A\lambda_{0}}$ $s_{2A\lambda_{0}}$
$s_{4A\lambda_{0}}$ $s_{8A\lambda_{0}}$ $s_{16A\lambda_{0}}$ $k_{0}h$ - m3-γ ∘
∘ ∘ ∘ ∘ ∘ - 2 1.00$\times 10^{-6}$ 2.87 0.79 0.55 0.37 0.24 0.13 0.31
1.00$\times 10^{-5}$ 9.00 2.49 1.73 1.17 0.75 0.40 0.97 2.00$\times 10^{-5}$
12.62 3.52 2.44 1.66 1.06 0.56 1.37 3.00$\times 10^{-5}$ 15.33 4.31 2.99 2.03
1.30 0.69 1.68 4.00$\times 10^{-5}$ 17.57 4.97 3.45 2.34 1.50 0.80 1.94 1.5
1.47$\times 10^{-7}$ 5.34 0.79 0.47 0.27 0.15 0.07 0.22 1.47$\times 10^{-6}$
16.44 2.49 1.47 0.86 0.48 0.23 0.71 2.93$\times 10^{-6}$ 22.65 3.52 2.08 1.21
0.68 0.32 1.00 4.40$\times 10^{-6}$ 27.07 4.31 2.55 1.49 0.83 0.39 1.23
5.86$\times 10^{-6}$ 30.54 4.97 2.94 1.72 0.96 0.45 1.42 2.5 5.90$\times
10^{-6}$ 1.61 0.79 0.63 0.48 0.35 0.21 0.44 5.92$\times 10^{-5}$ 5.09 2.49
1.99 1.53 1.11 0.65 1.39 1.18$\times 10^{-4}$ 7.18 3.52 2.81 2.17 1.57 0.93
1.96 1.78$\times 10^{-4}$ 8.77 4.31 3.44 2.66 1.92 1.13 2.41 2.37$\times
10^{-4}$ 10.11 4.97 3.97 3.07 2.21 1.31 2.78

Table 1: Roughness parameters used in the simulations. All parameters are listed for $f_{0}=100$kHz and $c=1500$ m/s. Simulations at 10 kHz use the same dimensionless mean square parameters. Units for each parameter are given in the second row. The rms slope is reported using the upper limit computed with the sampling interval for the rough surface, as well as using the acoustic resolution. The parameter $A\approx 2.12$. $\Delta X/\lambda_{0}$ | $N_{p}$ | Err (%) | Err (dB)
---|---|---|---
1 | 4275 | 1.53 | 0.066
2 | 2137 | 2.16 | 0.093
4 | 1069 | 3.06 | 0.13
8 | 534 | 4.33 | 0.18
16 | 267 | 6.12 | 0.26
Table 2: Number of independent samples, $N_{p}$ for the different resolutions
at small grazing angles, and relative uncertainty associated with the finite
ensemble assuming intensity is exponentially distributed. Uncertainty is
reported in terms of percent and decibels. The frequency-domain simulations
used 5000 surface realizations, with an uncertainty of 1.14%, or 0.061 dB.

###  VI.3 Sampling parameters

The sampling interval, $\delta x$ was specified so that the errors caused by
discretization of the integral equation were not observable to within Monte-
Carlo error. Since very large numbers of independent samples were used, these
estimates of the broadband scattering cross section had uncertainty of about
0.07-0.3 dB. With this small uncertainty, a noticeable bias in the scattering
strength occurred when $\delta x/\lambda_{0}>1/12$, due to discretization
error. A conservative value of $\delta x/\lambda_{0}=1/15$ was therefore used
to minimize this bias.

The focus in this work is moderate to low grazing angles. The lower limit to
the grazing angles that can be reliably estimated in these numerical
simulations is set by the surface length, which in turn is set by the memory
limitations and acceptable number of CPU hours used. These latter constraints
limited us to $g=250\lambda_{0}$. At the center frequency, the relative
angular width for that value of $g$ was about 3.5% at 10∘ grazing angle. Lower
frequencies are more problematic, since, for a constant value of $g$,
decreasing the frequency will increase the angular width of the field. At the
lower 6 dB down point of the largest bandwidth signal used (equivalent to
$0.844f_{0}$), the angular width for $g=250\lambda_{0}$ was approximately 4.2∘
at a nominal angle of $10^{\circ}$ grazing angle. With these angular widths,
$10^{\circ}$ was taken to be an acceptable, if conservative, lower limit to
the grazing angles that can be reliably estimated in this work. In choosing
the surface length, $g$ was increased to $400\lambda_{0}$ without any change
in the behavior of $\sigma_{bb}$ or $SI$ above 10∘. Grazing angles less than
$10^{\circ}$ likely require fast approximate methods to solve the integral
equation, such as the fast multipole method Liu (2009), or hierarchical
matrices Hackbusch (2015).

The rough surface length, $L$ was set to $4.5g$ to allow the incident field
taper to decay sufficiently at the edges of the computational domain. Near the
edges, 2.5% of the time-domain samples were discarded on each side to remove
edge effects (e.g. the large peak at -11 ms in the blue curve in Fig. 2(b)),
so that only 95% of the full time series available was used.

In the broadband simulations, four roughness realizations were used, and 5000
realizations used in the frequency domain. For the broadband case, the time-
domain response for each angle contains many independent samples, $N_{p}$
(including multiple realizations). At small grazing angles for all ensembles,
$N_{p}$ is $4\times 0.95L/(\Delta X)$, where the factor of $0.95$ results from
discarding samples near the edges of the surface. Assuming Gaussian pressure
statistics (or exponentially distributed intensity statistics), one can define
a standard error $Err_{\mathrm{std}}=1/\sqrt{N_{p}}$ which characterizes the
relative uncertainty of the broadband scattering cross section estimate. The
decibel representation is $Err_{\mathrm{dB}}=10\log_{10}(1+1/\sqrt{N_{p}})$.
However, the scattered intensity is manifestly not exponentially distributed
for large values of spectral strength, small pulse length, and small grazing
angles. A more realistic characterization of the uncertainty is to divide the
standard deviation of the intensity by the mean intensity and $\sqrt{N_{p}}$,
amounting to $\sqrt{SI/N_{p}}$. Since the scintillation index is shown in the
next section to vary with roughness parameters and grazing angle, it is
difficult to give an overall characterization of the uncertainty. Therefore,
uncertainty calculated using Gaussian statistics are given in Table 2, with
the smallest error of 1.5% or 0.066 dB, and largest uncertainty of 6.1% or
0.26 dB. If the uncertainty for the non-Gaussian cases are desired, then the
scintillation index plots presented in Figs. (5-7) can be used along with
$N_{p}$ to provide this quantity.

The frequencies required for the largest bandwidth simulation $\Delta\tau
f_{0}=2$, spanned approximately $0.2f_{0}$ to $1.8f_{0}$. This computational
bandwidth was about seven times the largest 3 dB bandwidth used for the
Gaussian spectrum. Using a frequency spacing of $\delta f=c/(3L)$, the number
of frequencies per simulation was approximately 6000. For the proportional
bandwidths studied in this work, these surface parameters resulted in surfaces
with $N=16,875$ points. The matrix $V$ resulting from this discretization
required 4.3 GB of memory storage for double precision complex numbers.
Simulations were performed on the Hamming high-performance computing cluster
at the Naval Postgraduate School.

##  VII Results

Figure 3: (color online) Scattering Strength comparison for $f_{0}$= 100 kHz
and $\gamma=2$. Each color represents scattering strength estimated using
different bandwidths, as indicated by the caption. The dot-dashed line
indicates the frequency-domain estimate. Five different values of spectral
strength are plotted on separate lines, with scattering strength monotonically
increasing with spectral strength (i.e. the smallest spectral strength has the
lowest scattering strength). Values for the spectral strength are given in
Table 1. Note that the random realizations are generated only with inner
cutoff scale of $\ell_{0}=\lambda_{0}/6$. The time-domain results are
difficult to distinguish, and the ratio for the time-domain and frequency
results are shown in Fig. 4. Figure 4: (color online) Scattering Strength
ratio for $f_{0}$=100 kHz and $\gamma=2$ Figure 5: (color online)
Scintillation Index for $f_{0}$=100 kHz and $\gamma=2$. Note that the vertical
axis is different for each subfigure. Figure 6: (color online) Scintillation
Index for $f_{0}=$100 kHz and $\gamma=1.5$. Note that the vertical axis is
different for each subfigure. Figure 7: (color online) Scintillation index for
$f_{0}$=100 kHz and $\gamma=2.5$. Note that the vertical axis is different for
each subfigure.

Results for backscattering strength as a function of grazing angle,
$\theta_{i}$, resolution $\Delta X$, and spectral strength, $w$, are presented
in Fig. 3 for $f_{0}=100$ kHz and $\gamma=2$. The numerical method can handle
bistatic geometries (with differing $\theta_{i}$ and $\theta_{s}$), but only
monostatic geometries are presented here. Broadband scattering strength,
$\sigma_{bb}$ as well as the frequency domain version, $\sigma_{f}$ is plotted
on the vertical axis, with grazing angle on the horizontal axis. The angles
have a lower limit of 10∘ grazing angle due to finite surface length, and an
upper limit of 70∘, due to the difficulty in estimating the broadband cross
section near vertical using broadband pulses Hefner (2015); Hellequin _et al._
(2003). Each resolution is plotted as a different color and line style, with
the frequency domain version as a black dashed-dot line. Results for different
spectral strengths are on the same figure since they are well separated from
one another, with the smallest $w$ corresponding to the lowest scattering
strength. For each value of the spectral strength, the frequency-domain
scattering cross section is indistinguishable from the broadband version.

To compare these results more closely, the broadband cross section is divided
by the frequency-domain cross section and the dB value taken. This quantity,
which is called the scattering strength difference, is plotted in Fig. 4. At
the largest spectral strength, some systematic oscillations as a function of
$\theta_{i}$ are present, but cannot be easily disentangled from the rapid
Monte-Carlo fluctuations. Other than that case, all differences appear to be
random. The standard deviation of the ratio $\sigma_{bb}/\sigma_{f}$ across
angles between 10 and 70 degrees grazing is about 4% for $\Delta
X/\lambda_{0}=1$, and about 10% for $\Delta X/\lambda_{0}=16$. Note that since
this ratio is between two random variables, the uncertainty may be higher than
the theoretical uncertainty for $\sigma_{bb}$ or $\sigma_{f}$ alone. These
numbers are relatively consistent across all spectral strengths. Larger
uncertainty for long pulses is a consequence of having fewer independent
samples per roughness realization. The uncertainty for long pulse lengths is
consistent with the theoretical uncertainty shown in Table 2, but the
uncertainty for short pulse lengths is higher by about factor of 2.5. This
increased uncertainty is likely a result of the departure from exponential
intensity statistics for short pulses. Broadband scattering strength and the
frequency domain scattering strength are indistinguishable, to within the
Monte-Carlo error of the numerical simulations. This result confirms that the
broadband scattering strength is a robust metric across different resolutions
for scattering from a power-law seafloor with spectral exponent $\gamma=2$.

The frequency domain scattering cross section exhibits power-law frequency
dependence, as seen by the nearly $f^{1}$ trend in Fig. 2(a). When
incoherently averaged in the frequency domain, weighted by $S(f)$, the
frequency dependence in the cross section for $\gamma\in[1,3]$, results in a
negligible bandwidth dependence, under 4 percent. If $\sigma_{f}$ exhibited
very strong frequency dependence, e.g. with sharp peaks seen in scattering
from layered seafloors (e.g. Jackson and Ivakin (1998); Jackson and Olson
(2020)), then more significant bandwidth dependence may be observable.

The scintillation index for this case is plotted in Fig. 5 for broadband
signals as well as the single frequency case. Each spectral strength is
plotted in its own subfigure, and each resolution has its own line within the
subfigure. The vertical axis is $SI$, and the horizontal axis is grazing angle
in degrees. The narrowband result is approximately unity for the entire
angular domain shown, for all spectral strengths. This is the expected result
if the central-limit theorem is employed. For the broadband signals, there is
a profound dependence on resolution, with the scintillation index increasing
as the resolution cell becomes small. This behavior can be seen in the example
realization shown in Fig. 2(c), in which the intensity peaks become higher as
$\Delta X$ becomes small. Additionally, holding resolution constant, the
scintillation index increases as the grazing angle becomes small,
monotonically for this case. For most broadband cases, $SI$ asymptotically
approaches unity as grazing increases to its upper limit.

However, for the highest resolution cases, $\Delta X/\lambda_{0}=(1,2)$, and
the largest spectral strengths, this high angle asymptote is greater than one,
indicating that for all angles examined here, scattered complex pressure
magnitude is non-Rayleigh. In Lyons _et al._ (2016) a K distribution was
required to describe the pdf of the scattered field at moderate grazing
angles, which agrees with this result. The Monte-Carlo fluctuations are
significantly less than the difference between the high-angle $SI$ asymptote
and unity, indicating that this is a statistically significant finding.

Lupien (1999) also observed non-Rayleigh scattering for broadband scattering
from rough surfaces with a power-law exponent of $\gamma=3$, but statistical
tests barely rejected the Rayleigh distribution. That analysis did not remove
the effect of the Gaussian taper, so the conclusions are not comparable to the
present work.

Broadband and single-frequency scattering strength and $SI$ were also computed
for 10 kHz, $\gamma=2$, and spectral strengths that were ten times the value
in the previous section. The relative resolution, $\Delta X/\lambda_{0}$ was
held constant, but consequently the resolution $\Delta X$ was a factor of 10
larger. These parameters were chosen such that the dimensionless second-order
quantities were the same as in the previous 100 kHz simulations. This scaling
is only true for $\gamma=2$, and would be modified if $\gamma$ were a
different value. It was found that the scattering strengths, scattering
strength dB error, and the scintillation index were the same for the 10 kHz
and 100 kHz cases, to within Monte-Carlo fluctuation. This set of simulations
was performed to verify that characterizing the simulations non-dimensionally
was valid. Since plots for the 10 kHz case do not add significantly new
information, they are not shown here. These results indicate that departure
from Rayleigh statistics is not isolated to very high-frequency imaging
systems, and may occur in lower-frequency sonar systems as well, so long as
the seafloor has the appropriate roughness parameters.

The spectral exponent was changed to $\gamma=1.5$ to examine the effect of
changing the shape of the power spectrum. New values of $w$ were used, as
specified in Table 1. Again, $\sigma_{bb}$ and $\sigma_{f}$ are the same.
Scattering strength comparisons and the scattering strength ratio are not
shown. The scintillation index is plotted in Fig. 6, and is seen to depend on
angle, resolution, and spectral strength, as in the previous case. The rough
interfaces have the same effective large-scale slope for $\gamma=2$ and
$\gamma=1.5$ at $\Delta X/\lambda_{0}=1$. $SI$ qualitatively depends similarly
on spectral strength and resolution, as compared to the $\gamma=2$ case.
Quantitatively, the values of $SI$ are slightly different than the $\gamma=2$
case, but are similar.

Finally, the spectral exponent was changed to $\gamma=2.5$. The values of
spectral strength can be found in Table 1. Again, the broadband scattering
cross section, and the frequency-domain version computed at the center-
frequency were the same to within the Monte-Carlo error of the simulations,
and are not shown here. The scintillation index is plotted in Fig. 7, and
again depends on resolution, spectral strength and grazing angle. As the
spectral strength is increased in the same proportions as the earlier plots
(the second through fifth spectral strengths are 10, 20, 30, and 40 times the
smallest spectral strength respectively), the scintillation index increases
much more rapidly than either the $\gamma=2$ and $\gamma=1.5$ cases. For the
four largest values of spectral strength, the $SI$ is elevated at the higher
angles as well. The smallest resolution cases also have elevated $SI$ for the
entire angular domain, and the four highest spectral strength values.

##  VIII Discussion

In the results presented in Sec. VII, it was shown that the broadband
scattering strength is indistinguishable from frequency-domain scattering
strength, and is independent of pulse length. This conclusion is not
surprising. However, as the pulse length changes, the properties of the
ensemble used to estimate scattering strength changes as well. It is
encouraging to see that although the ensemble is changing with respect to
resolution (i.e the rough patch within a resolution cell is different for each
resolution), $\sigma_{bb}$ is invariant to pulse length. It is expected that
this result holds for 3D environments as well if the roughness is isotropic,
although numerical tests using this geometry are needed. Based on these
results, high-resolution systems should be able to reliably estimate
scattering strength, and this work also confirms that it may be a stable
quantity to use for seafloor remote sensing. However, for highly anisotropic,
non-stationary scenarios, such as those studied by Olson _et al._ (2016,
2019); Lyons _et al._ (2010), the measured scattering strength may depend on
pulse length.

The scintillation index (also called structure Wang and Bovik (2002),
lacunarity Williams (2015), or contrast Marston and Plotnick (2015)) was shown
to be is highly dependent on all the parameters studied: resolution, grazing
angle, spectral strength, and spectral exponent. For moderate to low grazing
angles, $SI$ monotonically increases as grazing angle decreases, resolution
cell decreases, and spectral strength increases. SI is similar for
$\gamma=1.5$ and $\gamma=2$, but is larger for $\gamma=2.5$, for the values of
$w$ chosen for these numerical experiments. Contrary to scattering strength,
$SI$, and therefore the scattering process in general, is fundamentally
different in the frequency and time domains for broadband pulses.

In Lyons _et al._ (2016), it was hypothesized that the physical cause of
heavy-tailed statistics in high resolution sonar imagery was local tilting of
the seafloor due to roughness wavelengths larger than the acoustic resolution.
The scattered pressure amplitude (envelope) was modeled as a product between a
random variable due to sub-resolution roughness, and a random variable that
took into account the effect of tilting by longer wavelengths – the small-
scale scattering strength evaluated at the nominal grazing angle modified by
the local slope. Local tilting (due to large-wavelength roughness components)
modulates the Rayleigh-distributed field (due to small-wavelength components)
and causes the $SI$ to be greater than unity. An rms slope with an upper
cutoff related to the acoustic resolution was used as the input parameter to a
simplified version of the composite roughness model McDaniel and Gorman
(1983), which was then used to compute the scintillation index.

Here, the role of slopes at scales at or larger than the pulse length on the
intensity fluctuations is investigated. A representative function is used to
model the effect of slope modulation that is inspired by the composite
roughness model. The composite roughness approximation McDaniel and Gorman
(1983) uses separate scattering models on the small and large-scale surfaces.
The validity of the composite model is subject to the validity of the
perturbation approximation applied to the small-scale surface, and the
Kirchhoff approximation applied to the large-scale surface. The perturbation
approximation is valid when $k_{0}h_{s}$ is small Thorsos and Jackson (1989),
where $h_{s}$ is the rms roughness of the small-scale roughness. This
parameter can be found by integrating the roughness spectrum over wavenumbers
with magnitudes between $\pi/(A\Delta X)$ and $\pi/\ell_{0}$. In this work,
$k_{0}h_{s}$ was less than 2.0 for all simulations. Since perturbation theory
must be applicable on the small-scale surface, $k_{0}h_{s}<1$ was required.
This is possible if the highest spectral strength for each value of $\gamma$
is ignored, and the restriction $\Delta X/\lambda_{0}\leq 8$ is made. We found
that scattering strength was well modeled by perturbation theory for a few
cases where $k_{0}h>1$ (including all the wavenumber components), but exclude
these cases since the validity region of this model for general power-law
roughness spectra has nost been previously studied.

In McDaniel and Gorman (1983); Jackson _et al._ (1986b), the validity of the
Kirchhoff approximation was stated to be a function of the rms radius of
curvature of the large-scale surface. However, Thorsos (1988), found that the
radius of curvature is not important for the Kirchhoff approximation, but that
the characteristic length of the surface, $L_{c}$, must be large compared to
the wavelength. In the separation of scales performed here, the characteristic
length of the large scale surface is approximately equal to $\Delta X$, which
is always $\lambda_{0}$ or greater. Therefore $k_{0}L_{c}\geq 2\pi$, which is
typically acceptable for the Kirchhoff approximation Thorsos (1988). At low
grazing angles, which are most relevant here, the Kirchhoff approximation
fails due to the presence of multiple scattering, which is important when
$\theta_{i}$ and $\theta_{s}$ are each less than or equal to $2\epsilon$,
twice the large-scale rms slope angle Thorsos (1988); Thorsos and Jackson
(1991). Shadowing is another source of inaccuracy, which occurs when
$\theta_{i}$ and $\theta_{s}$ are less than or equal to $\epsilon$ Thorsos
(1988). The large-scale slope angles can be seen in Table 1. For the rms slope
angles that are based on the acoustic resolution as the upper cutoff, the
global maximum is about 5∘. Since this number is comparable to the smallest
grazing angles examined here, both multiple scattering and shadowing may be
present in these cases.

To explore the local tilting hypothesis, the relationship between intensity as
a function of space $I_{\Delta X}(x)$, by which is $I(t^{\prime})$ for a given
resolution mapped to $x$ through $x=-ct^{\prime}/(2\cos\theta)$ for
backscattering, with the large-scale surface slope at position $x$ is
investigated. $I_{\Delta X}(x)$ should decrease if the large-scale slope at
$x$ is positive, and increase if negative. Consequently there should be a
statistical correlation between $I_{\Delta X}(x)$, the intensity at a given
resolution at position $x$, and $-s_{\Delta X}(x)$, the negative of the slope
field low-pass filtered to remove wavelengths shorter than $\Delta X$. More
specifically, the scattered intensity from the integral equations can be
compared to the intensity produced by the effect of tilting in the composite
roughness approximation. A simplified version of this approximation is used,
which is explained in the following paragraphs.

The composite model requires two types of averages, one over the small-scale
and one over the large-scale. The large-scale average is a simple average of
the small-scale scattering cross section over the pdf of the large-scale
slopes. The small scale average results in the scattering cross section from
the small-scale roughness. Using the results of McDaniel and Gorman (1983),
the composite model may be written as (after converting between differing
conventions of the roughness power spectrum, and adapting to 1D roughness used
here),

$\displaystyle\sigma_{CR}$
$\displaystyle=4k^{3}\langle\sin(\theta_{i}-\epsilon_{L})^{4}W_{s}(\Delta
K_{mod})\rangle_{L}$ (44) $\displaystyle\Delta K_{mod}$ $\displaystyle=\Delta
K+\frac{df_{L}}{dx}\Delta k_{z}$ (45)

where $\langle\cdot\rangle_{L}$ denotes averaging over the large-scale slope,
and $\epsilon_{L}=\tan^{-1}(df_{L}/dx)$ is the large-scale slope angle at
position $x$. The power spectrum of the small-scale roughness is denoted
$W_{s}(\Delta K_{mod})$, where the modified Bragg wavenumber, $\Delta K_{mod}$
is specialized to the backscattering case and takes into account the effect of
tilting on the Bragg components. The difference between the vertical component
of the scattered and incident vertical wavenumber is $\Delta
k_{z}=2k_{0}\sin\theta_{i}$, and the horizontal wavenumber difference is
$\Delta K=2k_{0}\cos\theta_{i}$. Although not present in most applications of
the composite roughness model (e.g. Kur’yanov (1963); Bachmann (1973);
McDaniel and Gorman (1983)), the modulation of the Bragg wavenumber due to
large-scale slopes may be significant. The form of the Bragg modulations is
given in an intermediate result (Eq. (22)) of McDaniel and Gorman (1983), and
adapted to the present notation. This model was formulated in the frequency
domain, and is not directly applicable to these broadband simulations.
However, in the broadband case, the effect of tilting would be preserved, even
if the scattered power from the small-scale surface is altered.

Because of this simplification and the broadband nature of the simulations,
this model is called the slope modulation model (SMM) instead of the two-
scale, or composite roughness model. Terms that are constant as a function of
the large-scale slope were ignored in the SMM. Averaging over the large-scale
slopes was not performed in the model, to compare fluctuations produced by
tilting to the fluctuations from the integral equation results for an given
ensemble. The intensity fluctuations caused by tilting can be written as,

$\displaystyle\begin{split}I^{SMM}_{{\Delta
X}}(x)=\sin^{4}&\left(\theta_{i}-\tan^{-1}\left[\frac{df_{\Delta
X}}{dx}(x)\right]\right)\\\ \times&W(\Delta K+\Delta k_{z}\,\frac{df_{\Delta
X}}{dx}(x))\end{split}$ (46)

where the subscript $\Delta X$ on the intensity and slope explicitly denotes
that a Gaussian function with 3 dB width $\Delta X$ was used to filter the
large-scale slopes. The filtered slope is defined by

$\displaystyle\frac{df_{\Delta
X_{surf}}}{dx}(x)=\int\limits_{-\infty}^{\infty}\frac{df}{dx}(x^{\prime}-x)e^{-x^{\prime
2}\ln(4)/(\Delta X_{surf})^{2}}\mathrm{d}x^{\prime}\,.$ (47)

Eq. (46) has been averaged over the small-scale roughness, but not the large-
scale slopes. It is important to note that in the integral equation
simulations, it is impossible to perform this partial averaging. Thus, taking
the variance of Eq. (46) does not result in the scintillation index of the
scattered pressure, but only the variance of slope-induced fluctuations.
Taking the limit of zero large-scale rms slope in the previous equation
results in no slope-induced fluctuations, and fluctuations in the total
scattered field would have a Rayleigh distribution, with a scintillation index
of unity. It is strongly emphasized that this model is only a representative
function of the effect of slopes at the scale of the pulse resolution, and is
not an adequate model for the scattered intensity. This inadequacy is due to
the broadband nature of these simulations, and the formulation of the
composite model in the frequency domain. The integral equation results are
considered to be the accurate “ground truth” here.

Although not shown in (46), $I^{SMM}_{{\Delta X}}(x)=0$ if the argument to the
sine function is less than zero, to include a rudimentary form of shadowing.
This form of shadowing has a small effect, and the results presented below are
essentially unchanged if this simple form of shadowing is left out. Non-local
shadowing may have a significant effect, but is not examined here. The
simplifications used here are employed because a rigorous formulation of the
composite roughness approximation for very broadband signals is not available
in the literature. Development of such a model, and comparison to results
using the broadband integral equation technique developed here are both
fruitful areas for future work.

Similarity between the intensity from the numerical solutions and model
predictions can be quantified in a crude but straightforward manner using the
Pearson product-moment correlation coefficient, $\rho$, defined for random
variables $U$ and $V$.

$\displaystyle\rho\left(U,V\right)=\frac{\langle\left(U-\langle
U\rangle\right)\left(V-\langle
V\rangle\right)\rangle}{\sqrt{\langle\left(U-\langle
U\rangle\right)^{2}\rangle\langle\left(V-\langle V\rangle\right)^{2}\rangle}}$
(48)

This coefficient quantifies the linear variation between a dependent and an
independent variable. The correlation coefficient is used here because a
rigorous application of the composite roughness approximation was not used,
and absolute intensity cannot be compared directly.

The local tilting hypothesis is tested by forming the correlation coefficient
between the scattered intensity from numerical simulations, and
$I^{SMM}_{{\Delta X}}(X)$. Both the acoustic resolution and the surface filter
size are varied. The surface filter scale that maximizes the correlation for
each acoustic resolution is estimated, and this process is repeated for all
acoustic resolutions. The acoustic resolution is $\Delta X$, and the surface
filter size is $\Delta X_{\textrm{surf}}$. The $\Delta X_{\textrm{surf}}$ that
maximizes $\rho$ is denoted $\Delta X_{\textrm{max}}$. If $\Delta
X_{\textrm{max}}$ varies in proportion to $\Delta X$, then it may be concluded
1) that slope modulation is responsible in part for the intensity
fluctuations, and 2) that slopes at (or larger than) the scale of the acoustic
resolution are responsible in part for the fluctuations. Note that the surface
filters used here are zero-phase, acausal filters (Eq. (47)), and are applied
in the same way as the Fourier synthesis used to obtain the time-domain
scattered pressure, Eq. (23).

Figure 8: Correlation coefficient as a function of the acoustic resolution,
$\Delta X$, and the surface filter size, $\Delta X_{surf}$. The Gaussian pulse
shape has been used as the surface filter here, as specified by Eq. 47. Figure
9: (color online) Surface slope filter scale that maximizes the product moment
cross correlation coefficient for each acoustic resolution. Grazing angle has
been held constant at $20^{\circ}$, and each line represents a different
spectral strength. Each subplot contains a different spectral exponent with
(a)$\gamma=1.5$, (b)$\gamma=2$, and (c),$\gamma=2.5$. Dashed black lines have
a slope of unity and intercept of zero for reference. Note that the Gaussian
pulse shape, Eq. (47) has been used as the surface filter.

The correlation coefficient for parameters $\gamma=2$,
$\theta_{i}=20^{\circ}$, and $w=3\times 10^{-5}$m is plotted in Fig. 8.
$\Delta X$ is on the horizontal axis, $\Delta X_{\textrm{surf}}$ is on the
vertical axis, and $\rho$ is denoted by grayscale. Holding $\Delta X$
constant, there is a distinct peak in $\rho$ as a function of $\Delta
X_{\textrm{surf}}$. The peak value of $\rho$ increases as the acoustic
resolution becomes small, indicating that the role of slope modulation is
greater for smaller resolution cells. Additionally, the peak location in
$\Delta X_{\textrm{surf}}$ varies with $\Delta X$, indicating the local
tilting hypothesis may be correct. This plot has a similar structure for other
$\theta_{i}$, $\gamma$, and $w$.

$\Delta X_{\textrm{max}}$ as a function of $\Delta X$ is plotted in Fig. 9.
$\theta_{i}$ is constant at $20^{\circ}$ grazing angle, and each spectral
strength is plotted as its own line. The different values of $\gamma$ appear
in subfigures. For each spectral strength and exponent, $\Delta
X_{\textrm{surf}}$ varies monotonically with the acoustic resolution, except
for the smallest spectral strength for $\gamma=2.5$. A line with unit slope
and zero intercept is also plotted for reference. For many cases, the slope of
these lines is approximately unity when the acoustic resolution is small. Some
of the lines retain that slope when the acoustic resolution is larger, but
others taper off and have a smaller slope. For example the smallest roughness
cases in $\gamma=2$ and $\gamma=2.5$ curve downards to an approximately
constant function at large $\Delta X/\lambda_{0}$. For these two curves in
particular, the scintillation index was approximately unity, and therefore
slope modulation was not required to explain their behavior. Other cases where
the slope of these lines becomes small for large $\Delta X/\lambda_{0}$, but
result in in $SI>1$, indicate that slopes at scales smaller than the pulse
resolution are the most important. Future work modeling broadband scattering
from power law surfaces is needed to explain this behavior.

The monotonic dependence of $\Delta X_{\textrm{max}}$ on $\Delta X$ indicates
that the large-scale slope near the acoustic resolution accounts for a
significant part of the intensity fluctuations for short pulses. The unit
slope for some of these cases indicates that the 3 dB width of the point
spread function is sometimes the appropriate scale for calculating the large-
scale slope, and is a good estimate for all cases. For longer pulses, the
lines in Fig. 9 have slopes generally less than unity, meaning that the
horizontal roughness components that are important for slope modulation are
not exactly the same as the 3 dB width of the incident intensity, but are
proportionally slightly smaller. This departure from unit slope may indicate
that the slope modulation model is not sufficient to explain the data, and
other effects, such as multiple scattering or shadowing may be present. The
same analysis was performed for the large-scale surface heights, and no
systematic trends were found. Therefore, the large-scale rough interface at a
single point is uncorrelated with the scattered intensity mapped to that
point.

To understand what kind of intensity fluctuations can be caused by this model
when for the large-scale rms slope values in these numerical simulation, the
model for $SI$ from Lyons _et al._ (2016) was applied to the case of
$w=3\times 10^{-5}$ m, $\gamma=2$, and $\Delta X=2\lambda_{0}$. This model
accounts for fluctuations due to large-scale tilting, as well as Rayleigh-
distributed scattering from the small-scale roughness to produce the
scintillation index. In the model of Lyons _et al._ (2016), $SI>1$ if the
large-scale slope variance is greater than 0, and the small-scale scattering
cross section is neither constant nor linear in $\theta_{i}$. The rms slope
was taken from Table 1, and was $2.99^{\circ}\approx 3^{\circ}$
($s_{2A\lambda_{0}}$ in Table 1). A model-data comparison is shown in Fig. 10.
The model follows the same trends as the integral equations. It is slightly
higher than 1 at high grazing angles, and increases as the grazing angle
becomes small. However, the model $SI$ increases less rapidly than the
integral equation simulation at low grazing angles. Given the evidence from
Figs. 8, and 9, as well as evidence that the composite roughness model
produces non-Rayleigh behavior that is similar to the integral equations it
may be concluded that the tilting hypothesis is responsible in part for the
intensity fluctuations. Note that both the $\sin^{4}$ term and Bragg
modulation term from Eq. (46) were included to account for fluctuations due to
slopes. Neglecting the Bragg term decreases $SI$ by about 5% at 10∘, and about
20% at 70∘ grazing angle. Including modifications to the Bragg wavenumber is
responsible for the good agreement at high and moderate grazing angles.

Since the model of Lyons _et al._ (2016) underestimates $SI$, especially for
small grazing angles, the tilting hypothesis cannot be fully responsible for
heavy-tailed scattering seen here. We suspect multiple scattering and
shadowing may be responsible for the model-simulation mismatch at low angles,
since the large-scale rms slope is 3∘ for the surface in Fig. 10, and is
comparable to the smaller range of grazing angles Thorsos (1988); Thorsos and
Jackson (1991). Shadowing may increase $SI$, since it increases the degree of
intensity fluctuations. Multiple scattering may move $SI$ towards unity, since
scattering from multiple locations on the rough interface may increase the
effective ensonified length, and contribute an additive Gaussian component to
the scattered pressure.

Figure 10: Comparison of the $SI$ from the integral equation results with the
model of Lyons _et al._ (2016), for the parameters $w=3\times 10^{-5}$m,
$\gamma=2$, and $\Delta X=2\lambda_{0}$.

A physically accurate theoretical model that includes all of these effects is
evidently required to predict the scintillation index at low grazing angles.
Such a model is, at this time, not available, and is a fruitful opportunity
for future research. The independence of the broadband scattering cross
section on resolution established here imposes a useful constraint. Any
theoretical model that breaks apart the solution of the exact integral
equation into physically interpretable phenomena (such as tilting, shadowing,
or multiple scattering) must also satisfy this constraint.

##  IX Conclusion

In this work Fourier synthesis combined with numerical solution of the
Helmholtz integral equation was used to analyze the scattered field in terms
of the broadband scattering cross section and scintillation index. The
dependence of these two quantities on acoustic resolution was examined to
understand the effects that contemporary high-resolution acoustic imaging
systems have on quantitative measurements of seafloor scattering. For power-
law surfaces, the scattering strength is independent of pulse length, which
indicates it is a stable quantity to use across measurement systems with
different geometries (at least for the kinds of rough surfaces examined here).
The scintillation index depends strongly on pulse length, which indicates that
the scattering processes is fundamentally different in the time and frequency
domains, and that further research is needed to understand or predict
intensity fluctuations in high resolution broadband sonar systems.

Although simulations were performed in two dimensions, these results may be
present in three dimensions as well, if the roughness spectrum is assumed to
be isotropic. The exact values of the scintillation index will be different
for 3D environments, as the rms slope is calculated differently, and out of
plane effects may be important.

Heavy-tailed, or non-Rayleigh scattering is commonly observed in scattering
measurements and is usually attributed to non-stationary, or patchy
environments Abraham and Lyons (2004); Lyons _et al._ (2009). Heavy-tailed
statistics have been observed in seemingly homogeneous seafloors by Lyons _et
al._ (2016), and these numerical simulations have verified that statistically
homogeneous surfaces can produce heavy-tailed statistics when interrogated by
a broadband high-resolution system. The slope modulation model used in Lyons
_et al._ (2016) was investigated, and found that it is in part responsible for
intensity fluctuations. Other sources of fluctuations, such as shadowing and
multiple scattering, were found to be necessary to accurately model the
scintillation index. Further research is required to better understand this
aspect of seafloor scattering.

A few consequences are noted for heavy-tailed statistics arising in high-
resolution systems in homogeneous roughness environments. Heavy-tailed
statistics are a significant source of false alarms in acoustic target
detection systems. Since scintillation index increases at low angles, long
range systems may suffer from decreased performance. The benefit of high
resolution systems, more pixels per target, may be offset due to the increased
false alarms.

There may be some benefits of resolution dependence of the pdf of the
scattered field. Some autofocus algorithms for synthetic aperture systems
(e.g. Blacknell _et al._ (1992); Callow (2003); Marston and Plotnick (2015)),
use the scintillation index, or contrast as their cost function. If the
acoustic field is entirely due to point scatterers, as is commonly assumed
Brown _et al._ (2017), $SI$ will be unity for all resolutions. If an autofocus
algorithm is applied, and the point spread function of the imaging algorithm
becomes smaller, then the field will still be Rayleigh and contrast will not
increase. Therefore, an autofocus algorithm based on $SI$ or lacunarity will
not be sensitive to the focus, unless there are discrete scatterers in the
scene. However, as shown here, $SI$ is a strong function of resolution for a
statistically homogeneous power law surfaces, especially at low grazing
angles. Improving the focus at low grazing angles will lead to an increase in
$SI$, and thus the autofocus algorithm will be more sensitive to the actual
degree of focus for seemingly featureless seafloors.

Median filters are often used to “remove” the speckle or intensity
fluctuations from acoustic or electromagnetic images before use in remote
sensing or target detection algorithms (e.g. Kwon and Nasrabadi (2005);
Williams (2015, 2018); Galusha _et al._ (2018) and references therein).
Although the broadband scattering cross section, which uses the arithmetic
mean of the intensity, is insensitive to resolution, the median is not, since
the median is highly dependent on the probability density function. Thus
either the broadband cross section can be used as a resolution independent
quantity that has a large variance, or its variance can be reduced with the
consequence that the pixel intensity will no longer be directly related to
broadband scattering strength.

###### Acknowledgements.

The authors thank the US Office of Naval Research for financial support of
this work through grants N00014-13-1-0056, N00014-18-WX-00769,
N00014-19-WX-00427, and N00014-16-1-2268, as well as the Research Initiation
Program of the Naval Postgraduate School. Numerical simulations were performed
on the Hamming high-performance computing cluster at Naval Postgraduate
School. The authors thank the anonymous reviewers for their comments on
previous versions of the manuscript.

## Appendix A

In this appendix, the steps are given to generate random realizations of a
given roughness power spectrum. Here, the method of Thorsos (1988) is
followed, specifically Eqs. (17-18). The spatial frequency is $u=K/(2\pi)$.
The sampling of this vector is set by the surface length, $L$, and sampling
interval, $\delta x$. The number of points, $N$, is $L/\delta x$, and the
spatial frequency spacing is $\delta u=1/L$. A spatial frequency vector is
specified, $u_{m}=m\delta u$ for $m\in[-N/2+1,N/2]$. Then, a sampled power
spectrum vector is created, $W_{m}=W(2\pi u_{m})$. To create a random
realization from this power spectrum, a sampled version of a randomized
complex amplitude spectrum, $F$, is defined by

$\displaystyle F_{m}=\sqrt{2\pi W_{j}L}\mu_{m}$ (49)

where $\mu_{m}$ are elements of a complex random vector with unit variance. To
generate a real vector of heights, $f_{n}$, the amplitude spectrum $F$ must
have conjugate symmetry. This property is achieved by setting

$\displaystyle\mu_{m}=\begin{cases}\frac{N(0,1)+iN(0,1)}{\sqrt{2}}&m\in[1,N/2-1]\\\
\mu_{-m}^{\ast}&m<0\\\ N(0,1)&m=0,N/2\end{cases}\,,$ (50)

where $N(0,1)$ is an independent random draw from the normal distribution with
zero mean, and unit variance. The first case creates random vector for
positive spatial frequencies, the second case enforces conjugate symmetry. The
third case stipulates that points in the sampled spatial frequency vector that
do not have symmetric pairs (the point at zero frequency, and the single
Nyquist point) are real. This method is valid for even values of $N$. For odd
$N$, $m\in[-(N-1)/2,(N-1)/2]$, and the lower Nyquist point, $m=-(N-1)/2$ must
be the complex conjugate of the point $M=(N-1)/2$, and only the case $m=0$
should be real.

The rough interface $f(x)$ is found through the inverse discrete Fourier
transform. For even $N$, this is

$\displaystyle f_{n}=f(n\delta x)=\delta
u\sum\limits_{m=-N/2+1}^{N/2}F_{m}e^{i2\pi\frac{mn}{N}}\,.$ (51)

For odd $N$, the sum is the same, except that the limits are from $m=-(N-1)/2$
to $(N-1)/2$. In practice, this sum is approximated using the fast Fourier
transform algorithm. The MATLAB computer language was used to generate the
rough surfaces, which itself uses the FFTW library (fastest Fourier transform
in the west) Frigo and Johnson (2005). This library requires that the input
spatial frequency spectrum be contained in a vector starting with the $m=0$
point, followed by the points with $m>0$, followed by the points with $m<0$.
The FFTW library also divides by $N$ when performing the inverse FFT, so the
output must be multiplied by $N$, and by $\delta u$ to form the sum in Eq.
(51).

## Appendix B

Here, the effective energy flux is calculated for the incident field used in
this work. Effects of the $w_{t}$ term in the extended Gaussian beam are not
included, since analytically tractable results are not available for that
case. Although the full extended Gaussian beam was used as the incident field
in the numerical simulations, the error introduced by omitting it in
calculating the energy flux was insignificant, as the parameter
$(k_{0}g\sin\theta_{i})^{-1}$ was always small, with a maximum of 3.7$\times
10^{-3}$.

Setting $w_{t}=0$, the incident pressure for a broadband, tapered pulse is
given by

$\displaystyle\begin{split}p_{i}(\mathbf{r},t)=p_{0}&\psi(\mathbf{r})\\\
\times&\exp\left(i\omega_{0}(t-t_{i}(\mathbf{r})\right)\\\
\times&\exp\left(-(t-t_{i}(\mathbf{r}))^{2}/\tau^{2}\right)\,,\end{split}$
(52)

where

$\displaystyle
t_{i}(\mathbf{r})=\frac{\mathbf{r}\cdot\hat{\mathbf{k}}_{i}}{c}=\frac{-x\cos\theta_{i}-z\sin\theta_{i}}{c},$
(53)

and

$\displaystyle\psi(\mathbf{r})=\exp\left(-(x-z\cot\theta_{i})^{2}/g^{2}\right)$
(54)

is the spatial taper. The vertical component of the acoustic particle velocity
can be computed as (Pierce, 1994, p. 19)

$\displaystyle v_{z}(\mathbf{r},t)$ $\displaystyle=\frac{\partial}{\partial
z}\Phi(\mathbf{r},t),$ (55)

and the acoustic pressure as

$\displaystyle p(\mathbf{r},t)$
$\displaystyle=-\rho_{0}\frac{\partial}{\partial t}\Phi(\mathbf{r},t),$ (56)

where $\Phi$ is the acoustic velocity potential. Using Eq. (52) in (56), the
velocity potential can be computed using Eqs. 3.322(1-2) from Gradshteyn and
Rhyzik (2007),

$\displaystyle\begin{split}\Phi(\textbf{r},t)&=\frac{-p_{0}}{\rho_{0}}\psi(\mathbf{r})M(t-t_{i}(\mathbf{r}),f_{0},\tau)\end{split}$
(57)

where

$\displaystyle
M(t,f_{0},\tau)=\frac{\sqrt{\pi}\tau}{2}e^{-\pi^{2}f_{0}^{2}\tau^{2}}\operatorname{erf}\left[t/\tau-i\pi
f_{0}\tau\right],$ (58)

and $\operatorname{erf}\left[z\right]$ is the error function (Abramowitz and
Stegun, 1972, p. 297). If the opposite time convention is used, then the minus
sign inside the error function argument must be changed to a plus. In the
integration used to compute Eq. (57), an integration constant independent of
time is required, but is assumed to be zero here, since derivatives are taken
to compute $p$ and $v_{z}$.

Inserting Eq. (57) into Eq. (55), one obtains

$\displaystyle\begin{split}v_{iz}(\mathrm{r},t&)=-\frac{p_{0}}{\rho_{0}}\left[\frac{2\cot\theta_{i}}{g^{2}}(x-z\cot\theta_{i})\psi(\mathbf{r})M(t-t_{i}(\mathbf{r}),f_{0},\tau)\right.\\\
&\left.+\frac{\sin\theta_{i}\psi(\mathbf{r})}{c}\exp\left(i\omega_{0}(t-t_{i}(\mathbf{r}))-(t-t_{i}(\mathbf{r}))^{2}/\tau^{2}\right)\right]\end{split}$
(59)

To compute the effective energy flux passing through the $z=0$ plane, Eq. (28)
is computed using Eqs. (52) and (59). The time variable is set to
$t=R/c+t^{\prime}$, to match the convention in Eq. (23). Making these
substitutions,

$\displaystyle
E_{f}^{\prime}=\frac{|p_{0}|^{2}}{2\rho_{0}}\left(\frac{2\cot\theta_{i}}{g^{2}}I_{1}+\frac{\sin\theta_{i}}{c}I_{2}\right)\,.$
(60)

The two integrals are defined as

$\displaystyle\begin{split}I_{1}&=\frac{\sqrt{\pi}\tau}{2}e^{-(\omega_{0}\tau/2)^{2}}\\\
&\times
Re\biggl{[}\int\limits_{-\infty}^{\infty}x\,\mathrm{erf}\left(\frac{t^{\prime}}{\tau}+\frac{x\nu_{x}}{c\tau}+i\frac{\omega_{0}\tau}{2}\right)\biggr{.}\\\
&\times\biggl{.}e^{-2x^{2}/g^{2}-\left(\frac{t^{\prime}}{\tau}+\frac{x\nu_{x}}{c\tau}\right)^{2}+i\omega_{0}\tau\left(\frac{t^{\prime}}{\tau}+\frac{x\nu_{x}}{c\tau}\right)}\,\mathrm{d}x\biggr{]}\end{split}$
(61) $\displaystyle I_{2}$
$\displaystyle=\int\limits_{-\infty}^{\infty}e^{-2x^{2}/g^{2}-2(t^{\prime}+x\nu_{x}/c)^{2}/\tau^{2}}\,\mathrm{d}x,$
(62)

where the identity $\mathrm{erf}^{\ast}(z)=\mathrm{erf}(z^{\ast})$ was used.

The integral $I_{1}$ may be cast into a generic form by completing the square
in the exponent, resulting in

$\displaystyle
I_{1}=\frac{\sqrt{\pi}\tau}{2}Re\left[e^{-\zeta}\int\limits_{-\infty}^{\infty}x\,\mathrm{erf}(ax+b)e^{-(\alpha
x+\beta)^{2}}\,\mathrm{d}x\right]\,,$ (63)

with

$\displaystyle a$ $\displaystyle=\nu_{x}/(c\tau)$ (64) $\displaystyle b$
$\displaystyle=t^{\prime}/\tau+i\frac{\omega_{0}\tau}{2}$ (65)
$\displaystyle\alpha$
$\displaystyle=\sqrt{\frac{2}{g^{2}}+\left(\frac{\nu_{x}}{c\tau}\right)^{2}}$
(66) $\displaystyle\beta$ $\displaystyle=\frac{ab^{\ast}}{\alpha}$ (67)
$\displaystyle\zeta$
$\displaystyle=(b^{\ast})^{2}+2\left(\frac{\omega_{0}\tau}{2}\right)^{2}-\left(\frac{ab^{\ast}}{\alpha}\right)^{2}\,.$
(68)

This integral can be transformed via integration by parts, $\int
u\,\mathrm{d}v=uv|_{-\infty}^{\infty}-\int v\,\mathrm{d}u$, using

$\displaystyle u$ $\displaystyle=\mathrm{erf}(ax+b)$ (69) $\displaystyle dv$
$\displaystyle=xe^{-(\alpha x+\beta)^{2}}dx$ (70) $\displaystyle du$
$\displaystyle=\frac{2a}{\sqrt{\pi}}e^{-(ax+b)^{2}}dx$ (71) $\displaystyle v$
$\displaystyle=(2\alpha^{2})^{-1}\left(-e^{-(\alpha
x+\beta)^{2}}-\beta\sqrt{\pi}\mathrm{erf}(\alpha x+\beta)\right).$ (72)

Using the large-argument expansion of the error function, the product of the
error function terms in $uv$ evaluated at $+\infty$ is 1, and is 1 at
$x=-\infty$, so the difference is zero. The exponential in $v$ evaluated at
$\pm\infty$ is zero. Thus the $uv$ term evaluates to zero.

Turning to $-\int v\,\mathrm{d}u$,

$\displaystyle\begin{split}-\int
v\,\mathrm{d}u=&\frac{a}{\sqrt{\pi}\alpha^{2}}\int\limits_{-\infty}^{\infty}e^{-(ax+b)^{2}}\\\
\times&\left(e^{-(\alpha x+\beta)^{2}}+\beta\sqrt{\pi}\mathrm{erf}(\alpha
x+\beta)\right)\,\mathrm{d}x\,.\end{split}$ (73)

There are two integrals in the previous equation:

$\displaystyle I_{11}$
$\displaystyle=\int\limits_{-\infty}^{\infty}e^{-(ax+b)^{2}-(\alpha
x+\beta)^{2}}\,\mathrm{d}x$ (74) $\displaystyle I_{12}$
$\displaystyle=\beta\sqrt{\pi}\int\limits_{-\infty}^{\infty}\mathrm{erf}(\alpha
x+\beta)e^{-(ax+b)^{2}}\,\mathrm{d}x\,.$ (75)

$I_{11}$ straightforward to compute, and is

$\displaystyle
I_{11}=\frac{\sqrt{\pi}}{\sqrt{a^{2}+\alpha^{2}}}e^{-(b\alpha-a\beta)^{2}/(a^{2}+\alpha^{2})}\,,$
(76)

with the restriction that both $a$ and $\alpha$ are greater than zero (which
is satisfied here).

$I_{12}$ can be solved using the following tabulated integral, which is Eq.
(13) of section 4.3 in Ng and Geller (1969),

$\displaystyle\int\limits_{-\infty}^{\infty}\mathrm{erf}(y)e^{-(zy+w)^{2}}\,\mathrm{d}y=-\frac{\sqrt{\pi}}{z}\mathrm{erf}\left(\frac{w}{\sqrt{z^{2}+1}}\right).$
(77)

Using Eq. (77), in combination with the substitutions,

$\displaystyle y$ $\displaystyle=\alpha x+\beta$ (78) $\displaystyle z$
$\displaystyle=\frac{a}{\alpha}$ (79) $\displaystyle w$
$\displaystyle=b-\frac{\beta a}{\alpha}$ (80)

the integral can be calculated as

$\displaystyle I_{12}$
$\displaystyle=-\frac{\beta\pi}{a}\mathrm{erf}\left(\frac{b\alpha-a\beta}{\sqrt{\alpha^{2}+a^{2}}}\right)\,.$
(81)

Here, $Re[z^{2}]=Re[(a/\alpha)^{2}]$ must be greater than zero, which is
satisfied here. Note that $y$ is complex, so the integration limits in Eq.
(77) should have the imaginary part of $\beta$ added to both. In effect, this
changes the integral from over the real line, to a line with a constant
imaginary part, parallel to the real line. This change has no effect on the
integral for the parameters in this work, as it has been verified using
numerical quadrature that Eq. (81) produces the correct result.

Putting these results together, the integral $I_{1}$ can be calculated as

$\displaystyle I_{1}=\tau\frac{a}{2\alpha^{2}}\begin{aligned}
Re&\left\\{e^{-\zeta}\left[\frac{\sqrt{\pi}}{\sqrt{a^{2}+\alpha^{2}}}e^{-(b\alpha-a\beta)^{2}/(a^{2}+\alpha^{2})}\right.\right.\\\
-&\left.\left.\frac{\pi\beta}{a}\mathrm{erf}\left(\frac{b\alpha-a\beta}{\sqrt{\alpha^{2}+a^{2}}}\right)\right]\right\\}\,.\end{aligned}$
(82)

Turning to $I_{2}$, it can be computed as

$\displaystyle
I_{2}=\sqrt{\frac{\pi}{2}}\frac{e^{-\frac{2(ct^{\prime}/\nu_{x})^{2}}{g^{2}+(c\tau/\nu_{x})^{2}}}}{\sqrt{g^{-2}+(c\tau/\nu_{x})^{-2}}}\,.$
(83)

Using (82) and (83) in (60), the energy flux can be written as

$\displaystyle\begin{split}E_{f}^{\prime}&=\frac{|p_{0}|^{2}}{2\rho_{0}}\left(\frac{\sin\theta_{i}}{c}\sqrt{\frac{\pi}{2}}\frac{e^{-\frac{2(ct^{\prime}/\nu_{x})^{2}}{g^{2}+(c\tau/\nu_{x})^{2}}}}{\sqrt{g^{-2}+(c\tau/\nu_{x})^{-2}}}\right.+\frac{2\cot\theta_{i}}{g^{2}}\frac{\tau\nu_{x}}{2c\tau}\frac{1}{2/g^{2}+(\nu_{x}/(c\tau))^{2}}\mathrm{Re}\left\\{\vphantom{\frac{\sqrt{\pi}}{\sqrt{2/g^{2}+(\nu_{x}/(c\tau))^{2}}}}\right.\\\
\times&\exp\left[\left(\left(\frac{\nu_{x}}{c\tau\sqrt{2/g^{2}+(\nu_{x}/(c\tau))^{2}}}\right)^{2}-1\right)\left(\left(\frac{t^{\prime}}{\tau}\right)^{2}-i\omega_{0}\tau-\left(\frac{\omega_{0}\tau}{2}\right)^{2}\right)-2\left(\frac{\omega_{0}\tau}{2}\right)^{2}\right]\\\
\times&\left(\frac{\sqrt{\pi}}{\sqrt{2}\sqrt{1/g^{2}+(\nu_{x}/(c\tau))^{2}}}\exp\left[-\left(\frac{b\alpha-a\beta}{\sqrt{a^{2}+\alpha^{2}}}\right)^{2}\right]\right.-\left.\left.\left.\frac{\pi(t^{\prime}/\tau-i\omega_{0}\tau/2)}{\sqrt{2/g^{2}+(\nu_{x}/(c\tau))^{2}}}\operatorname{erf}\left[\frac{b\alpha-a\beta}{\sqrt{a^{2}+\alpha^{2}}}\right]\vphantom{\frac{\sqrt{\pi}}{\sqrt{2/g^{2}+(\nu_{x}/(c\tau))^{2}}}}\right)\right\\}\right)\,,\end{split}$
(84)

where the factor

$\displaystyle\frac{b\alpha-a\beta}{\sqrt{a^{2}+\alpha^{2}}}$
$\displaystyle=\frac{\left(t^{\prime}/\tau+i\omega_{0}\tau/2\right)\sqrt{2/g^{2}+(\nu_{x}/(c\tau))^{2}}-(\nu_{x}/(c\tau))^{2}\frac{t^{\prime}/\tau-i\omega_{0}\tau/2}{\sqrt{2/g^{2}+(\nu_{x}/(c\tau))^{2}}}}{\sqrt{2}\sqrt{1/g^{2}+(\nu_{x}/(c\tau))^{2}}}$
(85)

is not explicitly included in Eq. (84) for space reasons.

It is advantageous to find a simplification for the error function with a
complex argument appearing in (84). Many software packages only implement the
error function with a real argument (e.g. the basic Matlab installation, the
GNU C math library, or the GNU C++ cmath library), and an approximation of
this special function in terms of elementary functions makes implementing the
broadband scattering cross section much simpler numerically. The approximation
derived below is very accurate for all parameters investigated here, and has a
more succinct representation than Eq. (84).

To see what kind of approximation is appropriate, it is necessary to know if
the error function argument is large, or small. To obtain an estimate of the
order of magnitude, for this problem at moderate and low grazing angles, $g\gg
c\tau/\nu_{x}$. Therefore, $\alpha\approx a$, and $\alpha^{2}+a^{2}\approx
2a^{2}$. The argument of the error function is $q$, and combined with the
previous approximation, can be simplified to

$\displaystyle
q\approx\frac{b-b^{\ast}}{\sqrt{2}}=\frac{i\omega_{0}\tau}{\sqrt{2}}$ (86)

The largest value of $f_{0}\tau$ examined in this work is 1.7, resulting in
$q\approx 7.5i$. Therefore the magnitude of $q$ is large, and its real part is
small compared to its imaginary part. The large-argument asymptotic series
approximation is

$\displaystyle\mathrm{erf}(q)=1-\frac{e^{-q^{2}}}{q\sqrt{\pi}}\left(1-\frac{1}{2q^{2}}+\frac{3}{4q^{4}}+\dots\right).$
(87)

For $f_{0}\tau=1.7$, $\exp(-q^{2})\approx 10^{24}$, so additive factor of 1
outside the parentheses can be safely neglected. The terms inside the
parentheses become progressively smaller. The second term is about $10^{-2}$,
and the third term is about $10^{-4}$, for $f_{0}\tau=1.7$, so only retaining
the first term is appropriate. If longer pulse lengths are used, this
approximation improves, and if shorter pulse lengths are used, then the
approximation is poorer. If extremely short pulses are used, then more terms
could be kept in the approximation.

Keeping only the largest terms in the power series of $\mathrm{erf}(q)$, the
integral $I_{1}$ is approximately

$\displaystyle I_{1}$
$\displaystyle\approx\sqrt{\frac{\pi}{2}}\frac{e^{-\frac{2(ct^{\prime}/\nu_{x})^{2}}{g^{2}+(c\tau/\nu_{x})^{2}}}}{\sqrt{a^{2}+\alpha^{2}}}\frac{\tau
a\sqrt{2}}{2\alpha^{2}}\left(1+Re\left\\{\frac{\beta}{a}\frac{\alpha^{2}+a^{2}}{b\alpha-a\beta}\right\\}\right).$
(88)

where $e^{-\zeta}$ has been combined with the other exponential terms.

This approximation for $I_{1}$ can be used to obtain a good approximation to
the effective energy flux

$\displaystyle\begin{split}E^{\prime}_{f}\approx&\frac{|p_{0}|^{2}}{2\rho_{0}c}L_{eff}e^{-\frac{2(ct^{\prime}/\nu_{x})^{2}}{g^{2}+(c\tau/\nu_{x})^{2}}}\sin\theta_{i}\\\
&\times\left[1+\frac{\cot\theta_{i}}{\sin\theta_{i}}\frac{\nu_{x}}{g^{2}\alpha^{2}}\left(1+Re\left\\{\frac{\beta(\alpha^{2}+a^{2})}{a(b\alpha-a\beta)}\right\\}\right)\right]\end{split}$
(89)

where $L_{eff}$ is defined in the main text, in Eq. (34). Substituting in the
definitions of $a$, $b$, $\alpha$ and $\beta$, and taking the real part
results in Eqs. (33) though (39) in the main text.

## References

  * Abraham and Lyons (2002) Abraham, D. A., and Lyons, A. P. (2002). “Novel physical interpretations of K-distributed reverberation,” IEEE J. Ocean Eng. 27, 800–813.
  * Abraham and Lyons (2004) Abraham, D. A., and Lyons, A. P. (2004). “Reverberation envelope statistics and their dependence on sonar bandwidth and scattering patch size,” IEEE J. Ocean Eng. 29, 126–137.
  * Abramowitz and Stegun (1972) Abramowitz, M., and Stegun, I. A. (1972). _Handbook of Mathematical Functions_ (Dover, Mineola, NY).
  * Anderson _et al._ (1999) Anderson, E., Bai, Z., Bischof, C., Blackford, S., Demmel, J., Dongarra, J., Croz, J. D., Greenbaum, A., Hammarling, S., McKenney, A., and Sorensen, D. (1999). _LAPACK Users’ Guide_ , 3rd ed. (SIAM, Philadelphia,PA).
  * Atkinson (1997) Atkinson, K. (1997). _The Numerical Solution of Integral Equations of the Second Kind_ (Cambridge University Press, Cambridge, UK).
  * Bachmann (1973) Bachmann, W. (1973). “A theoretical model for the backscattering strength of a composite‐roughness sea surface,” J. Acoust. Soc. Am. 54(3), 712–716, http://dx.doi.org/10.1121/1.1913652.
  * Bellettini and Pinto (2009) Bellettini, A., and Pinto, M. (2009). “Design and experimental results of a 300-khz synthetic aperture sonar optimized for shallow-water operations,” IEEE J. Ocean Eng. 34(3), 285–293.
  * Blacknell _et al._ (1992) Blacknell, D., Blake, A. P., Oliver, C. J., and White, R. G. (1992). “A comparison of SAR multilook registration and contrast optimisation autofocus algorithms applied to real SAR data,” in _92 International Conference on Radar_ , pp. 363–366.
  * Brown _et al._ (2017) Brown, D. C., Johnson, S. F., and Olson, D. R. (2017). “A point-based scattering model for the incoherent component of the scattered field,” J. Acoust. Soc. Am. 141(3), EL210–EL215, 10.1121/1.4976584.
  * Callow (2003) Callow, H. J. (2003). “Signal processing for synthetic aperture sonar image enhancement,” Ph.D. thesis, University of Canterbury, Christchurch, New Zealand.
  * Dillon (2018) Dillon, J. (2018). “Real-time interferometric SAS processing with ultra-low power consumption,” in _OCEANS 2018 MTS/IEEE Charleston_ , pp. 1–6, 10.1109/OCEANS.2018.8604512.
  * Fossum _et al._ (2008) Fossum, T. G., Sæbø, T. O., Langli, B., Callow, H., and Hansen, R. E. (2008). “HISAS 1030 - High resolution interferometric synthetic aperture sonar,” in _Proc. Can. Hydrog. Conf._ , Victoria BC, Canada.
  * Frigo and Johnson (2005) Frigo, M., and Johnson, S. G. (2005). “The design and implementation of FFTW3,” Proceedings of the IEEE 93(2), 216–231 special issue on “Program Generation, Optimization, and Platform Adaptation”.
  * Fritsch and Carlson (1980) Fritsch, F. N., and Carlson, R. E. (1980). “Monotone piecewise cubic interpolation,” SIAM Journal on Numerical Analysis 17(2), 238–246, 10.1137/0717021.
  * Galusha _et al._ (2018) Galusha, A. P., Keller, J. M., Zare, A., and Galusha, G. (2018). “A fast target detection algorithm for underwater synthetic aperture sonar imagery,” in _Detection and Sensing of Mines, Explosive Objects, and Obscured Targets XXIII_ , edited by J. C. Isaacs and S. S. Bishop, SPIE, 10.1117/12.2304976.
  * Gradshteyn and Rhyzik (2007) Gradshteyn, I., and Rhyzik, I. (2007). _Table of Integrals, Series, and Products_ , Seventh ed. (Elsevier Academic Press, Burlington, MA).
  * Hackbusch (2015) Hackbusch, W. (2015). _Hierarchical Matrices: Algorithms and Analysis_ , 1st ed. (Springer Publishing Company, Incorporated).
  * Harris (1978) Harris, F. J. (1978). “On the use of windows for harmonic analysis with the discrete Fourier transform,” Proceedings of the IEEE 66(1), 51–83, 10.1109/proc.1978.10837.
  * Hefner (2015) Hefner, B. T. (2015). “Inversion of high frequency acoustic data for sediment properties needed for the detection and classification of UXOs,” Final Report, https://www.serdp-estcp.org/content/download/34593/333838/file/MR-2229-FR.pdf.
  * Hellequin _et al._ (2003) Hellequin, L., Boucher, J. M., and Lurton, X. (2003). “Processing of high-frequency multibeam echo sounder data for seafloor characterization,” IEEE J. Ocean Eng. 28(1), 78–89, 10.1109/JOE.2002.808205.
  * Ishimaru (1978) Ishimaru, A. (1978). _Wave Propagation and Scattering in Random Media_ , 2 (Academic Press).
  * Jackson and Olson (2020) Jackson, D., and Olson, D. R. (2020). “The small-slope approximation for layered, fluid seafloors,” The Journal of the Acoustical Society of America 147(1), 56–73, 10.1121/10.0000470, 10.1121/10.0000470.
  * Jackson _et al._ (1986a) Jackson, D. R., Baird, A. M., Crisp, J. J., and Thomson, P. A. G. (1986a). “High-frequency bottom backscatter measurements in shallow water,” J. Acoust. Soc. Am. 80(4), 1188–1199, http://dx.doi.org/10.1121/1.393809.
  * Jackson and Ivakin (1998) Jackson, D. R., and Ivakin, A. N. (1998). “Scattering from elastic sea beds: First-order theory,” J. Acoust. Soc. Am. 103, 336 – 345.
  * Jackson and Richardson (2007) Jackson, D. R., and Richardson, M. D. (2007). _High-Frequency Seafloor Acoustics_ , 1st ed. (Springer, New York, NY).
  * Jackson _et al._ (1986b) Jackson, D. R., Winebrenner, D. P., and Ishimaru, A. (1986b). “Application of the composite roughness model to high-frequency bottom backscattering,” J. Acoust. Soc. Am. 79, 1410–1422.
  * Jakeman (1980) Jakeman, E. (1980). “On the statistics of K-distributed noise,” J. Phys. A: Math. Gen. 13, 31–48.
  * Kur’yanov (1963) Kur’yanov, B. (1963). “The scattering of sound at a rough surface with two types of irregularity,” Sov. Phys. Acoust. 8(3), 252–257.
  * Kwon and Nasrabadi (2005) Kwon, H., and Nasrabadi, N. (2005). “Kernel RX-algorithm: a nonlinear anomaly detector for hyperspectral imagery,” IEEE Transactions on Geoscience and Remote Sensing 43(2), 388–397, 10.1109/tgrs.2004.841487.
  * Landau and Lifshitz (1987) Landau, L. D., and Lifshitz, E. M. (1987). _Fluid Mechanics_ , 2 ed. (Butterworth Heinemann, New York).
  * Li and Johnson (2017) Li, H., and Johnson, J. T. (2017). “On the amplitude distributions of bistatic scattered fields from rough surfaces,” IEEE Transactions on Geoscience and Remote Sensing 55(12), 6883–6892, 10.1109/tgrs.2017.2735862.
  * Liu (2009) Liu, Y., ed. (2009). _Fast Multipole Boundary Element Method: Thory and Applications in Engineering_ (Cambridge University Press, Cambridge, UK).
  * Lupien (1999) Lupien, V. (1999). “The role of scale structure in scattering from random rough surfaces,” J. Acoust. Soc. Am. 105, 2187–2202.
  * Lyons _et al._ (2016) Lyons, A., Olson, D. R., and Hansen, R. (2016). “Quantifying the effect of random seafloor roughness on high-frequency synthetic aperture sonar image statistics,” in _Acoustic and Environmental Variability, Fluctuations, and Coherence_ , Institute of Acoustics, Cambridge, UK, Vol. 38, pp. 151–158.
  * Lyons _et al._ (2010) Lyons, A. P., Abraham, D. A., and Johnson, S. F. (2010). “Modeling the effect of seafloor ripples on synthetic aperture sonar speckle statistics,” IEEE J. Ocean Eng. 35(2), 242–249, 10.1109/JOE.2009.2039656.
  * Lyons _et al._ (2009) Lyons, A. P., Johnson, S. F., Abraham, D. A., and Pouliquen, E. (2009). “High-frequency scattered envelope statistics of patchy seafloors,” IEEE J. Ocean Eng. 34, 451–458.
  * Lysanov (1973) Lysanov, Y. P. (1973). “A property of the scattering coeffciient in the Fraunhofer zone,” Sov. Phys. Acoust. 18, 505 – 506.
  * Marston and Plotnick (2015) Marston, T. M., and Plotnick, D. S. (2015). “Semiparametric statistical stripmap synthetic aperture autofocusing,” IEEE Transactions on Geoscience and Remote Sensing 53(4), 2086–2095, 10.1109/TGRS.2014.2353515.
  * McDaniel and Gorman (1983) McDaniel, S. T., and Gorman, A. D. (1983). “An examination of the composite-roughness scattering model,” J. Acoust. Soc. Am. 73, 1476–1486.
  * Moccia and Renga (2011) Moccia, A., and Renga, A. (2011). “Spatial resolution of bistatic synthetic aperture radar: Impact of acquisition geometry on imaging performance,” IEEE Trans. Geosci. Remote Sens. 49(10), 3487–3503, 10.1109/TGRS.2011.2115250.
  * Ng and Geller (1969) Ng, E. W., and Geller, M. (1969). “A table of integrals of the error functions,” J. Research of the National Bureau of Standards – B. Mathematical Sciences 73B(1), 1–20.
  * Olson _et al._ (2019) Olson, D. R., Lyons, A. P., Abraham, D. A., and Sæbø, T. O. (2019). “Scattering statistics of rock outcrops: Model-data comparisons and Bayesian inference using mixture distributions,” J. Acoust. Soc. Am. 145(2), 761–774, 10.1121/1.5089892.
  * Olson _et al._ (2016) Olson, D. R., Lyons, A. P., and Sæbø, T. O. (2016). “Measurements of high-frequency acoustic scattering from glacially eroded rock outcrops,” J. Acoust. Soc. Am. 139(4), 1833–1847, 10.1121/1.4945589.
  * Pierce (1994) Pierce, A. D. (1994). _Acoustics: An introduction to its physical principles and applications_ (the Acoustical Society of America).
  * Pinto (2011) Pinto, M. (2011). “Interferometric synthetic aperture sonar design optimized for high area coverage shallow water bathymetric survey,” in _Proc. 4th Int. Conf Exhib. Underwater Ac. Meas._ , Kos, Greece, pp. 505–512.
  * Sauter and Schwab (2011) Sauter, S. A., and Schwab, C. (2011). _Boundary Element Methods_ (Springer-Verlag, Heidelberg).
  * Sternlicht _et al._ (2016) Sternlicht, D. D., Fernandez, J. E., Prater, J. L., Weaver, J. N., Isaacs, J. C., Montgomery, T. C., Loeffler, C. M., and Purcell, M. (2016). “Advanced sonar technologies for high clearance rate mine countermeasures,” in _OCEANS 2016 MTS/IEEE Monterey_ , pp. 1–8, 10.1109/OCEANS.2016.7761133.
  * Tatarski (1961) Tatarski, V. I. (1961). _Wave Propagation in a Turbulent Medim_ (Dover Reprint, New York).
  * Thorsos (1988) Thorsos, E. I. (1988). “The validity of the Kirchhoff approximation for rough surface scattering using a Gaussian roughness spectrum,” J. Acoust. Soc. Am. 83, 78–92.
  * Thorsos and Jackson (1989) Thorsos, E. I., and Jackson, D. R. (1989). “The validity of the perturbation approximation for rough surface scattering using a Gaussian roughness spectrum,” J. Acoust. Soc. Am. 86, 261–277.
  * Thorsos and Jackson (1991) Thorsos, E. I., and Jackson, D. R. (1991). “Studies of scattering theory using numerical methods,” Waves Random Media 1(3), S165–S190, 10.1088/0959-7174/1/3/014.
  * Toporkov _et al._ (1998) Toporkov, J., Marchand, R., and Brown, G. (1998). “On the discretization of the integral equation describing scattering by rough conducting surfaces,” Antennas and Propagation, IEEE Transactions on 46(1), 150 –161, 10.1109/8.655462.
  * Urick (1954) Urick, R. J. (1954). “The backscattering of sound from a harbor bottom,” J. Acoust. Soc. Am. 26, 231–235.
  * Urick (1983) Urick, R. J. (1983). _Principles of Underwater Sound_ , 3rd ed. (Peninsula, Los Altol Hills, California).
  * Valenzuela and Liang (1971) Valenzuela, and Liang (1971). “On the statistics of sea clutter,” US Naval Research Laboratory Report 7349\.
  * Wang and Bovik (2002) Wang, Z., and Bovik, A. C. (2002). “A universal image quality index,” IEEE Signal Processing Letters 9(3), 81–84, 10.1109/97.995823.
  * Williams (2015) Williams, D. P. (2015). “Fast unsupervised seafloor characterization in sonar imagery using lacunarity,” IEEE Transactions on Geoscience and Remote Sensing 53(11), 6022–6034, 10.1109/TGRS.2015.2431322.
  * Williams (2018) Williams, D. P. (2018). “The Mondrian detection algorithm for sonar imagery,” IEEE Transactions on Geoscience and Remote Sensing 56(2), 1091–1102, 10.1109/TGRS.2017.2758808.
  * Williams _et al._ (2002) Williams, K., Jackson, D., Thorsos, E., Tang, D., and Briggs, K. (2002). “Acoustic backscattering experiments in a well characterized sand sediment: Data/model comparisons using sediment fluid and Biot models,” IEEE J. Ocean Eng. 27, 376–387.
  * Winebrenner and Ishimaru (1986) Winebrenner, D., and Ishimaru, A. (1986). “On the far-field approximation for scattering from randomly rough surfaces,” Antennas and Propagation, IEEE Transactions on 34(6), 847–849, 10.1109/TAP.1986.1143901.
  * Wu (2000) Wu, T. W., ed. (2000). _Boundary Element Acoustics: Fundamentals and Computer Codes_ (WIT Press, Southampton, UK).
  * Yang and McDaniel (1991) Yang, C. C., and McDaniel, S. T. (1991). “Fourth moments of acoustic waves forward scattered by a rough ocean surface,” Waves in Random Media 1(4), 419–439, 10.1088/0959-7174/1/4/012.