XX Month, XXXX XX Month, XXXX XX Month, XXXX XX Month, XXXX XX Month, XXXX
OJIM.2022.1234567

Corresponding Author: M. Muzeau (e-mail: max.muzeau@ onera.fr).

# Deep Learning-Based Anomaly Detection in Synthetic Aperture Radar Imaging

Max Muzeau1,2  Chengfang Ren2  Member  IEEE  Sébastien Angelliaume2  Mihai
Datcu3  Fellow  IEEE  and Jean-Philippe Ovarlez1,2  Member  IEEE  SONDRA,
CentraleSupélec, Université Paris-Saclay, France  DEMR, ONERA, Université
Paris-Saclay, France  University Politehnica of Bucharest (UPB), Romania and
German Aerospace Center (DLR), Germany

###### Abstract

In this paper, we proposed to investigate unsupervised anomaly detection in
Synthetic Aperture Radar (SAR) images. Our approach considers anomalies as
abnormal patterns that deviate from their surroundings but without any prior
knowledge of their characteristics. In the literature, most model-based
algorithms face three main issues. First, the speckle noise corrupts the image
and potentially leads to numerous false detections. Second, statistical
approaches may exhibit deficiencies in modeling spatial correlation in SAR
images. Finally, neural networks based on supervised learning approaches are
not recommended due to the lack of annotated SAR data, notably for the class
of abnormal patterns. Our proposed method aims to address these issues through
a self-supervised algorithm. The speckle is first removed through the deep
learning SAR2SAR algorithm. Then, an adversarial autoencoder is trained to
reconstruct an anomaly-free SAR image. Finally, a change detection processing
step is applied between the input and the output to detect anomalies.
Experiments are performed to show the advantages of our method compared to the
conventional Reed-Xiaoli algorithm, highlighting the importance of an
efficient despeckling pre-processing step.

adversarial autoencoder, anomaly detection, deep-learning, despeckling, SAR,
self-supervised

##  1 INTRODUCTION

Anomaly detection is one of the most critical issues in multidimensional
imaging, mainly in hyperspectral [1, 2] and medical imaging . Even if we have
no prior information about a target or background signature, anomalies
generally differ from surrounding pixels due to their dissimilar signatures.
Anomaly detection for radar and SAR imaging aims to discover abnormal patterns
hidden in multidimensional radar signals and images. Such anomalies could be
man-made changes in a specific location or natural processes affecting a
particular area. These anomalies can characterize, for example, several
potential applications: Oil slick detection [3], turbulent ship wake [4],
levee anomaly [5] or archaeology [6]. They can also be related to any change
detection in time-series SAR images. This research field is essential in data
mining for quickly isolating irregular or suspicious segments in large amounts
of the database. Many anomaly detection schemes have been proposed in the
literature [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]. Among them, the
unsupervised methods are the most interesting since they are widely applicable
and do not require labeling the data.

Deep learning techniques for anomaly detection are often based on an encoding-
decoding network architecture to learn healthy data or images that do not
contain anomalies, such as references [14, 15]. Another approach is based on
the one-class SVM method in the latent space [17]. Anomaly detection in SAR
using the difference between the input and the reconstructed image obtained
through the autoencoder scheme is also proposed in [16]. Still, it could
suffer from the SAR speckle noise.

Our approach consists in mitigating this annoying speckle before training an
Adversarial Autoencoder (AAE) from these despeckled data. It does not need to
separate abnormal SAR images from the training dataset.We then compare our
method with the conventional RX anomaly detector [7] by analyzing their global
performance, and we highlight the importance of the despeckling process.

This paper is organized as follows. After some introduction in Section I,
Section II introduces the context of SAR imaging. The proposed methodology
being based on despeckling of SAR images, Section III describes the filtering
process. Section IV presents the self-supervised neural architecture that aims
to reconstruct the final SAR image without anomalies. Section V will define
the last phase of the processing based on the change detection step and will
discuss the quality of the anomaly detection procedure. The proposed
methodology is applied to experimental polarimetric SAR images. Final Section
VI gives some conclusions and perspectives.

Figure 1: Despeckling method with the architecture described in [22]. First,
we decompose a polarimetric SAR image in four different intensity images
(discarding the phase). Each polarization is then denoised separately, and the
polarimetric SAR image is finally reconstructed.

##  2 SYNTHETIC APERTURE RADAR

Airborne and spaceborne SAR aims to provide images of Earth’s surface at radar
frequency bands. Contrary to optical imaging, they can work day and night
using an active radar that transmits and receive pulses in the scanning area.
They offer an opportunity to monitor changes and anomalies on the Earth’s
surface. Their principle is to combine multiple received signals coherently to
simulate an antenna with a larger aperture. This procedure allows building a
complex-valued image of the terrain with very high range and azimuth
resolutions, independently of the distance between the radar and the imaged
area. Electromagnetic waves can also be polarized at emission and reception.
Horizontal and vertical polarizations are generally used, resulting in four
coherent SAR channels ($HH,HV,VH,VV$). $HV$ corresponds to a horizontal
polarization for the emission and vertical polarization for the reception. The
information provided by polarimetry is crucial for most geoscience
applications from SAR remote sensing data [25].

###  2.1 Speckle

The main issue in SAR imaging concerns the corruption of the backscattered
signal within a resolution cell by a multiplicative noise called speckle. This
is due to the coherent summation of multiple scatterers in one resolution
cell, which may cause destructive or constructive inferences. This phenomenon
disturbs the exploitation of radar data for detection and geoscience
applications. Much work in the literature has proposed reducing its effect
(e.g. [18, 19, 20, 21, 22, 23, 24]) while preserving the resolution as much as
possible.

A model of the speckle has been defined by Goodman [26]. In the case of a
single look complex (SLC) image, each pixel power or intensity $I$ may be
distributed according to an exponential distribution:

$p(I|R)=\frac{1}{R}\,\exp{\left(-\frac{I}{R}\right)}\,,$ (1)

where $R$ denotes the mean reflectivity level. A useful general statistical
model representation is the so-called multiplicative noise model:

$I=R\,S\,,$ (2)

where $S$ characterizes the exponentially distributed speckle and $R$ is a
scalar positive random variable characterizing the texture.

To change the multiplicative noise model into an additive one, a log
transformation is often used, which leads to a new distribution:

$p(y|x)=\exp{(y-x)}\,\exp(-\exp{(y-x)})\,,$ (3)

where $y=\log(I)$ and $x=\log(R)$. One crucial detail is the spatial
correlation between the pixels which is not always taken into account in the
above model. This can be mainly explained by a possible apodization of the
images or some oversampling lower than the spatial resolution during the image
synthesis. We can directly work on full-resolution data thanks to recent
developments in deep learning despeckling algorithms that know how to keep the
spatial resolutions and preserve details such as lines and small bright
targets (boats or vehicles, for example). They, therefore, render detection
tasks and false alarm regulation easier on speckle-free data.

###  2.2 Dataset

In this paper, the analyzed dataset is composed of a time series of X-band SAR
images (each of size $4800\times 30000\times 4$) acquired by SETHI, the
airborne instrument developed by ONERA [27, 28]. The resolutions of these
images are about $20$ cm in both azimuth and range domains for the four
polarization channels ($HH,HV,VH,VV$).

In the monostatic case, the channels $VH$ and $HV$ are often averaged because
they contain the same information. The averaging decreases the speckle impact
without degrading the resolution (so-called reciprocity principle [29]). The
resulting three channels $[HH,\frac{1}{2}(HV+VH),VV]$ are then thresholded
separately using a value $\lambda(I)$ chosen as:

$\lambda(I)=\mu_{I}+3\,\sigma_{I}\,,$ (4)

where $\mu_{I}$ and $\sigma_{I}$ are respectively the estimated mean and the
standard deviation of the image $I$. The threshold $\lambda(I)$ is the same
for all SAR images to have a consistent visualization, as shown in Fig. 1.

##  3 Despeckling

To reduce the speckle impact, the deep-learning despeckling algorithm SAR2SAR
[22] has been used. As discussed previously, this method can estimate the
spatial correlation between pixels and thus works with full-resolution images.

The training phase is based on a method called Noise2Noise first developed in
[30]. It shows that there is no need for a ground truth to train a denoising
deep neural network, but only for multiple images of the same area with
different noise realizations. A common strategy for estimating the true
unknown image is to find, for the quadratic $L^{2}$ loss function, an image
$\bf{z}$ that has the smallest average deviation from the measurements:

$\operatorname*{argmin}_{\mathbf{z}}\mathbb{E}_{\mathbf{y}}\big{[}\left\lVert\mathbf{z}-\mathbf{y}\right\rVert^{2}\big{]}=\mathbb{E}_{\mathbf{y}}[\mathbf{y}]\,,$
(5)

where $\mathbf{y}$ represents a set of observations. If the sample noise is
additive, centered, and independent in all the observations, estimating
$\mathbb{E}_{\mathbf{y}}[\mathbf{y}]$ is the same as predicting noise-free
observations.

The same principle has been applied and refined for SAR images in the
algorithm [22]: let $\mathbf{y}_{1}$ and $\mathbf{y}_{2}$ be two independent
realizations of identically distributed random variables and let a denoising
network $f_{\mathbf{\theta}}(.)$ be parameterized by $\mathbf{\theta}$. A loss
$L_{speck}$ based on the negative log-likelihood of log-intensity images can
be used according to the distribution defined in Eq. (3) and leads to:

$L_{speck}=\sum_{k}f_{\mathbf{\theta}}([\mathbf{y}_{1}]_{k})-[\mathbf{y}_{2}]_{k}+\exp\left([\mathbf{y}_{2}]_{k}-f_{\mathbf{\theta}}([\mathbf{y}_{1}]_{k}\right)$
(6)

The loss is a sum of all partial loss for each pixel represented with the
index $k$. The use of Eq. (6), instead of a mean square error, allows a faster
convergence.

The other difference between SAR2SAR and Noise2Noise is the training phase
which is decomposed into three steps:

  * •

Phase A: The training step is done with a speckle-free image $\mathbf{y}$. Two
fake noisy image are created such that
$\mathbf{y}_{1}=\mathbf{y}+\mathbf{s}_{1}$ and
$\mathbf{y}_{2}=\mathbf{y}+\mathbf{s}_{2}$ where $\mathbf{s}_{1}$ and
$\mathbf{s}_{2}$ are both distributed according the distribution from Eq. (3).
The network is learning how to denoise data without any knowledge of spatial
correlation.

  * •

Phase B: The training step is performed with experimental SAR images and a
change compensation pre-processing: two same areas acquired at a different
time $\mathbf{y}_{t_{n}}$ and $\mathbf{y}_{t_{m}}$ are used. To compensate for
changes that could have occurred between the acquisition at two dates, an
estimation of the reflectivity is done beforehand (on sub-sampled images to
remove spatial correlation). This leads to $\mathbf{y}_{1}=\mathbf{y}_{t_{n}}$
and
$\mathbf{y}_{2}=\mathbf{y}_{t_{m}}-f_{\mathbf{\theta}}(\mathbf{y}_{t_{m}}^{sub})+f_{\mathbf{\theta}}(\mathbf{y}_{t_{n}}^{sub})$
where the $sub$ exponent characterizes sub-sampled data. This finally allows
$\mathbf{y}_{2}$ to have the speckle of $\mathbf{y}_{t_{m}}$ and the
reflectivity of $\mathbf{y}_{t_{n}}$.

  * •

Phase C: The same operation is done, but this time, the network has learned
how to model spatial correlation, so the images are not sub-sampled to
estimate a reflectivity beforehand. This gives
$\mathbf{y}_{1}=\mathbf{y}_{t_{n}}$ and
$\mathbf{y}_{2}=\mathbf{y}_{t_{m}}-f_{\mathbf{\theta}}(\mathbf{y}_{t_{m}})+f_{\mathbf{\theta}}(\mathbf{y}_{t_{n}})$.

##  4 Image reconstruction

Once the pre-processing step has been performed, the goal is to highlight
anomalies. This is obtained in an unsupervised manner, and the network has no
information about what should be an anomaly and what shouldn’t. This makes the
task harder than, for example, with a supervised convolution network. But, if
we can reach the state-of-the-art with an unsupervised algorithm, one of the
most significant drawbacks of artificial intelligence, corresponding to the
need for a high-quality labeled dataset will be overcome.

###  4.1 Proposed architecture

We use an AAE [31] with convolution layers to deal with this problem. To
understand the architecture, it is first necessary to introduce what is a
Generative Adversarial Network (GAN) [32]. We will use the notations
$\mathcal{E}(.)$, $\mathcal{D}(.)$ and $\mathcal{D}_{c}(.)$ to define the
encoder, the decoder, and the discriminator respectively. A description of the
network is illustrated in Fig. 4.

####  4.1.1 Generative adversarial networks

A GAN is a self-supervised deep learning algorithm that was originally used to
generate synthesized images based on what it saw. A generator and a
discriminator are trained jointly against each other. The generator’s goal is
to create the most realistic image possible based on a source vector of low
dimension. Usually, this vector is distributed according to a multivariate
Normal distribution. The generated image should fool the discriminator into
thinking there is no difference between real and fake images. In opposition,
the discriminator’s goal is to be able to distinguish fake generated data from
real ones. In our architecture, the role of the GAN is to generate a latent
vector that is distributed according to a Normal distribution representative
of a patch of polarimetric SAR speckle-free image.

####  4.1.2 Adversarial Autoencoder

An AAE is composed of two things:

  * •

A generator, which is an encoder and a decoder, one after another,

  * •

A discriminator that is placed between the encoder and the decoder. His goal
is to ensure that the latent space follows a normal distribution.

Our input is the logarithm of a denoised polarimetric SAR patch
$\mathbf{X}\in\mathbb{R}^{h\times w\times c}$ of size $h\times w$ and depth
$c$, where $h$ and $w$ are respectively the numbers of pixels in azimuth and
range, and where $c$ is the polarization channel. The logarithm operation is
used to reduce the dynamic of the data (for example, there can be a difference
of $10^{3}$ between the amplitude of a strong scatterer and the amplitude of a
pixel located in a vegetation area). After passing through the encoder, we get
a latent vector $\mathbf{z}=\mathcal{E}(\mathbf{X})$ with $\mathbf{z}\sim
p(\mathbf{z})$ and where $p(.)$ is the a priori distribution of our encoder.
Based on this vector, the decoder will then make an estimation of our input
patch $\hat{\mathbf{X}}=\mathcal{D}(\mathbf{z})$. In our architecture, we can
observe a GAN present with the generator $\mathcal{E}(.)$ and the
discriminator $\mathcal{D}_{c}(.)$ whose purpose is to differentiate
$\mathbf{z}_{real}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$ from
$\mathbf{z}_{fake}\sim p(\mathbf{z})$. Here, the latent space is distributed
according to a reduced and centered Normal distribution, but the Uniform
distribution on [0,1] could also have given similar results. One restriction
is that it could not have been a long tail distribution because the
information would not have been compacted in a small enough space. We then
train the AAE in two successive phases:

  * 1

Reconstruction error: this characterizes a loss that ensures to have a low
pixel-per-pixel error. We minimize $L_{rec}$, which is a $L^{1}$ norm. It has
been preferred instead of a $L^{2}$ norm because we do not want to heavily
penalize the network when the difference between $\mathbf{X}$ and
$\hat{\mathbf{X}}$ is huge.

$L_{rec}=\frac{1}{h\,w\,c}\displaystyle\sum\limits_{i,j,k}\big{\lVert}\mathbf{X}_{i,j,k}-\hat{\mathbf{X}}_{i,j,k}\big{\rVert}_{1}\,.$
(7)

  * 2

Regularization error: such a loss allows us to control the distribution
$p(\mathbf{z})$. It also gives a better reconstruction according to [17, 31].
To do so, the weights of the encoder $\mathcal{E}$ and the discriminator
$\mathcal{D}_{c}$ are estimated according to:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\underset{\mathcal{E}}{\min}\,\underset{\mathcal{D}_{c}}{\max}\,L_{lat}\\!\\!\\!$
$\displaystyle\\!\\!\widehat{=}\\!\\!$
$\displaystyle\\!\\!\mathbb{E}_{\mathbf{z}_{real}\sim\mathcal{N}(\mathbf{0},\mathbf{I})}\left[\log(\mathcal{D}_{c}(\mathbf{z}_{real})\right]$
(8) $\displaystyle+$
$\displaystyle\\!\\!\mathbb{E}_{\mathcal{E}(\mathbf{X})\sim
p(z)}\left[\log\left(1-\mathcal{D}_{c}(\mathcal{E}(\mathbf{X})\right)\right]\,.$

####  4.1.3 Link between anomaly detection and AAE

It is not apparent to see the link between this architecture and anomaly
detection. The goal of an AAE for this application will be to make an accurate
estimation for Normal distributed data and a bad estimation for abnormal data.
This will be based on the assumption that the network is not powerful enough
to reconstruct every aspect of the image. Only the recurrent patterns will be
remembered. By definition, a spatial anomaly occurs rarely compared to the
rest of the data. If all these assumptions are valid, only rare patterns will
be seen as anomalies, which is exactly what we want.

Figure 2: Overview of our proposed architecture for SAR anomaly detection.
$X_{noisy}$ corresponds to the amplitude of SAR image, $X$ is the speckle-free
version, and $\hat{X}$ is its reconstruction using an adversarial autoencoder.
Finally, the anomaly map $A^{X}$ is obtained through a change detection
algorithm between $X$ and $\hat{X}$

##  5 Change detection method

Once the network has delivered an estimation of the input image, the goal is
to detect changes between $\mathbf{X}$ and $\hat{\mathbf{X}}$. This detection
will be represented in an anomaly map $A^{\mathbf{X}}\in[0,1]^{h\times w}$.
The closer the value is to one, the more the pixel is likely to be an anomaly.
The global flow of the detection procedure is shown in Fig. 2.

###  5.1 Problem formulation

There are multiple approaches to detect changes between two or more images
[33]. As described in [33], pixel-level comparisons are widely used in SAR
imaging community. The problem is formulated within the set of the following
two hypotheses:

$\left\\{\begin{array}[]{ll}H_{0}:\boldsymbol{\theta}_{1}=\boldsymbol{\theta}_{2}\text{
(no anomaly),}\\\
H_{1}:\boldsymbol{\theta}_{1}\neq\boldsymbol{\theta}_{2}\text{
(anomaly),}\end{array}\right.$ (9)

where $\boldsymbol{\theta}_{1}$ and $\boldsymbol{\theta}_{2}$ are vectors of
the estimated parameters of the distribution used to model pixel values.
Methods relying on the hypothesis test on the statistics of the image are
difficult to exploit. The non-Gaussianity, the heterogeneity of the SAR
images, or their complex-valued nature make this derivation very difficult.

When pixels have been transformed into log intensity for better contrast and
when the speckle effect has been reduced beforehand, the final distributions
of the proposed anomaly detection tests are generally unknown. Hence, to be
able to compare all the proposed strategies, a visualization process is used
to threshold each map according to the following equation:

$\tilde{x}=\min(x,t)\,.$ (10)

For a given percentage $p$, the threshold $t$ is fixed such that $p\%$ of the
pixels are above it. This ensures to have the same number of pixels of value
$t$ for each result. The dynamic is compressed in a way that allows us to see
at the same time the anomalies and the background. The Probability of False
Alarm (PFA) is here, for convenience, characterized by this value $p$ in the
sense that $p$ and PFA are equal if all the pixels in the map correspond to
$H_{0}$ hypothesis. They will have the same significance in the sequel.

One way to detect changes between two SAR images consists of testing the
equality between the two estimated covariance matrices of the corresponding
pixel under test for each pixel. This can be made statistically if the
knowledge of the data statistic is known [34] or through matrix distances
[35].

Figure 3: Comparison of anomaly maps clipped at the highest $2\%$ values:
(left) Input image, (middle) $L_{1}$ loss between $X$ and $\hat{X}$ and
(right) $A^{\mathbf{X}}_{E}$.

###  5.2 Covariance estimation

Covariance matrices are estimated locally around the pixel under test. This is
useful to strongly reduce the noise in the anomaly map compared to a standard
$L^{1}$ loss, as it is shown in Fig. 3. Indeed, to estimate a covariance
matrix, we use a sliding window represented by a boxcar $\mathcal{B}_{k,l}$
where $k,l$ is the coordinate of its center. For multivariate Gaussian
distribution $\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$, the Maximum
Likelihood Estimators of the mean vector $\boldsymbol{\mu}$ and the covariance
matrix $\boldsymbol{\Sigma}$ lead to the well known Sample Mean Vector (SMV)
and the Sample Covariance Matrix (SCM) which are defined as:

$\hat{\boldsymbol{\mu}}^{\mathbf{X}}_{k,l}=\displaystyle\frac{1}{|\mathcal{B}_{k,l}|}\sum\limits_{i,j\in\mathcal{B}_{k,l}}\mathbf{X}_{i,j}\,,$
(11)
$\hat{\boldsymbol{\Sigma}}^{\mathbf{X}}_{k,l}=\frac{1}{\left|\mathcal{B}_{k,l}\right|}\sum\limits_{i,j\in\mathcal{B}_{k,l}}\left(\mathbf{X}_{i,j}-\hat{\boldsymbol{\mu}}^{\mathbf{X}_{k,l}}\right)\left(\mathbf{X}_{i,j}-\hat{\boldsymbol{\mu}}^{\mathbf{X}_{k,l}}\right)^{T}\,,$
(12)

where $\hat{\boldsymbol{\mu}}^{\mathbf{X}_{k,l}}$ is the estimate of the mean
vector $\boldsymbol{\mu}$, where
$\hat{\boldsymbol{\Sigma}}_{k,l}^{\mathbf{X}}$ is the estimate of the
covariance matrix $\boldsymbol{\Sigma}$ associated to the boxcar
$\mathcal{B}_{k,l}$ in the image $\mathbf{X}$ and where $|.|$ denotes the
cardinal operator.

Figure 4: Architecture of the Adversarial Autoencoder. The numbers above
blocks represent the number of input channels for the encoder and the number
of output channels for the decoder and the discriminator (disc). The
description of each block is done with Pytorch notations.

###  5.3 Distance metric

There are multiple possibilities to compute a distance between two matrices
[35, 36]. A common way to do this consists in computing the square of the
Frobenius norm of the difference between the two matrices.

$A^{\mathbf{X}}_{E}(k,l)=\left\|\hat{\boldsymbol{\Sigma}}_{k,l}^{\mathbf{X}}-\hat{\boldsymbol{\Sigma}}_{k,l}^{\hat{\mathbf{X}}}\right\|_{F}^{2}\,.$
(13)

Other methods, based on the Frobenius norm, could also be applied on log-
matrices or root matrices. In our detection case, a lot of importance is given
to the intensity difference between $\hat{\boldsymbol{\Sigma}}^{\mathbf{X}}$
and $\hat{\boldsymbol{\Sigma}}^{\hat{\mathbf{X}}}$. The Euclidean metric (13)
highlights the difference in intensity between the two covariance matrices
characterizing each pixel and its reconstructed value while preserving a low
PFA. The pseudo-code for the proposed anomaly detection method is detailed in
Algorithm 1.

# SAR2SAR: despeckling network

# AAE: reconstruction network

# k: semi kernel size

# scm: compute the sample covariance matrix

# norm: Frobenius norm

for x_noisy in eval_data:

x = SAR2SAR(x_noisy)

x_rec = AAE(x)

(_,m,n) = x.shape # channels,height,width

ano_map = zeros(m,n)

# Compute anomaly score for each pixel

for i in range(m):

for j in range(n):

# Crop to compute scm

x_c = x[:,i-k:i+k+1,j-k:j+k+1]

x_r_c = x_rec[:,i-k:i+k+1,j-k:j+k+1]

c, c_r = scm(x_c),scm(x_r_c)

ano_map[i,j] = norm(c - c_r).pow(2)

m,M = min(ano_map), max(ano_map)

ano_map = (ano_map - m) / (M - m)

save(ano_map)

Algorithm 1 Anomaly detection pseudocode in python style

##  6 Experiences and analysis

In this section, we experiment with the proposed algorithm on the ONERA SETHI
SAR dataset. First, the analysis of the despeckling network is presented. We
then evaluate the change detection method and compare it with a standard
approach.

###  6.1 Despeckling quality

The proposed algorithm is applied to images that are decomposed into patches
of size $256\times 256$.

For phase A, we first need to compute a speckle-free image by averaging co-
registered images acquired at different dates. Every polarization has been
averaged, and we then used the algorithm MuLoG [21] to remove the remaining
speckle. It gives us training data of 3052 patches grouped in batches of size
32, which leads to a total of 955 batches. The network is then trained for 20
epochs.

Phases B and C use images from the dataset previously described, and the pre-
processing step is the one described in Phase B and C, Section 3. There are
four piles of two images, one for each polarization, and each image is
decomposed into 33047 patches with a batch size of 32, which leads to a total
of 1032 batches. The network is trained for ten epochs in phases B and C.

Figure 5 presents the result of the despeckling process for $HH$ polarization,
and the corresponding ratio image is defined as:

$ratio=I_{speckle}/I_{despeckling}\,,$ (14)

where $I_{speckle}$ and $I_{despeckling}$ represent respectively the intensity
of the original image and the speckle-free image.

This test image, relatively heterogeneous, comprises vegetation, roads, thin
lines, and a strong scatterer annotated as a vehicle. The image dynamic has
not been altered, and because there is no structure in the ratio image, we may
confidently affirm that our despeckling network has succeeded in only removing
the speckle.

Figure 5: Comparison between (top left): $HH$ SAR image, (top-right): its
speckle-free version. (bottom-left): the intensity ratio of these images.
(bottom-right): superposition of the histogram of the ratio image and the
theoretical Exponential distribution with parameter $\lambda=1$.

The speckle on its own should follow an exponential distribution with its
parameter equal to one. In the histogram in Fig. 5, we found that, as
expected, the actual distribution follows an exponential probability density
function.

###  6.2 Anomaly detection

To assess the quality of our results, they are compared with the so-called
Reed-Xiaoli detector [7] on complex-valued SAR images. The algorithm consists
in estimating locally, around a pixel characterized by its polarimetric
response $\bf{x}$, the associated mean vector $\hat{\boldsymbol{\mu}}$ and the
covariance matrix $\hat{\boldsymbol{\Sigma}}$ of its surrounding background
and then to test if this pixel under test is belonging, or not, to this
background. The anomaly score is computed through the well-known Mahalanobis
distance:

$RX(\bf{x})=\left(\bf{x}-\hat{\boldsymbol{\mu}}\right)^{H}\,\hat{\boldsymbol{\Sigma}}^{-1}\,\left(\bf{x}-\hat{\boldsymbol{\mu}}\right)\,,$
(15)

where $H$ denotes the transpose conjugate operator. The parameters are locally
estimated through the Gaussian Sample Mean Vector and Sample Covariance Matrix
estimates. An exclusion window prevents the use of anomalous data in parameter
estimation (guard cells).

Additional comparisons will be made using our proposed AAE algorithm without
the despeckling pre-processing and the $L^{1}$ norm between $\mathbf{X}$ and
$\hat{\mathbf{X}}$.

To run the proposed algorithm, we first need to train the AAE on the dataset
described in 2.2.2

####  6.2.1 AAE training and architecture

The training process is self-supervised. The training dataset has no labels
and potentially can contain anomalies. In this case, the evaluation and
training datasets can be mixed because the training process requires only
unlabeled data. In the same way as for the despeckling algorithm, a log
transformation is used, and the input data are normalized between 0 and 1.
These images are decomposed into patches of size $64\times 64$. The sliding
window has a stride of 16. This leads to a total of 163785 patches grouped in
batches of size 128, so we get 1279 batches for one training epoch with a
total of 20 epochs.

The architecture of the network is illustrated in Fig. 4. It has been designed
according to [37]. To update the weights, we use the optimizer Adam [38] with
a cyclical learning rate [39] that goes linearly from $10^{-3}$ to $10^{-2}$.
It takes 2558 batches to go from one value to another. A complete cycle will
take four epochs. This method helps to have a robust training phase that will
converge even if we do not know the perfect learning rate for the network.

####  6.2.2 Evaluation dataset

The proposed method has been qualitatively evaluated on a known area of the
dataset described in Section 2.2.2. Since there is no ground truth available,
the boundary between the anomaly and the background is not apparent, which
makes annotation almost impossible in many cases. The area is composed of
known anomalies (vehicles) intentionally placed by the ONERA team during a
measurement campaign.

(a)

(b)

(c)

Noisy

Denoised

Reconstruction

Difference

Figure 6: Evolution of the different steps of the process for three different
areas of the SAR image. (left): (a) SAR image containing a fence and vehicles.
(b): SAR image containing vehicles and trihedral sidelobes corresponding to
purple lines. (c): SAR image containing three vehicles.

HH

Intensity

Position

HV

Intensity

Position

Figure 7: Comparisons of horizontal profiles on a strong scatterer of the
image (a) in Figure 6 in $HH$ and $HV$ polarization.

####  6.2.3 Qualitative results

The evaluation of noisy images is displayed in Fig. 6 with their denoised and
reconstructed versions. The map named Difference represents the absolute value
between Denoised and Reconstruction maps. We can remark that the reconstructed
images are blurry, but all the essential structures are kept.

The despeckling pre-processing is not removing any structure, even the
smallest one, like the green line in the bottom right of the image (b). In all
the reconstructed images, the intensity of each point-like anomaly is greatly
reduced. Figure 7 illustrates this fact and shows a comparison of horizontal
profiles at the location characterizing a strong scatterer of image (a). Red
dotted lines correspond to the visualization threshold defined in Eq. (4).
This intensity attenuation is fundamental to the algorithm because the change
detection is only applied between the denoised and the reconstructed image in
order to highlight potential anomalies.

Figure 9 presents a comparison of anomaly detection results obtained with
different methods. All the results are here clipped at the highest $1\%$
values, see Eq 10.

The anomaly maps obtained with the methods $A^{\mathbf{X_{noisy}}}_{E}$, the
Reed-Xiaoli RX detector and $L^{1}(\mathbf{X},\hat{\mathbf{X}})$ have a
similar drawback characterizing the presence of a large number of false
detections. The proposed method $A^{\mathbf{X}}_{E}$ is shown to have better
performance for the same PFA.

input

label

RX

$A^{\mathbf{X_{noisy}}}_{E}$

$L^{1}(\mathbf{X},\hat{\mathbf{X}})$

$A^{\mathbf{X}}_{E}$

Figure 8: Comparison between different anomaly maps on a real SAR image with
fake anomalies and its corresponding label.

####  6.2.4 Quantitative results

For quantitative evaluation, we have embedded synthetic test patterns with
different intensity levels in a true anomaly-free crop SAR image.

The fake anomaly map results are represented in Fig. 8. The values are based
on true intensities characterizing the dataset with a decreasing significance
from top to bottom. In this figure, anomaly detection maps are clipped at
$10\%$ of the higher value because some anomalies have a value close to the
background. This makes the detection harder; thus, we reduce the dynamic range
to have a better representation. One advantage of autoencoder-based detection
is detecting abnormal areas of low intensities (left cross), contrary to the
RX detector. The proposed $A^{\mathbf{X}}_{E}$ method outperforms all the
other ones and even for the detection of low-intensity patterns. This is
illustrated in Fig. 10, which shows the overall performance obtained in terms
of Area Under the Curve (AUC). The proposed $A^{\mathbf{X}}_{E}$ still has the
better performance.

$A^{\mathbf{X}}_{E}$

RX

$L^{1}(\mathbf{X},\hat{\mathbf{X}})$

$A^{\mathbf{X_{noisy}}}_{E}$

Figure 9: Comparison between different anomaly maps obtained from the images
described in Fig. 6.

##  7 Conclusion

In this article, we propose a novel anomaly detection algorithm. It is mainly
based on the use of adversarial autoencoders. They have not been widely used
in the SAR community because of the speckle noise, which dramatically degrades
algorithm performance. Thanks to the recent advances in deep learning
despeckling algorithms specifically developed for SAR images, we can now
efficiently develop new algorithms based on these methods to enhance anomaly
detection performance.

Our proposed algorithm outperforms the conventional Reed-Xiaoli method since
it can detect abnormal areas with low-intensity values. Because this self-
supervised strategy does not require labeled data, it can easily be extended
to another type of data as long as the anomaly quantity remains negligible.

The perspectives envisaged concern the improvement of change detection
techniques that can help in giving better detection performance as well as in
regulating the False Alarm Rate.

00.10.20.30.40.50.60.70.80.9100.10.20.30.40.50.60.70.80.91False positive
rateTrue positive rate$A^{\mathbf{X}}_{E}$ \-
AUC:0.8774$L^{1}(\mathbf{X},\hat{\mathbf{X}})$ \-
AUC:0.8099$A^{\mathbf{X}_{noisy}}_{E}$ \- AUC:0.7751RX - AUC:0.7531 Figure 10:
ROC curve and value of the AUC for the anomaly detection algorithms displayed
Fig. 8.

## References

  * [1] D. W. J. Stein, S. G. Beaven, L. E. Hoff, E. M. Winter, A. P. Schaum, and A. D. Stocker, ”Anomaly detection from hyperspectral imagery,” in IEEE Signal Processing Magazine, 19(1):58-69, Jan. 2002.
  * [2] S. Matteoli, M. Diani and G. Corsini, ”A tutorial overview of anomaly detection in hyperspectral images,” in IEEE Aerospace and Electronic Systems Magazine, 25(7):5-28, July 2010.
  * [3] A. Alpers, B. Holt, K. Zeng, Oil spill detection by imaging radars: Challenges and pitfalls, Remote Sensing of Environment, 201:133-147, 2017.
  * [4] M. D. Graziano, M. Grasso, M. D’Errico, ”Performance Analysis of Ship Wake Detection on Sentinel-1 SAR Images,” Remote Sensing, 9(11):1107, 2017.
  * [5] W. D. Fisher, T. K. Camp, V. V. Krzhizhanovskaya, ”Anomaly detection in Earth dam and levee passive seismic data using support vector machines and automatic feature selection,” Journal of Computational Science, 20:143-153, 2017.
  * [6] I. Scollar, A. Tabbagh, A. Hesse, and I. Herzog, ”Archaeological prospecting and remote sensing,” Cambridge University Press, 1990.
  * [7] I. S. Reed and X. Yu, ”Adaptive multiple-band CFAR detection of an optical pattern with unknown spectral distribution,” IEEE Transactions on Acoustics, Speech, and Signal Processing, 38:1760–1770, 1990.
  * [8] J. Frontera-Pons, M. A. Veganzones, S. Velasco-Forero, F. Pascal, J.-P. Ovarlez, J. Chanussot, ”Robust Anomaly Detection in Hyperspectral Imaging,” IEEE International Geoscience and Remote Sensing Symposium, Quebec, Canada, 2014.
  * [9] M. A. Veganzones, J. Frontera-Pons, F. Pascal, J.-P. Ovarlez, J. Chanussot, ”Binary Partition Trees-Based Robust Adaptive Hyperspectral RX Anomaly Detection,” IEEE International Conference on Image Processing, Paris, 2014.
  * [10] E. Terreaux, J.-P. Ovarlez and F. Pascal, ”Anomaly Detection and Estimation in Hyperspectral Imaging using RMT tools,” IEEE Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), Cancun, Mexico, 2015
  * [11] J. Frontera-Pons, M. A. Veganzones, F. Pascal and J.-P. Ovarlez, ”Hyperspectral Anomaly Detectors using Robust Estimators,” IEEE Journal of Selected Topics in Applied Earth Observation and Remote Sensing (JSTARS), 9(2):720-731, 2016.
  * [12] A. W. Bitar, L.-F. Cheong, J.-P. Ovarlez, ”Sparse and Low-Rank Matrix Decomposition for Automatic Target Detection in Hyperspectral Imagery,” IEEE Transactions on Geoscience and Remote Sensing, 57(8):5239-5251, August 2019.
  * [13] Y. Haitman, I. Berkovich, S. Havivi, S. Maman, D. G. Blumberg and S. R. Rotman, ”Machine Learning for Detecting Anomalies in SAR Data,” IEEE International Conference on Microwaves, Antennas, Communications and Electronic Systems (COMCAS), Tel-Aviv, Israel, pp. 1-5, 2019.
  * [14] S. Akçay, A. Atapour-Abarghouei, and T. P. Breckon, ”SkipGANomaly: Skip Connected and Adversarially Trained Encoder-Decoder Anomaly Detection,” International Joint Conference on Neural Networks, pp. 1–8, July 2019.
  * [15] T. Schlegl, P. Seeböck, S. M. Waldstein, G. Langs, and U. Schmidt-Erfurth, ”f-AnoGAN: Fast unsupervised anomaly detection with generative adversarial networks,” Medical Image Analysis, 54:30–44, May 2019.
  * [16] S. Sinha, S. Giffard-Roisin, F. Karbou, M. Deschatres, A. Karas, N. Eckert, C. Coléou, and C. Monteleoni, ”Variational Autoencoder Anomaly Detection of Avalanche Deposits in Satellite SAR Imagery,” 10th International Conference on Climate Informatics, pp. 113–119, ACM, Sept. 2020.
  * [17] S. Mabu, S. Hirata, and T. Kuremoto, ”Anomaly Detection Using Convolutional Adversarial Autoencoder and One-class SVM for Landslide Area Detection from Synthetic Aperture Radar Images,” Journal of Robotics, Networking and Artificial Life, 8(2):139, 2021.
  * [18] C.-A. Deledalle, L. Denis, S. Tabti, and F. Tupin, ”Iterative Weighted Maximum Likelihood Denoising With Probabilistic Patch-Based Weights,” IEEE Transactions on Image Processing, 18(12):2661–2672, 2009.
  * [19] C.-A. Deledalle, L. Denis, S. Tabti, and F. Tupin, ”NL-InSAR: Nonlocal Interferogram Estimation,” IEEE Transactions on Geoscience and Remote Sensing, 49(4):1441–1452, 2011.
  * [20] C.-A. Deledalle, L. Denis, F. Tupin, A. Reigber, and M. Jäger, ”NL-SAR: A Unified Nonlocal Framework for Resolution-Preserving (Pol)(In)SAR Denoising,” IEEE Transactions on Geoscience and Remote Sensing, 53(4):2021–2038, 2015.
  * [21] C.-A. Deledalle, L. Denis, S. Tabti, and F. Tupin, ”MuLoG, or how to apply Gaussian denoisers to multi-channel SAR speckle reduction?,” IEEE Transactions on Image Processing,” 26(9):4389–4403, 2017.
  * [22] E. Dalsasso, L. Denis, and F. Tupin, ”SAR2SAR: A semi-supervised despeckling algorithm for SAR images,” IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 14:4321–4329, 2021.
  * [23] E. Dalsasso, L. Denis, M. Muzeau, and F. Tupin, ”Self-supervised training strategies for SAR image despeckling with Deep Neural Networks,” 14th European Conference on Synthetic Aperture Radar (EUSAR), Leipzig, Germany, 2022.
  * [24] E. Dalsasso, L. Denis, and F. Tupin, ”As if by magic: Self-supervised training of deep despeckling networks with MERLIN,” IEEE Transactions on Geoscience and Remote Sensing, 60:1–13, 2022.
  * [25] J. S. Lee E. and Pottier, Polarimetric Radar Imaging: From Basics to Applications, CRC Press: Boca Raton, FL, 2009.
  * [26] J. W. Goodman, ”Some fundamental properties of speckle,” Journal of the Optical Society of America, 66(11):1145–1150, Nov. 1976.
  * [27] R. Baqué, P. Dreuillet, and H. Oriot, ”SETHI: Review of 10 years of development and experimentation of the remote sensing platform,” 2019 International Radar Conference (RADAR), pp. 1–5, 2019.
  * [28] S. Angelliaume, X. Ceamanos, F. Viallefont-Robinet, R. Baqué, P. Déliot, and V. Miegebielle, “Hyperspectral and radar airborne imagery over controlled release of oil at sea,” Sensors, 17(8):1772, Aug. 2017.
  * [29] L. Pallotta, ”Reciprocity Evaluation in Heterogeneous Polarimetric SAR Images,” IEEE Geoscience and Remote Sensing Letters, 19:1-5, 2022.
  * [30] J. Lehtinen, J. Munkberg, J. Hasselgren, S. Laine, T. Karras, M. Aittala, and T. Aila, ”Noise2Noise: Learning Image Restoration without Clean Data,” Proceedings of the 35th International Conference on Machine Learning, PMLR 80:2965-2974, 2018.
  * [31] A. Makhzani, J. Shlens, N. Jaitly, I. Goodfellow, and B. Frey, ”Adversarial Autoencoders,” International Conference on Learning Representations, San Juan, Puerto Rico, May 2016.
  * [32] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. WardeFarley, S. Ozair, A. Courville, and Y. Bengio, ”Generative adversarial networks,” Communications of the ACM, 63:139-144, Oct. 2020.
  * [33] A. Mian, G. Ginolhac, J.-P. Ovarlez, A. Breloy, F. Pascal, ”An Overview of Covariance-based Change Detection Methodologies in Multivariate SAR Image Time Series,” Chapter 3 in Change Detection and Image Time Series Analysis, Unsupervised Methods, Vol. 1, ISTE/Wiley Encyclopedia of Sciences - Remote Sensing Imagery, A. M. Atto, F. Bovolo and L. Bruzzone (eds.), 2021.
  * [34] L. M. Novak, ”Coherent Change Detection for Multi-Polarization SAR,” Conference Record of the Thirty-Ninth Asilomar Conference on Signals, Systems and Computers, pp. 568-573, 2005.
  * [35] W. Förstner and B. Moonen, A Metric for Covariance Matrices, Book chapter in: Grafarend, E.W., Krumm, F.W., Schwarze, V.S. (eds) Geodesy-The Challenge of the 3rd Millennium. Springer, Berlin, Heidelberg, 2003.
  * [36] I. L. Dryden, A. Koloydenko, and D. Zhou, ”Non-Euclidean Statistics for Covariance Matrices, with Applications to Diffusion Tensor Imaging.” The Annals of Applied Statistics, 3(3):1102–23, 2009.
  * [37] A. Radford, L. Metz, and S. Chintala, “Unsupervised Representation Learning with Deep Convolutional Generative Adversarial Networks,” International Conference on Learning Representations, San Juan, Puerto Rico, 2016.
  * [38] D. P. Kingma and J. Ba, “Adam: A Method for Stochastic Optimization,” International Conference on Learning Representations, San Diego, CA, USA, May 2017.
  * [39] L. N. Smith, ”Cyclical Learning Rates for Training Neural Networks,” arXiv:1506.01186, 2015.


# Regularization of the Hill four-body problem with oblate bodies

Edward Belbruno  Yeshiva University, Department of Mathematical Sciences, New
York, NY 10016, USA<EMAIL_ADDRESS>,  Marian Gidea  Yeshiva
University, Department of Mathematical Sciences, New York, NY 10016, USA
<EMAIL_ADDRESS>and  Wai-Ting Lam  Florida Atlantic University, Department
of Mathematical Sciences, Boca Raton, FL 33431, USA<EMAIL_ADDRESS>
###### Abstract.

We consider the Hill four-body problem where three oblate, massive bodies form
a relative equilibrium triangular configuration, and the fourth, infinitesimal
body orbits in a neighborhood of the smallest of the three massive bodies. We
regularize collisions between the infinitesimal body and the smallest massive
body, via McGehee coordinate transformation. We describe the corresponding
collision manifold and show that it undergoes a bifurcation when the
oblateness coefficient of the small massive body passes through the zero
value.

E.B. and M.G. were partially supported by NSF grant DMS-1814543.

W-T.L. was partially supported by NSF grant DMS-1814543 and DMS-2138090.

##  1\. Introduction

We consider the Hill approximation of the circular restricted four-body
problem with oblate bodies, on the motion of an infinitesimal body under the
gravitational influence of three massive bodies of oblate shapes; the three
bodies are assumed to be in a relative equilibrium triangular configuration,
and the motion of the infinitesimal body is assumed to take place in a
neighborhood of the smallest of the three bodies, which we think of as an
asteroid. See [BGCG+20]. The resulting gravitational field in the Hill
approximation contains a non-Newtonian term which depends on the oblateness
coefficient of the asteroid. We use McGehee coordinates to regularize
collisions between the infinitesimal body and the asteroid, which amounts to
blowing up the collision set to a manifold that captures the dynamics in the
singular limit. (Note that, due to the non-Newtonian term in the potential,
the Levi-Civita regularization does not apply to this setting.) We describe
the collision manifold and the regularized dynamics in a neighborhood of it.
We show that each collision solution is branch regularizable, and each
extension of a collision solution is a reflection. We also show that the
collision manifold is not block regularizable. Moreover, we show that the
collision manifold undergoes a double saddle-node bifurcation as the
oblateness coefficient of the asteroid passes through the zero value. When the
shape of the asteroid becomes prolate, no collisions between the infinitesimal
body and the asteroid are possible.

The four body system that we consider here can be viewed as a model for the
Sun-Jupiter-Hektor-Skamandrios system; Hektor is a Juipter’s trojan asteroid,
and Skamandrios is a moonlet of Hektor. Hektor’s shape can be approximated by
a dumb-bell figure and has one of the largest oblateness coefficients amongst
objects of similar size in the solar system [Des15]. The moonlet Skamandrios
appears to have a complicated orbit, which is close to 1:10 and 2:21
orbit/spin resonances; a small change could potentially eject the moonlet or
make it collide with the asteroid [MDCR+14]. This justifies our interest in
understanding collision orbits.

McGehee coordinate transformation was introduced in [McG81] to regularize
collisions in a central force field of the form $U({\bf x})=|{\bf
x}|^{-\alpha}$, where ${\bf x}\in\mathbb{R}^{2}$ and $\alpha>0$. He also
introduced the concept of branch regularization. A solution is branch
regularizable if it has a unique real analytic extension past the collision.
Branch regularization concerns the extension of individual solutions. The
concept of block regularization considers collective extensions of solutions;
it was introduced by Easton in [Eas71] who referred to it as ‘regularization
by surgery’. A flow is called block regularizable if it is diffeomorphic to
the trivial parallel flow in a deleted neighborhood of the collision set.

McGehee regularization can also be applied to quasi-homogeneous central force
fields of the form $U({\bf x})=\gamma_{1}|{\bf
x}|^{-\alpha_{1}}+\gamma_{1}|{\bf x}|^{-\alpha_{2}}$, with
$\gamma_{1},\gamma_{2},\alpha_{1},\alpha_{2}>0$; see [SSM00]. Belbruno used
McGehee transformation to regularize collisions with a black hole in [BP11] in
order to establish the relationship between the null geodesic structure of the
Schwarzschild black hole solution, and the corresponding inverse-cubic
Newtonian central force problem. Belbruno and collaborators also used the
McGehee transformation to study the regularizability of the big bang
singularity, including the case when random perturbations modeled by Brownian
motion are present in the system [Bel13, XB14, BX18]. Other applications of
related interest include [DMS00, GM14, ElB09, ORS22].

A contribution of our work is that we perform McGehee regularization of
collisions in a four-body problem (rather than in a central force field),
where the non-Newtonian part of the gravitational potential is owed to the
shape of the body. As a matter of fact, our work assumes a more general
setting, of a Hill four-body problem with a general quasi-homogenous
potential, which includes the oblateness effect as a particular case. Another
contribution is that we perform a bifurcation analysis as the oblateness
coefficient varies, with the surprising conclusion that collisions cease to
occur as we switch from oblate to prolate shape.

##  2\. Background

###  2.1. Hill four-body problem with oblate bodies

In this section we describe the Hill approximation of the circular restricted
four-body problem with oblate massive bodies. This problem concerns the
dynamics of an infinitesimal body (particle) moving in a plane under the
gravitational influence of three oblate bodies of masses $m_{1}>m_{2}>m_{3}$,
but without influencing their motion. We refer to these three bodies as
primary, secondary, and tertiary, respectively. We express the gravitational
potential of each body in terms of spherical harmonics truncated up to second
order zonal harmonic, that is,

(2.1)
$V_{i}(x_{1},x_{2},x_{3})=\frac{m_{i}}{r}+\frac{m_{i}}{r}\left(\frac{R_{i}}{r}\right)^{2}\left(\frac{C^{i}_{20}}{2}\right)\left(3\left(\frac{x_{3}}{r}\right)^{2}-1\right)$

where $r=(x_{1}^{2}+x_{2}^{2}+x_{3}^{2})^{1/2}$, $R_{i}$ is the average radius
of the $i$-th body, and the gravitational constant is normalized to $1$. The
dimensionless quantity $C^{i}_{20}$ is the coefficient of the zonal harmonic
of order $2$, with $C^{i}_{20}<0$ for an oblate body, $C^{i}_{20}=0$ for a
spherical body, and $C^{i}_{20}>0$ for a prolate body. Further, we denote
$C_{i}=C^{i}_{20}R_{i}^{2}/2$.

For the circular restricted four-body problem, the assumption is that the
three massive bodies are in a relative equilibrium configuration, that is,
they move on circular orbits around their center of mass while preserving
their mutual distances constant over time. In the case when the bodies have no
oblateness, the only non-collinear relative equilibrium configuration is the
Lagrangian equilateral triangle. When the bodies are oblate, the gravitational
potential is no longer Newtonian, and the relative equilibrium is no longer an
equilateral triangle. It has been shown in [BGCG+20] that there is a unique
relative equilibrium which is a scalene triangle. Such triangle has the
property that the body with the larger $C_{i}$ is opposite to the longer side
of the triangle. We normalize the units of distance so that the distance
between $m_{1}$ and $m_{2}$ is set to $1$, and we let $u_{1}$ be the distance
from $m_{1}$ to $m_{3}$, and $u_{2}$ be the distance from $m_{2}$ to $m_{3}$.
See Fig. 1. The sides $u_{1}$ and $u_{2}$ are uniquely determined by the
implicit equations

(2.2)
$1-3C_{12}=\frac{1}{u_{1}^{3}}-\frac{3C_{13}}{u_{1}^{5}}=\frac{1}{u_{2}^{3}}-\frac{3C_{23}}{u_{2}^{5}},$

where we denote $C_{ij}=C_{i}+C_{j}$.

Figure 1. Scalene triangle relative equilibrium.

Given such a relative equilibrium configuration, the motion of the particle in
a vicinity of the tertiary is described by the Hamiltonian of the circular
restricted four-body problem (see, e.g., [BGCG+20]). However, the
corresponding Hamiltonian equations are difficult to treat analytically.
Therefore we consider below the Hill approximation of the circular restricted
four-body problem. This is derived by rescaling the distances by a factor of
$m_{3}^{1/3}$, writing the associated Hamiltonian in the rescaled coordinates
as a power series in $m_{3}^{1/3}$, and neglecting all the terms of order
$O(m_{3}^{1/3})$ in the expansion. The oblateness coefficient $C_{i}$ also
gets rescaled to $c_{i}=m_{3}^{-2/3}C^{i}_{20}R_{i}^{2}/2$. This procedure
yields an approximation of the motion of the particle in an
$O(m_{3}^{1/3})$-neighborhood of the tertiary, while the primary and the
secondary are ‘sent to infinity’. We obtain a much simpler Hamiltonian then
the one for the circular restricted four-body problem, for which the
contribution of the primary and of the secondary to the gravitational
potential is given by a quadratic polynomial. Specifically, the Hamiltonian of
the Hill four-body problem relative to some convenient co-rotating frame is
given by

(2.3)
$\begin{split}H=&\frac{1}{2}(y_{1}^{2}+y_{2}^{2}+y_{3}^{2})+x_{2}y_{1}-x_{1}y_{2}\\\
&\quad+\left(\frac{1-\lambda_{2}}{2}\right)x_{1}^{2}+\left(\frac{1-\lambda_{1}}{2}\right)x_{2}^{2}+\frac{1}{2}\left(\frac{(1-\mu)}{u_{1}^{3}}+\frac{\mu}{u_{2}^{3}}\right)x_{3}^{2}\\\
&\quad-\left(\frac{(1-\mu)c_{1}}{u_{1}^{3}}\right)\left(3\left(\frac{x_{3}}{u_{1}}\right)^{2}-1\right)-\left(\frac{\mu
c_{2}}{u_{2}^{3}}\right)\left(3\left(\frac{x_{3}}{u_{2}}\right)^{2}-1\right)\\\
&\quad-\frac{1}{(x_{1}^{2}+x_{2}^{2}+x_{3}^{2})^{1/2}}-\frac{c_{3}}{(x_{1}^{2}+x_{2}^{2}+x_{3}^{2})^{3/2}}\left(\frac{3x_{3}^{2}}{x_{1}^{2}+x_{2}^{2}+x_{3}^{2}}-1\right).\end{split}$

where $\mu=\frac{m_{2}}{m_{1}+m_{2}}$, and $\lambda_{1}$ and $\lambda_{2}$ are
given by the following formulas

(2.4)
$\begin{split}\lambda_{1}=&\frac{1}{2}\left(2-\frac{2(1-\mu)}{u_{1}^{5}}-\frac{2\mu}{u_{2}^{5}}+\frac{3(1-\mu)}{u_{1}^{3}}+\frac{3\mu}{u_{2}^{3}}-\frac{3}{u_{1}^{3}u_{2}^{3}}\sqrt{\Delta}\right),\\\
\lambda_{2}=&\frac{1}{2}\left(2-\frac{2(1-\mu)}{u_{1}^{5}}-\frac{2\mu}{u_{2}^{5}}+\frac{3(1-\mu)}{u_{1}^{3}}+\frac{3\mu}{u_{2}^{3}}+\frac{3}{u_{1}^{3}u_{2}^{3}}\sqrt{\Delta}\right),\end{split}$

where

$\Delta=(\mu
u_{1}^{3}+(1-\mu)u_{2}^{3})^{2}-\mu(1-\mu)u_{1}u_{2}\left(-u_{1}^{4}-u_{2}^{4}+2u_{1}^{2}+2u_{2}^{2}+2u_{1}^{2}u_{2}^{2}-1\right).$

When we restrict to the planar problem ($x_{3}=0$) the Hamiltonian becomes

(2.5)
$\begin{split}H=&\frac{1}{2}(y_{1}^{2}+y_{2}^{2})+x_{2}y_{1}-x_{1}y_{2}+\left(\frac{1-\lambda_{2}}{2}\right)x_{1}^{2}+\left(\frac{1-\lambda_{1}}{2}\right)x_{2}^{2}\\\
&\quad-\frac{1}{(x_{1}^{2}+x_{2}^{2})^{1/2}}+\frac{c_{3}}{(x_{1}^{2}+x_{2}^{2})^{3/2}},\end{split}$

where the constant terms $\frac{(1-\mu)c_{1}}{u_{1}^{3}}$ and $\frac{\mu
c_{2}}{u_{2}^{3}}$ were dropped, as they do not appear in the Hamiltonian
equations. We note that in the planar problem the oblateness of the primary
and the secondary plays no role.

We denote by $\mathbf{M}_{h}$ the $3$-dimensional energy manifold

(2.6) $\mathbf{M}_{h}=\\{H=h\\}.$

######  Remark 2.1.

An example of a system that can be modeled by the Hill four-body problem is
the Sun-Jupiter-Hektor-Skamandrios system [BGCG+20]. Hektor is a Jupiter
Trojan, which is approximately located at Lagrangian point $L_{4}$ of the Sun-
Jupiter system, thus forming an approximate triangular relative equilibrium
configuration with Sun and Jupiter. Hektor is the biggest Jupiter Trojan and
has one of the largest values of the oblateness coefficients among the objects
of its size in in the Solar system. Hektor’s moonlet, Skamandrios, can be
viewed as the fourth, infinitesimal body. In this case, the constants that
appear in (2.3) are $c_{3}=-1.327161\times 10^{-7}$, $\mu=0.0009533386$,
$u_{1}=1-5.94154\times 10^{-11}$, $u_{2}=1-1.99318\times 10^{-12}$,
$\lambda_{1}=0.002144$, and $\lambda_{2}=2.997855$.

###  2.2. Branch and block regularization

We give a brief review of branch and block regularization following [McG81].

For a differential equation

(2.7) $\dot{\mathbf{x}}=F(\mathbf{x})$

with $F$ a real analytic vector field on some open set
$U\subseteq\mathbb{R}^{n}$, and $\dot{}=\frac{d}{dt}$. The standard existence
and uniqueness theorem for ODE’s gives for each initial condition
$\mathbf{x}(0)\in U$ a unique, real analytic solution $\mathbf{x}(t)$ defined
on a maximal interval $(t^{-},t^{+})$ with $-\infty\leq
t^{-}<0<t^{+}\leq+\infty$. Solutions for which $-\infty<t^{-}$ or
$t^{+}<+\infty$ are said to have a singularity at $t^{*}=t^{-}$ or
$t^{*}=t^{+}$, respectively.

We briefly describe the concept of branch regularization.

If $\mathbf{x}_{1}(t)$ and $\mathbf{x}_{2}(t)$ are solutions of (2.7), with
$\mathbf{x}_{1}$ ending in a singularity at time $t^{*}$ and $\mathbf{x}_{2}$
beginning in a singularity at $t^{*}$, and there exists a multivalued analytic
complex function having a branch at $t^{*}$ and extending both
$\mathbf{x}_{1}$ and $\mathbf{x}_{2}$ when we regard the time $t$ as complex,
then $\mathbf{x}_{1}$, $\mathbf{x}_{2}$ are said to be branch extensions of
one another at $t^{*}$.

A solution $\mathbf{x}(t)$ of equation(2.7) with a singularity at $t^{*}$ is
said to be branch regularizable at $t^{*}$ if it has a unique branch extension
at $t^{*}$. The extension is called a ‘reflection’ if the velocity vector
reverses direction at collision, and is called a ‘transmission’ if the
direction of the velocity vector is preserved at collision.

The equation (2.7) is said to be branch regularizable if every solution is
branch regularizable at every singularity.

Now, consider the motion of a single particle in a potential field given by

(2.8) $U(\mathbf{x})=|\mathbf{x}|^{-\alpha}\textrm{ with
$\mathbf{x}\in\mathbb{R}^{2}$.}$

The equation of motion is given by the second order equation

$\ddot{\mathbf{x}}=\nabla U(\mathbf{x}),$

or equivalently, by the first order system

$\left\\{\begin{split}\dot{\mathbf{x}}=&\mathbf{y},\\\
\dot{\mathbf{y}}=&-\alpha|\mathbf{x}|^{-\alpha-2}\mathbf{x}.\end{split}\right.$

Let $\beta=\frac{\alpha}{2}$ and
$\gamma=\frac{1}{1+\beta}=\frac{2}{2+\alpha}$.

We recall the following result from [McG81]:

######  Theorem 2.2.

A collision solution for the potential (2.8) is branch regularizable if and
only if $\gamma=\frac{p}{q}$ with $p<q$ positive integers,
$\textrm{gcd}(p,q)=1$, and $q$ odd.

Moreover, if $p$ is even the extension solution is a ‘reflection’, and when
$p$ is odd the extension solution is a ‘transmission’.

In [SSM00] this result has been extended for quasi-homogeneous potentials of
the form

(2.9)
$U(\mathbf{x})=\gamma_{1}|\mathbf{x}|^{-\alpha_{1}}+\gamma_{2}|\mathbf{x}|^{-\alpha_{2}}\textrm{
with }\mathbf{x}\in\mathbb{R}^{2},$

where $\gamma_{1},\gamma_{2}>0$, $\alpha_{2}>\alpha_{1}>0$.

######  Theorem 2.3.

A collision solution for the potential (2.9) is branch regularizable if and
only if both

$\frac{2}{2+\max(\alpha_{1},\alpha_{2})},\textrm{ and
}\frac{\min(\alpha_{1},\alpha_{2})}{2+\max(\alpha_{1},\alpha_{2})}$

are of the form $\frac{p}{q}$ with $p<q$ positive integers,
$\textrm{gcd}(p,q)=1$, and $q$ odd.

We now describe the concept of block regularization. Denote by
$\phi^{t}=\phi(\cdot,t)$ the flow of (2.7).

A compact invariant set $\mathbf{N}\subseteq U$ is called isolated if there
exists an open set $V\subseteq U$ – referred to as an isolating neighborhood –
such that $\mathbf{N}\subset V$ is the maximal invariant subset of $V$.

Let $\mathbf{B}\subseteq U$ be a compact set with non-empty interior, and
assume that the boundary $\mathbf{b}=\partial\mathbf{B}$ of $\mathbf{B}$ is a
smooth submanifold. Define

$\begin{split}\mathbf{b}^{+}=&\\{\mathbf{x}\in\mathbf{b}\,|\,\phi(\mathbf{x},(-\varepsilon,0))\cap\mathbf{B}=\emptyset,\textrm{
for some }\varepsilon>0\\},\\\
\mathbf{b}^{-}=&\\{\mathbf{x}\in\mathbf{b}\,|\,\phi(\mathbf{x},(0,\varepsilon))\cap\mathbf{B}=\emptyset,\textrm{
for some }\varepsilon>0\\},\\\
\mathbf{t}=&\\{\mathbf{x}\in\mathbf{b}\,|\,\dot{\phi}(\mathbf{x},0)\textrm{ is
tangent to }\mathbf{b}\\}.\end{split}$

The set $\mathbf{B}$ is called an isolating block if
$\mathbf{b}^{+}\cap\mathbf{b}^{-}=\mathbf{t}$.

If $\mathbf{N}$ is an isolated invariant set, we say that $\mathbf{B}$
isolates $\mathbf{N}$ if the interior set $\textrm{Int}(\mathbf{B})$ of
$\mathbf{B}$ is an isolating neighborhood for $\mathbf{N}$. For every isolated
invariant set $\mathbf{N}$ there exists an isolating block which isolates
$\mathbf{N}$. If $\mathbf{B}$ is an isolating block, then there exists an
isolated invariant set $\mathbf{N}$ (possibly empty) which is isolated by
$\mathbf{B}$. See [CE71].

The asymptotic sets to $\mathbf{N}$ are defined by

$\begin{split}\mathbf{a}^{+}=&\\{\mathbf{x}\in\mathbf{b}^{+}\,|\,\phi(\mathbf{x},(0,+\infty))\subset\mathbf{B}\\}\\\
\mathbf{a}^{-}=&\\{\mathbf{x}\in\mathbf{b}^{-}\,|\,\phi(\mathbf{x},(-\infty,0))\subset\mathbf{B}\\}\\\
\end{split}$

The map across the block is defined as

$\begin{split}\Phi:\mathbf{b}^{+}\setminus\mathbf{a}^{+}\to\mathbf{b}^{-}\setminus\mathbf{a}^{-},\\\
\Phi(\mathbf{x})=\phi(\mathbf{x},T(\mathbf{x})),\end{split}$

where $T(\mathbf{x})=\inf\\{t>0\,|\,\phi(\mathbf{x},t)\not\in\mathbf{B}\\}$ is
the time spent inside the block.

If $\mathbf{B}$ is an isolating block, then the application $\Phi$ is a
diffeomorphism. See [CE71].

An isolating block $\mathbf{B}$ is said to be trivializable if the map $\Phi$
extends uniquely to a diffeomorphism from $\mathbf{b}^{+}$ to
$\mathbf{b}^{-}$.

The theory of isolating blocks can be applied to singularities by essentially
replacing the role of an isolated invariant set $\mathbf{N}$ as above with the
set of singularities, as we shall see below.

In Section 4 we will see that, going through regularized coordinates and time
rescaling, the set of singularities for (2.3), which consists of the origin,
gets transformed into an invariant set, which is in fact a manifold (referred
to as a collision manifold).

Let $F(\mathbf{x})$ be a vector field defined on $U\setminus\mathbf{N}$, where
$\mathbf{N}$ is a compact set representing the singularities of the vector
field. Let $\mathbf{B}\subseteq U$ be compact set with non-empty interior,
such that $\mathbf{b}=\partial\mathbf{B}$ is a smooth submanifold, and with
$\mathbf{b}\cap\mathbf{N}=\emptyset$. Define the subsets
$\mathbf{b}^{+},\mathbf{b}^{-}\subset\mathbf{b}$ in the same way as above.
Under these conditions, the definition of an isolating block is the same as
before.

The orbit through a point $\mathbf{x}$ is defined by

$O(\mathbf{x})=\\{\phi(\mathbf{x},t)\,|\,\phi(\mathbf{x},t)\textrm{ is defined
}\\}.$

That is, there are no invariant sets in $\mathbf{B}$.

An isolating block $\mathbf{B}$ is said to isolate the singularity set
$\mathbf{N}$ if $\mathbf{N}\subset\textrm{Int}(\mathbf{B})$ and if
$O(\mathbf{x})\not\subset\mathbf{B}$ for all
$\mathbf{x}\in\mathbf{B}\setminus\mathbf{N}$.

The asymptotic sets to $\mathbf{N}$ are defined by

$\begin{split}\mathbf{a}^{+}=&\\{\mathbf{x}\in\mathbf{b}^{+}\,|\,\phi(\mathbf{x},t)\in\mathbf{B}\textrm{
for all $t\geq 0$ for which $\phi(\mathbf{x},t)$ is defined}\\}\\\
\mathbf{a}^{-}=&\\{\mathbf{x}\in\mathbf{b}^{-}\,|\,\phi(\mathbf{x},t)\in\mathbf{B}\textrm{
for all $t\leq 0$ for which $\phi(\mathbf{x},t)$ is defined}\\}.\end{split}$

We define the map across the block
$\Phi:\mathbf{b}^{+}\setminus\mathbf{a}^{+}\to\mathbf{b}^{-}\setminus\mathbf{a}^{-}$
as before.

The singularity set $\mathbf{N}$ is said to be block regularizable if there
exists a trivializable block $\mathbf{B}$ which isolates $\mathbf{N}$ .

Regarding block regularization we recall the following result from [McG81]:

######  Theorem 2.4.

A collision set for the potential (2.8) is block regularizable if and only if
$\beta=1-\frac{1}{n}$ for $n$ positive integer.

##  3\. Main result

######  Theorem 3.1.

For the system (2.5) with oblate tertiary, i.e., $c_{3}<0$, each collision
solution is branch regularizable, and each extension of a collision solution
is a reflection. The collision manifold is not block regularizable.

At $c_{3}=0$ the reduced system of equations associated to the collision
manifold undergoes a double saddle-node bifurcation. For $c_{3}=0$, the
collision manifold is branch and block regularizable.

For the system (2.5) with a prolate tertiary, i.e., $c_{3}>0$, there are no
collisions.

The collision manifold and the corresponding reduced system of equations are
described in Section 5. The proof of this result is provided in the following
sections.

We note here that the bifurcation observed for this system is very different
from the one described by [McG81] for the system (2.8), which occurs when the
parameter $\alpha$ is varied.

##  4\. McGehee transformation

We rewrite the Hamiltonian (2.5) in a simpler form

(4.1)
$\begin{split}H=&\frac{1}{2}(y_{1}^{2}+y_{2}^{2})+x_{2}y_{1}-x_{1}y_{2}+Ax_{1}^{2}+Bx_{2}^{2}\\\
&\quad-\frac{1}{(x_{1}^{2}+x_{2}^{2})^{\nu/2}}-\frac{c}{(x_{1}^{2}+x_{2}^{2})^{\alpha/2}},\end{split}$

where $\dot{}=\frac{d}{dt}$, and $1\leq\nu<\alpha$. The corresponding
potential is quasi-homogeneous.

In the case of the potential (2.5), we have $\nu=1$, $\alpha=3$,
$A=\frac{1-\lambda_{2}}{2}$, $B=\frac{1-\lambda_{1}}{2}$, and $c=-c_{3}$.

We identify $x,y\in\mathbb{R}^{2}$ with the complex numbers $x_{1}+ix_{2}$,
$y_{1}+iy_{2}$, respectively. The corresponding Hamilton equations are

(4.2) $\begin{split}\dot{x}=&\frac{\partial H}{\partial y}=y-ix,\\\
\dot{y}=&-\frac{\partial H}{\partial x}=-\frac{\nu
x}{|x|^{\nu+2}}-\frac{\alpha cx}{|x|^{\alpha+2}}-iy-Tx,\end{split}$

where $T$ is the real-linear transformation given by
$T(x_{1}+ix_{2})=2Ax_{1}+2Bx_{2}i$, and $|x|=(x_{1}^{2}+x_{2}^{2})^{1/2}$.

We perform a coordinate change to new real coordinates $(r,\theta,v,w)$, with
$r>0$ and $\theta\in\mathbb{T}^{1}$, defined as follows

(4.3) $\begin{split}x=&r^{\gamma}e^{i\theta},\\\
y=&r^{-\gamma\beta}(v+iw)e^{i\theta},\end{split}$

where

$\beta=\frac{\alpha}{2},\textrm{ and
}\gamma=\frac{1}{\beta+1}=\frac{2}{\alpha+2}.$

Writing (4.3) in terms of components we have

(4.4) $\begin{split}x_{1}=&r^{\gamma}\cos\theta,\\\
x_{2}=&r^{\gamma}\sin\theta,\\\
y_{1}=&r^{-\gamma\beta}(v\cos\theta-w\sin\theta),\\\
y_{2}=&r^{-\gamma\beta}(v\sin\theta+w\cos\theta).\end{split}$

The new coordinates $(r,\theta,v,w)$ in terms of the old coordinates
$(x_{1},x_{2},y_{1},y_{2})$ are given by

(4.5) $\begin{split}r=&|x|^{\frac{1}{\gamma}},\\\ \theta=&\arg(x),\\\
v=&r^{\gamma\beta}(y_{1}\cos\theta+y_{2}\sin\theta),\\\
w=&r^{\gamma\beta}(-y_{1}\sin\theta+y_{2}\cos\theta).\end{split}$

The new coordinates are known as the McGehee coordinates [McG81]. We rewrite
the Hamiltonian equations (4.2) in the new coordinates and equate the real and
imaginary parts on the two sides. From

$\begin{split}\dot{x}=&\left[\gamma
r^{\gamma-1}\dot{r}+ir^{\gamma}\dot{\theta}\right]e^{i\theta},\\\
\frac{\partial H}{\partial
y}=&\left[r^{-\beta\gamma}(v+iw)-ir^{\gamma}\right]e^{i\theta},\\\
\dot{y}=&\left[-\beta\gamma
r^{-\beta\gamma-1}\dot{r}(v+iw)+r^{-\beta\gamma}(\dot{v}+i\dot{w})+r^{-\beta\gamma}(v+iw)i\dot{\theta}\right]e^{i\theta}\\\
=&\left[-\beta
r^{-\beta\gamma-1}v(v+iw)+r^{-\beta\gamma}(\dot{v}+i\dot{w})+r^{-\beta\gamma}(v+iw)i(r^{-1}w-1)\right]e^{i\theta},\\\
-\frac{\partial H}{\partial x}=&\left[-\frac{\nu
r^{\gamma}}{r^{\gamma(\nu+2)}}-\frac{\alpha
cr^{\gamma}}{r^{\gamma(\alpha+2)}}-i\frac{v+iw}{r^{\gamma\beta}}\right]e^{i\theta}-2Ar^{\gamma}\cos\theta-i2Br^{\gamma}\sin\theta\\\
&=\left[-\nu r^{\gamma(-1-\nu)}-\alpha
cr^{\gamma(-1-\alpha)}-i(v+iw)r^{-\gamma\beta}\right.\\\
&\quad\left.+2r^{\gamma}(-A\cos^{2}\theta-B\sin^{2}\theta)+i2r^{\gamma}(A-B)\sin\theta\cos\theta\right]e^{i\theta}\end{split}$

we obtain

(4.6) $\begin{split}\dot{r}=&(\beta+1)v,\\\ \dot{\theta}=&r^{-1}w-1,\\\
\dot{v}=&\frac{\beta v^{2}+w^{2}-\alpha
c}{r}-\frac{\nu}{r^{\gamma(\nu+2)-1}}-2Ar\cos^{2}\theta-2Br\sin^{2}\theta,\\\
\dot{w}=&\frac{(\beta-1)vw}{r}+2(A-B)r\sin\theta\cos\theta.\end{split}$

In the above, after equating $\dot{x}=\frac{\partial H}{\partial y}$ we obtain
$\dot{r}=(\beta+1)v$ and $\dot{\theta}=r^{-1}w-1$, which we substitute in the
equation for $\dot{y}$. We also use that $\frac{1-\gamma}{\gamma}=\beta$,
$-\gamma\beta=\gamma-1$, and $\alpha=2\beta$. The fact that $T$ is real-linear
transformation but not-complex linear is taken into account when factoring out
$e^{i\theta}$ in the equation for $-\frac{\partial H}{\partial x}$ by
expressing $Tx=(Txe^{-i\theta})e^{i\theta}$.

The equations (4.6) have a singularity at $r=0$. We remove the singularity by
introducing a new time parameter $\tau$ given by

(4.7) $dt=rd\tau.$

The equations (4.6) expressed in terms of the new time $\tau$ become

(4.8) $\begin{split}r^{\prime}=&(\beta+1)vr,\\\ \theta^{\prime}=&w-r,\\\
v^{\prime}=&(\beta v^{2}+w^{2}-\alpha c)-\nu
r^{2-\gamma(\nu+2)}-2Ar^{2}\cos^{2}\theta-2Br^{2}\sin^{2}\theta,\\\
w^{\prime}=&(\beta-1)vw+2(A-B)r\sin\theta\cos\theta,\end{split}$

where ${}^{\prime}=\frac{d}{d\tau}$. Since $\nu<\alpha$, we have that
$2-\gamma(\nu+2)>0$. Thus, the obtained differential equations have no
singularity at $r=0$; the singularity has been ‘removed’. We also note that
the terms $\nu
r^{2-\gamma(\nu+2)}-2Ar^{2}\cos^{2}\theta-2Br^{2}\sin^{2}\theta$ and
$2(A-B)r\sin\theta\cos\theta$ tend to $0$ as $r\to 0$, so they can be
neglected for $r$ sufficiently small.

The energy condition $H=h$ in the new coordinates, when we use (4.4), becomes

(4.9) $\begin{split}&\frac{1}{2}r^{-2\gamma\beta}(v^{2}+w^{2})\\\
&+r^{\gamma-\gamma\beta}\sin\theta(v\cos\theta-w\sin\theta)-r^{\gamma-\gamma\beta}\cos\theta(w\cos\theta+v\sin\theta)\\\
&+Ar^{2\gamma}\cos^{2}\theta+Br^{2\gamma}\sin^{2}\theta-r^{-\gamma\nu}-cr^{-\gamma\alpha}=h,\end{split}$

which, after multiplying both sides by $r^{2-2\gamma}=r^{2\gamma\beta}$ yields

(4.10)
$\begin{split}\frac{v^{2}+w^{2}-2c}{2}-rw+r^{2}(A\cos^{2}\theta+B\sin^{2}\theta)-r^{2-\gamma(\nu+2)}=r^{2-2\gamma}\,h.\end{split}$

Thus, the energy manifold $\mathbf{M}_{h}$ in the coordinates $(r,\theta,v,w)$
is defined by (4.10).

When $r=0$ the energy condition (4.10) reduces to

(4.11) $\begin{split}v^{2}+w^{2}-2c=0.\end{split}$

######  Remark 4.1.

In the case of the potential (2.8), one obtains a system of $4$-equations,
similar to (4.8), however that system is partially decoupled – the first two
equations are determined by the last two equations. Also, the energy manifold
$\mathbf{M}_{h}$ projects onto $\\{v^{2}+w^{2}<1\\}$ when $h<0$, onto
$\\{v^{2}+w^{2}>1\\}$ when $h>0$, and onto $\\{v^{2}+w^{2}=1\\}$ when $h=0$.
See[McG81].

In our case the system (4.8) is fully coupled, and there is no obvious
relation between the regions bounded by $v^{2}+w^{2}=2c$ in the $(v,w)$-plane
and the energy $h$.

Writing the energy condition (4.10) as

$\frac{v^{2}+w^{2}-2c}{2}=r^{2-2\gamma}\,h+rw-r^{2}(A\cos^{2}\theta+B\sin^{2}\theta)+r^{2-\gamma(\nu+2)},$

by neglecting the $O(r^{2})$ terms we can obtain that the points in
$\mathbf{M}_{h}$ with $h>0$, $w>0$ and $r$ sufficiently small project onto
$\\{v^{2}+w^{2}>2c\\}$.

Also, provided that $\alpha>2(\nu+1)$, we can neglect $r^{2-\gamma(\nu+2)}$
and obtain that the points in $\mathbf{M}_{h}$ with $h<0$, $w<0$ and $r$
sufficiently small project onto $\\{v^{2}+w^{2}<2c\\}$.

##  5\. Collision manifold

From (4.11), the intersection between the energy manifold $\mathbf{M}_{h}$ and
the $3$-dimensional hyperplane

$\mathbf{Z}=\\{r=0\\}$

is a $2$-dimensional manifold corresponding to collisions

(5.1) $\mathbf{N}_{h}=\\{H=h\\}\cap\\{r=0\\}.$

It is referred to as the _collision manifold_. Thus, from (4.11) we obtain
that

(5.2)
$\mathbf{N}_{h}=\\{(r,\theta,v,w)\,|\,r=0,\theta\in\mathbb{T}^{1},v^{2}+w^{2}=2c\\},$

so the collision manifold is a $2$-dimensional torus provided $c>0$. Note that
it is in fact independent of the energy level $h$.

The collision manifold is an isolated invariant set for the flow of (4.8). If
a trajectory approaches the singularity, i.e., $r\to 0$ as $t\to\pm t^{*}$,
then in the $(r,\theta,v,w)$ coordinates the trajectory approaches the
collision manifold $\mathbf{N}_{h}$ as $\tau(t)\to\mp\infty$. The argument is
the same as in [McG81].

Since $r^{\prime}=0$ when $r=0$, it follows that the set $\mathbf{Z}$ is
invariant under the solutions of the system (4.8). Thus, we can consider the
restriction of (4.8) to $\mathbf{Z}$, which is given by

(5.3) $\begin{split}\theta^{\prime}=&w,\\\ v^{\prime}=&\beta
v^{2}+w^{2}-\alpha c,\\\ w^{\prime}=&(\beta-1)vw.\end{split}$

The dynamics on $\mathbf{Z}$ is the skew product between the dynamics in the
variables $(v,w)$ and the dynamics in $\theta$. The solution of the equation
in $\theta$ is determined by the solutions of the $(v,w)$-subsystem, which is
independent of $\theta$. We refer to the $(v,w)$-subsystem of (5.3) as the
_reduced system associated to the collision manifold_.

Define

(5.4) $K=|w|^{\alpha}|v^{2}+w^{2}-2c|^{1-\beta}.$

We claim that $K$ is an integral of motion for the $(v,w)$-subsystem of (5.3).
Indeed, using (5.3) we obtain

$\begin{split}K^{\prime}=&\alpha|w|^{\alpha-1}w^{\prime}|v^{2}+w^{2}-2c|^{1-\beta}+|w|^{\alpha}|v^{2}+w^{2}-2c|^{-\beta}(2vv^{\prime}+2ww^{\prime})\\\
=&\alpha|w|^{\alpha-1}(\beta-1)vw|v^{2}+w^{2}-2c|^{1-\beta}\\\
&+|w|^{\alpha}|v^{2}+w^{2}-2c|^{-\beta}(2v(\beta v^{2}+w^{2}-\alpha
c)+2w(\beta-1)vw)\\\
=&|w|^{\alpha-1}|v^{2}+w^{2}-2c|^{-\beta}(\beta-1)vw\left[\alpha(v^{2}+w^{2}-2c)\right.\\\
&\left.\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad-2(\beta
v^{2}+w^{2}-\alpha c)-2(\beta-1)w^{2}\right]\\\ =&0.\end{split}$

By (5.2), the collision manifold $\mathbf{N}_{h}$ intersects the $(v,w)$ plane
along the $0$-level set of the integral $K$.

We now describe the geometry of the $(v,w)$-subsystem. The equilibrium points
are $S_{\pm}=(0,\pm\sqrt{\alpha c})$ and $Q_{\pm}=(\pm\sqrt{2c},0)$. The
circle

$\mathbf{C}=\\{v^{2}+w^{2}=2c\\}$

is invariant under the flow of the subsystem, and passes through the points
$Q_{\pm}$.

The eigenvalues of the linearized system at $S_{\pm}$ are

$\pm\sqrt{2(\beta-1)\alpha c},$

and since one is positive and the other is negative, both points are saddle
points. The eigenvalues of the linearized system at $Q_{\pm}$ are

$\pm 2\beta\sqrt{2c},\quad\pm(\beta-1)\sqrt{2c}.$

Both eigenvalues at $Q_{+}$ are positive hence this is a source. Both
eigenvalues at $Q_{-}$ are negative hence this is a sink. The line $w=0$ is
also invariant under the flow, where $v^{\prime}<0$ for $|v|<\sqrt{2c}$ and
$v^{\prime}>0$ for $|v|>\sqrt{2c}$. The phase portrait is shown in Fig. 2.

We obtain the existence of several types of orbits:

  * •

orbits beginning and ending in collision orbits

  * •

orbits that only begin or only end in collision

  * •

swing-by orbits, passing near the singularity

  * •

hyperbolic periodic orbits possessing asymptotic orbits (stable and unstable
manifolds in $\mathbf{M}_{h}$) than begin or end in collision

$\begin{array}[]{cc}\includegraphics[width=216.81pt]{contour_c_1.pdf}&\includegraphics[width=216.81pt]{contour_c_0_1.pdf}\end{array}$

Figure 2. Left: $c=1$. Right: $c=0.1$.

Recall that for the system 2.5 we have $\nu=1$, $\alpha=3$,
$\beta=\frac{3}{2}$, and $\gamma=\frac{2}{5}$. By Theorem 2.3 it follows that
each collision solution is branch regularizable. Since $p=2$ is even, each
extension solution is a ‘reflection’.

By examining Fig. 2 we observe that the collision manifold $\mathbf{N}_{h}$ is
not an isolated invariant set, and therefore it is not block regularizable.
This agrees with the case of the potential (2.8), where for $\beta\geq 1$ the
collision manifold $\mathbf{N}$ is not an isolated set.

As $c\to 0$, the two saddles, the source, and the sink disappear through a
double saddle-node bifurcation. See Fig. 3.

For $c=0$ the collision manifold is reduced to a point, and it is both branch
and block regularizable.

We now discuss the case when $c<0$. This describes a situation when the
tertiary is a prolate body, In this case the set of $(v,w)$ with
$v^{2}+w^{2}=2c$ is the empty set. Thus the collision set $\mathbf{N}_{h}$ is
empty. Then the $(v,w)$-subsystem

(5.5) $\begin{split}v^{\prime}=&\beta v^{2}+w^{2}-\alpha c,\\\
w^{\prime}=&(\beta-1)vw,\end{split}$

has the property that $v^{\prime}>0$. The phase portrait is shown in Fig. 3.
In this case there are no collisions.

$\begin{array}[]{cc}\includegraphics[width=216.81pt]{contour_c_0.pdf}&\includegraphics[width=216.81pt]{contour_c_neg.pdf}\end{array}$

Figure 3. Left: $c=0$. Right: $c=-1$.

The physical interpretation is the following. Denoting $c=-\tilde{c}$ where
$\tilde{c}>0$, the Hamiltonian (2.5) becomes

(5.6)
$\begin{split}H=&\frac{1}{2}(y_{1}^{2}+y_{2}^{2})+x_{2}y_{1}-x_{1}y_{2}+Ax_{1}^{2}+Bx_{2}^{2}\\\
&\quad-\frac{1}{|x|^{\nu/2}}+\frac{\tilde{c}}{|x|^{\alpha/2}}.\end{split}$

The term $-\frac{1}{|x|^{\nu/2}}$ in the potential corresponds to an
attractive force, and the term $\frac{\tilde{c}}{|x|^{\alpha/2}}$ corresponds
to a repulsive force. When the particle approaches the tertiary, since
$1\leq\nu<\alpha$ the repulsive force becomes dominating, preventing
collisions between the particle and the tertiary to occur. This situation is
also described in [Saa74].

######  Remark 5.1.

There are several moons in the Solar System that are considered to be
approximately prolate spheroids in shape, for example, Uranus’ moons Cordelia,
Cressida, Desdemona, Juliet, Ophelia, and Rosalind.

######  Remark 5.2.

In the case when $c=0$, the term $-\frac{c}{|x|^{\alpha}}$ in (4.1) vanishes.
Then one can perform the coordinate change (4.3) with $\gamma=\frac{2}{2+\nu}$
and $\beta=\frac{\nu}{2}$, as in [McG81]. The resulting collision set is a
circle as in Fig. 4, and is both branch and block regularizable. It is
interesting that this coordinate change leads to a different collision
manifold from the one above, but nevertheless its branch and block
regularization properties are the same.

Figure 4. Collision manifold for $c=0$ with the coordinate change (4.3) with
$\gamma=\frac{2}{2+\nu}$ and $\beta=\frac{\nu}{2}$.

##  6\. Conclusions

In this paper we provide an explicit McGehee coordinate transformation to
regularize collision in the planar Hill four-body problem with oblate bodies.
This transformation can be used to understand the behavior of collision and
near-collision orbits. In particular, our formulas can be implemented in
numerical integrators to compute orbits that pass close to an oblate Jupiter’s
trojan asteroid.

We also describe the collision manifold and show that it undergoes a
bifurcation as the oblateness coefficient of the asteroid passes through the
zero value. It is interesting to note that when the oblateness approaches zero
(and hence the gravitational potential becomes Newtonian), the limiting
collision manifold that we obtain is not the same as the collision manifold
obtained by applying the McGehee coordinate transformation to the Newtonian
potential. It would be interesting to see if there is a McGehee-type
coordinate transformation for which the limiting collision manifold is the
same as in the Newtonian case. Another interesting problem would be to extend
these results to the spatial Hill four-body problem with oblate bodies.

## References

  * [Bel13] Edward Belbruno.  On the regularizability of the big bang singularity.  Celestial Mechanics and Dynamical Astronomy, 115(1):21–34, 2013\.
  * [BGCG+20] Jaime Burgos-García, Alessandra Celletti, Catalin Gales, Marian Gidea, and Wai-Ting Lam.  Hill Four-Body Problem with Oblate Bodies: An Application to the Sun–Jupiter–Hektor–Skamandrios System.  Journal of Nonlinear Science, pages 1–46, 2020.
  * [BP11] Edward Belbruno and Frans Pretorius.  A dynamical system’s approach to Schwarzschild null geodesics.  Classical and Quantum Gravity, 28(19):195007, 2011.
  * [BX18] Edward Belbruno and BingKan Xue.  Regularization of the big bang singularity with random perturbations.  Classical and Quantum Gravity, 35(6):065013, 2018.
  * [CE71] Charles Conley and Robert Easton.  Isolated invariant sets and isolating blocks.  Transactions of the American Mathematical Society, 158(1):35–61, 1971.
  * [Des15] Pascal Descamps.  Dumb-bell-shaped equilibrium figures for fiducial contact-binary asteroids and EKBOs.  Icarus, 245:64 – 79, 2015.
  * [DMS00] Florin Diacu, Vasile Mioc, and Cristina Stoica.  Phase-space structure and regularization of Manev-type problems.  Nonlinear Analysis: Theory, Methods & Applications, 41(7-8):1029–1055, 2000.
  * [Eas71] Robert Easton.  Regularization of vector fields by surgery.  Journal of Differential Equations, 10(1):92–99, 1971.
  * [ElB09] Mohamed Sami ElBialy.  Collective branch regularization of simultaneous binary collisions in the 3D N-body problem.  Journal of mathematical physics, 50(5):052702, 2009.
  * [GM14] Pablo Galindo and Marc Mars.  McGehee regularization of general SO (3)-invariant potentials and applications to stationary and spherically symmetric spacetimes.  Classical and Quantum Gravity, 31(24):245008, 2014.
  * [McG81] Richard McGehee.  Double collisions for a classical particle system with nongravitational interactions.  Commentarii Mathematici Helvetici, 56(1):524–557, 1981.
  * [MDCR+14] F. Marchis, J. Durech, J. Castillo-Rogez, F. Vachier, M. Cuk, J. Berthier, M.H. Wong, P. Kalas, G. Duchene, M.A. Van Dam, et al.  The puzzling mutual orbit of the binary Trojan asteroid (624) Hektor.  The Astrophysical journal letters, 783(2):L37, 2014.
  * [ORS22] M. Ollé, Ó. Rodríguez, and J. Soler.  Study of the ejection/collision orbits in the spatial RTBP using the McGehee regularization.  Communications in Nonlinear Science and Numerical Simulation, 111:106410, 2022.
  * [Saa74] Donald G Saari.  Regularization and the artificial Earth satellite problem.  Celestial mechanics, 9(1):55–72, 1974.
  * [SSM00] Gheorghe Stoica, Cristina Stoica, and Vasile Mioc.  Branch regularization of quasihomogeneous functions.  Revue Roumaine de Mathematiques Pures et Appliquees, 45(5):897–906, 2000.
  * [XB14] BingKan Xue and Edward Belbruno.  Regularization of the big bang singularity with a time varying equation of state $w>1$.  Classical and Quantum Gravity, 31(16):165002, 2014.