# Microlensing signatures of extended dark objects using machine learning

Miguel Crispim Romao<EMAIL_ADDRESS>Institute for Particle Physics
Phenomenology, Department of Physics, Durham University, Durham DH1 3LE, U.K.
Djuna Croon<EMAIL_ADDRESS>Institute for Particle Physics
Phenomenology, Department of Physics, Durham University, Durham DH1 3LE, U.K.

###### Abstract

This paper presents a machine learning-based method for the detection of the
unique gravitational microlensing signatures of extended dark objects, such as
boson stars, axion miniclusters and subhalos. We adapt MicroLIA, a machine
learning-based package tailored to handle the challenges posed by low-cadence
data in microlensing surveys. Using realistic observational timestamps, our
models are trained on simulated light curves to distinguish between
microlensing by point-like and extended lenses, as well as from other object
classes which give a variable magnitude. We show that boson stars, examples of
objects with a relatively flat mass distribution, can be confidently
identified for $0.8\lesssim r/r_{E}\lesssim 3$. Intriguingly, we also find
that more sharply peaked structures, such as NFW-subhalos, can be distinctly
recognized from point-lenses under regular observation cadence. Our findings
significantly advance the potential of microlensing data in uncovering the
elusive nature of extended dark objects. The code and dataset used are also
provided.

††preprint: IPPP/23/

##  I Introduction

Macroscopic dark matter candidates, with masses ranging from large asteroids
($\sim 10^{-15}\rm M_{\odot}$) to stars ($\sim\rm M_{\odot}$), offer
compelling alternatives to the traditional particle-based theories. These
celestial objects, potentially formed in the early universe, are primarily
detectable via their gravitational effects: gravitational lensing (e.g. Alcock
_et al._ (1998); Niikura  _et al._ (2019a, b); Green and Kavanagh (2021)) and
gravitational waves (e.g. Bird  _et al._ (2016); Ali-Haïmoud  _et al._ (2017);
Kavanagh  _et al._ (2018); Bertone  _et al._ (2020); Abbott  _et al._ (2019);
Chen and Huang (2020); Franciolini (2021); Croon  _et al._ (2023)). Indeed,
gravitational microlensing is one of the most important ways of probing
compact objects such as “machos” or primordial black holes (PBHs), a dark
matter (DM) candidate consisting of compact objects formed in the early
Universe. Through surveys of a range of sources, microlensing of such point-
like lenses has been used to constrain the fraction of DM such objects can
comprise in a wide range of masses (see e.g. Green and Kavanagh (2021)).

It has also been proposed that dark matter can instead be comprised of
extended objects, such as boson stars (e.g. Ruffini and Bonazzola (1969);
Gleiser (1988); Schunck and Mielke (2003); Croon  _et al._ (2019)), axion
miniclusters Fairbairn  _et al._ (2017), and subhalos Erickcek and Sigurdson
(2011); Barenboim and Rasero (2014); Fan  _et al._ (2014); Dror  _et al._
(2018); Graham  _et al._ (2016). Like compact objects, such objects can also
bend the light of distant stars. Whether this effect can be probed by a
microlensing survey depends on the comparison between the object radius and
the Einstein radius – the characteristic length scale which is a function of
the mass of the dark matter lens and the distance to the light source. The
effectiveness of (micro-)lensing then depends on the size of the object
compared to the Einstein radius: dilute dark objects, which are transparent to
light, are ineffective lenses. Using conservative assumptions about the number
of events observed, Refs.Croon  _et al._ (2020a, b) derived modified
constraints on extended dark matter objects.

Interestingly, structures with radii close to the Einstein radius may give
distinct microlensing signatures. How the mass is spread within these
structures affects the lensing effect. This was demonstrated explicitly for
microlensing of various dark matter structures in Croon  _et al._ (2020b); Bai
_et al._ (2020); Croon  _et al._ (2020a). In objects with a flatter density
profile, such as boson stars, the microlensing magnification time series can
deviate significantly from that expected from a point-like lens such as a PBH,
for example featuring caustics crossings as can be seen in Fig. 1. These
distinguishable features of extended dark lenses can in principle be used to
make a positive discovery.

Figure 1: The 10 most distinctive Boson Star light curves, using the dataset
generated with OGLE-II timestamps..

In this work we take the first step towards a positive discovery of a
microlensing signature by an extended dark object. To this end we develop an
analysis pipeline to search for signatures of extended dark matter objects in
time series data. Specifically, we train a histogram-based gradient boosted
classifier to identify such features, and describe how these observations can
be used to search for dark objects in a range of experiments.

This paper is organised as follows. We first review microlensing by extended
lenses and estimate the sensitivity of the OGLE survey to the modified light
curves of boson stars. We then describe the generation of our dataset and
methodology, followed by an analysis of our results and a discussion.

##  II Microlensing signatures of extended lenses

Let us first review some basics of gravitational microlensing, largely
following the treatment in Ref. Narayan and Bartelmann (1996). We will need to
define some important parameters in the lens setup; a depiction of the
geometry can be found in Croon  _et al._ (2020a). The observer-lens, lens-
source, and observer-source distances are denoted $D_{\rm L}$, $D_{\rm S}$,
and $D_{\rm LS}$ respectively. From the perspective of the observer, the lens
center subtends angles of $\beta$ and $\theta_{i}$ with the source and images
of the source respectively. As $D_{\rm L}$, $D_{\rm S}$, and $D_{\rm LS}$ are
much larger than all other scales in the problem, lensing calculations can be
simplified using the small angle approximation. In this approximation, the
deflections $\alpha=4GM/(c^{2}\xi)$ only occur when starlight encounters the
“lens plane” perpendicular to the observer-source axis.

Assuming that the lens is spherically symmetric with density distribution
$\rho(r)$ (and total mass $M=4\pi\int_{0}^{\infty}drr^{2}\rho(r)$), the
lensing equation may be written

$\beta=\theta-\frac{\theta_{\rm E}^{2}}{\theta}\frac{M(\theta)}{M}~{},$ (1)

with the surface mass density projected onto the lens plane as

$\displaystyle M(\theta)$ $\displaystyle=$ $\displaystyle 2\pi D_{\rm
L}^{2}\int_{0}^{\theta}d\theta^{\prime}\theta^{\prime}\Sigma(\theta^{\prime})~{},$
$\displaystyle\Sigma(\theta)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}dz\,\rho\left(\sqrt{D_{\rm
L}^{2}\theta^{2}+z^{2}}\right)~{},$ (2)

and where Einstein (1936)

$\theta_{\rm E}\equiv\sqrt{\frac{4GM}{c^{2}}\frac{D_{\rm LS}}{D_{\rm L}D_{\rm
S}}}~{},$ (3)

is the point-like Einstein angle, the value of $\theta$ for a point-like lens
($M(\theta)\rightarrow M$) at zero impact parameter ($\beta=0$). This in turn
defines the point-like Einstein radius $r_{\rm E}\equiv D_{\rm L}\theta_{\rm
E}$ on the lens plane.

As the only way in which the total lens mass $M$ and the distances $D_{\rm
L},D_{\rm LS},D_{\rm S}$ enter the problem is through their contributions to
the Einstein radius $r_{\rm E}$, it is convenient to express all angles in
units of $\theta_{\rm E}$ and all distances $r_{\rm E}$. Thus, we define
$u\equiv\beta/\theta_{\rm E}=D_{\rm L}\beta/r_{\rm E}$, and
$\tau\equiv\theta/\theta_{\rm E}=D_{\rm L}\theta/r_{\rm E}$ (note that the
latter was named $t$ in Croon  _et al._ (2020b, a)), which allows us to
rewrite (1) as

$u=\tau-\frac{m(\tau)}{\tau},$ (4)

where we have also defined $m(\tau)\equiv M(\theta_{\rm E}\tau)/M$ which
describes the distribution of the lens mass projected onto the lens plane,

$m(\tau)=\frac{\int_{0}^{\tau}d\sigma\sigma\int^{\infty}_{0}d\lambda\,\rho(r_{\rm
E}\sqrt{\sigma^{2}+\lambda^{2}})}{\int_{0}^{\infty}d\gamma\gamma^{2}\rho(r_{\rm
E}\gamma)}~{},$ (5)

where $\rho$ is the density distribution of the lens.

The lensing equation (4) can be used to find the position(s) of the images
$\theta_{i}$ given a position of the lens $\beta$. As the images subtend
different solid angles than the unlensed source, microlensing alters the
observed flux. The magnification is the ratio of the angular extent of an
image to that of the source:

$\displaystyle\begin{split}\mu=&\sum_{i}\left|\frac{\theta_{i}}{\beta}\frac{d\theta_{i}}{d\beta}\right|=\sum_{i}\left|\frac{\tau_{i}}{u}\frac{d\tau_{i}}{du}\right|\\\
=&\sum_{i}\left|1-\frac{m(\tau_{i})}{\tau_{i}^{2}}\right|^{-1}\left|1+\frac{m(\tau_{i})}{\tau_{i}^{2}}-\frac{1}{\tau_{i}}\frac{dm(\tau_{i})}{d\tau_{i}}\right|^{-1}~{}.\end{split}$
(6)

The light curve as a function of time $t$ for a lens with velocity $v$ and
minimum impact parameter $\xi_{\rm min}$ can now be found through $\beta
r_{\rm E}=\sqrt{\xi^{2}_{\rm min}+v^{2}t^{2}}$. Whether a microlensing event
is observable depends on the minimum detectable magnification for a given
telescope, as well as the range of cadences at the microlensing survey, which
sets the transit timescales to which it is sensitive. Typically, a transit is
counted as a lensing event if $\mu>1.34$, which occurs for a point-like lens
for impact parameter (in units of the Einstein radius) $u\equiv\xi/r_{\rm
E}=\beta/\theta_{\rm E}=1$

In Croon  _et al._ (2020a, b) the microlensing efficiency of an extended lens
compared to that of a point-like lens was defined as the maximum impact
parameter $u_{1.34}$ for which a threshold magnification is produced:
$\mu_{\rm tot}(u\leq u_{1.34})\geq 1.34~{}.$ For a point-like lens,
$u_{1.34}=1$ by definition, and the naive expectation for extended lenses
would be that $u_{1.34}\leq 1$. Remarkably, this is not necessarily the case:
in particular for lenses with a reasonably flat density profile, $u_{1.34}$
can be larger than one. Given $u_{1.34}$, the microlensing differential event
rate for a single source with respect to the typical event timescale $t_{E}$
and $x=D_{\rm L}/D_{\rm S}$, observed by a particular experiment, can be
calculated as

$\frac{d^{2}\Gamma}{dxdt_{\rm E}}=\varepsilon(t_{\rm E})\frac{2D_{\rm
S}}{v_{0}^{2}M}f_{\rm DM}\rho_{\rm DM}(x)v_{\rm E}^{4}(x)e^{-v_{\rm
E}^{2}(x)/v^{2}_{0}}~{},$ (7)

where $\varepsilon(t_{\rm E})$ is the efficiency of telescopic detection,
$v_{\rm E}(x)\equiv 2u_{1.34}(x)r_{\rm E}(x)/t_{\rm E}$, $v_{0}=220$ km/s is
the dark matter circular speed in the galaxy, and $\rho_{\rm DM}(x)$ is the DM
density projected onto the line of sight, for example following an isothermal
profile in the Milky Way galaxy,

$\begin{split}\rho_{\rm DM}(r)&=\frac{\rho_{\rm s}}{1+(r/r_{\rm
s})^{2}}~{},\\\ r&\equiv\sqrt{R^{2}_{\rm Sol}-2xR_{\rm Sol}D_{\rm
S}\cos\ell\cos b+x^{2}D_{\rm S}^{2}}~{},\end{split}$

with $R_{\rm Sol}$ = 8.5 kpc, $\rho_{\rm s}=1.39$ GeV/cm3, and $r_{\rm
s}=4.38$ kpc Cirelli  _et al._ (2011); $\ell$ and $b$ are the longitude and
latitude of the source in galactic co-ordinates. The total number of events is
then given by

$N_{\rm events}=N_{\star}T_{\rm obs}\int^{1}_{0}dx\int^{t_{\rm E,max}}_{t_{\rm
E,min}}dt_{\rm E}\frac{d^{2}\Gamma}{dxdt_{\rm E}}~{},$ (8)

where $N_{\star}$ is the number of observed sources, $T_{\rm obs}$ is the
total observation time, and $t_{\rm E,min}$ ($t_{\rm E,max}$) is the minimum
(maximum) timescale of an event.

Comparing $N_{\rm events}$ with the number of observed events in microlensing
surveys, Croon  _et al._ (2020a, b) set constraints on extended lenses.
However, it is important to point out that the surveys identify microlensing
events based on a comparison with the point-like magnification light curve.
For extended lenses, in particular lenses with $\tau_{m}\equiv r_{\rm
lens}/r_{E}\sim 1$, this is not a good approximation. In Fig. 2 we estimate
the range of lens masses and radial sizes for which we expect significant
deviations from the point-like light curve, for the example of a boson star
observed in the OGLE-IV survey (5 year dataset) Niikura  _et al._ (2019b),
anticipating that the non-compactness of the lens can be resolved for
$0.8<\tau_{m}<3$ (which we will verify in the next section). Here we follow
the mass profile of boson stars $m(t)$ outlined in the appendix of Croon  _et
al._ (2020a). As in this work, the sensitivity is computed using the
Poissonian $90\%$ confidence limit based on a signal comprised of dark matter
and astrophysical foregrounds (see Niikura  _et al._ (2019b)). We note that as
the OGLE collaboration did not search for the particular microlensing light
curves predicted by boson stars, this is an estimate only. Comparing Fig. 2
with the sensitivity of the OGLE-IV survey to boson stars in Croon  _et al._
(2020a), we note that for a given lens size, the lighter masses lead to
distinguishable features in the light curve, as can be expected from the
dependence of the Einstein radius on $M$.

Figure 2: Sensitivity of the OGLE survey to a modified lensing light curve of
a boson star. Here we anticipate that a boson star can be distinguished from a
point-like lens for $0.8<\tau_{m}<3$.

Because the Einstein radius varies along the line of sight to the source, for
a lens of a particular size $r_{\rm lens}$, a range of $\tau_{m}(x)$ is
relevant, where $x=D_{\rm L}/D_{\rm S}$. One might wonder what the
distribution in $\tau_{m}$ is. This is given by ${\rm d}N/{\rm
d}\tau_{m}=({\rm d}N/{\rm d}x)({\rm d}x/{\rm d}\tau_{m})=f_{\rm DM}\rho_{\rm
DM}(x)\;{\rm d}x/{\rm d}\tau_{m}$. For the OGLE sources, we find that this
distribution peaks for $\tau_{m}(x=0.5)$, and rapidly falls like
$\tau_{m}^{-3}$ away from it. Thus, we expect that given a particular lens
mass, it is in practise possible to interpret a measurement of $\tau_{m}$
directly in terms of a lens size.

##  III Dataset Generation and Methodology

In this work we adapt MicroLIA, a tool developed and detailed in Godines  _et
al._ (2019). The classifier developed in this paper is a machine learning
model utilising the Random Forest algorithm, specifically designed for the
detection of microlensing events in astronomical surveys. It is tailored to
handle the challenges posed by low-cadence data, which typically suffer from
irregular signal sampling and thus lower signal-to-noise ratios, making
microlensing event detection more difficult. MicroLIA distinguishes between
microlensing and other variable star events using 148 features derived from
the light curve and its derivative time series. The classifier categorises
events into classes like microlensing, eclipsing binaries, and regular
variable stars, focussing on accurate identification of microlensing amidst
these.

In this work, we extend the scope of MicroLIA by including extended dark
matter objects which act like non-point-like microlensing lenses, such as
Boson Stars (BS) and Navarro-Frenk-White subhalos (NFW). For this, we
generated microlensing light curves for these extended objects using their
mass profile, $m(\tau)$, which have previously been calculated in Croon  _et
al._ (2020a). From these mass profiles, we fit an interpolating function, and
for each value of the impact parameter, $u$, we solve the microlensing
equation 4 to obtain the total magnification to the light curve, 6.

The observed impact parameters are closely related to the cadence of a survey.
For a lens with a characteristic timescale defined by the crossing of the
Einstein radius $t_{E}{\equiv R_{E}/v}$ and a minimum impact parameter $u_{0}$
producing a magnification peak occurring at $t_{0}$, the values of the
observed impact parameters are

$u=\sqrt{u_{0}^{2}+\left(\frac{t-t_{0}}{t_{E}}\right)^{2}}\ ,$ (9)

where $t$ is the survey time. Therefore, the cadence at which the survey
collects observational data will have a significant impact in the observed
light curve. For this study, we considered two possible cases for the data
collection timestamps. In the first case we used OGLE-II timestamps (but not
the light curves nor their errors), which are input into MicroLIA and randomly
sampled when simulating individual light curves. We note that the regularity
of the cadence in OGLE-II is not significantly different from other iterations
of the survey. In the second case we considered a _perfectly regular daily
cadence_ , where timestamps are all equally spaced by the same interval. As we
will see below, this case demonstrates some identification opportunities that
are obscured by irregular cadences, such as the case with OGLE-II timestamps.
We will refer to these to cases as OGLE-II Timestamps and Regular Daily
Cadence, respectively.

Both datasets were generated by simulating 100 000 light curves for each of
the six classes: Cataclysmic Variables (CV), RR Lyrae and Cepheid Variables
(VARIABLE), Mira long-period variables (LPV), Point-like Microlensing (ML),
Boson Stars (BS), and NFW Subhalos (NFW). The CV, VARIABLE, LPV, and ML light
curves were generated using MicroLIA’s simulation, and BS and NFW were
generated using our own simulation. We applied the same selection criteria for
the three microlensing source events (ML, BS, NFW). The criteria are the same
as MicroLIA’s, but we further demanded that the observed magnification be at
least $1.34$, a common criterion imposed by microlensing surveys. The light
curves were simulated with a minimum magnitude of 15, a maximum magnitude of
20, and with Gaussian noise.

The extended microlensing sources, BS and NFW, have mass profiles which depend
on the parameter $\tau_{m}$, which follows some distribution depending on the
prevalence of such objects. For this study, we sampled
$\tau_{m}\sim\mathcal{U}(0.5,5)$ logarithmically, since, as discussed in the
previous section, this distribution is strongly peaked at the $\tau_{m}$
reached at $x=0.5$, and the largest deviations in the light curves are
expected for $\tau_{m}\sim 1$. Finally, we sampled the minimal impact
parameter, $u_{0}$, differently for each of the three microlensing sources

$u_{0}\sim\begin{cases}\mathcal{U}(0,1)&\text{for ML}\\\
\mathcal{U}(0,1.5)&\text{for BS}\\\ \mathcal{U}(0,1.1)&\text{for NFW}\
,\end{cases}$ (10)

where we emphasise that these are used during the simulation step, i.e. before
the selection criteria, including $\mu\geq 1.34$, are imposed.

In total, we generated 600 000 curves, of which only a few
$\mathcal{O}(0.1\%)$ with missing values for some of the features were
dropped. The dataset was then split into train, validation, and test subsets
with proportions 0.5:0.25:0.25, respectively corresponding to around 300 000,
150 000, and 150 000 light curves, with the six classes being equally
represented in each of the sets. For each light curve, we computed 74 features
using the light curve time series, in addition to the same 74 using the
derivative of the light curve time series, to a total of 148 features. The
computed features relate to statistics of the time series of the light curve,
as well as other time series quantities. See Godines  _et al._ (2019) and
MicroLIA’s API reference for the full description and reference of each
feature.111https://microlia.readthedocs.io/en/latest/autoapi/MicroLIA/features/
The training set was used to conduct exploratory data analysis and to train
machine learning models. The validation set was used for model selection and
comparison, and to produce preliminary analysis plots and statistics. The test
set was used to produce the final analysis plots and statistics, presented in
the next section. The full dataset can be downloaded here Crispim Romão and
Croon (2024).

For the multiclassification task, we trained a histogram-based gradient
boosted classifier, which is a member of the broader machine learning
algorithm family of gradient boosting machines (GBM). Like Random Forests
(RF), implemented in MicroLIA, GBM are ensemble models that leverage the power
of multiple weaker estimators, usually small trees, to produce a strong
estimator. In RF the trees are independent of each other and the final
prediction is obtained by the ensemble average over all the predictions from
the trees in the forest. In GBM the trees are not independent but sequentially
trained as to improve on the previous iteration. Schematically, consider
$F(X)$ to be the output of the GBM as a function of the data, $X$. In the
first iteration, we want to train a simpler estimator, $F_{1}(X)$, to match a
target label, $y$. Since the estimator is simple, it will not be very accurate
and the prediction will have a certain residual error, $y-F_{1}(X)$. In the
second iteration, we train another weak estimator, $h_{1}(X)$ not on the
desired output, but on residual error of the previous step to create a new,
improved, estimator $F_{2}(X)=F_{1}(X)+h_{1}(X)$, since in the ideal case
$h_{1}(X)=y-F_{1}(X)$ and we would have the desired prediction. However, each
step will have a residual error $y-F_{m}(X)$, and so the process can be
repeated until the desired accuracy (or maximum number of iterations) is
reached. The final GBM output will then be a function of all the steps,
$F(X)=\sum_{m}F_{m}(X)+h_{m}(X)$. This presentation is schematic and can be
generalised to any problem with a differentiable loss function.222The term
_boosting_ arises from iteratively boosting the performance using the output
of another weak learner. The term _gradient_ comes from the observation that
fitting successive weak estimators on the residuals of the previous iteration
is equivalent to a gradient descent step in a function space spanned by the
weak learners. GBM have been known in the literature to be powerful estimators
for tabular data, which we explore in this work. The usual implementation of
GBM uses simple trees at each step. However, this requires sorting the data at
each iteration, making tree-based GBM computationally heavy for datasets
larger than a few tens of thousands examples, which is our case. To mitigate
this, a histogram-base variant has been developed that sorts and binarises
(i.e., assign each data point to the bins of the histogram) the data once,
producing orders of magnitude speed improvements in training and prediction.
In this work, we used the histogram-based GBM implemented by `scikit-learn`’s
`HistGradientBoostingClassifier` (HGBC).

In a preliminary study, we compared HGBC against other GBM implementations and
`scikit-learn`’s RF, observing all GBM to outperform the RF when computing the
Area Under the Curve (AUC) of the Receiver Operating Characteristic (ROC)
evaluated on the validation set. We then tuned the hyperparameters of the
HGBC, observing minimal improvements of its discriminating performance, for
which reason we decided to keep the default hyperparamters for the rest of the
study presented herein.333The performance improvements only impacted at most
the third significant digit of the validation ROC AUC, which is arguably at
the level of statistical fluctiation of the dataset itself, and therefore
meaningless.

##  IV Analysis

In this section, we present the multiclassification analysis on the two
generated sets: the first generated with the OGLE-II timestamps, and the
second generated with ideal regular daily cadence. The purpose of performing
two analyses is to assess whether and how the discrimination power is affected
by the cadence.

###  IV.1 OGLE-II Timestamps

In fig. 3 we present the confusion matrix of the six-way (All vs. All)
multiclassification performed by the HGBC. We observe that the non-
microlensing events – CV, LPV, VARIABLE – have minimal to nonexisting overlap
with any class. Conversely, the HGBC prediction for the microlensing events –
ML, BS, and NFW – has significant overlap. However, we observe that BS suffer
from considerable less contamination from ML and NFW events than these two do
of each other. This suggests that BS events are the easiest to isolate, and
therefore to detect, of the three microlensing cases.

Figure 3: Confusion matrix for the six-way All vs. All multiclassification
performed by the HGBC using the dataset generated with OGLE-II timestamps. The
entries are rounded to three significant digits.

Since fig. 3 suggests that BS events are easy to isolate, it is important to
understand the nature of the BS we identify. All BS follow the same mass
profile, $m(t)$, which can be wider or narrower, depending on the parameter
$\tau_{m}$, associated with the radius of the BS. In fig. 4 we show how the
probability of a BS being identified as such by the HGBC depends on
$\tau_{m}$. We find that there is a _sweet spot_ at $\tau_{m}\simeq 2$ to
maximise the correct identification of BS, with some some BS light curves
being classified with high confidence, i.e. with $P(y=BS|X)\simeq 1$. This
also motivates a posterior to our choice of interval $0.8<\tau_{m}<3$ in fig.
2. In fig. 1 we present the 10 BS light curves with the highest value of
$P(y=BS|X)$, and we observe that they all exhibit the three peak magnification
profile produced by caustics that one expects from BS sources.

Figure 4: Probability of a BS light curve being correctly identified as one by
the HGBC versus the Boson Star $\tau_{m}$ parameter, using the datset
generated with OGLE-II timestamps.

Our analysis so far has shown that the discrimination between microlensing
sources (ML, BS, NFW) and other sources of light curves (CV, LPV, VARIABLE) is
a relatively easy task when OGLE-II cadence is observed, with almost
nonexisting overlap between these two broad classes. As such, we now focus on
the difficult three-way classification focused on discriminating between the
three microlensing sources. In fig. 5 we present the ROC and the AUC of each
class (versus the other two) for a HGBC trained on this three-way
multiclassification problem.444The AUC of a ROC curve can then be seen as an
average performance over all possible thresholds set on the predicted
probability used to predict the class. The advantage of analysing a ROC curve,
and its area, over a classification matrix is that the ROC curve captures the
classification performance over all possible output threshold cuts, while the
confusion matrix only shows the true and false positives of the assigned
predicted class as the one that has maximal probability. For example, in fig.
5 we can see that BS light curves can be isolated with high purity (i.e. with
False Positive Rate around $0.0$), while that is not possible for neither ML
nor NFW light curves, again reinforcing that BS light curves are easier to
identify than the remaining microlensing sources due to their unique three-
peak profile. Furthermore, fig. 5 suggests that point-like microlensing light
curves are only slightly easier to identify than NFW ones, since their ROC
curve has a True Positive Rate greater than that of the NFW for all values of
the False Positive Rates. This is also evidenced by the ML AUC, $0.65$, which
is greater than that of the NFW, $0.62$.

Figure 5: ROC curves and their areas obtained from a HGBC trained on the
three-way multiclassification task using the dataset generated with OGLE-II
timestamps.

Having been able to isolate the BS light curves from other microlensing
sources, we now turn to interpreting our results. Unfortunately, although HBGC
are incredibly powerful estimators they do no provide a clear way of
interpreting their prediction, a common challenge when using machine learning
in highly dimensional multivariate analyses. In this study, we implemented a
backward Sequential Feature Selection (SFS) loop to assess which of the 148
features are relevant for this classification task. The backward SFS loop
starts by fitting an HGBC using all 148 features and evaluating its
performance on the validation set. Then we train 148 HGBC on all the possible
148 subsets comprised of 147 features (i.e., with one feature removed) and
evaluate their performance on the validation set. Keeping only the best set of
147 features, the same step can be done for the 146 subsets and so one, for a
total of 148 iterations totaling 11026 HGBC trained on 11026 subset of
features, which is a feasible study given the increased training speed offered
by the HGBC.555Noting that a set of 148 elements has $2^{148}-1$ non-empty
subsets, 11026 is far fewer than the total possible combinations of features.
In fig. 6 we show the performance of the HGBC on the three-way microlensing
focused multiclassification for each of the possible sources and their
geometric mean. We see that only around 25 features derived from the light
curve time series and its derivatives, out of the 148, are relevant for this
task.666It is important to point out that while this methodology is aimed at
removing uninformative features, it will also remove redundant, i.e. highly
correlated, features. Therefore, the final set of features is only unique up
to fluctuations of the data and the stochastic initialisation of the HGBC. To
mitigate this, we fixed the seed of the run, but modifications of the training
and validation set could still hold a different set of features. We also note
that the discrimination performance degrades for all the three classes at the
same stage, suggesting that the final 25 features are relevant for all three
classes. The list of the final 25 features by survival ranking is presented in
table 1.

Figure 6: Backward SFS (features reduced from right to left) for a HGBC
trained on the three-way multiclassification task using the dataset generated
with OGLE-II timestamps.

For the dataset using OGLE-II timestamps, which is the subject of the analysis
in this section, we find that the _complexity_ of the time series is the most
important feature to separate the light curves from the three microlensing
sources, surviving the whole backward selection loop until there is only one.
In fig. 7 we show the distribution of the time series complexity for BS
against the $\tau_{m}$, where we witness a boost in the time series complexity
for $\tau_{m}\simeq 2$, the same sweet spot identified above in fig. 4,
suggesting that the complexity of the time series of a BS light curve is a
driving feature for it to be correctly identified as a BS. The same is not
observed for the NFW events, which we do not show for the sake of presentation
tidiness.

Figure 7: Complexity vs $\tau_{m}$ scatter for BS events with OGLE-II
timestamps.

In fig. 1 we show the 10 most distinctive BS light curves, i.e. the 10 BS
light curves with highest $P(y=BS|X)$ as identified by the HGBC trained on the
three-way multiclassification task on the OGLE-II dataset. These light curves
all exhibit very clear three magnitude peaks, a hallmark feature arising from
the caustics produced by the extended nature of the lens.

###  IV.2 Regular Daily Cadence Timestamps

So far our analysis has focused on light curves simulated using OGLE-II
timestamps. This reflects well the sensitivity to extended objects in a
realistic microlensing survey, with its given observational constraints.
However, one might wonder to what extent the conclusions previously drawn are
sensitive to the observation cadence details, especially its irregularity. To
address this, we conduct the analysis with light curves simulated with regular
daily cadence, i.e. where observations are taken exactly 24 hours apart.

We begin with the six-way All vs All multiclassification task. In fig. 8 we
show the confusion matrix obtained using the HGBC. Although similar to its
counterpart with OGLE-II timestamps, fig. 3, it has some noticeable
differences. First, we notice that the contamination in BS positive
predictions has decreased from $\mathcal{O}(10\%)$ down to $\mathcal{O}(1\%)$,
suggesting an improved capacity of the HGBC in producing a pure sample of BS
light curves. Secondly, we notice that the cross misclassification between NFW
and ML light curves has also decreased, improving upon the OGLE-II timestamps
case.

Figure 8: Confusion matrix for the six-way All vs. All multiclassification
performed by the HGBC using the datset generated with regular daily cadence
timestamps. The entries are rounded to three significant digits.

In fig. 9 we revisit the probability of a BS light curve being classified as
such versus the BS $\tau_{m}$, but this time for the regular cadence case. We
see that the previously identified  _sweet spot_ at $\tau_{m}\simeq 2$ has now
almost no missclassifications into other classes, showing how a regular
cadence can lead to an even better identification of BS light curves; this is
due to a better resolution of the light curve being able to capture the three
magnification peaks.

Figure 9: Probability of a BS light curve being correctly identified as one by
the HGBC versus the Boson Star $\tau_{m}$ parameter, using the datset
generated with regular cadence timestamps.

As before, we now focus on the three-way multiclassification problem focused
on the three microlensing sources. The ROC curves for each of the three
classes and their respective AUC can be seen in fig. 10. Again, we see that
there are some noticeable differences to the OGLE-II case shown in fig. 5. The
first thing to notice is that all AUC have increased compared to the OGLE-II
timestamps case, indicating an easier discrimination when using regular daily
cadence timestamps. Next, we see how the ROC curve for the BS is significantly
more peaked for small False Positive Rate, providing further evidence that the
BS light curves are better classified with regular daily cadence. More
interestingly, however, is how the ML and NFW ROC curves no longer follow the
same trend as in the OGLE-II case. More precisely, we can observe that they
cross, as before they did not. This happens halfway through the curve, with
the NFW ROC curve having the _upper hand_ over the ML ROC curve for lower
values of the False Positive Rate, with an added feature that the ROC curve
for the NFW can have non-vanishing True Positive Rate at low False Positive
Rate values. This implies that with regular daily cadence it is possible to
isolate NFW light curves with little to no contamination of any other class,
something that was impossible to achieve using the OGLE-II cadence.

Figure 10: ROC curves and their areas obtained from a HGBC trained on the
three-way multiclassification task using the datset generated with regular
daily cadence timestamps.

The previous result points at the possibility in principle of completely
isolating NFW light curves. It is then important to understand the nature of
the NFW light curves that we can isolate. In fig. 11 we show how the
probability for an NFW light curve to be correctly identified as such can vary
with the NFW $\tau_{m}$ parameter. Although less noticeable than what we saw
before for the BS light curves, we can observe a _sweet spot_ for correct
classification at $1\lesssim\tau_{m}\lesssim 2$, with some light curves being
assigned $P(y=NFW|X)\simeq 1$. We notice that the equivalent scatter for the
OGLE-II timestamps, not shown here to declutter the presentation, does not
exhibit this pattern, suggesting that there is important information in the
light curve that can only be obtained with regular cadence.

Figure 11: Probability of a NFW light curve being correctly identified as one
by the HGBC versus the NFW $\tau_{m}$ parameter, using the dataset generated
with regular daily cadence timestamps.

Contrary to BS light curves, NFW light curves do not have a clear visual
profile compared to point-like microlensing sources. In fig. 12 we compare the
100 most easily identifiable point-like microlensing curves, i.e. those with
the highest $P(y=ML|X)$, against the 100 most easily identifiable NFW curves,
i.e. those with the highest $P(y=NFW|X)$, where the probabilities are obtained
from the HGBC trained on the three-way multiclass classification on the
regular cadence timestamps dataset. We normalised the magnitudes by min-maxing
their values to fall under the $[0,1]$ interval to aid visual comparison.
Although the difference is very nuanced, the NFW curves tend to be _narrower_
than the ML curves. This subtlety explains why regular cadence is so important
for identifying NFW light curves, as one needs a better _resolution_ of the
light curve profile to be able to identify this nuance, which proved
impossible when using the OGLE-II timestamps.

Figure 12: Comparison between the 100 most easily identifiable ML curves and
the 100 most easily identifiable NFW curves.

Finally, we perform a backward SFS for the regular daily cadence dataset using
a three-way HGBC focused on the microlensing classes. The discrimination
performance versus the number of features is shown in fig. 13, where we
observe a similar trend than that already discussed for the OGLE-II timestamps
dataset with the ROC AUC degrading significantly below 25 features. However,
it is worth noting that the performance degrades first for the NFW and ML
cases, with the ROC AUC associated with the classification of the BS light
curves staying at its maximum value until around 15. This suggests that there
are less relevant features to distinguish BS light curves from the other
classes than there are to correctly identify the others. In table 1 we present
the top 25 features, and we observe considerable overlap with the OGLE-II
study. Perhaps curiously, we count 10 of the features to have been computed on
the derivative of the time series, whereas for the OGLE-II the number is six,
possibly hinting at the importance of the derivative of the time series in
identifying the nuanced differences in shape between the ML and the NFW light
curve magnification peak.

Figure 13: Backward SFS (features reduced from right to left) for a HGBC
trained on the three-way multiclassification task using the dataset generated
with regular daily cadence timestamps.

##  V Discussion

In this work we have studied the observation of microlensing signatures due to
extended dark objects in time series data, and its distinction from other
signals, most importantly point-like lenses. We have focused on two lens
profiles, each exemplifying a class of objects: the profile of an NFW-subhalo
is more peaked, whereas boson stars have a more diffuse profile. As expected,
the boson stars can be confidently distinguished from point-like lenses, for
$0.8\lesssim\tau_{m}(=r_{\rm lens}/r_{E})\lesssim 3$. This is owed to the
characteristic caustics in the light curves for these objects.

As an exercise, we studied how the regularity of the time series cadence
affects the confidence of the detections. Interestingly, for regular daily
cadence we also find confident detections of NFW-subhalos. These detections
occur for $0.9\lesssim\tau_{m}\lesssim 4$, despite the fact that no caustics
are observed in the light curves. Though the regularity of (daily)
observations is dependent on many conditions that are beyond the observer’s
control, this is an interesting observation which warrants investigation for
other microlensing scenarios and targets.

Stellar binaries or exoplanets orbiting a lensing body may give rise to
perturbations and caustics in the light curve. Unlike for extended dark
objects, these are asymmetric or one-sided features (see e.g. Khakpash  _et
al._ (2021)). Microlensing has allowed for the discovery of over 100
exoplanets, particularly of near-terrestial size Exo . In this work, we have
not considered these signatures, which may give rise to confusion with boson
star identifications with low significance. We will leave such an analysis for
future work.

We performed the analysis in this work by extending the MicroLIA algorithm,
which utilises 148 features derived from the light curves in a single band to
distinguish between objects. Our SFS analysis showed that of these, only 25
were needed for the optimal classification. For our regular cadence analysis,
further features could help distinguish between NFW-subhalos and point-like
lenses. In future work, we will study the potential benefits of learning
directly on the light curves.

We note that in our analysis we have assumed a Gaussian noise model, which may
not be realistic. In future work towards the application of our methodology,
further noise models should be considered. These could lead to larger
misclassification between classes of events.

Microlensing surveys typically only release data on candidate events, which
have passed through a strict selection in which our extended dark object light
curves were likely cut. Exceptions include the UKIRT Microlensing Survey UKI ;
Shvartzvald  _et al._ (2017) and the VISTA Variables in the Via Lactea Survey
Minniti  _et al._ (2010), which we plan to analyse in future work. Upcoming
microlensing opportunities include the launch of infrared astronomy
experiments in the mid-2020s. In particular, the Nancy Grace Roman Space
Telescope (previously WFIRST; a space mission by NASA) has as a key objective
to discover exoplanets through the microlensing technique Spergel  _et al._
(2015), but will also probe primordial black holes Fardeen  _et al._ (2023)
and extended dark objects.

Finally, for some surveys, the finite extend of the source becomes important.
This is the case, for example, for the Subaru-hsc survey of M31 Niikura  _et
al._ (2019c): because of its sensitivity to small transit times (and hence
small Einstein radii), the angular extent of source stars corresponds to a
distance at the lens larger than the Einstein radius, suppressing the
magnification relative to point-like sources Witt and Mao (1994); Montero-
Camacho  _et al._ (2019); Smyth  _et al._ (2020); Croon  _et al._ (2019). We
leave an analysis of the magnification curves with finite source effects for
future work.

## Software

For the computation of the mass profile of the extended sources and the
sensitivity estimation, we used ` Mathematica ` version 12.

The dataset generation, machine learning training, and the final analysis were
performed in `python 3.9.18`, making use of several packages, of which we make
especial note:

  * •

Light curve simulation and time series feature extraction was performed using
MicroLIA. We used the code directly from the MicroLIA github repository, as it
includes considerable changes to the code provided with the packaged version
`2.6.0`. More precisely, the version used makes proper use of the derivative
time series, as well as its errors as statistical weights. These two fixes
were contributions of this work.

  * •

`scikit-learn 1.3.2` for HGBC implementation.

  * •

`mlxtend 0.23.0` for the backwards SFS.

  * •

`SciencePlots 2.1.1` for plotting.

The full list of packages and their versions can be found in the
`requirements.txt` file included in the code repository hosting the code used
in this work.777https://gitlab.com/miguel.romao/microlensing-extended-objects-
machine-learning

## Acknowledgements

We thank Daniel Godines for useful discussions and help regarding MicroLIA.
MCR and DC are supported by the STFC under Grant No. ST/T001011/1.

##  Appendix A Backward Sequential Feature Selection Results

In table 1 we collect the top 25 features obtained by the backward Sequential
Feature Selection (SFS), for both the OGLE-II and the regular daily cadence
timestamps datasets. The details of the features can be seen in Godines  _et
al._ (2019) and in MicrioLIA’s code. The feature rank is assocaited with how
long it survived the SFS loop. E.g., a rank $1$ means that the feature is
present in all subsets down to the last subset of a single feature, whereas a
rank $n$ means that the feature is present up until the subset of size $n$,
but is then removed and is not present in smaller subsets. This set of
features is only unique up to fluctuations of the data and the stochastic
variability associated to HGBC training. In principle, any feature presented
in this table could be replaced by a highly correlated counterpart, and as
such the list presented herein should be interpreted as the maximal set of
mutually uninformative features obtained using our training and validation
sets.

Feature Rank | OGLE-II | Daily Regular Cadence
---|---|---
1 | complexity | stetsonJ (derivative)
2 | median buffer range | half mag amplitude ratio
3 | ratio recurring points | FluxPercentileRatioMid80
4 | amplitude | median buffer range
5 | medianAbsDev | count below
6 | FluxPercentileRatioMid80 | medianAbsDev
7 | count below | check max last loc
8 | cusum | longest strike below
9 | shapiro wilk | complexity
10 | mean change | PercentDifferenceFluxPercentile
11 | longest strike above (derivative) | longest strike above
12 | mean second derivative | shapiro wilk (derivative)
13 | number of crossings | FluxPercentileRatioMid65
14 | LinearTrend | number of crossings
15 | longest strike below (derivative) | number cwt peaks
16 | sample entropy | number cwt peaks (derivative)
17 | FluxPercentileRatioMid65 | time reversal asymmetry (derivative)
18 | mean change(derivative) | FluxPercentileRatioMid35
19 | longest strike above | peak detection(derivative)
20 | mean n abs max (derivative) | mean n abs max (derivative)
21 | FluxPercentileRatioMid20 | longest strike above (derivative)
22 | quantile (derivative) | check min last loc
23 | check min last loc | quantile (derivative)
24 | Gskew | sample entropy (derivative)
25 | check max last loc (derivative) | below3 (derivative)
Table 1: Top 25 features, ranked by _survivability_ of the backward Sequential
Feature Selection loop. Features computed in the derivative of the time series
are identified as such. See Godines  _et al._ (2019) for details about the
features.

## References

  * Alcock  _et al._ (1998) C. Alcock, R. A. Allsman, D. Alves, R. Ansari, É. Aubourg, T. S. Axelrod, P. Bareyre, J.-P. Beaulieu, A. C. Becker, D. P. Bennett, S. Brehin, F. Cavalier, S. Char, K. H. Cook, R. Ferlet, J. Fernandez, K. C. Freeman, K. Griest, P. Grison, M. Gros, C. Gry, J. Guibert, M. Lachièze-Rey, B. Laurent, M. J. Lehner, É. Lesquoy, C. Magneville, S. L. Marshall, É. Maurice, A. Milsztajn, D. Minniti, M. Moniez, O. Moreau, L. Moscoso, N. Palanque-Delabrouille, B. A. Peterson, M. R. Pratt, L. Prévôt, F. Queinnec, P. J. Quinn, C. Renault, J. Rich, M. Spiro, C. W. Stubbs, W. Sutherland, A. Tomaney, T. Vandehei, A. Vidal-Madjar, L. Vigroux, and S. Zylberajch, The Astrophysical Journal 499, L9 (1998).
  * Niikura  _et al._ (2019a) H. Niikura _et al._ , Nature Astron. 3, 524 (2019a), arXiv:1701.02151 [astro-ph.CO] .
  * Niikura  _et al._ (2019b) H. Niikura, M. Takada, S. Yokoyama, T. Sumi, and S. Masaki, Phys. Rev. D 99, 083503 (2019b), arXiv:1901.07120 [astro-ph.CO] .
  * Green and Kavanagh (2021) A. M. Green and B. J. Kavanagh, J. Phys. G 48, 043001 (2021), arXiv:2007.10722 [astro-ph.CO] .
  * Bird  _et al._ (2016) S. Bird, I. Cholis, J. B. Muñoz, Y. Ali-Haïmoud, M. Kamionkowski, E. D. Kovetz, A. Raccanelli, and A. G. Riess, Phys. Rev. Lett. 116, 201301 (2016), arXiv:1603.00464 [astro-ph.CO] .
  * Ali-Haïmoud  _et al._ (2017) Y. Ali-Haïmoud, E. D. Kovetz, and M. Kamionkowski, Phys. Rev. D 96, 123523 (2017), arXiv:1709.06576 [astro-ph.CO] .
  * Kavanagh  _et al._ (2018) B. J. Kavanagh, D. Gaggero, and G. Bertone, Phys. Rev. D 98, 023536 (2018), arXiv:1805.09034 [astro-ph.CO] .
  * Bertone  _et al._ (2020) G. Bertone _et al._ , SciPost Phys. Core 3, 007 (2020), arXiv:1907.10610 [astro-ph.CO] .
  * Abbott  _et al._ (2019) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 123, 161102 (2019), arXiv:1904.08976 [astro-ph.CO] .
  * Chen and Huang (2020) Z.-C. Chen and Q.-G. Huang, JCAP 08, 039 (2020), arXiv:1904.02396 [astro-ph.CO] .
  * Franciolini (2021) G. Franciolini,  _Primordial Black Holes: from Theory to Gravitational Wave Observations_ , Ph.D. thesis, Geneva U., Dept. Theor. Phys. (2021), arXiv:2110.06815 [astro-ph.CO] .
  * Croon  _et al._ (2023) D. Croon, S. Ipek, and D. McKeen, Phys. Rev. D 107, 063012 (2023), arXiv:2205.15396 [astro-ph.CO] .
  * Ruffini and Bonazzola (1969) R. Ruffini and S. Bonazzola, Phys. Rev. 187, 1767 (1969).
  * Gleiser (1988) M. Gleiser, Phys. Rev. D38, 2376 (1988), [Erratum: Phys. Rev.D39,no.4,1257(1989)].
  * Schunck and Mielke (2003) F. E. Schunck and E. W. Mielke, Classical and Quantum Gravity 20, R301–R356 (2003).
  * Croon  _et al._ (2019) D. Croon, J. Fan, and C. Sun, JCAP 1904, 008 (2019), arXiv:1810.01420 [hep-ph] .
  * Fairbairn  _et al._ (2017) M. Fairbairn, D. J. E. Marsh, and J. Quevillon, Phys. Rev. Lett. 119, 021101 (2017), arXiv:1701.04787 [astro-ph.CO] .
  * Erickcek and Sigurdson (2011) A. L. Erickcek and K. Sigurdson, Phys. Rev. D 84, 083503 (2011).
  * Barenboim and Rasero (2014) G. Barenboim and J. Rasero, JHEP 04, 138 (2014), arXiv:1311.4034 [hep-ph] .
  * Fan  _et al._ (2014) J. Fan, O. Özsoy, and S. Watson, Phys. Rev. D 90, 043536 (2014).
  * Dror  _et al._ (2018) J. A. Dror, E. Kuflik, B. Melcher, and S. Watson, Phys. Rev. D 97, 063524 (2018).
  * Graham  _et al._ (2016) P. W. Graham, J. Mardon, and S. Rajendran, Phys. Rev. D 93, 103520 (2016).
  * Croon  _et al._ (2020a) D. Croon, D. McKeen, and N. Raj, Phys. Rev. D 101, 083013 (2020a), arXiv:2002.08962 [astro-ph.CO] .
  * Croon  _et al._ (2020b) D. Croon, D. McKeen, N. Raj, and Z. Wang, Phys. Rev. D 102, 083021 (2020b), arXiv:2007.12697 [astro-ph.CO] .
  * Bai  _et al._ (2020) Y. Bai, A. J. Long, and S. Lu, JCAP 09, 044 (2020), arXiv:2003.13182 [astro-ph.CO] .
  * Narayan and Bartelmann (1996) R. Narayan and M. Bartelmann, in  _13th Jerusalem Winter School in Theoretical Physics: Formation of Structure in the Universe Jerusalem, Israel, 27 December 1995 - 5 January 1996_ (1996) arXiv:astro-ph/9606001 [astro-ph] .
  * Einstein (1936) A. Einstein, Science 84, 506 (1936).
  * Cirelli  _et al._ (2011) M. Cirelli, G. Corcella, A. Hektor, G. Hütsi, M. Kadastik, P. Panci, M. Raidal, F. Sala, and A. Strumia, Journal of Cosmology and Astroparticle Physics 2011, 051–051 (2011).
  * Godines  _et al._ (2019) D. Godines, E. Bachelet, G. Narayan, and R. Street, Astronomy and Computing 28, 100298 (2019).
  * Crispim Romão and Croon (2024) M. Crispim Romão and D. Croon, “Light curves for variable, point-like microlensing, and extended objects microlensing sources with regular cadence and OGLE-II timestamps cadence.”  (2024).
  * Khakpash  _et al._ (2021) S. Khakpash, J. Pepper, M. Penny, B. S. Gaudi, and R. Street, The Astronomical Journal 161, 132 (2021).
  * (32) “Nasa exoplanet archive,” https://exoplanetarchive.ipac.caltech.edu/, accessed: 25-01-2024.
  * (33) “Ukirt microlensing survey information,” https://exoplanetarchive.ipac.caltech.edu/docs/UKIRTMission.html, accessed: 25-01-2024.
  * Shvartzvald  _et al._ (2017) Y. Shvartzvald, G. Bryden, A. Gould, C. Henderson, S. Howell, and C. Beichman, The Astronomical Journal 153, 61 (2017).
  * Minniti  _et al._ (2010) D. Minniti, P. Lucas, J. Emerson, R. Saito, M. Hempel, P. Pietrukowicz, A. Ahumada, M. Alonso, J. Alonso-Garcia, J. I. Arias, _et al._ , New Astronomy 15, 433 (2010).
  * Spergel  _et al._ (2015) D. Spergel _et al._ , (2015), arXiv:1503.03757 [astro-ph.IM] .
  * Fardeen  _et al._ (2023) J. Fardeen, P. McGill, S. E. Perkins, W. A. Dawson, N. S. Abrams, J. R. Lu, M.-F. Ho, and S. Bird, (2023), arXiv:2312.13249 [astro-ph.GA] .
  * Niikura  _et al._ (2019c) H. Niikura _et al._ , Nat. Astron. 3, 524 (2019c), arXiv:1701.02151 [astro-ph.CO] .
  * Witt and Mao (1994) H. J. Witt and S. Mao, ApJ 430, 505 (1994).
  * Montero-Camacho  _et al._ (2019) P. Montero-Camacho, X. Fang, G. Vasquez, M. Silva, and C. M. Hirata, JCAP 08, 031 (2019), arXiv:1906.05950 [astro-ph.CO] .
  * Smyth  _et al._ (2020) N. Smyth, S. Profumo, S. English, T. Jeltema, K. McKinnon, and P. Guhathakurta, Phys. Rev. D 101, 063005 (2020), arXiv:1910.01285 [astro-ph.CO] .


# Classification of Large-Scale High-Resolution SAR Images with Deep Transfer
Learning

Zhongling Huang, Corneliu Octavian Dumitru, Zongxu Pan, Bin Lei, and Mihai
Datcu, This work was supported by the University of Chinese Academy of
Sciences (UCAS) Joint PhD Training Program scholarship. (Corresponding author:
Mihai Datcu.)Z. Huang, Z. Pan and B. Lei are with the Aerospace Information
Research Institute, Chinese Academy of Sciences, Beijing 100094, China, School
of Electronic, Electrical and Communication Engineering, University of Chinese
Academy of Sciences, Huairou District, Beijing 101408, China, and also with
Key Laboratory of Technology in Geo-spatial Information Processing and
Application System, Chinese Academy of Sciences, Beijing, China. (e-mail:
huangzhongling15@mails.ucas.ac.cn).Z. Huang, C.O. Dumitru and M. Datcu are
with the EO Data Science Department, German Aerospace Center (DLR), 82234
Wessling, Germany (e-mail<EMAIL_ADDRESS>Mihai.Datcu@dlr.de).

###### Abstract

The classification of large-scale high-resolution SAR land cover images
acquired by satellites is a challenging task, facing several difficulties such
as semantic annotation with expertise, changing data characteristics due to
varying imaging parameters or regional target area differences, and complex
scattering mechanisms being different from optical imaging. Given a large-
scale SAR land cover dataset collected from TerraSAR-X images with a
hierarchical three-level annotation of 150 categories and comprising more than
100,000 patches, three main challenges in automatically interpreting SAR
images of highly imbalanced classes, geographic diversity, and label noise are
addressed. In this letter, a deep transfer learning method is proposed based
on a similarly annotated optical land cover dataset (NWPU-RESISC45). Besides,
a top-2 smooth loss function with cost-sensitive parameters was introduced to
tackle the label noise and imbalanced classes’ problems. The proposed method
shows high efficiency in transferring information from a similarly annotated
remote sensing dataset, a robust performance on highly imbalanced classes, and
is alleviating the over-fitting problem caused by label noise. What’s more,
the learned deep model has a good generalization for other SAR-specific tasks,
such as MSTAR target recognition with a state-of-the-art classification
accuracy of 99.46%.

###### Index Terms:

TerraSAR-X, high-resolution SAR images, land cover classification, transfer
learning, label noise

##  I Introduction

Synthetic Aperture Radar (SAR) land use and land cover (LULC) classification
is an important step in SAR image interpretation. Recognizing different land
covers in SAR images requires an expertise in interpreting the complex
scattering characteristics under various circumstances such as different
wavelengths, incidence angles and target areas which makes the annotation of
large-scale SAR datasets time consuming. With the development of deep learning
algorithms, hierarchical features can be extracted automatically. Unsupervised
learning approaches, such as deep belief networks (DBNs) [1], and auto-
encoders (AEs) [2, 3], are widely applied to SAR images for classification.
With more detailed texture information, high-resolution (HR) PolSAR images are
used to learn spatial and polarimetric features jointly via deep convolutional
neural networks (DCNNs) for supervised classification which achieves state-of-
the-art performance [4].

However, many SAR LULC classification algorithms mainly focus on specific
tasks which may not generalize well on other SAR tasks, for example,
classifying vegetation and agricultural land covers in farmland areas [5], or
training an urban-specific network with two scene images in urban areas [6].
Although a big volume of SAR data all over the world is acquired every day, it
is challenging to automatically interpret land covers in SAR images with such
diversity and to obtain a well-generalized model for SAR image understanding.

In this letter, a global dataset of multi-level annotated HR SAR images [7]
has been experimented. These SAR images vary in geographic locations,
incidence angles, radiometric calibration factors, and ascending or descending
satellite orbit branches. With the difficulties of highly imbalanced
categories, geographic diversity and noisy annotations, our contributions are
as follows:

(1) A very deep residual network with transitive transfer learning method from
natural images to remote sensing and then to SAR data is used to understand a
large amount of land cover SAR images from all over the world.

(2) A cost-sensitive top-2 smooth combined loss function with cross-entropy
loss is proposed to alleviate the imbalanced problem and model bias to
mislabeled training samples.

(3) The proposed network is proved to have a good generalization on other SAR
datasets and tasks.

In the following, we will firstly introduce the TerraSAR-X (TSX) annotated
land cover dataset in Section II. Our transfer learning based method, as well
as the experimental results will be presented in Section III and Section IV,
followed by a conclusion in Section V.

Figure 1: Nine kinds of land cover of our annotated level-1 TSX dataset with
imbalanced classes. Statistically, 40 categories in NWPU-RESISC45 are
semantically similar to the TSX dataset. The number in brackets indicates the
assigned label of TSX, which can be indexed in TABLE A.I., Appendix [7]. Due
to the small number of patches in ”Military facilities” and indistinct label
”Unclassified”, only seven classes are applied in our experiments.

##  II Datasets

The HR SAR annotated land cover dataset [7] exploited in this paper has been
collected from the X-band TSX instrument [8]; we used horizontally polarized
(HH), multi-look ground range detected (MGD) products which were
radiometrically enhanced. The selected SAR images were taken in High
Resolution Spotlight mode with a resolution of 2.9 m, acquired with an
incidence angle between 20∘ and 50∘, in both descending and ascending pass
directions.

This annotated semantic dataset covers urban and non-urban areas from all over
the world, grouped into 46 collections in total. With a pixel spacing of 1.25
m, each SAR image was tiled into non-overlapping patches with a size of
160$\times$160 pixels, resulting in 200$\times$200 m2 on the ground so that a
patch covered a typical object on the ground. A single label was assigned to
each patch based on the dominating class. With three-level hierarchical
annotation, only level-1 labels are used in this paper.

Automatically interpreting the worldwide large-scale SAR dataset is
challenging due to the following reasons:

####  II-1 Highly imbalanced classes

As shown in Fig. 1, the imbalanced classes result in a big difference of data
volume for each category in level-1 annotation case, mainly due to the fact
that the dataset focuses on urban areas. This imbalance problem appears more
serious in level-2 and level-3 annotation.

####  II-2 Geographic diversity

Collected world-widely, different architectural styles of cities bring out a
large geographic diversity of images. The typical regional characteristics
result in high intra-class diversity and inter-class similarity as shown in
Fig. 2. Besides, the acquisition period of these images ranges from 2007 to
2012. During this period of time, some non-constant parameters such as the
radiometric calibration factor, will result in different characteristics of
the SAR images.

Figure 2: The high intra-class diversity of ”low-density residential areas” in
two different cities. Figure 3: Examples of label noise and model bias due to
the noisy labels.

####  II-3 Label noise

By covering an area of 200$\times$200 m2 on the ground, sometimes a patch
contains more than one land cover type. Even though the annotation was based
on the dominating class, there still existed some incomplete annotations or
some semantic biases [9] which may lead the trained model to over-fitting when
using these misleading samples. Fig. 3 gives an example of label noise in our
dataset.

Due to the low number of only 285 patches from two collections of ”Military
facilities” and the unlabeled patches of ”Unclassified”, we only used the
remaining seven land cover classes in this paper.

Another remote sensing LULC dataset used in this paper is NWPU-RESISC45 [10],
which has a land cover annotation similar to the seven classes of TSX dataset
used in our experiment, shown in Fig. 1, containing 31,500 images covering 45
scene classes.

##  III Proposed Method

Facing the difficulties mentioned above, we propose a deep transfer learning
based method in this part to tackle the SAR land cover classification issue.

Figure 4: Residual network architecture applied in our experiments.

###  III-A Transfer Learning with Deep Residual Networks

ResNet-18 [11] shown in Fig. 4, containing four residual blocks repeated
twice, with two convolution layers in each block, followed by an average
pooling layer and a fully-connected layer related to the class number of the
task, was applied in our experiments. Such a very deep well-trained neural
network (with a data volume of about 43.7 MB) has a great ability in
extracting hierarchical features but is difficult to train from scratch.
Transfer learning provides a better initialization to the parameters in deep
layers, making the optimization easier and speeding up the training process.
The ImageNet pre-trained model is popular as a transferred source. However,
the large difference lying between SAR land cover patches and natural object
images may harm the transferability, making it necessary to find a bridge to
narrow the gap [12].

Rather than using the ImageNet pre-trained model as most remote sensing
applications do, we proposed a transitive transfer learning method with the
help of the NWPU-RESISC45 dataset [10] to build a bridge linking the natural
and SAR land cover images. ResNet-18 pre-trained on ImageNet [11] was applied
as our initial model and then the convolution layers were fine-tuned with the
NWPU-RESISC45 dataset, obtaining the remote sensing enhanced model as the
transfer source for our applications.

###  III-B Imbalanced Classification

When we trained a model with a dataset of highly imbalanced classes, we
applied some techniques such as sampling and cost-sensitive learning [13]. For
evaluation, we considered some metrics beyond accuracy, like precision, recall
and F1-score.

Each time when we input a batch of samples to the network, over-sampling of
the minority classes and under-sampling of the majority classes were applied
to ensure a balanced distribution in that batch. Moreover, the training
samples were weighted according to the corresponding proportions of the
dataset in order to penalize the misclassification of the minority classes
more heavily than for the majority classes.

###  III-C Top-2 Smooth Loss Function

Here, we introduced the top-2 smooth loss function [14] to reduce the label
noise in our dataset. With the classification of seven types of land cover in
level-1, we mainly focused on the top-2 predictions
$\mathrm{\mathbf{\bar{y}}}\in\mathbb{R}^{2}$. We define the label space
$\mathcal{Y}=\\{1,2,...,n\\}$ and the ground truth label $y\in\mathcal{Y}$,
the top-2 predictions $\mathrm{\mathbf{\bar{y}}}\in\mathcal{Y}^{(2)}$. The
subset of $\mathcal{Y}^{(2)}$ where the ground truth $y$ is included, is
written as $\mathcal{Y}_{y}^{(2)}$. The top-2 smooth loss function with
temperature can be defined as:

$\begin{split}L_{2,\tau}(\mathbf{s},y)=&\tau\log[\sum\limits_{\mathrm{\mathbf{\bar{y}}}\in\mathcal{Y}^{(2)}}\exp(\frac{1}{\tau}(\Delta(\mathrm{\mathbf{\bar{y}}},y)+\frac{1}{2}\sum\limits_{j\in\mathrm{\mathbf{\bar{y}}}}s_{j}))]\\\
&-\tau\log[\sum\limits_{\mathrm{\mathbf{\bar{y}}}\in\mathcal{Y}_{y}^{(2)}}\exp(\frac{1}{2\tau}\sum\limits_{j\in\mathrm{\mathbf{\bar{y}}}}s_{j})].\end{split}$
(1)

where $\mathbf{s}\in\mathbb{R}^{n}$ denotes a vector of scores of $n$ labels,
$\tau$ the temperature parameter and
$\Delta(\mathrm{\mathbf{\bar{y}}},y)\triangleq\mathbbm{1}(y\notin\mathbf{\bar{y}})$
denotes the margin. The top-2 smooth loss function takes two most likely
predictions into consideration, making it possible to prevent over-fitting
with label noise when learning only with the cross-entropy loss of a single
label.

In the proposed method, the top-2 smooth loss was combined with the cross-
entropy loss with a trade-off $\lambda$ rather than totally replacing it as
[14] did, because $k=2$ is not much smaller than $n=7$ in our case. The top-2
loss was used as a constraint to alleviate model bias towards the single
misleading label when there was more than one kind of land cover on the
ground. Finally, adding an L2 regularization term and the loss weight $w_{y}$,
the combined loss function can be written as:

$\begin{split}L_{2,\tau,\mu,\lambda}(\mathbf{s},y)=&(1-\lambda)(-\sum\limits_{i}s_{i}\log\frac{\exp(s_{i})}{\sum\limits_{j=1}^{C}\exp(s_{j})})\\\
&+w_{y}\lambda
L_{2,\tau}(\mathbf{s},y)+\mu\sum\limits_{i}\|\beta_{i}\|^{2},\end{split}$ (2)

where $\beta_{i}$ denotes the weights in each layer, and $\lambda$, $\tau$ and
$\mu$ are three hyper-parameters which control the trade-off between the top-2
and cross-entropy losses, temperature tuning, and the penalty of the
regularization term. The cost weight $w_{y}$ related to the class $y$ is given
by:

$w_{y}=\frac{1}{n-1}(1-\frac{N_{y}}{\sum_{i=1}^{n}{N_{i}}}),$ (3)

where $n$ denotes the class number and $N_{i}$ denotes the number of samples
in the $i^{th}$ class.

In addition, we use bagging [15], training five sub-models with re-sampling
data and voting the predictions of all sub-models in order to make the trained
model more robust and general.

##  IV Experimental Results and Discussions

We firstly introduce our experimental settings in Section IV-A. In Section
IV-B and IV-C we analyze the effectiveness of different transfer sources and
the combined top-2 smooth loss. Finally, we use three examples to demonstrate
the good generalization of the proposed SAR specific deep model in IV-D.

###  IV-A Experimental Settings

On the one hand, we randomly chose 200 patches in each category from
collection 01 to 36 as Test Data 1 and the remaining as Training Data. On the
other hand, patches of the full scene images from another seven collections
were used as Test Data 2 which had a different distribution with Training Data
and Test Data 1. The details of training and test data settings are presented
in Table I.

TABLE I: Training and Test Dataset Category | Settlements |  | Industrial
---
areas
| Public
---
transportations
| Agricultural
---
land
| Natural
---
vegetation
| Bare
---
ground
| Water
---
bodies
Total
Training Data | 24,930 | 2,979 | 4,485 | 6,029 | 4,911 | 2,240 | 6,826 | 52,400
Test Data 1 | 200 | 200 | 200 | 200 | 200 | 200 | 200 | 1,400
Test Data 2 | 3,155 | 563 | 928 | 925 | 1,383 | 1,319 | 2,133 | 10,406

The TSX detected images are represented by 16 bits spanning a large dynamic
range. In our experiments, the 16 bit data after log transform was fed into
the network to prevent the gradients from vanishing during training. When
transferring the pre-trained model, only the last fully-connected layer was
adapted to fit the new task.

The hyper-parameters of $\lambda$, $\tau$ and $\mu$ which control the trade-
off between the top-2 and cross-entropy losses, the temperature tuning, and
the penalty of the regularization term, were empirically set to 0.2, 1.0 and
0.25, respectively. Adam optimization algorithm [16] was applied and the
initial learning rate was set to $10^{-5}$ when using the top-2 loss, and to
$10^{-4}$ otherwise.

###  IV-B The Transitive Transferring Model

We compared the results of overall classification accuracy in Test Data 1 and
the average stopping time of training a ResNet-18 network for the following
cases to demonstrate the effectiveness of transferring from proper source data
to SAR land cover classification:

####  IV-B1 RD-$\mathbf{TSX}$

By randomly initializing all layers of ResNet-18, the network was trained with
the TSX dataset from scratch.

####  IV-B2 TL-$\mathbf{ImageNet}$2$\mathbf{TSX}$

By using the pre-trained model of ImageNet, the network was fine-tuned with
the TSX dataset with only the last layer related to the classes trained from
scratch.

####  IV-B3 TL-$\mathbf{ImageNet}$2$\mathbf{RS}$2$\mathbf{TSX}$

The hierarchical convolution layers of ImageNet pre-trained model were fine-
tuned with NWPU-RESISC45 optical remote sensing images for 45 types of land
cover classification and then transferred to the TSX dataset.

Figure 5: The lines denote the loss values on training and validation data
before early stopping, and the bars denote the validation accuracy of three
experimental settings.

Fig. 5 illustrates that training such a very deep network like ResNet-18 with
the TSX dataset from scratch is difficult and time-consuming. Consequently, it
is necessary to find a better initialization of millions of parameters.
Although the ImageNet pre-trained model has been widely used for remote
sensing applications, it is controversial to directly transfer the natural
images to SAR data due to the huge differences between them. Here, we took
advantage of the remote sensing images with similar categories to our SAR land
cover dataset to narrow the gap between natural images and SAR data. We see a
performance improvement both in test accuracy and training time.

###  IV-C Combined with Top-2 Loss

Besides applying some strategies such as Dropout and L2 regularization to
avoid over-fitting during training, we found that label noise may play an
important role with a negative influence on the trained model. To address this
prediction bias issue, we introduced a top-2 smooth loss function to constrain
the cross-entropy from introducing a bias to the noisy labels. We use
different metrics of top-1, top-2 overall accuracy, and also F1-score to
evaluate the results which is defined by $F_{1}=2/(p^{-1}+r^{-1})$, where
$p=\frac{TP}{TP+FP}$, $r=\frac{TP}{TP+FN}$, and $TP$, $FP$, $FN$ denote the
true positive, false positive and false negative samples, respectively.

TABLE II: Classification Result on Test Data 1 and Test Data 2 |  Eval Loss | ce loss | top-2 loss | combined loss
---|---|---|---|---
Test Data 1 | F1 score | 0.8318 | 0.7954 | 0.8327
top-1 acc | 83.14% | 79.5% | 83.21%
top-2 acc | 95.36% | 93.29% | 95.00%
Test Data 2 | F1 score | 0.709 | 0.684 | 0.72
top-1 acc | 74.09% | 73.24% | 76.78%
top-2 acc | 91.26% | 91.25% | 92.62%
training time | 16 epochs | 6 epochs | 9 epochs

The top-2 smooth loss focuses on two most likely land cover classes, thus
penalizing the predictions only when neither of the top-2 likely labels are
the ground truth. Table II shows the training time of using different loss
functions and the test results of F1-score, top-1, and top-2 accuracies both
on Test Data 1 and 2. The top-2 accuracy is the percentage of time that the
classifier gives the correct class in the top-2 probabilities. We found that
training with top-2 smooth loss only would lead to a fast convergence of the
deep network to speed up training. However, the top-2 smooth loss was too
relaxing for the optimization resulting in much smaller gradients after a
period of time and weights being difficult to update. As a consequence, we
combined the top-2 smooth loss with the cross-entropy as a constraint to
achieve a faster convergence with nine epochs compared with 16 epochs of
learning with cross-entropy. With the improvement in F1-score and accuracy,
together with the obvious speed-up in training time, the combined loss
performed considerably better.

Figure 6: The classification results on seven new collections of Test Data 2,
learning with cross-entropy, and the combined loss shown as blue/red dots,
respectively.

###  IV-D Generalization to Other Data and Tasks

We could demonstrate the good generalization of our SAR deep ResNet model by
testing two tasks: (1) land cover classification on new collections of Test
Data 2; (2) MSTAR [17] target recognition without data augmentation, as shown
in Table III.

TABLE III: Transfer to New Tasks on SAR Data Dataset | Task | Classes | Transfer Result
---|---|---|---
Test Data 2 | scene classification | 7 | 76.78%
MSTAR [17] | target recognition | 10 | 99.46%

Fig. 6 shows the test results of each new collection in Test Data 2 which has
a different distribution with the training data. In most new collections, the
performance of the combined loss was better than cross-entropy (2.69% higher
on average) which demonstrates that the trained model with combined loss
reduces over-fitting due to label noise and improves generalization on new
dataset.

Furthermore, when we transferred the proposed model to an MSTAR [17] 10-class
recognition task, we achieved a state-of-the-art overall accuracy of 99.46%,
compared with 98.02% using an ImageNet pre-trained model, and 99.09% by other
CNN-based methods [18].

##  V Conclusion

In this letter, a deep transfer learning method combining the top-2 smooth
loss was proposed to solve the land cover classification problem for a large-
scale HR TSX dataset, with highly imbalanced classes, geographic diversity and
noisy labels. Rather than applying the ImageNet pre-trained model of ResNet-18
to SAR images directly, we used an optical remote sensing land cover dataset
to narrow the gap between SAR and natural images which results in a
significant improvement on feature transferability. The proposed combined loss
function is successful in reducing model bias to noisy labels and our proposed
deep ResNet model shows a good generalization which allows it to be exploited
for other SAR-specific tasks.

## Acknowledgment

We thank the TerraSAR-X Science Service System for the provision of images
(Proposals MTH-1118 and LAN-3156), and University of Chinese Academy of
Sciences (UCAS) Joint PhD Training Program for the provision of a scholarship.
We would also thank Dr. Gottfried Schwarz for improving this letter.

## References

  * [1] Q. Lv, Y. Dou, X. Niu, J. Xu, J. Xu, and F. Xia, “Urban land use and land cover classification using remotely sensed SAR data through deep belief networks,” _Journal of Sensors_ , vol. 2015, 2015.
  * [2] L. Zhang, W. Ma, and D. Zhang, “Stacked sparse autoencoder in PolSAR data classification using local spatial information,” _IEEE Geosci. Remote Sens. Lett._ , vol. 13, no. 9, pp. 1359–1363, Sept 2016.
  * [3] Z. Huang, Z. Pan, and B. Lei, “Transfer learning with deep convolutional neural network for SAR target classification with limited labeled data,” _Remote Sens._ , vol. 9, no. 9, 2017.
  * [4] W. Wu, H. Li, L. Zhang, X. Li, and H. Guo, “High-resolution PolSAR scene classification with pretrained deep convnets and manifold polarimetric parameters,” _IEEE Trans. on Geosci. Remote Sens._ , vol. 56, no. 10, pp. 6159–6168, Oct 2018.
  * [5] C. Hütt, W. Koppe, Y. Miao, and G. Bareth, “Best accuracy land use/land cover (LULC) classification to derive crop types using multitemporal, multisensor, and multi-polarization SAR satellite images,” _Remote Sens._ , vol. 8, no. 8, 2016.
  * [6] J. Zhao, W. Guo, B. Liu, Z. Zhang, W. Yu, and S. Cui, “Preliminary exploration of SAR image land cover classification with noisy labels,” in _Proc. IEEE International Geoscience and Remote Sensing Symposium (IGARSS)_ , Jul 2017, pp. 3274–3277.
  * [7] C. O. Dumitru, G. Schwarz, and M. Datcu, “Land cover semantic annotation derived from high-resolution SAR images,” _IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens._ , vol. 9, no. 6, pp. 2215–2232, 2016.
  * [8] “TerraSAR-X basic product specification document, issue 1.9,” https://sss.terrasar-x.dlr.de/docs/TX-GS-DD-3302.pdf, accessed December 10, 2018.
  * [9] A. M. M. de Oca, R. Bahmanyar, N. Nistor, and M. Datcu, “Earth observation image semantic bias: A collaborative user annotation approach,” _IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens._ , vol. 10, no. 6, pp. 2462–2477, Jun 2017.
  * [10] G. Cheng, J. Han, and X. Lu, “Remote sensing image scene classification: Benchmark and state of the art,” _Proc. IEEE_ , vol. 105, no. 10, pp. 1865–1883, Oct 2017.
  * [11] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proc. The IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , Jun 2016.
  * [12] Z. Huang, Z. Pan, and B. Lei, “What, where, and how to transfer in sar target recognition based on deep cnns,” _IEEE Trans. on Geosci. Remote Sens._ , pp. 1–13, 2019.
  * [13] M. Kukar, I. Kononenko _et al._ , “Cost-sensitive learning with neural networks.” in _ECAI_ , 1998, pp. 445–449.
  * [14] L. Berrada, A. Zisserman, and M. P. Kumar, “Smooth loss functions for deep top-k classification,” in _Proc. IEEE International Conference on Learning Representations (ICLR)_ , 2018.
  * [15] E. Bauer and R. Kohavi, “An empirical comparison of voting classification algorithms: Bagging, boosting, and variants,” _Mach. Learn._ , vol. 36, no. 1, pp. 105–139, Jul 1999.
  * [16] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _3rd International Conference on Learning Representations, ICLR 2015_ , 2015\.
  * [17] “The air force moving and stationary target recognition database,” https://www.sdms.afrl.af.mil/datasets/mstar/, accessed 3 February, 2016\.
  * [18] S. Chen, H. Wang, F. Xu, and Y.-Q. Jin, “Target classification using the deep convolutional networks for SAR images,” _IEEE Trans. Geosci. Remote Sens._ , vol. 54, no. 8, pp. 4806–4817, 2016.