# TetraScatt model: Born approximation for the estimation of acoustic
dispersion of fluid-like objects of arbitrary geometries.

Edmundo F. Lavia(1,2) Guadalupe Cascallares(3,4) Juan D. Gonzalez(1,2)
1Acoustic Propagation Department. Argentinian Navy Research Office (DIIV)
Laprida 555 Vicente Lopez, (1638) Buenos Aires, Argentina
2UNIDEF (National Council of Scientific and Technical Research – Ministry of
Defense)
3National Institute for Fisheries Research and Development (INIDEP)
Paseo Victoria Ocampo 1, Escollera Norte, Mar del Plata, Argentina.
4CONICET.

((10 May 2022))

###### Abstract

Modelling the acoustic scattering response due to penetrable objects of
arbitrary shapes, such as those of many marine organisms, can be
computationally intensive, often requiring high-performance computing
equipment when considering a completely general situation. However, when the
physical properties (sound speed and density) of the scatterer object under
consideration are similar to those of the surrounding medium, the Born
approximation provides a computationally efficient way to calculate the
scattering. For simple geometrical shapes like spheres and spheroids, the
acoustic scattering in the far field can be evaluated through the Born
approximation recipe as a formula which has been historically employed to
predict the response of weakly scattering organisms, such as zooplankton.
Moreover, the Born approximation has been extended to bodies whose geometry
can be described as a collection of non-circular rings centred on a smooth
curve. In this work, we have developed a numerical approach to calculate the
far-field backscattering by arbitrary 3D objects under the Born approximation.
The object’s geometry is represented by a volumetric mesh composed of
tetrahedrons, and the computation is efficiently performed through analytical
3D integration, yielding a solution expressed in terms of elementary functions
for each tetrahedron. The method’s correctness has been successfully validated
against benchmark solutions. Additionally, we present acoustic scattering
results for species with complex shapes. To enable other researchers to use
and validate the method a computational package named tetrascatt implemented
in the R programming language was developed and published in the CRAN
(Comprehensive R Archive Network).

###### keywords:

Born approximation; weakly scattering; arbitrary geometry.

##  1 Introduction

Acoustic sampling has proven to be a useful tool to estimate biomass in
fishery acoustics. The physical and mathematical theory on which it is based
was established by the pioneering works of some of most important researchers
from XIX century and consequently methods and techniques aimed to solve the
acoustic scattering problem when certain simplifying assumptions are fulfilled
are bonded to such big names as Helmholtz, Kirchhoff or Born.

When the shape of the body which produces the scattering can be approximated
by a simple geometrical shape (sphere, spheroid or cylinder, for instance) the
previously mentioned approaches even provide an analytical solution. However,
the bodies of interest in fisheries acoustics are usually not simple
geometries, so it is necessary fall to numerical techniques.

Accurate computation of scattering from an arbitrary shaped obstacle can be a
challenging task, as it requires solving the exterior Helmholtz equation in
three dimensions under appropriate boundary conditions. When the scattering
body can be considered as an impenetrable one, a reasonable approximation may
be the Kirchhoff approximation [1]. It gives valid solutions for (mainly)
convex bodies in the high frequency regime, and its application avoids solving
many of the problems that occur when calculating the exact solution, such as
the resolution of high-dimensional matrix systems and the computation of
singular integrals.

When the scattering objects are penetrable, but their physical properties do
not vary significantly with respect to the properties of the medium where the
wave propagates, an approach known as the Born Approximation (BA) can be made
[2]. The conditions of this approximation are what is known in the literature
as weakly scattering. The key assumption is that, up to first order, the
perturbation induced for the impinging wave in the interior of this weakly
scattering object is neglected. As a result, the acoustic field within the
object is assumed to be identical to the incident field. This allows us to
solve the scattering through an explicit evaluation of an integral, without
the need to solve an integral equation or its discretized version.

Improved versions of the original BA are the Distorted Wave Born Approximation
(DWBA), and also its phase compensated version [3, 4] which are of widespread
use in many works regarding scattering by weakly organisms as zooplankton. The
first models focused on marine organisms were based on canonical shapes [5].
Then different schemes allows to deal with complex bodies in an approximated
way, see [6, 4] and references therein. Further, use of the Born approximation
was extended to bodies whose geometry can be described as a collection of non-
circular rings centred on a smooth curve [7]. Such approaches are not
applicable to incidences far from the broadside incidence.

A work addressing the BA for arbitrary shaped geometries is presented in the
Ref. 8 where the pressure field is calculated resorting to a fundamental
theorem on diffraction tomography which involves Fourier transforms.

In this work, a numerical approach to calculate the far field backscattering
by arbitrary three-dimensional objects under the Born approximation is
developed and computationally implemented in the R programming language. The
model is based on direct (analytic) integration over a discretized version of
the scattering volume, which only requires elementary techniques avoiding thus
the use of quadrature rules. In this way a fast and general approach is
obtained. A computational code is made publicly available as an R package
named tetrascatt[9] in the CRAN (Comprehensive R Archive Network) repository.
The TetraScatt model is validated using as benchmarks exact solutions of the
Born approximation and exact solutions of general scattering problems
evaluated in weakly conditions. Then the code is used to explore the validity
range of the Born approximation itself.

This work is organized as follows, in the Section 2 a brief summary of the BA
is provided followed by a description of the mathematical and numerical
implementation of the model. Section 3 shows a set of tests against benchmark
solutions which gives a validation for the model. In Section 4 numerical
experiments are conducted to explore the range of validity where the BA gives
predictions with an acceptable error margin. Exact solutions of general
scattering are employed for this task. As an application of the model, a
backscattering calculation from a native Argentinian species of crustacean is
presented in Section 5. Finally, conclusions of the work are drawn in Section
6.

##  2 Model description

###  2.1 Scattering by the Born approximation

The acoustic scattering problem we considered here is depicted in Fig. 1. A
plane harmonic wave of frequency $f$ propagating in a medium of sound speed
and density $c_{0},\rho_{0}$ is incident on a scattering object with internal
properties $c_{1},\rho_{1}$. Two wavenumbers $k_{i}=2\pi f/c_{i}$ (with
$i=0,1$) are thus defined.

Figure 1: Geometry of the scattering by an arbitrary body of volume $V$.

The backscattering far-field form function $f_{\infty}$ under the Born
approximation[2] can be expressed as the volume integral

$f_{\infty}^{\text{BA}}=\frac{k_{0}^{2}}{4\pi}\>\int_{V}\>(\gamma_{\kappa}-\gamma_{\rho})\>e^{2ik_{0}\hat{\bf{k}}\cdot\bf{x}}\>dV(\bf{x}),$
(1)

where the factors $\gamma$ are related to the density $\rho$ and
compressibility $\kappa=1/(\rho c)^{2}$ of the two media and can be expressed
in terms of the ratios $g\equiv\rho_{1}/\rho_{0}$ and $h\equiv c_{1}/c_{0}$,
according to

$\gamma_{\rho}=\frac{\rho_{1}-\rho_{0}}{\rho_{1}}=\frac{g-1}{g}\qquad\qquad\gamma_{\kappa}=\frac{\kappa_{1}-\kappa_{0}}{\kappa_{0}}=\frac{1-gh^{2}}{gh^{2}}.$

Since the Born approximation assumes similar properties for the exterior media
“0” and interior media “1”, it follows that $g\approx 1$ and $h\approx 1$. A
modified version of this approach is the distorted wave Born approximation
(DWBA)[3] which replaces the wave vector $\mathbf{k}_{0}$ in both the phase of
the internal field and the Green function111After performing the far-field
limit. by $\mathbf{k}_{1}$, thus representing a distorted wave.

This modification provides, in the backscattering calculations, a slightly
improvement over the original expression. For homogeneous media the DWBA
results in

$f_{\infty}^{\text{DWBA}}=\frac{k_{1}^{2}}{4\pi}\>(h^{2}\gamma_{\kappa}-\gamma_{\rho})\>\int_{V}\>e^{2ik_{1}\hat{\bf{k}}\cdot\bf{x}}\>dV(\bf{x}).$
(2)

For completeness, it is worth mentioning that exists at least one alternative
realization[4] of the DWBA which in the far-field expression carry out the
change $k_{0}\to k_{1}$ in all the places where it appears. This alternative
expression differs from Eq. (2) in the absence of the $h^{2}$ factor affecting
the $\gamma_{\rho}$ term.

###  2.2 TetraScatt formulation

The TetraScatt model calculates the DWBA backscattering due a fluid volume $V$
by discretizing it in a collection $V=\cup_{j=1}^{N}\>T_{j}$ of $N$ tetrahedra
$T_{j}$ and then evaluating the integration defined in (2) as the sum of
integrals over each tetrahedron. Then, the integral of Eq. (2) results in

$f_{\infty}=\sum_{j=1}^{N}\>\frac{k_{1}^{2}}{4\pi}\>(h^{2}\gamma_{\kappa}-\gamma_{\rho})\>\int_{T_{j}}\>e^{2ik_{1}\hat{\bf{k}}\cdot\bf{x}}\>dV(\bf{x}),$
(3)

where $\bf{x}$ are the coordinates of the $T_{j}$ tetrahedron.

The arbitrary chosen vertex $\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D}$ of
each tetrahedron (see Figure 2, where a generic tetrahedron is illustrated)
define oriented segments
$\mathbf{AB}=\mathbf{B}-\mathbf{A},\mathbf{AC}=\mathbf{C}-\mathbf{A},\mathbf{AD}=\mathbf{D}-\mathbf{A}$
which can be used for its parameterization.

Figure 2: Arbitrary tetrahedron of vertex
$\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D}$ which define, with respect to
vertex $\mathbf{A}$, the oriented segments
$\mathbf{AB},\mathbf{AC},\mathbf{AD}$.

Then, integration over the $j$-th tetrahedron with vertex
$\mathbf{A}_{j},\mathbf{B}_{j},\mathbf{C}_{j},\mathbf{D}_{j}$ is carried out
using the parameterization

$\mathbf{x}:R\to
T_{j}\subseteq\mathbb{R}^{3},\quad\mathbf{x}(u,v,w)=\mathbf{A}_{j}+w\>\mathbf{AB}_{j}+v\>\mathbf{AC}_{j}+u\mathbf{AD}_{j},$
(4)

where $R\subset\mathbb{R}^{3}$ is defined as

$R=\\{(u,v,w):0\leq u+v+w\leq 1\quad\text{with}\quad u,v,w\geq 0\\},$ (5)

which describes the j-th tetrahedron referred to the vertex $\mathbf{A}$,
arbitrarily chosen as origin of the local coordinates.

The integration over a single tetrahedron appearing in Eq. (3) in terms of the
parametrization (4) reads

$\int\limits_{T_{j}}\>e^{2ik_{1}\hat{\mathbf{k}}\cdot\mathbf{x}}\>dV(\mathbf{x})=J\int\limits_{R}e^{2ik_{1}\hat{\mathbf{k}}\cdot(\mathbf{A}_{j}+w\mathbf{AB}_{j}+v\mathbf{AC}_{j}+u\mathbf{AD}_{j})}dV=JF(\mathbf{A}_{j},\mathbf{B}_{j},\mathbf{C}_{j},\mathbf{D}_{j},\hat{\bf{k}},k_{1})$
(6)

where $J=J(\mathbf{A}_{j},\mathbf{B}_{j},\mathbf{C}_{j},\mathbf{D}_{j})$ is
the Jacobian accompanying the volume differential, which is constant within
the tetrahedron $T_{j}$, and $F$ is an elementary function based on
polynomials an exponential functions only, whose detailed expression is given
in the Appendix A.

The discretization of a volume $V$ in tetrahedra can be carried out in a CAD
software as SALOME[10] which will result in a volumetric mesh (a structure
containing information about the vertex, the tetrahedra which constitutes the
volume and the plane triangles which build the surface of the enclosed
volume).

Because a tetrahedron has sides which are plane triangles any curvature in the
scattering body surface has to be represented in an approximated way through a
discretization that includes enough elements to ensure that curvature is
captured by them.

In summary, the TetraScatt model takes as input: a mesh representing the
geometry of the scattering body, the physical parameters of the involved media
and the scattering parameters.

##  3 Model validation

We perform two types of verification for the TetraScatt model: (a) To ensure
the accuracy of the integral computation implementation, we compare the
model’s output against exact DWBA integral solutions. This comparison is
conducted for cases where the DWBA integral expressed in Eq. (2) has a closed-
form expression, which happens for example when the scatterer’s geometry is an
sphere or an spheroid. (b) In order to assess the validity of the presented
model as a solution of the acoustic scattering problem for arbitrary shaped
objects we compare our results with a well established benchmark solutions for
penetrable acoustic scattering [11]. Specifically, we utilize the Boundary
Element Method (BEM) solution [12] and the method based on spheroidal wave
functions [13] which can solve the scattering by prolate and oblate spheroids.
While the BEM solution is not constrained by material properties or
geometrical shapes of the scatterers, making it applicable beyond weak
scattering scenarios, it is both computationally and theoretically more
intricate. As a test of rigorous scattering modelling without any assumptions
we need to carry on these tests in the $g,h$ domain where the Born
prescription gives a negligible error, so that any departure can be attributed
to the TetraScatt formulation and not to a breakdown of the Born approach
itself.

The first two subsections (3.1 and 3.2) deal with analytical DWBA solutions
(the (a)-type verifications), while subsection 3.3 presents a comparison with
BEM results for weakly scattering by a copepod (a (b)-type verification). For
all cases, we chose media parameters with values used in Ref. 11, as listed in
the Table 1. The surrounding media properties are those of seawater.

Parameter |  $c$ (m/s) |  $\rho$ (g/cm3) |  $c_{1}$ (m/s) |  $\rho_{1}$ (g/cm3)
---|---|---|---|---
Value | 1477.40 | 1026.80 | 1480.30 | 1028.90
|  |  |  |
Table 1: Material properties (density $\rho$ and sound speed $c$) of the
targets and surrounding medium.

A comment on the nature of the exact benchmark solutions used in this work is
relevant here. The far-field scattering by a fluid sphere, the classical
partial-wave solution of Ref. 14, which also can be found in Ref. 11, is
constituted by a modal series solution whose coefficients has an explicit form
given in terms of spherical Bessel functions. The exact far-field solution for
prolate and oblate spheroids is far more involved than the spherical one since
the computation of the spheroidal wave functions is a more difficult task than
the Bessel function computation. Furthermore, the expansion coefficients has
not analytic expression and are to be founded by solving a family of recursive
matrix systems. Bessel function evaluation is integrated in almost any
computational environment aimed to scientific numerical usage. On the
contrary, spheroidal function evaluation is restricted to a few, specific,
implementations. In this work we adopt the model proposed in Ref. 15, which is
based on arbitrary-precision arithmetic for computing the modal series of
prolate spheroidal wave functions, aimed to a parallel computation environment
of high performance computing. The BEM requires a mesh representation of the
scatterer surface, usually in the form of a collection of $N$ planar
triangles. The solution of the fluid scattering problem case results in a
matrix system whose size is $4N^{2}$. Given that the length $\ell$ of the
segments which constitutes each triangle in this mesh is required to fulfill a
condition $\ell\leq\lambda/5$, where $\lambda$ is the minimum implied
wavelength
($\lambda=\mathrm{Min}\left\\{\lambda_{0}=c_{0}/f,\lambda_{1}=c_{1}/f\right\\})$),
the size $N$ increases with an increase in the frequency $f$.

###  3.1 Sphere

Computing DWBA for a sphere of radius $a$ through Equation (2) gives[3]

$f_{\infty}=k_{1}^{2}a^{3}\>(h^{2}\gamma_{\kappa}-\gamma_{\rho})\>\frac{j_{1}(2k_{1}a)}{2k_{1}a},$
(7)

being $j_{1}$ the spherical Bessel function of first order.

For the comparison with the TetraScatt model, two meshes with $N=350$ and
$N=9521$ tetrahedra representing a sphere of $a=0.01$ m have been used. Two
tests has been conducted. In first place an evaluation of the backscattering
$|f_{\infty}|$ for the frequency range $(1,400)$ kHz was carried out. The
three resulting curves are shown in Fig. 3. The coarse mesh $N=350$ shows
increasing departure from the exact solution of Eq. (7) whereas the finer one
is in good agreement along the entire frequency range. As expected, as long as
the number of tetrahedron in the mesh increases the error respect the exact
solution decreases.

Figure 3: Absolute value $|f_{\infty}|$ for the backscattering by an sphere
under the DWBA according to the exact formula (7) and two executions of
TetraScatt for meshes with the indicated number $N$ of tetrahedra.

In second place, an evaluation of the relative error
$E=|f_{\infty}-f_{\infty}^{\text{ex}}|/|f_{\infty}^{\text{ex}}|$ for the range
$0.9\leq g\leq 1.1$ and $0.9\leq h\leq 1.1$ at fixed frequency $f=$ 300 kHz
was carried out, being the results exhibited in the plots at top of the Fig. 4
as heatmaps. The colorbar scale was saturated at 1 in both plots to facilitate
a direct comparison. The behaviour of error is evidently influenced by $h$ as
indicated by the translational symmetry observed in both plots along the
$g$-direction. There are two critical values $h=0.956$ and $h=1.083$ where
local maxima is attained for all $g$-values. They are associated to nulls
offset as depicted in the plot at the bottom of the figure.

Figure 4: Relative error
$|f_{\infty}-f_{\infty}^{\text{ex}}|/|f_{\infty}^{\text{ex}}|$ for the
backscattering with TetraScatt with respect to the exact formula (7) executed
on the $N=350$ mesh (top left) and on the $N=9521$ mesh (top right) at $f=300$
kHz. Both plots were saturated at 1 to provide a comparison on equal footing.
The two error peaks seen at $h=0.956$ and $h=1.083$ are related to nulls on
the acoustic response as can be verified by the plot (bottom) at $g=0.96$. The
point $g=h=1$ (zero error by definition there) has been removed from both top
plots.

###  3.2 Prolate spheroid

For a prolate spheroid [3] of minor and major semi-axis $a,b$, respectively,
the DWBA has an especially simple expression in the mutually perpendicular
cases of broadside and end-on incidences. In these cases it results in

$f_{\infty}^{\text{DWBA}}=\begin{cases}\displaystyle
k_{1}^{2}a^{3}e\>(h^{2}\gamma_{\kappa}-\gamma_{\rho})\>\frac{j_{1}(2k_{1}a)}{2k_{1}a},\qquad\text{Broadside
incidence}\\\ \\\ \displaystyle
k_{1}^{2}a^{3}e\>(h^{2}\gamma_{\kappa}-\gamma_{\rho})\>\frac{j_{1}(2k_{1}b)}{2k_{1}b},\qquad\text{End-
on incidence}\end{cases}$ (8)

where $e=b/a>1$ is the spheroid aspect-ratio. The sphere solution is recovered
when $e\to 1$ (see Eq. (7)). For the comparison against the model two meshes
representing an spheroid with $a=0.01$ m and $b=0.07$ m of $N=5391$ and
$N=25574$ tetrahedra have been used. The Fig. 5 shows evaluation of the
backscattering $|f_{\infty}|$ for the ranges $(1,400)$ and $(1,100)$ kHz in
broadside and end-on incidences, respectively (inserts at the right of each
figure illustrates the incidence direction).

Figure 5: Absolute value $|f_{\infty}|$ for the backscattering by an spheroid
under the DWBA according to the exact formula (8) and two executions of
TetraScatt with the indicated meshes. Top: broadside incidence. Bottom: end-on
incidence.

The results exhibit qualitatively similar behaviour to the sphere case: an
error that increases with frequency, and finer meshes providing better
agreement than coarser ones. When compared with the sphere, the spheroid
requires meshes with a higher number of tetrahedra for TetraScatt, attributed
to its larger surface area and the strong curvature imposed by the high aspect
ratio $e=7$ (refered as 1:7).

###  3.3 Weakly scattering by a copepod

The two examples previously used as a validation of the model show that if the
number of tetrahedral elements in the volumetric mesh is high enough to ensure
a proper external curvature representation the error can be negligible.

In order to evaluate the TetrasScatt model with a complex scattering object, a
profile representative of a copepod available in the shape catalog inside
ZooScatR package [7] was used. Using this profile as a curved-axis of symmetry
a scattering body can be defined. Since there is not analytical solution for
this kind of geometry, we use the BEM implementation of Ref. 12 as a
benchmark.

For the TetrasScatt model two meshes of $N=2110$ and $N=95291$ were drawn from
the geometry shown in the insert in the plot on Fig. 6 where also are
indicated the incidence direction and a reference scale. Absolute value
$|f_{\infty}|$ as a function of frequency from 1 to 400 kHz is showed in the
main window of the figure for the model and the BEM benchmark. The physical
parameters used were the indicated in Table 1. For the BEM evaluation a mesh
with $N=6608$ triangular elements was enough to evaluate the entire frequency
interval.

Figure 6: Absolute value $|f_{\infty}|$ for the backscattering by a idealized
marine organism: a copepod. The BEM solution was taken as the benchmark
whereas for the TetraScatt execution two meshes of the indicated tetrahedral
elementre were used. The insert in the main window shows the mesh, the
reference system and the incidence direction.

##  4 Validity range of the DWBA for $g$ and $h$

There is general agreement between researchers about the validity of the BA
and the DWBA; when the contrast parameters $g,h$ differ even slightly from
unity the error increases significantly. Further, this behaviour is uneven in
$g$ and $h$ and non-monotonous.

In fisheries acoustics, as well as in other marine applications, sound
scattering by an object is analyzed in the logarithmic scale using the target
strength (TS) parameter which can be expressed as

$\text{TS}=10\log_{10}\left(|f_{\infty}|^{2}\right)\text{ dB re 1 m${}^{2}$}.$
(9)

In addition, since sea measurements are affected by a large number of factors
upon which researchers have little control, TS scattering predictions with a
certain amount of error are allowed. In this vein a more appropriate, and
permissive, quantification for the model error can be defined[11] as

$\varepsilon=\left|\text{TS}^{\text{model}}-\text{TS}^{\text{benchmark}}\right|.$
(10)

In the next sub-sections predictions of the TetraScatt model for the range
$0.9\leq g\leq 1.1$ and $0.9\leq h\leq 1.1$ will be compared with appropriate
benchmark solutions and errors according to the Eq. (10). Two benchmark
solutions of the general, i.e. not weakly, fluid acoustic scattering problem
were used: the farfield scattering by a fluid sphere, the classical partial-
wave solution[14, 11] and the scattering by a prolate fluid spheroid[13] under
the end-on incidence.

The sphere and spheroid keep the given values at sections 3.1 and 3.2,
respectively. To ensure a proper curvature representation in both the sphere
and spheroid, meshes of an appropriate number of tetrahedral elements were
used. This is crucial in the spheroid case since end-on incidence is the
computationally hardest direction when high aspect ratios are considered.

###  4.1 Sphere

Heatmaps of absolute error in TS according to Eq. (10) are shown in the panels
of Fig. 7 for the frequencies $f=$ 38, 120, 200 and 333 kHz in the range
(0.9,1.1) for both $g$ and $h$. The TetraScatt predictions were computed with
a $N=$ 9521 tetrahedra mesh. The central point $(g,h)=(1,1)$ is always
excluded since it represents a no-scattering condition.

Colours were saturated at $\varepsilon=$ 5 dB to clearly define zones of
unacceptable error (more than 5 dB). The pattern is clearly not trivial,
gaining in complexity with increasing frequency and there are curves of high
error related to nulls or resonances in the fluid sphere scattering. In the
context of the DWBA for homogeneous spheres peaks error appears over the curve
$g=g(h)$ in the $(g,h)$ space domain, corresponding to regions where the pre-
factor $(\gamma_{\kappa}h^{2}-\gamma_{\rho})$ vanishes.

Figure 7: Absolute error
$|\text{TS}(f_{\infty})-\text{TS}(f_{\infty}^{\text{ex}})|$ for the
backscattering with TetraScatt executed on a spherical mesh of $N=9521$
elements with respect to the exact solution of a fluid sphere. Results for
$f=$ 38, 120, 200 and 333 kHz are shown in the top left, top right, bottom
left and bottom right panels, respectively. All plots were saturated at 5 dB
to provide a comparison on equal footing.

The translational symmetry in the direction of parameter $g$, present at
$f=38$ kHz is already broken at $f=120$ kHz.

###  4.2 Prolate spheroid in end-on incidence

The previous comparison were made for the sphere which is an isotropic object
with constant curvature. To quantify the curvature influence on the scattering
by the BA, two prolate spheroids of different aspect-ratios $e$ under end-on
incidence were considered.

Heatmaps of absolute error in TS for a prolate spheroid of 1:3 and 1:7 both
evaluated at $f=$ 38 kHz are shown in the left and right panels of Fig. 8,
respectively, in the range (0.9,1.1) for both $g$ and $h$. As observed in the
previous subsection central point $(g,h)=(1,1)$ is excluded. Spheroidal meshes
of $N=$ 25574 for TetrasScatt evaluation of the 1:3 and 1:7 ratio spheroids
were used for both cases. Compared with the $f=$ 38 kHz sphere is noticeable
as the pattern is more complex. The low error region is more restricted now,
there are more regions where the saturated 5 dB error are reached.

Figure 8: Absolute error
$|\text{TS}(f_{\infty})-\text{TS}(f_{\infty}^{\text{ex}})|$ for the end-on
incidence backscattering at $f=$ 38 kHz with TetraScatt executed on a
spheroidal mesh of $N=25574$ elements for 1:3 prolate spheroid (left panel)
and 1:7 spheroid (right panel). All plots were saturated at 5 dB to provide a
comparison on equal footing.

As is known $g$ and $h$ should be close to $1$ to use Born approach but in
practical applications it is worthwhile to investigate to what extents the
approximation remains valid when these parameters deviates from the unitary
value. In the context of acoustic backscattering applied to aquatic ecosystem
research, validity of these approximations for $g,h$ has been previously
discussed. References 3, 4, 11 suggest that the approximations can be applied
for values of both $g$ and $h$ ranging from $1$ to $1.06$, approximately. In
Fig. 8, a solid black line square $[1,1.05]\times[1,1.05]$ is included for
comparison purposes. It can be observed that the error can reach values up to
4 dB there, indicating that this solution may introduce significant errors
when applied to an elongated organism, as the represented by our prolate
spheroid in the end-on incidence.

It can be observed (also in the plot for the sphere case in the Fig. 7) that
regardless of the scatterer geometry or frequency, a significant error in the
diagonal $g+h\approx 2$ appears. The reason for this behavior is explained by
the factor $(\gamma_{\kappa}h^{2}-\gamma_{\rho})$ in Eq. (2). When this
quantity is equal to zero the DWBA predicts $|f_{\infty}|=0$, independently of
the volume integral value. By considering $g=1+\epsilon_{g}$ and
$h=1+\epsilon_{h}$ where $\epsilon_{g},\epsilon_{h}\ll 1$ and neglecting
second-order terms in the Taylor expansion about $g=1,g=1$ is easy to show
that

$\gamma_{\kappa}h^{2}-\gamma_{\rho}\approx-\frac{2}{g}\>(g+h-2),$

so that modeled $f_{\infty}$ and TS values are poorly estimated in the region
$2/g\>(g+h-2)\approx 0$, specially when relative error is considered. It is
important to highlight the difference between the four quadrants in the $g,h$
space: it is fortunate that most of the modelling in the aquatic ecosystem
research correspond to the first quadrant, i.e., parameters $g,h>1$, where the
Born approach for spheres and spheroids at least shows a better behaviour
regarding absolute error.

##  5 Application to Argentinian native species: peisos

The sergestid Peisos petrunkevitchi is a crustacean with abundance on the
Argentinean sea of the Patagonia. There are part of the diet of the hake and
other fish species. The Fig. 9 shows a typical profile of an adult specimen
(top view). Assuming that physical properties of Peisos are in the validity
region of the Born approximation a representative volumetric mesh was built
(see bottom view of Fig. 9) and the backscattering according to TetraScatt and
BEM were calculated.

Figure 9: Typical Peisos profile (top) and representative mesh with simplified
surfaces and legs and tail removed (bottom).

Both the TetraScatt and BEM methods requires meshes of the scatterer’s
geometry, but their idiosyncrasies differ significantly. In TetraScatt, the
mesh is volumetric, and the solution is calculated through simple integral
evaluations. On the other hand, the BEM method employs a surface mesh
representing the body, with the solution calculated from a matrix system
solving procedure. The elements of this matrix are the outcome of integrals
computed using quadrature methods. BEM scattering calculations are
significantly more time-consuming and require greater computational resources.
However, the advantage lies in its ability to handle various contrast values,
unlike the Born approximation.

Since the $\rho,c$ values appropriate for Peisos are currently not known, the
calculations were carried out by using the Table 1 values. The backscattering
TS vs frequency pattern is shown in Fig. 10, where a good agreement between
TetraScatt with a mesh of $16938$ tetrahedral elements and BEM with a mesh of
$5306$ planar triangular elements is achieved.

Figure 10: TS value for the backscattering by the mesh of figure 9 (bottom)
representing a typical Peiso organism. The incidence is illustrated within the
insert on the figure.

For a quick comparison in computation resources used by BEM and TetraScatt it
is interesting comment about the computation times. In this case, for the
frequency interval considered in the Figure 10, which correspond to 389
values, TetraScatt took 3 seconds running in a single core of an Intel
i7-12700 whereas BEM took one hour for the same task in a parallel execution
on a 130 cores mixed cluster environment (Xeon E5-2620 and Xeon Gold 6240R).

##  6 Conclusions

A model based on the DWBA that describes the backscattering of weakly
scattering objects was developed, computationally implemented in the R
language and uploaded to a public repository. The presented TetraScatt model
has been validated against benchmark solutions based on analytical integration
and on general scattering methods (BEM and a previously developed scattering
spheroidal code based on spheroidal wave function calculation, which are not
restricted to a weakly condition). A study on the range of validity of the
DWBA with respect to the contrast parameters $g,h$ in terms of an error in dB
showed a complicated pattern which hampers to obtain an easy prescription for
the safe application of the model regarding an error threshold.

Since the integration is exact in each tetrahedron any error in the scattering
response (assuming that the weakly condition is fulfilled) must come from the
inaccurate representation of the surface curvature of the bodies whose
scattering is being calculated. In any case larger volumetric meshes required
to overcome this curvature issue doesn’t handicap the method noticeably
because the computational task to be done (integration reduces to a formula)
is elementary and its complexity doesn’t rise with an increase of frequency.
Instead, BEM formulation requires the solving of a matrix system whose size is
related to the number of elements of the mesh and the elements themselves are
evaluated through integration by quadrature rules which are increasingly
expensive with the frequency. Off course, BEM is an exact method in the sense
that no approximations or limitations are assumed in the scattering problem
besides a proper relationship between the triangle’s segment length and the
wavelength of the incident field.

The TetraScatt only requires a volumetric mesh in the Gamma Mesh Format (GMF)
to be able to compute the scattering without any kind of limitations about the
geometrical shape of the object. This last feature represent an add to the
tools which are available to the underwater acoustic community and it is based
on remarkably simple ideas on modelling and requires only conventional
computational resources to be executed.

## References

  * [1] E. Lavia, J. D. Gonzalez and S. Blanc, Modelling High-Frequency Backscattering from a Mesh of Curved Surfaces Using Kirchhoff Approximation, Journal of Theoretical and Computational Acoustics 26 (2018) 1–17.
  * [2] P. M. Morse and K. U. Ingard, Theoretical Acoustics (Princeton University Press, 1986).
  * [3] D. Chu and Z. Ye, A phase-compensated distorted wave born approximation representation of the bistatic sscattering by weakly scattering objects: Application to zooplankton, The Journal of the Acoustical Society of America 106 (1999) 1732–1743.
  * [4] B. A. Jones, A. C. Lavery and T. K. Stanton, Use of the distorted wave born approximation to predict scattering by inhomogeneous objects: Application to squid, The Journal of the Acoustical Society of America 125 (2009) 73–88.
  * [5] T. K. Stanton, D. Chu, P. H. Wiebe and C. S. Clay, Average echoes from randomly oriented random-length finite cylinders: Zooplankton models, The Journal of the Acoustical Society of America 94 (1993) 3463–3472.
  * [6] T. K. Stanton and D. Chu, Review and recommendations for the modelling of acoustic scattering by fluid-like elongated zooplankton: euphausiids and copepods, ICES Journal of Marine Science 57 (2000) 793–807.
  * [7] S. Gastauer, D. Chu and M. J. Cox, Zooscatr-an r package for modelling the scattering properties of weak scattering targets using the distorted wave born approximation, The Journal of the Acoustical Society of America 145 (2019) EL102–EL108.
  * [8] E. H. Pees, An efficient method for computing backscattering from born objects of arbitrary shape, The Journal of the Acoustical Society of America 129 (2011) 3470–3474.
  * [9] J. D. Gonzalez, E. F. Lavia and G. Cascallares, tetrascatt: Acoustic Scattering for Complex Shapes by Using the DWBA, 2023, r package version 0.1.0.
  * [10] A. Ribes and C. Caremoli, SALOME platform component model for numerical simulation, in Computer Software and Applications Conference, 2007. COMPSAC 2007. 31st Annual International (IEEE, 2007), volume 2, pp. 553–564.
  * [11] J. M. Jech, J. K. Horne, D. Chu, D. A. Demer, D. T. Francis, N. Gorska, B. Jones, A. C. Lavery, T. K. Stanton and G. J. Macaulay, Comparisons among ten models of acoustic backscattering used in aquatic ecosystem research, The Journal of the Acoustical Society of America 138 (2015) 3742–3764.
  * [12] J. D. Gonzalez, E. F. Lavia, S. Blanc, M. Maas and A. Madirolas, Boundary element method to analyze acoustic scattering from a coupled swimbladder-fish body configuration, Journal of Sound and Vibration 486 (2020) 115609.
  * [13] J. D. González, E. F. Lavia and S. Blanc, A computational method to calculate the exact solution for acoustic scattering by fluid spheroids, Acta Acustica united with Acustica 102 (2016) 1061–1071.
  * [14] V. C. Anderson, Sound scattering from a fluid sphere, The Journal of the Acoustical Society of America 22 (1950) 426–431.
  * [15] E. F. Lavia, Acoustic scattering by two fluid confocal prolate spheroids, European Journal of Mechanics-B/Fluids 88 (2021) 47–57.

##  Appendix A Function $F$

Considering the following notational conventions

$\displaystyle kab=k_{1}\hat{\bf{k}}\cdot\mathbf{AB},\quad
kac=k_{1}\hat{\bf{k}}\cdot\mathbf{AC},\quad
kad=k_{1}\hat{\bf{k}}\cdot\mathbf{AD}$ $\displaystyle
kbc=k_{1}\hat{\bf{k}}\cdot\mathbf{BC},\quad
kbd=k_{1}\hat{\bf{k}}\cdot\mathbf{BD},\quad
kdc=k_{1}\hat{\bf{k}}\cdot\mathbf{DC},$ $\displaystyle
ka=k_{1}\hat{\bf{k}}\cdot\mathbf{A},\quad
kb=k_{1}\hat{\bf{k}}\cdot\mathbf{B},\quad
kc=k_{1}\hat{\bf{k}}\cdot\mathbf{C},\quad
kd=k_{1}\hat{\bf{k}}\cdot\mathbf{D},$

where it is supposed that $kad\neq 0$, (otherwise it is always possible
rearrange the four vertices in a way that this condition is fullfiled), then
the function
$F(\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D},\hat{\bf{k}},k_{1})$ which is
results from the integration prescribed by the Eq. (6) is given by

$F=\begin{cases}\displaystyle\>-\frac{1}{4kab}\left[\frac{e^{i(kd+kc)}\mathrm{sinc}(kdc)-e^{i(kb+kc)}\mathrm{sinc}(kbc)}{kbd}\>+\right.\\\
\\\ \displaystyle\hskip
140.00021pt\left.\frac{e^{i(ka+kc)}\mathrm{sinc}(kac)-e^{i(kd+kc)}\mathrm{sinc}(kdc)}{kad}\right],\qquad|kbd|\geq\varepsilon\\\
\\\
\displaystyle\>-\frac{1}{4kab}\left[e^{i(kb+kc)}\>\Gamma(kbd,kbc)\>+\>\frac{e^{i(ka+kc)}\mathrm{sinc}(kac)-e^{i(kd+kc)}\mathrm{sinc}(kdc)}{kad}\right],\qquad|kbd|<\varepsilon\\\
\end{cases}$

where $\varepsilon=1$E-8 is a fixed parameter and the $\Gamma(x,y)$ function
is defined according to

$\Gamma(x,y)=\frac{e^{ix}\mathrm{sinc}(y-x)-\mathrm{sinc}(y)}{x}.$

To ensure proper behaviour when $x,y\leq\varepsilon$ appropriate Taylor
expansions are used in the numerical implementation. The sinc function is
defined as $\mathrm{sinc}(x)=\sin(x)/x$.


# Near-optimal Deep Reinforcement Learning Policies from Data for Zone
Temperature Control

Loris Di Natale, Bratislav Svetozarevic, Philipp Heer, and Colin N. Jones This
research was supported by the Swiss National Science Foundation under NCCR
Automation, grant agreement 51NF40_180545.L. Di Natale, B. Svetozarevic, P.
Heer are with the Urban Energy Systems Laboratory, Swiss Federal Laboratories
for Materials Science and Technology (Empa), 8600 Dübendorf, Switzerland
{loris.dinatale, bratislav.svetozarevic<EMAIL_ADDRESS>Di Natale
and C.N. Jones are with the Laboratoire d’Automatique, Swiss Federal Institute
of Technology Lausanne (EPFL), 1015 Lausanne, Switzerland<EMAIL_ADDRESS>
###### Abstract

Replacing poorly performing existing controllers with smarter solutions will
decrease the energy intensity of the building sector. Recently, controllers
based on Deep Reinforcement Learning (DRL) have been shown to be more
effective than conventional baselines. However, since the optimal solution is
usually unknown, it is still unclear if DRL agents are attaining near-optimal
performance in general or if there is still a large gap to bridge.

In this paper, we investigate the performance of DRL agents compared to the
theoretically optimal solution. To that end, we leverage Physically Consistent
Neural Networks (PCNNs) as simulation environments, for which optimal control
inputs are easy to compute. Furthermore, PCNNs solely rely on data to be
trained, avoiding the difficult physics-based modeling phase, while retaining
physical consistency. Our results hint that DRL agents not only clearly
outperform conventional rule-based controllers, they furthermore attain near-
optimal performance.

##  I Introduction

In the European Union, buildings are responsible for $36\%$ of the energy-
related greenhouse gas emissions, consuming $40\%$ of the end-use energy [1].
Almost two-thirds of the latter can be linked to space heating and cooling
[2], which calls for solutions to decrease the energy intensity of the sector.
Building occupants play a major role through the generated heat gains and the
constraints that arise from their comfort requirements [3]. The latter
transform the building control problem into a multi-objective optimization,
where one tries to minimize energy consumption while maximizing the comfort of
the users. To attain this objective, one solution is to replace poorly
performing conventional zone temperature Rule-Based Controllers (RBCs) with
more advanced methods, such as Deep Reinforcement Learning (DRL) policies [4].
Indeed, numerous studies have already shown how DRL agents are able to improve
upon the performance of classical RBCs in terms of energy consumption, cost
savings, and/or comfort satisfaction, e.g. in [5, 6, 7, 8, 9].

Because of the slow thermal dynamics of buildings and the data-inefficiency of
DRL algorithms, months of data are routinely needed to converge to a good
solution, e.g. in [6, 10]. This is infeasible in real case studies and
researchers thus have to rely on simulation environments. Many of them, e.g.
[6, 8, 9], are trained in complex physics-based simulation environments, such
as EnergyPlus [11]. Setting up such environments can however be a very
difficult and lengthy procedure [12], and it furthermore results in nonlinear
models. On the other hand, only historical data is required to fit black-box
models, such as in [5], but one cannot ensure that such models are physically
accurate and they remain highly nonlinear. Consequently, studies relying on
such complex models usually do not provide upper bounds on the performance of
DRL agents since the theoretical optimal solution (i.e. knowing the noise and
disturbances) is often impossible to compute.

To the best of the authors’ knowledge, works comparing DRL agents to the
optimal performance rely on relatively simple models built from first
principles, such as [10, 13]. In this paper, we consider predefined dynamic
comfort bounds, which is an improvement over previously studied fixed bounds.
Moreover, [10] trains discrete DRL agents and [13] assumes access to solar
photovoltaic electricity generation and battery storage. Furthermore, in this
paper, we use Physically Consistent Neural Networks (PCNNs) [14] as simulation
environments instead of first principles models. PCNNs are trained from
historical data, hence bypassing the complex design stage of physics-based
methods, but remain physically consistent with respect to control inputs.
However, they introduce nonlinearities in the model through the use of Neural
Networks (NNs). Nonetheless, we can compute super-optimal control inputs, i.e.
assuming access to a perfect oracle of the future, with Linear Programs (LPs),
which are easy to solve. This allows us to assess the performance of DRL
policies in various settings in feasible time.

To summarize our contribution, we analyze the performance of DRL agents for
zone temperature control compared to industrial rule-based baselines and to
the super-optimal solution, using PCNNs as simulation environments. Our
pipeline solely relies on past historical data to fit the PCNNs, avoiding the
difficult physics-based modeling procedure while retaining physical
consistency. The agents are trained in continuous action spaces to maintain
the temperature of the zone between predefined dynamic comfort bounds. We then
investigate the sensitivity of the obtained policies to the choice of random
seed and reward function, analyzing the gap with the theoretically optimal
solutions in each setting.

The rest of the paper is structured as follows. Section II briefly introduces
DRL, PCNNs, and the case study. We then detail our methods in Section III,
present and discuss the results in Section IV, and conclude the paper with
Section V.

##  II Background

###  II-A Deep Reinforcement Learning

Reinforcement Learning (RL) problems are usually formulated as Markov Decision
Processes (MDPs), which are represented by tuples $<S,A,P,\rho_{0},r,\gamma>$,
where $S$ is the state space, $A$ the action space, $P=P(s^{\prime}|s,a)$ the
probability of transitioning from state $s\in S$ to $s^{\prime}\in S$ when
action $a\in A$ is taken, $\rho_{0}$ the initial state distribution,
$r=r(s,a)$ the reward function, and $0<\gamma<1$ the discount factor. At each
time step $t$, the agent observes the state $s_{t}$ and chooses an action
$a_{t}$. The environment then steps to the next state $s_{t+1}$, provides the
reward $r(s_{t},a_{t})$ and $s_{t+1}$ to the agent, and the process repeats
itself until the end of an episode after some time $H$.

The goal of RL algorithms is to find the policy $\pi(a|s)$ that maximizes the
expected return $R_{t_{0}}$, i.e. the expected sum of discounted rewards, over
an episode starting at time $t_{0}$:

$\displaystyle
J(\pi)=\mathbb{E}_{s_{t_{0}}\sim\rho_{0},a\sim\pi}[R_{t_{0}}]=\mathbb{E}_{s_{t_{0}}\sim\rho_{0},a\sim\pi}[\sum_{t=t_{0}}^{t_{0}+H}\gamma^{t-t_{0}}r(s_{t},a_{t})]$
(1)

To that end, many algorithms rely on learning an approximation of this
objective, named the Q-function, which estimates the expected return the agent
will receive if it takes action $a$ in state $s$ and then follows the current
policy $\pi$:

$\displaystyle Q^{\pi}(s,a)=E_{a\sim\pi}[R_{t_{0}}|s_{t_{0}}=s,a_{t_{0}}=a]$
(2)

The Q-function is often referred to as the critic, while the policy function
is called the actor. In DRL, one parametrizes both the critic and the actor
with NNs, i.e. $\pi_{\theta}(a|s)$ and $Q^{\pi}_{\phi}(s,a)$, and optimizes
over $\\{\theta,\phi\\}$. Among the countless improvements and techniques
presented in various contributions, we want to point out the influence of
target networks [15], which are used in this work. The idea is to keep a copy
of the actor and the critic in memory and only update these copies slowly to
decrease the usual overestimation bias of Q-values and stabilize the learning
process.

###  II-B Physically Consistent Neural Networks

While classical NNs have been proven to be very effective function
approximators in various fields, they remain agnostic to the underlying
physical laws when modeling physical systems and might fail to grasp
fundamental principles [14]. For example, when modeling the temperature in a
zone, turning the heating on should always have the effect of increasing the
temperature, as DRL agents could otherwise learn spurious behaviors. To still
leverage the expressiveness of NNs while retaining physical consistency with
respect to some inputs, PCNNs were proposed and analyzed in [14], and we refer
the reader to the original paper for the details. The main feature of PCNNs is
that they are physically consistent with respect to power inputs, external
temperatures, and temperatures in neighboring zones by design. In particular,
this means that applying heating or cooling power to the zone will have the
expected impact on the change in temperature, ensuring that DRL agents learn
meaningful behaviors.

In brief, the core idea of PCNNs is to treat part of the inputs – the ones
that require physical consistency – in parallel of the main recurrent NN
pipeline, in a linear module $E$. Heavily inspired from classical physics-
based resistance-capacitance models to ensure consistency, this module
accumulates or dissipates energy at each time step to capture the impact of
the heating/cooling power $u$, the outside temperature $T^{out}$, and the
temperature in the neighboring zone $T^{neigh}$ on the temperature of the
modeled zone $T$. Simultaneously, a recursive NN pipeline $D$ processes the
other inputs $x$ to capture unforced dynamics, i.e. the evolution of the
temperature if no power is applied and heat losses are neglected. The final
temperature prediction $T$ is then the sum of both components, the unforced
dynamics and the accumulated energy since the beginning of the prediction.
Altogether, we can write PCNN predictions as follows:

$\displaystyle D_{k+1}$ $\displaystyle=D_{k}+f(D_{k},x_{k})$ (3)
$\displaystyle E_{k+1}$
$\displaystyle=E_{k}+au_{k}-b(T_{k}-T^{out}_{k})-c(T_{k}-T^{neigh}_{k})$ (4)
$\displaystyle T_{k+1}$ $\displaystyle=D_{k+1}+E_{k+1},$ (5)

where $f$ is any nonlinear function, typically comprised of recurrent NNs, and
$a$, $b$, $c$ are parameters learned simultaneously with $f$ through
backpropagation. While this represents the situation during the heating
season, the cooling case is analogous, using another parameter $d$ instead of
$a$ to capture the efficiency of the cooling system instead of the heating
one. Note that we define heating as positive and cooling as negative power
throughout this paper, so that cooling decreases the energy accumulated in
$E$, as expected.

###  II-C Case study

The data used to train the PCNNs and DRL agents is taken from the Urban Mining
and Recycling (UMAR) unit, an apartment in the NEST building, located in
Duebendorf, Switzerland, and pictured in Fig. 1 [16]. It consists of two
bedrooms separated by a living room, all powered by heating/cooling panels
letting hot/cold water flow through their ceiling. More than three years of
data is available and was processed as explained in [14], including the
normalization procedure.

Figure 1: NEST building, Duebendorf, and the UMAR unit circled in white ©
Zooey Braun, Stuttgart.

##  III Methods

###  III-A Problem setting

The task of the agents is to control the heating or cooling power of one
bedroom in UMAR at each time step of $15\text{\,}\mathrm{min}$, i.e.
$a_{t}=u_{t}$, to minimize the energy consumption while respecting predefined
dynamic comfort bounds. While we require the temperature to stay between
$23\text{\,}\mathrm{\SIUnitSymbolCelsius}$ and
$24\text{\,}\mathrm{\SIUnitSymbolCelsius}$ during the night, i.e. from
$20\text{\,}\mathrm{h}$ to $8\text{\,}\mathrm{h}$, it can be relaxed during
the day, when the bedroom is generally unoccupied. Specifically, we relax the
lower bound to $21\text{\,}\mathrm{\SIUnitSymbolCelsius}$ during the day in
the heating season, and the upper bound is raised to
$26\text{\,}\mathrm{\SIUnitSymbolCelsius}$ from $8\text{\,}\mathrm{h}$ to
$20\text{\,}\mathrm{h}$ during the cooling season.

The state $s$ observed by the agents at each time step is composed of the
input features of the PCNN, i.e. zone temperatures (of the controlled and
neighboring room), ambient conditions (temperature and solar irradiation), and
time information (sine and cosine functions of the month of the year and time
of the day, and the day of the week) [14]. Additionally, agents know whether
they can heat or cool and the current temperature comfort bounds for the
occupants. To have more expressive policies aware of the evolution of the
environment in time, we also add $12$ autoregressive terms, with
$15\text{\,}\mathrm{min}$ time steps, of the zone temperatures and ambient
conditions, so that agents know the state of these variables during the last
three hours when taking decisions.

To train and evaluate the agents, we create $12\text{\,}\mathrm{h}$ to
$72\text{\,}\mathrm{h}$-long sequences of data with no missing values,
hereafter referred to as trajectories. These sequences, possibly overlapping
each hour, are then separated into a training and a validation data set. To
learn policies that are robust to measurement noise, we add a small
independent Gaussian noise on top of the zone temperature measurement.

The reward function is defined as the negative weighted sum of comfort
violations, i.e. how far from the designed bounds the temperature inside the
zone is, and energy consumption. Mathematically, we thus have:

$\displaystyle r(s_{t},a_{t})$
$\displaystyle=-\max{\\{L_{t}-T_{t},0\\}}-\max{\\{T_{t}-U_{t},0\\}}$
$\displaystyle\ \quad-\lambda|a_{t}|,$ (6)

where $L$ and $U$ represent the lower and upper comfort bounds, respectively,
and $\lambda$ the weighting factor. As a rule of thumb, we designed $\lambda$
such that the agent receives the same penalty for using a power of
$1\text{\,}\mathrm{kW}$ and for being
$0.5\text{\,}\mathrm{\SIUnitSymbolCelsius}$ outside of the comfort bounds.

###  III-B Algorithm

In this work, we rely on the Twin Delayed Deep Deterministic policy gradient
algorithm (TD3) [17], an actor-critic method improving upon the classical Deep
Deterministic Policy Gradient algorithm (DDPG) [18]. The main issue with DDPG
is that the critic network has a tendency to sharply overestimate Q-values,
which is then exploited by agents and can lead to poorly performing policies
[17]. To correct this bias, TD3 introduces three main improvements over DDPG:

  * •

Inspired from the success of Double Q-learning [19], two critic networks are
learned in parallel, and the smallest of the two approximated Q-values is used
to limit overestimation.

  * •

To avoid instability arising from fast-changing Q-functions, the actor and the
targets networks are updated less frequently than the critics.

  * •

To reduce the ability of the policy to exploit overestimated Q-values, noise
is added to the action chosen by the actor before it is evaluated by the
critics.

In our implementations, all the NNs have three hidden layers of $512$ neurons
and we use a slightly modified version of the Adam optimizer [20] with a
learning rate of $10^{-4}$.

###  III-C Performance assessment

To analyze the performance of the DRL agents, we compare them to two classical
RBC baselines and the optimal solution. Baseline 1 is tracking a reference
$0.5\text{\,}\mathrm{\SIUnitSymbolCelsius}$ away from the bound, turning the
heating on and off as soon as the target temperature is met. Baseline 2 is a
classical rule-based controller with a one-degree hysteresis, i.e. it starts
heating at full power when the temperature reaches the lower bound and until
one degree has been gained. During the cooling season, similar strategies are
used, with the controllers starting to cool once the upper bound or the
reference temperature half a degree from it is reached.

To compute the super-optimal control inputs over an episode, knowing all the
external conditions and the measurement noise $N$ over the horizon, we need to
solve the following LP from time $t_{0}$, where the objective function is
designed to match the reward of the agents:

$\displaystyle\min_{u_{0},...,u_{H-1}}\quad\sum_{k=0}^{H-1}$
$\displaystyle\tilde{\lambda}u_{k}+\epsilon^{L}_{k+1}+\epsilon^{U}_{k+1}$ (7)
$\displaystyle s.t.\qquad\ E_{0}$ $\displaystyle=E(t_{0})$ (8) $\displaystyle
E_{k+1}$ $\displaystyle=E_{k}+au_{k}-b(T_{k}-T^{out}_{k})$
$\displaystyle\qquad\quad-c(T_{k}-T^{neigh}_{k})$ (9) $\displaystyle T_{k+1}$
$\displaystyle=D(t_{k+1})+E_{k+1}+N(t_{k+1})$ (10) $\displaystyle-T_{k+1}$
$\displaystyle\leq-L_{k+1}+\epsilon^{L}_{k+1}$ (11) $\displaystyle T_{k+1}$
$\displaystyle\leq U_{k+1}+\epsilon^{U}_{k+1}$ (12) $\displaystyle u^{low}$
$\displaystyle\leq u_{k}\leq u^{high},$ (13)

where all the constraints hold for $k=0,...,H-1$. Note that we use
$\tilde{\lambda}$ here, which equals $\lambda$ in heating cases and $-\lambda$
during the cooling season, so that it always penalizes the absolute value of
the power used, as in the reward function (6). The key property of PCNNs
rendering this optimization feasible is the fact that the highly nonlinear
unforced dynamics $D$ are independent of the control inputs $u$, which appear
linearly in the temperature predictions of PCNNs in (3)-(5). $D$ can hence be
computed a priori for the entire horizon (knowing all external conditions) and
fed into the optimization procedure as an external variable to compute the
temperature evolution in (10). This LP can be solved very efficiently with
common tools, which allows us to compute the optimal solution for thousands of
trajectories in feasible time. The code and data can be found here:
https://gitlab.nccr-automation.ch/loris.dinatale/NoDRL.

##  IV Results and Discussion

###  IV-A Performance analysis

Since DRL policies are notoriously sensitive to the random seed used [21], we
trained $10$ agents with different seeds to compare their performance. In Fig.
2, we plot the convergence rate of the agents in green, with the median in
bold, and zoom in the bottom plot to get a clearer picture. Here, one epoch
represents $5000$ time steps in the environment, i.e. slightly over $50$ days
of data, after which the agent is evaluated on $50$ trajectories from the
validation set. The maximal performance attained by each agent of the training
horizon is reported in dashed lines, and we included the baseline and optimal
performance for reference. One can see the reward obtained by each agent
fluctuating a lot, and most agents attained a maximal reward between $-2.25$
and $-2.3$. All DRL agents present comparable learning patterns: they converge
to policies that perform similarly to the baselines after $35$-$40$ epochs and
then consistently outperform them, with a decrease in the variance of the
obtained rewards. However, we can see a few exceptions along the training
pattern, with some agents’ performance plummetting for a few epochs but
recovering quickly. It is also noteworthy that one of the agents performed
significantly worse than the other ones, only reaching a best reward of
$-2.51$. This proves that DRL agents are also sensitive to the choice of
random seed in our setting, and one should hence always perform several runs
to rule out the possibility of the random seed impacting the quality of the
results.

Figure 2: Convergence rate of DRL agents with different random seeds in green,
with the median in bold and the maximal reward attained by each agent in
dashed lines. For reference, the baseline and optimal performance are also
shown in dashed lines.

For each agent, we selected the best policy obtained during the training
phase, i.e. the one that achieved the maximum reward on the $50$ test
trajectories, and analyzed its average performance on almost $2000$ $3$-day
long trajectories from the validation set in terms of rewards, energy
consumption, and comfort violations. In the top plot of Fig. 3, one can
observe the median results achieved by the different random seeds, where this
metric was preferred over the mean to decrease the impact of the poorly
performing seed. Note that the random seed also impacts the performance of the
baselines and the optimal solution because it changes the noise added to the
temperature measurements. However, very little variance could be observed
between different runs, the influence of the Gaussian noise cancels out over
the entire validation data set and was hence not reported. The bottom plot
shows the performance gap between the optimal one and the other controllers.
Note that we subtracted unavoidable penalties from all statistics: since each
episode is initialized from a situation observed in the historical data, the
temperature in the zone might be initialized out of the defined comfort
bounds, leading to unavoidable penalties for any controller (see Fig. 4). We
removed them to give a better picture of the difference between each
controller.

Figure 3: Comparison between the baselines, the average performance of three
DRL agents, and the optimal solution on the validation set. The top plot
presents the average performance of each controller after subtraction of the
unavoidable penalties, while the bottom one show the gap in performance with
respect to the optimal trajectories for each controller.

Remarkably, with the chosen parameters, agents were able to converge to a
near-optimal solution. They found a relatively similar trade-off between
energy consumption and comfort violations as the optimal solution, consuming
roughly the same amount of energy as the optimal solution at the cost of
slightly increased comfort violations of a little over $20\%$. This is,
however, still much better than what the baselines are able to do, both in
terms of energy savings and comfort improvements. Alltogether, these results
confirm that DRL agents are able to converge to policies that not only clearly
beat classical controllers but actually attain near-optimal performance.

Figure 4: Example of the behavior of a DRL agent over $3$ days, compared to
the baselines and the optimal solution. The top plot shows the room
temperature, the bottom one the power input.

In general, all agents were able to capture the expected and desired
preheating/precooling behavior. Indeed, they always take action earlier than
the other controllers, especially the rule-based ones, to anticipate
constraint tightenings, as for example pictured in Fig. 4 (without Gaussian
noise for clarity). They usually heat or cool the zone until the temperature
is slightly above or below the comfort bound and then stop. This often results
in the temperature reaching the bound again just before the constraints are
relaxed. Overall, this strategy is not far from the optimal one, which starts
to act later before the constraints tighten and uses the full available power
to reach the narrow comfort bounds faster. It then takes advantage of its full
knowledge of the environment to input just enough energy in the system for the
temperature to stay exactly at the desired limit and avoid comfort penalties
unless, for example, solar heat gains are expected to create more violations
later (e.g. Fig. 4). Due to their reactive nature, both baselines cannot
anticipate constraint tightenings and relaxations, acting too fast or too
late. Furthermore, they do not anticipate any external input, such as the
solar irradiation, for example leading to overheating behavior.

In the particular case of Fig. 4, the agent consumed slightly less energy than
the optimal solution, $10.29\text{\,}\mathrm{kW}\text{\,}\mathrm{h}$ against
$10.66\text{\,}\mathrm{kW}\text{\,}\mathrm{h}$, but had slightly more comfort
violations, $13.2\text{\,}\mathrm{K}\text{\,}\mathrm{h}$ against
$12.4\text{\,}\mathrm{K}\text{\,}\mathrm{h}$. The difference can be noticed
when the comfort bounds were relaxed: the agent stopped heating slightly
earlier to save energy and anticipate potential solar gains, while the optimal
solution kept heating to stay longer above the lower bound. In general, DRL
agents have a tendency to converge to risk-averse polices because of the noisy
observations returned by the environment. Indeed, they usually slightly
overheat/overcool the zone when the constraints are tightened to end up
further away from the bounds and hence avoid comfort penalties arising from
the noisy temperature measurement jumping outside of the comfort zone. On the
other hand, when computing the optimal trajectory, the noise is assumed to be
known, such that the optimization knows when to input power to the zone to
always stay in bounds despite the impact of the noise.

###  IV-B Sensitivity to the weighting factor

To assess the impact of the weighting factor $\lambda$ on the policy found by
the DRL agents, we performed a sensitivity analysis to investigate the range
of values that lead to near-optimal solutions and the different trade-offs
reached by the different policies. To that end, we trained agents in the same
environment, with a fixed random seed, but multiplying or dividing $\lambda$
by increasing powers of two to reflect situations where more and more
importance is put on decreasing the energy consumption or the amount of
comfort violations, respectively. We can then analyze the Pareto frontiers
arising from the trade-off between both objectives for the DRL policies and
the optimal solutions in Fig. 5, where the weighting factor was decreased from
$4\lambda$ to $\frac{1}{16}\lambda$ from the left to the right along the
dashed Pareto frontiers. The performance of the baselines was also plotted for
reference.

Figure 5: Sensitivity analysis of DRL agents and the optimal solution to
different weighting factors in (6), decreasing it from left to right along the
dashed lines. The trade-off obtained by both baselines is also plotted for
reference, as well as the results obtained with different seeds in Section
IV-A in shaded green.

As can be observed, the agent is generally able to strike a trade-off similar
to the optimal one, usually consuming roughly the same amount of energy at the
cost of slightly more comfort violations, as long as the weighting factor is
not too big. Once $\lambda$ is multiplied by $4$, the agent struggles to find
an interesting solution: it uses very little energy but does not improve the
comfort of the occupants much compared to the baselines (top left green
point). When we tried to increase $\lambda$ by a factor of $8$, the agent
quickly converged to a very poor policy that never uses energy at all. On the
other hand, decreasing the weighting factor impacts the quality of the
solution less: policies consume slightly more energy each time, slowly
reducing the amount of comfort violations.

In this analysis, it is important to remember that the choices of random seed
again impacts the results, which could explain why the trade-off obtained by
the agent with a weighting factor of $\frac{1}{8}\lambda$ is slightly higher
than the Pareto front. This is confirmed by the shaded green markers showing
the performance of the agents discussed in Section IV-A, where we see the
impact of different random seeds for the same weighting factor. All random
seeds usually lead to similar solutions, but we can also observe the outlier
pointed out previously (top left shaded green marker). Interestingly, even
though it obtained worse rewards than all the other agents, we can see its
behavior still lying near the Pareto front. Hence, it would seem that this
agent simply converged to a different solution that was not optimal in this
situation but might be the expected behavior under different circumstances.
Note that the choice of random seed also impacts the performance of the
baselines and the optimal solution, but we again found that the differences
were small and hence only plotted one solution for clarity.

To summarize, one should be careful with the design of the weighting factor,
as it might impact the quality of the solution. However, our tests showed that
a wide range of values, from $2\lambda$ to $\frac{1}{16}\lambda$ in our case,
could be selected, reflecting different preferences for the occupants, and
still lead to near-optimal behaviors. In general, choosing a value that is too
large seems to be more problematic than the contrary. Finally, while the
random seed seems to impact in which region of the Pareto frontier the
solution converges, our experiments provide evidence that DRL agents are
always converging to a region near the Pareto frontier of all DRL policies.
However, the latter is diverging from the optimal Pareto frontier when agents
try to use too little energy, and it thus seems safer to choose small values
for $\lambda$.

##  V Conclusion

In this work, leveraging the expressiveness and physical consistency of PCNNs,
we trained DRL agents to control the temperature of a building zone, balancing
energy consumption and comfort violations. Thanks to the structure of PCNNs,
we were able to compute the theoretical optimal control inputs over each
episode, which showed that DRL agents were not only clearly outperforming
rule-based baselines, they were also able to obtain near-optimal policies. We
analyzed the impact of the random seed and of the weighting factor balancing
energy consumption and comfort violations on the obtained solutions. The
former was shown to sometimes lead to lower quality solutions in terms of the
obtained rewards, but still lying near the Pareto frontier of both objectives.
For the latter, we found a large range of possible values leading to near-
optimal solutions, until the factor is too big.

Interestingly, the entire pipeline to design and fit PCNNs and DRL agents only
relies on historical data, avoiding the difficult physics-based design phase.
However, DRL agents still take years of artificial data to converge to
interesting policies. In future works, we thus want to address this data-
inefficiency issue and create faster converging DRL agents through physics-
inspired structures. We also plan to scale up our analysis to larger and more
complex systems, and to compare DRL agents to Model Predictive Controllers,
both in simulation but also on real buildings.

##  VI Acknowledgments

The authors gratefully acknowledge the support of the Swiss National Science
Foundation.

## References

  * [1] European Commission (EC), “Factsheet: The energy performance of buildings directive,” 2019, accessed: 01.02.2022. [Online]. Available: https://ec.europa.eu/energy/sites/ener/files/documents/buildings\\_performance\\_factsheet.pdf
  * [2] Eurostat, statistics explained, “Energy consumption in households,” 2020, accessed: 2021.05.28. [Online]. Available: https://ec.europa.eu/eurostat/statistics-explained/index.php?title=Energy\\_consumption\\_in\\_households
  * [3] F. Oldewurtel, D. Sturzenegger, and M. Morari, “Importance of occupancy information for building climate control,” _Applied energy_ , vol. 101, pp. 521–532, 2013.
  * [4] B. Svetozarevic, C. Baumann, S. Muntwiler, L. Di Natale, M. N. Zeilinger, and P. Heer, “Data-driven control of room temperature and bidirectional EV charging using deep reinforcement learning: Simulations and experiments,” _Applied Energy_ , vol. 307, p. 118127, 2022.
  * [5] Z. Zou, X. Yu, and S. Ergan, “Towards optimal control of air handling units using deep reinforcement learning and recurrent neural network,” _Building and Environment_ , vol. 168, p. 106535, 2020.
  * [6] S. Touzani, A. K. Prakash, Z. Wang, S. Agarwal, M. Pritoni, M. Kiran, R. Brown, and J. Granderson, “Controlling distributed energy resources via deep reinforcement learning for load flexibility and energy efficiency,” _Applied Energy_ , vol. 304, p. 117733, 2021.
  * [7] S. Brandi, M. Fiorentini, and A. Capozzoli, “Comparison of online and offline deep reinforcement learning with model predictive control for thermal energy management,” _Automation in Construction_ , vol. 135, p. 104128, 2022.
  * [8] B. Chen, P. Donti, K. Baker, J. Z. Kolter, and M. Berges, “Enforcing Policy Feasibility Constraints through Differentiable Projection for Energy Optimization,” _arXiv preprint arXiv:2105.08881_ , 2021.
  * [9] W. Valladares, M. Galindo, J. Gutierrez, W.-C. Wu, K.-K. Liao, J.-C. Liao, K.-C. Lu, and C.-C. Wang, “Energy optimization associated with thermal comfort and indoor air control via a deep reinforcement learning algorithm,” _Building and Environment_ , vol. 155, pp. 105–117, 2019.
  * [10] A. Nagy, H. Kazmi, F. Cheaib, and J. Driesen, “Deep reinforcement learning for optimal control of space heating,” _arXiv preprint arXiv:1805.03777_ , 2018.
  * [11] D. B. Crawley, L. K. Lawrie, F. C. Winkelmann, W. F. Buhl, Y. J. Huang, C. O. Pedersen, R. K. Strand, R. J. Liesen, D. E. Fisher, M. J. Witte _et al._ , “EnergyPlus: creating a new-generation building energy simulation program,” _Energy and buildings_ , vol. 33, no. 4, pp. 319–331, https://doi.org/10.1016/S0378--7788(00)00\,114--6, 2001.
  * [12] Z. Zhang, A. Chong, Y. Pan, C. Zhang, and K. P. Lam, “Whole building energy model for HVAC optimal control: A practical framework based on deep reinforcement learning,” _Energy and Buildings_ , vol. 199, pp. 472–490, 2019.
  * [13] C. Huang, H. Zhang, L. Wang, X. Luo, and Y. Song, “Mixed Deep Reinforcement Learning Considering Discrete-Continuous Hybrid Action Space for Smart Home Energy Management,” _Journal of Modern Power Systems and Clean Energy_ , 2022.
  * [14] L. Di Natale, B. Svetozarevic, P. Heer, and C. N. Jones, “Physically Consistent Neural Networks for building thermal modeling: theory and analysis,” _arXiv preprint arXiv:2112.03212_ , 2021.
  * [15] H. Van Hasselt, A. Guez, and D. Silver, “Deep reinforcement learning with double Q-learning,” in _Proceedings of the AAAI conference on artificial intelligence_ , vol. 30, no. 1, 2016.
  * [16] Empa, “NEST,” 2021, accessed: 04.10.2021. [Online]. Available: https://www.empa.ch/web/nest/overview
  * [17] S. Fujimoto, H. Hoof, and D. Meger, “Addressing function approximation error in actor-critic methods,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 1587–1596.
  * [18] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning,” _arXiv preprint arXiv:1509.02971_ , 2015.
  * [19] H. Hasselt, “Double Q-learning,” _Advances in neural information processing systems_ , vol. 23, pp. 2613–2621, 2010.
  * [20] L. Liu, H. Jiang, P. He, W. Chen, X. Liu, J. Gao, and J. Han, “On the Variance of the Adaptive Learning Rate and Beyond,” in _Proceedings of the Eighth International Conference on Learning Representations (ICLR 2020)_ , April 2020.
  * [21] P. Henderson, R. Islam, P. Bachman, J. Pineau, D. Precup, and D. Meger, “Deep reinforcement learning that matters,” in _Proceedings of the AAAI conference on artificial intelligence_ , vol. 32, no. 1, 2018.


# SaRNet: A Dataset for Deep Learning Assisted Search and Rescue with
Satellite Imagery

Michael Thoreau  Department of Electrical and Computer Engineering
New York University
<EMAIL_ADDRESS>Frazer Wilson  No Affiliation
<EMAIL_ADDRESS>
###### Abstract

Access to high resolution satellite imagery has dramatically increased in
recent years as several new constellations have entered service. High revisit
frequencies as well as improved resolution has widened the use cases of
satellite imagery to areas such as humanitarian relief and even Search and
Rescue (SaR). We propose a novel remote sensing object detection dataset for
deep learning assisted SaR. This dataset contains only small objects that have
been identified as potential targets as part of a live SaR response. We
evaluate the application of popular object detection models to this dataset as
a baseline to inform further research. We also propose a novel object
detection metric, specifically designed to be used in a deep learning assisted
SaR setting.

###### Index Terms:

Satellite, Remote, Search, Object Detection

##  I Introduction

In-domain datasets are currently indispensible for applying machine learning
models to real-world problems. One such problem is the search for missing
persons in remote locations, where access can be difficult and timing is
critical. So far, datasets for visual search and rescue (SaR) have mostly
contained images taken by UAVs or light aircraft. This data cannot be simply
translated to a satellite imagery setting, where there are currently no SaR
datasets that we know of, due to greater diversity in viewpoints and
relatively large target sizes. Modern high-resolution satellite constellations
that can be tasked with imaging almost anywhere on the planet in a small
number of hours, might soon enable a powerful complement to aerial searches,
particularly in combination with recent advances in deep learning.

We propose a novel object detection dataset, collected in a live search
setting, that we use to demonstrate the concept of deep learning assisted SaR.
The dataset was created during the search for a missing paraglider pilot, lost
in a remote and mountainous area of the western United States. Over 500
volunteers labelled potential targets in high-resolution images using axis-
aligned bounding boxes. The true target, as seen inset in figure 1, was found
after a three week search and the labels generated for potential targets were
saved. These images and annotations were post-processed to form this dataset
of 2552 images.

Figure 1: Inset - The paraglider wing as it was found. Main - The wing as
detected by our prototype system.

Search and rescue via satellite imagery is a challenging application for off-
the-shelf deep learning methods. Metrics typically used to evaluate a models
performance on datasets such as MS-COCO [6] are very informative when the
ground truth is fairly irrefutable and labelling is consistent, but have some
issues when labels are noisy. Systems that use human verification as part of
target acquisition can also generally tolerate a lower precision. We propose a
new metric that is better suited to deep learning assisted SaR, that provides
an intuitive way to choose a detection threshold for a given set of test
images. We will evaluate a number of popular object detection models using
this new metric on our dataset. Our contributions are as follows:

  * •

We present a novel dataset for satellite imagery based SaR, as well as;

  * •

A novel and specially informative metric for object detection in a SaR
setting, and;

  * •

We perform an comparative study of popular object detection models trained and
tested on this dataset.

Figure 2: Representative training samples, selected at random.

##  II Dataset Details

This dataset contains 2552 images with a total of 4206 axis aligned bounding
boxes of a single ‘target’ class. Volunteers were instructed to label anything
they think could be the missing paraglider wing and were provided with
examples of similar objects visible in the source data. The wing is shown
inset in figure 1 as it was found after a three week search. A total of
approximately 5000 annotations were originally generated, however bounding
boxes with a height or width greater than 20 meters were discarded, along with
the corresponding images. A mosaic of 20 labelled targets, selected at random,
can be seen in figure 2. Each 1000x1000 image in the dataset is a jpeg
rendered at a high quality factor corresponding to a 500m x 500m tile of
satellite data, with a pixel pitch of 0.5m. The images have had all geographic
and vendor information removed. The scale of the targets can be as small as
3-4 pixels in some cases, which is a particular challenge for machine learning
models, as we will discuss in a later section. The bounding box labels in this
dataset are generally slightly oversized relative to the target outlines. The
distribution of the length of the longest side of the boxes in the training
set can be seen in figure 3.

Figure 3: Distribution of the longest side of bounding boxes in the training
set.

In contrast to most object detection datasets, including those in the remote
sensing domain, the targets in this dataset cannot be considered a strong
ground truth, however we will argue that they can still be used to train a
useful object detector. Objects that might have been considered by labelers to
be a potential target may have varied over time and between volunteers. This
noise in the labels can be seen qualitatively in figure 2, where there is some
variance in the obviousness of the targets. For a dataset with a weaker ground
truth, metrics based directly on false positives and false negatives are not
as informative, as we will discuss in section V.

Examples have been split in to 70%, 20%, 10% divisions for training,
validation, and testing respectively, with annotations following the MS-COCO
[6] convention. This dataset will be available online with example code:
https://github.com/michaelthoreau/SearchAndRescueNet.

##  III Deep Learning Assisted SaR

It is important to consider how the outputs of any model might be used in a
wider system. Let us consider the case of the search [1] that created this
dataset, where targets were identified in a two step process. First,
volunteers labeled potential targets in image tiles without contextual
information. These labels were then verified by other volunteers with
additional information such as historical satellite imagery for comparison. In
this environment, a high number of false positives from the first step could
be tolerated, as they would be filtered in the second step, with only a small
number of potential targets being forwarded on to ground/air search teams.
Anecdotally, we found that verification of the potential targets (second
stage) was less tiring for volunteers than actively searching for targets in
image tiles (first stage).

We propose that an object detection model could be used to replace or assist
humans in the first stage of this SaR pipeline. In this proposed system,
humans would initially label data in a new target domain, and a model would be
trained to provide detections as search areas expanded. These detections could
be used either as proposals in the first stage, highlighting potential targets
and reducing strain on volunteers, or directly as potential targets for
verification in the second stage. In both cases, search areas could be covered
faster while keeping a human in the loop and reducing strain on resources.

##  IV Prior Works

A primary challenge for designing, evaluating, and deploying models to assist
in SaR missions is the lack of applicable datasets. One popular dataset [12]
for remote sensing object detection has bounding box annotations for a variety
of object categories, with the smallest being ‘small vehicle’. However this
dataset contains aerial imagery with a finer pixel pitch than the proposed
imagery. Fewer datasets exist for satellite based remote sensing data, and
none currently available contain objects as small as those in the proposed
dataset.

Detection of objects in remote imagery is a fairly well studied field, with
state of the art methods [8] using region based Convolutional Neural Networks
(CNNs) to detect objects such as aircraft and oil tanks with a high degree of
accuracy. These objects however have a scale in the order of 100s of pixels
across, making the features learnt on common pre-training tasks such as
ImageNet [3] very applicable and making transfer learning possible. One of the
closest approaches [10] to our problem achieves robust detection of small
objects with highly variable backgrounds, however this domain has a much more
variable object scale as well as significant viewpoint changes.

The use of deep learning to assist in SaR operations has been explored a
number of times [13, 2, 4] but not as far as we can tell with satellite
imagery. The authors of [13] discuss how the human eye has incredible power to
use context to discern true from false targets but is slow to scan images and
can quickly become fatigued. This group [13] also describes how small targets
can be detected faster than by the human eye in some circumstances. We propose
that an object detector can assist in SaR as a first step in a machine-human
process, with humans verifying potential targets.

##  V Object Detection Metrics

A key metric used in object detection is mean Average Precision (mAP) which is
calculated as the area under the curve when precision is plotted against
recall for all classes. Informed by the precision vs recall curve,
practitioners applying object detection models choose a threshold that
corresponds to a point on the precision recall curve that they deem most
reasonable. Due to the nature of our dataset that does not have a strong
ground truth, precision is not directly informative of the performance of the
model on the task. In the appendix, figure 5 shows some ‘false positives’ that
degrade the apparent performance of the model in standard metrics but appear
to be reasonable detections. We also found that when applying object detection
in an SaR setting, we were setting the detection threshold based on the
perceived density of detections. We propose a metric that can be directly
related to the time cost of verifying candidate detections. We will define the
detection density as the number of detections per square km. For a set of
detections with confidence values $\mathcal{Y}=\\{y_{1},y_{2},...,y_{M}\\}$
estimated for a given image $x_{i}$ in a set
$\mathcal{X}=\\{x_{1},x_{2},...,x_{N}\\}$ the detection density $D(\tau)$ can
be found as:

$D(\tau)=\sum_{i=1}^{N}\frac{|y_{i}>\tau|}{\textit{Area(}x_{i}\textit{)}}$ (1)

Where $|y_{i}>\tau|$ denotes the number of detections corresponding to a given
image with a confidence greater than $\tau$.

Similar to mAP, we propose a metric that considers the average performance
over a range of thresholds. In this case we will consider the recall of the
detector at various detection densities between 0 and some reasonable number
per square km, which can be chosen based on the dataset. We found that 20
detections per square km was a reasonable maximum density in this case. Recall
was included in the metric as we found it remained fairly relevant in the
presence of label noise, likely due to the labels being noisy but generally
conservative. We call this metric the Average Recall-Density to 20 or AR-d20.
The metric can be calculated as:

$\text{AR-d20}=\frac{1}{20}\sum_{k=0}^{20}\textit{Recall(}\tau_{k}\textit{)}$
(2)

Where $\tau_{k}$ is the lowest threshold that corresponds to a density lower
than $k$ detections per square km and $\textit{Recall(}\tau_{k}\textit{)}$ is
the recall for a given threshold. Crucially, this metric may inform the use of
a particular model based on the human resources available for considering
detections. In practice, the recall and detection densities can be estimated
on a small amount of labeled test data in the target domain.

##  VI Baseline Results

To acquire a baseline solution, we fine tune a number of popular object
detection models on the training data and evaluate the performance of each
model on the test set using our new metric. The framework used to evaluate
these models is Faster R-CNN[9] in which we consider three different feature
extractors available in the Detectron2 [11] ‘model zoo’. As the targets are
very small and the bounding box labels are relatively oversized, we relax the
Intersection over Union (IoU) threshold for detections to $0.1$ to avoid the
elimination of fair detections. We train all three models with a learning rate
of 0.0001 for 5000 iterations. The optimiser was Adam, and each minibatch
contained 4 images. All models were pre-trained on MS-COCO. We evaluate the
baseline solution on the AR-d20 metric and present the results in table I.

Figure 4: Density of detections and recall of 3 different Faster R-CNN models
fine-tuned on this dataset.

The highest scoring model ‘R_50_FPN’ includes a 50 layer resnet [5] feature
extractor. This model, which includes a Feature Pyramid Network [7], had an
average recall of 41.8% for detection densities between 0 and 20 detections
per square km. Figure 4 shows detection density and recall plotted for all
three models at a range of confidence thresholds. The FPN model has a
significantly higher recall than the other models but a similar detection
density across most thresholds, making it a compelling choice. Qualitative
detection results from applying this model to the validation set can be seen
in figure 6 in the appendix.

model | ARd-20
---|---
faster_rcnn_R_50_FPN_3x | 41.82
faster_rcnn_R_50_C4_3x | 28.88
faster_rcnn_R_50_DC5_3x | 35.46
TABLE I: Average Recall-density scores for three popular models.

##  VII Conclusion

We presented a novel dataset that we hope will inform future research into
satellite imagery based SaR. We introduced a new metric that may assist in the
application of object detectors to SaR problems and presented a baseline model
trained on this dataset to demonstrate the deep learning assisted SaR concept.
We believe that satellite based SaR is an emerging field that may someday be
used to save lives and bring closure to families of missing persons. We firmly
believe in the applicability of this technology to the search for missing
aircraft, watercraft, and a variety of other targets we cannot yet imagine.

## Acknowledgment

The authors would like to thank Planet Labs, Airbus Defence and Space, and
Maxar Technologies for providing satellite imagery. The authors would also
like to thank the 500+ volunteers who labelled the data, as well as all those
who contributed to the physical air and land search.

##  Appendix A Supplementary figures

Figure 5: A selection of ‘False Positives’ on the test set. These are objects
that were detected by the model but missed by volunteers.

## References

  * [1] L. Binding.  Body of well-known paraglider found in us mountains - weeks after he went missing.  Sky News, Sep 2020.
  * [2] G. Castellano, C. Castiello, C. Mencar, and G. Vessio.  Preliminary evaluation of tinyyolo on a new dataset for search-and-rescue with drones.  In 2020 7th International Conference on Soft Computing Machine Intelligence (ISCMI), pages 163–166, 2020.
  * [3] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei.  Imagenet: A large-scale hierarchical image database.  In 2009 IEEE conference on computer vision and pattern recognition, pages 248–255. Ieee, 2009.
  * [4] S. Gotovac, D. Zelenika, Z. Marusic, and D. Božić-Štulić.  Visual-based person detection for search-and-rescue with uas: Humans vs. machine learning algorithm.  Remote Sensing, 12, 10 2020.
  * [5] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition, 2015.
  * [6] T. Lin, M. Maire, S. J. Belongie, L. D. Bourdev, R. B. Girshick, J. Hays, P. Perona, D. Ramanan, P. Dollár, and C. L. Zitnick.  Microsoft COCO: common objects in context.  CoRR, abs/1405.0312, 2014.
  * [7] T.-Y. Lin, P. Dollár, R. Girshick, K. He, B. Hariharan, and S. Belongie.  Feature pyramid networks for object detection, 2017.
  * [8] Y. Long, Y. Gong, Z. Xiao, and Q. Liu.  Accurate object localization in remote sensing images based on convolutional neural networks.  IEEE Transactions on Geoscience and Remote Sensing, 55(5):2486–2498, 2017.
  * [9] S. Ren, K. He, R. B. Girshick, and J. Sun.  Faster R-CNN: towards real-time object detection with region proposal networks.  CoRR, abs/1506.01497, 2015.
  * [10] M. Schembri and D. Seychell.  Small object detection in highly variable backgrounds.  In 2019 11th International Symposium on Image and Signal Processing and Analysis (ISPA), pages 32–37, 2019.
  * [11] Y. Wu, A. Kirillov, F. Massa, W.-Y. Lo, and R. Girshick.  Detectron2.  https://github.com/facebookresearch/detectron2, 2019.
  * [12] G. Xia, X. Bai, J. Ding, Z. Zhu, S. J. Belongie, J. Luo, M. Datcu, M. Pelillo, and L. Zhang.  DOTA: A large-scale dataset for object detection in aerial images.  CoRR, abs/1711.10398, 2017.
  * [13] K. Yun, L. Nguyen, T. Nguyen, D. Kim, S. Eldin, A. Huyen, T. Lu, and E. Chow.  Small target detection for search and rescue operations using distributed deep learning and synthetic data generation.  CoRR, abs/1904.11619, 2019.

Figure 6: Example detections from the best performing model on the validation
set. Each tile in the mosaic is a fixed size and is centred on the detection.
Some examples are shown multiple times due to overlapping bounding boxes of
very different sizes.