# Joint Ranging and Phase Offset Estimation for Multiple Drones using ADS-B
Signatures

Mostafa Mohammadkarimi, , Geert Leus, ,
and Raj Thilak Rajan The authors are with the Faculty of Electrical
Engineering, Mathematics and Computer Science, Delft University of Technology,
2628 CD Delft, The Netherlands (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>R.T.Rajan@tudelft.nl). This work is partially funded by
the European Leadership Joint Undertaking (ECSEL JU), under grant agreement No
876019, and the ADACORSA project - “Airborne Data Collection on Resilient
System Architectures.” (https://adacorsa.eu/).

###### Abstract

A new method for joint ranging and Phase Offset (PO) estimation of multiple
drones/aircrafts is proposed in this paper. The proposed method employs the
superimposed uncoordinated Automatic Dependent Surveillance–Broadcast (ADS-B)
packets broadcasted by drones/aircrafts for joint range and PO estimation. It
jointly estimates range and PO prior to ADS-B packet decoding; thus, it can
improve air safety when packet decoding is infeasible due to packet collision.
Moreover, it enables coherent detection of ADS-B packets, which can result in
more reliable multiple target tracking in aviation systems using cooperative
sensors for detect and avoid (DAA). By minimizing the Kullback–Leibler
Divergence (KLD) statistical distance measure, we show that the received
complex baseband signal coming from $K$ uncoordinated drones/aircrafts
corrupted by Additive White Gaussian Noise (AWGN) at a single antenna receiver
can be approximated by an independent and identically distributed (i.i.d.)
Gaussian Mixture (GM) with $2^{K}$ mixture components in the two-dimensional
(2D) plane. While direct joint Maximum Likelihood Estimation (MLE) of range
and PO from the derived GM Probability Density Function (PDF) leads to an
intractable maximization, our proposed method employs the
Expectation–Maximization (EM) algorithm to estimate the modes of the 2D
Gaussian mixture followed by a reordering estimation technique through
combinatorial optimization to estimate range and PO. An extension to a
multiple antenna receiver is also investigated in this paper. While the
proposed estimator can estimate the range of multiple drones/aircrafts with a
single receive antenna, a larger number of drones/aircrafts can be supported
with higher accuracy by the use of multiple antennas at the receiver. The
effectiveness of the proposed estimator is supported by simulation results. We
show that the proposed estimator can jointly estimate the range of three
drones/aircrafts accurately.

###### Index Terms:

Range estimation, phase offset, cooperative navigation,
expectation–maximization (EM), Gaussian mixture (GM), ADS-B, multiple receive
antennas, detect and avoid (DAA).

##  I Introduction

Automatic Dependent Surveillance–Broadcast (ADS-B) is one of the two Automatic
Dependent Surveillance (ADS) systems that tracks aircrafts without employing
radar. It is intended to improve traffic surveillance capabilities by sharing
accurate aircraft position information between pilots and air traffic
controllers. In the ADS-B system, aircrafts regularly and asynchronously
broadcast their real-time position information, velocity, and identification
to no specific receiver using a transponder, typically combined with a Global
Positioning System (GPS), to transmit highly accurate positional information
to traffic controllers and directly to other aircrafts. This transmission is
known as ADS-B Out and its accuracy and update rate are much greater than
conventional primary radar surveillance. ADS-B In is the reception of the
ADS-B packet by an aircraft [1, 2, 3].

The ADS-B system can be considered as a promising solution for autonomous
navigation of drones in future drone technology. It can be employed in detect
and avoid (DAA) systems to improve air safety. In order to avoid aviation
accidents, each drone needs to be aware of the position and speed of the
surrounding drones so that it can keep a safe separation distance with the
other drones. This safe separation can be achieved by a cooperative sensor
system, such as the ADS-B system. In this solution, drones are equipped with
GPS, an Inertial Measurement Unit (IMU), and a miniaturized transponder and
they broadcast their real-time position information, which can be employed by
the surrounding drones or Ground Controllers (GCs) to maintain a safe
operation distance of drones at low altitude and congested airspace. The
design of communication schemes employing Wi-Fi and other wireless
technologies for autonomous navigation based on a cooperative sensor system
was investigated in [4] and [5].

One of the main challenges in the employment of a cooperative sensor system,
such as ADS-B, for future drone technology is packet collisions due to a
larger number of drones compared to aircrafts in the airspace. As the number
of drones in the airspace increases, the probability of packet collision
increases. The ADS-B system in its current form cannot handle packet
collision; hence, a large number of packets are lost. Packet loss means less
information and more uncertainty for the surrounding drones, resulting in less
air safety.

While the performance of ADS-B for a single aircraft has been widely studied
in the literature [1], ADS-B for multiple aircrafts in the presence of packet
collisions and intentional/unintentional interference has not been thoroughly
investigated. Coherent and non-coherent detection of ADS-B packets for a
single aircraft was studied in [6] and closed form expressions for the packet
error rate were derived. The challenge of preamble detection for ADS-B for a
single aircraft at low Signal-to-Noise Ratios (SNRs) was studied in [7], where
the authors developed a constant false alarm rate detection algorithm.

A signal separation algorithm based on time-domain cancellation for two ADS-B
packets was proposed in [8]. This is the first attempt to handle packet
collisions in ADS-B. To mitigate the interference caused by other ADS-B
signals, the idea of sectorized antennas and blind adaptive beamforming was
proposed in [9]. Security issues for the ADS-B system were studied in [10, 11,
12]. A packet loss tolerant security framework for the ADS-B system was
investigated in [10]. A complete cryptographic solution to ensure the privacy
and integrity of ADS-B messages was proposed in [11]. Security issues in ADS-B
are studied in [12]. Recently, the ADS-B signal has also been employed as a
Signal-of-Opportunity (SoOp) to develop alternative positioning, navigation,
and timing systems. A synchronization method for large-scale IoT networks
using ADS-B packets of multiple aircrafts was proposed in [13]. Wide-area
multilateration was developed in [14] that implements a time difference of
arrival (TDOA) positioning by employing ADS-B packets. The effects of jamming
on ADS-B receivers were analyzed and discussed using real data in [15]. The
authors reported a strong reduction of the coverage, a reduced detection
probability, and a possible overload of the processor due to jamming. In [16],
the authors used different machine learning-based algorithms, including KNN,
Naive Bayes, and logistic regression to secure the ADS-B system against
deliberate interference, i.e., jamming.

Information extraction from collided ADS-B packets, such as range, velocity,
Angle of Arrival (AoA), etc. is of interest for drone/aircraft avionics,
specifically for DAA. To the best of the authors’ knowledge, joint estimation
of range and PO for multiple drones/aircrafts using collided ADS-B packets has
not been investigated yet. Existing blind source separation algorithms suffer
from sign ambiguity and require the number of receive antennas not to be lower
than the number of drones/aircrafts [17, 18, 19]. In this paper, we derive the
Maximum Likelihood (ML) cost function for joint range and PO estimation of
multiple drones/aircrafts. Then, a simple solution based on the EM algorithm
is proposed. Our proposed estimator enables avionic systems employing
cooperative sensors to obtain range information of the surrounding
drones/aircrafts while ADS-B packet decoding is not feasible due to collision.
Furthermore, the estimated range and PO can be employed for coherent detection
of ADS-B packets, which offers higher detection performance. Our proposed
estimator can estimate the range of multiple drones/aircrafts with a single
receive antennas.

(a) Range estimation in a GC.

(b) Range range estimation in a flying drone.

Figure 1: Range estimation using the asynchronous ADS-B In signatures of the
drones at the receiver.

###  I-A Contributions

The contributions of the paper can be categorized as:

  * •

By minimizing the KLD statistical distance measure, we show that the received
complex baseband signal coming from $K$ drones/aircrafts asynchronously
broadcasting ADS-B packets can be approximated by an i.i.d. GM random variable
with $2^{K}$ mixture components in the 2D plane, which is independent of the
arrival time of the ADS-B packets at the receiver.

  * •

By using the approximate PDF, the ML cost function for the joint range and PO
estimation of multiple drones/aircrafts is derived.

  * •

A low complexity suboptimal joint range and PO estimator by employing EM is
proposed.

  * •

We extend the derived maximum likelihood cost function and the proposed EM-
based joint estimator to the case of multiple receive antennas.

  * •

To remove the ambiguity in the order of the estimated modes of the 2D GM in
the proposed EM-based joint estimator, we introduce a permutation-based
combinatorial optimization procedure in order to reorder the estimated modes.
Then, a low complexity solution for the reordering is proposed.

  * •

A further performance improvement for multiple receive antennas by using
outlier detection algorithms is proposed.

  * •

The proposed algorithm can jointly estimate range and PO of multiple
drones/aircrafts. Moreover, it is applicable to both near and far field
ranging.

###  I-B Notations

The identity matrix, all-zero vector, and all-one vector of length $N$ are
denoted by ${\bf{I}}_{N}$, ${\bf{0}}_{N}$, and ${\bf{1}}_{N}$, respectively.
Throughout the paper, $(\cdot)^{*}$, $(\cdot)^{T}$, and $(\cdot)^{\rm{H}}$
show the complex conjugate, transpose, and Hermitian transpose, respectively.
Also, $|\cdot|$, $\lfloor\cdot\rfloor$, $*$, and $\otimes$ represent the
absolute value operator, the floor function (greatest integer value), linear
convolution, and Kronecker product, respectively. $\mathbb{E}\\{\cdot\\}$ is
the statistical expectation, $\hat{{x}}$ is an estimate of $x$. The complex
Gaussian distribution with mean vector $\bm{\mu}$ and covariance matrix
$\bm{\Sigma}$ is denoted by
${\cal{C}\cal{N}}\big{(}\bm{\mu},\bf{\Sigma}\big{)}$. The continuous uniform
distribution between $a$ and $b$ and the discrete uniform distribution between
$N_{1}$ and $N_{2}$ are denoted by $\mathcal{U}_{\rm c}[a,b]$ and
$\mathcal{U}_{\rm d}[N_{1},N_{2}]$, respectively.

The remainder of the paper is organized as follows. Section II introduces the
system model. Section III describes the GM distribution approximation by
minimizing the KLD. In Section IV, the maximum likelihood cost function and
the EM-based joint ranging and PO estimation algorithm are analytically
derived. Reordering estimation for the proposed EM-based joint estimator
through permutation-based combinatorial optimization is investigated in
Section V. Joint estimation by taking the advantage of diversity gain through
multiple receive antennas is discussed in Section VI. Simulation results are
provided in Section VII, and conclusions are drawn in Section VIII.

Figure 2: An ADS-B packet is composed of a preamble and data in the PPM form.

##  II System Model

We consider $K$ drones broadcasting ADS-B packets through their transponders
to the GC111The GC can be a simple Software-defined Radio (SDR) receiver.
Number of drones/aircrafts, $K$, can be determined by employing model-order
selection techniques [20, 21] prior to joint ranging and PO estimation.  and
also directly to other flying drones/aircrafts. Typical scenarios for signal
reception at the GC and the flying drone are shown in Fig. 1. The packets are
transmitted at $1090$ Mhz and use PPM at a rate of $1$ Mbit per second. A mode
S ADS-B packet consists of a preamble followed by a data block as shown in
Fig. 2. The preamble is used to synchronize transmission timing between the
drones and the receiver. In the data field, the first 5 bits are used to
indicate the type of data since an ADS-B system does not carry all the
necessary and critical information, such as speed, location, altitude and so
on in a single packet.

The next 3-bit capability field serves as an additional identifier. The
subsequent $24$ bits in the data block represent the International Civil
Aviation Organization (ICAO) code, which can uniquely identify one
drone/aircraft. We use the ICAO code to classify the information from
different drones. ADS-B adds flight state information into the next 56-bit
data field. The last field is a $24$-bit parity check for receivers to
validate the correctness of the received packets.

(a) Observation window for joint range and PO estimation.

(b) Equivalent observation window for joint range and PO estimation.

Figure 3: ADS-B packets are received with different delays at the receiver.

It is assumed that the drones/aircrafts asynchronously broadcast their ADS-B
packets every $T_{\rm P}$ seconds. We consider an observation window of length
$T_{\rm{w}}=T_{\rm P}$ for parameter estimation at the receiver (Fig. 3(a)).
To make the joint ranging and PO estimation independent of the arrival time of
the ADS-B packets at the receiver, we approximate the received samples in the
observation interval by an i.i.d. complex random variable as it will be
explained in Section III. Hence, without loss of generality, we can consider
the ADS-B packet reception in Fig. 3(b) to simplify the derivation of the
joint range and PO estimator. In this case, it is assumed that the ADS-B
packet of the $k$th drone/aircraft with a packet length of $T_{\rm{A}}$ is
received at the receiver with time delay $\tau_{k}\in[0,\tau_{\rm max}]$ in
the timing reference of the receiver, where $\tau_{k}$ is unknown and random
in each observation window of length $T_{\rm{w}}=T_{\rm P}$, and $\tau_{\rm
max}=T_{\rm P}-T_{\rm A}$ is the maximum time delay of a packet.

The received complex baseband signal at the ADS-B receiver can be expressed as

$\displaystyle
y(t)=\sum_{k=1}^{K}\sqrt{P_{k}L_{k}}x_{k}(t-\tau_{k})e^{j(2\pi\Delta
f_{k}t+\theta_{k})}+w(t),$ (1)

where $t\in[0,T_{\rm w}]$, and where $P_{k}$, $x_{k}(t)$, $w(t)$, and $L_{k}$
denote, the transmit power by the $k$th drone at time $t$, the transmit packet
by the $k$th drone, the additive noise with Power Spectral Density (PSD)
$N_{0}$ over the frequency $f\in[-B/2,B/2]$ ($B$ is the bandwidth of the band-
pass filter at the receiver), and the path loss between the $k$th
drone/aircraft and the ADS-B receiver, respectively. For free-space path loss,
we have

$\displaystyle L_{k}\triangleq\Big{(}\frac{\lambda_{\rm{c}}}{4\pi
r_{k}}\Big{)}^{2},$ (2)

where $r_{k}$ is the range between the $k$th drone and the receiver,
$\lambda_{\rm{c}}\triangleq c/f_{\rm c}$ is the wavelength of the carrier
wave, $c$ denotes the speed of light, and $f_{\rm c}$ represents the carrier
frequency. For the ADS-B system, since $f_{\rm c}=1090$MHz, we have
$\lambda_{\rm{c}}\approx 0.2752$ m. In (1), $\Delta f_{k}$ and $\theta_{k}$
further denote the Carrier Frequency Offset (CFO) and PO of the $k$th drone in
the observation window. In this paper, we consider that $\Delta f_{k}T_{\rm
w}\ll 1$ and $\exp(j2\pi\Delta f_{k}t)\approx 1$ for $t\in[0,T_{\rm w}]$,
$k=1,2,\cdots,K$; thus, we can write

$\displaystyle
y(t)\approx\sum_{k=1}^{K}\sqrt{P_{k}L_{k}}x_{k}(t-\tau_{k})e^{j\theta_{k}}+w(t).$
(3)

Assumption 1: We assume that $P_{k}$, $k=1,2,\ldots,K$, is known at the
receiver, and that $P_{1}L_{1}>P_{2}L_{2}>\ldots,P_{K}L_{K}$.

Since the bit-rate for ADS-B systems is $1{\rm M\ bit}/{\rm s}$, a sampling
rate of $f_{s}=\frac{1}{T_{\rm{s}}}=2{\rm M\ sample/s}$ is sufficient to
capture the bit transitions, detect the preamble and decode packets. A typical
sampled ADS-B signal is shown in Fig. 4 for visualization. The discrete-time
received baseband signal after sampling, i.e., $y_{n}\triangleq
y(nT_{\rm{s}})$, $n=0,1,\ldots,N$, where $N=239+M$,222A sampling rate of
$2{\rm M\ sample/s}$ results in $240$ samples per ADS-B packet.
$M\triangleq\lfloor\frac{\tau_{\rm max}}{T_{\rm s}}\rfloor$, can be written in
vector form as

$\displaystyle{\bf y}=\sum_{k=1}^{K}h_{k}{\bf x}_{k}+{\bf
w}=\sum_{k=1}^{K}{\bf z}_{k}+{\bf w}={\bf g}+{\bf w},$ (4)

where ${\bf g}\triangleq\sum_{k=1}^{K}h_{k}{\bf x}_{k}=\sum_{k=1}^{K}{\bf
z}_{k}$, ${\bf z}_{k}=h_{k}{\bf x}_{k}$,

$\displaystyle{\bf y}$ $\displaystyle\triangleq\big{[}y_{0}\ y_{1}\ \dots\
y_{N}\big{]}^{T},$ (5a) $\displaystyle{\bf g}$
$\displaystyle\triangleq\big{[}g_{0}\ g_{1}\ \dots\ g_{N}\big{]}^{T},$ (5b)
$\displaystyle{\bf w}$ $\displaystyle\triangleq\big{[}w_{0}\ w_{1}\ \dots\
w_{N}\big{]}^{T},$ (5c) $\displaystyle{\bf z}_{k}$
$\displaystyle\triangleq\big{[}z_{k,0}\ z_{k,1}\ \cdots\ z_{k,N}\big{]}^{T},$
(5d)

$w_{n}\triangleq w(nT_{\rm{s}})$, $h_{k}\triangleq\beta_{k}e^{j\theta_{k}}$,
$\beta_{k}\triangleq\sqrt{P_{k}L_{k}}$, and

$\displaystyle{\bf x}_{k}$ $\displaystyle\triangleq\big{[}x_{k,0}\ x_{k,1}\
\cdots\ x_{k,N}\big{]}^{T}$ (6)
$\displaystyle\triangleq\big{[}{\bf{0}}_{m_{k}}^{T}\ {\bf s}^{T}\ {\bf
d}_{k}^{T}\ {\bf{0}}_{M-{m_{k}}}^{T}\big{]}^{T}.$

In (6), $x_{k,n}\triangleq x_{k}(nT_{\rm{s}}-\tau_{k})$, ${\bf
d}_{k}\in\\{0,1\\}$, is the data vector of the $k$th drone with a length of
$224$ symbols, ${\bf s}$ is the preamble vector given by ${\bf s}=\big{[}1\ 0\
1\ 0\ 0\ 0\ 0\ 0\ 0\ 1\ 0\ 1\ 0\ 0\ 0\ 0\big{]}^{T},$ and
$M\triangleq\lfloor\frac{\tau_{\rm max}}{T_{\rm s}}\rfloor$ and
$m_{k}\triangleq\lfloor\frac{\tau_{k}}{T_{\rm s}}\rfloor$ denote the maximum
possible integer delay for a drone/aircraft and the integer delay of the $k$th
drone/aircraft, respectively. The integer delays $m_{k}$, $k=1,2,\cdots,K$,
are unknown at the receiver and their values change from one observation
window to another due to mobility of the drones/aircrafts. The vector ${\bf
w}\triangleq[w_{0}\ w_{1}\ \ldots\ w_{N}]^{T}$ in (4) denotes the zero-mean
additive white Gaussian noise (AWGN) with covariance matrix
$\mathbb{E}\big{\\{}{\bf w}{\bf w}^{T}\big{\\}}=\sigma_{\rm w}^{2}{\bf
I}=N_{0}B{\bf I}$.

Let us define hypothesis $H_{m}^{k}$ as follows

$\displaystyle H_{m}^{k}:\,\,\,\,\,\ {\bf x}_{k}=\big{[}{\bf{0}}_{m}^{T},{\bf
s}^{T},{\bf d}_{k}^{T},{\bf{0}}_{M-{m}}^{T}\big{]}^{T},$ (7)

which represents the ADS-B packet of the $k$th drone arriving at the receiver
with integer delay $m_{k}=m\in\\{0,1,\ldots,M\\}$.

Figure 4: Amplitude of the noisy ADS-B packet after sampling.

##  III Distribution Approximation

To remove the dependency of joint ranging and PO estimation from the unknown
arrival time of the ADS-B packets at the receiver, we approximate the received
noisy samples in the observation interval by i.i.d. complex random variables.
In this paper, we show that an i.i.d. two-dimensional (2D) GM model can be
used to model the received baseband superimposed and noise corrupted signal at
the receiver.

######  Theorem 1.

By maximizing the KLD criterion, the elements of the ADS-B packet of the $k$th
drone, i.e., ${\bf x}_{k}=[{x}_{k,0}\ {x}_{k,1}\ \ldots\
{x}_{k,N}]^{T}=\big{[}{\bf{0}}_{m_{k}}^{T}\ {\bf s}^{T}\ {\bf d}_{k}^{T}\
{\bf{0}}_{M-{m_{k}}}^{T}\big{]}^{T}$ can be approximated by an i.i.d. random
variable that are Bernoulli distributed with Probability Mass Function (PMF)

$\displaystyle q\big{(}x;p\big{)}=\begin{cases}p&\text{if }x=0,\\\
1-p&\text{if }x=1,\end{cases}$ (8)

where

$\displaystyle p=\frac{M+124}{M+240}.$ (9)

###### Proof.

See Appendix A. ∎

As seen in Theorem 1, the approximated PMF $q$ does not depend on $H_{m}^{k}$.

Since ${\bf z}_{k}=h_{k}{\bf x}_{k}$ is the scaled version of ${\bf x}_{k}$
and the parameter of the Bernoulli distribution, $p$, is independent of
$m_{k}$, the elements of the complex vector ${\bf z}_{k}$ can be approximated
by i.i.d. complex random variables $Z_{k}$ with PDF
$f_{{Z}_{k}}\big{(}z;p,h_{k}\big{)}$ as follows

$\displaystyle f_{{Z}_{k}}\big{(}z;p,h_{k}\big{)}=p\delta_{\rm
c}(z)+(1-p)\delta_{\rm c}(z-h_{k}),$ (10)

where $\delta_{\rm c}(z)$, $z=z_{\rm r}+jz_{\rm I}\in{\mathbb{C}}$, is the
complex Delta function and is defined as

$\displaystyle\delta_{\rm c}(z)\triangleq\delta(z_{\rm r})\delta(z_{\rm I}),$
(11)

where $\delta(t)$, $t\in{\mathbb{R}}$, is the Dirac Delta function.

We consider that $g_{i}\triangleq\sum_{k=1}^{K}z_{k,i}\sim G$ and $z_{k,i}\sim
Z_{k}$ given in (10), where the symbol $\sim$ denotes distributed according
to. The PDF of the sum of independent random variables is obtained as the
convolution of the PDFs. For the complex random variable,
$G=\sum_{k=1}^{K}Z_{k}$, by employing the multi-binomial theorem [22], we can
obtain the PDF of $G$ as

$\displaystyle f_{G}$
$\displaystyle(g;p,{\bf{h}})=\sum_{v_{1}=0}^{1}\cdots\sum_{v_{K}=0}^{1}\bigg{[}p^{\sum_{k=1}^{K}v_{k}}(1-p)^{K-{\sum_{k=1}^{K}v_{k}}}$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\
\times\delta_{\rm c}\Big{(}g-\sum_{k=1}^{K}(1-v_{k})h_{k}\Big{)}\bigg{]},$
(12)

where ${\bf{h}}\triangleq[h_{1},h_{2},\cdots,h_{K}]^{T}$.

###### Proof.

See Appendix B. ∎

For the circularly symmetric complex Gaussian noise vector, ${\bf
w}\triangleq\big{[}w_{0}\ w_{1}\ \dots\ w_{N}\big{]}^{T}$, the PDF of the
random variable $W$ associated with the noise elements is expressed as

$\displaystyle f_{{W}}({w};\sigma_{\rm
w}^{2})\triangleq\frac{1}{{\pi\sigma_{\rm{w}}^{2}}}\exp{\bigg{(}\frac{-|w|^{2}}{\sigma_{\rm{w}}^{2}}\bigg{)}},$
(13)

where $w\in\mathbb{C}$. From (4), we have $y_{n}=g_{n}+w_{n}$,
$n=0,1,\dots,N$, where $g_{n}\sim G$ and $w_{n}\sim W$. Since $G$ and $W$ are
independent complex random variables, the PDF of $Y=G+W$ is obtained by the
linear convolution of the PDFs in (12) and (13), which results in

$\displaystyle
f_{Y}({{y}};p,{\bm{\beta}},{\bm{\theta}},\sigma_{\rm{w}}^{2})=\sum_{a=0}^{2^{K}-1}\frac{\xi_{a}}{{\pi\sigma_{\rm{w}}^{2}}}\mathcal{C}\mathcal{N}(y;\mu_{a},\sigma_{\rm{w}}^{2})$
(14)
$\displaystyle=\sum_{a=0}^{2^{K}-1}\frac{\xi_{a}}{{\pi\sigma_{\rm{w}}^{2}}}\exp{\Bigg{(}-\frac{|y-\mu_{a}|^{2}}{\sigma_{\rm{w}}^{2}}\Bigg{)}},$

where

$\displaystyle{\bm{\beta}}$ $\displaystyle\triangleq\Big{[}\beta_{1}\
\beta_{2}\ \ldots\ \beta_{K}\Big{]}^{T},$ (15a) $\displaystyle{\bm{\theta}}$
$\displaystyle\triangleq\big{[}\theta_{1}\ \theta_{2}\ \ldots\
\theta_{K}\big{]}^{T},$ (15b) $\displaystyle\xi_{a}$ $\displaystyle\triangleq
p^{\sum_{k=1}^{K}b_{k}}(1-p)^{K-{\sum_{k=1}^{K}b_{k}}},$ (16)

and

$\displaystyle{{\mu}}_{a}\triangleq\sum_{k=1}^{K}(1-b_{k})h_{k}=\sum_{k=1}^{K}(1-b_{k})\beta_{k}\exp(j\theta_{k}).$
(17)

with $b_{i}$ the $i$th bit in the binary representation of $a$ as

$\displaystyle a$
$\displaystyle=(b_{K},b_{K-1},\ldots,b_{1})_{2},\,\,\,\,\,\,\,\,\,\,\,\,\,\
b_{i}\in\\{0,1\\},$ (18)

and $a=0,1,\dots,2^{K}-1$.

As seen, $f_{Y}({{y}};p,{\bm{\beta}},{\bm{\theta}},\sigma_{\rm{w}}^{2})$
represents a 2D GM, where its modes are located at the delta functions given
in (12). Since $\theta_{1},\theta_{2},\ldots,\theta_{K}$, are independent
uniform random variables in the range of $[0,2\pi)$, and
$\beta_{1}>\beta_{2}>\ldots>\beta_{K}$, with probability of almost one, the
number of distinct mixtures is $2^{K}$. Fig. 5 illustrates the constellation
of the received signal for $K=3$ drones and maximum integer delay $M=20$ at
$20$ dB Signal-to-noise Ratio (SNR).

Figure 5: The constellation of the received signal for $K=3$ drones and
maximum integer delay $M=20$ at $20$ dB SNR. The values are multiplied by
$10^{6}$.

##  IV Joint Range and PO Estimation

The MLE for the vector parameters $[{\bm{\beta}}^{T}\ \bm{\theta}^{T}]^{T}$
given observation vector ${\bf{y}}=[{y}_{0}\ {y}_{2}\ \ldots\ {y}_{N}]^{T}$ is
expressed as

$\displaystyle\\{\hat{\bm{\beta}},\hat{\bm{\theta}}\\}=\operatorname*{arg\,max}_{{\bm{\beta}},{\bm{\theta}}}\sum_{n=0}^{N}\ln
f_{{Y}}({y}_{n};p,{\bm{\beta}},\bm{\theta},\sigma_{\rm{w}}^{2}).$ (19)

The maximization problem in (19) cannot be analytically solved in a trackable
manner. An alternative simple solution is to employ the EM algorithm to
estimate the $2^{K}$ modes of the GM components; then, we can decouple the
desired parameters, i.e., ${\bm{\beta}}$ and ${\bm{\theta}}$ from the
estimated modes.

Let $\bm{\mu}\triangleq[\mu_{0}\ \mu_{1}\ \ldots\ \mu_{2^{K}-1}]^{T}$ denote
the mode vector of the GM, where $\mu_{a}$, $a=0,1,\ldots,2^{K}-1$, is given
by (17). For $K=2$, we have

$\displaystyle\bm{\mu}\triangleq[\mu_{0}\ \mu_{1}\ \mu_{2}\ \mu_{3}]^{T}$ (20)
$\displaystyle=[\beta_{1}\exp(j\theta_{1})+\beta_{2}\exp(j\theta_{2})\,\,\
\beta_{2}\exp(j\theta_{2})\,\,\ \beta_{1}\exp(j\theta_{1})\,\,\ 0]^{T},$

where $\beta_{1}>\beta_{2}$, and for $K=3$, we can write

$\displaystyle\bm{\mu}$ $\displaystyle\triangleq[\mu_{0}\ \mu_{1}\ \mu_{2}\
\mu_{3}\ \mu_{4}\ \mu_{5}\ \mu_{6}\ \mu_{7}]^{T}$ (21)
$\displaystyle=\big{[}\beta_{1}\exp(j\theta_{1})+\beta_{2}\exp(j\theta_{2})+\beta_{3}\exp(j\theta_{3})$
$\displaystyle\,\,\,\,\,\,\
\beta_{2}\exp(j\theta_{2})+\beta_{3}\exp(j\theta_{3})\,\,\
\beta_{1}\exp(j\theta_{1})+\beta_{3}\exp(j\theta_{3})$
$\displaystyle\,\,\,\,\,\,\ \beta_{3}\exp(j\theta_{3})\,\,\
\beta_{1}\exp(j\theta_{1})+\beta_{2}\exp(j\theta_{2})\,\,\
\beta_{2}\exp(j\theta_{2})$ $\displaystyle\,\,\,\,\,\,\
\beta_{1}\exp(j\theta_{1})\,\,\ 0\big{]}^{T},$

with $\beta_{1}>\beta_{2}>\beta_{3}$. For $K=4$, we can write

$\displaystyle\bm{\mu}$ $\displaystyle\triangleq[\mu_{0}\ \mu_{1}\ \mu_{2}\
\mu_{3}\ \mu_{4}\ \mu_{5}\ \mu_{6}\ \mu_{7}\ \ldots\ \mu_{13}\ \mu_{14}\
\mu_{15}]^{T}$ (22)
$\displaystyle=\big{[}\beta_{1}\exp(j\theta_{1})+\beta_{2}\exp(j\theta_{2})+\beta_{3}\exp(j\theta_{3})+\beta_{4}\exp(j\theta_{4})$
$\displaystyle\,\,\,\,\,\,\
\beta_{2}\exp(j\theta_{2})+\beta_{3}\exp(j\theta_{3})+\beta_{4}\exp(j\theta_{4})$
$\displaystyle\,\,\,\,\,\,\
\beta_{1}\exp(j\theta_{1})+\beta_{3}\exp(j\theta_{3})+\beta_{4}\exp(j\theta_{4})$
$\displaystyle\,\,\,\,\,\,\
\beta_{3}\exp(j\theta_{3})+\beta_{4}\exp(j\theta_{4})$
$\displaystyle\,\,\,\,\,\,\
\beta_{1}\exp(j\theta_{1})+\beta_{2}\exp(j\theta_{2})+\beta_{4}\exp(j\theta_{4})$
$\displaystyle\,\,\,\,\,\,\
\beta_{2}\exp(j\theta_{2})+\beta_{4}\exp(j\theta_{4})$
$\displaystyle\,\,\,\,\,\,\
\beta_{1}\exp(j\theta_{1})+\beta_{4}\exp(j\theta_{4})\,\,\,\,\
\beta_{4}\exp(j\theta_{4})$ $\displaystyle\,\,\,\,\,\,\
\beta_{1}\exp(j\theta_{1})+\beta_{2}\exp(j\theta_{2})+\beta_{3}\exp(j\theta_{3})$
$\displaystyle\,\,\,\,\,\,\
\beta_{2}\exp(j\theta_{2})+\beta_{3}\exp(j\theta_{3})\,\,\,\,\
\beta_{1}\exp(j\theta_{1})+\beta_{3}\exp(j\theta_{3})$
$\displaystyle\,\,\,\,\,\,\ \beta_{3}\exp(j\theta_{3})\,\,\,\,\
\beta_{1}\exp(j\theta_{1})+\beta_{2}\exp(j\theta_{2})\,\,\,\,\
\beta_{2}\exp(j\theta_{2})$ $\displaystyle\,\,\,\,\,\,\
\beta_{1}\exp(j\theta_{1})\,\,\,\,\ 0\big{]}^{T},$

where $\beta_{1}>\beta_{2}>\beta_{3}>\beta_{4}$.

Let us define the discrete function $\chi_{q}(n)$ as

$\displaystyle\chi_{q}(n):\\{1,2,\ldots,l\\}\longrightarrow\\{1,2,\ldots,l\\},$
(23)

where for $n_{1}\neq n_{2}$, $\chi_{q}(n_{1})\neq\chi_{q}(n_{2})$. There are
$Q_{l}\triangleq l!$ unique functions in the form of (23), where $!$ denotes
the factorial function. Using (23), $Q_{l}$ permutation matrices of size
$l\times l$ can be defined as

$\displaystyle{{\bf\Lambda}_{q}}=\left[\begin{array}[]{l}{{\bf
e}_{{\chi_{q}}(1)}}\\\ {{\bf e}_{{\chi_{q}}(2)}}\\\ \,\,\,\,\vdots\\\ {{\bf
e}_{{\chi_{q}}({l})}}\end{array}\right],$ (28)

where ${\bf{e}}_{\ell}$, $\ell=1,2,\ldots,l$, denote the standard basis
vectors of length $l$ with a $1$ in the $\ell$th coordinate and 0’s elsewhere.
The set composed of all permutations of vector ${\bf a}=[a_{1}\ a_{2}\ \ldots\
a_{l}]^{T}$ is given by

$\displaystyle\mathcal{F}_{\bf
a}^{Q_{l}}\triangleq\Big{\\{}{{{\bf{\Lambda}}_{1}}{\bf{a}},{{\bf{\Lambda}}_{2}}{\bf{a}},\ldots,{{\bf{\Lambda}}_{Q_{l}}}{\bf{a}}}\Big{\\}},$
(29)

where $Q_{l}=l!$.

The EM algorithm estimates the permuted mode vector ${\bm{\eta}}=[\eta_{0}\
\eta_{1}\
\ldots\eta_{2^{K}-1}]^{T}\in\mathcal{F}_{\bm{\mu}}^{Q_{2^{K}}}\subset\mathbb{C}^{2^{K}}$,
where $Q_{2^{K}}=2^{K}!$ and ${\bm{\mu}}\triangleq[{{\mu}}_{0}\ {{\mu}}_{1}\
\ldots\ {{\mu}}_{2^{K}-1}]^{T}$. The EM algorithm defines a latent random
vector ${\bf u}\triangleq[u_{0}\ u_{1}\ \ldots\ u_{N}]^{T}$ that determines
the GM component from which the observation originates, i.e.,
$f_{{Y}|U}({{y}}_{n}|u_{n}=a;p,{\bm{\beta}},{\bm{\theta}},\sigma_{\rm{w}}^{2})\sim\mathcal{C}\mathcal{N}(y_{n};\mu_{a},\sigma_{\rm{w}}^{2})$,
where $P_{U}(u_{n}=a)=\xi_{a}$ for $n=0,1,\ldots,N$ and
$a=0,1,\ldots,2^{K}-1$. The EM algorithm iteratively maximizes the expected
value of the complete-data log-likelihood function to estimate the permuted
mode vector ${\bm{\eta}}=[\eta_{0}\ \eta_{1}\
\ldots\eta_{2^{K}-1}]^{T}\in\mathcal{F}_{\bm{\mu}}^{Q_{2^{K}}}\subset\mathbb{C}^{2^{K}}$
of the GM as follows [23, 24, 25]

$\displaystyle\hat{\bm{\eta}}^{(t+1)}=\operatorname*{arg\,max}_{\bm{\eta}}Q({\bm{\eta}}|{\bm{\eta}}^{(t)}),$
(30)

where ${\bm{\eta}}^{(0)}$ is the initialization vector,

$\displaystyle
Q({\bm{\eta}}|{\bm{\eta}}^{(t)})=\mathbb{E}_{{\bm{U}}|{\bm{Y}},{\bm{\eta}}^{(t)}}\Big{\\{}\ln
f_{{\bm{Y}},{\bm{U}}}({\bf y},{\bf
u};p,{\bm{\eta}},\sigma_{\rm{w}}^{2})\Big{\\}}$ (31)
$\displaystyle=\sum_{n=0}^{N}\sum_{a=0}^{2^{K}-1}\lambda_{a,n}^{(t)}\Bigg{(}\ln\frac{\xi_{a}}{\pi\sigma_{\rm
w}^{2}}-\frac{|y_{n}-\eta_{a}|^{2}}{\sigma_{\rm{w}}^{2}}\Bigg{)},$

with

$\displaystyle\lambda_{a,n}^{(t)}$
$\displaystyle=P_{U|Y}\big{(}u_{n}=a|{y}_{n};{\bm{\eta}}^{(t)}\big{)}$ (32)
$\displaystyle=\frac{P_{U}(u_{n}=a)f_{Y|U}\big{(}{y}_{n}|u_{n}=a;{\bm{\eta}}^{(t)}\big{)}}{f_{Y}(y_{n})}$
$\displaystyle=\frac{\xi_{a}\mathcal{C}\mathcal{N}\big{(}y_{n};\eta_{a}^{(t)},\sigma_{\rm{w}}^{2})}{\sum_{q=0}^{2^{K}-1}\xi_{q}\mathcal{C}\mathcal{N}(y_{n};\eta_{q}^{(t)},\sigma_{\rm{w}}^{2})},$

and the complete-data likelihood function is given by

$\displaystyle f_{{\bm{Y},\bm{U}}}({\bf y},{\bf
u};p,\bm{\eta},\sigma_{\rm{w}}^{2})=\prod_{n=0}^{N}\prod_{a=0}^{2^{K}-1}\big{(}\xi_{a}\mathcal{C}\mathcal{N}(y_{n};\eta_{a},\sigma_{\rm{w}}^{2})\big{)}^{\mathbb{I}\\{u_{n}=a\\}}.$
(33)

In (33), $\mathbb{I}\\{\cdot\\}$ denotes the indicator function, and $\xi_{a}$
is a function of $p$ and is given in (16). The EM algorithm at the $(t+1)$th
iteration estimates the vector ${\bm{\eta}}^{(t+1)}=[\eta_{0}^{(t+1)}\
\eta_{1}^{(t+1)}\ \ldots\eta_{2^{K}-1}^{(t+1)}]^{T}$ which is a permuted
version of the vector $\bm{\mu}$. The order of ${\bm{\eta}}^{(t+1)}$ depends
on the initialization of the EM algorithm, i.e., ${\bm{\eta}}^{(0)}$. By
solving the maximization problem in (30), the elements of
${\bm{\eta}}^{(t+1)}$ are iteratively updated as follows

$\displaystyle\eta_{a}^{(t+1)}=\frac{\sum_{n=0}^{N}\lambda_{a,n}^{(t)}y_{n}}{\sum_{n=0}^{N}\lambda_{a,n}^{(t)}},$
(34)

for $a=0,1,\dots,2^{K}-1$, where the convergence condition for the EM
algorithm is $\|{\bm{\eta}}^{(t+1)}-{\bm{\eta}}^{(t)}\|<\epsilon$, with
$\epsilon$ a preset threshold. We denote $\hat{\bm{\eta}}={\bm{\eta}}^{(t+1)}$
when the EM algorithm converges at the $(t+1)$th iteration.

The EM algorithm may converge to a local maximum of the observed data
likelihood function, depending on the starting value. A variety of heuristic
or metaheuristic approaches exist to escape a local maximum, such as random-
restart hill climbing where the EM algorithm starts with different random
initial estimates. The recursive expression in (34) is very straightforward
and can be easily computed for multiple initial points. Hence, EM with random-
restart hill climbing can be employed for the problem at hand. As mentioned
earlier, there is an order ambiguity in the vector $\hat{\bm{\eta}}$; thus,
reordering estimation is required to resolve this ambiguity. In the next
section, we propose different solutions based on permutation- based
combinatorial optimization for one to one mapping between the elements of
$\hat{\bm{\eta}}$ and ${\bm{\mu}}$.

##  V Reordering Estimation

The goal of reordering estimation is to change the order of the elements in
$\hat{\bm{\eta}}\triangleq[\hat{{\eta}}_{0}\ \hat{{\eta}}_{1}\ \ldots\
\hat{{\eta}}_{2^{K}-1}]^{T}$, obtained by the EM algorithm, to achieve a new
vector $\hat{\bm{\mu}}\triangleq[\hat{{\mu}}_{0}\ \hat{{\mu}}_{1}\ \ldots\
\hat{{\mu}}_{2^{K}-1}]^{T}$ that corresponds to an estimate of
${\bm{\mu}}\triangleq[{{\mu}}_{0}\ {{\mu}}_{1}\ \ldots\
{{\mu}}_{2^{K}-1}]^{T}$. Let us denote $\hat{{\eta}}_{i}$ as the element of
$\hat{\bm{\eta}}$ that corresponds to ${\mu}_{2^{K}-1}=0$. The index $i$ can
be estimated as

$\displaystyle|\hat{\eta}_{i}|<\min\Big{\\{}|\hat{\eta}_{0}|,|\hat{\eta}_{1}|,\ldots,|\hat{\eta}_{i-1}|,|\hat{\eta}_{i+1}|,\ldots,|\hat{\eta}_{2^{K}-1}|\Big{\\}}.$
(35)

Accordingly, we have $\hat{{\mu}}_{2^{K}-1}=\hat{\eta}_{i}$. Let us now define

$\displaystyle{\hat{\bm{\eta}}}_{i}\triangleq[\hat{{\eta}}_{0}\
\hat{{\eta}}_{1}\ \ldots\ \hat{\eta}_{i-1}\ {\hat{\eta}}_{i+1}\ \ldots\
\hat{{\eta}}_{2^{K}-1}]^{T},$ (36)

and

$\displaystyle{\hat{\mathcal{A}}}_{l}$
$\displaystyle\triangleq\Big{\\{}[\phi_{1}\ \phi_{2}\ \ldots\
\phi_{l}]^{T}\big{|}\forall{d}\in\\{1,2,\ldots,l\\},$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\
\phi_{d}\in\\{\hat{{\eta}}_{0},\hat{{\eta}}_{1},\ldots\
\hat{\eta}_{i-1},{\hat{\eta}}_{i+1}\ \ldots\ \hat{{\eta}}_{2^{K}-1}\\},$ (37)
$\displaystyle\,\,\,\,\,\,\,\,\,\,\
|\phi_{1}|>|\phi_{2}|>\ldots|\phi_{l}|\Big{\\}}.$

Different reordering estimation methods for $\hat{\bm{\beta}}$ and
$\hat{\bm{\theta}}$ can be considered. The Least Squares (LS) reordering
estimation for $\hat{\bm{\beta}}$ and $\hat{\bm{\theta}}$ form
$\hat{\bm{\eta}}_{i}$ is given by

$\displaystyle\hat{\bf
h}=\hat{\bm{\beta}}e^{j\hat{\bm{\theta}}}=\hat{\bf\Lambda}{\bf
A}\hat{\bf\Phi},$ (38)

where

$\displaystyle\\{\hat{\bf\Lambda},\hat{{\bf\Phi}}\\}$
$\displaystyle=\operatorname*{arg\,min}_{{\bf\Lambda},{{\bf\Phi}}}\big{\|}{\bf\Lambda}{\bf
A}{\bf\Phi}-\hat{\bm{\eta}}_{i}\big{\|}_{2},$ (39)
$\displaystyle\textrm{s.t.}\,\
{\bf\Phi}\triangleq[\Phi_{1},\Phi_{2},\ldots,\Phi_{K}]^{T}\in{\hat{\mathcal{A}}}_{K}$
$\displaystyle\,\,\,\,\,\,\,\
{\bf\Lambda}\in\\{{\bf\Lambda}_{1},{\bf\Lambda}_{2},\ldots,{\bf\Lambda}_{(2^{K-1})!}\\}$

where ${\bf\Lambda}_{i}$ is a permutation matrix of size $2^{K-1}\times
2^{K-1}$ given in (28) for $l=2^{K-1}$, and ${\bf A}$ is the $2^{K-1}\times K$
matrix given by

$\displaystyle{\bf
A}\triangleq\left[\begin{array}[]{l}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\
{{\bf{e}}_{1}}\\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\ {{\bf{e}}_{2}}\\\
\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\ \vdots\\\
\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\ {{\bf{e}}_{K}}\\\
\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{\bf{e}}_{1}}+{{\bf{e}}_{2}}\\\
\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{\bf{e}}_{1}}+{{\bf{e}}_{3}}\\\
\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\vdots\\\
{{\bf{e}}_{1}}+{{\bf{e}}_{2}}+\ldots+{{\bf{e}}_{K}}\end{array}\right],$ (48)

where ${\bf{e}}_{j}$, $j=1,2,\ldots,K$, denote the standard basis vectors of
length $K$.

An alternative minimization problem with lower computational complexity for
vector reordering is given by

$\displaystyle\\{\hat{\bf\Lambda},\hat{{\bf\Phi}}\\}$
$\displaystyle=\operatorname*{arg\,min}_{{\bf\Lambda},{{\bf\Phi}}}\big{\|}{\bf\Lambda}{\bf
A}{\bf\Phi}-\hat{\bm{\eta}}_{i}\big{\|}_{2},$ (49)
$\displaystyle\textrm{s.t.}\,\
{\bf\Phi}\triangleq[\Phi_{1},\Phi_{2},\ldots,\Phi_{K}]^{T}\in\mathbb{C}^{K}$
$\displaystyle\,\,\,\,\,\,\,\
{\bf\Lambda}\in\\{{\bf\Lambda}_{1},{\bf\Lambda}_{2},\ldots,{\bf\Lambda}_{(2^{K-1})!}\\}$

By employing the solution of the unconstrained LS minimization for a linear
observation model [26] and the fact that
${\bf\Lambda}^{T}{\bf\Lambda}={\bf\Lambda}{\bf\Lambda}^{T}={\bf I}_{2^{K}-1}$,
we can easily write

$\displaystyle{{\bf\Phi}}=({\bf A}^{T}{\bf A})^{-1}{\bf
A}^{T}{\bf\Lambda}^{T}\hat{\bm{\eta}}_{i},$ (50)

and thus we can formulate the following problem

$\displaystyle\hat{\bf\Lambda}$
$\displaystyle=\operatorname*{arg\,min}_{{\bf\Lambda}}\hat{\bm{\eta}}_{i}^{T}(\hat{\bm{\eta}}_{i}-{\bf\Lambda}{\bf
A}{{\bf\Phi}})$ (51) $\displaystyle\textrm{s.t.}\,\ {{\bf\Phi}}=({\bf
A}^{T}{\bf A})^{-1}{\bf A}^{T}{\bf\Lambda}^{T}\hat{\bm{\eta}}_{i}$
$\displaystyle\,\,\,\,\,\,\,\ |\Phi_{1}|>|\Phi_{2}|>\ldots>|\Phi_{K}|$
$\displaystyle\,\,\,\,\,\,\,\
{\bf\Lambda}\in\\{{\bf\Lambda}_{1},{\bf\Lambda}_{2},\ldots,{\bf\Lambda}_{(2^{K-1})!}\\}.$

Our simulation experiments show that for low and moderate SNRs, the vector
reordering based on the combinatorial optimization in (49) outperforms the one
in (39) in terms of estimation error. While the vector reordering based on the
minimization formulation in (49) offers a lower computational complexity
compared to the minimization in (39), the computational complexity of both
methods is still high for $K>4$.

An alternative low complexity solution is to define a combinatorial
optimization with lower cardinality. By selecting all $Q$-combination of the
set $\\{\hat{{\eta}}_{0}\ \hat{{\eta}}_{1}\ \ldots\ \hat{\eta}_{i-1}\
{\hat{\eta}}_{i+1}\ \ldots\ \hat{{\eta}}_{2^{K}-1}\\}$, we can define a
$C_{Q}^{2^{K}-1}$-cardinality combinatorial optimization problem that
minimizes a linear/non-linear combinations of the elements in the set in such
a way that one or more elements of the set can be unambiguously assigned to
the elements of $\bm{\mu}$.

As an example, for $K=4$ drones/aircrafts, our goal is to estimate
${\mu}_{14}=\beta_{1}\exp(j\theta_{1})$,
${\mu}_{13}=\beta_{2}\exp(j\theta_{2})$,
${\mu}_{11}=\beta_{3}\exp(j\theta_{3})$, and
${\mu}_{7}=\beta_{4}\exp(j\theta_{4})$, where
$\beta_{1}>\beta_{2}>\beta_{3}>\beta_{4}$ (refer to (22)). We can easily show
that there is a unique solution for the linear equation

$\displaystyle z_{0}+z_{1}+z_{2}+z_{3}-z_{4}=0,$ (52)

where $z_{i}\in\\{\mu_{0}\ \mu_{1}\ \ldots\ \mu_{14}\\}$ and
$|z_{1}|>|z_{2}|>|z_{3}|>|z_{4}|$. This unique solution is
${\mu}_{14}+{\mu}_{13}+{\mu}_{11}+{\mu}_{7}-{\mu}_{0}=0$, where
${\mu}_{0}=\sum_{i=1}^{4}\beta_{i}\exp(j\theta_{i})$,
${\mu}_{14}=\beta_{1}\exp(j\theta_{1})$,
${\mu}_{13}=\beta_{2}\exp(j\theta_{2})$,
${\mu}_{11}=\beta_{3}\exp(j\theta_{3})$, and
${\mu}_{7}=\beta_{4}\exp(j\theta_{4})$ for
$\beta_{1}>\beta_{2}>\beta_{3}>\beta_{4}$. By taking this into account, we can
define the following combinatorial optimization for reordering estimation

$\displaystyle\hat{\bm{\Phi}}$
$\displaystyle=\operatorname*{arg\,min}_{\bm{\Phi}}\ \big{|}{{\bf
v}^{T}}{\bm{\Phi}}\big{|}$ (53) $\displaystyle\textrm{s.t.}\,\
{\bm{\Phi}}\triangleq[\phi_{1}\ \phi_{2}\ \ldots\
\phi_{5}]^{T}\in\hat{\mathcal{A}}_{5}$

where $\hat{\bm{\Phi}}\triangleq[\hat{\phi}_{1}\ \hat{\phi_{2}}\ \ldots\
\hat{\phi_{5}}]^{T}$, the set $\hat{\mathcal{A}}_{5}$ is defined in (37) for
$l=5$, and ${\bf v}\triangleq[v_{1}\ v_{2}\ v_{3}\ v_{4}\ v_{5}]^{T}=[1\ 1\ 1\
1\ -1]^{T}$. It is obvious that (53) represents a
$C_{5}^{15}=3003$-cardinality combinatorial optimization problem. By solving
the combinatorial optimization in (53), we obtain

$\displaystyle\hat{\bf
h}=\hat{\bm{\beta}}e^{j\hat{\bm{\theta}}}=[\hat{\phi}_{1}\ \hat{\phi_{2}}\
\hat{\phi}_{3}\ \hat{\phi_{4}}]^{T}.$ (54)

Different combinatorial optimization problems can be defined for reordering
estimation. For $K=4$, let us consider the following linear equation

$\displaystyle 7\sum_{i=0}^{3}z_{i}-\sum_{i=4}^{14}z_{i}=0,$ (55)

where $z_{i}\in\\{\mu_{0}\ \mu_{1}\ \ldots\ \mu_{14}\\}$ and
$|z_{1}|>|z_{2}|>|z_{3}|>|z_{4}|$. Since
$\beta_{1}>\beta_{2}>\beta_{3}>\beta_{4}$ (refer to assumption 1), we can show
that the solution of (55) is given by
$7\mu_{14}+7\mu_{13}+7\mu_{11}+7\mu_{7}-{\bf c}^{T}{\bf\Lambda}^{T}=0$, where
${\bf c}\triangleq[\mu_{0}\ \mu_{1}\ \mu_{2}\ \mu_{3}\ \mu_{4}\ \mu_{5}\
\mu_{6}\ \mu_{8}\ \mu_{9}\ \mu_{10}\ \mu_{12}]^{T}$ and ${\bf\Lambda}$ is a
$11\times 11$ permutation matrix.333The second summation is independent of the
permutation of $\bf c$. Taking this equation into account, the combinatorial
optimization for reordering estimation can be expressed as

$\displaystyle\hat{\bm{\Phi}}$
$\displaystyle=\operatorname*{arg\,min}_{{\bm{\Phi}}}\
\Big{|}7{\bm{\Phi}}^{T}{\bf 1}_{4}-\alpha\Big{|},$ (56)
$\displaystyle\textrm{s.t.}\,\ {\bm{\Phi}}\triangleq[\phi_{1}\ \phi_{2}\
\phi_{3}\ \phi_{4}]^{T}\in\hat{\mathcal{A}}_{4}$ $\displaystyle\,\,\,\,\,\,\,\
\alpha=\sum_{\kappa\in{\mathcal{S}}}\kappa$

where ${\mathcal{S}}\triangleq\\{\hat{{\eta}}_{0},\hat{{\eta}}_{1},\ldots\
\hat{\eta}_{i-1},{\hat{\eta}}_{i+1}\ \ldots\
\hat{{\eta}}_{2^{K}-1}\\}-\\{\phi_{1},\phi_{2},\phi_{3},\phi_{4}\\}$,
$\hat{\bm{\Phi}}\triangleq[\hat{\phi}_{1}\ \hat{\phi}_{2}\ \hat{\phi}_{3}\
\hat{\phi}_{4}]^{T}$, ${\bf 1}_{4}=[1\ 1\ 1\ 1]^{T}$, and the set
$\hat{\mathcal{A}}_{4}$ is defined in (37) for $l=4$. By using the solution of
the minimization in (56), we can write $\hat{\bf
h}=\hat{\bm{\beta}}e^{j\hat{\bm{\theta}}}=[\hat{\phi}_{1}\ \hat{\phi_{2}}\
\hat{\phi}_{3}\ \hat{\phi_{4}}]^{T}$.

It should be mentioned that other combinatorial optimizations can also be
defined for reordering estimation. Ambiguity removal through multiple
combinatorial optimization problems can also be considered.

###  V-A Joint Range and PO Estimation

For joint range and PO estimation of $K$ drones, the modes
$\mu_{a_{k}}\triangleq\sqrt{P_{k}L_{k}}e^{j\theta_{k}}$, $k=1,2,\ldots,K$, are
needed to be estimated, where

$\displaystyle a_{k}=\sum_{\begin{subarray}{c}n=0\\\ n\neq
k-1\end{subarray}}^{K-1}2^{n}.$ (57)

Let $\hat{\mu}_{a_{k}}$, $k=1,2,\ldots,K$, denote the estimated modes after
reordering estimation. By using (2), the range and PO for the $k$th drone are
estimated as

$\displaystyle{\hat{r}}_{k}=\frac{\lambda_{\rm
c}\sqrt{P_{k}}}{\big{|}4\pi\hat{\mu}_{a_{k}}\big{|}},$ (58)

and

$\displaystyle\hat{\theta}_{k}=\begin{cases}\tan^{-1}\frac{\Im\\{\hat{\mu}_{a_{k}}\\}}{\Re\\{\hat{\mu}_{a_{k}}\\}},\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\
{\Re\\{\hat{\mu}_{a_{k}}\\}}\geq 0\\\
\tan^{-1}\frac{\Im\\{\hat{\mu}_{a_{k}}\\}}{\Re\\{\hat{\mu}_{a_{k}}\\}}+\pi,\,\,\,\,\,\,\,\,\
{\Re\\{\hat{\mu}_{a_{k}}\\}}<0\end{cases}.$ (59)

where $\Im\\{\cdot\\}$ and $\Re\\{\cdot\\}$ denotes the real and imaginary
operators, respectively. The proposed EM-based joint ranging and PO estimation
is summarized in Algorithm 1, where $T_{\rm EM}$ denotes the number of
iterations of the EM algorithm.444The pseudocode has been written for the
reordering methods in (39) and (49). Similar pseudocode can be written for
other reordering estimation methods.

Algorithm 1 : EM-based joint ranging and PO estimation

1:Input: $P_{1}>P_{2}>\ldots>P_{K}$ and $\bf y$

2:Output: $\hat{r}_{k}$ and $\hat{\theta}_{k}$ for $k=1,2,\ldots,K$

3:Initialize ${\bm{\eta}}^{(0)}$

4:for $t=0,1,\ldots,T_{\rm EM}-1$ do

5: Use (34) to obtain ${\bm{\eta}}^{(t+1)}\triangleq[{{\eta}}_{0}^{(t+1)}\
{{\eta}}_{1}^{(t+1)}\ \ldots\ {{\eta}}_{2^{K}-1}^{(t+1)}]^{T}$

6: if $\|{\bm{\eta}}^{(t+1)}-{\bm{\eta}}^{(t)}\|<\epsilon$ or $t=T_{\rm EM}-1$

7: $\hat{\bm{\eta}}={\bm{\eta}}^{(t+1)}$ and $t=T_{\rm EM}$

8: end if

9:end for

10:Use $\eqref{first_elm}$ to estimate the element of $\hat{\bm{\eta}}$ that
corresponds to ${\mu}_{2^{K}-1}=0$ and obtain ${\hat{\bm{\eta}}}_{i}$ in (36).

11:Solve the combinatorial optimization in (39) or (49) to obtain
$\\{\hat{\bf\Lambda},\hat{{\bf\Phi}}\\}$

12:Estimate $\hat{\bf h}=[\hat{\mu}_{a_{1}}\ \hat{\mu}_{a_{2}}\ldots\
\hat{\mu}_{a_{K}}]$ by using (38) for $a_{k}$, $k=1,2,\ldots,K$, given in
(57).

13:Estimate $r_{k}$ and $\theta_{k}$ by employing (58) and (59).

##  VI Multiple Antennas Receiver

In this section, we extend the derived maximum likelihood cost function and
the proposed EM-based joint range and PO estimators to the case of multiple
receive antennas.

###  VI-A Maximum Likelihood and EM Estimators

We consider that the receiver is equipped with $N_{\rm r}$ receive antennas.
We assume spatially uncorrelated Gaussian noise. With the assumption that the
time delay between receive antennas is negligible, and the path loss between
the $k$th drone and $\ell$th receive antenna is the same for all receive
antennas, i.e, $L_{\ell,k}=L_{k}$, $k=1,2,\ldots,K$, $\ell=1,2,\ldots,N_{\rm
r}$, the received complex baseband signal at the multiple-receive antennas can
be expressed as

$\displaystyle{\bf Y}={\bf H}{\bf X}+{\bf W},$ (60)

where ${\bf X}\triangleq[{\bf x}_{1}\ {\bf x}_{2}\ \ldots\ {\bf
x}_{K}]^{T}\in\mathbb{C}^{K\times(N+1)}$, ${\bf Y}\triangleq[{\bf y}_{0}\ {\bf
y}_{1}\ldots\ {\bf y}_{N}]\in\mathbb{C}^{N_{\rm{r}}\times(N+1)}$, ${\bf
x}_{k}$ is given by (6), and ${\bf y}_{n}\triangleq[y_{1,n}\ y_{2,n}\ \ldots\
y_{N_{\rm r},n}]^{T}$ denotes the received vector at time index $n$. In (60),
the matrices ${\bf{H}}\in\mathbb{C}^{N_{\rm{r}}\times K}$ and
${\bf{W}}\in\mathbb{C}^{N_{\rm{r}}\times(N+1)}$ are given as

$\displaystyle{\bf{H}}=\left[\begin{array}[]{l}{\bf{h}}_{1}^{T}\\\
{\bf{h}}_{2}^{T}\\\ \,\,\vdots\\\
{\bf{h}}_{{N_{\rm{r}}}}^{T}\end{array}\right],\,\,\,\,\,\,\,\,\,\
{\bf{W}}=\left[\begin{array}[]{l}{\bf{w}}_{1}^{T}\\\ {\bf{w}}_{2}^{T}\\\
\,\,\vdots\\\ {\bf{w}}_{{N_{\rm{r}}}}^{T}\end{array}\right],$ (69)

where

$\displaystyle{\bf{h}}_{\ell}$ $\displaystyle\triangleq\big{[}h_{\ell,1}\
h_{\ell,2}\ \ldots\ h_{\ell,K}\big{]}^{T}$ (70)
$\displaystyle=\big{[}\beta_{1}\exp(j\theta_{\ell,1})\ \ldots\
\beta_{N}\exp(j\theta_{\ell,K})\big{]}^{T},$

$\beta_{k}=\beta_{\ell,k}=\sqrt{P_{k}L_{k}}$, ${\bf
w}_{\ell}\triangleq[w_{\ell,0}\ w_{\ell,1}\ \ldots\ w_{\ell,N}]^{T}$ with
$w_{\ell,n}\sim\mathcal{C}\mathcal{N}(0,\sigma_{\rm{w}}^{2})$ the complex
Gaussian noise at the $\ell$th receive antenna at time index $n$. As seen,
while $\beta_{1,k}=\beta_{2,k}=\cdots=\beta_{N_{\rm{r}},k}=\beta_{k}$, the
phases $\theta_{1,k},\theta_{2,k},\ldots,\theta_{N_{\rm{r}},k}$ are
independent random values in $[0\ 2\pi)$.

The joint PDF of the elements of ${\bf Y}$ is given by

$\displaystyle f_{\bf
Y}({\bf{Y}};p,{\bm{\beta}},{\bm{\Theta}},\sigma_{\rm{w}}^{2})=\prod_{\ell=1}^{N_{\rm
r}}\prod_{n=0}^{N}f_{Y}({{y}_{\ell,n}};p,{\bm{\beta}},{\bm{\theta}}_{\ell},\sigma_{\rm{w}}^{2}),$
(71)

where ${\bm{\beta}}\triangleq[\beta_{1}\ \beta_{2}\ \ldots\ \beta_{K}]^{T}$,
${\bm{\Theta}}\triangleq\big{[}{\bm{\theta}}_{1}^{T}\ {\bm{\theta}}_{2}^{T}\
\ldots\ {\bm{\theta}}_{N_{\rm r}}^{T}\big{]}^{T}$,
${\bm{\theta}}_{\ell}\triangleq\big{[}\theta_{\ell,1}\ \theta_{\ell,2}\
\ldots\ \theta_{\ell,K}\big{]}^{T}$, and
$f_{Y}({{y}};p,{\bm{\beta}},{\bm{\theta}},\sigma_{\rm{w}}^{2})$ is given in
(14). We can easily write the direct MLE of the parameter vector $\bm{\beta}$
and ${\bm{\Theta}}$ as

$\displaystyle\\{\hat{\bm{\beta}},\hat{\bm{\Theta}}\\}=\operatorname*{arg\,max}_{{\bm{\beta}},{\bm{\Theta}}}\sum_{\ell=1}^{N_{\rm
r}}\sum_{n=0}^{N}\ln
f_{{Y}}\big{(}{y}_{\ell,n};p,{\bm{\beta}},\bm{\theta}_{\ell},\sigma_{\rm{w}}^{2}\big{)}.$
(72)

Solving the maximization problem in (72) is not straightforward. Hence, the
direct MLE cannot be obtained. However, indirect estimation can be obtained by
estimating the modes of the GM similar to the single receive antenna.

Analogous to the single receive antenna scenario, we can employ the EM
algorithm for estimating the modes of the GM. For multiple receive antennas,
the EM algorithm defines an identical latent random vector ${\bf
u}\triangleq[u_{0}\ u_{1}\ \ldots\ u_{N}]^{T}$ for all receive antennas. This
random vector determines the GM component from which the observation
originates, i.e.,

$\displaystyle
f_{Y}({{y}_{\ell,n}}|u_{n}=a;p,{\bm{\beta}},{\bm{\theta}}_{\ell},\sigma_{\rm{w}}^{2})\sim\mathcal{C}\mathcal{N}(y_{\ell,n};\mu_{\ell,a},\sigma_{\rm{w}}^{2}),$
(73)

where $n=0,1,\ldots,N$, $\ell=1,2,\ldots,{N_{\rm{r}}}$, $a=0,1,\dots,2^{K}-1$,
${{\mu}}_{\ell,a}\triangleq\sum_{k=1}^{K}(1-b_{k})h_{\ell,k}=\sum_{k=1}^{K}(1-b_{k})\beta_{k}\exp(j\theta_{\ell,k})$,
$a=(b_{K},b_{K-1},\ldots,b_{1})_{2},b_{i}\in\\{0,1\\}$, and
$P_{U}(u_{n}=a)=\xi_{a}$.

The EM algorithm iteratively maximizes the expected value of the complete-data
log-likelihood function to estimate the vector
${\bm{\Gamma}}\triangleq[{\bm{\eta}}_{1}^{T}\ {\bm{\eta}}_{2}^{T}\ \ldots\
{\bm{\eta}}_{N_{\rm r}}^{T}]^{T}$ as

$\displaystyle{\bm{\Gamma}}^{(t+1)}=\operatorname*{arg\,max}_{{\bm{\Gamma}}}C\big{(}{\bm{\Gamma}}|{\bm{\Gamma}}^{(t)}\big{)},$
(74)

where the vector ${\bm{\eta}}_{\ell}\triangleq[\eta_{\ell,0}\ \eta_{\ell,1}\
\ldots\ \eta_{\ell,2^{K}-1}]^{T}$ denotes the $2^{K}$ modes of the 2D GM at
the $l$th receive antenna, ${\bm{\Gamma}}^{(0)}$ is the initialization vector,

$\displaystyle
C({\bm{\Gamma}}|{\bm{\Gamma}}^{(t)})=\mathbb{E}_{{\bm{U}}|{\bm{Y}};{\bm{\Gamma}}^{(t)}}\Big{\\{}\ln
f_{{\bm{Y}},{\bm{U}}}({\bf Y},{\bf
u};p,{\bm{\Gamma}},\sigma_{\rm{w}}^{2})\Big{\\}}$ (75)
$\displaystyle=\sum_{\ell=1}^{N_{\rm
r}}\sum_{n=0}^{N}\sum_{a=0}^{2^{K}-1}\delta_{a,n}^{(t)}\Bigg{(}\ln\frac{\xi_{a}}{\pi\sigma_{\rm
w}^{2}}-\frac{|y_{\ell,n}-\eta_{\ell,n}|^{2}}{\sigma_{\rm{w}}^{2}}\Bigg{)},$
$\displaystyle\delta_{a,n}^{(t)}$
$\displaystyle=P_{U|{{\bm{Y}}}}\big{(}u_{n}=a|{\bf{y}}_{n};p,{\bm{\Gamma}}^{(t)},\sigma_{\rm
w}^{2}\big{)}$ (76)
$\displaystyle=\frac{P_{U}(u_{n}=a)f_{{\bm{Y}}|{U}}\big{(}{\bf{y}}_{n}|u_{n}=a;p,{\bm{\Gamma}}^{(t)},\sigma_{\rm{w}}^{2}\big{)}}{\sum_{a=0}^{2^{K}-1}f_{\bm{Y},U}({\bf
y}_{n},u_{n}=a;p,{\bm{\Gamma}}^{(t)},\sigma_{\rm w}^{2})}$
$\displaystyle=\frac{\xi_{a}\prod_{\ell=1}^{N_{\rm{r}}}\mathcal{C}\mathcal{N}\big{(}y_{\ell,n};\eta_{\ell,a}^{(t)},\sigma_{\rm{w}}^{2})}{\sum_{a=0}^{2^{K}-1}\xi_{a}\prod_{\ell=1}^{N_{\rm
r}}\mathcal{C}\mathcal{N}(y_{\ell,n};\eta_{\ell,a}^{(t)},\sigma_{\rm{w}}^{2})},$

and the joint complete-data likelihood function is given by

$\displaystyle f_{{\bm{Y}},{\bm{U}}}({\bf Y},$ $\displaystyle{\bf
u};p,\bm{\Gamma},\sigma_{\rm{w}}^{2})$ (77)
$\displaystyle=\prod_{\ell=1}^{N_{\rm{r}}}\prod_{n=0}^{N}\prod_{a=0}^{2^{K}-1}\big{(}\xi_{a}\mathcal{C}\mathcal{N}(y_{\ell,n};\mu_{\ell,a},\sigma_{\rm{w}}^{2})\big{)}^{\mathbb{I}\\{u_{n}=a\\}}.$

By solving the maximization problem in (74), the elements of
${\bm{\eta}}^{(t+1)}$ are iteratively updated as follows

$\displaystyle\eta_{\ell,a}^{(t+1)}=\frac{\sum_{n=0}^{N}\delta_{a,n}^{(t)}y_{\ell,n}}{\sum_{n=0}^{N}\delta_{a,n}^{(t)}}.$
(78)

The convergence condition for the EM algorithm is
$\|{\bm{\Upsilon}}^{(t+1)}-{\bm{\Upsilon}}^{(t)}\|<\epsilon$, where $\epsilon$
is a preset threshold. We denote $\hat{\eta}_{\ell,a}={\eta}_{\ell,a}^{(t+1)}$
when the EM algorithm converges at the $(t+1)$th iteration.

While the EM algorithm estimates $N_{\rm r}2^{K}$ parameters, only $N_{\rm
r}K$ parameters are used for joint ranging and PO estimation of $K$ drones.
These $N_{\rm r}K$ modes are
${\mu}_{\ell,{a_{1}}},{\mu}_{\ell,a_{2}},\ldots,{\mu}_{\ell,a_{K}}$, where
$a_{k}$ is defined in (57). After the EM convergence, each receive antenna
independently applies estimation mapping as explained in section V. Let
$\hat{\mu}_{\ell,{a_{1}}},\hat{\mu}_{\ell,a_{2}},\ldots,\hat{\mu}_{\ell,a_{K}}$
denote the estimated and reordered modes at the $\ell$th receive antenna. By
averaging, then we can write

$\displaystyle|\hat{\mu}_{{a_{k}}}|=\frac{1}{N_{\rm r}}\sum_{\ell=1}^{N_{\rm
r}}|\hat{\mu}_{\ell,{a_{k}}}|\propto\frac{1}{\hat{r}_{k}},$ (79)

where $k=1,2,\ldots,K$. By substituting (79) into (58), we can estimate the
range of the $k$th drone. The PO for each drone-receive antenna is obtained by
replacing
$\hat{\mu}_{\ell,{a_{1}}},\hat{\mu}_{\ell,a_{2}},\ldots,\hat{\mu}_{\ell,a_{K}}$,
$\ell=1,2,\cdots,N_{\rm r}$ into (59).

It should be mentioned that the EM algorithm can also be implemented
independently at each receive antennas. Then, after reordering estimation, the
ranges are obtained by averaging. This results in a lower complexity solution.

###  VI-B Outlier Detection

For multiple receive antennas at the receiver, it is possible that one or
multiple outliers appear in the estimated sequence
$|\hat{\mu}_{\ell,{a_{1}}}|,|\hat{\mu}_{\ell,a_{1}}|,\ldots,|\hat{\mu}_{\ell,a_{K}}|$,
where $a_{k}$ is given in (57). These outliers can increase the estimation
error of $|\hat{\mu}_{{a_{k}}}|$ in (79) that is used for range estimation. To
remove the effect of outliers, we can use outlier detection algorithms, such
as the median absolute deviation (MAD) [27]. In this case, the averaging is
performed over the receive antennas without outliers. This results in a
significant performance improvement in range estimation.

##  VII Simulation

In this section, we examine the performance of the proposed EM-based joint
range and PO estimator through several simulation experiments to confirm the
effectiveness of the proposed algorithm.

###  VII-A Simulation Setup

Unless otherwise mentioned, we considered the following scenarios:

1) $K=3$ drones with ranges $r_{1}\in\mathcal{U}_{\rm c}[500,1000]$,
$r_{2}\in\mathcal{U}_{\rm c}[1500,2000]$, and $r_{3}\in\mathcal{U}_{\rm
c}[2500,3000]$.

2) $K=2$ drones with ranges $r_{1}\in\mathcal{U}_{\rm c}[500,1500]$ and
$r_{2}\in\mathcal{U}_{\rm c}[2000,3000]$.

3) $K=1$ drones with range $r_{1}\in\mathcal{U}_{\rm c}[500,3000]$.

Figure 6: Range tracking by the proposed EM-based estimator for $K=3$ drones
at $\gamma=20$ dB SNR. The dashed and solid lines denote the true and
estimated range, respectively.

We considered free space path loss model in (2) and $\lambda_{\rm{c}}\approx
0.2752$ m for ADS-B systems. The number of receive antennas was set to $N_{\rm
r}=5$. The discrete time delay of each drone for each ADS-B packet
transmission is assumed to be modeled by a discrete i.i.d. random variable
with uniform distribution as $m_{1},m_{2},\ldots,m_{K}\in\mathcal{U}_{\rm
d}[0,M]$, where $M=20$. We also assumed that the PO between the $k$th drone
and the $\ell$th receive antenna is modeled by a continuous i.i.d. random
variable with uniform distribution as $\theta_{\ell,k}\in\mathcal{U}_{\rm
c}[0,2\pi)$, $k=1,2,\ldots,K$, $\ell=1,2,\ldots,N_{\rm r}$. The average SNR at
each receive antenna, $\gamma$, is defined as

$\displaystyle\gamma\triangleq\frac{\sum_{k=1}^{K}L_{k}^{\rm
av}P_{k}}{\sigma_{\rm w}^{2}},\,\,\,\,\,\,\,\,\,\ L_{k}^{\rm
av}\triangleq\Big{(}\frac{\lambda_{\rm{c}}}{4\pi(\mathbb{E}\\{r_{k}\\})^{2}}\Big{)}^{2}.$

The noise power in dBm is considered $\sigma_{\rm w}^{2}=-174+\log_{10}(B)$
for each receive antennas. For $B=4$ MHz, the noise power in dBm is -104.58.
Without loss of generality, we considered that $P_{1}=P_{2}=\dots=P_{K}$. The
performance of the proposed EM-based joint range and PO estimator was
evaluated in terms of $1-P_{{\rm{out}},{r}}$ and
$1-P_{{\rm{out},{\rm\theta}}}$ where

$\displaystyle
P_{{\rm{out}},r}\triangleq\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}\bigg{\\{}\frac{|\hat{r}_{k}-r_{k}|}{r_{k}}>\alpha_{r}\bigg{\\}},$

and

$\displaystyle
P_{{\rm{out}},\theta}\triangleq\frac{1}{KN_{\rm{r}}}\sum_{\ell=1}^{N_{\rm
r}}\sum_{k=1}^{K}\mathbb{P}\bigg{\\{}\frac{|\hat{\theta}_{\ell,k}-\theta_{\ell,k}|}{\theta_{\ell,k}}>\alpha_{\theta}\bigg{\\}},$

with $\hat{r}_{k}$ as the estimated range of the $k$th drone and
$\hat{\theta}_{\ell,k}$ as the PO of the $k$th drone at the $\ell$th receive
antenna. The EM algorithm was run with $100$ different random initial values
by using the $k$-means$++$ initialization. The number of Monte Carlo runs is
$10^{4}$. We considered the minimization problem in (49) for reordering
estimation.

###  VII-B Simulation Results

Fig. 6 shows range tracking for $K=3$ drones at $\gamma=20$ dB SNR. It is
assumed that the ADS-B packets are periodically broadcasted by the drones
every $T_{\rm p}$ seconds without any coordination between them. The range and
PO variations versus index of the ADS-B packet, $n$, for the three drones are
modeled as $r_{1}(n)=750+250\cos(0.1\pi n)$, $r_{2}(n)=1750+250\cos(0.05\pi
n+2\pi/4)$, $r_{3}(n)=2750+250\cos(0.2\pi n-\pi/3)$,
$\theta_{\ell,k}(n)\in\mathcal{U}_{\rm c}[0,2\pi)$, $k=1,2,3$,
$\ell=1,2,\ldots,5$, where $n=1,2,\ldots,100$,
$\mathbb{E}\\{\theta_{\ell_{1},k_{1}}(n_{1})\theta_{\ell_{2},k_{2}}(n_{2})\\}=(\pi^{2}/3)\delta[\ell_{1}-\ell_{2}]\delta[k_{1}-k_{2}]\delta[n_{1}-n_{2}]$.
As seen, the proposed EM-based estimator can accurately track the range of
$K=3$ drones. While the estimation error for the closest drone to the receiver
is lower, i.e.,
$\mathbb{E}\\{|{r}_{1}-\hat{r}_{1}|\\}^{2}<\mathbb{E}\\{|{r}_{2}-\hat{r}_{2}|\\}^{2}<\mathbb{E}\\{|{r}_{3}-\hat{r}_{3}|\\}^{2}$,
$1-P_{{\rm{out}},r}$ is almost the same for all drones.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 7: The performance of the proposed joint EM-based estimator versus SNR,
$\gamma$, for different number of drones $K$, $\alpha_{\rm r}$, and
$\alpha_{\rm\theta}$, for $N_{\rm r}=5$ receive antennas.

Fig. 7 illustrates $1-P_{{\rm{out}},r}$ and $1-P_{{\rm{out}},\theta}$ of the
proposed EM-based joint estimator versus SNR, $\gamma$, in dB for different
number of drones $K=1,2,3$, $\alpha_{\rm r}$ and $\alpha_{\theta}$. As
expected, by increasing $\gamma$, $1-P_{{\rm{out}},r}$ and
$1-P_{{\rm{out}},\theta}$ increase since the modes of the 2D GM are more
accurately estimated. However, the accuracy of the reordering estimation
through the LS minimization in (49) decreases at higher values of SNRs because
the constraint of the minimization is ${\bf\Phi}\in\mathbb{C}^{K}$, which is
suitable for low and moderate values of SNRs. The reordering estimation via
the minimization in (49) becomes less effective at higher SNRs as can be seen
in Figs. 7(f) and 7(e), where $1-P_{{\rm{out}},r}$ and
$1-P_{{\rm{out}},\theta}$ increase. This problem can be solved by employing
the reordering estimation through the minimization in (39) at the expense of a
higher computational complexity. In this case, $1-P_{{\rm{out}},r}$ and
$1-P_{{\rm{out}},\theta}$ monotonically increase with SNR. Moreover, as seen
in Fig. 7, the larger number of drones, higher estimation error. The reason is
that the number of GM components to be estimated by the EM algorithm
exponentially increases with $K$, i.e., $2^{K}$; however, the number of
observation samples remains fixed.

Figure 8: The effect of the number of receive antennas, $N_{\rm r}$, on
$1-P_{{\rm{out}},r}$ versus $\alpha_{\rm r}$ for $K=2$ drones at $\gamma=20$
dB SNR.

(a)

(b)

Figure 9: The effect of maximum delay, $M$, on the performance of the proposed
EM-based joint estimator for $N_{\rm{r}}=5$ and $K=2$ drones at $\gamma=20$ dB
SNR.

Fig. 8 illustrates the effect of the number of receive antennas on the range
estimation versus $\alpha_{\rm r}$ for $K=2$ drones at $\gamma=20$ dB SNR. As
seen, the larger number of receive antennas, $N_{\rm r}$, higher
$1-P_{{\rm{out}},r}$. Also, the rate of performance improvement decreases as
the number of receive increases. It should be mentioned that increasing number
of receive antennas does not affect $1-P_{{\rm{out}},\theta}$ because spatial
diversity is not employed for the estimation of the $KN_{\rm r}$ independent
POs. However, the range estimation takes the advantage of spatial diversity.
The interesting property of the proposed algorithm is that it can jointly
estimate the range and PO of multiple drones/aircrafts with a single receive
antenna.

(a)

(b)

Figure 10: Performance comparison of the proposed EM-based joint estimator
with the efficient estimator [26] versus $\alpha_{\rm r}$ and
$\alpha_{\rm\theta}$ for $N_{\rm{r}}=1$ and $M=20$ at $\gamma=20$ dB SNR.

Fig. 9 illustrates the effect of the maximum time delay, $M$, on the
performance of the proposed EM-based joint estimator for $K=2$
drones/aircrafts and $N_{\rm r}=5$ receive antennas for different values of
$\alpha_{\rm r}$ and $\alpha_{\rm\theta}$ at $\gamma=20$ dB SNR. As seen, the
proposed EM-based estimator is robust to $M$ for different values of
$\alpha_{\rm r}$ and $\alpha_{\rm\theta}$. Moreover, as seen, the larger $M$,
more accurate ranging and PO estimation because the probability of packet
collision and thus interference decreases.

In Fig. 10, we compare the performance of the proposed EM-based joint
estimator with the efficient estimator [26].555An efficient estimator is an
unbiased estimator that attains the the Cramer-Rao Lower Bound (CRLB) [28].
For the efficient estimator, the transmit symbols and delay of the
drones/aircrafts are assumed to be known a priori at the receiver. As seen,
there is a small gap between the performance of the proposed estimator and the
efficient estimator while the transmit symbols and delay of the
drones/aircraft are unknown in the EM-based estimator. It should be noted that
the efficient estimator represents the upper bound for the mean squared error
(MSE) of the unbised estimators. However, the proposed EM-based estimator is
not an unbiased estimator, and also $1-\alpha_{\rm r}$ and
$1-\alpha_{\rm\theta}$ are the metrics for performance evaluation. Thus, the
EM-based estimator can outperform the efficient estimator in terms of
$1-\alpha_{\rm r}$ and $1-\alpha_{\rm\theta}$.

##  VIII Conclusion

In this paper, we showed that the lost ADS-B packets due to packet collisions
can be employed to jointly estimate range and PO of multiple drones/aircrafts
in the airspace. This enables drones to maintain safe operation distance in
the congested airspace, where ADS-B packet decoding is impossible due to
packet collisions. To achieve this, we derived the maximum likelihood and EM-
based joint range and PO estimators using an approximate PDF obtained by KLD
minimization. The proposed estimators consider uncoordinated and asynchronous
ADS-B packet transmission. For joint range and PO estimation, a priori
knowledge or estimation of the drones’ time delay is not required. Simulation
results showed that the EM-based joint estimator can estimate the range of
multiple drones/aircrafts with a single receive antenna. Performance
improvement by employing multiple receive antennas is obtained.

##  Appendix A Proof of Theorem 1

The vector ${\bf x}_{k}$ given hypothesis $H_{m}^{k}$ in (7) is composed of
$4$ ones and $12$ zeros in ${\bf s}$. The number of zeros and ones in the data
field ${\bf d}_{k}$ is $112$ due to Manchester encoding. Moreover, $M$ zeros
are added irrespective to the hypothesis $H_{m}^{k}$ because of the maximum
integer delay. Hence, the total number of zeros and ones in ${\bf x}_{k}$ are
${M+124}$ and ${116}$, respectively.

Let $\mathcal{X}_{m}^{k}$ denote all possible vector for ${\bf x}_{k}$ given
hypothesis $H_{m}^{k}$. The cardinality of $\mathcal{X}_{m}^{k}$, i.e.,
$|\mathcal{X}_{m}^{k}|=2^{112}$ since the randomness in ${\bf x}_{k}$ because
of the data field ${\bf d}_{k}$ results in $2^{112}$ different Manchester
encoded sequences. Thus; the joint PMF of $x_{k,0},x_{k,1},\cdots,x_{k,N}$ is
given by

$\displaystyle f({\bf x}_{k})=f\big{(}$ $\displaystyle
x_{k,0},x_{k,1},\ldots,x_{k,N}|H_{m}^{k}\big{)}$ (80)
$\displaystyle=\begin{cases}\frac{1}{2^{112}}\,\,\,\,\,\,\,\,\ {\bf
x}_{k}\in\mathcal{X}_{m}^{k},\\\ 0\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\ {\bf
x}_{k}\notin\mathcal{X}_{m}^{k}.\end{cases}$

Let us consider the following approximation for the joint PMF in (80).

$\displaystyle f({\bf x}_{k})=f\big{(}x_{k,0},x_{k,1},$
$\displaystyle\ldots,x_{k,N}|H_{m}^{k}\big{)}$ (81)
$\displaystyle\approx\prod_{n=0}^{N}g(x_{k,n}|H_{m}^{k})$

The KLD for $f$ and $g$ in (81) is given by

$\displaystyle D$
$\displaystyle(f\|g)=\sum_{\mathcal{X}_{m}^{k}}\bigg{[}f\big{(}x_{k,0},x_{k,1},\ldots,x_{k,N}|H_{m}^{k}\big{)}$
$\displaystyle\,\,\,\
\times\ln\frac{f\big{(}x_{k,0},x_{k,1},\cdots,x_{k,N}|H_{m}^{k}\big{)}}{\prod_{n=0}^{N}g(x_{k,n}|H_{m}^{k})}\bigg{]}.$
(82)

Because $x_{k,n}\in\\{0,1\\}$ for $n=0,1,\ldots,N$, the KLD is minimized for
Bernoulli distribution. The PMF of Bernoulli distribution is expressed as

$\displaystyle g(i|H_{m}^{k})=\begin{cases}p&\text{if }i=0,\\\ 1-p&\text{if
}i=1,\end{cases}$ (83)

where $p\in[0,1]$.

By substituting the PMF in (80) and (83) into (82), and by using the fact that
$\sum_{n=0}^{N}x_{k,n}=116$ for ${\bf x}_{k}\in\mathcal{X}_{m}^{k}$, we can
write the KLD in (82) as

$\displaystyle D(f\|g)$
$\displaystyle=\frac{1}{2^{112}}\sum_{\mathcal{X}_{m}^{k}}\ln\bigg{[}\frac{1}{2^{212}p^{M+124}(1-p)^{116}}\bigg{]}$
(84) $\displaystyle=-\ln\Big{[}{2^{212}p^{M+124}(1-p)^{116}}\Big{]}.$

To obtain $p$, we need to minimize $D(f\|g)$. Because $\ln(\cdot)$ is a
monotonically increasing function, $D(f\|g)$ is minimized when
$p^{M+124}(1-p)^{116}$, $p\in[0,1]$, is maximized.

By taking the derivative with respect to $p$ and setting it to zero, we obtain

$\displaystyle\diff{\big{(}p^{M+124}(1-p)^{116}\big{)}}$
$\displaystyle={(M+124)}p^{M+123}(1-p)^{116}$ (85)
$\displaystyle-116(1-p)^{115}p^{M+124}=0.$

By solving (85), we obtain (9).

##  Appendix B Proof of the PDF of $G=\sum_{k=1}^{K}Z_{k}$

Let us consider $G=\sum_{k=1}^{K}Z_{k}$. Since $Z_{k}$, $k=1,2,\cdots,K$, are
independent random variables, the PDF of $G$ is given by

$\displaystyle f_{G}(g;p,{\bf{h}})$
$\displaystyle=f_{{Z}_{1}}\big{(}z;p,h_{1}\big{)}*f_{{Z}_{2}}\big{(}z;p,h_{2}\big{)}$
(86) $\displaystyle\,\,\,\,\,\ *\ldots*f_{{Z}_{N}}\big{(}z;p,h_{K}\big{)},$

where $z\triangleq z_{\rm r}+iz_{\rm I}\in\mathbb{C}$,
${\bf{h}}\triangleq[h_{1},h_{2},\cdots,h_{K}]^{T}$,
$h_{k}\triangleq{(h_{k})}_{\rm r}+i{(h_{k})}_{\rm I}\in\mathbb{C}$, and
$f_{{Z}_{k}}\big{(}z;p,h_{k}\big{)}$ is given in (10). The two-dimensional
Laplace transform of $f_{{Z}_{k}}\big{(}z;p,h_{k}\big{)}$ is expressed as

$\displaystyle
F_{Z_{k}}\big{(}s_{1},s_{2};p,h_{k}\big{)}=\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}f_{{Z}_{k}}\big{(}z;p,h_{k}\big{)}e^{-s_{1}{z_{\rm
r}}-s_{2}{z}_{\rm I}}{\rm d}{z}_{\rm r}{\rm d}z_{\rm I}$
$\displaystyle=\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}p\delta(z_{\rm
r})\delta({z}_{\rm I})e^{-s_{1}{z_{\rm r}}-s_{2}{z_{\rm I}}}{\rm d}{z}_{\rm
r}{\rm d}z_{\rm I}$
$\displaystyle+\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}(1-p)\delta({z_{\rm
r}}-{(h_{k})}_{\rm r})\delta(z_{\rm I}-{(h_{k})}_{\rm I})e^{-s_{1}{z_{\rm
r}}-s_{2}{z_{\rm I}}}{\rm d}{z}_{\rm r}{\rm d}z_{\rm I}$
$\displaystyle=p+(1-p)e^{-({(h_{k)}}_{\rm r}s_{1}+{(h_{k})}_{\rm I}s_{2})}.$

Using the fact that the linear convolution is equivalent to multiplication in
the Laplace domain, we can write

$\displaystyle F_{G}\big{(}s_{1},s_{2};p,{\bf h}\big{)}$
$\displaystyle=\prod_{k=1}^{K}F_{Z_{k}}\big{(}s_{1},s_{2};p,h_{k}\big{)}$ (87)
$\displaystyle=\prod_{k=1}^{K}\Big{[}p+(1-p)e^{-({(h_{k})}_{\rm
r}s_{1}+{(h_{k})}_{\rm I}s_{2})}\Big{]},$

where $F_{G}\big{(}s_{1},s_{2};p,{\bf h}\big{)}$ denotes the two-dimensional
Laplace transform of $f_{G}(g;p,{\bf{h}})$. Let us consider the multi-binomial
theorem as follows [22]

$\displaystyle\prod_{i=1}^{d}(a_{i}+b_{i})^{n_{i}}$ (88)
$\displaystyle=\sum_{v_{1}=0}^{n_{1}}\cdots\sum_{v_{d}=0}^{n_{d}}\binom{n_{1}}{v_{1}}a_{1}^{v_{1}}b_{1}^{n_{1}-v_{1}}\cdots\binom{n_{d}}{v_{d}}a_{d}^{v_{d}}b_{d}^{n_{d}-v_{d}}.$

By substituting $a_{i}\triangleq p$, $b_{i}\triangleq(1-p)e^{-({(h_{i})}_{\rm
r}s_{1}+{(h_{i})}_{\rm I}s_{2})}$, $d=K$, and $n_{i}=1$, $i=1,2,\ldots,d$,
into (88), we obtain

$\displaystyle F_{G}$ $\displaystyle\big{(}s_{1},s_{2};p,{\bf
h}\big{)}=\sum_{v_{1}=0}^{1}\cdots\sum_{v_{K}=0}^{1}\bigg{[}p^{\sum_{k=1}^{K}v_{k}}(1-p)^{K-{\sum_{k=1}^{K}v_{k}}}$
$\displaystyle\times\exp{\Big{[}-\sum_{k=1}^{K}(1-v_{k})({({{h}}_{k})}_{\rm
r}s_{1}+{(h_{k})}_{\rm I}s_{2})\Big{]}}\bigg{]}.$ (89)

The two-dimensional inverse Laplace transform of the exponential term in (89)
is given by

$\displaystyle\mathcal{L}^{-1}\bigg{\\{}\exp{\Big{[}-\sum_{k=1}^{K}(1-v_{k})({(h_{k})}_{\rm
r}s_{1}+{(h_{k})}_{\rm I}s_{2})\Big{]}}\bigg{\\}}$ (90)
$\displaystyle=\delta\Big{(}s_{\rm r}-\sum_{k=1}^{K}(1-v_{k}){(h_{k})}_{\rm
r}\Big{)}\delta\Big{(}s_{\rm I}-\sum_{k=1}^{K}(1-v_{k}){(h_{k})}_{\rm
I}\Big{)}$ $\displaystyle=\delta_{\rm
c}\Big{(}s-\sum_{k=1}^{K}(1-v_{k})h_{k}\Big{)},$

where $s=s_{\rm r}+is_{\rm I}\in{\mathbb{C}}$. By taking two-dimensional
inverse Laplace transform from both sides of (89) and then employing (90), we
obtain $f_{G}(g;p,{\bf{h}})$ as in (12).

## References

  * [1] M. Strohmeier, M. Schäfer, V. Lenders, and I. Martinovic, “Realities and challenges of nextgen air traffic management: the case of ADS-B,” _IEEE Commun. Mag._ , vol. 52, no. 5, pp. 111–118, May 2014.
  * [2] E. Tohidi, A. Hariri, H. Behroozi, M. M. Nayebi, G. Leus, and A. P. Petropulu, “Compressed-domain detection and estimation for colocated MIMO radar,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 56, no. 6, pp. 4504–4518, May 2020.
  * [3] Y. Kim, J.-Y. Jo, and S. Lee, “ADS-B vulnerabilities and a security solution with a timestamp,” _IEEE Aerosp. Electron. Syst. Mag._ , vol. 32, no. 11, pp. 52–61, Nov. 2017.
  * [4] F. Minucci, E. Vinogradov, and S. Pollin, “Avoiding collisions at any (low) cost: ADS-B like position broadcast for UAVs,” _IEEE Access_ , vol. 8, pp. 121 843–121 857, July 2020.
  * [5] E. Vinogradov, F. Minucci, and S. Pollin, “Wireless communication for safe UAVs: From long-range deconfliction to short-range collision avoidance,” _IEEE Veh. Technol. Mag._ , vol. 15, no. 2, pp. 88–95, Apr. 2020.
  * [6] A. Abdulaziz, A. S. Yaro, A. A. Adam, M. T. Kabir, and H. B. Salau, “Optimum receiver for decoding automatic dependent surveillance broadcast (ADS-B) signals,” _American Journal of Signal Processing_ , vol. 5, no. 2, pp. 23–31, 2015.
  * [7] P. Ren, J. Wang, H. Song, P. Zhang, and T. Ban, “A novel multi-criteria preamble detection algorithm for ADS-B signals,” _IEEE Access_ , vol. 7, pp. 97 319–97 332, July 2019.
  * [8] K. Li, J. Kang, H. Ren, and Q. Wu, “A reliable separation algorithm of ADS-B signal based on time domain,” _IEEE Access_ , vol. 9, pp. 88 019–88 026, 2021.
  * [9] W. Wang, R. Wu, and J. Liang, “ADS-B signal separation based on blind adaptive beamforming,” _IEEE Trans. Veh. Technol._ , vol. 68, no. 7, pp. 6547–6556, 2019.
  * [10] S. Sciancalepore and R. Di Pietro, “SOS: Standard-compliant and packet loss tolerant security framework for ADS-B communications,” _IEEE Trans. Depend. Sec. Comput._ , vol. 18, no. 4, pp. 1681–1698, 2019.
  * [11] H. Yang, Q. Zhou, M. Yao, R. Lu, H. Li, and X. Zhang, “A practical and compatible cryptographic solution to ADS-B security,” _IEEE Internet Things J._ , vol. 6, no. 2, pp. 3322–3334, 2018.
  * [12] Z. Wu, T. Shang, and A. Guo, “Security issues in automatic dependent surveillance-broadcast (ADS-B): A survey,” _IEEE Access_ , vol. 8, pp. 122 147–122 167, July 2020.
  * [13] S. Zhu, X. Zheng, L. Liu, and H. Ma, “AirSync: Time synchronization for large-scale iot networks using aircraft signals,” _IEEE Trans. Mobile Comput._ , pp. 1681–1698, Apr. 2021.
  * [14] S.-L. Jheng, S.-S. Jan, Y.-H. Chen, and S. Lo, “1090 MHz ADS-B-based wide area multilateration system for alternative positioning navigation and timing,” _IEEE Sensors J._ , vol. 20, no. 16, pp. 9490–9501, Apr. 2020\.
  * [15] M. Leonardi, E. Piracci, and G. Galati, “ADS-B jamming mitigation: A solution based on a multichannel receiver,” _IEEE Aerosp. Electron. Syst. Mag._ , vol. 32, no. 11, pp. 44–51, Nov. 2017.
  * [16] S. Khan, J. Thorn, A. Wahlgren, and A. Gurtov, “Intrusion detection in automatic dependent surveillance-broadcast (ADS-B) with machine learning,” in _Proc. IEEE DASC_ , Oct. 2021, pp. 1–10.
  * [17] Z. Luo, C. Li, and L. Zhu, “A comprehensive survey on blind source separation for wireless adaptive processing: Principles, perspectives, challenges and new research directions,” _IEEE Access_ , vol. 6, pp. 66 685–66 708, 2018\.
  * [18] G. Chabriel, M. Kleinsteuber, E. Moreau, H. Shen, P. Tichavsky, and A. Yeredor, “Joint matrices decompositions and blind source separation: A survey of methods, identification, and applications,” _IEEE Signal Processing Magazine_ , vol. 31, no. 3, pp. 34–43, 2014.
  * [19] G. R. Naik and D. K. Kumar, “An overview of independent component analysis and its applications,” _Informatica_ , vol. 35, no. 1, 2011.
  * [20] M. Shi, Y. Bar-Ness, and W. Su, “Adaptive estimation of the number of transmit antennas,” in _Proc. IEEE MILCOM_ , Feb. 2007, pp. 1–5.
  * [21] M. Mohammadkarimi, E. Karami, O. A. Dobre, and M. Z. Win, “Number of transmit antennas detection using time-diversity of the fading channel,” _IEEE Trans. Signal Process._ , vol. 65, no. 15, pp. 4031–4046, Mar. 2017.
  * [22] C. Morris, “Central limit theorems for multinomial sums,” _The Annals of Statistics_ , pp. 165–188, Jan. 1975.
  * [23] G. J. McLachlan, S. X. Lee, and S. I. Rathnayake, “Finite mixture models,” _Annual review of statistics and its application_ , vol. 6, pp. 355–378, 2019\.
  * [24] S. Theodoridis, _Machine Learning: a Bayesian and Optimization Perspective_. Academic press, 2015.
  * [25] T. K. Moon, “The expectation-maximization algorithm,” _IEEE Signal Process. Mag._ , vol. 13, no. 6, pp. 47–60, Nov. 1996.
  * [26] S. M. Kay, _Fundamentals of Statistical Signal Processing: Estimation Theory_. Prentice-Hall, Inc., 1993.
  * [27] V. Chandola, A. Banerjee, and V. Kumar, “Anomaly detection: A survey,” _ACM Computing Surveys_ , vol. 41, p. 15, 2009.
  * [28] Y. Wang, G. Leus, and A.-J. van der Veen, “Cramer-Rao bound for range estimation,” in _Proc. IEEE ICASP_ , 2009, pp. 3301–3304.

  *[PO]: Phase Offset
  *[ADS-B]: Automatic Dependent Surveillance–Broadcast
  *[KLD]: Kullback–Leibler Divergence
  *[AWGN]: Additive White Gaussian Noise
  *[i.i.d.]: independent and identically distributed
  *[GM]: Gaussian Mixture
  *[MLE]: Maximum Likelihood Estimation
  *[PDF]: Probability Density Function
  *[EM]: Expectation–Maximization
  *[GPS]: Global Positioning System
  *[IMU]: Inertial Measurement Unit
  *[SNRs]: Signal-to-Noise Ratios
  *[SoOp]: Signal-of-Opportunity
  *[ML]: Maximum Likelihood
  *[GC]: Ground Controller
  *[PPM]: Pulse Position Modulation
  *[SDR]: Software-defined Radio
  *[PSD]: Power Spectral Density
  *[CFO]: Carrier Frequency Offset
  *[PMF]: Probability Mass Function
  *[SNR]: Signal-to-noise Ratio
  *[LS]: Least Squares


# Forged Image Detection using SOTA Image Classification Deep Learning Methods
for Image Forensics with Error Level Analysis

Raunak Joshi
Department of IT
University of Mumbai
Mumbai, 400032, India
<EMAIL_ADDRESS>
&Abhishek Gupta
Department of EXTC
University of Mumbai
Mumbai, 400032, India
<EMAIL_ADDRESS>
&Nandan Kanvinde
Department of MCA
TIMSCDR
Mumbai, 400101, India
<EMAIL_ADDRESS>
&Pandharinath Ghonge
Department of EXTC
St. John College of Engineering and Management
Palghar - 401404, India
<EMAIL_ADDRESS>

###### Abstract

The advancement in the area of computer vision has been brought using deep
learning mechanisms. Image Forensics is one of the major areas of computer
vision application. Forgery of images is sub-category of image forensics and
can be detected using Error Level Analysis. Using such images as an input,
this can turn out to be a binary classification problem which can be leveraged
using variations of convolutional neural networks. In this paper we perform
transfer learning with state-of-the-art image classification models over error
level analysis induced CASIA ITDE v.2 dataset. The algorithms used are VGG-19,
Inception-V3, ResNet-152-V2, XceptionNet and EfficientNet-V2L with their
respective methodologies and results.

_K_ eywords Image Classification $\cdot$ Image Forensics $\cdot$ Error Level
Analysis $\cdot$ Forged Image Detection $\cdot$ Transfer Learning

##  1 Introduction

The area of computer vision [1] has excelled in terms of innovation and
performance delivered by leveraging Deep Learning [2]. The various tasks of
computer vision are classification, object detection [3], object counting [4],
image segmentation [5, 6] and many more. Classification [7] can be termed as
one of the most primordial tasks in computer vision. The task of
classification is identification of an entity or object by prognosticating its
appropriate label is done effectively using Convolutional Neural Networks [8]
abbreviated as CNN. The standard CNN gone under massive improvements in the
latter period when ImageNet [9] Large Scale Visual Recognition Challenge
(ILSVRC) [10] came into existence yielding many models that are currently
state-of-the-art deep learning models for image classification. These models
were later refined to a greater extent and turned out to be a base for many
new models effectively. One sub-category in the area of classification is
image forgery [11], where the deep learning can be used with binary
classification [12, 13, 14] problem. The task will be simple as classifying
whether the image is original or fake, hence leading to binary classification
problem.

The image forgery can be termed as a sub-category of image forensics [15]
where there are many methods where Error Level Analysis [16] is one of the
commonly used techniques. The image when tampered performs compression which
is used as the main point for Error Level Analysis abbreviated as ELA. The
compression is obtained at 95% of the input image and saved, which is later
evaluated with the original image for differences. Images have grid squares
which are thoroughly evaluated by ELA which yield errors due to successive
saving of the image after editing. The forged section of the image is bright
after ELA operation whereas non-tampered area is darker than all the parts of
the image. The section highlights with lighter shades where tampering is done
specifically. This data can be labelled and given to CNN for classification
where the ELA filters over an image are input. Standard CNN application [17]
is very obvious and lack the necessary performance where the transfer learning
[18] can be performed with state-of-the-art models. The varied popular models
can be used, viz. VGG [19], ResNet [20, 21], InceptionNet [22], XceptionNet
[23] and EfficientNet [24]. The appropriate data for such an application is
CASIA dataset, where the CASIA ground-truth [25] is the first dataset that has
distinct 8 classes of images. The tampering over this dataset was performed
which yielded CASIA ITDE v.2 [26] dataset having 2 distinct classes, authentic
images and forged images. This paper is compendious comparison of state-of-
the-art models using transfer learning over the CASIA ITDE v.2, which will be
intricately explained in further sections of this paper.

##  2 Methodology

###  2.1 VGG-19

The VGG [19] is one of the primarily used deep learning computer vision
classification based algorithm which has 2 variants, 16 layered known as
VGG-16 and 19 layered known as VGG-19. The network consists of blocks, where
each block has 2-dimensional convolutional layers and 2-dimensional max-
pooling [27] layers. The VGG-19 has total 5 blocks, where first 2 blocks have
2 convolutional layers followed by one max-pooling layer. The remaining 3
blocks have 4 convolutional layers followed by max-pooling layer. Since the
model was originally trained on imagenet [9] data which has 1000 classes, we
are using the same trained weights of the imagenet for transfer learning. We
ignore the output layer of model and append 2-dimensional global average
pooling [28] layer followed by 2 layers. The first layer has 1024 hidden
neurons with ReLU [29] activation function, followed by the output layer that
has 2 hidden neurons with softmax [30] activation function. Since we are
dealing with binary classification problem, we could use 1 hidden neuron with
sigmoid [31] activation function, but we wanted the output probabilities for
inference examination in a through manner. The network for backward
propagation [32] uses binary cross-entropy loss function [33] and Adam [34]
optimizer for loss optimization. The network has learned 20,551,746 total
parameters which are approximately, 20.55 million parameters. If we remove the
personally initialized layers the VGG-19 itself consists of roughly 20.02
million parameters.

###  2.2 Inception-V3

The computationally efficient network was definitely the InceptionNet [22]
consisting of 22 layers. The inception module was the main point of focus for
the network that later made it a state-of-the-art model. The inception modules
force the network to train many different filters in parallel fashion with
concatenation for depth-wise aspect. This network over the period of time had
many improvements, where batch normalization layers were used which is known
as Inception-V2 [35] network. Later the network was made deeper with 42 layers
which still turned out to be computationally efficient as compared to VGG.
This modified network was Inception-V3 [36] that we have used, where we
stripped of the output layers and applied 2-dimensional global average pooling
[28] followed by one 1024 hidden neuron ReLU [29] activation layer and one 2
hidden neuron softmax [30] activation output layer. The loss function is
binary cross-entropy and loss optimizer is Adam. The network has total
23,903,010 parameters, which approximately are around 23.90 million. The
trainable parameters are 23,868,578 and non-trainable parameters are 34,432.
21,802,784 are the sole parameters of the Inception-V3 which are roughly 21.80
million respectively.

###  2.3 ResNet

The ResNet [20] is termed as one of the most prominently used network for
image classification tasks. The main problem that ResNet emphasized on was the
degradation problem with convergence of the network. This is basically
saturation of the accuracy with increasing depth of the network. For this the
residual block[21] was developed which uses a concept known as skip connection
which abruptly skips the layers in between while back-propagation [32] to
avoid the saturation in loss optimization. The network we have used is ResNet
with 152 layers. Now it has 2 versions, where version 1 uses non-linearity in
the last layer and version 2 has removed it. We have used the ResNet 152 v2
[37], where we have removed the output layer and added 2-dimensional global
average pooling [28] with 1024 hidden neuron ReLU [29] activation layer
followed by 2 hidden neuron softmax [30] activation output layer. The loss
function and optimizer is same as the earlier networks and it learns
60,431,874 parameters which are approximately 60.4 million parameters. After
removing the custom added layers, the ResNet alone has 58,331,648 parameters.

###  2.4 XceptionNet

The XceptionNet [23] is 71 layers deep extension network of Inception-V3 [36]
that replaces the inception modules with depth-wise separable convolutions.
The depth of the input is not only covered in this network but also spatial
dimensions which makes it special which also covers the kernels that cannot be
separated into further sections. We remove the output layer and add custom
layers where first 2-dimensional global average pooling is applied. Followed
by it is one layer with 1024 hidden neurons and ReLU activation function.
Finally one output layer with 2 hidden neurons and softmax activation function
is used. The loss function used is binary cross-entropy and loss optimizer
used is Adam. The network learns total of 22,961,706 parameters which are
approximately 22.96 million. The non-trainable parameters are 54,528 only and
22,907,178 are trainable. The parameters excluding the custom defined layers
are 20,861,480, which roughly are around 20.86 million respectively.

###  2.5 EfficientNet-V2L

The need of efficiency in training and computation of the networks has risen
in a broader perspective for which EfficientNet [24] drives as a good example.
The scaling of the base features of the input like width, height and depth is
done arbitrarily which has been made uniform with EfficientNet. The compound
scaling is the main feature and many variation of it were made where
advancement also reached the version 2. We have used EfficientNet-V2L [38]
which is large network with over 110 million parameters. The network we have
trained has 119,060,642 parameters which are around 119 million roughly where
the non-trainable parameters are around 500k and rest are trainable. The
network without our layers is 117,746,848 parameters which is approximately
117.74 million parameters. The last layer of the EfficientNet-V2L is replaced
and 2-dimensional global average pooling is applied. Followed by it is one
dense layer with 1024 hidden neurons and ReLU activation function. The final
output layer is applied with 2 hidden neurons and softmax activation function.
The loss function used is binary cross-entropy and adam loss optimizer for
loss optimization.

##  3 Results

###  3.1 Training and Validation Accuracy

The below figure 1 gives a general overview of the training and validation
accuracy for all the models used in the implementation. The highest accuracy
in the last epoch of training accuracy for VGG-19 is 95.46% and validation
accuracy is 93.51% which is very effectively generalized. Similarly the
training accuracy for Inception-V3 is 97.35% whereas the validation accuracy
is 89.67% which makes model less generalized but a better training accuracy
than VGG-19. The training accuracy for ResNet152-V2 is 99.6% and validation
accuracy is 90.75%. The training accuracy for XceptionNet is 97.71% and
validation accuracy is 90.87%. Finally the training accuracy for
EfficientNet-V2L is 96.36% and validation accuracy is 93.15% making it
generalized and effective as compared to other networks yet this is not
sufficient metrics for stating an inference.

$0$$2$$4$$6$$8$$0.7$$0.8$$0.9$EpochsTraining AccuracyVGG-19Inception-V3ResNet-152-V2XceptionNetEfficientNet-V2L | $0$$2$$4$$6$$8$$0.6$$0.7$$0.8$$0.9$EpochsValidation AccuracyVGG-19Inception-V3ResNet-152-V2XceptionNetEfficientNet-V2L
---|---
Figure 1: Training Accuracy and Validation Accuracy

###  3.2 Training and Validation Loss

The figure 2 gives a graphical depiction of the models converging throughout
the epochs. The training loss for VGG-19 is 12.25% in the last epochs whereas
the validation loss is 18.72%. The training loss for Inception-V3 is 8.8% and
validation loss is 28.91%. The training loss for ResNet-152V2 is 2.39% while
the validation loss is 30.03% effectively. The training loss for XceptionNet
is 5.98% whereas the validation loss is 28.48%. Finally the training loss for
EfficientNet-V2L is 9.67% and validation loss is 21.05% respectively. The
generalization can be seen better with VGG-19 but highest training loss
convergence can be seen with ResNet-152V2. EfficientNet-V2L attains a middle
ground effectively.

$0$$2$$4$$6$$8$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$EpochsTraining LossVGG-19Inception-V3ResNet-152-V2XceptionNetEfficientNet-V2L | $0$$2$$4$$6$$8$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$EpochsValidation LossVGG-19Inception-V3ResNet-152-V2XceptionNetEfficientNet-V2L
---|---
Figure 2: Training Loss and Validation Loss

###  3.3 Precision

The precision [39] in simple terms is positive sample for classification in
model accuracy. The better explanation can be given as number of positive
samples classified with respect to total number of samples. The table 1 below
gives an explanation of precision scores for all the respective models.

Table 1: Precision Score Model | Precision Score
---|---
VGG-19 | 98.25%
Inception-V3 | 94.11%
ResNet-152-V2 | 90.38%
XceptionNet | 94.61%
EfficientNet-V2L | 85.20%

The precision obtained for VGG-19 is highest which is an indication of higher
number of positive samples classified as per the total samples effectively.
The precision alone though is not the perfect metric for deriving the final
inference.

###  3.4 Recall

The recall [39] is number of positive samples classified appropriately with
respect to total number of positive samples. The table 2 gives the recall
scores for all the algorithms used.

Table 2: Recall Score Model | Recall Score
---|---
VGG-19 | 79.34%
Inception-V3 | 82.62%
ResNet-152-V2 | 88.26%
XceptionNet | 86.61%
EfficientNet-V2L | 94.60%

The Recall score is highest for EfficientNet-V2L which is was directly lowest
in terms of precision. This is the reason recall also cannot be used for final
inference judgement.

###  3.5 RoC and AUC

The Receiver operating characteristic [40] curve is specifically for binary
classification that gives the performance of classification at all the
thresholds ranging from 0 to 1. It uses false and true positive rates for
evaluation of the RoC curve which the abbreviation basically. The area under
curve is also its crucial aspect where the entire area covering the RoC and
signifies the separability of classes. It is difficult to make a judgement
from figure 3 since the predictions are merging, for which zooming into the
depth can help derive some intuition.

$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$0.2$$0.4$$0.6$$0.8$$1$False Positive RateTrue
Positive RateVGG-19Inception-V3ResNet-152-V2XceptionNetEfficientNet-V2L Figure
3: RoC and AUC $0.2$$0.4$$0.6$$0.8$$1$$0.9$$0.92$$0.94$$0.96$$0.98$$1$False
Positive RateTrue Positive
RateVGG-19Inception-V3ResNet-152-V2XceptionNetEfficientNet-V2L Figure 4:
Zoomed RoC with AUC
$0.2$$0.25$$0.3$$0.35$$0.4$$0.45$$0.5$$0.94$$0.95$$0.95$$0.96$$0.96$$0.97$$0.97$False
Positive RateTrue Positive
RateVGG-19Inception-V3ResNet-152-V2XceptionNetEfficientNet-V2L Figure 5: Super
Zoomed RoC with AUC in range of 94% to 97% Y-axis

The figure 4 shows the pivots every algorithm made after a stagnant set of
predictions. All the algorithms eventually reach the peak that is 1 value. The
EfficientNet-V2L has it around 96% as seen in figure 5 clearly. Big leap for
Inception-V3 is at 95% approximately. For ResNet-152-V2 it is somewhere above
95% approximately. The VGG-19 turns just before 96% and XceptionNet clearly at
96%. These values are an indication of the learning procedures adopted for
every classification threshold and area under their curve.

###  3.6 F Measure

All the metrics are very contradictory for making a judgement for the best
algorithm. The precision indicates that VGG-19 is better whereas the
EfficientNet-V2L is highest for recall. These both inferences are very far
away from each other and proper metric is required. The averaging system for
this is required where the F Measure also known as F-Score [39] is used.

Table 3: F-Measure for all the used Models Model | F-Measure
---|---
VGG-19 | 87.79%
Inception-V3 | 88%
ResNet-152-V2 | 89.31%
XceptionNet | 90.44%
EfficientNet-V2L | 89.65%

The F measure has 1 as the $\beta$ value where the precision and recall are in
the balanced state for metric measuring. If the $\beta$ value is 0.5 or 2, the
metric inclines towards precision and recall respectively. These values can be
represented with a table for comparison. The table 3 is the depiction of F
measure for all the models used where XceptionNet outperforms all the models.
Now this is most different observation and least expected. Since one relies on
half measures, the results can turn out to be very varied but the least
expected always turns out to be correct.

##  4 Conclusion

The area of image forgery is a potential application which in this paper has
been leveraged using state-of-the-art deep learning image classification
models with transfer learning. The dataset we used was CASIA ITDE v.2 which is
a binary classification dataset that differentiates between authentic and
tampered images. After performing the error level analysis from image
forensics, the variety of famous models were used. Making an inference for the
best model was a difficult task indeed due to results section variability yet
the right judgement was attained. This paper does open potential thought for
transfer learning models and wide variety of applications in the area of deep
learning for computer vision. The limitations of the paper are only bounded
and entitled to subtle differences between various state-of-the-art models
that are very well explained throughout the result section of the paper. The
distinction and clarity of reaching an inference in the paper is bounded to
the methods of proving a point and we hope with our best belief that some more
extensions with more better results can be obtained in future work.

## References

  * [1] Victor Wiley and Thomas Lucas.  Computer vision and image processing: A paper review.  International Journal of Artificial Intelligence Research, 2018\.
  * [2] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  nature, 521(7553):436–444, 2015.
  * [3] Zhong-Qiu Zhao, Peng Zheng, Shou-tao Xu, and Xindong Wu.  Object detection with deep learning: A review.  IEEE transactions on neural networks and learning systems, 30(11):3212–3232, 2019.
  * [4] Shuyue Guan and Murray Loew.  Understanding the ability of deep neural networks to count connected components in images.  In 2020 IEEE Applied Imagery Pattern Recognition Workshop (AIPR), pages 1–7. IEEE, 2020.
  * [5] Shervin Minaee, Yuri Y Boykov, Fatih Porikli, Antonio J Plaza, Nasser Kehtarnavaz, and Demetri Terzopoulos.  Image segmentation using deep learning: A survey.  IEEE transactions on pattern analysis and machine intelligence, 2021\.
  * [6] Raunak M Joshi and Deven Shah.  Refactoring faces under bounding box using instance segmentation algorithms in deep learning for replacement of editing tools.  In Intelligent Computing and Networking, pages 236–247. Springer, 2022.
  * [7] Robin Cormack.  A review of classification.  1971\.
  * [8] Yann LeCun, Yoshua Bengio, et al.  Convolutional networks for images, speech, and time series.
  * [9] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei.  Imagenet: A large-scale hierarchical image database.  In 2009 IEEE Conference on Computer Vision and Pattern Recognition, pages 248–255, 2009.
  * [10] Olga Russakovsky, Jia Deng, Hao Su, Jonathan Krause, Sanjeev Satheesh, Sean Ma, Zhiheng Huang, Andrej Karpathy, Aditya Khosla, Michael Bernstein, et al.  Imagenet large scale visual recognition challenge.  International journal of computer vision, 115(3):211–252, 2015\.
  * [11] Navpreet Kaur Gill, Ruhi Garg, and Er. Amit Doegar.  A review paper on digital image forgery detection techniques.  In 2017 8th International Conference on Computing, Communication and Networking Technologies (ICCCNT), pages 1–7, 2017.
  * [12] Roshan Kumari and Saurabh Kr Srivastava.  Machine learning: A review on binary classification.  International Journal of Computer Applications, 160(7), 2017.
  * [13] Nandan Kanvinde, Abhishek Gupta, and Raunak Joshi.  Binary classification for high dimensional data using supervised non-parametric ensemble method.  arXiv preprint arXiv:2202.07779, 2022.
  * [14] Abhishek Gupta, Sarvesh Thustu, Riti Thakor, Saniya Patil, Raunak Joshi, and Ronald Melvin Laban.  Prediction of maneuvering status for aerial vehicles using supervised learning methods.  arXiv preprint arXiv:2206.10303, 2022.
  * [15] Alessandro Piva.  An overview on image forensics.  International Scholarly Research Notices, 2013:1–22, 2013.
  * [16] Nor Bakiah Abd Warif, Mohd. Yamani Idna Idris, Ainuddin Wahid Abdul Wahab, and Rosli Salleh.  An evaluation of error level analysis in image forensics.  In 2015 5th IEEE International Conference on System Engineering and Technology (ICSET), pages 23–28, 2015.
  * [17] Abhishek Gupta, Raunak Joshi, and Ronald Laban.  Detection of tool based edited images from error level analysis and convolutional neural network.  arXiv preprint arXiv:2204.09075, 2022.
  * [18] Fuzhen Zhuang, Zhiyuan Qi, Keyu Duan, Dongbo Xi, Yongchun Zhu, Hengshu Zhu, Hui Xiong, and Qing He.  A comprehensive survey on transfer learning.  Proceedings of the IEEE, 109(1):43–76, 2020.
  * [19] Karen Simonyan and Andrew Zisserman.  Very deep convolutional networks for large-scale image recognition.  arXiv preprint arXiv:1409.1556, 2014.
  * [20] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [21] Abhishek Gupta, Sruthi Nair, Raunak Joshi, and Vidya Chitre.  Residual-concatenate neural network with deep regularization layers for binary classification.  In 2022 6th International Conference on Intelligent Computing and Control Systems (ICICCS), pages 1018–1022. IEEE, 2022.
  * [22] Christian Szegedy, Wei Liu, Yangqing Jia, Pierre Sermanet, Scott Reed, Dragomir Anguelov, Dumitru Erhan, Vincent Vanhoucke, and Andrew Rabinovich.  Going deeper with convolutions.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 1–9, 2015.
  * [23] François Chollet.  Xception: Deep learning with depthwise separable convolutions.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 1251–1258, 2017.
  * [24] Mingxing Tan and Quoc Le.  Efficientnet: Rethinking model scaling for convolutional neural networks.  In International conference on machine learning, pages 6105–6114. PMLR, 2019.
  * [25] Nam Thanh Pham, Jong-Weon Lee, Goo-Rak Kwon, and Chun-Su Park.  Hybrid image-retrieval method for image-splicing validation.  Symmetry, 11(1):83, 2019.
  * [26] Jing Dong, Wei Wang, and Tieniu Tan.  Casia image tampering detection evaluation database.  In 2013 IEEE China Summit and International Conference on Signal and Information Processing, pages 422–426, 2013.
  * [27] Hossein Gholamalinezhad and Hossein Khosravi.  Pooling methods in deep neural networks, a review.  arXiv preprint arXiv:2009.07485, 2020.
  * [28] Min Lin, Qiang Chen, and Shuicheng Yan.  Network in network.  arXiv preprint arXiv:1312.4400, 2013.
  * [29] Abien Fred Agarap.  Deep learning using rectified linear units (relu).  arXiv preprint arXiv:1803.08375, 2018.
  * [30] John S. Bridle.  Training stochastic model recognition algorithms as networks can lead to maximum mutual information estimation of parameters.  In Proceedings of the 2nd International Conference on Neural Information Processing Systems, NIPS’89, page 211–217, Cambridge, MA, USA, 1989\. MIT Press.
  * [31] Ali A. Minai and Ronald D. Williams.  Original contribution: On the derivatives of the sigmoid.  Neural Netw., 6(6):845–853, jun 1993.
  * [32] David E Rumelhart, Geoffrey E Hinton, and Ronald J Williams.  Learning internal representations by error propagation.  Technical report, California Univ San Diego La Jolla Inst for Cognitive Science, 1985.
  * [33] Usha Ruby.  Binary cross entropy with deep learning technique for image classification.  2020\.
  * [34] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [35] Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In International conference on machine learning, pages 448–456. PMLR, 2015.
  * [36] Christian Szegedy, Vincent Vanhoucke, Sergey Ioffe, Jon Shlens, and Zbigniew Wojna.  Rethinking the inception architecture for computer vision.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 2818–2826, 2016.
  * [37] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Identity mappings in deep residual networks.  In European conference on computer vision, pages 630–645. Springer, 2016.
  * [38] Mingxing Tan and Quoc Le.  Efficientnetv2: Smaller models and faster training.  In International Conference on Machine Learning, pages 10096–10106. PMLR, 2021.
  * [39] David MW Powers.  Evaluation: from precision, recall and f-measure to roc, informedness, markedness and correlation.  arXiv preprint arXiv:2010.16061, 2020.
  * [40] Tom Fawcett.  An introduction to roc analysis.  Pattern Recognit. Lett., 27:861–874, 2006.