Matrix CFAR detectors based on symmetrized Kullback–Leibler and total Kullback–Leibler divergences

Digital Signal Processing 69 (2017) 106–116

Contents lists available at ScienceDirect

Digital Signal Processing www.elsevier.com/locate/dsp

Matrix CFAR detectors based on symmetrized Kullback–Leibler and total Kullback–Leibler divergences Xiaoqiang Hua ∗ , Yongqiang Cheng, Hongqiang Wang, Yuliang Qin, Yubo Li, Wenpeng Zhang a r t i c l e

i n f o

Article history: Available online 30 June 2017 Keywords: Target detection Matrix CFAR sKL mean sKL median tKL t center

a b s t r a c t Target detection in clutter is a fundamental problem in radar signal processing. When the received radar signal contains only few pulses, it is difficult to achieve a satisfactory performance using the traditional detection algorithm. In recent times, a generalized constant false alarm rate (CFAR) detector on the Riemannian manifold of Hermitian positive-definite (HPD) matrix was proposed. The employment of this detector, which compares the Riemannian distance between the covariance matrix of the cell under test (CUT) and an average matrix of reference cells with a given threshold, has significantly improved the detection performance. However, the application of this detector in real scenarios is still limited by two problems; it is computationally expensive and the detection performance is not very good since the Riemannian distance is utilized. In this paper, the symmetrized Kullback–Leibler (sKL) and the total Kullback–Leibler (tKL) divergences, instead of the Riemannian distance, are used as dissimilarity measures in the matrix CFAR detector. According to sKL and tKL divergences, three average matrices, the sKL mean, the sKL median, and the tKL t center, are derived. Furthermore, the relationship between the detection performance and the anisotropy of the distance measure used in the matrix CFAR detector is explored. Numerical experiments and real radar sea clutter data are given to confirm the superiority of the proposed algorithms in terms of the computational complexity and the detection performance. © 2017 Elsevier Inc. All rights reserved.

1. Introduction Target detection in a clutter is of great importance in radar signal processing when only few pulses are available. However, the classical fast Fourier transform (FFT) based constant false alarm rate (CFAR) detection algorithms [1] cannot achieve satisfactory performance owing to the poor Doppler resolution as well as the energy spread of the Doppler filter banks. To circumvent these drawbacks, many strategies are conceived to cope with such situations. For instance, in [2–6], the authors exploit the priori information about the surrounding environment to achieve significant performance improvements. Another example is provided in [7], wherein the Bayesian approach is employed to assume a suitable distribution about the unknown clutter covariance matrix, and similar methods are found in [8–10]. In recent times, a generalized CFAR technique on the Riemannian manifold of the Hermitian positive-definite (HPD) matrix, referred to as the Riemannian distance-based matrix CFAR, was proposed by F. Barbaresco [11–13]. In this algorithm, the received radar complex signal z (z = { z0 , z1 , . . . , zn−1 } represents the set of

*

Corresponding author. E-mail address: [email protected] (X. Hua).

http://dx.doi.org/10.1016/j.dsp.2017.06.019 1051-2004/© 2017 Elsevier Inc. All rights reserved.

n pulses) in each cell in one coherent processing interval (CPI) is modeled as a complex circular multivariate Gaussian distribution with zero mean, which can be represented by its covariance matrix R = E[ zz H ]. The detection procedure of the Riemannian distancebased matrix CFAR detector can be formulated as follows: For each cell under test (CUT), the Riemannian distance between the covariance matrix of CUT and the Riemannian mean or median matrix of reference cells around the CUT is computed; if this distance is greater than a given threshold determined by the Monte Carlo method in order to maintain the false alarm constant, then one can conclude that there is a target at the location of the CUT. It has been proved that the Riemannian distance-based matrix CFAR detector has better detection performance than the classical FFT-CFAR detection algorithm [12,13]. However, there are two shortcomings in the Riemannian distance-based matrix CFAR detector: the computational cost is very high; the detection performance is not very good. To meliorate the above deficiencies, a straight conception is to explore some alternative distance measures on the Riemannian manifold instead of the Riemannian distance. Simultaneously, these distance measures have lower computational complexity and better detection performance than the Riemannian distance used in the matrix CFAR detector. In addition to the Riemannian distance,

X. Hua et al. / Digital Signal Processing 69 (2017) 106–116

many of the measure structures can be given on the Riemannian manifold of HPD matrices. For instance, the square loss function has been applied to regression analysis; the Bhattacharyya divergence is exploited for diffusion tensor magnetic resonance image (DT-MRI) segmentation [14,15]; and the Kullback–Leibler (KL) divergence [16] has been used for measuring the dissimilarity between two probability density functions. The KL divergence is an information measure already widely used in many detection problems, even involving radar signals. In the multi-target recognition application, the mutual information and the KL divergence are used for the multiple-input multiple-output (MIMO) radar optimal waveform design [17]. In [18], the authors employ the KL divergence to measure the difference between the probability densities of the observations under two alternative hypotheses. Moreover, the difference is applied to address the problem of space–time code design for a MIMO radar detection system. In the synthetic aperture radar (SAR) image change detection application, the problem of change detection on SAR images acquired at different dates, is addressed via the non-parametric estimation [19] and the parametric estimation [20,21] of the KL divergence is used as a distance between statistical distributions of local features. Moreover, the KL divergence has been applied to anomaly detection with the wall human detection [22] and the monitoring of large-scale technical systems [23], Stereo matching [24], incipient fault detection [25], and the optimal distribution approximation [26]. In radar target detection applications, our previous work [27] and Zhao [28] have explored the matrix CFAR detection method based on the KL divergence. Numerical experiments [27,28] and X-band radar clutter [12] are given to confirm that the detection performance of KL divergence-based matrix CFAR detector outperforms the Riemannian distance-based matrix CFAR detector, and both these two matrix CFAR detection algorithms have better detection performance than the classical FFF-CFAR. Among many modern radar applications, a class of detector is based on the Cell-Averaging (CA) CFAR technique. The cellaveraging process is an example of what is known as a sliding window detector. It is based upon the properties of the Gaussian processes instead of the Neyman–Pearson Lemma. The matrix CFAR detector [12,27] and the classical CA-CFAR detector [1] are of similar schemes under the constant false alarm rate formulation. The slow-time dimension of the received clutter data in each cell is modeled and represented as a HPD matrix. This HPD matrix denotes the correlation or power level between multiple pulses, and is estimated by the pulses data according to its correlation coefficient. The matrix CFAR detector differs from the CA-CFAR detector in three ways: (1) the observation data in each cell is a HPD matrix, and not the original pulse data; (2) the distance measures used in matrix CFAR detector is Riemannian distance or divergence measure, and not the Euclidean distance; and (3) the averaging process in matrix CFAR detector is geometric mean of a set of HPD matrices, not the arithmetic mean of scalar number. These differences imply that the matrix CFAR detector performs on the HPD matrix space, in other words, the different geometry considered in detection. In this paper, we develop the detection algorithms in the framework of the matrix CFAR detector based on the symmetrized Kullback–Leibler (sKL) divergence [29] and the total Kullback– Leibler (tKL) divergence [30]. According to sKL and tKL divergences, three average matrices, the sKL mean [31], median, and the tKL t center [30], are derived and explored to replace the Riemannian mean and median. Furthermore, we explore the relationship between the detection performance of the matrix CFAR detectors based on different distance measures and the anisotropy of the distance measure. We show that different distance measures have different anisotropy, and the distance measure, which has a better anisotropy, has a better detection performance. These results

107

would have been proven by numerical experiments and real sea clutter data. The rest of this paper is organized as follows: Section 2 gives a brief description of the Riemannian distance-based matrix CFAR detector; the detector proposed in this paper is detailed in Section 3; results obtained from simulated data are presented in Section 4; and Section 5 concludes our work. 1.1. Notation A lot of notations are adopted as follows. We use math italic for scalars x, uppercase bold for matrices A, and lowercase bold for vectors x. The conjugate transpose operator is denoted by the symbol (·) H . tr(·) and det(·) are the trace and the determinant of the square matrix argument, respectively. I denotes the identity matrix, and Cn , H(n) are the sets of n-dimensional vectors of complex numbers and of n × n Hermitian matrices, respectively. The Frobenius norm of the matrix A is denoted by  A  F . For any A ∈ H(n), A > 0 means that A is a HPD matrix, and denoted by P(n). Finally, E(·) denotes statistical expectation. 2. Riemannian distance-based matrix CFAR detector 2.1. Signal model and signal manifold For the radar received complex clutter data z = { z0 , z1 , . . . , zn−1 } in each cell in one CPI, where n is the length of pulses, assuming z is a complex circular multivariate Gaussian distribution, z ∼ C N (0, R ), with zero mean and covariance matrix R,

p(z| R ) =



1

π n det( R )

exp − z H R −1 z



(1)

with the covariance matrix R given by

⎡ 



r0 ⎢ r1 ⎢

R = E zz H = ⎢

. ⎣ ..

r n −1

⎤ · · · r¯n−1 · · · r¯n−2 ⎥ ⎥ .. ⎥ , .. .. . . . ⎦ · · · r1 r0 r¯1 r0

rk = E[ zi z¯ i +k ], 0 ≤ k ≤ n − 1, 0 ≤ i ≤ n − 1

(2)

where rk = E[ zn z¯ n+k ] is called the correlation coefficient and z¯ denotes the complex conjugate of z. R is a Toeplitz HPD matrix with R H = R. The estimation of covariance matrix R is called the direct data domain method, which has been used in array signal processing [32] and space time adaptive processing [33]. It is well known that the stationary Gaussian processes have both ergocity and strict stationarity. According to the ergodicity, the correlation coefficient rk of data z can be calculated by averaging it over time instead of its statistical expectation E[ zn z¯ n+k ], as

rˆk =

1 n−k

n− 1−k

z j z¯ j +k ,

0≤k≤n−1

(3)

j =0

The pulse data in each cell in one CPI is modeled by Eqs. (1) and (2), and are represented by an HPD matrix. This matrix stands for the correlation or power level between multiple pulses. Through parameterization by the HPD matrix, the received clutter data in each cell in one CPI z = { z1 , z2 , . . . , zn } can be mapped into an n dimensional parameter space,

ψ : Cn → P(n),

z → A ∈ P(n)

(4)

Here, P(n) forms a differentiable Riemannian manifold [34]. In the following, we present the distance measure and its mean and, median on manifold.

108

X. Hua et al. / Digital Signal Processing 69 (2017) 106–116

2.2. Riemannian mean and median For the two given models C N (0, R 1 ) and C N (0, R 2 ), the Riemannian distance between the two HPD matrices R 1 , and R 2 is given as [12]



−1/2 

2

−1/2

d2 ( R 1 , R 2 ) = log R 1

R2 R1

F

=

n

log2 (λk )

(5)

k =1

−1/2

−1/2

where λk is the kth eigenvalue of R 1 R2 R1 , log(·) is the logarithmic map on the Riemannian manifold of HPD matrices. With respect to the Riemannian distance, the mean or median of a set of matrices determined by the distance metric used in the detector has been widely studied in the last few decades [35–38]. M. Moakher has done pioneering work [35] to define the metricbased mean for more than two positive definite matrices associated with the Riemannian distance (or divergence). X. Pennec [39] defined the geometric mean, and obtained the geometric mean using a gradient descent algorithm. It is proved that the geometric mean for a set of HPD matrices exist and is unique. A simple and fast convergence algorithm to compute the geometric mean is proposed using the fixed-point algorithm [40]. The Riemannian mean of N given HPD matrices { R 1 , R 2 , . . . , R N } is defined as the minimal value of the sum of the squared distances to the given matrices R i ,

¯ = arg min R

N

R

d2R ( R , R i )

(6)

i =1

where d R (·,·)represents the Riemannian distance. The existence and uniqueness of the solution in Eq. (6) was proven in [41]. For a set of HPD covariance matrices { R 1 , R 2 , . . . , R N }, the Riemannian mean can be derived using the fixed-point algorithm or the gradient descent algorithm [39]. The iteration formula is [12],





¯ n+1 = R¯ n1/2 exp −εn R

N





¯ n−1/2 R k R¯ n−1/2 log R



¯n R

1/2

(7)

k =1

where εn is the step size, it varies with n or a constant. For a set of HPD matrices { R 1 , R 2 , . . . , R N }, the Riemannian median is defined as the minimal value of the sum of the distances to the given matrices R i ,

ˆ = arg min R

N

R

dR (R, Ri )

(8)

i =1

Similar to Eq. (6), the Riemannian median matrix can be obtained using the fixed-point algorithm as the following iteration [13], 1/2

ˆ n+1 = Rˆ n R 

 

exp tn

k∈G R˜

−1/2

ˆn Cn = log R

−1/2 

ˆn Rk R

,



Cn

Cn  F

1/2

ˆn , R

Fig. 1. Riemannian distance-based matrix CFAR detector.

Then, the distance between the covariance matrix R D of the CUT ¯ of the reference cells around and the mean or median matrix R the CUT is calculated. Finally, the detection is performed by com¯ for a given threshold γ . paring the distance between R D and R The rule of target detection can be formulated as follows,

ˆ) d( R D , R

target present

≷

n

3. Proposed matrix CFAR detectors In this section, the matrix CFAR detectors based on the sKL and tKL divergences are presented in detail. The sKL and tKL divergences are used as a replacement of the Riemannian distance to measure the dissimilarity between the two HPD matrices on the Riemannian manifold of the HPD matrix. Based on the sKL and tKL divergences, the sKL mean and median, and the tKL t center are derived. Fortunately, there are closed-form expressions for the sKL mean and the tKL t center in the HPD matrix space. Moreover, we analyze the anisotropy of these average matrices compared to that of the Riemannian distance. In the performance assessment part, we will prove the coherence between the detection performance and the anisotropy of distance measure. 3.1. sKL and tKL divergences 3.1.1. sKL divergence In information geometry, the dissimilarity between two probability distributions p and q on a space Ω is usually measured by the KL divergence (or relative entropy) [43],



(9)

where tn is the step size, and the computations of the Riemannian mean and median are referred to [37,42] for more details. 2.3. Detection According to the Riemannian distance and its corresponding mean and median presented in the above section, the Riemannian distance-based matrix CFAR detector can be illustrated as shown in Fig. 1. The pulse data in each cell in one CPI is a HPD matrix estimated by the sample data z according to its correlation coefficient.

(10)

ˆ ) is the test statistic, and d(·,·) denotes the Riewhere d( R D , R mannian distance. Since there is no analytical expression for the threshold γ , the threshold γ is derived by the Monte Carlo method in order to maintain the false alarm constant. Moreover, the different measures exploited in the detector have different detection performance. The detection performance is related to anisotropy of the geometric measure, and this will be illustrated in numerical experiments.

n

ˆ n} G R˜ = {k/ R k = R

γ

target absent

KL( p , q) =

p (x) log

p ( x) q ( x)

dx

(11)

Ω

For two zero-mean complex Gaussian distributions,

p ( x| R 1 ) = q ( x| R 2 ) =

1

π n det( R 1 ) 1

π

n det( R

2)









1 exp −x H R − 1 x 1 exp −x H R − 2 x

(12)

The KL divergence between them whose covariant matrices are R 1 and R 2 gives rise to the KL divergence for the two HPD matrices R 1 and R 2 [43],

X. Hua et al. / Digital Signal Processing 69 (2017) 106–116

109

Fig. 2. Difference between the tBD and BD divergences.







1 −1 KL( R 1 , R 2 ) = tr R − 2 R 1 − I − log det R 2 R 1



(13)

Here, a fact must be emphasized that the KL divergence does not define a distance as it is neither symmetric with respect to its two arguments nor does it satisfy the triangle inequality. Its symmetrized form can be defined as [29],

sKL( R 1 , R 2 ) =

1



KL( R 1 , R 2 ) + KL( R 2 , R 1 ) 2  1  1 −1 = tr R − 2 R 1 + R 1 R 2 − 2I 2

(14)

From (16) it becomes clear that the sKL divergence behaves as the square of a distance; however, it is not a distance, as the triangle inequality does not hold. 3.1.2. tKL divergence The tKL divergence is a special case of the total Bregman divergence (tBD), which was proposed by B.C. Vemuri [30], and it was a robust divergence measure. The tBD d f associated with a strictly convex and differentiable function f defined on a convex set X between points x, y ∈ X is defined as [30],

d f (x, y ) =

f (x) − f ( y ) − x − y , ∇ f ( y )



1 + ∇ f ( y )2

(15)

where ∇ f ( y ) is the gradient of f at y and ·,· is the inner product on the space X . Let p, q be the two  probability distributions, the tBD becomes tKL when f ( p ) = p (x) log p (x)dx. The difference between the tBD and the Bregman divergence (BD) is that tBD measures the orthogonal distance between the value of a convex and a differentiable function f at the first argument and its tangent at the second argument, while the BD is a measure of the distance that is perpendicular to the horizontal axis. This can be observed from Fig. 2. If the coordinate system rotates an angle, BD(x, y ) changes, while the tBD(x, y ) does not. This means that tBD is invariant to linear transformation and that it is a statistically robust divergence. For the given two complex Gaussian distributions p and q with zero-mean, the tKL divergence between them gives rise to the tKL divergence for the two HPD matrices R 1 and R 2 [30], 1 −1 log(det( R − 1 R 2 )) + tr ( R 2 R 1 ) − n tKL( R 1 , R 2 ) =  (log(det R 2 ))2 2π ) 2 c+ − n(1+log log(det R 2 ) 4 2

Fig. 3. Plots of the iso-surface of different distance measures.

vectors span an n-dimensional isotropic unit sphere. However, this sphere might not be isotropous at a location on the Riemannian manifold, since the curvatures along the directions of the unit basis vectors of this location are not identical. In order to illustrate the difference of anisotropy of the aforementioned distance measures, we provide some visual results of the Riemannian distance, sKL divergence, and tKL divergence in the case of tree-dimensional HPD matrix. These anisotropy isosurfaces, centered at the identity matrix with radius r = 0.1, 0.5, and 1, are shown in Fig. 3. Given the identity matrix I and a three dimensional matrix P , the distance measures can be expressed as,

d R ( I , P ) =  log P  F = r

1

sKL( I , P ) =

2



P 1/2 − P −1/2

=r F

log(det( P )) + tr( P −1 ) − n

tKL( I , P ) =  =r (log(det( P )))2 n(1+log 2π ) 2 c+ − log ( det ( P )) 4 2 (16)

where c is (3n/4) + (n2 log 2π /2) + ((n log 2π )2 /4), and n stands for the dimension of matrix R. 3.1.3. Anisotropy of distance measures The anisotropy of the distance measure indicates the difference between the normalized basis vectors on the Riemannian manifold. In the n-dimensional Euclidean space, the normalized basis

(17) It can be observed from Fig. 3 that the anisotropy of various measures are different at point I on manifold P(3). Moreover, the iso-surface of each measure is an anisotropy sphere. This also implies that the curvatures along the different directions are not identical. The curvature usually reflects the structure of manifold. In numerical experiments, we will give a quantitative illustration about the difference of anisotropy at different locations.

110

X. Hua et al. / Digital Signal Processing 69 (2017) 106–116

3.2. Divergence-based mean and median 3.2.1. sKL mean and median According to the sKL divergence, similar to the definition of (4) and (6), we can define the sKL mean and median. Definition 1. The sKL mean, relative to sKL divergence, of N HPD matrices { R 1 , R 2 , ..., R N } is defined as

¯ sKL = arg min R

N

R ∈P(n)

sKL( R , R i )

(18)

i =1

Definition 2. The sKL median, associated with the sKL divergence, of N HPD matrices { R 1 , R 2 , ..., R N } is defined as

ˆ = arg min R

N 

R ∈P(n)

∗

R = sKL( R , R i )

(19)

i =1

Thus, the gradient of f ( R ) =

N √ i =1

N

i =1 sKL( R ,

R i ) and g ( R ) =

sKL( R , R i ) at R can be given by,



∇ f (R) = R2

N

 R k−1

−

N

Ri

N −1 1 − R −2 R i + R i

2N

(20)

sKL( R , R i )

i =1

The equation ∇ f ( R ) = 0 and ∇ g ( R ) = 0 both have unique solutions on the space of HPD matrices. Its existence and uniqueness can be proved using Banach’s fixed-point theorem [44]. The sKL mean and median can be obtained using the fixed-point algorithm,

 ¯ = R

N 1

N

Ri

i =1

 ˆ k +1 = R



N

N 1 −1 Rk N k =1

Ri



ˆ k, Ri) sKL( R i =1

−1 1/2

N

−1

Rj

(21)

It can be noted from (23) that the sKL mean has a closedform expression while the solution of sKL median is the iteration equation. An interesting fact should be emphasized that the sKL mean is the geometric mean of the arithmetic mean and the harmonic mean of the N given  NHPD matrices { R 1 , R 2 , . . . , R N } [29]. Moreover, the minimize i =1 KL( R i , R ) is the arithmetic mean

N

i =1 KL( R ,

R i ) is the har-

3.2.2. tKL t center The tKL t center defined by B. Vemuri [31] is efficient to compute, which is a suitable characteristic. In many applications, the tKL t center acts as a robust cluster center for clustering and classification owing to its remarkable anisotropy. Given a set of HPD matrices { R 1 , R 2 , . . . , R N }, the tKL t center is defined as the following minimization,

R ∗ = arg min R ∈P( N )

N

tKL( R , R i )

(22)

i =1

∗  NTo find R , one can take the derivative of f ( R ) = i =1 tKL( R , R i ) with respect to R, and set it to 0. Then the t

center R ∗ can be given by [31],

−1

wi Ri

μi j μj

wi = 

,

i





μi = 2 c +

(log(det R i ))2 4

−

n(1 + log 2π ) 2

 −1 log(det R i ) (23)

trix R. It can be noted from (23) that the weight w i is different from each other, and it is related to its matrix closely. Another significant property of the t center is that the weight is inversely proportional to the value of the gradient of the convex function ∇ f ( R ) used in defining the divergence. As noise data and outliers have a greater magnitude of the gradient, the t center is robust to clutter. This is important for our proposed matrix CFAR detector. As the divergence measure and average matrix are presented above, our proposed matrix CFAR detectors can be derived as illustrated in Fig. 4. Our proposed detectors are more competitive than the Riemannian distance-based matrix CFAR detector in terms of both the computational complexity and the detection performance. 4. Performance analysis

−1 1/2

ˆ k, R j) sKL( R j =1

of { R 1 , R 2 , . . . , R N }, and the minimize monic mean of { R 1 , R 2 , . . . , R N }.

−1

where w i is the weight of the ith matrix R i , and c is(3n/4) + (n2 log 2π /2) + ((n log 2π )2 /4), n denotes the dimension of ma-

i =1

k =1

∇ g( R ) =

Fig. 4. Proposed matrix CFAR detector.



In this section, numerical experiments and real sea clutter data are performed to evaluate the detection performance of the proposed matrix CFAR detectors in comparison with the Riemannian distance-based matrix CFAR detector. 4.1. Numerical experiments The numerical examples are obtained by means of the standard Monte Carlo counting techniques, since the analytical expressions for detection probability ( P d ) and probability of false alarm ( P f a ) are not available. The data samples are from N = 7 received pulses. A total of 50 range cells are considered, while M = 16 range cells are used for averaging. One moving target with a Doppler frequency f d = 3.83 Hz is located at the 13th range cell. We use K distribution with the shape parameter ν = 1, scale parameter μ = 0.5 to simulate the clutter, which can properly fit the sea clutter data [45]. The amplitude-pdf of K-distribution is

√

2 ν /μ

p (x) = ν −1 2 Γ (ν )



2ν

μ



ν x

K ν −1

2ν

μ

 x ,

x≥0

(24)

where Γ (·) is the gamma function, and K ν −1 (·) is the modified Bessel function of the second kind with order ν − 1. Here, the anisotropy index at the location R, related to the geometric measure, is defined to be the distance to its closest isotropic matrix as L = min d( R , α I ), where α is a scalar parameter. The normalized anisotropic index N L, defined as anisotropy

X. Hua et al. / Digital Signal Processing 69 (2017) 106–116

111

Fig. 5. Plot of normalized anisotropy index in 50 range cells.

factor L divided by the maximum value of 50, is plotted in Fig. 5. It can be noted from Fig. 5 that different locations have different anisotropy, and different measures at a location have different anisotropy. The normalized anisotropic indexes in the range cells with the target are much larger than that in the range cells without the target in these measures. Moreover, the tKL divergence has the best anisotropy, followed by the sKL divergence, while the Riemannian distance has the worst anisotropy. The pulse data in each cell in one CPI is modeled as a circular complex Gaussian distribution, and a HPD matrix is the new observation data for each cell. The statistics in each range cell, which are the distance measures between the covariance matrix of CUT and its corresponding mean or median matrix, are plotted in Fig. 6 when Signal to Clutter Rate (SCR) equals to 0 dB, 5 dB, and 10 dB. It can be seen from Fig. 6 that the statistics in the range cells with the target are much larger than that in the range cells without the target in these detectors. This means that the target can be detected by all the five algorithms correctly. However, the normalized statistics of these algorithms in the range cells without the target are quite different. It is clear that the statistics of the range cells without the target for sKL mean, sKL median are much lower than that of the Riemannian mean and median; and the tKL t center have the lowest statistics in the range cells without the target. This indicates that both the sKL mean and sKL median have better detection performance than the Riemannian mean and median-based matrix CFAR; and the tKL t center has the best detection performance. It can be also noted that the statistics of the range cells with the target for the tKL t center-based matrix CFAR are the largest followed by the sKL mean and sKL median-based matrix CFAR when compared to the range cells without the target. This phenomenon is consistent with the anisotropy related to the Riemannian distance, sKL divergence, and tKL divergence plotted in Fig. 5. Furthermore, it can be observed from Fig. 6(a) that the statistics of the sKL mean and Riemannian mean detectors in some range cells are very close. The reason is that their detection probabilities are very close to each other and are very low in the SCR = 0 dB scenario. This can be found from Fig. 7 in the next. Monte Carlo simulations are carried out to give a more accurate comparison of the detection performance between these algorithms. For the underlying detectors, the detection thresholds are not only related to the background clutter power, but are also associated with the threshold multiplier, which usually does not exist in a closed-form expression. In consequence, it is difficult to determine the detection threshold analytically. In the simulation, the detection threshold is obtained via Monte Carlo runs. The test statistics in the absence of a target is calculated via 106 Monte

Carlo runs, where the detection threshold is determined according to the given probability of false alarm P f a . To ensure the detection probability is accurately estimated, we repeat the simulations 200 times to estimate the detection probability for different values of SCR, keeping the noise power constant and varying the value of the signal energy. The detection probability is estimated by the relative frequencies. As the adaptive matched filter (AMF) is often used for the target detection in clutter, our proposed detectors are also compared to the AMF algorithm with the covariance matrix estimated by Eq. (2). Fig. 7 plots the probability of detection versus SCR under different probabilities of false alarm. The SCR varies from −5 to 15 dB with an interval of 1 dB. The P f a for these plots are 10−5 , and 10−4 respectively. It can be noted from Fig. 7 that the tKL t center detector has the best detection performance followed by the sKL divergence-based matrix CFAR detector. These detectors have better detection performance than the Riemannian distance-based matrix CFAR detector. Moreover, the detection performances of these matrix CFAR detectors are better than that of the AMF detection algorithm. The P d of tKL t center detector has 4–5 dB higher than the Riemannian mean and median detector. The P d of the sKL mean and median detectors is 1–2 dB higher than that of the Riemannian mean and median detectors, respectively. These results prove the superiority of our proposed matrix CFAR detectors. The computational complexity of these average matrices derived from different distance measures is analyzed by counting the number of multiplications. Specifically, the Riemannian mean and median are more demanding, as it involves an iterative procedure and the number of multiplications are (3kn3 + kn2 ) pk and (3kn3 + 2kn2 ) pk , respectively, where k denotes the number of HPD matrices used for averaging, pk is the number of iteration, and n denotes the length of the pulse data in the range cell. The sKL median also requires an iterative procedure, but it involves simple matrix inversions. The number of multiplications are ((k + 3)n3 + 2kn2 ) pk . The sKL mean and tKL t center have the lowest computational cost as they have closed-form expressions. The number of multiplications of the sKL mean and tKL t center are (k + 3)n3 and (k + 1)n3 + 2kn2 , respectively. In order to compare the computational efficiency between these matrix CFAR detectors, Table 1 shows the time consumption during the calculation of average matrices with different number of matrices (M) used for averaging and different number of pulses (N). It is clear from Table 1 that the time taken for sKL Mean, sKL Median, and tKL t center are less than the Riemannian Mean and Median.

112

X. Hua et al. / Digital Signal Processing 69 (2017) 106–116

Fig. 6. Plots of the statistics in each cell with SCR = 0, 5, 10 dB.

4.2. Real sea clutter data The measured data we use from the McMaster University IPIX radar, which was collected at the Osborne Head Gunnery Range (OHGR), Dartmouth, Nova Scotia, Canada [46]. We exploit the file 19931107_135603_starea.cdf, 19931107_141630_starea.cdf, 19931107_145028_starea.cdf, and 19931108_213827_starea.cdf to test the detection performance of Matrix CFAR detectors. The complex data is composed of 131072 pulses and 14 range cells. The

average SCR varies in the range bin of 0–6 dB. The target parameters are detailed in Table 2. The table lists the range cell where the target is strongest (“primary”), and neighboring range cells where the target may also be visible (“secondary”). The weak static target is a spherical block of Styrofoam, wrapped with wire mesh. The observation data that was previously processed by ipixload.m consists of two steps: Remove mean and standard deviation from I and Q channels separately and remove the phase imbalance. Then, the preprocessed data in each cell is modeled and

X. Hua et al. / Digital Signal Processing 69 (2017) 106–116

113

Fig. 7. P d versus SCR with different detectors.

Table 1 Time taken of different average matrices. Average matrix

Time (ms)

M

30

N

7

14

21

7

14

21

7

14

21

Riemannian Mean Riemannian Median sKL Mean sKL Median tKL t center

547 590 4.2 89 4.5

1852 3020 6 160 6.5

3149 8230 6.1 265 7.2

703 753 4.6 120 4.9

2270 3829 6.2 1919 6.8

4202 8089 6.7 2438 8.4

1019 889 5.5 157 5.8

3603 7275 6.4 265 7.9

6086 17772 9.1 401 11.6

40

60

Table 2 Real target parameters with IPIX radar. File name

Primary range cell

Secondary range cells

Range cells

19931107_135603_starea 19931107_141630_starea 19931107_145028_starea 19931108_213827_starea

9 9 8 7

8–11 8–11 7–9 6–8

14 14 14 14

described as a HPD matrix. By comparing the statistics of these matrix CFAR detectors, we can investigate in which range cell targets are visible in Fig. 8. It is clear from Fig. 8 that the primary target cells can be well detected by these detectors. However, the normalized statistics of range cells without the target are quite different. The tKL t center detector has the lowest statistic, and the statistics of the sKL divergence-based matrix CFAR detectors is slightly smaller than

that of the Riemannian distance-based matrix CFAR detectors in the range cells without the target. This also verifies that the tKL divergence-based matrix CFAR detector has the best detection capability, followed by the sKL divergence-based matrix CFAR detectors, which are better than the Riemannian distance-based matrix CFAR detectors. In some cases, the sKL Median performs worse than the Riemannian median. A high number of outliers or discontinuity in the measured data might be the reason.

114

X. Hua et al. / Digital Signal Processing 69 (2017) 106–116

Fig. 8. Plots of statistics with mean, median detectors in range cells.

5. Conclusion In this paper, we deal with the problem of target detection in the clutter, when the received radar signal contains only few pulses. Matrix CFAR detectors based on sKL and tKL divergences are proposed. The advantages of the proposed matrix CFAR detectors are twofold: (1) computational complexity of the proposed detectors is much lower than the Riemannian distancebased matrix CFAR detector. As there are explicit expressions for the sKL mean and tKL t center, the calculation of the sKL median

merely involves some matrix inversions; (2) the proposed detectors have a higher P d than the Riemannian distance-based matrix CFAR detector. At the analysis stage, we have assessed the detection performance of the matrix CFAR detector and the anisotropy of measures by means of numerical experiments and real sea clutter. The results have shown that the better the anisotropy of a measure, the better the detection performance it exhibits. A possible future research track might concern the performance analysis of the matrix CFAR detectors from the manifold structure.

X. Hua et al. / Digital Signal Processing 69 (2017) 106–116

115

Fig. 8. (continued)

References [1] M.A. Richards, Fundamentals of Radar Signal Processing, second edition, McGraw-Hill, 2014. [2] F. Gini, M. Rangaswamy, Knowledge Based Radar Detection, Tracking and Classification, Adaptive and Learning Systems for Signal Processing, Communications and Control Series, Wiley-Interscience, 2008. [3] W.L. Melvin, J.R. Guerci, Knowledge-aided signal processing: a new paradigm for radar and other advanced sensors, IEEE Trans. Aerosp. Electron. Syst. 42 (2006) 983–996. [4] C. Hao, D. Orlando, X. Ma, C. Hou, Persymmetric Rao and Wald tests for partially homogeneous environment, IEEE Signal Process. Lett. 19 (2012) 587–590. [5] A. Aubry, A. De Maio, L. Pallotta, A. Farina, Radar detection of distributed targets in homogeneous interference whose inverse covariance structure is defined via unitary invariant functions, IEEE Trans. Signal Process. 61 (2013) 4949–4961. [6] A.D. Maio, D. Orlando, C. Hao, G. Foglia, Adaptive detection of point-like targets in spectrally symmetric interference, IEEE Trans. Signal Process. 64 (2015) 3207–3220. [7] A. De Maio, A. Farina, G. Foglia, Design and experimental validation of knowledge-based constant false alarm rate detectors, IET Radar Sonar Navig. 1 (2007) 308–316. [8] P. Wang, H. Li, B. Himed, Knowledge-aided parametric tests for multichannel adaptive signal detection, IEEE Trans. Signal Process. 59 (2011) 5970–5982. [9] A. De Maio, A. Farina, G. Foglia, Knowledge-aided Bayesian radar detectors & their application to live data, IEEE Trans. Aerosp. Electron. Syst. 46 (2010) 170–183. [10] A. Aubry, V. Carotenuto, A.D. Maio, G. Foglia, Exploiting multiple a priori spectral models for adaptive radar detection, IET Radar Sonar Navig. 8 (2014) 695–707. [11] F. Barbaresco, New foundation of radar Doppler signal processing based on advanced differential geometry of symmetric spaces: Doppler matrix CFAR and radar application, in: International Radar Conference, Radar ’09, Bordeaux, France, 2009. [12] J. Lapuyade-Lahorgue, F. Barbaresco, Radar detection using Siegel distance between autoregressive processes, application to HF and X-band radar, in: 2008 IEEE Radar Conference, 2008, pp. 1–6. [13] F. Barbaresco, Robust statistical radar processing in Fréchet metric space: OSHDR-CFAR and OS-STAP processing in Siegel homogeneous bounded domains, in: 12th International Radar Symposium, IRS 2011, 2011, pp. 639–644. [14] M. Charfi, Z. Chebbi, M. Moakher, B.C. Vemuri, Using the Bhattacharyya mean for the filtering and clustering of positive-definite matrices, in: F. Nielsen, F. Barbaresco (Eds.), Geometric Science of Information: First International Conference. Proceedings, GSI 2013, Paris, France, August 28–30, 2013, Springer, Berlin, Heidelberg, 2013, pp. 551–558. [15] M. Charfi, Z. Chebbi, M. Moakher, B.C. Vemuri, Bhattacharyya median of symmetric positive-definite matrices and application to the denoising of diffusiontensor fields, in: 2013 IEEE 10th International Symposium on Biomedical Imaging, 2013, pp. 1227–1230. [16] S. Kullback, R.A. Leibler, On information and sufficiency, Ann. Math. Stat. 22 (1951) 79–86. [17] L. Wang, K.K. Wong, H. Wang, Y. Qin, MIMO radar adaptive waveform design for extended target recognition, Int. J. Distrib. Sens. Netw. 2015 (2015) 1–11. [18] E. Grossi, M. Lops, Kullback–Leibler divergence region in MIMO radar detection problems, in: International Conference on Information Fusion, 2012, pp. 896–901.

[19] S. Cui, C. Luo, Feature-based non-parametric estimation of Kullback–Leibler divergence for SAR image change detection, Remote Sens. Lett. 7 (2016) 1102–1111. [20] Q. Xu, L.J. Karam, Change detection on SAR images by a parametric estimation of the KL-divergence between Gaussian mixture models, in: IEEE International Conference on Acoustics, Speech and Signal Processing, 2013, pp. 2109–2113. [21] J. Inglada, Change detection on SAR images by using a parametric estimation of the Kullback–Leibler divergence, in: IEEE International Geoscience and Remote Sensing Symposium, IGARSS ’03 2003, in: Proceedings, vol. 4106, 2003, pp. 4104–4106. [22] W. Wang, B. Zhang, D. Wang, Y. Jiang, S. Qin, L. Xue, Anomaly detection based on probability density function with Kullback–Leibler divergence, Signal Process. 126 (2016) 12–17. [23] J. Zeng, U. Kruger, J. Geluk, X. Wang, L. Xie, Detecting abnormal situations using the Kullback–Leibler divergence, Automatica 50 (2014) 2777–2786. [24] J. Liang, S.J. Maybank, Y. Zhang, Stereo matching-based definition of saliency via sample-based Kullback–Leibler divergence estimation, Mach. Vis. Appl. 26 (2015) 607–618. [25] A. Youssef, C. Delpha, D. Diallo, An optimal fault detection threshold for early detection using Kullback–Leibler divergence for unknown distribution data, Signal Process. 120 (2015) 266–279. [26] G.V. Weinberg, Kullback–Leibler divergence and the Pareto–Exponential approximation, SpringerPlus 5 (2016) 1–9. [27] Y. Cheng, X. Hua, H. Wang, Y. Qin, X. Li, The geometry of signal detection with applications to radar signal processing, Entropy 18 (2016) 381. [28] X. Zhao, S. Wang, An improved matrix CFAR detection method base on KL divergence, J. Electron. Inform. Technol. 38 (2016) 934–940. [29] M. Moakher, P.G. Batchelor, Symmetric positive-definite matrices: from geometry to applications and visualization, in: J. Weickert, H. Hagen (Eds.), Visualization and Processing of Tensor Fields, Springer, Berlin, Heidelberg, 2006, pp. 285–298. [30] B.C. Vemuri, M. Liu, S.I. Amari, F. Nielsen, Total Bregman divergence and its applications to DTI analysis, IEEE Trans. Med. Imaging 30 (2011) 475–483. [31] W. Zhizhou, B.C. Vemuri, DTI segmentation using an information theoretic tensor dissimilarity measure, IEEE Trans. Med. Imaging 24 (2005) 1267–1277. [32] W. Choi, T.K. Sarkar, H. Wang, E.L. Mokole, Adaptive processing using real weights based on a direct data domain least squares approach, IEEE Trans. Antennas Propag. 54 (2006) 182–191. [33] D. Cristallini, W. Burger, A robust direct data domain approach for STAP, IEEE Trans. Signal Process. 60 (2012) 1283–1294. [34] R. Bhatia, Positive Definite Matrices, Princeton University Press, Princeton, NJ, 2007. [35] M. Moakher, A differential geometric approach to the geometric mean of symmetric positive-definite matrices, SIAM J. Matrix Anal. Appl. 26 (2005) 735–747. [36] D. Petz, Means of positive matrices: geometry and a conjecture, Ann. Math. Inform. 32 (2005) 129–139. [37] P.T. Fletcher, S. Venkatasubramanian, S. Joshi, Robust statistics on Riemannian manifolds via the geometric median, in: 2008 IEEE Conference on Computer Vision and Pattern Recognition, CVPR, 2008, pp. 1–8. [38] L. Yang, Riemannian median and its estimation, LMS J. Comput. Math. 13 (2009) 461–479. [39] X. Pennec, Intrinsic statistics on Riemannian manifolds: basic tools for geometric measurements, J. Math. Imaging Vis. 25 (2006) 127–154. [40] M. Moakher, On the averaging of symmetric positive-definite tensors, J. Elast. 82 (2006) 273–296.

116

X. Hua et al. / Digital Signal Processing 69 (2017) 106–116

[41] R. Bhatia, J. Holbrook, Riemannian geometry and matrix geometric means, Linear Algebra Appl. 413 (2006) 594–618. [42] D.A. Bini, B. Iannazzo, Computing the Karcher mean of symmetric positive definite matrices, Linear Algebra Appl. 438 (2013) 1700–1710. [43] S.I. Amari, Differential-Geometrical Methods in Statistics, Springer-Verlag, 1985. [44] Z. Chebbi, M. Moakher, Means of Hermitian positive-definite matrices based on the log-determinant α -divergence function, Linear Algebra Appl. 436 (2012) 1872–1889. [45] M. Greco, P. Stinco, F. Gini, M. Rangaswamy, Impact of sea clutter nonstationarity on disturbance covariance matrix estimation and CFAR detector performance, IEEE Trans. Aerosp. Electron. Syst. 46 (2010) 1502–1513. [46] S. Haykin, C. Krasnor, T.J. Nohara, B.W. Currie, A coherent dual-polarized radar for studying the ocean environment, IEEE Trans. Geosci. Remote Sens. 29 (1991) 189–191.

Xiaoqiang Hua received the B.S. in surveying and mapping engineering from Wuhan University, Wuhan, China, in 2008, and the M.S. degrees in photogrammetry and remote sensing from National University of Defense Technology, Changsha, China, in 2012, respectively. He is currently pursuing the Ph.D. degree from the National University of Defense Technology. His research interests lie in the areas of information geometry, statistical signal processing, photogrammetry and remote sensing, image processing and computer vision. Yongqiang Cheng received the M.S. and Ph.D. degrees in information and communication engineering from National University of Defense Technology, Changsha, China, in 2007 and 2009, respectively. He is a Professor

of the National University of Defense Technology. His research interests lie in the areas of information geometry and statistical signal processing. Hongqiang Wang is a Professor of the National University of Defense Technology. His research interests lie in the areas of statistical signal processing, and radar imaging. Yuliang Qin is a Professor of the National University of Defense Technology. His research interests lie in the areas of statistical signal processing, and radar imaging. Yubo Li received the B.S. in University of Electronic Science and Technology, Chengdu, China, in 2006, and the M.S. degrees in National University of Defense Technology, Changsha, China, in 2010, respectively. He is currently pursuing the Ph.D. degree from the National University of Defense Technology. His research interests lie in the areas of information geometry, statistical signal processing, target tracking and computer technology. Wenpeng Zhang received the B.S. and M.S. degrees in information and communication engineering from National University of Defense Technology, Changsha, China, in 2008 and 2012, respectively. He is currently pursuing the Ph.D. degree from the National University of Defense Technology. His research interests lie in the areas of micro-motion target detection, GMTI and statistical signal processing.