doi:10.1155/2008/786136 Research Article Reiterative Robust Adaptive Thresholding for Nonhomogeneity Detection in Non-Gaussian Noise

A robust and data-dependent adaptive thresholding algorithm for nonhomogeneity detection in non-Gaussian interference is addressed. The algorithm is to be used as a preprocessing technique to select a set of homogeneous data from a bulk of nonhomogeneous compound-Gaussian secondary data employed for adaptive radar. An iterative version of the algorithm is also suggested in situations of multiple outliers in the secondary data. Performance analysis is conducted with simulated data as well as with real sea clutter data.


INTRODUCTION
The deleterious impact of nonhomogeneous secondary data in covariance estimation for adaptive radar systems has been widely reported [1][2][3][4][5][6]. Typically, the unknown interference covariance matrix is estimated from a set of identically and independently distributed (iid) target-free data, which is representative for the interference statistics in a cell under test (CUT). Frequently, the secondary data is subject to contamination by discrete scatterers or interfering targets (outliers). In both events, the training data becomes nonhomogeneous, leading to a nonrepresentative set for the interference in the CUT. Estimates of the covariance matrix from nonhomogeneous training data result in under-nulled clutter. Consequently, constant false alarm rate (CFAR) and detection performance are greatly degraded. This has led to the development of improved training data selection techniques, seeking to discard bins that are nonhomogeneous from the bulk of training data, and using the resulting outlier-free data for covariance matrix estimation.
Recently, several algorithms for outlier removal have been proposed [7][8][9]. The works of [7,9] addressed the use of the nonhomogeneity detector (NHD) based on the generalized inner product (GIP) measure involving Gaussian interference scenarios. However, in most practical situations, the Gaussian model is no longer valid. In particular, for high-resolution radars operating at low grazing angles, or in a maritime environment, a satisfactory fit of the clutter amplitude probability density function (apdf) can be achieved through families of non-Gaussian distributions [10,11]. For this non-Gaussian interference, the corresponding nonhomogeneity detection problem has received limited attention. This is due to the fact that there is no unique model for representing the joint probability density function (pdf) of a set of K-correlated non-Gaussian random variables. In [12], the interference has been modeled by a spherically invariant random process (SIRP) with the pdf of the texture assumed to be known a priori. Also, in [12] the expectation maximization (EM) algorithm is used to estimate the interference covariance matrix from representative training data, generally those in the neighborhood of the CUT, sharing the covariance structure of the CUT. The estimated covariance matrix is then used in a scale invariant test for nonhomogeneity detection, where a data-independent fixed threshold for excision is numerically calculated with respect to a type-I error criterion. The type-I error, denoted P e in this paper, is the probability of incorrectly excising homogeneous data. In real-life situations, in addition to the fact that the texture PDF is unknown, the methods cited above suffer from masking effects. This problem appears when, for example, one or multiple interfering targets are present in the vicinity of the CUT.
In the process of outlier removal, the most important issue is the setting of the threshold for excision. A survey 2 EURASIP Journal on Advances in Signal Processing of previous works in this area reveals two basic approaches. First, a prescribed number of bins corresponding to the outlier values may be removed (see [8] and references therein). This leads to an improved result, but it is not clear how to optimally choose the number of bins to be removed. In a second approach, the threshold setting is based on the knowledge of the pdf of the test statistic used for the NHD. This method works well if the assumed pdf accurately represents the data. However, the presence of nonhomogeneities can significantly alter this pdf, making the set threshold inaccurate for practical use.
This paper presents a data-dependent adaptive thresholding algorithm that can be used for nonhomogeneity detection in compound-Gaussian noise. The proposed algorithm does not require the knowledge of the texture pdf, the noise, and the number of targets to be removed. We use the bootstrap [13] to estimate the pdf of the test statistic from which the threshold for excision is set according to a prescribed type-I error P e = α. For multiple targets, an iterative version of the algorithm is proposed. The test has to be applied to nonhomogeneous secondary data in order to select a homogeneous set. The latter is to be used for estimating the noise covariance matrix in an adaptive radar detector.
The paper is organized as follows: Section 2 presents some preliminaries. In Section 3, we present the NHD test statistic based on the bootstrap principle. Performance analysis with simulated as well as real sea clutter data is shown and discussed in Section 4. A comparison with other existing methods is also included. Conclusion is given in Section 5.

PRELIMINARIES
Let z i , i = 1, . . . , K, denote K identically and independently distributed (iid) complex vectors from secondary data. Under the compound-Gaussian model, z i = [z1 z 2 · · · z N ] T can be written as [10,11] with x i ∼CN(0, R) is the speckle component following a complex normal pdf with zero mean and unknown covariance matrix R, τ i ∼ f τi (τ i ), is a positive random variable called the texture component and f τ1 ( is the texture pdf, which is unknown in this paper. One of the most popular high-resolution radar clutter models is the K-distribution model which is obtained when assuming that the texture is Gamma distributed, that is, with mean value E(τ) = μ, which is also the average clutter power, and shape parameter ν. The lower the value of ν, the spikier is the clutter. The product model of (1) corresponds to the case of texture completely correlated in azimuth. This model accurately describes the scattering mechanism for observation time intervals in the order of the coherent processing interval of the radar system [14].
As mentioned before, an adaptive radar uses the K training vectors, supposed to be homogeneous, to estimate the covariance matrix of the noise in the cell under test [15]. If outliers resembling a target of interest are present in the data, a preprocessing step is needed to excise them [16]. This can be cast in the following statistical hypothesis test: H 0 : data homogeneous, H 1 : data non-homogeneous. ( Most of the methods used in the Gaussian case for this purpose rely on the GIP given by where is the sample covariance matrix which is also the maximum likelihood estimate (ML) of R calculated from range bins surrounding the CUT [15]. The cell bins corresponding to P i > γ are declared to be outliers and discarded. The remaining ones are then used to estimate the covariance matrix. The threshold γ is set according to a certain criterion (see [8,9] and references therein for more details). For the compound-Gaussian case, the GIP statistic cannot be directly applied since in this case, the sample covariance matrix R S does not estimate correctly the true covariance matrix. Reference [17] contains relevant methods to estimate the covariance matrix of heavy tailed noise.
Our aim is to detect the presence of interfering targets in the secondary data. To this end, denote by p i the Ndimensional complex vector representing these outliers and z i the received signal from the ith range cell. The test for nonhomogeneity can then be recast in the following test: Bearing in mind that we are concerned with training data containing interfering targets, which share the same steering vector as that of the desired target, we can model p i as: p i = α i s, where s is the steering vector of interest. The test statistic, which is widely used for Gaussian as well as non-Gaussian noise for this kind of problems, is given by the adaptive coherence estimator (ACE) [12,[18][19][20] where R is an appropriate estimate of the covariance matrix.
Since the texture PDF is unknown, we should use an estimate of R, which is robust or at least independent of this PDF.

Robust estimation of the covariance matrix
In [18], a robust estimator of the covariance matrix, known as the normalized sample covariance matrix (NSCM), was proposed and given by where τ j is the estimate of the texture component in the jth bin based on the method of moments given by Inserting (9) in (8), we get The choice of the estimator in (10) is motivated by its easy implementation. Other estimators of R can be used, such as in [17,21,22], but they require much more computations.
It is interesting to note that the test statistic of (7) with R given by the NSCM of (10) is independent of the texture PDF. This is of great importance in practical situations where this PDF is generally unknown and even if it were, the corresponding parameters are always unknown and must be estimated from the data. This requires a large number of samples which is not always available. One notes that for the obtained test statistic there is no closed-form expression for its pdf. This means that we cannot set the excision threshold analytically. We will use therefore an alternative way based on the boostrap to set this threshold.

Bootstrap-based detection
In the next section, we will use the bootstrap to detect nonhomogeneities in the secondary data. We therefore give a brief review of the bootstrap-based detection. The bootstrap is a general tool to solve difficult statistical problems. It is extremely valuable in situations of unknown distributions and where data sizes are too small to invoke asymptotic results. Two principles are at its core: substituting estimates for unknowns and replacing tedious mathematical analysis with computer simulations. Figure 1 gives the basic idea of the general bootstrap detection procedure. Essentially, the bootstrap provides a data-dependent threshold for the detector. This threshold is computed from a large number of realizations of T(x * ) based on data resampled (bootstrap data) from one set of observations (original data). The resampling procedure is usually done by using a pseudorandom number generator to draw a new random sample X * of the same size as X, with replacement, from the original set X. To ensure a type-I error P e = α, the threshold must be set as the (1 − α) quantile of the empirical distribution of the test T(x * ).
Compute threshold . . Figure 1: General bootstrap detection procedure for a nominal P e = α and B bootstrap resamples.
The bootstrap-based methods are generally computationally very expensive, but in an area of exponentially declining computational costs, computer-intensive methods such as the bootstrap are becoming increasingly attractive. These methods have found wide applications in diverse fields such as engineering, astronomy, biology, genetics, economics, and finance [23][24][25][26][27]. It is not surprising that in the near future, the bootstrap will gain considerable attention. In a recently published text book [24], we found this comment "The revolution [bootstrap and jacknife procedures] has not yet become widespread . . .. Nevertheless, the revolution is happening and you should be prepared for its eventual results."

THE NHD TEST
Our NHD method compares the test statistic (7) to a threshold γ. If T(z i ) > γ, the ith cell is declared to be nonhomogeneous and eliminated. Thus the remaining bins are used to get a new estimate of R. The method is repeated to remove possible outliers in the secondary data set. At the end of the iterative process, the survival set is declared homogeneous and can be used by any adaptive radar detector.
We will present now the method to set the excision threshold γ. If the true or an asymptotic distribution of the test statistic is known, it is relatively simple to set the threshold according to the prescribed P e (analytically or numerically). But this is not the case here since the distribution of the test is unknown and an asymptotic one is not available because the number of secondary data is limited. In such situations, the bootstrap is well motivated [13]. At each iteration, the distribution of the test statistic is estimated using the bootstrap [25] and γ is set to be equal to the (1 − α) quantile, where α is the desired type-I error P e . The algorithm is as shown in Algorithm 1.
At the end of the algorithm, the survival set is decided to be homogeneous.

EURASIP Journal on Advances in Signal Processing
Step 0. Data collection: Get the matrix of secondary data Z = [ z 1 z 2 · · · z K ] from K range bins surrounding the CUT.
Step 1. Resampling: After centering the data, resample with replacement to obtain Step 2. Bootstrap test statistic: Compute (a) The sample bootstrap normalized covariance matrix using (10) as follows: (7) with R * replacing R. Step 3. Repetition: Repeat Steps 1 and 2, B times to obtain T * The number of iterations N it depends on the convergence criterion to be used. One can set a priori a small number of iterations. Another way to proceed is to stop the iterative process if after a given number of successive iterations, no outliers are detected. In the next section we will investigate the performance of these two methods in homogeneous clutter as well as in nonhomogeneous clutter.

PERFORMANCE ANALYSIS
In this section, we analyze the performance of our proposed algorithm in simulated K-distributed clutter, which is well motivated for sea clutter [28], as well as in real sea clutter collected by IPIX radar in February 1998. To the best of our knowledge, this is the first time the bootstrap is tested with real sea clutter.

Simulated clutter
We first consider the case of simulated clutter. The secondary data is generated according to the product model (1). The texture follows the pdf as given in (2). The covariance matrix of the speckle component is generated according to [R] i, j = ρ |i− j| . The one-lag correlation coefficient ρ is 0.9 (common value for sea clutter [17]). The steering vector is taken to be We consider the homogeneous as well as the nonhomogeneous secondary data (presence of interfering targets) case. At the first stage, we set the number of iterations N it = 1 for both the homogeneous and nonhomogeneous situation. Later on, we will present an automatic iterative version of the algorithm.

The homogeneous case
In all the following, when not stated otherwise, the number of bootstrap resamples is B = 100, [29]. Figure 2 shows the test statistic versus range in a homogeneous situation (no interfering targets in the secondary cells) for spiky clutter (ν = 0.5) with the following test parameters: μ = 1, N = 8, K = 32, N it = 1, and P e = 0.025. We see that the values of the test statistic do not exceed the threshold, confirming homogeneity of the data. Figure 3 shows also a homogeneous situation for a very spiky clutter (ν = 0.1) with the same test parameters as before. Even if the clutter is very spiky, the test does not exceed the threshold. Figure 4 depicts the case of large data samples with P e = 0.01. One notes that, for a large number of secondary data, it is necessary to choose low values of P e , because high values of P e (corresponding to lower thresholds) lead to erroneously discarded homogeneous data. This situation can be avoided in practice if the number of secondary range bins is not taken too large so as to avoid nonhomogeneities. Experimental tests have demonstrated that the number of homogeneous secondary data is often limited [30]. In a practical radar system, P e is a design parameter which is set by the operator according to a certain tolerable level of false rejecting homogenous data. It is desirable to minimize this error since we do not wish to edit out range bins that do not contain interfering targets. On the other hand, relatively high values of P e allow us to detect outliers with low power level. Therefore, the design of P e is always a compromise between all these constraints. Figure 5 plots the bootstrap and Monte Carlo cumulative density functions (CDF) of the test statistic for different values of N and K. One notes the negligible difference between the Monte Carlo CDF and the estimated CDF using the bootstrap.
To investigate the impact on P e , we set a prescribed value of P e0 = 0.025 and run 5000 Monte Carlo simulations on homogeneous data and estimate the actual probability of error P e . Tables 1, 2, 3, and 4 show the ratio P e /P e0 for different values of N, K, and ν. We have considered several values of the shape parameter ν, starting from low values corresponding to a spiky clutter to high values corresponding A. Younsi et al.   to a less spiky clutter. It is worth observing that in a real-live situation, the shape parameter varies from range cell to range cell [28,31,32] leading to secondary data with different measure of clutter spikiness. From the results, one can see that the error ratio is approximately constant for a given N to K ratio as a function of shape parameter.

The nonhomogeneous case
We now consider the nonhomogeneous case. We inject for example an interfering target in the range bin number K/2, and assess the capability of the algorithm to excise it (this   cell acts here as a CUT). We use the range bins surrounding this CUT (which are homogeneous in this case) to set the adaptive threshold according to the algorithm. Figure 6 depicts the result where it is apparent that the algorithm detects the nonhomogeneity.
The situation of Figure 6 is repeated 1000 times for N = 8 and different values of K. Table 5 shows the number of times the CUT exceeds the threshold. Similar results were also obtained with different values of N and K but not reported here.
Until now, we have considered that the range bins surrounding the cell containing the nonhomogeneity used to   To illustrate this phenomenon, we consider the same situation as in Table 5, except that we inject randomly in one or more range bins surrounding the CUT another interfering target (outliers). Let n i denote the number of injected outliers. Column 3 of Table 6 shows the performance of the NHD detector for different situations with N = 8, μ = 1, and ν = 0.5. The detection of the interfering target in the CUT is seriously affected by the contamination number n i , compared to Table 5.   To alleviate this problem, a heuristic solution can be considered by increasing the number of iterations N it . But the question is how to choose this number or equivalently, what is the convergence criterion? Here we adopt an automatic stopping rule as follows: if after two successive iterations no outliers are detected, the algorithm is stopped. To check the performance of this method, we consider again the same situation as before. The results obtained are shown in column  4 of Table 6. We see that the detection of the interfering target in the CUT is improved considerably.
To complete the performance analysis of this automatic NHD, we must look at its influence on the probability of error. As in Table 1, we run 5000 Monte Carlo simulations in homogeneous and nonhomogeneous clutter and estimate the actual probability of error. The ratio P e /P e0 is shown in Table 7.
In Figure 7, we plot the NHD test used in [2,3,7] which compares the GIP statistic z H i R −1 S z i to its theoretically mean value N. In Figure 8, we plot the NHD test proposed in [9], comparing the normalized GIP z H i R −1 S z i /K with the threshold-setting determined according to [9, equation (4.2)]. Here R S is simply the sample covariance matrix. It is evident that, for almost range bins, the GIP and the normalized GIP statistics exceed the threshold, leading to the declaration of nonhomogeneity, when in fact, the data is homogeneous. For the same data set, the proposed test statistic using the bootstrap and ACE decides for homogeneity of the set under test (see Figure 9). To confirm the results of these Figures, we conduct a Monte carlo analysis with the GIP and normalized GIP statistics. We found P e /P e0 = 15.5 for the GIP and P e /P e0 = 8.8 for the normalized GIP. We observe a drastic degradation of P e with the two methods allowing us to claim that the GIP and the normalized GIP statistics cannot be used in non-Gaussian clutter.

Real sea clutter data
This section is devoted to the performance assessment of the algorithm in the presence of real sea clutter collected by the X-band McMaster university, Canada, IPIX radar in February 1998 [33]. Table 8 shows the specifications of the data files which are exploited here. The data is first preprocessed in order to remove the DC offset and the phase imbalance due to hardware imperfections and then stored in an N t × N c complex matrix. A detailed statistical analysis   of adopted real data has been conducted in [31] for the file 19980204 − 220849 and [32] for the file 19980223 − 165836. The results have shown that the former data set follows a generalized K-distribution model and the later follows a Kdistribution.
The procedure used to assess the performances is now described. We consider an N × (K + 1) data window, where N is the number of pulses and K the number of secondary cells, the primary cell denoted by P c is set in the middle of the window. The data window is slid in space from range bin to range bin and in time from N samples to the next N samples until the end of the data set. The total number of To evaluate the actual probability of error P e , we perform the bootstrap-based nonhomogeneity detector for each data window and count the total number of times the test exceeds the threshold, say, N total . The actual P e is then P e = N total /N trials .

The homogeneous case
All the subsequent analyses assume that the steering vector where T is the pulse repetition time and f d the Doppler frequency which is chosen in order to coincide with the peak of the clutter power spectral density (PSD) of the considered data set. This is tantamount to considering the worst case of a target embedded in deep clutter. We set the prescribed probability of error to P e0 = 0.025. The results  of the experiment are illustrated in Table 9 for different polarization of the data set. The horizontally polarized data shows a little increase in the P e because of the spiky nature of the data on this polarization. The iterative algorithm using the automatic stopping rule shows also a little increase in the P e , this is due to the fact that in the homogeneous case, increasing the number of iterations N it , we can discard some cells even if they are homogeneous. This reduces the number of homogenous secondary data to be used in covariance estimation leading to an increased P e . One notes that P e is approximately constant and very close to the prescribed value for both the iterative and noniterative algorithm. Figure 10 depicts the test statistic versus range in homogeneous real sea clutter.

The Nonhomogeneous case
We now inject an interfering target p i = α i s with an interference to noise ratio defined as INR = α 2 i /E[z 2 ] in the primary cell P c . Figure 11 shows the test statistic versus range where it is clear that the test in the primary cell exceeds the threshold. In Figure 12, we plot the probability of correctly rejecting (P r ) this outlier, versus INR, following the same procedure outlined before to estimate P e .
The situation of multiple outliers in the secondary data is depicted in Figure 13, in addition to the interfering target in P c , another interfering one with the same INR A. Younsi et al.  is injected randomly in the other cells. The figure shows the performance of the algorithm using only one iteration (N it = 1) and the algorithm using the automatic stopping rule. One notes that with only one iteration the probability of correctly rejecting the primary interfering target is highly degraded (masking effect) while the iterative algorithm improves considerably P r . Figure 14 shows another situation where two interfering targets are present in the secondary cells, here also the automatic algorithm gives the better P r .
In all the previous analysis, we have neglected the effect of the receiver noise. This is often a realistic assumption since in general, for radar systems with reliable detection performance, the clutter power is much greater than the power of the thermal noise. As an example, we have estimated the  clutter-to-thermal noise ratio (CNR) of the real data used here and we have found CNR 6 dB and CNR 20 dB for the file 19980204 − 220849 and the file 19980223 − 165836, respectively. To investigate the effect of thermal noise on the performance of our algorithm, we add an artificial noise to the real data and estimate the actual probability of error. The additional thermal noise is generated according to the complex normal model CN(0, σ 2 I). Table 10 shows the ratio P e /P e0 versus the clutter-to-noise ratio CNR = μ/σ 2 . We can conclude that for CNR ≥ −5 dB the algorithm maintains the probability of error approximately constant. This means that the proposed algorithm is robust even in the presence of high-level thermal noise. The degradation of P e grows with the power of the additional noise.  The impact of thermal noise on the probability of correctly reject interfering targets is depicted in Figure 15 for CNR = 0 dB and CNR = −5 dB. For hight INRs, the effect of thermal noise is insignificant. Elsewhere, the degradation of P r is acceptable for CNR = 0 dB since the loss induced on the INR is less than 2 dB but not for the case of CNR = −5 dB.
Finally, an other important point to consider is the implementation issue. As all computer-intensive methods, the bootstrap-based methods suffer from computational complexity. The algorithm proposed here suffers also from this problem since it requires an inversion of a data matrix many times. It involves O(BN it KN 2 ) floating point operations, where we used the usual Landau notation O(n) which means that an algorithm is O(n) if its implementation requires a number of floating points operations (flops) proportional to n. However, for some kinds of radar systems, the parameters N and K are not very large, leading to a lowdimension matrix to be inverted. N is the number of hits during the time on target and it is dependant on the radar parameters such as the scan rate, the bandwidth and the PRF. As concerning the number of training data, in [34] a number K = 5N was considered and in [35] a rule of thumb K ≥ 2N was proposed. Also, some simplifications can be obtained in practice by considering, for example, application of parallel processing techniques since the supercomputing industry is undergoing a significant resurgence in the development of these techniques and in the development of processors with very high speed. In addition to that, a suitable choice of the matrix inversion algorithm (QR factorization techniques are often employed in the matrix inversion) can significantly reduce the computational burden.

CONCLUSION
We have presented a data adaptive thresholding algorithm that can be used as a preprocessing stage for outlier removal from secondary data, in an adaptive radar which uses an estimate of the interference covariance matrix, from range bins surrounding the cell under test. The noise is supposed to follow the compound-Gaussian model. The proposed algorithm does not require the knowledge of the pdf of the spiky component of the clutter. It automatically adapts the threshold according to the data it collects by means of bootstrap estimation of the pdf of the test statistic. The algorithm is tested with very spiky K-distributed simulated clutter as well as with real sea clutter. In both cases, the algorithm gives good results. For a critical situation with multiple outliers in secondary data, an automatic iterative version of the algorithm is suggested. The robustness of the algorithm in the presence of thermal noise was also considered and we found that the algorithm is robust even in the presence of high-level thermal noise.