Dual-Channel Particle Filter Based Track-Before-Detect for Monopulse Radar

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.


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 noiselimited 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][5][6].Particle filter based TBD (PF-TBD) was introduced by [7] and extended by [8][9][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 Mathematical Problems in Engineering 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.

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   = 2V  /, where  is the wavelength and V  is the radial velocity.
Although the dynamic model can be as general as x  =  −1 (x −1 , k −1 ) for a particle implementation, where k −1 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 V (1)   .The target echo amplitude and monopulse ratio are modeled as random walk processes with process noises V (2)   and V (3)  , respectively.The process noises V (1)   , V (2)   , and V (3)   are mutually independent, zero mean white noise with variances  2 (1) ,  2 (2) , and  2 (3) , 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   ∈ {0, 1}.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  = 0, 1, . . .,  − 1 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 ,  Σ/Δ () = 0 for  >  − 1 (also known as zero padding),   is the windowing function, and  = 0, 1, . . .,   − 1 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 , 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 2 2  .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, z  Σ/Δ = {  Σ/Δ (  : (  −   − 1))}, where   = ceil(2V  /(  )), V  is the horizontal velocity of the missile, and   = 1/(    ) 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  0 (⋅) 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 (x  ) 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 Z  = {z  Σ , z  Δ ,  = 1, . . ., }.It is computational complex to calculate the |  (   , )| for bins in (x  ) in real time applications.The contribution of x  to bin  in (x  ) (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 ℎ(x  , ) is where  = ∑ −1 =0   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  problem, we present an approximation approach which is calculation-free and more precise.Note that }| can be expressed as a function  w () with a parameter w = {  } and a variable  = |   −   |.Because the windowing function w can be taken as known a priori, we can quantize  w () into a number of points (e.g.,  w (  / app ),  = 0, 1, . . .,  app − 1 for  ∈ [0,   ), where  app 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 ⌊ + 0.5⌋ rounds  to the nearest integer.The result of this approximation is also presented in Figure 2.

Recursive Bayesian Filtering Procedure
The posterior probability of target existence    ≜ {  = 1 | Z  } and x  are estimated recursively by a Bayesian method as follows.Given the joint posterior PDF at frame  − 1, (x −1 ,  −1 | Z −1 ) and the latest measurement Z  , the goal is to construct the joint posterior PDF at frame , (x  ,   | Z  ).   and x  are then estimated using (x  ,   = 1 | Z  ).
Prediction.Prediction of   is given by If   = 0, x  is undefined and no prediction of it is needed.If   = 1, the prediction step of x  can be expressed as where The transitional density (x  | x −1 ,   = 1,  −1 = 1) is defined by the target dynamic model (2).The PDF   (x  ) denotes the initial target density on its appearance.
Update.The update equation using Bayes' rule is given by where the prediction density (x  ,   = 1 | Z −1 ) is given by ( 17), the normalizing constant in the denominator is where the likelihood function (z  Σ/Δ | x  ,   = 1) is given by (12).

Estimate. 𝑃 𝑘
is estimated by taking marginal of (x  ,   = 1 | Z  ) as follows: Using expected a posterior (EAP) estimator, x  is estimated by

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 resamplemove 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    , > 1, 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  = 0 and generate   samples {  0 } Step 2 (update).In the SIR filter, the prior PDF (  |   −1 ) 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 (  Σ (),   Δ () | x   ) as follows: From ( 10) and (11), (  Σ/Δ () | x   ) can be simplified as where ] and replace them using systematic resampling algorithm [18].The weights of the new samples are not required since they are all equal to 1/  .
Step 4 (MCMC move).Generate a new set of samples from [{   }   =1 ] 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    = 1.
Step 5 (state estimation).Estimate the posterior probability of target existence    by If P  exceeds a certain threshold Th ∈ (0, 1), 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 Mathematical Problems in Engineering  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 [ min ,  max ], where  min and  max 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   (x  | z  )), the information should be used rather than the uniform one to improve the performance.
The bins in (x  ) should be selected carefully, one practical choice is (x  ) = { 0 − , . . .,  0 − 1,  0 ,  0 + 1, . . .,  0 + }, where  0 is the bin nearest to the predicted x   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  = 1 in this application.

Experiment 1:
Stationary Scenario.The radar parameters are set as follows: the wavelength is  = 3 cm, the PRI is   = 4 s, and the number of pulses per CPI is  = 1000.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  2  (1) = 0.01 ⋅   ,  2  (2) = 0.001, and  2 (3) = 0.01 (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  0  = 0.1, the threshold  1 = 0.32, and each bin of the point spread function is quantized into  app = 64 points.The birth density   (x 0 | z 0 ) is selected as follows:  0 Σ ∼ (0.5, 3),  0 ∼ (0, 0.8), and  0  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. = 1 and 4000 particles are used.
Figure 3 shows the estimation result of the existence probability P  ; 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 Th = 0.6, 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 P  becomes stable.Now we evaluate the detection performance of the PF-TBD algorithm in the detection terminologies.We estimate 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 Th  .Setting  FA = 0.1 for both the PF-TBD and the classical detector (correspondingly, probability of false alarm for the classical detector in each single bin is 2.86 × 10 −4 and the threshold for the PF-TBD is Th = 0.45), 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.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  FA 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  FA .
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 P  becomes stable and then evaluate the detection and estimation performances in the stable region.

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  = 5000; thus, the CPI is 20 ms and there are 100 target presented frames.Because the number of bins in the clutterfree 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  2 (1) = 5 ⋅   ,  2 (2) = 0.05, and  2 (3) = 0.05.The birth density   (x 0 | z 0 ) is  0 Σ ∼ (1, 4),  0 ∼ (0.79,0.8), and  0  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 P  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 halfpower 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.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  > 0.1, 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.

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 sumonly methods.To further improve the detection performance, binary integration after the PF-TBD is proposed.Simulation

Figure 1 :
Figure 1: Amplitude of the sum and difference channels at different deviation angles.

Figure 3 :
Figure 3: Probability of target existence under different SNRs, asterisk signs ( * ) at the bottom indicate the presence of the target.

Figure 4 :
Figure4: Probability of detection.For single frame detection,  FA = 0.1.For binary integration (3-out-of-5), the  FA of classical detector is 0.02, while that of PF-TBD is 0 (no false alarm occurs in the 200 runs).

Figure 6 :
Figure 6: Probability of target existence (averaged by 100 runs), asterisk signs ( * ) at the bottom indicate the presence of the target.