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

Authors: Fei Cai; Hongqi Fan; Qiang Fu
Journal: Mathematical Problems in Engineering
 (2014)

---

## 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.


---

## Body

## 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.Figure 1
Amplitude of the sum and difference channels at different deviation angles.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 Section2 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(1)xk=[fdk,AΣk,γk]T,
where fdk, AΣk, and γk denote the Doppler frequency, echo amplitude of the sum channel, and monopulse ratio of the target in frame k, respectively. The Doppler frequency fd=2vr/λ, where λ is the wavelength and vr is the radial velocity.Although the dynamic model can be as general asxk=fk-1(xk-1,vk-1) for a particle implementation, where vk-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 vk(1). The target echo amplitude and monopulse ratio are modeled as random walk processes with process noises vk(2) and vk(3), respectively. The process noises vk(1), vk(2), and vk(3) are mutually independent, zero mean white noise with variances σ(1)2, σ(2)2, and σ(3)2, respectively. Thus, the system dynamic equation is
(2)xk=xk-1+T·vk-1,
where T is the CPI and vk-1=[vk-1(1),vk-1(2),vk-1(3)]T. This target model accommodates not only target maneuver but also fluctuations of the target intensity and the monopulse ratio.Target existence variableEk is modeled as a two-state Markov chain and Ek∈{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” (Pb) and “death” (Pd) as
(3)Pb≜P{Ek=1∣Ek-1=0}Pd≜P{Ek=0Ek-1=1}.
Thus the transitional probability matrix Π is given by
(4)Π=[1-PdPbPd1-Pb].


### 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 framek are denoted as sΣk(n) and sΔk(n) and given by
(5)sΣk(n)=AΣkexp⁡{j2π(fdknTr+ϕk)}+nΣk(n)(6)sΔk(n)=γkAΣkexp⁡{j2π(fdknTr+ϕk)}+nΔk(n)(7)≜AΔkexp⁡{j2π(fdknTr+ϕk)}+nΔk(n),
respectively, where AΔk is the amplitude of the difference channel, ϕk is some arbitrary phase, Tr is the pulse repetition interval (PRI), and n=0,1,…,N-1 is index of the sample in an CPI. The background thermal noises nΣk(n) and nΔk(n) are mutually independent, zero mean, and temporally white complex Gaussian processes with the same variance. The Doppler frequency fdk 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(8)yΣ/Δk(l)=∑n=0Nf-1‍sΣ/Δk(n)wnexp⁡{-j2π(nlNf)},
where the subscript Σ/Δ denotes sum channel Σ or difference channel Δ for simplification, Nf is the next power of two that is greater than or equal to N, sΣ/Δ(n)=0 for n>N-1 (also known as zero padding), wn is the windowing function, and l=0,1,…,Nf-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 k. Then the measurement can be modeled as
(9)zΣ/Δk(l)={|AΣ/ΔkBk(fdk,l)+uΣ/Δk(l)|Ek=1|uΣ/Δk(l)|Ek=0,
where |·| is the complex modulus, Bk(fdk,l)=exp⁡{j2πϕk}∑n=0N-1‍wnexp⁡{j2πfdknTr}exp⁡{-j2π(nl/Nf)}, and uΣk(l) and uΔk(l) are the background noises of the sum and difference channels, respectively, after the coherent integration. Because of linearity of the FFT, uΣk(l) and uΔk(l) are also zero mean i.i.d. complex Gaussian noise processes. We assume that they both have a variance 2σu2.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 framek; that is, zΣ/Δk={zΣ/Δk(lc:(Nf-lc-1))}, where lc=ceil(2vM/(λδf)), vM is the horizontal velocity of the missile, and δf=1/(TrNf) 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(10)p(zΣ/Δk(l)∣xk,Ek=1)=zΣ/Δk(l)σu2I0(AΣ/Δk∣Bk(fdk,l)∣zΣ/Δk(l)σu2)×exp⁡{-zΣ/Δk(l)2+(AΣ/Δk)2|Bk(fdk,l)|22σu2},
where I0(·) is the modified Bessel function of order zero. The likelihood when the target is absent has a Rayleigh distribution
(11)p(zΣ/Δk(l)∣Ek=0)=zΣ/Δk(l)σu2exp⁡{-zΣ/Δk(l)22σu2}.Because of the windowing before the FFT, the target (if present) power will spread into the bins in the vicinity of its location. LetC(xk) 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:
(12)p(zΣ/Δk∣xk,Ek=1)≈∏l∈C(xk)‍p(zΣ/Δk(l)∣xk,Ek=1)×∏l∉C(xk)‍p(zΣ/Δk(l)∣Ek=0)
and the likelihood function when the target is absent is
(13)p(zΣ/Δk∣Ek=0)=∏l=lcNf-lc-1‍p(zΣ/Δk(l)∣Ek=0).We denote the set of complete measurements up to framek as Zk={zΣi,zΔi,i=1,…,k}.It is computational complex to calculate the|Bk(fdk,l)| for bins in C(xk) in real time applications. The contribution of xk to bin l in C(xk) (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 h(xk,l) is
(14)h(xk,l)=AΣ/Δk|Bk(fdk,l)|≈AΣ/ΔkGexp⁡{-(lδf-fdk)22β2},
where G=∑n=0N-1‍wn 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 |Bk(fdk,l)|=|∑n=0N-1‍wnexp⁡{j2πnTr(fdk-lδf)}| can be expressed as a function gw(x) with a parameter w={wn} and a variable x=|fdk-lδf|. Because the windowing function w can be taken as known a priori, we can quantize gw(x) into a number of points (e.g., gw(kδf/Napp), k=0,1,…,Napp-1 for x∈[0,δf), where Napp 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,
(15)h(xk,l)≈AΣ/Δk|Bk(⌊|lδf-fdk|Nappδf+0.5⌋δfNapp,l)|,
where ⌊·⌋ is the floor function and ⌊x+0.5⌋ rounds x to the nearest integer. The result of this approximation is also presented in Figure 2.Figure 2
Comparison of different point spread function approximation methods. Hamming window is used,N=1024, Napp=64.


## 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(1)xk=[fdk,AΣk,γk]T,
where fdk, AΣk, and γk denote the Doppler frequency, echo amplitude of the sum channel, and monopulse ratio of the target in frame k, respectively. The Doppler frequency fd=2vr/λ, where λ is the wavelength and vr is the radial velocity.Although the dynamic model can be as general asxk=fk-1(xk-1,vk-1) for a particle implementation, where vk-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 vk(1). The target echo amplitude and monopulse ratio are modeled as random walk processes with process noises vk(2) and vk(3), respectively. The process noises vk(1), vk(2), and vk(3) are mutually independent, zero mean white noise with variances σ(1)2, σ(2)2, and σ(3)2, respectively. Thus, the system dynamic equation is
(2)xk=xk-1+T·vk-1,
where T is the CPI and vk-1=[vk-1(1),vk-1(2),vk-1(3)]T. This target model accommodates not only target maneuver but also fluctuations of the target intensity and the monopulse ratio.Target existence variableEk is modeled as a two-state Markov chain and Ek∈{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” (Pb) and “death” (Pd) as
(3)Pb≜P{Ek=1∣Ek-1=0}Pd≜P{Ek=0Ek-1=1}.
Thus the transitional probability matrix Π is given by
(4)Π=[1-PdPbPd1-Pb].


## 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 framek are denoted as sΣk(n) and sΔk(n) and given by
(5)sΣk(n)=AΣkexp⁡{j2π(fdknTr+ϕk)}+nΣk(n)(6)sΔk(n)=γkAΣkexp⁡{j2π(fdknTr+ϕk)}+nΔk(n)(7)≜AΔkexp⁡{j2π(fdknTr+ϕk)}+nΔk(n),
respectively, where AΔk is the amplitude of the difference channel, ϕk is some arbitrary phase, Tr is the pulse repetition interval (PRI), and n=0,1,…,N-1 is index of the sample in an CPI. The background thermal noises nΣk(n) and nΔk(n) are mutually independent, zero mean, and temporally white complex Gaussian processes with the same variance. The Doppler frequency fdk 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(8)yΣ/Δk(l)=∑n=0Nf-1‍sΣ/Δk(n)wnexp⁡{-j2π(nlNf)},
where the subscript Σ/Δ denotes sum channel Σ or difference channel Δ for simplification, Nf is the next power of two that is greater than or equal to N, sΣ/Δ(n)=0 for n>N-1 (also known as zero padding), wn is the windowing function, and l=0,1,…,Nf-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 k. Then the measurement can be modeled as
(9)zΣ/Δk(l)={|AΣ/ΔkBk(fdk,l)+uΣ/Δk(l)|Ek=1|uΣ/Δk(l)|Ek=0,
where |·| is the complex modulus, Bk(fdk,l)=exp⁡{j2πϕk}∑n=0N-1‍wnexp⁡{j2πfdknTr}exp⁡{-j2π(nl/Nf)}, and uΣk(l) and uΔk(l) are the background noises of the sum and difference channels, respectively, after the coherent integration. Because of linearity of the FFT, uΣk(l) and uΔk(l) are also zero mean i.i.d. complex Gaussian noise processes. We assume that they both have a variance 2σu2.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 framek; that is, zΣ/Δk={zΣ/Δk(lc:(Nf-lc-1))}, where lc=ceil(2vM/(λδf)), vM is the horizontal velocity of the missile, and δf=1/(TrNf) 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(10)p(zΣ/Δk(l)∣xk,Ek=1)=zΣ/Δk(l)σu2I0(AΣ/Δk∣Bk(fdk,l)∣zΣ/Δk(l)σu2)×exp⁡{-zΣ/Δk(l)2+(AΣ/Δk)2|Bk(fdk,l)|22σu2},
where I0(·) is the modified Bessel function of order zero. The likelihood when the target is absent has a Rayleigh distribution
(11)p(zΣ/Δk(l)∣Ek=0)=zΣ/Δk(l)σu2exp⁡{-zΣ/Δk(l)22σu2}.Because of the windowing before the FFT, the target (if present) power will spread into the bins in the vicinity of its location. LetC(xk) 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:
(12)p(zΣ/Δk∣xk,Ek=1)≈∏l∈C(xk)‍p(zΣ/Δk(l)∣xk,Ek=1)×∏l∉C(xk)‍p(zΣ/Δk(l)∣Ek=0)
and the likelihood function when the target is absent is
(13)p(zΣ/Δk∣Ek=0)=∏l=lcNf-lc-1‍p(zΣ/Δk(l)∣Ek=0).We denote the set of complete measurements up to framek as Zk={zΣi,zΔi,i=1,…,k}.It is computational complex to calculate the|Bk(fdk,l)| for bins in C(xk) in real time applications. The contribution of xk to bin l in C(xk) (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 h(xk,l) is
(14)h(xk,l)=AΣ/Δk|Bk(fdk,l)|≈AΣ/ΔkGexp⁡{-(lδf-fdk)22β2},
where G=∑n=0N-1‍wn 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 |Bk(fdk,l)|=|∑n=0N-1‍wnexp⁡{j2πnTr(fdk-lδf)}| can be expressed as a function gw(x) with a parameter w={wn} and a variable x=|fdk-lδf|. Because the windowing function w can be taken as known a priori, we can quantize gw(x) into a number of points (e.g., gw(kδf/Napp), k=0,1,…,Napp-1 for x∈[0,δf), where Napp 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,
(15)h(xk,l)≈AΣ/Δk|Bk(⌊|lδf-fdk|Nappδf+0.5⌋δfNapp,l)|,
where ⌊·⌋ is the floor function and ⌊x+0.5⌋ rounds x to the nearest integer. The result of this approximation is also presented in Figure 2.Figure 2
Comparison of different point spread function approximation methods. Hamming window is used,N=1024, Napp=64.


## 3. Recursive Bayesian Filtering Procedure

The posterior probability of target existencePEk≜P{Ek=1∣Zk} and xk are estimated recursively by a Bayesian method as follows. Given the joint posterior PDF at frame k-1, p(xk-1,Ek-1∣Zk-1) and the latest measurement Zk, the goal is to construct the joint posterior PDF at frame k, p(xk,Ek∣Zk). PEk and xk are then estimated using p(xk,Ek=1∣Zk).Prediction. Prediction of Ek is given by
(16)[P{Ek=1∣Zk-1}P{Ek=0∣Zk-1}]=Π[P{Ek-1=1∣Zk-1}P{Ek-1=0∣Zk-1}].
If Ek=0, xk is undefined and no prediction of it is needed. If Ek=1, the prediction step of xk can be expressed as
(17)p(xk,Ek=1∣Zk-1)=∫‍p(xk,Ek=1∣xk-1,Ek-1=1,Zk-1)di·p(xk-1,Ek-1=1∣Zk-1)dxk-1+∫‍p(xk,Ek=1∣xk-1,Ek-1=0,Zk-1)idiii·p(xk-1,Ek-1=0∣Zk-1)dxk-1,
where
(18)p(xk,Ek=1∣xk-1,Ek-1=1,Zk-1)=p(xk∣xk-1,Ek=1,Ek-1=1)P{Ek=1∣Ek-1=1}=p(xk∣xk-1,Ek=1,Ek-1=1)(1-Pd),p(xk,Ek=1∣xk-1,Ek-1=0,Zk-1)=p(xk∣xk-1,Ek=1,Ek-1=0)P{Ek=1∣Ek-1=0}=pb(xk)Pb.
The transitional density p(xk∣xk-1,Ek=1,Ek-1=1) is defined by the target dynamic model (2). The PDF pb(xk) denotes the initial target density on its appearance.Update. The update equation using Bayes’ rule is given by
(19)p(xk,Ek=1∣Zk)=p(zΣk,zΔk∣xk,Ek=1)p(xk,Ek=1∣Zk-1)p(zΣk,zΔk∣Zk-1),
where the prediction density p(xk,Ek=1∣Zk-1) is given by (17), the normalizing constant in the denominator is p(zΣk,zΔk∣Zk-1)=∫‍p(zΣk,zΔk∣x,Ek=1)p(x,Ek=1∣Zk-1)dx, and the likelihood function p(zΣk,zΔk∣xk,Ek=1) is
(20)p(zΣk,zΔk∣xk,Ek=1)=p(zΣk∣xk,Ek=1)p(zΔk∣xk,Ek=1),
where the likelihood function p(zΣ/Δk∣xk,Ek=1) is given by (12).Estimate. PEk is estimated by taking marginal of p(xk,Ek=1∣Zk) as follows:
(21)P^Ek=∫‍p(xk,Ek=1∣Zk)dxk.
Using expected a posterior (EAP) estimator, xk is estimated by
(22)x^k=∫xkp(xk,Ek=1∣Zk)dxkP^Ek.


## 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 qm(γk′∣γk) is defined, from which a sample is drawn for each particle after resampling. A monopulse ratio γk′ is obtained conditioned on the old monopulse ratio γk while keeping the other two states unchanged. Under the assumption that the proposal is symmetric, qm(γk′∣γk)=qm(γk∣γk′), the new particle is accepted or rejected on a test, formed by a ratio of likelihoods
(23)Tγ′,γ=L(zk∣fdk,AΣk,γk′)L(zk∣fdk,AΣk,γk).
If Tγ′,γ>1, then the new particle, with monopulse ratio γ′, is kept. Otherwise the new particle is kept in preference to the old particle only if U<Tγ′,γ, where U is a uniform random number between 0 and 1. The move operation is used twice in this application, firstly to the amplitude AΣk and then to the monopulse ratio γk. Truncated Gaussian distributions with different variances and means at AΣk and γk, respectively, are used as the proposal distributions.A detailed description of the TBD algorithm is given as follows.Initialization. Set k=0 and generate Ns samples {E0i}i=1Ns from PE0=P(E0=1). If E0i=1, generate x0i from the birth density qb(x0∣z0), or else, x0i is undefined.Then, given[{xk-1i}i=1Ns,zk] at each frame k, go from Steps 1 to 5.Step 1 (prediction).
Generate{Eki}i=1Ns on the basis of {Ek-1i}i=1Ns and Π. If Eki=0, xki is undefined. If Ek-1i=1 and Eki=1, predict xki according to (2). For the new born particles, that is, those with Ek-1i=0 and Eki=1, generate xki from the birth density qb(xk∣zk).Step 2 (update).
In the SIR filter, the prior PDFp(zk∣xk-1i) 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]:
(24)w~ki={∏l∈Ci(xk)‍L(zΣk(l),zΔk(l)∣xki)ifEki=11ifEki=0.
We simplify the likelihood L(zΣk(l),zΔk(l)∣xki) as follows:
(25)L(zΣk(l),zΔk(l)∣xki)≜p(zΣk(l),zΔk(l)∣xki,Ek=1)p(zΣk(l),zΔk(l)∣Ek=0)=p(zΣk(l)∣xki,Ek=1)p(zΔk(l)∣xki,Ek=1)p(zΣk(l)∣Ek=0)p(zΔk(l)∣Ek=0)=L(zΣk(l)∣xki)L(zΔk(l)∣xki).
From (10) and (11), L(zΣ/Δk(l)∣xki) can be simplified as
(26)L(zΣ/Δk(l)∣xki)=p(zΣ/Δk(l)∣xki,Ek=1)p(zΣ/Δk(l)∣Ek=0)=I0(|AΣ/Δ(k,i)Bk(fd(k,i),l)|zΣ/Δk(l)σu2)×exp⁡{-|AΣ/Δ(k,i)Bk(fd(k,i),l)|22σu2},
where AΔ(k,i) is calculated by AΔ(k,i)=γkiAΣ(k,i). Then get the normalized weights {wki}i=1Ns by wki=w~ki/∑i=1Ns‍w~ki.Step 3 (resample).
Generate a new set of samples[{Eki,xki}i=1Ns] from [{Eki,xki,wki}i=1Ns] and replace them using systematic resampling algorithm [18]. The weights of the new samples are not required since they are all equal to 1/Ns.Step 4 (MCMC move).
Generate a new set of samples from[{xki}i=1Ns] and replace them by move of AΣk using the Metropolis-Hasting method described above; do this again by move of γk. Note that this operation only changes the particles with Eki=1.Step 5 (state estimation).
Estimate the posterior probability of target existencePEk by
(27)P^Ek=∑i=1NsEkiNs.
If P^Ek exceeds a certain threshold Th∈(0,1), target presence is declared, and then the target state is estimated by
(28)x^k∣k=∑i=1NsxkiEki∑i=1NsEki.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,fdk, uniform samples are drawn from bins in the measurements which have amplitudes that exceed a predefined threshold. For echo amplitude Ak, the birth density is uniform over [Amin⁡,Amax⁡], where Amin⁡ and Amax⁡ are expected minimum and maximum intensity levels, respectively. For monopulse ratio, γk, 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 qb(xk∣zk)), the information should be used rather than the uniform one to improve the performance.The bins inC(xk) should be selected carefully, one practical choice is C(xk)={i0-p,…,i0-1,i0,i0+1,…,i0+p}, where i0 is the bin nearest to the predicted xki and p 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 p=1 in this application.


## 5. Experiments

### 5.1. Experiment 1: Stationary Scenario

The radar parameters are set as follows: the wavelength isλ=3 cm, the PRI is Tr=4 μs, and the number of pulses per CPI is N=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 σ(1)2=0.01·δf, σ(2)2=0.001, and σ(3)2=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”Pb and “death” Pd are both set as 0.05, the initial target existence probability is PE0=0.1, the threshold T1=0.32, and each bin of the point spread function is quantized into Napp=64 points. The birth density qb(x0∣z0) is selected as follows: AΣ0~U(0.5,3), γ0~U(0,0.8), and fd0 uniformly distributed in the clutter-free region. The variances of the proposal distributions in the MCMC move for AΣ and γ are 0.04 and 0.01, respectively. p=1 and 4000 particles are used.Figure3 shows the estimation result of the existence probability P^Ek; 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^Ek becomes stable.Probability of target existence under different SNRs, asterisk signs (*) at the bottom indicate the presence of the target.
(a)
Averaged by 200 runs
(b)
Single runNow we evaluate the detection performance of the PF-TBD algorithm in the detection terminologies. We estimate the probability of false alarmPFA using frames 1 to 10 of the 200 Monte Carlo runs, where no target is present. More explicitly,
(29)PFA=1200×10∑m=1200‍∑k=110‍(P^Ek,m≷01Th),
where m is the index of each Monte Carlo run. Similarly, PD is computed when the target is present. To see performance in the stable region as well as in the whole region, we estimate PD using frames 41 to 50 and frames 11 to 50, respectively. For example, PD using frames 41 to 50 is
(30)PD=1200×10∑m=1200‍∑k=4150‍(P^Ek,m≷01Th).
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 PFA=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 PD performances of them are shown in Figure 4(a). It can be observed that the PD 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.Probability of detection. For single frame detection,PFA=0.1. For binary integration (3-out-of-5), the PFA of classical detector is 0.02, while that of PF-TBD is 0 (no false alarm occurs in the 200 runs).
(a)
Single frame
(b)
Binary integrationTaking results of Figure4(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 PFA 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 PD 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 PD results are shown in Figure 4(b). We can see that the PD improvement over the classical detector is more compared with the single frame detection even under lower PFA.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 theP^Ek 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 Figure5 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 N=5000; 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 σ(1)2=5·δf, σ(2)2=0.05, and σ(3)2=0.05. The birth density qb(x0∣z0) is AΣ0~U(1,4), γ0~U(0.79,0.8), and fd0 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 γk 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.Figure 5
Missile and target trajectories. The “□” and “△” denote start and end of the trajectory, respectively.In Figure6, 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^Ek 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.Figure 6
Probability of target existence (averaged by 100 runs), asterisk signs (*) at the bottom indicate the presence of the target.Figures7(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.Estimation results ofAΣ, fd, and γ. The thick dashed lines show the mean value over 100 Monte Carlo runs. The thin dashed lines are mean ± one standard deviation.
(a)
A
Σ
(b)
f
d
(c)
Detail offd
(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 Figure1. 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.


## 5.1. Experiment 1: Stationary Scenario

The radar parameters are set as follows: the wavelength isλ=3 cm, the PRI is Tr=4 μs, and the number of pulses per CPI is N=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 σ(1)2=0.01·δf, σ(2)2=0.001, and σ(3)2=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”Pb and “death” Pd are both set as 0.05, the initial target existence probability is PE0=0.1, the threshold T1=0.32, and each bin of the point spread function is quantized into Napp=64 points. The birth density qb(x0∣z0) is selected as follows: AΣ0~U(0.5,3), γ0~U(0,0.8), and fd0 uniformly distributed in the clutter-free region. The variances of the proposal distributions in the MCMC move for AΣ and γ are 0.04 and 0.01, respectively. p=1 and 4000 particles are used.Figure3 shows the estimation result of the existence probability P^Ek; 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^Ek becomes stable.Probability of target existence under different SNRs, asterisk signs (*) at the bottom indicate the presence of the target.
(a)
Averaged by 200 runs
(b)
Single runNow we evaluate the detection performance of the PF-TBD algorithm in the detection terminologies. We estimate the probability of false alarmPFA using frames 1 to 10 of the 200 Monte Carlo runs, where no target is present. More explicitly,
(29)PFA=1200×10∑m=1200‍∑k=110‍(P^Ek,m≷01Th),
where m is the index of each Monte Carlo run. Similarly, PD is computed when the target is present. To see performance in the stable region as well as in the whole region, we estimate PD using frames 41 to 50 and frames 11 to 50, respectively. For example, PD using frames 41 to 50 is
(30)PD=1200×10∑m=1200‍∑k=4150‍(P^Ek,m≷01Th).
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 PFA=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 PD performances of them are shown in Figure 4(a). It can be observed that the PD 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.Probability of detection. For single frame detection,PFA=0.1. For binary integration (3-out-of-5), the PFA of classical detector is 0.02, while that of PF-TBD is 0 (no false alarm occurs in the 200 runs).
(a)
Single frame
(b)
Binary integrationTaking results of Figure4(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 PFA 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 PD 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 PD results are shown in Figure 4(b). We can see that the PD improvement over the classical detector is more compared with the single frame detection even under lower PFA.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 theP^Ek 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 Figure5 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 N=5000; 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 σ(1)2=5·δf, σ(2)2=0.05, and σ(3)2=0.05. The birth density qb(x0∣z0) is AΣ0~U(1,4), γ0~U(0.79,0.8), and fd0 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 γk 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.Figure 5
Missile and target trajectories. The “□” and “△” denote start and end of the trajectory, respectively.In Figure6, 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^Ek 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.Figure 6
Probability of target existence (averaged by 100 runs), asterisk signs (*) at the bottom indicate the presence of the target.Figures7(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.Estimation results ofAΣ, fd, and γ. The thick dashed lines show the mean value over 100 Monte Carlo runs. The thin dashed lines are mean ± one standard deviation.
(a)
A
Σ
(b)
f
d
(c)
Detail offd
(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 Figure1. 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.


## 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 Section4 has stated.



---

*Source: 750279-2014-08-18.xml*