figure table

# Super-resolving multiple scatterers detection in SAR Tomography assisted by
correlation information

Ahmad Naghavi  Isfahan University of Technology, Department of Electrical and
Computer Engineering, Isfahan, Iran  Mohammad Sadegh Fazel  Isfahan University
of Technology, Department of Electrical and Computer Engineering, Isfahan,
Iran  Mojtaba Beheshti  Research Institute for Information and Communication
Technology, Isfahan University of Technology, Isfahan, Iran  Ehsan Yazdian
Isfahan University of Technology, Department of Electrical and Computer
Engineering, Isfahan, Iran

###### Abstract

This paper proposes a method for detecting multiple scatterers (targets) in
the elevation direction for synthetic aperture radar (SAR) tomography. The
proposed method can resolve closely spaced targets through a two-step
procedure. In the first step, coarse detection is performed with a successive
cancellation scheme in which possible locations of targets are marked. Then,
in the second step, by searching in the reduced search space which is finely
gridded, the accurate location of the targets is found. For estimating the
actual number of targets, a model order selection scheme is used in two cases
of known and unknown noise variance. Also, by analytical investigation of the
probability of detection for the proposed method, the effect of the
influential parameters on the detection ability is explicitly demonstrated.
Compared to the super-resolution methods based on compressed sensing (CS), the
proposed method has a lower computational cost and higher estimation accuracy,
especially at low signal-to-noise ratio regime. Simulation results show the
superiority of the proposed method in terms of both 3D scatterer
reconstruction and detection ability.

###### keywords:

Synthetic aperture radar tomography (TomoSAR), nonlinear least square (NLS),
complexity reduction, detection, model order selection

*Ahmad Naghavi<EMAIL_ADDRESS>
##  1 Introduction

Synthetic aperture radar (SAR) facilitates a day-night all-weather two-
dimensional imaging for earth observation [1]. Space-borne missions such as
TerraSAR-X (TSX) [2] and COSMO-SkyMed [3] capture very high resolution (VHR)
data stacks. Using three-dimensional SAR imaging, the distribution of
scatterers in the elevation direction is also achieved. Synthetic aperture
radar interferometry (InSAR) is a powerful tool for height profile estimation
but it can not resolve the scatterers in height [4].

As a complementary method, SAR tomography (TomoSAR) is capable of resolving
the scatterers through a synthetic aperture in the elevation direction[5]. It
retrieves the location of the scatterers and their reflectivity by processing
a stack of range-azimuth focused images provided by captured data from
slightly different orbits [6]. Monitoring applications of SAR tomography are
in two main categories: natural land coverage and urban infrastructure.
Furthermore, by using multi-temporal SAR acquisition, TomoSAR has been
extended to higher dimensions. Four-dimensional (4D) tomography presents
deformation velocity of the scatterers in addition to elevation [7, 8]. By
five-dimensional (5D) tomography, thermal expansion of the scatterers is also
obtained [9, 10].

The main processing in TomoSAR is the inversion, which is the extraction of
the reflectivity vector of scatterers in elevation direction from the acquired
data. This leads to a spectral estimation problem which can be solved by a
parametric or non-parametric method. The first non-parametric solution is
beamforming [5]. Nonlinear least square (NLS), as a parametric method that is
accurate but is time-consuming [11]. Furthermore, the parametric methods such
as NLS and need knowledge about the number of targets [12] and hence, should
employ a model order selection scheme.

In addition to TomoSAR inversion methods using model order selection,
detection-based methods have been proposed for scatterers detection [13, 14].
Since there are some unknown parameters in the detection process, a common
solution is using the generalized likelihood ratio test (GLRT). In [13], GLRT
is used for single scatterer detection. For the double scatterers case, the
solution called sequential GLRT with cancellation (SGLRTC) is proposed which
is fast but has low resolution[14]. In [15], another GLRT-based method called
Sup-GLRT has been presented. This method is a high-resolution method and can
detect more than two targets at the expense of searching in larger
combinatorial spaces.

Compressed sensing (CS) as a sparse data recovery method, has been widely used
in signal processing [16]. To achieve super-resolution by CS, a limited number
of observations suffices, and knowledge of the number of targets is not
necessary. After compressed sampling, a reconstruction algorithm is needed to
recover the data. Among the several algorithms, $l_{1}$-norm minimization has
received the most attention [17][18]. Despite the high accuracy of
$l_{1}$-norm minimization algorithms, they are computationally complex and
suffer from outliers after reconstruction. [19][20]. In various recent studies
on TomoSAR, $l_{1}$-norm minimization has been applied to obtain height
profile of urban infrastructures [21, 22, 23]. In [23], the Scale-down by L1
norm Minimization, Model selection, and Estimation Reconstruction (SL1MMER)
method was proposed in which, by applying model order selection over the
downsized data resulted by CS reconstruction, the number of targets and their
locations are obtained.

While super-resolution capability is desirable in every imaging system, in
TomoSAR, it is more critical because of the layover phenomenon. When layover
occurs, multiple targets in elevation contribute to the same range-azimuth
cell. By utilizing tomography, layover can be resolved if the minimum distance
between the scatterers is more than one resolution unit in the elevation
direction.

In this paper, we propose a computationally-efficient super-resolution method
that accurately detects the scatterers through two consecutive steps. First,
in the coarse detection step, by exploiting correlation information, the
possible locations of targets are characterized. Then in the fine detection
step, a combinatorial search is applied for locating the targets. Also, in
this step, we used a model order selection scheme for the two cases of known
and unknown noise variance. Detection ability and estimation accuracy of the
proposed method are used as performance criteria and investigated regarding
the related factors. Moreover, by deriving an analytical expression for the
probability of detection of the proposed method, we investigate the impact of
the influential factors on the detection ability.

The rest of this paper is organized as follows. The system model appears in
Section 2. In Section 3, the proposed two-step detection method is presented.
The performance criteria of the proposed method are analyzed in Section 4.
Section 5 presents the results of TomoSAR experiments over simulated SAR data.
Section 6 states the advantages of the proposed method in computational cost.
Finally, section 7 concludes this work.

##  2 System Model

The multipass SAR geometry is shown in Fig. 1. In TomoSAR processing, first a
stack of complex datasets are acquired at $N$ orbit positions (orthogonal
baselines). Then, range-azimuth compression is applied on the raw data. It is
shown that (after some postprocessing), the focused complex image for a range-
azimuth pixel at $n$-th pass is [6]

$g_{n}=\underset{\Delta_{s}}{\int}\boldsymbol{\gamma}(s)\exp(j2{\pi}\xi_{n}s)ds,\hskip
4.26773ptn=0,\ldots,N-1$ (1)

where $\boldsymbol{\gamma}(s)$ is the reflectivity function along the
elevation direction $s$, $\xi_{n}={2b_{n}}/(\lambda R_{0})$ is the spatial
frequency, $b_{n}$ stands for the $n$-th orthogonal baseline relative to the
reference (master) acquisition, $\lambda$ is operating wavelength, $R_{0}$ is
the distance from the master baseline to the center of the ground scene and
$\Delta_{s}$ is the elevation extent. Equation (1) shows a spatial Fourier
transform relationship between the target reflectivity in the elevation
direction and the post-processed complex image at $n$-th pass. Discretizing
(1) along $s$ dimension yields the following linear model

Figure 1: Multipass SAR geometry in the range-elevation $(r,s)$ plane.

$\mathbf{g}=\mathbf{A}\boldmath{{\gamma}}+\mathbf{n}$ (2)

where $\mathbf{g}=\left[g_{1},g_{2},...,g_{N}\right]^{\mathrm{T}}$ is the
vector of the range-azimuth compressed images obtained at all passes and
$\mathrm{(.)^{T}}$ denotes the transpose operator. Let $M$ be the number of
samples (grids) in the elevation direction. $s_{m},m=1,2,...M$ is the
$m^{\prime}$th elevation samples and
$\boldsymbol{\gamma}=\left[\gamma_{1},\gamma_{2},\ldots,\gamma_{M}\right]^{\mathrm{T}}$
is the reflectivity vector. The matrix
$\mathbf{A}=[\boldsymbol{a}(s_{1}),\boldsymbol{a}(s_{2}),...,\boldsymbol{a}(s_{M})]$
is the $N\times{M}$ steering matrix where $\boldsymbol{a}(s_{m})$ is the m’th
steering vector whose $n$’th element is $\exp(j2\pi\xi_{n}s_{m})$ and
$j=\sqrt{-1}$. The noise term $\mathbf{n}$ is a zero-mean circular complex
Gaussian vector with covariance matrix
$E{\\{\mathbf{n}\mathbf{n}^{\mathrm{H}}}\\}=\sigma^{2}_{n}\mathbf{I}_{N}$,
where $\mathbf{I}_{N}$ is $N\times{N}$ identity matrix.

The basic tomographic inversion method is beamforming [5] which is actually
the correlation between $\mathbf{g}$ and the steering matrix $\mathbf{A}$.
Accordingly, Rayleigh resolution is defined as the minimum separation between
two point scatterers which can be discriminated by the beamforming method

$\rho_{s}=\frac{{\lambda}R_{0}}{2\Delta b}$ (3)

where $\Delta b$ is the total baseline extent.

##  3 Two-Step Multiple Scatterers Detection

In SAR tomography, by using an inversion method, all elements of the
reflectivity vector are estimated. However, for point scatterers, the
objective of TomoSAR is only the detection of non-zero elements of the
reflectivity vector. GLRT and SGLRTC are two detectors [13, 14] in which,
detection is based on beamforming method or equivalently correlation
information. Since the number of passes is usually small, the baseline extent
$\Delta b$ is limited, and as a result, the resolution of the SGLRTC detector
is low. A super-resolution method suitable for multiple scatterers detection
is the NLS method, which is based on maximum likelihood (ML) estimation [24].
The NLS method, though an accurate solution, has high computational complexity
due to the required exhaustive search [11].

To reduce the computational cost of NLS, a two-step method is proposed called
Correlation-Aided NLS (CA-NLS). In its first step, coarse detection is
performed through finding possible location of targets. Then by searching in
the reduced space obtained in the first step, fine detection is achieved. In
the following, the two-step method is presented in detail.

###  3.1 Coarse Detection

In the coarse detection step, we use the correlation information in an
iterative routine to characterize the possible location of the targets. To do
this, first the correlation between the received data $\mathbf{g}$ and the
steering matrix $\mathbf{A}$ is calculated. Then, the target is detected by
comparing the maximum power of the correlation vector (after scaling) with a
threshold as [25, 13]

$\max_{s_{m}}\frac{|\boldsymbol{a}^{H}(s_{m})\mathbf{g}|^{2}}{N\|\mathbf{g_{\perp}}\|^{2}}\underset{\mathcal{H}_{0}}{\overset{\mathcal{H}_{1}}{\gtrless}}T,\hskip
8.5359ptm=1,2,...M$ (4)

where $T$ is a threshold adjusted based on a constant false alarm rate (CFAR)
scheme and $\mathbf{g_{\perp}}$ is the residual component of $\mathbf{g}$ that
is orthogonal to $\boldsymbol{a}(s_{m})$

$\mathbf{g_{\perp}}=\mathbf{g}-\boldsymbol{a}(s_{m})[\boldsymbol{a}(s_{m})^{H}\boldsymbol{a}(s_{m})]^{-1}\boldsymbol{a}(s_{m})^{H}\mathbf{g}.$
(5)

For the multi-target case, the location of the detected target obtained by (4)
does not necessarily coincide with their real location. However, the vicinity
area of the detected peak can be characterized as the possible location of a
target. This area has the extent of $2\rho_{s}$ and is centered around the
detected peak.

The procedure of coarse detection is shown in Algorithm 1 which includes two
consecutive loops of iteration and decision. In the iteration loop, the index
of the dominant components is found using (4), and then, by canceling the
detected dominant component, the residual term is updated. This routine which
resembles orthogonal matching pursuit (OMP) [26], is successively done until
the maximum number of targets $k_{\mathrm{max}}$ are detected. The updating
parameters at $k$-th iteration are as follows: $\mathbf{r}_{k}$ is the
residual term, $p_{k}$ is the index of the peak point of correlation vector
between $\mathbf{r}_{k-1}$ and $\mathbf{A}$, $\Lambda_{k}$ is the vector of
detected peaks indexes, $\Gamma_{k}$ is the test statistics of (4) and
$supp_{k}$ is the partial support related to the detected peak. In the
decision loop, the derived test statistics are compared with the threshold $T$
in reverse order. If in the $k$’th loop $\Gamma_{k}$ exceeds the threshold,
$S$ which is the possible support of the targets will be the union of all $k$
detected partial supports.

Algorithm 1 Coarse Detection

1: $\mathbf{Input}:\mathbf{g},\mathbf{A},\rho_{s},k_{\mathrm{max}},T$

2: Initialization:
$k=0,\mathbf{r}_{0}=\mathbf{g},\Lambda_{0}=\emptyset,supp_{0}=\emptyset$

3: – Iteration loop

4: while $\hskip 4.26773pt(k\textless k_{\mathrm{max}})$ do

5: $k=k+1$

6:
$p_{k}=arg\max_{i}{\boldsymbol{|}\boldsymbol{a}^{H}(s_{i})\mathbf{r}_{k-1}\boldsymbol{|}},\hskip
4.26773pti=1,2,...,M$

7: $\Lambda_{k}=\Lambda_{k-1}\bigcup{p_{k}}$

8:
$\mathbf{r}_{k}=\mathbf{g}-\mathbf{A}_{\Lambda_{k}}(\mathbf{A}^{H}_{\Lambda_{k}}\mathbf{A}_{\Lambda_{k}})^{-1}\mathbf{A}^{H}_{\Lambda_{k}}\mathbf{g}$

9:
$\Gamma_{k}=\frac{\boldsymbol{|}\boldsymbol{a}^{H}(s_{p_{k}})\mathbf{r}_{k-1}\boldsymbol{|}^{2}}{N.\|\mathbf{r}_{k}\|^{2}}$

10:
$supp_{k}=\boldsymbol{[}s_{p_{k}}-\rho_{s},s_{p_{k}}+\rho_{s}\boldsymbol{]}$

11: end while

12: – Decision loop

13: while $k>0$ do

14: if $\Gamma_{k}>T$ then

15: break

16: else

17: $k=k-1$

18: end if

19: end while

20: $\mathbf{Output}$ : $S=\bigcup\limits^{k}_{i=0}{supp_{i}}$

For the special case of $k_{\mathrm{max}}=2$, the coarse detection algorithm
presented is similar to the SGLRTC method [14]. In our method, however, the
support information acquired during the coarse detection step, enables fast
detection of targets in the next step.

###  3.2 Fine detection

After identification of the possible location of targets, it is desired to
find the accurate location of targets. Assumming presence of $k$ targets,
$\Omega_{k}$ is defined as support of $k$ targets (position of nonzero
elements of reflectivity vector), $\boldsymbol{\gamma}_{\Omega_{k}}$ as the
related $k$ nonzero elements of the reflectivity vector and
$\mathbf{A}_{\Omega_{k}}$ is a $N\times{k}$ matrix composed of the $k$ columns
of matrix $\mathbf{A}$ whose indexes correspond to ${\Omega_{k}}$. The NLS
estimate of $\Omega_{k}$ and $\boldsymbol{\gamma}_{\Omega_{k}}$ are as follows
[24]

$\hat{\Omega}_{k}=arg\min_{\Omega_{k}\in{S}}\mathbf{g}^{H}\Pi^{\perp}_{\Omega_{k}}\mathbf{g}$
(6)
$\hat{\boldsymbol{\gamma}}_{\Omega_{k}}=(\mathbf{A}^{H}_{\Omega_{k}}\mathbf{A}_{\Omega_{k}})^{-1}\mathbf{A}^{H}_{\Omega_{k}}\mathbf{g}.$
(7)

where
$\Pi^{\perp}_{\Omega_{k}}=\mathbf{I}_{N}-\mathbf{A}_{\Omega_{k}}(\mathbf{A}^{H}_{\Omega_{k}}\mathbf{A}_{\Omega_{k}})^{-1}\mathbf{A}^{H}_{\Omega_{k}}$
is the orthogonal complement projection onto the subspace spanned by
$\mathbf{A}_{\Omega_{k}}$. The combinatorial search in (6) is done in $S$
which is the possible location of targets obtained in previous step. This
search space reduction makes accurate localization of targets feasible. For
determining the number of targets, the model order slection is required which
is done by optimizing a penalized likelihood criteria. Depending on noise
variance being known or unknown, the model can be selected in two ways that is
presented in the unified form as follows

$\hat{k}=arg\min_{k}\\{J_{k}\triangleq f(\varepsilon(k))+\mathcal{P}(k)\\}$
(8)

where $\varepsilon(k)$ and $f(x)$ are definde as

$\varepsilon(k)=\min_{\Omega_{k}\in{S}}\\{\mathbf{g}^{H}\Pi^{\perp}_{\Omega_{k}}\mathbf{g}\\}.$
(9) $f(x)=\left\\{\begin{array}[]{ll}\frac{x}{\sigma^{2}_{n}},\hskip
53.77568ptknown\hskip 5.69046pt\sigma^{2}_{n}\\\ \\\ N\ln(\frac{x}{N}),\hskip
23.61583ptunknown\hskip 5.69046pt\sigma^{2}_{n}.\end{array}\right.$ (10)

and $\mathcal{P}(k)$ is the penalty term. For the sake of brevity, the details
on derivation of (6)-(10) are refered in Appendix A. The flowchart of the
proposed method is shown in Fig. 2.

There are various penalty terms, among them, Akaike information criterion
(AIC) [27] and the Bayesian information criterion (BIC) [28] are commonly
used. AIC and BIC belong to the family of information-theoretic criteria with
the following penalty term [29]

$\mathcal{P}(k)=\eta.K$ (11)

where $\eta$ is the penalty factor and $K$ is the number of unknown
parameters. For AIC and BIC, the values of the penalty factor are
$\eta_{AIC}=1$ and $\eta_{BIC}=0.5\ln{N}$, respectively. Due to model
selection bias of AIC for finite samples, its corrected version (AICc) is also
suggested where

Figure 2: Flowchart of the proposed two-step detector (CA-NLS).

$\eta_{AICc}=\frac{N}{N-K-1}$ [30]. Since there are three unknown parameters
for every target (location, amplitude, and phase), the total number of unknown
parameters is $K=3k$.

Overall, in the fine detection step, first, a combinatorial search is
performed for each assumed $k$. Then, the model order is selected via
minimization of the cost function (8) over $k\in[0,k_{\mathrm{max}}]$.
However, if $J_{k}<J_{k+1}$, $k$ will be selected as the number of targets,
and no further search is required for the higher number of targets.

##  4 Performance criteria

In this section, the performance of the CA-NLS method in multiple scatterers’
detection is analyzed. Then, through the simulation of system model (2), the
efficiency of CA-NLS is compared with the other methods including SGLRTC,
SL1MMER, and Sup-GLRT. The first comparing method is a low-resolution method,
while the second and third are high-resolution ones. The performance of the
methods has been evaluated by two criteria: detection ability and estimation
accuracy.

In this section, we consider the presence of a maximum of two targets per a
range-azimuth cell ($k_{\mathrm{max}}=2$). SNR per scatterer is defined as the
square of the amplitude of a scatterer to the power of the noise. In double
scatterers scenario, amplitudes of the two scatterers are considered the same
and equal to $\sigma_{s}$. Thus, it is defined as
$SNR=\frac{\sigma^{2}_{s}}{\sigma^{2}_{n}}$. In addition, the distance of
targets is shown by $D_{s}$ and the normalized distance of targets $\alpha$ is
defined as $\alpha=\frac{D_{s}}{\rho_{s}}$. The parameters chosen for the
simulation are the number of passes $N=20$, baseline extent $\Delta b=903$ m
resulting in Rayleigh resolution $\rho_{s}=26$ m, elevation extent
$\Delta_{s}=360$ m, number of grids $M=234$ or equivalently 17 points per
$\rho_{s}$ and the threshold $T=0.8$ (to have probability of false alarm
$P_{FA}=10^{-3}$ for SGLRTC method).

###  4.1 Detection ability

For the detection study, $\mathcal{H}_{i}$ is defined as the hypothesis of the
presence of $i$ targets in the elevation direction and $\mathcal{D}_{k}$ as
decision about presence of $k$ targets. The probability of detection for the
resolution problem is defined as

$P_{D}=P(\mathcal{D}_{2}|\mathcal{H}_{2}).$ (12)

Let’s define $P_{D1}$ and $P_{D2}$ as the probabilities of detection of the
coarse and fine detectors, respectively. As a cascaded detector, it can be
shown that the total probability of detection will be the minimum of the two
probabilities of detection. Hence for an appropriate threshold-setting, it is
recommended to adjust $T$ in the first detector and penalty factor $\eta$ in
the second detector to ensure $P_{D1}>P_{D2}$. Otherwise, missed detection in
the coarse detection step may result in the null $S$ and the detection will
not proceed to the fine step.

To analytically assess the performance of the proposed two-step method,
$P_{D1}$ and $P_{D2}$ should be derived. For the sake of simplicity, we
emphasize on the derivation of $P_{D2}$ and leave the exact derivation of
$P_{D1}$ for future works. However, an upper-bound has been derived for
$P_{D1}$ in [14] which is valid only for faraway targets. In the rest of this
section, to consider only the detection ability of the fine detector, the
probability of detection of the coarse detector is assumed to be one i.e.
$P_{D1}=1$ which yields $P_{D}=P_{D2}$.

According to the definition of probability of detection in (12) and the model
order selection procedure in (8), the detection is successful when the cost
function of the two-targets case is the least one, i.e. $J_{2}<J_{1}$. By
using definition of $J_{n}$ in (8) for known noise variance case, we have

$\frac{{\|\mathbf{g}-\mathbf{A}_{\Omega_{2}}\gamma_{\Omega_{2}}\|}^{2}}{\sigma_{n}^{2}}+\mathcal{P}(2)<\frac{{\|\mathbf{g}-\mathbf{A}_{\Omega_{1}}\gamma_{\Omega_{1}}\|}^{2}}{\sigma_{n}^{2}}+\mathcal{P}(1).$
(13)

where each side of (13) composed of a likelihood term and a penalty term.

On the left-hand side of (13), in the likelihood term, the location of two
steering vectors is searched where $\mathbf{g}$ is composed of the noise plus
two steering vectors. Under good SNR conditions, the estimated location of
targets approximately coincides with their real value. Hence, after the
cancellation of steering vectors of the estimated targets, the likelihood term
on the left-hand side of (13) becomes the power of a pure Gaussian noise
scaled to the noise variance. Thus, the likelihood term on the left-hand side
of (13) will have a central chi-square distribution with $N$ degrees of
freedom

$\frac{1}{\sigma_{n}^{2}}{\|\mathbf{g}-\mathbf{A}_{\Omega_{2}}\gamma_{\Omega_{2}}\|}^{2}\approx\frac{1}{\sigma_{n}^{2}}{\|\mathbf{n}\|}^{2}\sim\chi_{N}^{2}.$
(14)

On the right hand side of (13), in the likelihood term, $\mathbf{g}$ is
composed of the noise plus two steering vectors which is subtracted from one
steering vector related to the falsely detected target. Consequently, the
likelihood term on the right hand side of (13) becomes the power of a pure
Gaussian noise plus an additional term scaled to the noise variance. As a
result, the likelihood term on right hand side of (13) will have a noncentral
chi-square distribution with $N$ degrees of freedom

$\frac{1}{\sigma_{n}^{2}}{\|\mathbf{g}-\mathbf{A}_{\Omega_{1}}\gamma_{\Omega_{1}}\|}^{2}\approx\frac{1}{\sigma_{n}^{2}}{\|\mathbf{n+r}\|}^{2}\sim\chi_{N}^{2}(\lambda_{r})$
(15)

where $\mathbf{r}$ is residual term added to Gaussian noise and $\lambda_{r}$
is noncentrality of the above chi-square distribution which are defined
respectively as follows

$\mathbf{r}=\sigma_{s}\boldsymbol{[}\boldsymbol{a}(s_{M_{1}})+e^{j\Delta{\phi}}\boldsymbol{a}(s_{M_{2}})-\mathbf{\hat{\gamma_{p}}}\boldsymbol{a}(s_{p})\boldsymbol{]}$
(16) $\lambda_{r}=\frac{\|\mathbf{r}\|^{2}}{\sigma_{n}^{2}}$ (17)

where $s_{M_{1}}$ and $s_{M_{2}}$ are real locations of targets,
$\Delta{\phi}$ is their phase difference, $s_{p}$ is the estimated single-
target location and $\mathbf{\hat{\gamma_{p}}}$ is the estimated reflectivity
of the single-target scaled to $\sigma_{s}$. By rearranging (13) and using
(14),(15), we obtain

$\frac{1}{\sigma_{n}^{2}}({{\|\mathbf{n+r}\|}^{2}}-{\|\mathbf{n}\|}^{2})>\mathcal{P}(2)-\mathcal{P}(1).$
(18)

According to the ITC penalty in (11), the above inequality can be rewritten as

$\frac{1}{\sigma_{n}^{2}}({{\|\mathbf{n+r}\|}^{2}}-{\|\mathbf{n}\|}^{2})>3\eta.$
(19)

The following proposition shows that the left hand side of (19) has a Gaussian
distribution.

###### Proposition 1

: Let $X$ be a random variable with central chi-square distribution
$\chi^{2}_{N}$ and $Y$ a random variable with noncentral chi-square
distribution $\chi^{2}_{N}(\lambda_{r})$ with $N$ degrees of freedom, both
built by a common circular Gaussian random vector $\mathbf{n}$. The
distribution of difference variable $Z=Y-X$ will be Gaussian with the mean and
variance values equal to $\lambda_{r}$ and $2\lambda_{r}$, respectively.

######  Proof 4.1.

See Appendix B.

As a result, $P_{D}$ can be derived as

$\displaystyle P_{D}$ $\displaystyle=P(Z>3\eta)$
$\displaystyle=\frac{1}{\sqrt{2\pi\sigma^{2}_{Z}}}\int_{3\eta}^{\infty}exp(-\frac{(Z-m_{Z})^{2}}{2\sigma^{2}_{Z}})dZ$
$\displaystyle=Q(\frac{3\eta-m_{Z}}{\sigma_{Z}})$
$\displaystyle=Q(\frac{3\eta}{\sqrt{2\lambda_{r}}}-\sqrt{\frac{\lambda_{r}}{2}})$
(20)

where $Q(x)=(1/\sqrt{2\pi})\int_{x}^{\infty}exp^{-t^{2}/2}dt$ is the
$Q$-function. Also, the noncentrality term in (17) can be expressed in terms
of other parameters as

$\lambda_{r}=N.SNR.\vartheta(\alpha,\Delta{\phi})$ (21)

where $\vartheta(.)$ is a function of normalized distance of two targets
$\alpha$ and their phase difference $\Delta{\phi}$. Derivation of $s_{p}$,
$\mathbf{\hat{\gamma_{p}}}$ and $\lambda_{r}$ is given in Appendix C.

Figure 3: $P_{D}$ versus SNR for CA-NLS method respecting the penalty term
used for model order selection. The normalized distance of targets is
$\alpha=\frac{1}{2}$.

Due to the complexity of $\vartheta(.)$, this two-dimensional function is
obtained by numerical calculation whose result is sketched in Fig. 12. It can
be observed that for each value of $\Delta{\phi}$, the function
$\vartheta(\alpha)$ is increasing with respect to $\alpha$ where
$0<\alpha<1.34$. This range of $\alpha$ is reasonable as the relationships led
to (21) are valid for closely-spaced targets. On the other hand, the effect of
$\Delta{\phi}$ on $\vartheta(.)$ is also demonstrated in Fig. 12. By looking
at the results, it is seen that the least probability of detection occurs in
$\Delta{\phi}=0$.

According to (21), by increasing $N.SNR$ and $\alpha$, the noncentrality
$\lambda_{r}$ increases. Since $Q(.)$ is a decreasing function with respect to
its argument, increasing of $\lambda_{r}$ leads to increase of $P_{D}$. As a
result, $N.SNR$ and $\alpha$ have increasing effects on $P_{D}$. The other
influential parameter on $P_{D}$ is the penalty factor $\eta$. From (21), it
is evident that increase of $\eta$ leads to decrease of $P_{D}$.

In the first numerical experiment, different cases of penalized cost function
(8) have been compared. These cases are related to the availability of the
information about the noise variance (known or unknown) and the penalty term
used (BIC, AIC, and AICc). In Fig. 3, $P_{D}$ of the different cases is
sketched versus SNR. Also for these cases, the average probability of false
detection $P_{FD}$ (i.e. $P(\mathcal{D}_{2}|\mathcal{H}_{1})$) is measured and
has been indicated in the legend of the figure.

Figure 4: $P_{D}$ as a function of SNR for SGLRTC (green circles), SL1MMER
(red squares), CA-NLS (blue stars) and analytical approximation (black
diamonds), for $N=20$ and normalized distance of targets $\alpha=\frac{1}{2}$.

It is observed that the probability of detection in the case of unknown noise
variance is slightly superior to the known noise variance one, but its average
$P_{FD}$ is also higher. Besides, for AIC, $P_{D}$ and average $P_{FD}$ are
higher comparing to BIC. For localization applications, when $P_{FD}$
increases, most of the single point scatterers are detected as two-point
scatterers which makes the final three-dimensional SAR image different from
reality. As a result, to maintain image quality, a trade-off between $P_{D}$
and $P_{FD}$ should be considered in designing the detector.

For different methods, $P_{D}$ is shown versus SNR and $\alpha$ in Figs. 4 and
5, respectively. Moreover, the analytical approximation of $P_{D}$ for CA-NLS
(mentioned in (20)-(21)) has been plotted. For a fair comparison with the
simulated methods using random $\Delta{\phi}$, the analytical function
$\vartheta(\alpha,\Delta{\phi})$ is integrated over
$\Delta{\phi}\in{[-\pi,\pi]}$ .To select the model order, the known noise
variance case with BIC penalty term is used. Also, for the SL1MMER method,
unconstraint form of Basis Pursuit DeNoising (BPDN) is implemented [17] as
follows

Figure 5: $P_{D}$ as a function of normalized distance of targets $\alpha$ for
SGLRTC (green circles), SL1MMER (red squares), CA-NLS (blue stars) and
analytical approximation (black diamonds), for $N=20$ and $SNR=9\hskip
2.84544ptdB$.

$\hat{\mathbf{\gamma}}=arg\min_{\mathbf{\gamma}}\\{\|\mathbf{g}-\mathbf{A}\mathbf{\gamma}\|_{2}^{2}+\lambda_{k}\|\gamma\|_{1}\\}$
(22)

where $\lambda_{k}$ is an adjusting parameter which trades off between
reconstruction fidelity and sparsity. The convex minimization problem in (22)
was solved by the Matlab CVX package [31].

As expected, the results show that the probability of detection rises by
increasing the SNR and distance of targets. It is observed that CA-NLS has the
highest probability of detection, SL1MMER is close to but lower than the CA-
NLS, and SGLRTC has the least $P_{D}$. Also, average $P_{FD}$ of the methods
is obtained as 0.002, 0.03, and 0.05 for SGLRTC, CA-NLS, and SL1MMER
respectively. These results show the advantage of CA-NLS to SL1MMER in having
both higher $P_{D}$ and lower $P_{FD}$. Furtherly, when SNR is high enough
($SNR>9$ dB), the analytical expression (20) is close to the numerical result
of CA-NLS.

Figure 6: RMSE of estimated locations of targets versus $SNR$ for SGLRTC
(green circles), SL1MMER (red squares), CA-NLS (blue stars) and square root of
CRLB (black diamonds), for $N=20$ and normalized target’s distance
$\alpha=\frac{1}{2}$.

###  4.2 Estimation accuracy

For evaluation of the estimation accuracy, the location of targets is
considered where the reflectivity is not of concern for localization purpose.
To have a criterion for the estimatin accuracy, Cramer Rao lower bound (CRLB)
is used as an asymptotic bound for variance of an unbiased estimator. In SAR
tomography, various studies have been conducted for obtaining CRLB [13, 23].
Associated CRLB for the one-target and two-target cases are respectively as
follows [32, 33]

$\hskip
31.2982ptCRLB_{1}=\frac{3}{2\pi^{2}}(\frac{\rho_{s}}{\sqrt{N.SNR}})^{2},\hskip
29.30629ptk=1\\\ $ (23) $\hskip 31.2982ptCRLB_{2}\approx
CRLB_{1}.\upzeta(\alpha,\Delta{\phi}),\hskip 27.03003ptk=2\\\ $ (24)

where $\upzeta(.)$ is the normalized CRLB which is a function of $\alpha$ and
$\Delta{\phi}$ with a value more than one such that $CRLB_{2}>CRLB_{1}$ [33,
23]. Also it is shown that the most estimation error occures at
$\Delta{\phi}=0$ [23]. Since multiple snapshots are available in the real
scenario, the phase difference of scatterers becomes random in nature. In this
case, a $\Delta{\phi}$-independent expression is derived for normalized CRLB
[32, 33]

$\upzeta(\alpha)\approx\max\\{\frac{15}{\pi^{2}}\alpha^{-2},1\\}$ (25)

Figure 7: RMSE of estimated locations of targets versus normalized distance of
targets $\alpha$ for SGLRTC (green circles), SL1MMER (red squares), CA-NLS
(blue stars) and square root of CRLB (black diamonds), for $N=20$ and
$SNR=9\hskip 2.84544ptdB$.

In Fig. 6, RMSE of different methods is sketched versus SNR. To measure the
estimation error, root mean square error (RMSE) of the location of point
targets is considered

$RMSE=\sqrt{\frac{1}{N_{Det}}\sum\limits_{i=1}^{N_{Det}}{\frac{(s_{M_{1}}-\hat{s}_{M_{1i}})^{2}+(s_{M_{2}}-\hat{s}_{M_{2i}})^{2}}{2}}}$
(26)

where $N_{Det}$ is the total number of the successful detections,
$s_{M_{1}},s_{M_{2}}$ are real positions of targets in elevation and
$\hat{s}_{M_{1i}},\hat{s}_{M_{2i}}$ are estimated positions of two targets at
$i$-th successful detection. It should be noted that in the following results,
RMSE of different methods together with the square root of CRLB is scaled to
Rayleigh resolution $\rho_{s}$. It is seen that CA-NLS has the least
estimation error and the error decreases by increasing SNR. For SGLRTC, it has
a nearly constant error versus SNR which is higher than that of other methods.
Similar to CA-NLS, the error of SL1MMER decreases by increasing SNR which is
close to that of CA-NLS. However, the SL1MMER method has a high estimation
error at lower SNR values. Hence, adequate SNR is required for using the
SL1MMER method, while CA-NLS performs well even at low SNR conditions.

In the next experiment, RMSE is measured versus normalized distance of targets
$\alpha$. As seen in Fig .7, the estimation error decreases by increasing the
distance of targets. The performance of CA-NLS is superior to those of SL1MMER
and SGLRTC. The error of CA-NLS and SL1MMER decreases by increasing $\alpha$
and approach zero, but for SGLRTC the error is nearly constant for
$\alpha<1.5$.

##  5 Experiments on simulated SAR data

In this section, the performance of different detection methods is measured
using simulated SAR data. To do this, two types of experiments have been
performed: layover separation and three-dimensional object reconstruction. Due
to working in the real scenario, when the noise variance is not known, the
unknown noise variance case is used in (8) for model selection.

###  5.1 Layover separation

In this experiment, two scatterers are considered, which are located on the
same range-azimuth cell, one on the ground and the other on the building
facade. The aim is to measure the separability of the targets when their
spacing varies from low to high. This type of evaluation has been previously
used in [21, 34]. The simulation parameters are shown in Table 1. For
preparing the tomographic stack, Range-Doppler Algorithm (RDA) has been used
for focusing SAR images. Next, coregistration and deramping have been applied
to the focused images.

Table 1: SAR Simulation Parameters. Parameter | Value
---|---
Platform velocity (m/s) | 200
Carrier frequency (GHz) | 4.5
Pulse repetition frequency (PRF) (Hz) | 300
Range distance to center of scene (km) | 20
Number of baselines (N) | 24
SNR (dB) | 9
Figure 8:  Estimated elevations of double scatterers located on ground and
facade at a linear increasing spacing. The sub-figures (a), (c) are for
SGLRTC, (b), (d) are for SL1MMER and (c), (e) are for CA-NLS. The top row
indicates targets with equal phase $\Delta{\phi}=0$ and the bottom row shows
targets with unequal phases $\Delta{\phi}\neq 0$. (Green and blue points
indicate detected double scatterers and the red points indicate unresolved
single scatterers).

The simulation result is shown in Fig. 8. The horizontal axis shows targets on
the ground, and the diagonal line belongs to the second targets on the facade.
Two parallel dashed lines (margin lines) show $\pm 3$ times of square root of
$CRLB_{1}$. The top subfigures are results of experiments done under
$\Delta{\phi}=0$ condition, and the bottom ones are the result of the
uniformly random phase difference ones. In this figure, the number of red
points indicates the number of missed detections. The deviation of detected
points (blue or green) from margin lines indicates the amount of estimation
error.

Figure 9: Three-dimensional reconstruction of a synthesized structure by
tomographic inversion, (a) Real points (b) SGLRTC (c) SL1MMER and (d) CA-NLS.
Figure 10: Side View of reconstructed images by TomoSAR methods, (a) Real
points in which double scatterers and triple scatterers zones are indicated
(b) SGLRTC (c) SL1MMER and (d) CA-NLS. Single, double and triple scatterer are
indicated by red, blue, and green dots, respectively.

The results of this experiment show the superiority of CA-NLS and SL1MMER
methods as they can resolve targets even at the lower distance of targets (the
area close to the origin). Also, SGLRTC has the highest probability of missed
detection and estimation error. By comparing the top and bottom sub-figures,
we conclude that the non-zero phase difference increases the probability of
detection (fewer red points) and decreases the estimation error. This verifies
the role of $\Delta{\phi}$ in $CRLB_{2}$ (24) and the analytical expression of
$P_{D}$ (20)-(21). Finally, the experiment took $0.3$, $13.8$, and $524$
seconds for SGLRTC, CA-NLS, and SL1MMER, respectively, which shows a
significant advantage of CA-NLS to SL1MMER in terms of complexity reduction.

###  5.2 3D target reconstruction

In this experiment, a three-part structure including facade, roof, and the
ground is considered. The 3D structure is shown in Fig. 9 together with the
results of TomoSAR reconstruction by SGLRTC, CA-NLS, and SL1MMER. The results
show that CA-NLS performs better in reconstruction relative to the SGLRTC,
especially at ground-facade and facade-roof conjunctions. In Fig. 10, side
view images of reconstructed structure is shown. Although by SL1MMER the
structure is formed, its parts are thickened and messy. In general, CA-NLS
shows the highest performance qualitatively.

To have a quantitative comparison, we analyzed the detection ability of the
methods. For constructing ground, facade, and roof, a total of 1905 point
targets were used, including 390 double scatterers (facade-ground) and 375
triple scatterers (roof-facade-ground) as indicated in Fig. 10. According to
the results, it is observed that for the existing 390 double scatterers, CA-
NLS, SL1MMER, and SGLRTC have detected 378, 244, and 326 double scatterers,
respectively. These statistics show that the double scatterers are better
detected by CA-NLS. Also, the total elapsed time of the experiment was $0.5$,
$290$, and $420$ seconds for SGLRTC, CA-NLS, and SL1MMER, respectively, which
shows the advantage of CA-NLS to SL1MMER in terms of computational cost.

Fig. 11 shows the reconstruction results for CA-NLS regarding the penalty
criteria used for model order selection. As the maximum number of targets is
three $(k_{\mathrm{Max}}=3)$, the model selection scheme should compromise
between the higher probability of detection and lower probability of false
detection. As seen in Fig. 11, the quality of the AICc scheme is superior to
that of AIC and BIC. The reason is the higher value of the penalty term for
AICc compared with AIC and BIC which leads to less model overfitting. In other
words, the probability of false detection, i.e.
$P(\mathcal{D}_{3}|\mathcal{H}_{2})$ for AICc is lower comparing to AIC and
BIC. Hence the parts of the structure in which two scatterers are present at a
range-azimuth pixel are better reconstructed by AICc, and consequently, a tidy
image has resulted. According to the results, for the existing 390 double
scatterers, BIC, AIC, and AICc have detected 278, 198, and 378 double targets,
respectively. Also, for the 375 triple scatterers existed, BIC, AIC, and AICc
have detected 474, 562, and 358 triple targets, respectively. It approves that
AICc is the most suitable penalty criteria for the CA-NLS method.

Figure 11: Comparison of model selection criteria used in CA-NLS for TomoSAR
3D reconstruction: (a) BIC (b) AIC and (c) AICc.

Although the effect of deformation and thermal expansion has not been
investigated in this paper, the performance of such exact modeling is
predictable. By adding more dimensions, the layovered scatterers can be better
resolved where the values of deformation velocity or thermal dilation of the
two targets are likely different [35]. As a result, in most cases for 3D
reconstruction, the fine detection step is not mandatory. The reason is the
two close point targets obtained the coarse detection step are separable due
to being located in different deformation or thermal bins. In [36], the
authors have used a non-accurate fast version of Sup-GLRT that improves the
detection performance by incorporating additional dimensions. However, some
considerations should be given to the threshold setting in the
multidimensional scenario.

##  6 Computational complexity

A clear advantage of the CA-NLS is the reduction of the combinatorial search
needed in the fine detection step. With the aid of coarse detection, the
number of points (grids) in which, the $k$ targets should be searched, reduces
from $M$ to $M_{r}$ and consequently, the number of likelihood function
evaluations decreases from $M\choose k$ to $M_{r}\choose k$.

In Table 2, the elapsed time of different methods is shown for various numbers
of grids. The normalized distance of targets is assumed $\alpha=\frac{1}{2}$.
Simulations are performed in MATLAB R2017a environment, using Intel Core (TM)
2 Duo 7700k, 4.27 GHz. Choosing a denser grid does not affect the SGLRTC
method but increases the time consumption of the Sup-GLRT method dramatically.
For CA-NLS, using the coarse detection technique makes the computational
complexity tolerable. Additionally, Sup-GLRT has another shortcoming that is
the need for $k_{\mathrm{max}}$-dimensional search at each step of sequential
detection, while for CA-NLS it is just necessary at the final step of model
order selection. In other words, in the no-target case, the two-targets search
(for finding $\Omega_{2}$) is not necessary for CA-NLS, while it is a must for
Sup-GLRT.

For the closely spaced targets, the computational complexity of SL1MMER is
more than that of CA-NLS. Also, for the SL1MMER method, applying compressed
sensing is always obligatory even in the absence of any coherent target(s),
where CA-NLS can easily reject the presence of target at coarse detection
step.

Table 2: Comparison of elapsed time (msec) for different methods. Number of grids | SGLRTC | CA-NLS | SL1MMER | Sup-GLRT
---|---|---|---|---
100 | 0.3 | 5 | 455 | 337
200 | 0.3 | 9 | 467 | 1938
300 | 0.3 | 16 | 484 | 3642

Beside search space reduction gained by CA-NLS, it has another advantage in
complexity reduction. For the case when the targets are far apart, it can be
shown that the fine detection step is not necessary. In other words, the
location of targets estimated in the fine detection step is identical to the
location of detected peaks acquired in the coarse detection step, as the
following proposition states

######  Proposition 1.

Assume $supp_{i},0<{i\leq{k}}$ are the $k$ partial supports detected in the
coarse detection step which are dis-contiguous i.e.
$supp_{i}\bigcap{supp_{j}}=\emptyset,0<{i,j\leq{k}},i\neq{j}$. Then, the
result of fine detection coincides with the associated peaks of detected
partial supports in the coarse detection, i.e.
$\Omega_{k}=\\{s_{p_{1}},s_{p_{2}},...,s_{p_{k}}\\}$.

######  Proof 6.1.

See Appendix D.

In other words the need for costly methods such as NLS and CS is eliminated
when the targets are widely separated. It makes the proposed method
computationally scalable regarding the spacing between the targets.

##  7 Conclusion

In this paper, a high-resolution target detection method has been addressed
for SAR tomography. The proposed two-step method utilizes the reduced search
space derived from the correlation information in the first step for finding
the accurate location of the targets in the second step. Also, an ITC based
model order selection scheme has been used in the second step in two cases of
known and unknown noise variance. The performance of the CA-NLS method was
analyzed and investigated regarding detection ability, estimation accuracy,
and computational cost. We have shown that our reduced-complexity method has
the least estimation error and possesses the highest probability of detection
with an acceptable probability of false detection. Additionally, we have
demonstrated that the phase difference of the two targets is an effective
parameter for detection ability other than $N$ and $SNR$. Furthermore, for far
enough targets, we have shown that coarse detection suffices for target
localization, and an accurate search is unnecessary. Overall, exploiting both
CFAR detection and penalized likelihood criterion in the proposed two-step
method gives a good compromise between $P_{D}$, $P_{FD}$, and the speed of
detection. Further works can be focused on reaching a closed-form analytical
solution for the case of closely-spaced targets, making the real-time
detection of all kinds of targets possible, regardless of their spacing.

##  Appendix A NLS estimation and model order selection

###  A.1 NLS estimation

Assuming the additive noise in (2) is complex Gaussian, the likelihood
function of the observed data in the presence of $k$ targets is defined as

$p(\mathbf{g}|\mathbf{\boldsymbol{\gamma}},k)=\frac{1}{\pi^{N}(\sigma_{n}^{2})^{N}}exp(-\frac{{\|\mathbf{g}-\mathbf{A}_{\Omega_{k}}\boldsymbol{\gamma}_{\Omega_{k}}\|}^{2}}{\sigma_{n}^{2}})$
(27)

The locations of the targets and their reflectivities are estimated by
maximizing the log-likelihood function as follows

$(\hat{\Omega}_{k},\hat{\boldsymbol{\gamma}}_{\Omega_{k}})=arg\max_{\Omega_{k}}\\{{\ln\,p(\mathbf{g}|\boldsymbol{\gamma},k)}\\}$
(28)

Rewriting (28) using (27) gives

$\displaystyle(\hat{\Omega}_{k},\hat{\boldsymbol{\gamma}}_{\Omega_{k}})$
$\displaystyle=arg\min_{\Omega_{k}}\\{-\ln{p(\mathbf{g}|\boldsymbol{\gamma},k)}\\}$
$\displaystyle=arg\min_{\Omega_{k}}\\{N\ln{2\pi}+N\ln{\sigma_{n}^{2}}+\frac{{\|\mathbf{g}-\mathbf{A}_{\Omega_{k}}\boldsymbol{\gamma}_{\Omega_{k}}\|}^{2}}{\sigma_{n}^{2}}\\}$
$\displaystyle=arg\min_{\Omega_{k}}\\{{\|\mathbf{g}-\mathbf{A}_{\Omega_{k}}\boldsymbol{\gamma}_{\Omega_{k}}\|}^{2}\\}$
(29)

which, for a given $\Omega_{k}$, results in the least square estimate of
reflectivity vector as [24]

$\hat{\boldsymbol{\gamma}}_{\Omega_{k}}=(\mathbf{A}^{H}_{\Omega_{k}}\mathbf{A}_{\Omega_{k}})^{-1}\mathbf{A}^{H}_{\Omega_{k}}\mathbf{g}.$
(30)

In order to find the location of targets,
$\hat{\boldsymbol{\gamma}}_{\Omega_{k}}$ is substituted in
${\|\mathbf{g}-\mathbf{A}_{\Omega_{k}}\boldsymbol{\gamma}_{\Omega_{k}}\|}^{2}$
and $\hat{\Omega}_{k}$ is obtained as

$\hat{\Omega}_{k}=arg\min_{\Omega_{k}}\mathbf{g}^{H}\Pi^{\perp}_{\Omega_{k}}\mathbf{g}$
(31)

###  A.2 Model order selection: known and unknown $\sigma_{n}^{2}$

To find the number of targets, the following penalized likelihood criterion
should be optimized

$\hat{k}=arg\min_{k}\\{-\ln{p(\mathbf{g}|\hat{\boldsymbol{\gamma}}_{\Omega_{k}},k)}+\mathcal{P}(k)\\}$
(32)

where $\mathcal{P}(k)$ is the penalty term. Substituting the likelihood
function (27) into (32) gives

$\hat{k}=arg\min_{k}\\{N\ln{2\pi}+N\ln{\sigma_{n}^{2}}+\frac{{\|\mathbf{g}-\mathbf{A}_{\Omega_{k}}\hat{\boldsymbol{\gamma}}_{\Omega_{k}}\|}^{2}}{\sigma_{n}^{2}}+\mathcal{P}(k)\\}$
(33)

The solution of the above problem depends on the information about the noise
variance $\sigma_{n}^{2}$. If the noise variance is known, the number of
targets is obtained by solving the following problem

$\hat{k}=arg\min_{k}\\{\frac{{\|\mathbf{g}-\mathbf{A}_{\Omega_{k}}\hat{\boldsymbol{\gamma}}_{\Omega_{k}}\|}^{2}}{\sigma_{n}^{2}}+\mathcal{P}(k)\\},\hskip
22.76228pt(Known\hskip 4.26773pt\sigma^{2}_{n}).$ (34)

When the noise variance is unknown, it should be estimated from the received
data. Using ML estimation we have

$\displaystyle\hat{\sigma}_{n}^{2}$
$\displaystyle=arg\min_{\sigma_{n}^{2}}\\{-\ln{p(\mathbf{g}|\sigma_{n}^{2},k)}\\}$
$\displaystyle=arg\min_{\sigma_{n}^{2}}\\{N\ln{2\pi}+N\ln{\sigma_{n}^{2}}+\frac{{\|\mathbf{g}-\mathbf{A}_{\Omega_{k}}\hat{\boldsymbol{\gamma}}_{\Omega_{k}}\|}^{2}}{\sigma_{n}^{2}}\\}$
(35)

The minimization problem (35) can be solved analytically by taking derivative
of its argument with respect to $\sigma_{n}^{2}$ and setting it to zero. The
estimation result is as follows [37, 38]

$\hat{\sigma}_{n}^{2}=\frac{1}{N}\|{\mathbf{g}-\mathbf{A}_{\Omega_{k}}{\hat{\boldsymbol{\gamma}}_{\Omega_{k}}}}\|^{2}.$
(36)

By substituting $\hat{\sigma}_{n}^{2}$ in (33), the third term of log-
likelihood function becomes constant, equal to $N$. Hence, after omitting the
constant terms of the likelihood function (33), the number of targets is
derived as

$\hat{k}=arg\min_{k}\\{N\ln(\frac{\|{\mathbf{g}-\mathbf{A}_{\Omega_{k}}}\hat{\boldsymbol{\gamma}}_{\Omega_{k}}\|^{2}}{N})+\mathcal{P}(k)\\},\hskip
9.95863pt(Unknown\hskip 4.26773pt\sigma^{2}_{n}).$ (37)

On the other side, at the NLS estimate of $\hat{\Omega}_{k}$ and
$\hat{\boldsymbol{\gamma}}_{\Omega_{k}}$ the
$\|{\mathbf{g}-\mathbf{A}_{\Omega_{k}}}\hat{\boldsymbol{\gamma}}_{\Omega_{k}}\|^{2}$
term in (34), (37) is minimized which is equal to

$\varepsilon(k)=\min_{\Omega_{k}}\\{\mathbf{g}^{H}\Pi^{\perp}_{\Omega_{k}}\mathbf{g}\\}.$
(38)

By considering (34), (37) and (38), proof of (8) is complete.

##  Appendix B Proof of Proposition 1

First, assume $\mathbf{n}$ is a zero-mean complex circular Gaussian random
vector $\mathcal{CN}(0,\sigma^{2}_{n}\mathbf{I}_{N})$, and $\mathbf{r}$ is a
$N\times 1$ complex vector. Next, let’s define $X$ and $Y$ as follows

$\displaystyle X=\frac{1}{\sigma_{n}^{2}}\|\mathbf{n}\|^{2}$ $\displaystyle
Y=\frac{1}{\sigma_{n}^{2}}\|\mathbf{n+r}\|^{2}$ (39)

which have standard central and noncentral chi-square distributions with $N$
degrees of freedom respectively. It is desired to obtain distribution of the
difference between above random variables

$Z=\frac{1}{\sigma_{n}^{2}}(\|\mathbf{n+r}\|^{2}-\|\mathbf{n}\|^{2}).$ (40)

and also to derive its statistical parameters.

By rewriting the above equation, we have

$\displaystyle Z$
$\displaystyle=\frac{1}{\sigma_{n}^{2}}((\mathbf{n+r})^{H}(\mathbf{n+r})-\mathbf{n}^{H}\mathbf{n})$
$\displaystyle=\frac{1}{\sigma_{n}^{2}}(\mathbf{n}^{H}\mathbf{n}+\mathbf{n}^{H}\mathbf{r}+\mathbf{r}^{H}\mathbf{n}+\mathbf{r}^{H}\mathbf{r}-\mathbf{n}^{H}\mathbf{n})$
$\displaystyle=\frac{1}{\sigma_{n}^{2}}(\mathbf{n}^{H}\mathbf{r}+\mathbf{r}^{H}\mathbf{n})+\frac{\mathbf{r}^{H}\mathbf{r}}{\sigma_{n}^{2}}$
$\displaystyle=\frac{2}{\sigma_{n}^{2}}\mathrm{Re}(\mathbf{r}^{H}\mathbf{n})+\frac{\|\mathbf{r}\|^{2}}{\sigma_{n}^{2}}.$
(41)

where $\mathrm{Re(.)}$ stands for real part. According to (17) for the
expression of the noncentrality term, we have

$Z=\frac{2}{\sigma_{n}^{2}}\mathrm{Re}(\mathbf{r}^{H}\mathbf{n})+\lambda_{r}.$
(42)

By expanding the real term of $\mathbf{r}^{H}\mathbf{n}$, we have

$\displaystyle Z$
$\displaystyle=\frac{2}{\sigma_{n}^{2}}\boldsymbol{[}\mathrm{Re}(\mathbf{r^{H}})\mathrm{Re}(\mathbf{n})-\mathrm{Im}(\mathbf{r^{H}})\mathrm{Im}(\mathbf{n})\boldsymbol{]}+\lambda_{r}$
$\displaystyle=\frac{2}{\sigma_{n}^{2}}\boldsymbol{[}\mathrm{Re}(\mathbf{r})\mathrm{Re}(\mathbf{n})+\mathrm{Im}(\mathbf{r})\mathrm{Im}(\mathbf{n})\boldsymbol{]}+\lambda_{r}$
(43)

where $\mathrm{Im}(.)$ stands for imaginary part.

Each of the two terms inside the bracket are linear combinations of zero mean
Gaussian random variables $\mathrm{Re}(\mathbf{n})$ and
$\mathrm{Im}(\mathbf{n})$, hence their summation has also zero mean Gaussian
distribution. As a result, $Z$ has normal distribution with the mean equal to
the noncentrality term

$m_{Z}=\lambda_{r}$ (44)

and the expression of the variance of $Z$ is

$\displaystyle\sigma^{2}_{Z}$
$\displaystyle=E\boldsymbol{[}(Z-m_{Z})^{2}\boldsymbol{]}$
$\displaystyle=E\\{\frac{2}{\sigma_{n}^{2}}\boldsymbol{[}\mathrm{Re}(\mathbf{r})\mathrm{Re}(\mathbf{n})+\mathrm{Im}(\mathbf{r})\mathrm{Im}(\mathbf{n})\boldsymbol{]}^{2}\\}$
$\displaystyle=\frac{4}{\sigma_{n}^{4}}\\{E\boldsymbol{[}(\mathrm{Re}(\mathbf{r})\mathrm{Re}(\mathbf{n}))^{2}\boldsymbol{]}+E\boldsymbol{[}(\mathrm{Im}(\mathbf{r})\mathrm{Im}(\mathbf{n}))^{2}\boldsymbol{]}$
$\displaystyle+2E\boldsymbol{[}\mathrm{Re}(\mathbf{r})\mathrm{Re}(\mathbf{n})\mathrm{Im}(\mathbf{r})\mathrm{Im}(\mathbf{n})\boldsymbol{]}\\}.$
(45)

Considering $\mathrm{Re}(\mathbf{n})$ and $\mathrm{Im}(\mathbf{n})$ that are
two uncorrelated Gausian random variables, cross-correlation of every linear
combination of them equals to zero. Consequently the third term of the (45)
becomes zero. So

$\sigma^{2}_{Z}=\frac{4}{\sigma_{n}^{4}}\\{E\boldsymbol{[}(\mathrm{Re}(\mathbf{r})\mathrm{Re}(\mathbf{n}))^{2}\boldsymbol{]}+E\boldsymbol{[}(\mathrm{Im}(\mathbf{r})\mathrm{Im}(\mathbf{n}))^{2}\boldsymbol{]}\\}.\\\
$ (46)

To calculate the remaining terms inside the brace, it is necessary to derive
variance of a linear combination of $N$ sample zero mean Gausian random
variables. Consider $\mathbf{x}=[x_{1},x_{2},...,x_{N}]^{\mathrm{T}}$ a vector
of $N$ independent zero-mean Gausian random variables
($E(x_{i}^{2})=\sigma_{x}^{2}$, $E(x_{i}x_{j})=0,i\neq{j}$) and the weighting
vector $\boldsymbol{a}=[a_{1},a_{2},...,a_{N}]$. Then we have

$\displaystyle E{(\boldsymbol{a}\mathbf{x})^{2}}$
$\displaystyle=E({\sum^{N}_{i=1}}a^{2}_{i}x^{2}_{i}+\sum^{N}_{i=1}\sum^{N}_{j=1,j\neq{i}}a_{i}x_{i}a_{j}x_{j})$
$\displaystyle=\sum^{N}_{i=1}a^{2}_{i}E(x^{2}_{i})+\sum^{N}_{i=1}\sum^{N}_{j=1,j\neq{i}}a_{i}a_{j}E(x_{i}x_{j})$
$\displaystyle=\sigma_{x}^{2}\sum^{N}_{i=1}a_{i}^{2}$
$\displaystyle=\sigma_{x}^{2}\|\boldsymbol{a}\|^{2}.$ (47)

According to the above result, the variance of $Z$ can be derived from (47) as

$\displaystyle\sigma^{2}_{Z}$
$\displaystyle=\frac{4}{\sigma_{n}^{4}}(\sigma_{nR}^{2}\|\mathrm{Re}(\mathbf{r})\|^{2}+\sigma_{nI}^{2}\|\mathrm{Im}(\mathbf{r})\|^{2})$
$\displaystyle=\frac{4}{\sigma_{n}^{4}}(\frac{\sigma_{n}^{2}}{2}\|\mathrm{Re}(\mathbf{r})\|^{2}+\frac{\sigma_{n}^{2}}{2}\|\mathrm{Im}(\mathbf{r})\|^{2})$
$\displaystyle=\frac{2}{\sigma_{n}^{2}}\|\mathbf{r}\|^{2}$
$\displaystyle=2\lambda_{r}$ (48)

where $\sigma_{nR}^{2}$ and $\sigma_{nI}^{2}$ are variance of real and
imaginary parts of the random vector $\mathbf{n}$ which because of their
independence, are equal to half the total variance $\sigma_{n}^{2}$.

Therefore, $Z$ is a Gaussian random variable with mean equal to the non-
centrality value and variance equal to two times of that value, and the proof
is complete.

##  Appendix C Derivation of $s_{p}$, $\hat{\boldmath{\gamma_{p}}}$ and
$\lambda_{r}$

To derive the expression of $\lambda_{r}$, the parameters $s_{p}$ and
$\hat{\boldmath{\gamma_{p}}}$ should be available. Initially, it is desired to
find the location of the falsely-detected target $s_{p}$ where two targets are
present at $s_{M_{1}}$ and $s_{M_{2}}$. According to (6), $s_{p}$ is obtained
by minimizing the following cost function

$\mathbf{J}={\mathbf{g}^{H}\boldsymbol{[}\mathbf{I}_{N}-\mathbf{A}_{\Omega_{1}}(\mathbf{{A^{H}}}_{\Omega_{1}}\mathbf{A}_{\Omega_{1}})^{-1}\mathbf{{A^{H}}}_{\Omega_{1}}\boldsymbol{]}\mathbf{g}}$
(49)

As ${\mathbf{g}}^{H}\mathbf{g}$ is a constant value, the minimization in (49)
is equivalent to maximization of the following cost function

$\mathbf{J}=\mathbf{g}^{H}\mathbf{A}_{\Omega_{1}}(\mathbf{{A^{H}}}_{\Omega_{1}}\mathbf{A}_{\Omega_{1}})^{-1}\mathbf{{A^{H}}}_{\Omega_{1}}\mathbf{g}.$
(50)

As $s_{p}$ is the maximizer the above cost function,
$\mathbf{A}_{\Omega_{1}}=\boldsymbol{[}\mathbf{\textit{a}}(s_{p})\boldsymbol{]}$.
Since steering vectors are coloumns of a Fourier matrix, we have
$\boldsymbol{a}^{H}(s_{p})\boldsymbol{a}(s_{p})=N$. Hence

$\mathbf{J}={\frac{1}{N}}{\boldsymbol{|}\boldsymbol{a}^{H}(s_{p})\mathbf{g}\boldsymbol{|}}^{2}.$
(51)

Neglecting the noise term, $\mathbf{g}$ is summation of two steering vectors
at $s_{M_{1}}$ and $s_{M_{2}}$ multiplied by $\sigma_{s}$. Hence $s_{p}$ is
the optimizer of following cost function

$\mathbf{J}={\frac{\sigma_{s}^{2}}{N}}{\boldsymbol{|}\boldsymbol{a}^{H}(s_{p})(\boldsymbol{a}(s_{M_{1}})+e^{j\Delta{\phi}}\boldsymbol{a}(s_{M_{2}}))\boldsymbol{|}}^{2}.$
(52)

By introducing the correlation between $m$th and $n$th steering vectors

$\mathcal{C}_{m,n}=\frac{\boldsymbol{a}^{H}(s_{m})\boldsymbol{a}(s_{n})}{N}$
(53)

the cost function (52) can be rewritten as

$\mathbf{J}={N\sigma_{s}^{2}}{\boldsymbol{|}\mathcal{C}_{p,M_{1}}+e^{j\Delta{\phi}}\mathcal{C}_{p,M_{2}}\boldsymbol{|}}^{2}.$
(54)

According to the definition of steering vectors, the value of correlation can
be obtained as follows

$\mathcal{C}_{m,n}={\frac{1}{N}}\sum^{N}_{i=1}exp(\frac{j4\pi(s_{n}-s_{m})b_{i}}{\lambda{R_{0}}}).$
(55)

Due to equi-spaced sampling in the baseline and elevation directions, we can
write $b_{i}=(\frac{i}{N-1})\Delta_{b}$ and
$s_{n}-s_{m}=\frac{(n-m)}{M-1}\Delta_{s}$. Then, by applying some mathematical
operations (summation of a geometric progression) on (55), the correlation
term can be expressed as the following digital sinc function

$\mathcal{C}_{m,n}=\exp(-jL(n-m))\frac{\sin(NL(n-m))}{N\sin(L(n-m))}$ (56)

where the factor $L$ is equal to

$L=\frac{2{\pi}\Delta s\Delta b}{(M-1)(N-1)\lambda{R_{0}}}.$ (57)

As is seen in (56), $\mathcal{C}_{m,n}$ is a function of the difference of
indexes of two steering vectors, and can be used as
$\mathcal{C}_{mn}=\mathcal{C}(n-m)$. Also conjugate symmetry property is held,
i.e. $\mathcal{C}(n-m)=\mathcal{C}^{*}(m-n)$.

For the correlation between $M_{1}$ and $M_{2}$ indexes related to the two
targets, $\mathcal{C}_{m,n}$ in (56) can be more simplified. Due to the
definition of the normalized distance of targets
$\alpha=(s_{M_{2}}-s_{M_{1}})/\rho_{s}$ and Rayleigh resolution in (3), the
following expression is resulted

$\mathcal{C}(M_{2}-M_{1})=\exp(-jL^{\prime}\alpha)\frac{\sin(NL^{\prime}\alpha)}{N\sin(L^{\prime}\alpha)}$
(58)

where the factor $L^{\prime}=\pi/(N-1)$ and is obtained with regard to the
definition of $\rho_{s}$ and $\alpha$.

Afterwards, by substituting (58) in (54) and through some inequality math
operations, it can be shown that $s_{p}$ the position of optimizer of (54) is
the average of $s_{M1}$ and $s_{M2}$

$s_{p}=\frac{s_{M_{1}}+s_{M_{2}}}{2}.$ (59)

Equivalently, due to correspondence between elevation samples and their
indexes $p=(M_{1}+M_{2})/2$. Now it is sufficient to get the least square
estimated value $\hat{\boldmath{\gamma_{p}}}$ at location
$(s_{M_{1}}+s_{M_{2}})/2$ with regard to (7) (after scaling to $\sigma_{s}$)

$\displaystyle\hat{\boldmath{\gamma_{p}}}$
$\displaystyle=\frac{1}{\sigma_{s}}.(a^{H}(s_{p})a(s_{p}))^{-1}a^{H}(s_{p})\mathbf{g}$
$\displaystyle=\frac{1}{N}a^{H}(s_{p})\frac{\mathbf{g}}{\sigma_{s}}$
$\displaystyle=\frac{1}{N}a(\frac{s_{M_{1}}+s_{M_{2}}}{2})^{H}(a(s_{M_{1}})+e^{j\Delta{\phi}}a(s_{M_{2}}))$
$\displaystyle=\mathcal{C}(M_{1}-\frac{M_{1}+M_{2}}{2})+e^{j\Delta{\phi}}\mathcal{C}(M_{2}-\frac{M_{1}+M_{2}}{2})$
$\displaystyle=\mathcal{C}(\frac{M_{1}-M_{2}}{2})+e^{j\Delta{\phi}}\mathcal{C}(\frac{M_{2}-M_{1}}{2})$

Regarding (16) and (17), the noncentrality term can be derived as follows

$\displaystyle\lambda_{r}=$
$\displaystyle\frac{\|\mathbf{r}\|^{2}}{\sigma_{n}^{2}}$ $\displaystyle=$
$\displaystyle\frac{\sigma_{s}^{2}}{\sigma_{n}^{2}}\|\boldsymbol{a}(s_{M_{1}})+e^{j\Delta{\phi}}\boldsymbol{a}(s_{M_{2}})-\hat{\boldmath{\gamma_{p}}}\boldsymbol{a}(s_{p})\|^{2}$
$\displaystyle=$ $\displaystyle
SNR\\{\|\boldsymbol{a}(s_{M_{1}}))\|^{2}+\|\boldsymbol{a}(s_{M_{2}}))\|^{2}+|\hat{\boldmath{\gamma_{p}}}|^{2}\|\boldsymbol{a}(s_{p}))\|^{2}$
$\displaystyle+2N.\mathrm{Re}\boldsymbol{[}e^{j\Delta{\phi}}\mathcal{C}(M_{2}-M_{1})-\hat{\boldmath{\gamma_{p}}}.\mathcal{C}(p-M_{1})$
$\displaystyle-e^{-j\Delta{\phi}}\hat{\boldmath{\gamma_{p}}}.\mathcal{C}(p-M_{2})\boldsymbol{]}\\}$
$\displaystyle=$ $\displaystyle
SNR\\{2N+|\hat{\boldmath{\gamma_{p}}}|^{2}N+2N.\mathrm{Re}\boldsymbol{[}e^{j\Delta{\phi}}\mathcal{C}(M_{2}-M_{1})$
$\displaystyle-\hat{\boldmath{\gamma_{p}}}.(\mathcal{C}(\frac{M_{2}-M_{1}}{2})+e^{-j\Delta{\phi}}\mathcal{C}(\frac{M_{1}-M_{2}}{2}))\boldsymbol{]}\\}$
$\displaystyle=$ $\displaystyle
N.SNR.\\{2+|\hat{\boldmath{\gamma_{p}}}|^{2}+2.\mathrm{Re}\boldsymbol{[}e^{j\Delta{\phi}}\mathcal{C}(M_{2}-M_{1})-|\hat{\boldmath{\gamma_{p}}}|^{2}\boldsymbol{]}\\}$
$\displaystyle=$ $\displaystyle
N.SNR.\\{2+2\mathrm{Re}\boldsymbol{[}e^{j\Delta{\phi}}\mathcal{C}(M_{2}-M_{1})\boldsymbol{]}-|\hat{\boldmath{\gamma_{p}}}|^{2}\\}.$

On the other hand, regarding real and imaginary parts of the correlation term
in (58), the second and third terms inside the brace of (LABEL:62) can be
written

$\mathrm{Re}\boldsymbol{[}e^{j\Delta{\phi}}\mathcal{C}(M_{2}-M_{1})\boldsymbol{]}=cos(\Delta{\phi})\cos(L^{\prime}\alpha)\frac{\sin(NL^{\prime}\alpha)}{N\sin(L^{\prime}\alpha)}+sin(\Delta{\phi})\frac{\sin(NL^{\prime}\alpha)}{N}.$
(62) $\displaystyle|\hat{\boldmath{\gamma_{p}}}|^{2}=$ $\displaystyle
2|\mathcal{C}(\frac{M_{1}-M_{2}}{2})|^{2}+2\mathrm{Re}.\boldsymbol{[}e^{-j\Delta{\phi}}\mathcal{C}^{2}(\frac{M_{1}-M_{2}}{2}){]}$
$\displaystyle=$ $\displaystyle
2\frac{\sin^{2}(\frac{NL^{\prime}\alpha}{2})}{N^{2}\sin^{2}(\frac{L^{\prime}\alpha}{2})}+2cos(L^{\prime}\alpha-\Delta{\phi})\frac{\sin^{2}(\frac{NL^{\prime}\alpha}{2})}{N^{2}\sin^{2}(\frac{L^{\prime}\alpha}{2})}$
$\displaystyle=$ $\displaystyle
4cos^{2}(\frac{L^{\prime}\alpha}{2}-\frac{\Delta{\phi}}{2})\frac{\sin^{2}(\frac{NL^{\prime}\alpha}{2})}{N^{2}\sin^{2}(\frac{L^{\prime}\alpha}{2})}$
(63)

which results in final expression of noncentrality term as

$\displaystyle\lambda_{r}=$ $\displaystyle
N.SNR.\\{2+2cos(\Delta{\phi})\cos(L^{\prime}\alpha)\frac{\sin(NL^{\prime}\alpha)}{N\sin(L^{\prime}\alpha)}+2sin(\Delta{\phi})\frac{\sin(NL^{\prime}\alpha)}{N}$
$\displaystyle-4cos^{2}(\frac{L^{\prime}\alpha}{2}-\frac{\Delta{\phi}}{2})\frac{\sin^{2}(\frac{NL^{\prime}\alpha}{2})}{N^{2}\sin^{2}(\frac{L^{\prime}\alpha}{2})}\\}$
$\displaystyle=$ $\displaystyle N.SNR.\vartheta(\alpha,\Delta{\phi}).$ (64)

Figure 12: Plot of $\vartheta(\alpha,\Delta{\phi})$.

And the proof is done. The plot of $\vartheta(\alpha,\Delta{\phi})$ is shown
in Fig. 12.

##  Appendix D Proof of Proposition 1

Based on the assumption, detection of $k$ points is desired in the $k$
detected partial supports. Regarding (6), location of targets are estimated as

$\hat{\Omega}_{k}=\arg\min_{\Omega_{k}}\\{\mathbf{g}^{H}[\mathbf{I}-\mathbf{A}_{\Omega_{k}}(\mathbf{{A^{H}}}_{\Omega_{k}}\mathbf{A}_{\Omega_{k}})^{-1}\mathbf{{A^{H}}}_{\Omega_{k}}]\mathbf{g}\\}$
(65)

Let’s define $\Omega_{k}=[m_{1},m_{2},...,m_{k}]$ as the $k$-element support
whose elements are in a range from 1 to $M$, and $\mathbf{A}_{\Omega_{k}}$ is
the related $k$-columns steering matrix
$\mathbf{A}_{\Omega_{k}}=[\boldsymbol{a}_{m_{1}},\boldsymbol{a}_{m_{2}},...,\boldsymbol{a}_{m_{k}}]$.

Regarding $\mathbf{g}^{H}\mathbf{g}$ which is constant, the minimization
problem in (65) is converted to the following maximization problem

$\hat{\Omega}_{k}=\arg\max_{\Omega_{k}}\\{\mathbf{g}^{H}\mathbf{A}_{\Omega_{k}}(\mathbf{{A^{H}}}_{\Omega_{k}}\mathbf{A}_{\Omega_{k}})^{-1}\mathbf{{A^{H}}}_{\Omega_{k}}\mathbf{g}\\}$
(66)

Simplifying the middle term of above argument gives

$\begin{split}\mathbf{{A^{H}}}_{\Omega_{k}}\mathbf{A}_{\Omega_{k}}&=\begin{bmatrix}\boldsymbol{a}_{m_{1}}^{H}\\\
\boldsymbol{a}_{m_{2}}^{H}\\\ \vdots\\\
\boldsymbol{a}_{m_{k}}^{H}\end{bmatrix}\begin{bmatrix}\boldsymbol{a}_{m_{1}}&\boldsymbol{a}_{m_{2}}&\ldots&\boldsymbol{a}_{m_{k}}\end{bmatrix}\\\
&=\begin{bmatrix}\boldsymbol{a}_{m_{1}}^{H}\boldsymbol{a}_{m_{1}}&\boldsymbol{a}_{m_{1}}^{H}\boldsymbol{a}_{m_{2}}&\ldots&\boldsymbol{a}_{m_{1}}^{H}\boldsymbol{a}_{m_{k}}\\\
\boldsymbol{a}_{m_{2}}^{H}\boldsymbol{a}_{m_{1}}&\boldsymbol{a}_{m_{2}}^{H}\boldsymbol{a}_{m_{2}}&\ldots&\boldsymbol{a}_{m_{2}}^{H}\boldsymbol{a}_{m_{k}}\\\
\ldots&\ldots&\ldots&\ldots\\\
\boldsymbol{a}_{m_{k}}^{H}\boldsymbol{a}_{m_{1}}&\boldsymbol{a}_{m_{k}}^{H}\boldsymbol{a}_{m_{2}}&\ldots&\boldsymbol{a}_{m_{k}}^{H}\boldsymbol{a}_{m_{k}}\end{bmatrix}\end{split}$
(67)

To search $k$ points inside $k$ distinct regions, every partial support should
include one point unless some of the partial supports do not contain any point
which is in contradiction with the definition of detected support. Hence,
every partial support includes an element of $\Omega_{k}$, i.e $m_{i}\subset
supp_{i},\forall i\in[1,k]$. On the other hand, regarding the steering matrix
which is a partial Fourier matrix, the columns whose spacing is more than
$2\rho_{s}$ are orthogonal. Consequently

$\boldsymbol{a}_{m_{i}}^{H}\boldsymbol{a}_{m_{j}}=\left\\{\begin{array}[]{ll}N,\hskip
25.6073pti=j\\\ \\\ 0,\hskip 28.45274pti\neq j\end{array}\right.$ (68)

Thus

$\begin{split}\mathbf{{A^{H}}}_{\Omega_{k}}\mathbf{A}_{\Omega_{k}}=\begin{bmatrix}N&0&\ldots&0\\\
0&N&\ldots&0\\\ \ldots&\ldots&\ldots&\ldots\\\
0&0&\ldots&N\end{bmatrix}&=N\mathbf{I}_{k}\end{split}$ (69)

By substituting (69) in the maximization problem (66), the argument will be

$\begin{split}\mathbf{g}^{H}\mathbf{A}_{\Omega_{k}}(N\mathbf{I}_{k})^{-1}\mathbf{{A^{H}}}_{\Omega_{k}}\mathbf{g}&=\frac{1}{N}\mathbf{g}^{H}\begin{bmatrix}\boldsymbol{a}_{m_{1}}&\boldsymbol{a}_{m_{2}}&\ldots&\boldsymbol{a}_{m_{k}}\end{bmatrix}\begin{bmatrix}\boldsymbol{a}_{m_{1}}^{H}\\\
\boldsymbol{a}_{m_{2}}^{H}\\\ \vdots\\\
\boldsymbol{a}_{m_{k}}^{H}\end{bmatrix}\mathbf{g}\\\
&=\frac{1}{N}\mathbf{g}^{H}[(\boldsymbol{a}_{m_{1}}\boldsymbol{a}_{m_{1}}^{H}+\boldsymbol{a}_{m_{2}}\boldsymbol{a}_{m_{2}}^{H}+\ldots+\boldsymbol{a}_{m_{k}}\boldsymbol{a}_{m_{k}}^{H})]\mathbf{g}\\\
&+\frac{1}{N}\mathbf{g}^{H}[\sum\limits_{i=1}^{N}\sum\limits_{j=1,j\neq{i}}^{N}\boldsymbol{a}_{m_{i}}\boldsymbol{a}_{m_{j}}^{H}]\mathbf{g}\end{split}$
(70)

Due to orthogonality property of steering vectors mentioned in (68), the
second term of (70) becomes zero which yields the following maximization

$\hat{\Omega}_{k}=\arg\max_{\Omega_{k}}\\{\frac{1}{N}\\{|\mathbf{g}^{H}\boldsymbol{a}_{m_{1}}|^{2}+|\mathbf{g}^{H}\boldsymbol{a}_{m_{2}}|^{2}+...+|\mathbf{g}^{H}\boldsymbol{a}_{m_{k}}|^{2}\\}\\}$
(71)

When the partial supports are dis-contiguous, the terms inside the argument of
(71) become independent. Hence, maximization of total expression is equivalent
to the maximization of individual terms as

$\hat{\Omega}_{k}=\bigcup\limits^{k}_{i=1}{\arg\max_{m_{i}\in{supp_{i}}}{|\mathbf{g}^{H}\boldsymbol{a}_{m_{i}}|^{2}}}$
(72)

On the other hand, according to Algorithm 1, the detected peaks in the course
detection step are obtained by maximization of the canceled terms (residuals).
Without loss of generality and for clarity, we consider $k=2$. In this case
the two peak points are detected as follows

$p_{1}=\arg\max_{i\in{supp_{1}}}|\mathbf{g}^{H}\boldsymbol{a}(s_{i})|^{2}$
(73)
$p_{2}=\arg\max_{i\in{supp_{2}}}|\mathbf{g}^{H}_{\perp}\boldsymbol{a}(s_{i})|^{2}$
(74)

Maximization of the argument of (73) for $p_{1}$ is equivalent to maximization
of the first term of (72). Hence $m_{1}=p_{1}$. About $p_{2}$, the
maximization in (74) can be simplified as

$\displaystyle p_{2}$
$\displaystyle=arg\max_{i\in{supp_{2}}}|\mathbf{g}^{H}_{\perp}\boldsymbol{a}(s_{i})|^{2}$
$\displaystyle=arg\max_{i\in{supp_{2}}}|\boldsymbol{a}^{H}(s_{i})\mathbf{g}_{\perp}|^{2}$
$\displaystyle=arg\max_{i\in{supp_{2}}}|\boldsymbol{a}^{H}(s_{i})(\mathbf{g}-\boldsymbol{a}(s_{p_{1}})[\boldsymbol{a}(s_{p_{1}})^{H}\boldsymbol{a}(s_{p_{1}})]^{-1}\boldsymbol{a}(s_{p_{1}})^{H}\mathbf{g})|^{2}$
$\displaystyle=arg\max_{i\in{supp_{2}}}|\boldsymbol{a}^{H}(s_{i})\mathbf{g}-\boldsymbol{a}^{H}(s_{i})\boldsymbol{a}(s_{p_{1}})[\boldsymbol{a}(s_{p_{1}})^{H}\boldsymbol{a}(s_{p_{1}})]^{-1}\boldsymbol{a}(s_{p_{1}})^{H}\mathbf{g})|^{2}$

As $s_{p_{1}}$ and $s_{p_{2}}$ are far away
($|s_{p_{2}}-s_{p_{1}}|>2\rho_{s}$) due to being in discontinuous partial
supports, their steering vectors are orthogonal i.e
$\boldsymbol{a}^{H}(s_{i})\boldsymbol{a}(s_{p_{1}})=0$. Consequently, the
second term inside $|.|$ in (LABEL:76) becomes zero and $p_{2}$ is obtained as

$p_{2}=arg\max_{i\in{supp_{2}}}|\mathbf{g}^{H}\boldsymbol{a}(s_{i})|^{2}$ (76)

which is equivalent to maximizing the second term of (72). Hence,
$m_{2}=p_{2}$ and the same logic can be applied for the rest of detected peaks
when $k>2$, i.e. $m_{i}=p_{i},1\leq{i\leq{k}}$, or equivalently

$\hat{\Omega}_{k}=\\{p_{1},p_{2},...,p_{k}\\}$ (77)

## References

  * [1] A. Moreira, P. Prats-Iraola, M. Younis, et al., “A tutorial on synthetic aperture radar,” IEEE Geoscience and remote sensing magazine 1(1), 6–43 (2013).
  * [2] M. Eineder, N. Adam, R. Bamler, et al., “Spaceborne spotlight sar interferometry with terrasar-x,” IEEE Transactions on Geoscience and Remote Sensing 47(5), 1524–1535 (2009).
  * [3] A. Budillon, G. Ferraioli, and G. Schirinzi, “Localization performance of multiple scatterers in compressive sampling sar tomography: Results on cosmo-skymed data,” IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 7(7), 2902–2910 (2014).
  * [4] R. Bamler and P. Hartl, “Synthetic aperture radar interferometry,” Inverse problems 14(4), R1 (1998).
  * [5] A. Reigber and A. Moreira, “First demonstration of airborne sar tomography using multibaseline l-band data,” IEEE Transactions on Geoscience and Remote Sensing 38(5), 2142–2152 (2000).
  * [6] G. Fornaro, F. Lombardini, and F. Serafino, “Three-dimensional multipass sar focusing: Experiments with long-term spaceborne data,” IEEE Transactions on Geoscience and Remote Sensing 43(4), 702–714 (2005).
  * [7] F. Lombardini, “Differential tomography: A new framework for sar interferometry,” IEEE Transactions on Geoscience and Remote Sensing 43(1), 37–44 (2005).
  * [8] G. Fornaro, D. Reale, and F. Serafino, “Four-dimensional sar imaging for height estimation and monitoring of single and double scatterers,” IEEE Transactions on Geoscience and Remote Sensing 47(1), 224–237 (2008).
  * [9] O. Monserrat, M. Crosetto, M. Cuevas, et al., “The thermal expansion component of persistent scatterer interferometry observations,” IEEE Geoscience and Remote Sensing Letters 8(5), 864–868 (2011).
  * [10] X. X. Zhu and R. Bamler, “Let’s do the time warp: Multicomponent nonlinear motion estimation in differential sar tomography,” IEEE Geoscience and Remote Sensing Letters 8(4), 735–739 (2011).
  * [11] X. X. Zhu and R. Bamler, “Very high resolution spaceborne sar tomography in urban environment,” IEEE Transactions on Geoscience and Remote Sensing 48(12), 4296–4308 (2010).
  * [12] P. Stoica, Introduction to spectral analysis, Prentice hall (1997).
  * [13] A. De Maio, G. Fornaro, and A. Pauciullo, “Detection of single scatterers in multidimensional sar imaging,” IEEE Transactions on Geoscience and Remote Sensing 47(7), 2284–2297 (2009).
  * [14] A. Pauciullo, D. Reale, A. De Maio, et al., “Detection of double scatterers in sar tomography,” IEEE Transactions on Geoscience and Remote Sensing 50(9), 3567–3586 (2012).
  * [15] A. Budillon and G. Schirinzi, “Glrt based on support estimation for multiple scatterers detection in sar tomography,” IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 9(3), 1086–1094 (2016).
  * [16] D. L. Donoho, “Compressed sensing,” IEEE Transactions on information theory 52(4), 1289–1306 (2006).
  * [17] S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomic decomposition by basis pursuit,” SIAM review 43(1), 129–159 (2001).
  * [18] R. Tibshirani, “Regression shrinkage and selection via the lasso,” Journal of the Royal Statistical Society: Series B (Methodological) 58(1), 267–288 (1996).
  * [19] Y. C. Eldar and G. Kutyniok, Compressed sensing: theory and applications, Cambridge university press (2012).
  * [20] D. L. Donoho, M. Elad, and V. N. Temlyakov, “Stable recovery of sparse overcomplete representations in the presence of noise,” IEEE Transactions on information theory 52(1), 6–18 (2005).
  * [21] X. X. Zhu and R. Bamler, “Tomographic sar inversion by $l\\_${$1$}-norm regularization—the compressive sensing approach,” IEEE Transactions on Geoscience and Remote Sensing 48(10), 3839–3846 (2010).
  * [22] A. Budillon, A. Evangelista, and G. Schirinzi, “Three-dimensional sar focusing from multipass signals using compressive sampling,” IEEE Transactions on Geoscience and Remote Sensing 49(1), 488–499 (2011).
  * [23] X. X. Zhu and R. Bamler, “Super-resolution power and robustness of compressive sensing for spectral estimation with application to spaceborne tomographic sar,” IEEE Transactions on Geoscience and Remote Sensing 50(1), 247–258 (2012).
  * [24] S. M. Kay, “Fundamentals of statistical signal processing: Estimation theory, ptr prentice hall,” Inc., Englewood Cliffs, NJ (1993).
  * [25] E. Conte, M. Lops, and G. Ricci, “Asymptotically optimum radar detection in compound-gaussian clutter,” IEEE Transactions on Aerospace and Electronic Systems 31(2), 617–625 (1995).
  * [26] J. A. Tropp and A. C. Gilbert, “Signal recovery from random measurements via orthogonal matching pursuit,” IEEE Transactions on information theory 53(12), 4655–4666 (2007).
  * [27] H. Akaike, “A new look at the statistical model identification,” in Selected Papers of Hirotugu Akaike, 215–222, Springer (1974).
  * [28] G. Schwarz et al., “Estimating the dimension of a model,” The annals of statistics 6(2), 461–464 (1978).
  * [29] M. Wax and T. Kailath, “Detection of signals by information theoretic criteria,” IEEE Transactions on Acoustics, Speech, and Signal Processing 33(2), 387–392 (1985).
  * [30] N. Sugiura, “Further analysts of the data by akaike’s information criterion and the finite corrections,” Communications in Statistics-Theory and Methods 7(1), 13–26 (1978).
  * [31] M. Grant and S. Boyd, “Cvx: Matlab software for disciplined convex programming, version 2.1,” (2014).
  * [32] H. B. Lee, “The cramér-rao bound on frequency estimates of signals closely spaced in frequency,” IEEE Transactions on Signal Processing 40(6), 1507–1517 (1992).
  * [33] D. N. Swingler, “Frequency estimation for closely spaced sinsoids: Simple approximations to the cramér-rao lower bound,” IEEE transactions on signal processing 41(1), 489 (1993).
  * [34] M. Schmitt and U. Stilla, “Maximum-likelihood-based approach for single-pass synthetic aperture radar tomography over urban areas,” IET Radar, Sonar & Navigation 8(9), 1145–1153 (2014).
  * [35] M. A. Siddique, U. Wegmüller, I. Hajnsek, et al., “Single-look sar tomography as an add-on to psi for improved deformation analysis in urban areas,” IEEE Transactions on Geoscience and Remote Sensing 54(10), 6119–6137 (2016).
  * [36] A. Budillon, A. C. Johnsy, and G. Schirinzi, “Extension of a fast glrt algorithm to 5d sar tomography of urban areas,” Remote Sensing 9(8), 844 (2017).
  * [37] S. Kay, “Conditional model order estimation,” IEEE Transactions on Signal Processing 49(9), 1910–1917 (2001).
  * [38] P. Stoica and Y. Selen, “Model-order selection: a review of information criterion rules,” IEEE Signal Processing Magazine 21(4), 36–47 (2004).

Ahmad Naghavi received his B.Sc. degree in electrical engineering from
Amirkabir University of Technology (Tehran Polytechnic), Tehran, Iran in 2001,
and the M.Sc. degree from the Ferdowsi University of Mashhad, Iran, in 2004.
He is currently a Ph.D. student in the Electrical Engineering department of
Isfahan University of Technology (IUT). His research interests include
statistical array signal processing, compressed sensing, and TomoSAR
processing.

Mohammad Sadegh Fazel received his B.Sc. degree in electrical engineering from
Isfahan University of Technology (IUT), Isfahan, Iran in 1996, and the M.Sc.
degree from the University of Tehran, Iran, in 1998, and the Ph.D. degree from
Institute for Communication Systems, University of Surrey, UK in 2010. Since
2011, he is with the Electrical Engineering department of Isfahan University
of Technology (IUT). His current research interests include 6G wireless
communications, massive MIMO, cooperative relaying, and spectral efficiency.

Mojtaba Beheshti received the B.Sc. degree from the Isfahan University of
Technology (IUT), Iran, the M.Sc. degree from the University of Tehran, and
the Ph.D. degree from IUT in 1996, 1999, and 2011 respectively, all in
electrical engineering. Currently, he is an assistant professor with the
Information and Communication Technology Institute of IUT. His research
interests include signal processing for synthetic aperture radar, multicarrier
and Massive MIMO systems.

Ehsan Yazdian received his B.Sc. degree in electrical engineering from Isfahan
University of Technology (IUT), Isfahan, Iran in 2004, and the M.Sc. and Ph.D.
degree in electrical engineering from the Sharif University of Technology,
Tehran, Iran, in 2006 and 2012, respectively. Since 2013, he is with the
Electrical Engineering department of Isfahan University of Technology (IUT).
His research interests are in the areas of statistical array signal
processing, wireless communications, digital communication systems, and
software-defined radio.