Robust Burg Estimation of Radar Scatter Matrix for Autoregressive
structured SIRV based on Fréchet medians
Alexis Decurninge1,* and Frédéric Barbaresco1

arXiv:1601.02804v4 [math.ST] 6 May 2016

1

Thales Air Systems, Voie Pierre Gilles de Gennes 91470 Limours, France

Abstract: We address the estimation of the scatter matrix of a scale mixture of Gaussian stationary
autoregressive vectors. This is equivalent to consider the estimation of a structured scatter matrix
of a Spherically Invariant Random Vector (SIRV) whose structure comes from an autoregressive
modelization. The Toeplitz structure representative of stationary models is a particular case for the
class of structures we consider.
For Gaussian autoregressive processes, Burg method is often used in case of stationarity for its
efficiency when few samples are available. Unfortunately, if we directly apply these methods to
estimate the common scatter matrix of N vectors coming from a non-Gaussian distribution, their
efficiency will strongly decrease. We propose then to adapt these methods to scale mixtures of
autoregressive vectors by changing the energy functional minimized in the Burg algorithm.
Moreover, we study several approaches of robust modification of the introduced Burg algorithms,
based on Fréchet medians defined for the Euclidean or the Poincaré metric, in presence of outliers
or contaminating distributions. The considered structured modelization is motivated by radar applications, the performances of our methods will then be compared to the very popular Fixed Point
estimator and OS-CFAR detector through radar simulated scenarios.

1.
1.1.

Motivations
Context

Real radar measurements of strong low grazing angle clutters such as ground or sea clutters
showed that these clutters should be described by non-Gaussian distributions, especially heavytailed [44, 48, 11, 49, 50]. The family of complex spherically invariant random vectors (SIRV),
a subfamily of the elliptically symmetric distributions [34] (which contains a lot of classical distributions such as multivariate Gaussian, multivariate Cauchy distributions and multivariate Kdistributions) is a useful generalization of Gaussian random vectors, inheriting of similar shape
and location parameters. This family has been often used to modelize such radar clutters (see e.g.
[19, 23, 24, 34]). In this paper, we propose estimators of the scatter matrix of a zero-mean SIRV
with a particular structure coming from an autoregressive (AR) modelization of the correlation
between coordinates of the vector. This structured model is natural to describe the temporal correlation between radar pulse responses when signal is stationary and was already considered for
example in [14, 33, 43].
In the following, the conjugate transpose operator applied on vectors or matrices is denoted by (.)∗
while the conjugate operator applied on a scalar z is denoted by z, vectors and matrices are denoted
a.s.
by bold letters, scalars by non-bold letters and → denotes the almost sure convergence.
1

Denote x = (x1 , ..., xd )T ∈ Cd a zero-mean SIRV representing for example the response of d
pulses for a radar spatial cell. x is then characterized by the existence of a Gaussian vector y of
d
covariance Σ and a positive non-Gaussian amplitude τ such that x = τ y; Σ is called the scatter
matrix [53]. Moreover, we suppose that the scatter matrix has the same structure as the covariance
matrix of a stationary AR vector of order M (i.e. y is the trace of a Gaussian AR process of order
M ; see Section 2).
Assuming that x1 , .., xN is an independently and identically distributed (iid) sample with a SIRV
distribution on Cd , the main focus of our study lies in the estimation of the constrained scatter matrix Σ of the underlying distribution. Within this framework, we consider two kinds of robustness
for the estimation of the scatter matrix :
• (R1) a robustness with respect to the distribution of the amplitude τ which is often heavytailed and will be considered as unknown.
• (R2) a robustness with respect to a contamination in the observed sample: a fraction of the
samples are outliers or come from a different distribution (see Section 4 for radar use cases
of contamination).
1.2.

Prior art

The estimation of covariance matrix of SIRV with or without structure has been the motivation
of many works in the past, especially for radar applications. In the Gaussian framework, taking
into account the structure of a covariance matrix has been shown to improve performance of both
estimation [41] and target detection [13].
Non-Gaussian models of low grazing angle clutters involving SIRV has been early proposed in order to improve the estimation of the scatter matrix for example through non-linear transformations
or cumulants as well as the performance of target detection in this framework [23, 25, 33]. These
approaches however often consider that the law of the non-Gaussian amplitude τ is known.
Maronna proposed a class of Huber-type M-estimators of the scatter matrix Σ that do not assume
this knowledge (therefore robust in the sense of (R1)). They are defined as solution of the equation
[32] :
N
1 X
Σ̂ =
u(x∗i Σ̂−1 xi )xi x∗i .
(1.1)
N i=1
The function u has to satisfy some conditions for the estimator to be defined and consistent. A
major drawback of these estimators is their non-invariance with respect to the distribution of the
amplitude. For this sake, Tyler [45] proposed the estimator satisfying
N
d X xi x∗i
Σ̂ =
.
N i=1 x∗i Σ̂−1 xi

(1.2)

The function u(x) = x1 does not satisfy the conditions of Maronna but Tyler has shown that the
estimator solution of Eq. (1.2) is well defined and consistent. It was furthermore shown to be
a maximum likelihood estimator for normalized samples kxx11 k , .., kxxNN k (often called multivariate
signs). Some authors rediscovered and studied this estimator in its complex version in the radar
context [18, 26, 36]. These estimators however do not consider any structure on the matrix Σ.

2

In case of stationary signals, we have to take into account a Toeplitz intrinsic structure for
the scatter matrix in the SIRV framework. It can be performed by computing the constrained
maximum likelihood of normalized samples (or other statistical criterion) in the space of positive
definite matrices with Toeplitz constraints or more specific structures (see [12, 40, 42, 35, 54, 51]).
Our approach is slightly different in the sense that we propose a reparametrization of the scatter
matrix and minimize a criterion based on the induced parameters.
Indeed, the autoregressive structure allows us to split the estimation of the matrix Σ of size d × d
into d estimations of Toeplitz matrices of size 2 × 2. This splitting corresponds to the so-called
“Burg technique” [15]. Instead of estimating the covariance of the raw sample x1 , ..., xN ∈ Cd ,
we iteratively define second-order samples in C2 whose theoretical covariance can be expressed in
function of Σ (see Section 3).
This technique was originally proposed in the context of stationary Gaussian AR time series. Note
that if we consider x as the trace of an AR process of order M < d − 1, we add more structure on the matrix Σ than the Toeplitz one. Actually, given the autocovariance taps E[x1 xk ] for
k = 1, ..., M with M ≤ d − 1, it is well known that the maximum entropy model pertaining to
the vector x = (x1 , ..., xd )T in Cd results as the complex Gaussian distribution in Cd , whose covariance coincides with the autocovariance of size d × d of the trace of a Gaussian AR process of
order M (see [7, 15, 37]).
The reparametrization of any Toeplitz covariance matrix by one real positive power parameter and
d − 1 complex coefficients (called reflection parameters; see Section 3) underlying the Burg technique was also denoted by Trench (see [43]). The Toeplitz structure is then a particular case for
the class of autoregressive structures when M = d − 1.

1.3.

Outline of the paper

The Burg technique was initially expressed for the estimation of the autoregressive model of one
process (hence one range case). However, it is possible to modify Burg estimates in order to combine multiple samples coming from different range cases (“segments”) in the Gaussian context
[27]. The first contribution of the following paper is the adaptation of this Multisegment Gaussian
Burg estimation, that we will call Normalized Burg, to the SIRV context.
An alternative to Multisegment estimations is to take the mean of models coming from each estimated segment [10]. Since reflection parameters do not depend on the amplitude realization (Burg
technique separates the power parameter from reflection parameters), the robustness with respect
to (R1) is ensured. However such estimators are not robust with respect to (R2) like Multisegment
(Gaussian and Normalized) estimates. We propose then to study a geometrical method consisting
in computing the median (instead of the mean) of AR parameters estimated from x1 , ..., xN in both
Euclidean and Riemannian context (see [2, 6, 52]) and a refinement consisting in a 2-step estimation through a selection of the “better” samples presented in [3, 4]. This will allow us to cumulate
robustness with respect to heavy-tailed amplitudes (R1) and contamination (R2).
The paper is organized as follows. Section 2 presents the mixture of AR vector model. We
introduce the Normalized Burg algorithm adapted to the case of SIRV models in Section 3. We
also present the average Burg estimators as well as its robust (with respect to (R2)) modifications
in this Section. Finally, we illustrate the performances of the introduced estimators through some
simulations of radar scenarios in Section 4.
3

2.

Mixtures of autoregressive processes

Let x ∈ Cd be the random vector coming from a mixture of stationary Gaussian AR random
vectors. Then, x is characterized by the existence of a scalar random variable τ > 0 and a scatter
matrix Σ such that :
x = τy
(2.1)
where y ∼ Nd (0, Σ) is a complex Gaussian vector (called speckle) of covariance matrix Σ independent of τ (called texture). Σ is then defined up to a multiplicative constant due to the presence
of τ (we can multiply Σ and divide τ by the same positive constant without changing the vector
x). We will then consider in the following the constraint Tr(Σ) = d (see [26]).
As y comes from a stationary Gaussian AR process of order M ≤ d − 1, if we note y =
(M )
(M )
(y1 , . . . , yd )T , there exist a1 , ..., aM ∈ C such that for 1 ≤ n ≤ d :
yn +

M
X

(M )

ai

yn−i = bn

(2.2)

i=1

where bn is a complex standard Gaussian variable independent of yn−1 , ..., yn−M and with the
convention y−i = 0 for all i ≥ 0.
We can remark that x is also an AR process with non-Gaussian innovations :
xn +

M
X

(M )

ai

xn−i = τ bn .

(2.3)

i=1

3.

Burg algorithms

3.1.

Multisegment Gaussian Burg method applied to Gaussian process

We first present the well-known Burg method for Gaussian AR vectors. All the definitions we
introduce for the process y are still valid for x.
Let define the autocorrelation function for t ≥ 0 : γ(t) = E[yn+t yn ] for any n such that the
expectation has a sense. γ is independent of n because of the stationarity of y. Moreover, the
stationarity condition can be summed up by Yule-Walker equation :


γ(0)
. . . γ(M − 1)  (M ) 


aM
γ(M )


γ(1)
. . . γ(M − 2)   . 

.
(3.1)

  ..  = −  ..  .
..
..
..


.
.
.
(M )
γ(1)
a1
γ(M − 1) . . .
γ(0)
(m)

Levinson algorithm inverts this equation by introducing the successive AR parameters (ak )1≤k≤m
of order 1 ≤ m ≤ M :
• Initialization : define P0 = γ(0) and


(1)

µ1 = a1 = − γ(1)
P0
.
P1 = P0 (1 − |µ1 |2 )
4

(3.2)

• For 1 ≤ m ≤ M − 1
P

γ(m+1)+ m
(m+1)
k=1 ak (m)γ(m+1−k)
µ
=
a
=
−

m+1
m+1
Pm


2

P
=
P
(1
−
|µ
|
)

m+1
m
m
 
  (m) 
 (m) 
(m+1)
am
a1
a1
.







.
.
.

.
..


 =  ..  + µm+1  . 



(m)
(m+1)
(m)
a1
am
am

(3.3)

This algorithm enhances the role of the parameters (µm )1≤m≤M , called reflection parameters, that
are sufficient with P0 to describe the AR vector y. Note furthermore that the condition P0 = 1 is
equivalent to the aforementioned condition Tr(Σ) = d.
Instead of estimating the covariance matrix directly from the samples which does not guarantee the
Toeplitz constraint, we estimate these reflection parameters adapted for AR random vectors (we
will then use the bijection given by equations (3.2) and (3.3) to recover an estimated covariance).
For this purpose, Burg proposed in the Gaussian framework to minimize an energy at each step
1≤m≤M:
d
X
(m)
|fm (n)|2 + |bm (n)|2
(3.4)
U
=
n=m+1

with fm and bm respectively the “forward” and “backward” errors defined for m + 1 ≤ n ≤ d :
(
P
(m)
a yn−k
fm (n) = yn + m
k=1
.
(3.5)
Pm k (m)
bm (n) = yn−m + k=1 ak yn−m+k
Note that the definition of the errors is still valid for m = 0. Thanks to Equation (3.3), we can state
for m + 2 ≤ n ≤ d :

fm+1 (n) = fm (n) + µm+1 bm (n − 1)
.
bm+1 (n) = bm (n − 1) + µm+1 fm (n)
Remark 1. When there is no prior information on the model order, namely M , we should take it
as high as possible, i.e. M = d − 1. However, when N is small, this choice could lead to a poor
estimation of the reflection parameters even if the “true” model order is low. A classical way to
solve this problem is to minimize the energy U (m) + γC (m) where C (m) corresponds to a spectral
smoothness of the AR process and γ tuned the compromise between regularization and estimation;
see [5] for details on the regularization of Gaussian Burg estimators and [22] for the regularized
version of Normalized Burg defined hereafter.
The estimation of the reflection parameters consists then in the solution of the minimization of
the empirical energy for a sample x1 , .., xN :
Û

(m)

=

N
d
X
X

|fi,m (n)|2 + |bi,m (n)|2

i=1 n=m+1

where, for each i, fi,m and bi,m are the forward and backward errors for the sample yi . Knowing
µ1 , ..., µm , we have a closed-form expression for the estimate of µm+1 (this estimator is called
Multisegment Gaussian Burg see [27]) :
PN Pd
(gauss)
i=1
n=m+2 fim (n)bim (n − 1)
µ̂m+1 = arg min Û (m+1) = −2 PN P
.
(3.6)
d
2
2
µm+1
i=1
n=m+2 |fim (n)| + |bim (n − 1)|
5

3.2.

Multisegment Normalized Burg method for non-Gaussian vectors

We now consider the AR vector x. The forward and backward errors are still defined by Equation
(3.5). The estimator defined by Equation (3.6) applied for x will suffer from the disparity of the
realizations of the scalar part τ which leads us to adapt the method by considering a different
energy independent of the realizations of the texture τ :
U

(m+1)

d
X
|fm+1 (n)|2 + |bm+1 (n)|2
=
.
2 + |b (n − 1)|2
|f
(n)|
m
m
n=m+2

(3.7)

The minimum of the empirical version of the previous energy is then :
µ̂m+1 = −

N
d
X
X
bi,m (n − 1)fi,m (n)
2
.
N (d − m − 1) i=1 n=m+2 |fi,m (n)|2 + |bi,m (n − 1)|2

(3.8)

The drawback is that µ̂m+1 is not consistent. We can however correct the asymptotic bias:
Proposition 1. For 1 ≤ m ≤ M and µ̂m defined by Eq. (3.8)
a.s.

µ̂m → B1 (|µm |)

µm
|µm |

(3.9)

with B1 defined for x > 0 by :
1 − x2
B1 (x) =
x




1
log(1 − x) − log(1 + x)
+
.
(3.10)
2x
1 − x2


fm (n)
Proof. This is an application of Theorem 1 of [9] applied for the vector
. We apply
bm (n − 1)
the law of large numbers for the empirical sum µ̂m by noting that for m ≥ 0 and m + 2 ≤ n ≤ d
(Prop. 1 of [21] applied for {k1 , .., km } = {1, ..., m})

E[|fm (n)|2 ] = E[|bm (n − 1)|2 ] = Pm
.
E[fm (n)bm (n − 1)|] = −Pm µm+1

The consistent version of (3.8) is then :
−1
µ̂(u)
m = B1 (|µ̂m |)

µ̂m
.
|µ̂m |

(3.11)

B1−1 is not explicit but can be pre-computed on a grid for a gain of time (see Fig. 1).
3.3.

Average Burg estimators

3.3.1. Euclidean Mean Burg:
combine two processes:

The Gaussian and Normalized Burg estimators presented above

• the iteration on the reflection parameters (with a propagation of errors from µ1 to µM )
6

Fig. 1. Bias function B1 .
• the averaging process with the N spatial range cases
The Gaussian and Normalized Burg both perform the spatial averaging in the temporal iteration
loop.
However, since spatial range cases may not be statistically homogeneous, it seems more robust to
treat these two processes separately in order that the error due to the presence of outliers are not
propagating in the iteration on the reflection parameters. For that purpose, since the realization of
the amplitude parameter τ is shared amongst one range case, a Gaussian Burg estimation of the
reflection parameters of each range case should be performed and a spatial “average” afterwards.
Since the number of pulses d may be small, the estimation of the reflection parameters for each
(i)
range case i (denoted hereafter by µ̂m ) may not be accurate, i.e. the variance of these local estimates will be high (it is comparable to the variance of xi x∗i as an estimate of the covariance) but
the spatial averaging counters this effect as detailed below.
With previous notations, a simple average estimator (defined e.g. in [10]) is then (the superscript
.(i) will refer to the i-th range case)
µ̂m+1

P
N
N
− dn=m+2 bi,m (n − 1)fi,m (n)
1 X
1 X (i)
=
:=
µ̂
P
N i=1 dn=m+2 1/2[|fi,m (n)|2 + |bi,m (n − 1)|2 ]
N i=1 m+1

(3.12)

Note that the iterative errors bi,m and fi,m are estimated only with the temporal samples xi of the
(i)
range case i. Hence, this is not a multisegment estimator since each reflection parameter µ̂m+1 is
estimated independently from the other range cases.
The average process will not affect the bias but will divide the variance by a factor N since
"
#
N
1 X (i)
µ̂ − µ = E[µ̂(1) − µ]
E
N i=1


2
N
X
1
1
and
E
µ̂(i) − µ  =
E[|µ̂(1) − µ|2 ].
N i=1
N
Remark that the regularization evoked in Remark 1 also induces a decreasing variance but an
increasing bias (since some a priori knowledge on the spectral smoothness is introduced). For N
7

large enough, it is therefore not necessary to introduce bias since the presence of multiple spatial
samples already decreases the variance of the estimators.
3.3.2. Poincaré Mean Burg: Eq. (3.12) corresponds to the Euclidean mean of Burg estimators
of each (short) time series xi . In the space of positive definite matrices, Aubry et al [3] showed
through simulations the superiority of non-Euclidean metric in terms of performance of the target detection. In the case of reflection parameters, a generalization to an arbitrary Riemannian
geometry is also possible
!
P
− dn=m+2 bi,m (n − 1)fi,m (n)
,
(3.13)
µ̂m+1 = mean Pd
2
2
n=m+2 1/2[|fi,m (n)| + |bi,m (n − 1)| ]
where we classically call mean the so-called Fréchet mean defined as a minimizer for a certain
distance d(., .)
(1)

(2)

mean(µ , µ , ..., µ

(N )

) = arg min

|µ|<1

N
X

d(µ, µ(i) )2 .

i=1

We will consider a Riemannian metric related to the statistical model parameterized by the reflection parameters. Indeed, a natural information geometry can be associated to any parametric
model by defining a Riemannian metric through the Fisher information matrix or its dual version
(see [2, 6, 7, 8, 52]), that we will consider here, defined by
X ∂φ
ds2 =
dwi dwj
∂wi ∂wj
i,j
where w = (P0 , µ1 , ..., µM )T and φ denotes the entropy of the Gaussian AR vector y
φ(P0 , µ1 , ..., µM ) = −

M
X

(M + 1 − k) log(1 − |µk |2 ) − (M + 1) log(πe.P0 )

k=1

which gives us:
2

ds = (M + 1)



dP0
P0

2
+

M
X

(M + 1 − k)

k=1

|dµk |2
.
(1 − |µk |2 )2

Fortunately, the geometry associated to the AR model reparametrized by reflection parameters is
simple in the sense that the geodesics do not have cross products which justifies the separation of
power and reflection parameters. For each of these latter, the natural metric is the Poincaré metric
in the unit disc D = {z ∈ C s.t. |z| < 1} :
ds2 =

|dz|2
.
(1 − |z|2 )2

which leads to the following distance function


1
1+δ
µ(1) − µ(2)
(1)
(2)
dP (µ , µ ) = log
with δ =
.
2
1−δ
1 − µ(1) µ(2)
8

(3.14)

Note that we derived the Poincaré metric from the entropy of a Gaussian distribution. We ignored
the non-Gaussian amplitude since we will not use in the following the power parameter but each
reflection parameter. Moreover, we consider each reflection parameters independently since we
do not want to mix the actual reflection parameters of the clutter (of low order) and the reflection
parameters (of high order) whose value should be estimated to 0.
If one wants to take into account the amplitude distribution in a SIRV setting, the related metric
will depend on the fixed amplitude distribution (see e.g. [16]).
The mean estimators have the robustness (R1) since they are independent with respect to the
texture τ but they fail to be robust with respect to (R2). This is the reason of our interest into the
replacement of the mean by a median.
3.3.3. Robust Euclidean and Poincaré Median Burg:
Riemannian manifold is usually defined by:
(1)

(2)

median(µ , µ , ..., µ

(N )

The Fréchet median of N points in a

) = arg min

|µ|<1

N
X

d(µ, µ(i) ).

i=1

The computations of means and medians for Poincaré metric presented above are available in
[2]. In the Euclidean framework, the median can be computed through Weiszfeld algorithm (see
[46] for a modification of Weiszfeld algorithm that is convergent for all initial point), which is
initialized with a point z0 ∈ C and for t ≥ 0
PN (i)
µ /|zt − µ(i) |
zt+1 = Pi=1
.
N
(i)
i=1 1/|zt − µ |
The Poincaré metric will favor estimates close to the center by penalizing the angle inhomo(1)
(N )
geneity of the reflection parameters of the local estimates µm , ..., µm . This is illustrated by Fig.
2 in Section 4.

3.4.

2-step procedures

The idea of 2-step procedure is to use a first estimation and to select the “best” samples in order to
remove spatial samples containing potential outliers. The choice of the discarded samples is made
according to the distance of the estimated reflection parameter of each range case to the robust
estimate
)
:N )
d(µ(1:N
, µ̂m ) ≤ · · · ≤ d(µ(N
, µ̂m ).
m
m
We then keep the only N/2 closest range cases that will be supposed to be statistically homogeneous. Unlike [3] where choice of such “secondary data” is made according to the so-called
generalized inner product (GIP) xi Σ̂−1 xi , our criterion is based on the Euclidean or Riemannian
distance between reflection parameter (see Algorithm 3). Indeed, the GIP criterion depends on the
power realizations and is then not robust with respect to (R1).

9

3.5.

Algorithm summary

Algorithm 1 Generalized Burg-Levinson algorithm
Aim : Estimation of the power and reflection parameters (P0 , µ1 , ..., µM )
Input : a sample of N vectors (x1 , .., xN ) in Cd , the order of the autoregressive process M
P Pd
2
P0 = N1d N
i=1
k=1 |xik |

fi,0 (n) = xin
For 1 ≤ i ≤ N and 1 ≤ n ≤ d ,
bi,0 (n) = xin
for m = 1...M
Estimation of µ̂m from fm−1 and bm−1 (through e.g. Normalized Burg estimator)
P
(1 −
|2 )Pm−1 
m =(m)
 |µ̂m
 (m−1) 
(m−1)
a1
am−1
a1
 ..  



..
.
..
 . =

 + µ̂m 
.
(m)
(m−1)
(m−1)
am−1
a1
am−1
(m)
am = µ̂m
Forward and backward errors for 1 ≤ i ≤ N and m + 1 ≤ n ≤ d ,
fi,m (n) = fi,m−1 (n) + µ̂m bi,m−1 (n − 1)
bi,m (n) = bi,m−1 (n − 1) + µ̂m fi,m−1 (n)
end
Algorithm 2 Normalized Burg
Aim : Estimation of the m-th coefficient of reflection µ̂m+1
Input : forward and backward errors fi,m (n) and bi,m (n)
PN Pd
bi,m (n−1)fi,m (n)
2
z = − N (d−m−1)
i=1
n=m+2 |fi,m (n)|2 +|bi,m (n−1)|2


log(1−x)−log(1+x)
1
1−x2
+ 1−x2
B1 = x 7→ x
2x
z
µ̂m+1 = B1−1 (|z|) |z|

Algorithm 3 2-step Median Burg
Aim : Estimation of the reflection parameters (µ1 , ..., µM )
Input : a sample of N vectors (x1 , .., xN ) in Cd , the order of the autoregressive process M
(i)
(i)
For each range case 1 ≤ i ≤ N compute the Gaussian Burg estimates (µ̂1 , ..., µ̂M ) (Eq. 3.6)
for m=1:M
(1)
(N )
µ̂0m = median(µ̂m , ..., µ̂m )
Order the reflection parameters with respect to the distance to µ̂0m
(1:N )
(N :N )
d(µ̂m , µ̂0m ) ≤ ... ≤ d(µ̂m , µ̂0m )
(1:N )
(N/2:N )
µ̂m = median(µ̂m , ..., µ̂m
)
end

10

4.

Simulations

4.1.

Simulated scenario

As an illustration of the performances of the defined algorithms, we modelize N range cells of a
clutter burst response z1 , ..., zN through independent realization of the following random vector
d

z i = xi + w i ,

(4.1)

where
d

• x = τ y ∈ Cd is a scale mixture of AR vectors; let us recall that τ is the texture and y the
speckle (see Section 2).
• w ∈ Cd is a white noise representing the thermal noise.
We choose a Weibull texture for the model of τ (considered for example in [17]) for its adequacy
with sea and ground clutters (see also [19] for a validation of Weibull distribution on real data). A
gamma-distributed texture corresponds to a K-distributed clutter which has been often proposed in
the literature [24, 48]. However, we choose the Weibull one in order to model heavy-tailed clutters.
We recall the expression of the density for a Weibull distribution :
for x ≥ 0,

fτ (x) =

ν  x ν−1 −(x/σ)ν
e
.
σ σ

(4.2)

The scale parameter σ is taken such that E[τ ] = σΓ(1 + 1/ν) is the desired clutter power whereas
ν (taken equal to 0.6 in order to modelize strong clutters) is the shape parameter representing the
disparity of the distribution. We take N = 64 samples and a speckle built from an AR vector of
order 1 of parameter µ1 and of dimension d. An AR(1) approximates a radar ground clutter or a
wind clutter with a single Doppler frequency.
The Riemannian mean error (RME) for NM C estimations is a natural error metric in the space of
positive definite matrices thanks to its affine invariance:
RME =

1

N
MC
X

NM C

i=1

log



Σ̂i

−1/2



Σ0 Σ̂i

−1/2 

(4.3)
F

where k.kF is the Frobenius norm. We will compare the following estimators of the scatter matrix
:
• (Multisegment) Gaussian Burg : given by Eq. (3.6).
• Fixed Point (FP) : the M-estimator proposed by Tyler [45].
• (Multisegment) Normalized Burg : estimator given by Equation (3.11).
• Euclidean/Poincaré Mean Burg : estimator given by Eq. (3.12) and (3.13).
• Euclidean/Poincaré Median Burg : estimator given in Section 3.3.3.
• 2-step Euclidean/Poincaré Median Burg : Algorithm 3 by using Poincaré or Euclidean
median and distances.
11

Table 1 Riemannian Mean Error for an AR(1) (µ1 = 0.9, d = 8)

ν = 0.1
ν = 0.5
ν=1
ν=2
ν=3
ν = 10

Normalized Burg
0.42
0.42
0.42
0.42
0.42
0.42

Gaussian Burg
3.88
1.49
0.73
0.48
0.41
0.36

• 2-step Fixed Point : 2-step procedure (Section 3.4) for Fixed Point algorithm with a selection
of secondary data performed according to the likelihood of normalized samples.
The order of the above Burg estimators is taken to be maximal, i.e. M = d − 1. Since the order
of the simulated AR vector is 1, this illustrates the robustness of the approach with respect to the
choice of the order.
The estimation and detection performances will be compared to the following approaches:
• a classical OS-CFAR used together with a Hamming window applied on outputs of Doppler
Filters Bank; see for example [38].
• Ideal detection: we assume that the scatter matrix is known and use it for the test of Section
4.3. This constitutes a best-case performance benchmark for detection performances.
4.2.

Estimation quality

4.2.1. Robustness with respect to non-Gaussian amplitude: Every tested estimator at the
exception of Gaussian Burg is independent to the amplitude realizations, then we present a comparison of the estimation quality between Gaussian Burg and Normalized Burg with respect to the
Weibull shape parameter in Table 1.
Note that the Gaussian distribution corresponds to the limit case ν → ∞. This is the reason of the
good behavior of Gaussian Burg estimates in Table 1 for large ν. Moreover, as expected, the Gaussian Burg is largely outperformed by its Normalized version when the texture is sub-exponential
(ν < 1), i.e. heavy-tailed.

4.2.2. Influence of the number of pulses per range case: In Table 2, the superiority in terms
of performance of the Normalized Burg with respect to Euclidean and Poincaré Mean Burg can
be explained by the bias of the latter that is important when d is small. Indeed, the bias of Mean
Burg algorithms does not depend on N contrary to Normalized Burg for which the performances
would have been the same if we had considered a single temporal sequence of length dN . The
same conclusion for the Multisegment estimate in the Gaussian case was given in [47].
The precision of Normalized Burg is deteriorated when d increases since the order of the estimated
autoregressive model (equal to d − 1) increases with d. For large d, it is then useful to consider
Mean Burg. However, since we are interested in contaminated scenarios where d is small, we will
restrict ourselves to the case d = 12 it in the following.
4.2.3. Estimation quality illustrated through a transition scenario: We summarize in Table
3 and 4 the estimation errors for two scenarios (µ1 = 0.9 and µ1 = 0.3).
12

Table 2 Riemannian Mean Error for an AR(1) (µ1 = 0.9, ν has no impact)

d=8
d = 16
d = 32
d = 64

Normalized Burg
0.42
0.44
0.47
0.50

Euclidean Mean Burg Poincaré Mean Burg
0.77
0.76
0.63
0.63
0.55
0.55
0.49
0.49

Table 3 Riemannian Mean Error for an AR(1) (µ1 = 0.9, d = 12); these errors are independent of the shape
parameter of the Weibull texture. We progressively increase the number of contaminating range cells (outliers).

Euclidean
Normalized Burg
Median Burg
0 outlier
5 outliers
10 outliers
20 outliers
30 outliers

0.42
1.08
2.03
3.17
3.77

0.57
0.64
0.77
1.10
1.96

2-step
Euclidean Median Burg
1.01
0.98
0.92
0.87
0.87

Poincaré
Median Burg
0.78
0.90
1.04
1.41
2.56

FP
0.70
0.78
1.30
2.70
3.71

Table 4 Riemannian Mean Error for an AR(1) (µ1 = 0.3, d = 12); these errors are independent of the shape
parameter of the Weibull texture. We progressively increase the number of contaminating range cells (outliers).

Euclidean
Normalized Burg
Median Burg
0 outlier
5 outliers
10 outliers
20 outliers
30 outliers

0.34
0.35
0.41
0.59
0.71

0.38
0.38
0.41
0.51
0.67

2-step
Euclidean Median Burg
0.90
0.85
0.85
0.85
0.97

Poincaré
Median Burg
0.30
0.33
0.37
0.49
0.64

FP
0.70
0.70
0.70
0.76
0.87

• In the non contaminated case, Normalized Burg show better accuracy than FP: taking into
account the Toeplitz structure of the scatter matrix then improves the estimation quality.
• When the spectrum is flatter (|µ1 | = 0.3), the Poincaré metric that favors the small coefficients
is slightly more efficient than the Euclidean one. Indeed, the reflection coefficients of high
order are then closer to 0. For reflection parameters of high modulus however, this behavior
makes Poincaré metric less efficient especially in the contaminated cases (see Fig. 2).
• The 2-step procedure drastically increases the estimation quality in the case of a strong contamination and a strong correlation (µ1 = 0.9). In that case, the outliers coming from the
“true” distribution are well separated from the correlation samples coming from the perturbing distribution and the 2-step procedure can then easily separate the two parts of the sample.
Otherwise, when |µ1 | is low, this separation is less clear.
• The surprising decrease of the error in the 2-step procedure when the amount of contamination
increases may be explained by the fact that the higher order reflection coefficient are better
and better estimated thanks to the diversity brought by the contamination.
13

Estimate of each range case
Euclidean median
Poincare median

1

0.5

0.5

0

0

−0.5

−0.5

−1

−1

−0.5

0

0.5

Estimate of each range case
Euclidean median
Poincare median

1

−1

1

(a) contamination by 30 range cases, µ1 = 0.9

−1

−0.5

0

0.5

1

(b) no contamination, µ1 = 0.9

Fig. 2. Estimated first coefficient of reflection for each range and their Riemannian and Euclidean
medians in case of a contamination by 30 range cases.

(a) Simulated spectra

(b) Normalized Burg

(c) Fixed Point

(d) OS-CFAR

(e) Poincaré Median Burg

(f) Euclidean Median Burg

(g) 2-step Poincaré Median Burg

(h) 2-step Euclidean Median Burg

(i) 2-step Fixed Point

Fig. 3. Estimated and simulated spectra for 100 range cases (The x-axis corresponds to the range
dimension whereas y-axis indicates the normalized frequency).

14

In Fig. 3, we illustrate the robustness of each estimator through a clutter transition. We considered a scenario where range cases 1 to 50 are simulated through an AR(1) of parameter µ1 = 0.9
and range cases 51 to 100 are simulated with an AR(1) of parameter µ1 = 0.9e0.3×2iπ (Fig. 3(a)).
For each test range case xi (i is represented in the x-axis), the represented Doppler spectrum results from the estimated covariance of the N = 64 neighbor cells xi−32 , ..., xi−1 , xi+1 , ..., xi+32 . In
order to control edge effects, we consider only available neighbor cells for 32 first and 32 last cells.
• For the non robust estimators (namely Normalized Burg and Fixed-Point) the estimated spectra have two frequencies for cases around the transition which is not the case for the other
estimators.
• For the robust estimators, the number of range cases where the estimated spectrum is not
accurate is respectively for OS-CFAR, Poincaré Median Burg, Euclidean Median Burg and
their respective 2 step versions of 16, 8, 8 , 8 and 2 cases. With this property, the detection
of a target with a normalized frequency of 0.3 is possible in an area close to the transition for
robust estimators.
• The 2-step procedure alone is not sufficient if the first estimation is not robust: this is illustrated by Fig. 3(i) where the secondary data selection after a first Fixed-Point estimation is
not sufficient to separate the two clutter frequencies.
• It can be observed that the spectra of OS-CFAR show a frequency resolution worse than its
competitors. This is due to the low number of pulses (d = 12) for each range cells which is
responsible for the low number of filters in OS-CFAR.
4.2.4. Estimation quality for a sea clutter scenario: In Fig. 4-5, we simulate a scenario
encountered when we face sea clutter, namely the position of the “peak” in the spectra (respectively
the spectral width for Fig. 5) of neighbor range cases is not stable and can be drifting (Fig. 4(a)).
Ideally, the estimated “mean” spectra should correspond to a mean behavior, spectrally speaking:
the position of the peak should be the mean of the neighbor peaks as well as the spectral width.
• In Fig. 4, since Fixed Point and Normalized Burg estimators take into account every neighbor
case with the same weights, the estimated spectrum is wider, this width representing the
incertitude on the position of the peak. On the other-hand, the median-based estimators are
only dependent on the considered Riemannian geometry in the space of reflection parameters.
Indeed, with our choice of geometry, the more diversity in the parameters, the lower the
absolute value of the median of these parameters and then the larger the spectrum. Since
the 2-step Euclidean Median estimator takes into account less neighbor cases, the diversity is
weaker and then, the accuracy of the estimated spectrum is higher. On the other hand, 2-step
Poincaré Median estimator is sensitive to highest order reflection parameters that are more
biased if we consider less range cases; this effect is moreover not sufficiently compensated
by the first reflection parameters.
• Similarly, in Fig. 5, robust estimators estimate more accurately the width of the “average”
spectrum while non-robust estimators under-estimate it.
4.3.

Detection quality

We will now compare the estimators through their detection performances (for the sake of clarity,
we restrain ourselves to Normalized Burg, 2-step Burg estimators, Fixed-Point and OS-CFAR). For
15

(a) Simulated spectra of neighbor range cases

Normalized Power (dB)

−10

Euclidean Median
2-step Poincaré Median
2-step Euclidean Median
Poincare Median Burg
Ideal

0
Normalized Power (dB)

Gaussian Burg
Normalized Burg
FP
Ideal

0

−20

−10

−20

−0.4

−0.2

0

0.2

−0.4

0.4

Normalized Frequency

−0.2

0

0.2

0.4

Normalized Frequency

(b) Estimated “mean” spectra of non-robust estimators

(c) Estimated “mean” spectra of robust estimators

Fig. 4. Simulated and estimated spectrum in a sea clutter typical scenario in a Gaussian context;
the simulated autoregressive process is an AR(1) with µ1 = 0.9.
that purpose, we assume that a cell under test is spatially surrounded by N neighbor cells sharing
the same distribution or not (depending on the scenario). We suppose that a target is present in the
cell under test. As for the target model, we consider small targets in the sense that they are present
in only one range cell which burst response is simulated by
d

z = αp + x + w,

(4.4)

where α represents the target power, p = (1, e2iπfD , ..., e2iπ(d−1)fD )T , fD the normalized Doppler
frequency and x and w are defined as in Eq. (4.1).
Several test statistics dedicated to target detection in a non-Gaussian environment have been proposed in the radar literature. A detector classically used is the GLRT (Generalized Likelihood
Ratio Test), also called ANMF (see e.g. [24]). A variant of the GLRT detector has been proposed
in [1] in the context of targets spread over several range cells present in an homogeneous or inhomogeneous clutter modelized as an autoregressive process.
We will use in the sequel the GLRT statistics and propose to compare its performances with a new
geometrical detector based on the geometry of reflection parameters.

16

(a) Simulated spectra of neighbor range cases

Gaussian Burg
Normalized Burg
FP
Ideal

−10

−20

−30

−0.4

−0.2

0

0.2

Euclidean Median
2-step Poincaré Median
2-step Euclidean Median
Poincare Median Burg
Ideal

0

Power (dB)

Power (dB)

0

−10

−20

−30

0.4

Normalized Frequency

−0.4

−0.2

0

0.2

0.4

Normalized Frequency

(b) Estimated “mean” spectra of non-robust estimators

(c) Estimated “mean” spectra of robust estimators

Fig. 5. Simulated and estimated spectrum in a sea clutter typical scenario in a Gaussian context;
the simulated autoregressive process is an AR(3).
4.3.1. GLRT detector: Denoting by Σ̂ one estimator of the scatter matrix of the neighbor cells,
the GLRT detector is defined by
GLRT(z) =

max
θ∈[−0.5;0.5[

|p(θ)∗ Σ̂−1 z|2
(z∗ Σ̂−1 z)(p(θ)∗ Σ̂−1 p(θ))

(4.5)

with p(θ) = (1, e2iπθ , ..., e2iπ(d−1)θ )T the steering vector and z ∈ Cd the data of the cell under test.
We compute the test threshold such that the probability of false alarm is set to 10−3 and compare
the probabilities of detection with the classical OS-CFAR test [38].
Figures 6-7 take into account a scenario where N −Nout neighbor cells are simulated through an
AR(1) of parameter µ1 (considered as clutter) and Nout contaminating cells are simulated through
AR(1) of parameter µ1 e0.3×2iπ of same power (considered as contaminating cells). Moreover, we
insert a target at different frequencies represented in the x-axis in the cell under test.
Figures 6-7 illustrate the Doppler resolution of the detectors with respect to the clutter and the
contaminating cells. A high frequency resolution is crucial in order to be able to detect small
targets (with moderate power) with velocity similar to the (ambient or contaminating) clutter’s.
17

Normalized Burg
2-step Euclidean Median
2-step Poincare Median
FP
Ideal
OS-CFAR

0.8

1
Probability of detection

Probability of detection

1

0.6
0.4
0.2
0

−0.4

−0.2

0

0.2

Normalized Burg
2-step Euclidean Median
2-step Poincare Median
FP
Ideal
OS-CFAR

0.8
0.6
0.4
0.2
0

0.4

−0.4

−0.2

Frequency

(a) No outliers (µ1 = 0.7, SCR=0dB)

Normalized Burg
2-step Euclidean Median
2-step Poincare Median
FP
Ideal
OS-CFAR

0.8

0.4
0.2

−0.4

−0.2

0

0.4

0.2

Normalized Burg
2-step Euclidean Median
2-step Poincare Median
FP
Ideal
GLRT detector
OS-CFAR

1

0.6

0

0.2

(b) 10 contaminating range cells with normalized
frequency 0.3 (µ1 = 0.7, SCR=0dB)

Probability of detection

Probability of detection

1

0
Frequency

0.8
0.6
0.4
0.2
0

0.4

Frequency

−20

0

20

SCR (dB)

(c) 10 contaminating range cells with normalized (d) Detection probability vs. Signal Clutter Ratio
frequency 0.3 (µ1 = 0.9, SCR=0dB)
(SCR) with a target of frequency 0.3 (µ1 = 0.7, no
outliers)

Fig. 6. Detection probability of a target vs. normalized frequency of the target or Signal-Clutter
Ratio (SCR) with a clutter-to-noise ratio CNR=40dB for a fixed PFA = 10−3 (GLRT detector).
We observe in Figures 6-7 that, since the Signal-Clutter Ratio is 0dB, the probability of detection
for frequencies close to the peak frequency of the clutter (i.e. 0) is close to 0. For normalized
frequencies close to 0.3, only the performances of FP and Normalized Burg estimators as well as
OS-CFAR detector decrease which illustrate the robustness of 2-step and Median procedures with
respect to outliers. With 10 contaminating range cells, the gap of detection probability can go up
to 30%.
We illustrate in Figure 6(d) and 7(d)the detection probability with respect to the Signal-Clutter
Ratio for a fixed target velocity. These two figures shows only slight differences between studied
estimators (at the exception of OS-CFAR) in the non-contaminated case since the estimators are
all independent with respect to the non-Gaussian texture.
Moreover, Burg estimators have a better frequency resolution with respect to the classical OS18

CFAR estimation for the same reason explained in Section 4.2.3.

4.4.

Comparison to a new geometrical detector

Instead of using the GLRT statistics, we consider the geometrical detector as alternative for Normalized Burg and 2-step Burg estimators :
AR(z) =

M
X

(M − k + 1)dP (µ̂k (z), µ̂k,amb )2

(4.6)

k=1

where (µ̂1,amb , ..., µ̂M,amb ) are the estimated reflection coefficients of the N surrounding cells and
(µ̂1 (z),..., µ̂M (z)) is a regularized estimation of the underlying autoregressive process of the cell
under test z ∈ Cd (see for example [5]). Let us recall that dP is the Riemannian distance in the
Poincaré disc considered in Section 3.3.2.
This detector does not directly provide an estimation of the normalized frequency of the detected
target. However, its performances, illustrated by Figure 7, show competitive results with respect
to the same scenario for the GLRT detector with a reduced complexity (neither an inversion of a
d × d matrix nor a max computation are needed for this detector).
4.5.

CFAR property of detectors

We illustrate the CFAR property of different detectors in Table 5 by computing the true probability
of false alarm (PFA) of the test statistics for a unique threshold. This property is crucial for the
detectors to be stable on real operations.
Since all proposed estimates are independent with respect to texture realizations (neglecting thermal noise), the CFAR property with respect to the texture distribution holds for both GLRT and
AR detectors.
Moreover, the GLRT detector coupled with Fixed Point estimator was shown to have CFAR property with respect to the clutter covariance (see e.g. [20] which also contains validation of this
property on real radar measurements); it is confirmed by Table 5 where only thermal noise have
a small impact on the PFA in the non-contaminated case. The presence of contaminating range
cells, however, breaks the CFAR property, especially for non-robust estimators (Normalized Burg
and FP).
Table shows that GLRT detector is more robust than AR detector with respect to CFAR property.
Let us note furthermore that the GLRT detector is more robust than AR detector with respect to
CFAR property. This could lead to further improvements of the latter.
5.

Conclusion

We have presented Burg methods for mixtures of autoregressive vectors which are a natural family
of distributions when we consider non-Gaussian stationary clutter. We proposed several estimators
independent of the non-Gaussian texture in this framework and studied their behavior especially
in terms of robustness with respect to contamination and efficiency. We can sum up these through
the following insights:
We have presented several Burg methods for mixtures of autoregressive vectors which are a
natural family of distributions when we consider non-Gaussian stationary clutter. We proposed
19

2-step Euclidean Median
2-step Poincare Median
Normalized Burg
OS-CFAR

1
Probability of detection

Probability of detection

1
0.8
0.6
0.4
0.2
0

−0.4

−0.2

0

2-step Euclidean Median
2-step Poincare Median
Normalized Burg
OS-CFAR

0.2

0.8
0.6
0.4
0.2
0

0.4

−0.4

−0.2

Frequency

1
Probability of detection

Probability of detection

2-step Euclidean Median
2-step Poincare Median
Normalized Burg
AR detector
OS-CFAR

0.6
0.4
0.2

−0.4

−0.2

0

0.4

(b) 10 contaminating range cells with normalized
frequency 0.3 (µ1 = 0.7, SCR=0dB)

2-step Euclidean Median
2-step Poincare Median
Normalized Burg
OS-CFAR

0.8

0

0.2

Frequency

(a) No outliers (µ1 = 0.7, SCR=0dB)

1

0

0.2

0.8
0.6
0.4
0.2
0

0.4

Frequency

−20

0

20

SCR (dB)

(c) 10 contaminating range cells with normalized (d) Detection probability vs. Signal Clutter Ratio
frequency 0.3 (µ1 = 0.9, SCR=0dB)
(SCR) with a target of frequency 0.3 (µ1 = 0.7, no
outliers)

Fig. 7. Detection probability of a target vs. normalized frequency of the target or Signal-Clutter
Ration (SCR) with a clutter-to-noise ratio CNR=40dB for a fixed PFA = 10−3 (AR detector).
several estimators independent of the non-Gaussian texture in this framework and studied their
behavior especially in terms of robustness with respect to contamination and efficiency. We can
sum up these through the following insights:
• Thanks to the exploitation of the Toeplitz structure, Burg methods complexity goes from
O(d3 ) (that corresponds to the complexity of the inversion of a matrix d × d) to O(M2 )
(recall that M stands for the order of the autoregressive model).
• It is useful to take into account the non-Gaussianity of the clutter for sub-exponential amplitude distributions.
• Considering medians instead of means in the Burg estimators furnishes a robustness for
medium contamination (10% to 30% outlier samples).
• 2-step procedures (consisting in a selection of secondary data) need a first estimation of the
20

Table 5 False alarm probability (to be multiplied by 10−3 ) for a fixed threshold computed
through Monte-Carlo simulations for the scenario with µ1 = 0.7 without contamination; x/y means PFA=x.10−3 for
the GLRT detector and PFA = y.10−3 for the AR detector (unavailable for FP estimator).

no contamination
(µ1 = 0.3)
no contamination
(µ1 = 0.7)
no contamination
(µ1 = 0.9)
20 contaminating range
cells (µ1 = 0.7)
20 contaminating range
cells (µ1 = 0.9)

Normalized Burg

2-step Euclidean
Median Burg

2-step Poincaré
Median Burg

FP

0.8/1.2

2.5/1.2

0.6/1.1

0.6/NA

1.0/1.0

1.0/1.0

1.0/1.0

1.0/NA

2.8/8.8

3.0/12.0

2.6/5.0

1.6/NA

5.8/15.2

3.2/0.40

2.2/0.60

3.2/NA

15.0/164.8

3.8/4.0

4.0/5.6

5.8/NA

scatter matrix of the clutter that is robust enough in order to be efficient. In that case,
high contamination (close to 50%) can be considered.
• Taking into account the Poincaré metric is efficient for reflection parameters of high order
(that should be close to 0) since they tends to under-estimate their modulus.
Future works will be devoted to the extension of these methods for non-stationary signal in
the burst [31, 39]. This work will be also extended for STAP (Space-Time Adaptive Processing)
based on OS-STAP algorithm described in [6] by computing mean and median Toeplitz-BlockToeplitz covariance matrices parameterized by Matrix-Valued Autoregressive model, and based on
numerical scheme described in [28, 29, 30].
Acknowledgment
The authors would like to thank the DGA/MRIS for its support as well as the referee for useful
comments.
6.

References

[1] Alfano, G., de Maio, A., Farina, A.: Model-based adaptive detection of range-spread targets,
IEE Proc.-Radar Sonar Navig, vol. 151, no. 1, pp.2-10, 2004.
[2] Arnaudon, M., Barbaresco, F., Yang, L.: Riemannian Medians and Means With Applications
to Radar Signal Processing, IEEE Journal of Selected Topics in Signal Processing, vol. 7, n.
4, August 2013.
[3] Aubry, A., De Maio, A., Pallotta, L., Farina, A.: Covariance matrix estimation via geometric
barycenters and its application to radar data selection, IET Radar, Sonar and Navigation, ,
vol. 7, no. 6, p. 600 - 614, 2013
[4] Aubry, A., De Maio, A., Pallotta, L., Farina, A.: Median matrices and their application to
radar training data selection, IET Radar, Sonar & Navigation, vol.8, no.4, p. 265 - 274, 2014
21

[5] Barbaresco, F.: Super resolution spectrum analysis regularization: Burg, capon and ago antagonistic algorithms in Proc. EUSIPCO’96, Trieste, pp. 2005-2008, Sept. 1996.
[6] Barbaresco, F.: Information Geometry of Covariance Matrix : Cartan-Siegel Homogeneous
Bounded Domains, Mostow/Berger fibration and Frechet median, in R. Bhatia & F.Nielsen
Ed., “Matrix Information Geometry”, Springer Lecture Notes in Mathematics, 2012.
[7] Barbaresco, F.: Koszul Information Geometry and Souriau Geometric Temperature/Capacity
of Lie Group Thermodynamics, Entropy, 16, pp 4521-4565., published in MDPI Book “Information, Entropy and Their Geometric Structures”, 2014.
[8] Barbaresco, F.: Symplectic Structure of Information Geometry: Fisher Metric and EulerPoincaré Equation of Souriau Lie Group Thermodynamics, Geometric Science of Information,
Vol. 9389 of the Springer series Lecture Notes in Computer Science, pp 529-540, 2015.
[9] Bausson, S., Pascal, F., Forster, P., Ovarlez, J.P., Larzabal, P.: First- and Second-Order Moments of the Normalized Sample Covariance Matrix of Spherically Invariant Random Vectors,
IEEE Signal Processing Letters, vol. 14, no. 6, 2007.
[10] Beex, A.A., Rahman, M.D.A.: On Averaging Burg Spectral, Estimators for Segments, IEEE
Trans. on Acoustics, speechand Signal Proc., vol. 34, no. 6,1986.
[11] Billingsley, J. B.: Ground clutter measurements for surface-sited radar, MIT, Tech. Rep. 780,
1993.
[12] Bini, D., Iannazzo, B., Jeuris, B., Vandebril, R.: Geometric means of structured matrices,
BIT, vol. 54, no. 1, pp. 55-83, 2014.
[13] Bose, S., Steinhardt, A.O.: A Maximal Invariant Framework for Adaptive Detection with
Structured and Unstructured Covariance Matrices, IEEE Trans. Signal Proc., vol. 43, no. 9,
1995
[14] Bouvier, C., Martinet, L., Favier, G., Artaud, M.: Simulation of radar sea clutter using autoregressive modelling and K-distribution IEEE International Radar Conference, pp. 425 430, 1995
[15] Burg, J. P.: Maximum entropy spectral analysis, Proceedings of the 37th Meeting of the
Society of Exploration Geophysicists, 1967, Reprinted in Modern Spectrum Analysis, D.G.
Childers, ed., IEEE Press, New York, 1978, pp. 34-41.
[16] Calvo, M., Oller, J.: A distance between elliptical distributions based in an embedding into
the Siegel group, Journal of Computational and Applied Mathematics, 145, pp. 319-334, 2002
[17] Chong, C.Y., Pascal, F., Ovarlez, J-P., Lesturgie, M.: MIMO Radar Detection in NonGaussian and Heterogeneous Clutter, IEEE Journal of selected topics in signal processing,
vol. 4, no. 1, pp 115-126, 2010.
[18] Conte, E., de Maio, A., Ricci, G.: Recursive Estimation of the Covariance Matrix of a
Compound-Gaussian Process and Its Application to Adaptive CFAR Detection, IEEE Trans.
Signal Processing, vol. 50, no. 8, 2002.
[19] Conte, E., de Maio, A., Galdi, C.: Statistical Analysis of Real Clutter at Different Range
Resolutions, IEEE Trans. Aerospace and Electronic Systems, vol. 40, no. 3, pp. 903-918, 2004.
22

[20] Conte, E., de Maio, A.: Mitigation Techniques for Non-Gaussian Sea Clutter, IEEE Journal
of Oceanic Engineering, vol. 29, no.2, pp. 284-302, 2004.
[21] Brockwell, P.J., Dalhaus, R.: Generalized Durbin-Levinson and Burg Algorithms, Journal of
Econometrics, 118, pp 129-149, 2003.
[22] Decurninge, A.: Univariate and multivariate quantiles, probabilistic and statistical approaches; radar applications PhD Thesis, Pierre and Marie Curie University, Paris, 2014.
[23] Farina, A., Russo, A., Scannapieco, F., Barbarossa, S.: Theory of radar detection in coherent
Weibull clutter, IEE Proc. Communications, Radar and Signal Processing, vol.134 , no.2, pp.
174-190, 1987
[24] Gini, F.: Sub-optimum coherent radar detection in a mixture of K-distributed and Gaussian
clutter, IEE Proc. Radar Sonar Navig., vol. 144, no. 1, pp. 39-48, 1997.
[25] Gini, F.: Cumulant-Based Adaptive Technique for Coherent Radar Detection in a Mixture of
K-Distributed Clutter and Gaussian Disturbance, IEEE Trans. Signal Processing, vol. 45, no.
6, pp. 1507-1519, 1997
[26] Gini, F., Greco, M.: Covariance matrix estimation for CFAR detection in correlated heavy
tailed clutter, Signal Processing, vol.82, no.12, pp. 1847-1859, 2002.
[27] Haykin, S., Currie, B.W., Kesler, S.B.: Maximum-entropy spectral analysis of radar clutter,
Proc. IEEE, vol. 70, 1982.
[28] Jeuris, B., Vandebril, R.: Averaging block-Toeplitz matrices with preservation of Toeplitz
block structure, SIAM Conference on Applied Linear Algebra (ALA), 2015.
[29] Jeuris, B., Vandebril, R.: The Kähler mean of Block-Toeplitz matrices with Toeplitz structured
block, TW Reports, TW660, 20, Department of Computer Science, KU Leuven, 2015
[30] Jeuris, B.: Riemannian Optimization for Averaging Positive Definite Matrices, PhD Thesis,
University of Leuven, 2015.
[31] Le Brigant, A., Barbaresco, F., Arnaudon, M.: Reparameterization Invariant Metric on the
Space of Curves, Geometric Science of Information, Vol. 9389 of the Springer series Lecture
Notes in Computer Science, pp 140-149, October 2015
[32] Maronna, R.A.: Robust M-estimators of multivariate location and scatter, The Annals of
Statistics, vol. 4, No 1, pp 51-67, 1976.
[33] Rangaswamy, M., Michels, J.H., Weiner, D.D.: Multichannel Detection for Correlated NonGaussian Random Processes Based on Innovations, IEEE Trans. Signal Processing, vol. 43,
no. 8, pp. 1915-1922, 1995
[34] Ollila, E., Tyler, D., Koivunen, V., Poor, V.: Complex Elliptically Symmetric Distributions :
Survey, New Results and Applications, IEEE Transactions on signal processing, vol. 60, no.
11, 2012.
[35] Pailloux, G.: Estimation Structurée de la Covariance du Bruit en Détection Adaptative PhD
Thesis, University of Paris Ouest Nanterre La Défense, 2010.

23

[36] Pascal, F., Forster, P., Ovarlez, J.P., Larzabal, P.: Performance Analysis of Covariance Matrix
Estimates in Impulsive Noise, IEEE Trans. on signal processing, vol. 56, no. 6, pp. 2206-2217,
2008.
[37] Pavon, M., Ferrante, A.: On the Geometry of Maximum Entropy Problems, SIAM Review,
Vol. 55, No. 3, 2013.
[38] Rohling, H.: Radar CFAR thresholding in clutter and multiple target situations, IEEE Transactions on Aerospace and Electronic Systems, no. 19, pp 608-621, 1983.
[39] Ruiz, M., Barbaresco, F.: Radar detection for non-stationary Doppler signal in one burst
based on Information Geometry distance between paths, International Radar Symposium,
IRS’15, Germany, 2015.
[40] Soloveychik, I., Wiesel, A.: Tyler’s Covariance Matrix Estimator in Elliptical Models With
Convex Structure, IEEE Trans. Signal Processing, vol. 62, no. 20, pp. 5251-5259 , 2014
[41] Li, H., Stoica, P., Li, J.: Computationally Efficient Maximum Likelihood Estimation of Structured Covariance Matrices, IEEE Trans. Signal Processing, vol. 47, no. 5, 1999
[42] Sun, Y., Babu, P., Palomar, D. P.: Robust Estimation of Structured Covariance Matrix for
Heavy-Tailed Elliptical Distributions, http://arxiv.org/abs/1506.05215
[43] Trench, W.: An algorithm for the inversion of finite Toeplitz matrices, J. Soc. Indust. Appli.
Math., Vol. 12, no. 3, pp 512-522, 1964.
[44] Trizna, D. B.: Statistics of low grazing angle radar sea scatter for moderate and fully developed ocean waves, IEEE Trans. Antennas Propag., vol. 39, no. 12, pp. 1681-1690, 1991
[45] Tyler, D.: A Distribution-Free M-Estimator of Multivariate Scatter, The Annals of Statistics,
Vol. 15, no. 1, pp 234-251, 1987.
[46] Vardi, Y., Zhang, C.: The multivariate L1-median and associated data depth, Proceedings
of the National Academy of Sciences of the United States of America, col. 97, no. 4, pp.
1423-1426 (electronic), 2000.
[47] de Waele, S., Broersen, P.: The Burg Algorithm for Segments, IEEE Trans. on Signal Proc.,
vol. 48, no. 10, pp. 2876-2880, 2000.
[48] Watts, S.: Radar detection prediction in sea clutter using the compound K-distribution model,
IEE Proc. F, vol. 132, no. 7, pp. 613-620, 1985.
[49] Watts, S., Rosenberg, L., Bocquet, S., Ritchie, M.: Doppler spectra of medium grazing angle
sea clutter; part 1: characterisation, in Radar, Sonar & Navigation, IET , vol.10, no.1, pp.2431, 1 2016.
[50] Watts, S., Rosenberg, L., Bocquet, S., Ritchie, M.: Doppler spectra of medium grazing angle
sea clutter; part 2: model assessment and simulation, in Radar, Sonar & Navigation, IET ,
vol.10, no.1, pp.32-42, 1 2016.
[51] Wiesel, A., Zhang, T.: Structured Robust Covariance Estimation, Foundations and Trends(r)
in Signal Processing, Vol. 8, No. 3, pp 127-216, 2015.

24

[52] Yang, L.: Medians of probability measures in Riemannian manifolds and applications to
radar target detection, PhD Thesis, University of Poitiers, 2012.
[53] Yao, K.: A representation theorem and its applications to spherically invariant random processes, IEEE Trans. Inf. Theory, vol. IT-19, no. 5, pp. 600-608, 1973.
[54] Zhang, T., Wiesel, A., Greco, M.: Multivariate generalized gaussian distribution: Convexity
and graphical models, IEEE Trans. Signal Process., vol. 61, no. 16, pp. 4141-4148, 2013.

25

