Radar Detection of Fluctuating Targets under Heavy-Tailed Clutter Using Track-Before-Detect

This paper considers the detection of fluctuating targets in heavy-tailed clutter through the use of dynamic programming based on track-before-detect (DP–TBD) in radar systems. The clutter is modeled in terms of K-distribution, which can be widely used to describe non-Gaussian clutter received from high-resolution radars and radars working at small grazing angles. Swerling type 1 is considered to describe the target fluctuation between scans. Conventional TBD techniques suffer from significant performance loss in heavy-tailed environments due to the more frequent occurrences of target-like outliers. In this paper, we resort to a DP–TBD algorithm based on prior information, which can enhance the detection performance by using the environment and target fluctuating information during the integration process of TBD. Under non-Gaussian background, the expressions of the likelihood ratio merit function for Swerling type 1 targets are derived first. However, the closed-form of the merit function is difficult to obtain. In order to reduce the complexity of evaluating the merit function and the computational load, an efficient approximation method as well as a two-stage detection approach is proposed and used in the integration process. Finally, several numerical simulations of the new strategy and the comparisons are presented to verify that the proposed algorithm can improve the detection performance, especially for fluctuating targets in heavy-tailed clutter.


Introduction
The detection of fluctuating targets with low signal-to-clutter ratio (SCR) is of significant importance in radar systems. Conventional detecting and tracking algorithms use thresholded detection as input. A target with low signal-to-clutter ratio is often lost due to information being irreversibly discarded after thresholding. Multi-frame integration is an effective strategy used in radar applications to detect dim targets by integrating signal returns over multiple consecutive scans. In the presence of a moving target, multi-frame integration requires track-before-detect (TBD) techniques to correctly correlate data over time.
Dynamic programming based on TBD (DP-TBD) is one of the TBD techniques [1,2], that has attracted extensive attention for the advantages of simplicity and needing less information. It transforms the integration into an optimal estimation of the physically admissible trajectory by maximal integration value of the merit function, which is a kind of multi-frame test statistic. DP-TBD can detect a target of arbitrary motion form and has been widely applied to several kinds of sensors [3,4]. In order to solve the problem of high-dimensional maximization under a multi-target environment, a novel partition method to cluster targets into well separated groups was proposed 1.
In order to limit complexity while still retaining the benefits of DP-TBD, we resort to a two-stage detection process with different resolution cells.

2.
For typical non-Gaussian distributed clutter (K-distribution) and a typical target amplitude fluctuation model (Swerling 1), the DP-TBD algorithm based on prior information is proposed. By using the likelihood ratio merit function in DP integration, the performance loss produced by the "heavy-tailed" clutter measurements can be reduced.

3.
An efficient but accurate approximation method is proposed to reduce the complexity of evaluating the merit function.
The remainder of this paper is organized as follows: Section 2 presents the notations and system models. In Section 3, a two-stage detection approach is proposed at first, and the expressions of the likelihood ratio merit function are derived in K-distributed clutter background for Swerling target of type 1; the implementation issues of the merit function are also discussed. Simulation results are showed by comparing different DP-TBD strategies in Sections 4 and 5 provides some conclusions.
Mathematical notations used in this paper are described as follows. s n is the target kinematic state at scan n; A n denotes an amplitude from the Swerling 1 target and c n (i, j) denotes the K-distributed clutter; a n is the measurement amplitude in K-distributed clutter background. I(s n ) is defined as the Sensors 2018, 18, 2241 3 of 13 merit function at scan n; V(s n−1 ) is defined as the maximal integration value of all the admissible trajectories; τ(s n ) is the collection of states at scan n for which transition to s n is possible; Ψ(s n ) is the retracing function, indicating the best state of the previous scan.

Kinematic Model
As shown in Figure 1, we assume that there is only one target in the surveillance region, whose kinematic state at scan n is denoted by the vector s n . The kinematic vector s n is specified by: where denotes matrix transpose, r n and θ n denote the range and azimuth measurement, respectively, R 2 denotes the two-dimensional state space, and N denotes the number of consecutive frames processed in a DP-TBD integration batch. The evolution of the target state is modeled by the linear process as: The term w n is the process noise, F is the transition matrix. Every real target must comply with some physical constraints on its kinematics, such as the maximum target velocity considered in this paper. The radial and tangential velocity can be calculated by two successive scans, which are given by: where T denotes the time interval between successive scans.

Kinematic Model
As shown in Figure 1, we assume that there is only one target in the surveillance region, whose kinematic state at scan n is denoted by the vector n s . The kinematic vector n s is specified by: where ' denotes matrix transpose, n r and n  denote the range and azimuth measurement, respectively, 2  denotes the two-dimensional state space, and N denotes the number of consecutive frames processed in a DP-TBD integration batch. The evolution of the target state is modeled by the linear process as: The term n w is the process noise, F is the transition matrix.
Every real target must comply with some physical constraints on its kinematics, such as the maximum target velocity considered in this paper. The radial and tangential velocity can be calculated by two successive scans, which are given by: sin n n n n n radial n n n n gential where T denotes the time interval between successive scans.

Measurement-Based Model
The measurement data consists of , at scan n can be expressed as [20]:

Measurement-Based Model
The measurement data consists of M r cells in the range-dimension and M θ cells in the azimuth-dimension. If no target exists (hypothesis H 0 ), the (i, j)th recorded resolution cell, z n (i, j) 1 ≤ i ≤ M r , 1 ≤ j ≤ M θ , at scan n can be expressed as [20]: while in the presence of a target (hypothesis H 1 ), the recorded resolution cell z n (i, j) can be expressed as: where A n denotes a complex fluctuated amplitude measurement from the target, and c n (i, j) denotes the K-distributed clutter, which is assumed in this paper. A Swerling 1 fluctuation model supposes that returned signal power per pulse is constant during a single scan, but fluctuates independently from scan to scan. The probability density function (PDF) of the Swerling 1 target amplitude A n is given by: with σ being the mean squared target amplitude.

K-Distributed Clutter Model
The K-distributed model is proposed as a model for radar clutter in this paper, which has the probability density function as: In formula, Γ(·) denotes the Gamma function and K α−1 (·) denotes the modified Bessel function of the second kind, a n is the measurement amplitude, β is scale parameter which describes the intensity of the clutter, and α is the shape parameter which determines the shape of the distribution function. For α → ∞ , the K-distribution turns into the Rayleigh distribution. PDFs of K-and Rayleighdistributions are shown in Figure 2a while the K-distributed clutter (α = 2, β = 1) is shown in Figure 2b.
where n A denotes a complex fluctuated amplitude measurement from the target, and   , n c i j denotes the K-distributed clutter, which is assumed in this paper. A Swerling 1 fluctuation model supposes that returned signal power per pulse is constant during a single scan, but fluctuates independently from scan to scan. The probability density function (PDF) of the Swerling 1 target amplitude n A is given by: with  being the mean squared target amplitude.

K-Distributed Clutter Model
The K-distributed model is proposed as a model for radar clutter in this paper, which has the probability density function as: In formula,     denotes the Gamma function and   1 K    denotes the modified Bessel function of the second kind, n a is the measurement amplitude,  is scale parameter which describes the intensity of the clutter, and  is the shape parameter which determines the shape of the distribution function. For    , the K-distribution turns into the Rayleigh distribution. PDFs of K-and Rayleigh-distributions are shown in Figure 2a   Meanwhile, K-distribution can also be viewed as a Rayleigh distribution modulated by a Gamma distribution for convenience:  Meanwhile, K-distribution can also be viewed as a Rayleigh distribution modulated by a Gamma distribution for convenience: p(a n ) = +∞ 0 p(a n |η )p(η)dη (8)  where p(a n |η ) = 2a n η exp − a n 2 η (9)

Development of the Proposed Strategies
The DP-TBD algorithm decomposes the integration among N successive scans into N sub-processes. The nth sub-process contains all the measurements up to scan n. The target can be detected and tracked by calculating the maximum of the energy integration value through a recursive model, which could be expressed as: where I(s n ) is defined as the merit function at scan n; V(s n−1 ) is defined as the maximal integration value of all the admissible trajectories; τ(s n ) is a collection of states at scan n for which a transition to s n is possible, and it can be obtained by the location and maximum velocity of the target; Ψ(s n ) is the retracing function, indicating the best state of the previous scan, which makes the integration value reach its maximum.
In summary, DP-TBD implements the equivalent of an exhaustive search in an efficient manner by enumerating and valuing all physical admissible state sequences, finally returning the state sequences whose final maximal integration value V(s N ) exceeds a given detection threshold γ, i.e.,: There are mainly two problems throughout the process. Firstly, the computational complexity of DP-TBD is unaffordable in the presence of a high-mobility target when the number of resolution elements is large. The discretization of state space is always based on the sensor's resolution so as to make full use of the measurements and achieve possibly accurate estimates. In this situation, strategies hardly lead to real-time implementable schemes, even resorting to a dynamic programming algorithm. In order to reduce the burden of computation, a two-stage detection approach is proposed in this work. Secondly, most of the previous work on DP-TBD assumed that the background model would be Rayleigh or Gaussian distribution with a known power. Such assumptions may not be adequate, as in the real world a more heavy-tailed background model than expected is often encountered. To improve the detection performance, we propose a novel DP-TBD algorithm based on the prior information to solve the aforementioned problem. In this paper, the merit function is set to be the likelihood ratio under both target-present hypothesis and null-target hypothesis in a surveillance region which is characterized by K-distributed background, and the simulated data would be tested for presenting the performance.

Two-Stage Detection Approach
Generally, DP-TBD is a grid-based method that estimates target trajectory by means of searching all the admissible paths in a discrete state space and the discretization of state space is based on the sensor's resolution. It is inefficient and computationally costly to carry out an accurate search (i.e., the search grid is exactly divided based on the sensor's resolution) since only a fraction of measurements are related to the actual target when the surveillance region is large. In order to reduce the computational load, while still retaining the benefits of TBD, here we resort to a two-stage detection approach which is illustrated in Figure 3. searching all the admissible paths in a discrete state space and the discretization of state space is based on the sensor's resolution. It is inefficient and computationally costly to carry out an accurate search (i.e., the search grid is exactly divided based on the sensor's resolution) since only a fraction of measurements are related to the actual target when the surveillance region is large. In order to reduce the computational load, while still retaining the benefits of TBD, here we resort to a two-stage detection approach which is illustrated in Figure 3. At stage 1, we first obtain the raw data at scan n, and roughly calculate the measurement n z  under the condition of low grid resolution. The target states are estimated by searching discrete grids with larger cell size based on the DP integration. After N times loop, the maximum of the At stage 1, we first obtain the raw data at scan n, and roughly calculate the measurement z n under the condition of low grid resolution. The target states are estimated by searching discrete grids with larger cell size based on the DP integration. After N times loop, the maximum of the energy integration value V(s N ) at scan N could be obtained by the process. For a single target model, the maximum integration value V(s N ) which exceeds detection threshold γ 1 is used to determine the existence of the target. If there is a target presented in the surveillance region, we could refine the target trajectory in stage 2.
In order to obtain a more accurate estimate, stage 2 is employed to recalculate the measurements under the high grid resolution condition. Once the maximum integration value V(s N ) exceeds the So the recovered trajectory estimate isŜ N = {ŝ 1 , . . . ,ŝ N }. The algorithmic description of the proposed two-stage TBD approach is shown in Table 1. Table 1. Algorithmic description of the proposed two-stage TBD.

Stage 1
Mearsurement: The surveillance region is divided into M r × M θ grid cells based on the resolution of the radar system, i.e., ∆r and ∆θ, where M r and M θ denote the number of cells in range and azimuth, respectively. To realize the target search with larger cell size, the state space is re-discretized by ∆r and ∆θ to obtain M r × M θ grid cells at first. As shown in Figure 4, all the measurements and DP integrations are processed in stage 1 based on the new state space, which may obtain a rough target trajectory by less computation. Then in stage 2, DP integration concentrates on the part of states which are indicated by stage 1. As calculations of less meaningful states could be avoided, the computational costs will become more reasonable. model, the maximum integration value   N V s which exceeds detection threshold 1  is used to determine the existence of the target. If there is a target presented in the surveillance region, we could refine the target trajectory in stage 2.
In order to obtain a more accurate estimate, stage 2 is employed to recalculate the measurements under the high grid resolution condition. Once the maximum integration value The surveillance region is divided into

Derivation and Implementation of the Merit Function
Combined (6), PDF for the Swerling 1 target in K-distributed clutter is given by: p(a n |η, s n ) = 2a n η + σ exp − a n 2 η + σ and p(a n |s n ) can be derived by marginalizing over η since η is random, i.e., p(a n |s n ) = ∞ 0 p(a n |η, s n )p(η)dη = ∞ 0 2a n η+σ exp − a n 2 η+σ η α−1 (16) where the integrand f (η) is given by: Substituting (7) and (16) into the expression of merit function I(s n ) at scan n, I(s n ) can be written as: I(s n ) = ln p(a n |s n ) p(a n ) = ln a n −α+1 Although the integrand f (η) in (17) has no closed-form solution, it can be evaluated with reasonable accuracy by using the trapezoidal rule, i.e., where f (η i ) is sample point drawn from the time interval δ i , δ i is a sampling interval which is short enough to cover the effective support of f (η), and N sa denotes the number of sample points. The sample points can be obtained by either deterministic sampling with a uniform grid or stochastic importance sampling. Since the integrand f (η) may tend quickly towards ∞ when η → 0 , while tending slowly towards 0 when η → ∞ . A reasonable approximation obtained by deterministic uniform grid sampling or stochastic importance sampling is difficult to carry out. A grid with a variable resolution method was proposed in [21] to approximate the merit function, which also leads to high computational complexity.
In order to reduce the complexity of approximation, we could possibly circumvent these problems by generating a lookup table offline with sample points using a uniform grid. The number of sample points with uniform grid is large enough to approximate the integrand f (η) accurately. Based on the lookup table, this calculating method trades little cost of precision and memory space for a great improvement on running speed in the calculation. The histograms of generation data and theory PDF are shown in Figure 5 for the Swerling targets type 1 with different parameters. According to Figure 5, we conclude that the approximation error is negligible.
Note that the K-distribution shape parameter α and the scale parameter β are supposed to be known in the derivation of merit function. In the case where the background is significantly heavy-tailed and the parameters are unknown, we should estimate the parameters first, which can be obtained through a numerical maximization of the likelihood function. Since the maximum likelihood techniques require numerical optimization routines and evaluation of Bessel functions, they are computationally intensive and, therefore, inappropriate for evaluation of large data sets. Abraham [22] recommended moment estimators based on the first and second moments, which can be used as our estimator in this work.
heavy-tailed and the parameters are unknown, we should estimate the parameters first, which can be obtained through a numerical maximization of the likelihood function. Since the maximum likelihood techniques require numerical optimization routines and evaluation of Bessel functions, they are computationally intensive and, therefore, inappropriate for evaluation of large data sets. Abraham [22] recommended moment estimators based on the first and second moments, which can be used as our estimator in this work.

Simulation
In this paper, the detection performances of conventional DP-TBD and the proposed strategy for a Swerling type 1 target are assessed. We assume that the measurement noise satisfies K-distribution, and each measurement frame consists of

Performance Analysis
We assess the performances of different strategies via the probability of track detection d P , which is a performance metric for both detection and tracking performance. d P is defined as the probability of the maximum integration value exceeding the detection threshold, and its final position is within a certain range of the actual target position. In addition, the root-mean square error (RMSE) on the estimation of the target position is also considered, which is defined as:

Simulation
In this paper, the detection performances of conventional DP-TBD and the proposed strategy for a Swerling type 1 target are assessed. We assume that the measurement noise satisfies K-distribution, and each measurement frame consists of M r × M θ = 180 × 90 resolution cells. The number of frames processed in a DP batch is N = 6, while the number of possible state transitions in a scan is q = 9. This scenario is run 1000 times for various SCR and shape parameters while the false alarm is fixed as P FA = 10 −3 .

Performance Analysis
We assess the performances of different strategies via the probability of track detection P d , which is a performance metric for both detection and tracking performance. P d is defined as the probability of the maximum integration value exceeding the detection threshold, and its final position is within a certain range of the actual target position. In addition, the root-mean square error (RMSE) on the estimation of the target position is also considered, which is defined as: where H 1 is the event that the target is confirmed, and e 2 (s n ) is the Euclidean distance between the true and estimated target position. Performance and RMSE comparison of the conventional DP method and the proposed method based on prior information is shown in Figure 6. For K-distributed clutter and a Swerling 1 target, the proposed method performs better than the traditional integration method. It can also be concluded that the proposed method, which is processed with only one stage (blue solid line) or two stages (red solid line), could achieve almost identical performance while the latter obtains further computational reduction. based on prior information is shown in Figure 6. For K-distributed clutter and a Swerling 1 target, the proposed method performs better than the traditional integration method. It can also be concluded that the proposed method, which is processed with only one stage (blue solid line) or two stages (red solid line), could achieve almost identical performance while the latter obtains further computational reduction. For different parameters of K-distributed clutter, detection performances are shown in Figure 7. With an increase of shape parameter α, both conventional DP and the proposed method on prior information achieve significant performance improvement. That is because when α is increasing, the K-distributed clutter is smoother and the frequency of target-like outliers is lower. Note that when α = 50, since K distribution almost degenerates to Rayleigh distribution in this case, the detection performances are nearly identical. For different parameters of K-distributed clutter, detection performances are shown in Figure 7. With an increase of shape parameter α, both conventional DP and the proposed method on prior information achieve significant performance improvement. That is because when α is increasing, the K-distributed clutter is smoother and the frequency of target-like outliers is lower. Note that when α = 50, since K distribution almost degenerates to Rayleigh distribution in this case, the detection performances are nearly identical. based on prior information is shown in Figure 6. For K-distributed clutter and a Swerling 1 target, the proposed method performs better than the traditional integration method. It can also be concluded that the proposed method, which is processed with only one stage (blue solid line) or two stages (red solid line), could achieve almost identical performance while the latter obtains further computational reduction. For different parameters of K-distributed clutter, detection performances are shown in Figure 7. With an increase of shape parameter α, both conventional DP and the proposed method on prior information achieve significant performance improvement. That is because when α is increasing, the K-distributed clutter is smoother and the frequency of target-like outliers is lower. Note that when α = 50, since K distribution almost degenerates to Rayleigh distribution in this case, the detection performances are nearly identical.

Computational Complexity Analysis
The complexity of the conventional DP-TBD method is are the number of range and azimuth resolution elements, q is the number of possible state transitions in a scan, and N is the number of the integration scans. In comparison with the conventional method, the two-stage detection approach schemed in Figure 3 has low computational

Computational Complexity Analysis
The complexity of the conventional DP-TBD method is O(M r M θ qN), where M r and M θ are the number of range and azimuth resolution elements, q is the number of possible state transitions in a scan, and N is the number of the integration scans. In comparison with the conventional method, the two-stage detection approach schemed in Figure 3 has low computational complexity.
For a different number of frames and different number of state transitions, detection performances are shown in Figure 8. With the increase of frames and state transitions, the proposed method achieves a performance improvement, but not very much. However, the computational cost increases rapidly.

Computational Complexity Analysis
The complexity of the conventional DP-TBD method is are the number of range and azimuth resolution elements, q is the number of possible state transitions in a scan, and N is the number of the integration scans. In comparison with the conventional method, the two-stage detection approach schemed in Figure 3 has low computational complexity.
For a different number of frames and different number of state transitions, detection performances are shown in Figure 8. With the increase of frames and state transitions, the proposed method achieves a performance improvement, but not very much. However, the computational cost increases rapidly. to realize the target search with larger cell size. Since the DP integration in stage 2 is concentrated on the part of states, which are indicated by stage 1, the computational cost is small enough not to care.
The computational cost of the strategies is listed in Table 2 for different parameters. It can be seen that the computational cost depends on the number of a possible state transition q and the resolution elements. For the same resolution cells, the scenario at q = 9 costs almost three times as long as the scenario q = 4. Meanwhile, the computational cost reduces rapidly as the number of The computational cost in stage 1 is O M r M θ qN , where M r × M θ denotes grid cells re-discretized by ∆r and ∆θ to realize the target search with larger cell size. Since the DP integration in stage 2 is concentrated on the part of states, which are indicated by stage 1, the computational cost is small enough not to care.
The computational cost of the strategies is listed in Table 2 for different parameters. It can be seen that the computational cost depends on the number of a possible state transition q and the resolution elements. For the same resolution cells, the scenario at q = 9 costs almost three times as long as the scenario q = 4. Meanwhile, the computational cost reduces rapidly as the number of resolution elements decreases. For example, when q = 4, the CPU times for M r × M θ = 180 × 90, M r × M θ = 90 × 45 and M r × M θ = 60 × 30 are 308 ms, 224 ms and 146 ms, respectively.

Conclusions
This paper has presented the systematic treatment of heavy-tailed clutter from a target detection and tracking perspective. Target detection in K-distributed clutter is more challenging than in Gaussianor Rayleigh-distributed clutter due to the higher likelihood of target-like outliers, especially for a fluctuating target. In this work, we dealt with the fluctuating target detection and tracking problem using a modified DP-TBD method. The contributions are as follows: first we have solved the target detection problem using two-stage detection architecture to avoid calculations of less meaningful states. Secondly, for a Swerling 1 target in a K-distributed background, the merit function was derived and implemented in the integration process of DP-TBD to enhance radar detection performance. In order to reduce the complexity of integral calculation, we also resorted to the trapezoidal rule with a generating lookup table.
Numerical analysis demonstrated that performance improvement could be applied via the proposed DP-TBD algorithm based on prior information, especially for heavy-tailed K-distributed clutter. Moreover, simulation results suggested that a trade-off between performance and computational complexity exists. Further research may investigate the performance of the proposed DP-TBD method experimentally. It may also be of interest to investigate other background models than the K-distribution.
Author Contributions: J.G. contributed to the original ideas and designed the simulations; J.G. and W.W. performed the simulations and analyzed the data; J.D. provided significant guidance in research planning and revised the manuscript. All authors contributed to and approved the written manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.