Abstract
A particle filter based track-before-detect (PF-TBD) algorithm is proposed for the monopulse high pulse repetition frequency (PRF) pulse Doppler radar. The actual measurement model is adopted, in which the range is highly ambiguous and the sum and difference channels exist in parallel. A quantization method is used to approximate the point spread function to reduce the computation load. The detection decisions of the PF-TBD are fed to a binary integrator to further improve the detection performance. Simulation results show that the proposed algorithm can detect and track the low SNR target efficiently. The detection performance is improved significantly for both the single frame and the multiframe detection compared with the classical detector. A performance comparison with the PF-TBD using sum channel only is also supplied.
1. Introduction
The developments of stealthy military aircraft and cruise missiles recently have emphasized the need for detection and tracking of low signal-to-noise ratio (SNR) targets. This need is especially urgent for a radar seeker because of its limited battery capacity and antenna size. High pulse repetition frequency (PRF) pulse Doppler is generally used in a radar seeker at early detection stage, which allows thermal noise-limited detection of targets with high radial velocities [1]. Noncoherent or binary integration is often used after the coherent processing to improve the detection performance. But the radar data rate and the unknown target motion have limited the coherent processing interval (CPI) and noncoherent/binary times. The azimuth and elevation are measured by monopulse generally, which is a widely used technique to provide accurate angle measurements in the tracking radar. A monopulse system for estimating one angle typically consists of two identical antennas, either separated by some distance (phase monopulse) or at the same phase center but with a squint angle (amplitude monopulse), whose outputs are summed up to produce a sum channel and are subtracted to yield the difference channel as shown in Figure 1. The angular information is contained in the monopulse ratio providing the function is reversible. Poor monopulse estimation performance under low SNR has also deteriorated the guidance performance.
Track-before-detect (TBD) is a simultaneous detection and tracking paradigm that uses unthresholded data or thresholded data with significantly lower thresholds than those used in conventional detectors and integrates them over time according to the target dynamic model to improve the sensitivity to low SNR targets. Typical TBD is implemented as a batch algorithm using the Hough transform [2] or dynamic programming [3]. The Hough transform TBD is suitable only for linear trajectories. The dynamic programming TBD is studied more for the radar application and is applied in pulse Doppler radars in [4–6]. Particle filter based TBD (PF-TBD) was introduced by [7] and extended by [8–10]. Compared to the typical methods, it is recursive and does not require discretization of the state space.
For simplicity, most researches on PF-TBD are based on grayscale-image-like measurements (e.g., [8, 10]). Boers and Driessen [9, 11] have studied PF-TBD on search radar measurements. A Rao-Blackwellised PF-TBD is proposed for over-the-horizon radar in [12]. Multisensor PF-TBD is studied for MIMO radar [13]. There is no open literature addressing PF-TBD on monopulse radar to the best of our knowledge. In monopulse radar systems, the sum and difference channels exist in parallel as Figure 1 has shown. A PF-TBD algorithm similar to [12] can be applied by using only the output of the sum channel as the measurements. The target Doppler and intensity are estimated by it and then the bearing and azimuth are estimated by classical monopulse methods (e.g., ML method proposed by [14]). But from Figure 1 we can see that amplitude of the difference channel is comparable to that of the sum channel when the target is not at the beam center, which often occurred in the target searching stage. So fusion of the sum and difference channels using Bayesian theory in the PF-TBD algorithm is possible to improve the detection performance as well as the monopulse estimation performance.
In this work, the target and measurement models of the monopulse high PRF pulse Doppler radar are constructed. Based on them, we derive a PF-TBD algorithm which can effectively detect and track the low SNR target. Its detection performance is compared with the classical detector, which shows that more than 7 dB gain in SNR can be attained. A quantization method of approximating the point spread function is proposed to reduce the computation load of the PF-TBD. Binary integration of the PF-TBD’s detection result is proposed to further improve the detection performance, which is shown to be very effective and not limited by the target maneuver.
The rest of the work is organized as follows. In Section 2 the target and sensor models are formulated. The recursive Bayesian TBD filter for this application is described in Section 3 and its PF implementation procedure is derived in Section 4. Two simulated examples are presented in Section 5, in which the detection and estimation performances of the proposed algorithm are evaluated in comparison with the classical method and the sum-only PF-TBD. Conclusions and future work are drawn in the last section.
2. Target and Measurement Models
2.1. Target Model
The high PRF can measure Doppler unambiguously, but it is highly ambiguous in range, which precludes the pulse delay ranging. The range information is not a must for a radar seeker, however, since the proportional navigation is commonly adopted. As a result, only the target Doppler is involved in the target state vector in this paper. The target azimuth and elevation are measured by monopulse. For the sake of brevity, only one difference channel (azimuth difference or elevation difference) is considered. Moreover, the unknown target echo amplitude is also incorporated to implement the PF-TBD algorithm. The target state vector is then defined as where , , and denote the Doppler frequency, echo amplitude of the sum channel, and monopulse ratio of the target in frame , respectively. The Doppler frequency , where is the wavelength and is the radial velocity.
Although the dynamic model can be as general as for a particle implementation, where is the process noise sequence, for simplicity we model the target motion relative to the radar as the nearly constant velocity model with a white acceleration noise . The target echo amplitude and monopulse ratio are modeled as random walk processes with process noises and , respectively. The process noises , , and are mutually independent, zero mean white noise with variances , , and , respectively. Thus, the system dynamic equation is where is the CPI and . This target model accommodates not only target maneuver but also fluctuations of the target intensity and the monopulse ratio.
Target existence variable is modeled as a two-state Markov chain and . Here 0 denotes the event that the target is absent, while 1 denotes the opposite [15]. Furthermore, we define the transitional probabilities of target “birth” () and “death” () as Thus the transitional probability matrix is given by
2.2. Measurement Model
We assume that the target is located in the clutter-free region; thus the clutter is not considered in the signal model. When the target is present, the received signal sequences at the video stage of the sum and difference channels in frame are denoted as and and given by respectively, where is the amplitude of the difference channel, is some arbitrary phase, is the pulse repetition interval (PRI), and is index of the sample in an CPI. The background thermal noises and are mutually independent, zero mean, and temporally white complex Gaussian processes with the same variance. The Doppler frequency is assumed to be constant within an CPI.
The coherent integrations of the sum and difference echoes are done via fast Fourier transform (FFT) independently. To reduce peak side-lobe levels, the signal sequences are windowed before the FFT. The result of the coherent integration is given by where the subscript denotes sum channel or difference channel for simplification, is the next power of two that is greater than or equal to , for (also known as zero padding), is the windowing function, and is the index of the frequency bin. The signal’s unknown phase component is useless, so the magnitude of the spectrum in each frequency bin forms the set of measurements in frame . Then the measurement can be modeled as where is the complex modulus, , and and are the background noises of the sum and difference channels, respectively, after the coherent integration. Because of linearity of the FFT, and are also zero mean i.i.d. complex Gaussian noise processes. We assume that they both have a variance .
As has been stated that not all the frequency bins of the FFT result are of interest, only bins in clutter-free region constitute the set of measurements at frame ; that is, , where , is the horizontal velocity of the missile, and is the Doppler bin size.
Following the model described above, the likelihood in each frequency bin when the target is present has a Ricean distribution where is the modified Bessel function of order zero. The likelihood when the target is absent has a Rayleigh distribution
Because of the windowing before the FFT, the target (if present) power will spread into the bins in the vicinity of its location. Let denote the bins affected by the target (i.e., the target’s effect on the other bins is negligible); then the likelihood function of the whole measurement set when the target is present can be approximated as follows: and the likelihood function when the target is absent is
We denote the set of complete measurements up to frame as .
It is computational complex to calculate the for bins in in real time applications. The contribution of to bin in (i.e., point spread function) is generally approximated by a Gaussian-like function (e.g., [7, 8] for optical sensor). Using the Gaussian approximation method, the point spread function is where is the coherent integration gain and is a parameter to be designed to better approximate the amount of blurring introduced by the FFT windowing functions. But this approximation is valid only within a limited range as Figure 2 has shown. To solve this problem, we present an approximation approach which is calculation-free and more precise. Note that can be expressed as a function with a parameter and a variable . Because the windowing function can be taken as known a priori, we can quantize into a number of points (e.g., , for , where is the number of points each bin is quantized into, and we can store them as a look-up table in the read-only memory (ROM). In real time operations, the value of the quantized point nearest to the true point is read from the ROM and used; that is, where is the floor function and rounds to the nearest integer. The result of this approximation is also presented in Figure 2.
3. Recursive Bayesian Filtering Procedure
The posterior probability of target existence and are estimated recursively by a Bayesian method as follows. Given the joint posterior PDF at frame , and the latest measurement , the goal is to construct the joint posterior PDF at frame , . and are then estimated using .
Prediction. Prediction of is given by If , is undefined and no prediction of it is needed. If , the prediction step of can be expressed as where The transitional density is defined by the target dynamic model (2). The PDF denotes the initial target density on its appearance.
Update. The update equation using Bayes’ rule is given by where the prediction density is given by (17), the normalizing constant in the denominator is , and the likelihood function is where the likelihood function is given by (12).
Estimate. is estimated by taking marginal of as follows: Using expected a posterior (EAP) estimator, is estimated by
4. Particle Filter Implementation
To implement the recursive Bayesian filtering procedure, a SIR particle filter based TBD algorithm described in [8] is adopted with some modifications. As the particle filter tends to suffer from a progressive degeneration as the sequence evolves, an MCMC step referred to as resample-move in [16] is employed after importance resampling, which adds diversity to the particles without altering the underlying distribution [10]. A Metropolis-Hasting resample-move method is used as described in [10, 17]. Taking move of the , for example, a proposal distribution is defined, from which a sample is drawn for each particle after resampling. A monopulse ratio is obtained conditioned on the old monopulse ratio while keeping the other two states unchanged. Under the assumption that the proposal is symmetric, , the new particle is accepted or rejected on a test, formed by a ratio of likelihoods If , then the new particle, with monopulse ratio , is kept. Otherwise the new particle is kept in preference to the old particle only if , where is a uniform random number between 0 and 1. The move operation is used twice in this application, firstly to the amplitude and then to the monopulse ratio . Truncated Gaussian distributions with different variances and means at and , respectively, are used as the proposal distributions.
A detailed description of the TBD algorithm is given as follows.
Initialization. Set and generate samples from . If , generate from the birth density , or else, is undefined.
Then, given at each frame , go from Steps 1 to 5.
Step 1 (prediction). Generate on the basis of and . If , is undefined. If and , predict according to (2). For the new born particles, that is, those with and , generate from the birth density .
Step 2 (update). In the SIR filter, the prior PDF is chosen to be the important density and, thus, unnormalized weights are proportional to the likelihood functions. Consequently, using the likelihood ratios as unnormalized weights will have no effect on the performance of the SIR filter. Thus the importance weights are calculated by the following [7]: We simplify the likelihood as follows: From (10) and (11), can be simplified as where is calculated by . Then get the normalized weights by .
Step 3 (resample). Generate a new set of samples from and replace them using systematic resampling algorithm [18]. The weights of the new samples are not required since they are all equal to .
Step 4 (MCMC move). Generate a new set of samples from and replace them by move of using the Metropolis-Hasting method described above; do this again by move of . Note that this operation only changes the particles with .
Step 5 (state estimation). Estimate the posterior probability of target existence by If exceeds a certain threshold , target presence is declared, and then the target state is estimated by
To be more specific, some application issues are discussed as follows.
If there is no additional information, the birth density should be a uniform density over the surveillance region. For example, for Doppler component, , uniform samples are drawn from bins in the measurements which have amplitudes that exceed a predefined threshold. For echo amplitude , the birth density is uniform over , where and are expected minimum and maximum intensity levels, respectively. For monopulse ratio, , we assume that the target only exists within the half-power beamwidth, and from Figure 1 we can get that takes value within [0, 0.8]; thus, we choose its birth density to be uniform within [0, 0.8]. If other information is available (e.g., angle, range, or Doppler information supplied by the carrier aircraft, which usually has a normal law of error distribution and can be easily sampled as , the information should be used rather than the uniform one to improve the performance.
The bins in should be selected carefully, one practical choice is , where is the bin nearest to the predicted and is a design parameter. Bins near the true Doppler position have comparatively higher amplitudes and can be beneficial to the performance, while the others will, on the contrary, deteriorate the performance because the signal amplitudes there are too low. As can be seen from Figure 2, the spread function for the points that are one bin away from the true position is below −20 dB; thus we choose in this application.
5. Experiments
5.1. Experiment 1: Stationary Scenario
The radar parameters are set as follows: the wavelength is cm, the PRI is μs, and the number of pulses per CPI is . Hamming FFT windowing function is used. The target SNR represents the envelope of the target return compared to that of just noise. The SNR is measured after the entire coherent process (losses caused by windowing and straddle effect are considered). The initial relative velocity between target and radar is 1900 m/s. The initial monopulse ratio is 0.2. There are 368 bins in the clutter-free region. The initial amplitudes for 3, 6, and 10 dB are 0.87, 1.23, and 1.95, respectively. The levels of process noise used in the target model are , , and (the SNR varies only marginally). The target is born at frame 11 and disappeared at frame 51.
The particle filter parameters are set as follows: the level of the process noise is perfectly matched to the simulated data, the probabilities of target “birth” and “death” are both set as 0.05, the initial target existence probability is , the threshold , and each bin of the point spread function is quantized into points. The birth density is selected as follows: , , and uniformly distributed in the clutter-free region. The variances of the proposal distributions in the MCMC move for and are 0.04 and 0.01, respectively. and 4000 particles are used.
Figure 3 shows the estimation result of the existence probability ; asterisk signs (*) at the bottom of the figure indicate the presence of the target. It can be seen that it is possible to detect target under an SNR as low as 3 dB. Setting the threshold , for example, we can see that the target can be detected after several frames’ accumulations. From Figure 3(b) we can see that the false alarms are isolated. Thus a binary integrator can be used to mitigate them and at the same time keep the successful detections, which are continuous after becomes stable.
(a) Averaged by 200 runs
(b) Single run
Now we evaluate the detection performance of the PF-TBD algorithm in the detection terminologies. We estimate the probability of false alarm using frames 1 to 10 of the 200 Monte Carlo runs, where no target is present. More explicitly, where is the index of each Monte Carlo run. Similarly, is computed when the target is present. To see performance in the stable region as well as in the whole region, we estimate using frames 41 to 50 and frames 11 to 50, respectively. For example, using frames 41 to 50 is For comparison, the classical detector is applied to the same data. Because the PF-TBD algorithm makes one decision in each frame, for a fairly comparison, the classical detector declares a detection once any bin in the clutter-free region exceeds the threshold . Setting for both the PF-TBD and the classical detector (correspondingly, probability of false alarm for the classical detector in each single bin is and the threshold for the PF-TBD is ), the performances of them are shown in Figure 4(a). It can be observed that the of PF-TBD in the stable region at 3 dB is better than that of the classical detector at 10 dB. Thus an SNR gain of up to 7 dB is obtained.
(a) Single frame
(b) Binary integration
Taking results of Figure 4(a) as the primary detection results, we apply the 3-out-of-5 binary integration strategy to both the PF-TBD and the classical detector. Once 3 or more frames of consecutive 5 frames pass the primary detection, a secondary detection is declared. The resulting of the classical detector is 0.02, while that of the PF-TBD is 0 (no false alarm occurs in the 200 runs), which has proved that the binary integration after the PF-TBD performs well at false alarm mitigation. The in binary integration is defined as the quotient of the number of secondary detections that have past the 3-out-of-5 logic divided by the total number of secondary detections. The results are shown in Figure 4(b). We can see that the improvement over the classical detector is more compared with the single frame detection even under lower .
Remark 1. As the number of Monte Carlo runs is comparatively small, these results are not intended to provide a performance assessment. More precise results can be attained by performing a large number of Monte Carlo simulations. Compared with the classical target detection problem, it seems more reasonable to define an index to describe the delay before the becomes stable and then evaluate the detection and estimation performances in the stable region.
5.2. Experiment 2: Maneuvering Target
Now we consider a real scenario on a 2D plane. As Figure 5 has shown, the missile performs a straight motion with its antenna direction 1 degree deviated off the south to the east side. After 10 noise only frames, the target enters the main beam of the seeker radar and performs a 2 s evasive maneuver. The trajectory of the target is generated by the simulation software JSBSim (http://jsbsim.sourceforge.net/). The target’s velocity is about 280 m/s and its normal acceleration during the maneuver is 6 g. The missile’s velocity is 1200 m/s and its monopulse sum and difference beam patterns are the same as those shown in Figure 1. The echo amplitude is inversely proportional to the square of the range between missile and target (the eclipsing effect and the target fluctuation are not considered). The radar parameters are the same as those in Experiment 1 except that ; thus, the CPI is 20 ms and there are 100 target presented frames. Because the number of bins in the clutter-free region is too large, only 200 bins (bins from 3100 to 3300) containing the target are used. The initial SNR is 6 dB. The levels of process noise used in the particle filter are set as , , and . The birth density is , , and uniformly distributed in the 200 bins. The other parameters of the particle filter are the same as Experiment 1. For comparison, the PF-TBD algorithm using sum channel only is also developed and tested using the same data. The sum-only PF-TBD is obtained through omitting the in the state vector and the filtering process. To distinguish them, the filter proposed in this paper is referred to as the dual-channel PF-TBD.
In Figure 6, the estimated probabilities of existence prove the effectiveness of the two filters in target detection. Note that the sum-only filter results in worse when the target is both absent and present, which means that its detection performance is worse than that of the dual-channel one. This is because the dual-channel PF-TBD benefits from the difference channel whose amplitude is high near the half-power point.
Figures 7(a)–7(c) present the state estimation results of the two filters. We can see that both of the two filters can successfully track in target maneuvering. The dual-channel filter has better Doppler estimation performance. Note that the target Doppler can travel across half the bin size per frame; the binary integration of the classical detector will fail while that of the PF-TBD is unaffected. As the sum-only filter does not output the monopulse estimation result, the monopulse estimation performance of the dual-channel PF-TBD is compared with the classical single frame monopulse estimation method as shown in Figure 7(d). To use the same a priori knowledge, the result of the classical method is constrained to be within (0, 0.8) and that is why its estimation result is biased. The classical method assumes index of the bin which contains the target is known while the PF-TBD does not use this information. In spite of this, the monopulse estimation performance of the dual-channel PF-TBD is better.
| (a) |
| (b) |
| (c) Detail of |
| (d) |
Remark 2. This example shows that the detection performance can be improved by using the difference channel when the target is near beam edge. When the target is at the beam center, however, the difference channel amplitude is approximately zero as can be seen from Figure 1. Then the detection performance may be deteriorated instead compared with the sum-only PF-TBD. In fact, through simulation we have found that when , detection performance of the dual-channel PF-TBD is better. In practical application, the two methods should be selected according to the scenario (e.g., whether there is precise angular targeting information), and the estimation performance should also be taken into account.
6. Conclusions and Future Work
Using PF-TBD in monopulse high PRF pulse Doppler radar to improve detection and estimation performances under low SNR is addressed in this paper. The target and measurement models are analyzed and defined for this application. Based on them, a PF-TBD algorithm with resample-move operations is developed. Extensive simulations have shown that the proposed algorithm can improve both the detection and estimation performances compared with the classical and sum-only methods. To further improve the detection performance, binary integration after the PF-TBD is proposed. Simulation result shows that it can effectively mitigate the false alarms in the PF-TBD detection result.
As a byproduct of the PF-TBD algorithm, the estimated amplitude can be used to predict range eclipsing and to estimate the SNR. Application of the PF-TBD requires exact knowledge of the thermal noise power, which can be estimated on-the-fly before the PF-TBD is enabled. For seekers incorporating multispectral sensors, targeting information (e.g., angular information of the target, probability of existence of target in the main beam) from other sensors like the infrared sensor or the passive radar can be fused easily as Section 4 has stated.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.