# An Adaptive PHD Filter for Multitarget Tracking with Multispectral Data Fusion

Authors: Guoliang Zhang; Chunling Yang; Yan Zhang
Journal: Journal of Spectroscopy
 (2015)

---

## Abstract

In order to improve the detection and tracking performance of multiple targets from IR multispectral image sequences, the approach based on spectral fusion algorithm and adaptive probability hypothesis density (PHD) filter is proposed. Firstly, the nonstationary adaptive suppression method is proposed to remove the background clutter. Based on the multispectral image sequence, the spectral fusion method is used to detect the abnormal targets. Spectral fusion produces the appropriate binary detection model and the computational probability of detection. Secondly, the particle filtering-based adaptive PHD algorithm is developed to detect and track multiple targets. This algorithm can deal with the nonlinear measurement on target state. In addition, the calculated probability of detection substitutes the fixed detection probability in PHD filter. Finally, the synthetic data sets based on various actual background images were utilized to validate the effectiveness of the detection approach. The results demonstrate that the proposed approach outperforms the conventional sequential PHD filtering in terms of detection and tracking performances.


---

## Body

## 1. Introduction

The conventional detection scheme of the long range point source target mainly focuses on employing the temporal features of target from broadband (single band) IR sensor [1, 2]. However, this methodology faces enormous challenges because of the latest decoys in modern electronic countermeasures. It is difficult to detect the decoys and targets by employing the temporal information when they have similar trajectories. It is noted that the spectral information, which is obtained through a large number of narrow contiguous spectral channels, denotes the intrinsic feature of target [3]. Therefore multispectral imaging sensors can portray more information than a broadband IR sensor. Multispectral IR image has been widely used to perform target detection by means of fusing the narrowband information [4–6]. However, the detection mechanism based on multispectral IR image only adopts the spectral data and cannot be used to detect moving targets with low signal-to-noise ratio (SNR). The fusion of the data in time, wavelength, or other observables should be explored to improve the ability to detect the targets with a low SNR.In addition, multiple target detection and tracking is a challenging problem that arises from heavy clutter, noise, and false alarm. The conventional approaches focused on the data association technique and the appropriate filtering methods (such as Kalman filtering). Typical data association algorithms include multiple hypothesis tracking (MHT), joint probabilistic data association (JPDA), and their own variants [7–9]. However, it has been demonstrated that data association-based tracking scheme is a NP-hard problem. Thus free data association tracking methods are hot topic in recent years. The Finite Random Set Statistics- (FRSST-) based PHD filter, which incorporates target birth, target spawn, and target death, has been extensively studied in multitarget detection and tracking applications [10–15]. But PHD filter is essentially approximation of the first moment of multitarget posterior density. This factor makes PHD filter performance degraded in low SNR environment. Moreover, owing to the difficulty of obtaining the suitable multisensor model for PHD filter, the preciously presented PHD filter employed the single sensor model in most cases. To solve this problem, Mahler proposed a sequential PHD filter to fuse measurements from multiple sensors. Some researchers also utilized the sequential PHD filter to track radar targets [16]. However, the sequential PHD filtering produces different updated PHD depending on the order of the multisensor update [17]. In the case that the probability of detection is not high for at least one sensor, the resulting intensity and estimates significantly rely on the order of the multisensor update. Therefore it is hard to determine the correct fusion order of multiple sensors.Though possessing attractive and promising application, the PHD filtering still has some problems not taken into consideration concretely. Actually, the PHD filter is mainly used in radar target area. The dynamics model and measurement model for radar target are mainly proposed based on active detection scheme. At present the PHD filter is rarely used for infrared targets detection. Two problems arise from IR target detection when using PHD filtering. Firstly the measurements which are produced by IR sensor are nonlinear and it is hard to formulate the measurement model by FRSST framework. Secondly, the low SNR targets immersed in IR images make the PHD filtering hard to implement.This paper, based on the PHD filtering, is aimed at detecting and tracking multiple targets from IR multispectral image sequences under diverse scenarios. The contribution of the paper can be summarized as follows.(1)
Constructing the general four-dimensional signal model for IR multispectral image sequences: the multitarget detection and tracking problem can be formulated in the adopted model which can be used in the airborne system or in the satellite-borne system.(2)
Proposing a method to suppress the background clutter noise: the proposed method differs from the traditional approaches on the aspect of removing the spatial correlation and temporal correlation in pixels simultaneously. The solution is put forward in accordance with the least square theorem. This method is also considered to be the preprocessing of multispectral image fusion.(3)
Developing the PHD filtering to detect and track multiple targets: an adaptive PHD (APHD) filtering algorithm is derived from using the conventional PHD method and the RX algorithm. The difference between the APHD filtering and the original PHD filtering is that the probability of detection and the probability of false alarm in our study are calculated based on the measured data. Therefore the presented approach may be more appropriate for the real physical model and the final results could have a higher degree of confidence. Additionally, the efficient particle filtering is also developed to implement the APHD filtering.The experimental results have shown that the proposed algorithms provide good detection and tracking capabilities, even for the targets with low signal-to-noise ratio. The paper is organized as follows. In Section2, the original PHD filtering is introduced. In Section 3, the general four-dimensional signal model and target imaging model are presented; meanwhile the adaptive suppression of background clutter is proposed. The multispectral data fusion is achieved by means of the RX algorithm. Section 4 derives the adaptive PHD filtering based on efficient particle filtering. The parameters in measurement update are described. In Section 5, the detailed performance results for the proposed detection and tracking algorithm are presented. The discussions and future work are stated in Section 6. The conclusions are given in Section 7.


## 2. The PHD Filtering

The PHD filter recursively estimates the first-order statistical moments of multitarget posterior density and is the suboptimal algorithm underlying the optimal Bayesian filtering framework. However, the PHD filter avoids the data association procedure and reduces the computational load. The PHD filter is described as follows by using Bayesian laws.The posterior intensity of FISST at framek-1 is denoted as Dk-1(x) and the predicted intensity is given by Bayesian laws. Consider(1)Dk∣k-1(x)=bk(x)+∫[ps,k(φ)fk∣k-1(x∣φ)+αk∣k-1(x∣φ)]Dk-1(x)dφ.When the new measurement zk is obtained, the intensity Dk∣k-1(x) is updated by(2)Dkx=∑zk∈ZkpD,k(x)l(z∣x)Dk∣k-1(x)κk(z)+∫pD,k(φ)l(z∣φ)Dk∣k-1(φ)dφ+1-pD,kxDk|k-1x,where bk(x) is the intensity of birth RFS at time k, ps,k(φ) is the keeping probability of target from time k-1 to k, and fk∣k-1(x∣φ) is the probability density of a target at time k, given that its previous state is φ. αk∣k-1(x∣φ) is the intensity of the RFS of targets spawned at time k by a target with previous state φ. pD,k(x) is the probability of detection given a state x at time k. l(z∣x) is the single-target measurement likelihood at time k. κk(z) is the intensity of clutter RFS at time k. It is usually defined as κk(z)=cλ(z), where c is subject to Poisson distribution and λ(z) is the probability density of uniform distribution. Zk is the available multitarget observations at time k. The integration of PHD over any region S from single-target state is the expected number of targets contained in S; that is, Nk∣k=∫SDk(x)dx. The intensity Dk(x) has many peak values (Nk∣k), and each peak value denotes the expectation state. EM algorithm and cluster algorithm can be used to extract multitarget states.It is noted that some assumptions are made in order to use the PHD filter. The sufficient conditions are as follows.Condition 1.
Targets evolve in time and generate measurements independently of one another.Condition 2.
The clutter RFS is Poisson and is independent of the measurements.Condition 3.
The predicted multitarget RFS is Poisson.In the real applications, all the three sufficient conditions should be taken into account. However, some studies neglected the third sufficient condition. The third sufficient condition is an approximation to update the predictedDk∣k-1(x) to obtain the posterior intensity Dk(x). This approximation is justifiable under the scenario that false alarm densities are small [10]. In this paper, the third sufficient condition is satisfied via suppression of background clutter and anomaly detection algorithm, which reduce false alarm densities.


## 3. Four-Dimensional Signal Model in Multispectral Image Sequence

### 3.1. Target State and Binary Detection Model

The multispectral images haveL bands. Each band consists of K frames. Each frame consists of R pixels in the x-dimension and C pixels in the y-dimension. The target state vector is xk=[rkvrkckvck]T, where (rk,ck) denotes the position and (vrk,vck) refers to the velocity of the target at frame k. The time evolution of target state can be modeled as(3)xk=fxk-1,nk-1,where f(·) is the nonlinear/linear state function, which depends on the target state xk-1 and the process noise nk-1. nk-1 is usually subject to Gaussian distribution with zero mean.It is assumed that each pixel only depends on one target at most. Thus the raw binary measurement signal model can be presented by the four-dimensional signal space. Conisder(4)H0:zrckl=brcklnotargetH1:zrckl=trckl+brckltarget,where zrckl is the received spectral irradiance at the position (rk,ck) of the image. brckl denotes the spectral irradiance of background and trckl represents the spectral irradiance of target.


### 3.2. Target Imaging Model

The target which is distant from the sensor covers one pixel or subpixel in the image. Since the IR imaging system is nonideal, the point source target has several pixels. The target imaging size is decided by the point spread function (PSF). PSF is usually approximated by a truncate 2D Gaussian function. PSF can be represented as(5)gr,c=exp⁡-r22σr2+c22σc2,where the two constants σr and σc are the spread parameters. Both can be obtained via experiments. In this paper we have σr=σc=1.


### 3.3. Adaptive Suppression of Background Clutter

It is necessary to remove clutter noise for the detection of dim point source target embedded in dense background clutter. The conventional methods are based on spatial domain, temporal domain, and frequency domain. However, the spatial correlation exists in the multispectral image sequences besides the temporal correlation. In order to remove background clutter, novel methods for suppressing clutter should be further explored.The spectral radiant intensity is defined asγrckl at wavelength l, and the mean of local spatial pixels is calculated by using the spatial filter. In order to avoid influencing the local spatial mean by the target pixels, it is supposed that the target covers a specific area in the image, which is shown in Figure 1.Figure 1
Pixels used for mean value estimation in spatial domain and temporal domain.In this paper, the spatial correlation and temporal correlation are removed simultaneously. Firstly, letI represent the multispectral image. Let γ(rk,ck) be a test pixel of interest in I. The PSF determines the covering area of target in I, which is shown in Figure 1. The corresponding local temporal mean is calculated by four nearest neighbors in temporal domain. The residual error after removing the spatial and temporal mean is presented as (6)γ~rckl=γrckl-ϕγ-rcklspatial-βγ-rckltemporal,where(7)γ-rcklspatial=∑γrcklspatial32,γ-rckltemporal=∑γrckltemporal4.It is denoted as(8)θ=ϕβ,H=γ-rc0lspatialγ-rc1lspatial⋯γ-rc(K-1)lspatialγ-rc0ltemporalγ-rc1ltemporal⋯γ-rc(K-1)ltemporalT,Λ=γrc0lγrc1l⋯γrcK-1lT.The solution of θ can be given by the least square theorem. It is given by (9)γ~rc0lγ~rc1l⋯γ~rcK-1lT=Λ-HHTH-1HTΛ.


### 3.4. Measurement Model Based on Multispectral Data Fusion

It is critical to acquire available measurements in the tracking system. The RX algorithm is utilized to produce measurements in this section. RX algorithm, as the benchmark method in multispectral applications, was proposed by Reed and Yu in 1990 [18]. It is an abnormal detection algorithm, supposing that the spectral distribution of target and the covariance of background clutter are unknown. The local nonstationary background mean should be removed when using RX to generate detections. This operation has been implemented in Section 3.3. It is supposed that we have J bands images; each residual pixel after background suppression is denoted as x(n)k=[x1(n)k,x2(n)k,…,xj(n)k]T. Thus sample matrix whose size is J×N at time k is defined as (10)Xk=x1k,x2k,…,xNk,where N is the number of samples, which come from the outside window of Figure 1. The binary detection model of RX algorithm is given by(11)H0:Xk=WknotargetH1:Xk=Wk+bkSktarget,where the hypothesis H0 shows that no target exists in the scene. Wk is the background clutter noise matrix. H1 indicates that the targets stay in the environment. bk=[b1k,b2k,…,bJk]T is a J vector of unknown signal intensities corresponding to the J scenes or channels. Sk=[s1k,s2k,…,sNk] is the signal pattern. Wk is subject to approximate Gaussian distribution. The RX detector is derived from the generalized likelihood ratio test (GLRT). Consider(12)r(Xk)=(XkST)XkXkT-1XkSTTSST≥r0,H1<r0,H0,where r(Xk) denotes the test statistics and r0 is the detection threshold. r0 can be calculated by the probability of detection (PD) and the probability of false alarm (PFA) according to the Neman-Pearson rules. The PD and PFA are given by (13)PD=∫r01f(r∣H1)dr,PFA=∫r01fr∣H0dr,where f(r∣H1) is a noncentral beta-distribution and f(r∣H0) is a standard beta-distribution. f(r∣H1) is represented as(14)fr∣H1=Γ(N/2)Γ((N-J)/2)Γ(J/2)1-r(N-J-2)/2r(J-2)/2·exp⁡⁡-a21F1(N2;J2;ar2)0<r<1,where F1(N/2;J/2;ar/2) is the confluent hypergeometric function and a denotes the generalized signal-to-noise ratio.f
(
r
∣
H
0
) is presented as(15)fr∣H0=Γ(N/2)Γ((N-J)/2)Γ(J/2)1-r(N-J-2)/2r(J-2)/2hhhhhhhhhhhhhhhhhhhhhhhhhhhh0<r<1.The data of spectral domain is fused by RX algorithm and the available measurements are produced. Then the fusion tracking can be implemented by using the test statistics of multiple frames. The measurements in each frame can be formulated by the random finite set.


## 3.1. Target State and Binary Detection Model

The multispectral images haveL bands. Each band consists of K frames. Each frame consists of R pixels in the x-dimension and C pixels in the y-dimension. The target state vector is xk=[rkvrkckvck]T, where (rk,ck) denotes the position and (vrk,vck) refers to the velocity of the target at frame k. The time evolution of target state can be modeled as(3)xk=fxk-1,nk-1,where f(·) is the nonlinear/linear state function, which depends on the target state xk-1 and the process noise nk-1. nk-1 is usually subject to Gaussian distribution with zero mean.It is assumed that each pixel only depends on one target at most. Thus the raw binary measurement signal model can be presented by the four-dimensional signal space. Conisder(4)H0:zrckl=brcklnotargetH1:zrckl=trckl+brckltarget,where zrckl is the received spectral irradiance at the position (rk,ck) of the image. brckl denotes the spectral irradiance of background and trckl represents the spectral irradiance of target.


## 3.2. Target Imaging Model

The target which is distant from the sensor covers one pixel or subpixel in the image. Since the IR imaging system is nonideal, the point source target has several pixels. The target imaging size is decided by the point spread function (PSF). PSF is usually approximated by a truncate 2D Gaussian function. PSF can be represented as(5)gr,c=exp⁡-r22σr2+c22σc2,where the two constants σr and σc are the spread parameters. Both can be obtained via experiments. In this paper we have σr=σc=1.


## 3.3. Adaptive Suppression of Background Clutter

It is necessary to remove clutter noise for the detection of dim point source target embedded in dense background clutter. The conventional methods are based on spatial domain, temporal domain, and frequency domain. However, the spatial correlation exists in the multispectral image sequences besides the temporal correlation. In order to remove background clutter, novel methods for suppressing clutter should be further explored.The spectral radiant intensity is defined asγrckl at wavelength l, and the mean of local spatial pixels is calculated by using the spatial filter. In order to avoid influencing the local spatial mean by the target pixels, it is supposed that the target covers a specific area in the image, which is shown in Figure 1.Figure 1
Pixels used for mean value estimation in spatial domain and temporal domain.In this paper, the spatial correlation and temporal correlation are removed simultaneously. Firstly, letI represent the multispectral image. Let γ(rk,ck) be a test pixel of interest in I. The PSF determines the covering area of target in I, which is shown in Figure 1. The corresponding local temporal mean is calculated by four nearest neighbors in temporal domain. The residual error after removing the spatial and temporal mean is presented as (6)γ~rckl=γrckl-ϕγ-rcklspatial-βγ-rckltemporal,where(7)γ-rcklspatial=∑γrcklspatial32,γ-rckltemporal=∑γrckltemporal4.It is denoted as(8)θ=ϕβ,H=γ-rc0lspatialγ-rc1lspatial⋯γ-rc(K-1)lspatialγ-rc0ltemporalγ-rc1ltemporal⋯γ-rc(K-1)ltemporalT,Λ=γrc0lγrc1l⋯γrcK-1lT.The solution of θ can be given by the least square theorem. It is given by (9)γ~rc0lγ~rc1l⋯γ~rcK-1lT=Λ-HHTH-1HTΛ.


## 3.4. Measurement Model Based on Multispectral Data Fusion

It is critical to acquire available measurements in the tracking system. The RX algorithm is utilized to produce measurements in this section. RX algorithm, as the benchmark method in multispectral applications, was proposed by Reed and Yu in 1990 [18]. It is an abnormal detection algorithm, supposing that the spectral distribution of target and the covariance of background clutter are unknown. The local nonstationary background mean should be removed when using RX to generate detections. This operation has been implemented in Section 3.3. It is supposed that we have J bands images; each residual pixel after background suppression is denoted as x(n)k=[x1(n)k,x2(n)k,…,xj(n)k]T. Thus sample matrix whose size is J×N at time k is defined as (10)Xk=x1k,x2k,…,xNk,where N is the number of samples, which come from the outside window of Figure 1. The binary detection model of RX algorithm is given by(11)H0:Xk=WknotargetH1:Xk=Wk+bkSktarget,where the hypothesis H0 shows that no target exists in the scene. Wk is the background clutter noise matrix. H1 indicates that the targets stay in the environment. bk=[b1k,b2k,…,bJk]T is a J vector of unknown signal intensities corresponding to the J scenes or channels. Sk=[s1k,s2k,…,sNk] is the signal pattern. Wk is subject to approximate Gaussian distribution. The RX detector is derived from the generalized likelihood ratio test (GLRT). Consider(12)r(Xk)=(XkST)XkXkT-1XkSTTSST≥r0,H1<r0,H0,where r(Xk) denotes the test statistics and r0 is the detection threshold. r0 can be calculated by the probability of detection (PD) and the probability of false alarm (PFA) according to the Neman-Pearson rules. The PD and PFA are given by (13)PD=∫r01f(r∣H1)dr,PFA=∫r01fr∣H0dr,where f(r∣H1) is a noncentral beta-distribution and f(r∣H0) is a standard beta-distribution. f(r∣H1) is represented as(14)fr∣H1=Γ(N/2)Γ((N-J)/2)Γ(J/2)1-r(N-J-2)/2r(J-2)/2·exp⁡⁡-a21F1(N2;J2;ar2)0<r<1,where F1(N/2;J/2;ar/2) is the confluent hypergeometric function and a denotes the generalized signal-to-noise ratio.f
(
r
∣
H
0
) is presented as(15)fr∣H0=Γ(N/2)Γ((N-J)/2)Γ(J/2)1-r(N-J-2)/2r(J-2)/2hhhhhhhhhhhhhhhhhhhhhhhhhhhh0<r<1.The data of spectral domain is fused by RX algorithm and the available measurements are produced. Then the fusion tracking can be implemented by using the test statistics of multiple frames. The measurements in each frame can be formulated by the random finite set.


## 4. Adaptive PHD Filtering Based on Particle Filtering

The PHD filter is used to track multiple targets by using the test statistics that were produced by RX. At present two methods are proposed to implement PHD filter: the Gaussian Mixture- (GM-) PHD and the Sequential Monte Carlo- (SMC-) PHD, that is, the particle filtering-based PHD [19, 20]. Because the target measurement from IR image is nonlinear, the particle filtering-based method is developed to implement PHD filter in this paper. However, the conventional approach based on particle filtering may produce sample degeneration and impoverishment. A new proposal distribution is proposed to deal with sample degeneration and impoverishment.The SMC approach is utilized to implement the algorithm which provides a mechanism to represent the posterior density by a set of random samples or particles. The posteriorDk-1∣k-1(xk-1∣Z1:k-1) can be represented by a set of particles {ωk-1(p),xk-1(p)}, p=1,2,…,Lk-1. Consider(16)Dk-1∣k-1(xk-1∣Z1:k-1)=∑p=1Lk-1ωk-1pδxk-1-xk-1p.


### 4.1. Time Update

The predicted particles are generated via the following distribution:(17)xk∣k-1(p)~qk(·∣xk(p)),p=1,2,…,Lk-1pk(·∣xk(p)),p=Lk-1+1,Lk-1+2,…,Lk-1+Jk,where qk(·∣xk(p)) and pk(·∣xk(p)) are a priori distribution. The predicted PHD can be expressed as (18)Dk∣k-1xk∣k-1∣Z1:k-1=∑p=1Lk-1+Jkωk∣k-1pδxk∣k-1-xk∣k-1p,where ωk∣k-1(p) is calculated according to (1):(19)ωk∣k-1(p)~psk∣k-1(xk∣k-1(p))fk(xk∣k-1(p)∣xk-1(p))+αk∣k-1(xk∣k-1(p)∣xk-1(p))qk(xk∣k-1(p)∣xk-1(p)),p=1,2,…,Lk-1bk(xk∣k-1(p))pk(xk∣k-1(p)∣xk-1(p)),p=Lk-1+1,Lk-1+2,…,Lk-1+Jk.


### 4.2. Measurement Update

When the available measurements are acquired at timek, the update PHD can be given by (2) and (18). Consider(20)Dk(xk∣Z1:k)=∑p=1Lk-1+Jkωk∗pδxk-xk∣k-1p,where ωk∗(p) is (21)ωk∗(p)=∑zk∈ZkPD(x)lk(zk∣xk∣k-1(p))ωk∣k-1(p)κk(z)+∑p=1Lk-1+JkpD,k(x)lk(zk∣xk∣k-1(p))ωk∣k-1(p)+1-PDxωk∣k-1p.The conventional PHD filtering always sets a fixed PD, which rarely exists in the real situation. In this paper, the specified PD is substituted by a computed PD according to (13). Thus the proposed algorithm is entitled by adaptive PHD (APHD). The likelihood function lk(zk∣xk|k-1(p)) in (21) is defined as (22)lkzk∣xk∣k-1p=∏i∈Ci(xk∣k-1p)∏j∈Cj(xk∣k-1p)lkzki,j∣xk∣k-1p,where Ci(xk∣k-1(p)) and Cj(xk∣k-1(p)) denote the covering area of the target in the image, which has been given in Figure 1. lk(zk(i,j)∣xk|k-1(p)) is the single pixel likelihood, which is expressed via (14) and (15). Consider(23)lkzki,j∣xk∣k-1p=f(zk(i,j)∣H1)f(zk(i,j)∣H0)=exp⁡-a21F1N2;J2;ar2.The expected number of targets is calculated by summing up the total weights; that is,(24)N^k=∑p=1Lk-1+Jkωk∗(p).Then resampling is performed to avoid the degeneracy. The updated posterior density at time stepk is presented as (25)Dk(xk∣Z1:k)=∑p=1Lkωkpδxk-xkp.In summary, the processing chain of the proposed approach is presented in Figure2.Figure 2
The processing chain of the proposed approach.


## 4.1. Time Update

The predicted particles are generated via the following distribution:(17)xk∣k-1(p)~qk(·∣xk(p)),p=1,2,…,Lk-1pk(·∣xk(p)),p=Lk-1+1,Lk-1+2,…,Lk-1+Jk,where qk(·∣xk(p)) and pk(·∣xk(p)) are a priori distribution. The predicted PHD can be expressed as (18)Dk∣k-1xk∣k-1∣Z1:k-1=∑p=1Lk-1+Jkωk∣k-1pδxk∣k-1-xk∣k-1p,where ωk∣k-1(p) is calculated according to (1):(19)ωk∣k-1(p)~psk∣k-1(xk∣k-1(p))fk(xk∣k-1(p)∣xk-1(p))+αk∣k-1(xk∣k-1(p)∣xk-1(p))qk(xk∣k-1(p)∣xk-1(p)),p=1,2,…,Lk-1bk(xk∣k-1(p))pk(xk∣k-1(p)∣xk-1(p)),p=Lk-1+1,Lk-1+2,…,Lk-1+Jk.


## 4.2. Measurement Update

When the available measurements are acquired at timek, the update PHD can be given by (2) and (18). Consider(20)Dk(xk∣Z1:k)=∑p=1Lk-1+Jkωk∗pδxk-xk∣k-1p,where ωk∗(p) is (21)ωk∗(p)=∑zk∈ZkPD(x)lk(zk∣xk∣k-1(p))ωk∣k-1(p)κk(z)+∑p=1Lk-1+JkpD,k(x)lk(zk∣xk∣k-1(p))ωk∣k-1(p)+1-PDxωk∣k-1p.The conventional PHD filtering always sets a fixed PD, which rarely exists in the real situation. In this paper, the specified PD is substituted by a computed PD according to (13). Thus the proposed algorithm is entitled by adaptive PHD (APHD). The likelihood function lk(zk∣xk|k-1(p)) in (21) is defined as (22)lkzk∣xk∣k-1p=∏i∈Ci(xk∣k-1p)∏j∈Cj(xk∣k-1p)lkzki,j∣xk∣k-1p,where Ci(xk∣k-1(p)) and Cj(xk∣k-1(p)) denote the covering area of the target in the image, which has been given in Figure 1. lk(zk(i,j)∣xk|k-1(p)) is the single pixel likelihood, which is expressed via (14) and (15). Consider(23)lkzki,j∣xk∣k-1p=f(zk(i,j)∣H1)f(zk(i,j)∣H0)=exp⁡-a21F1N2;J2;ar2.The expected number of targets is calculated by summing up the total weights; that is,(24)N^k=∑p=1Lk-1+Jkωk∗(p).Then resampling is performed to avoid the degeneracy. The updated posterior density at time stepk is presented as (25)Dk(xk∣Z1:k)=∑p=1Lkωkpδxk-xkp.In summary, the processing chain of the proposed approach is presented in Figure2.Figure 2
The processing chain of the proposed approach.


## 5. Simulation and Result Analysis

To validate the effectiveness of this algorithm, it is employed to detect and track three typical exoatmospheric targets. Warhead, fragment, and decoy are introduced in the multispectral infrared imagery sequence, which consist of 3 wavebands. The three infrared sensors include the long wave IR, mid-long wave IR, and short wave IR. Each waveband has 30 frames. And each frame size is180×240 pixels. Frame 1 from the image sequence at wavelength 4 μm is shown in Figure 3, which is the actual background with clutter noise.Figure 3
Frame 1 from the image sequence at wavelength 4μm.Warhead is introduced at position (25.6, 24.8) in frame 1 and remains within the surveillance region until the last time step. Fragment appears at position (28.3, 47.2) in frame 6 and died in frame 21. Decoy appears at position (56.7, 27.4) in frame 11 and disappeared in frame 26. This simulation does not include any spawning of new targets from existing ones. The IR intensity sequences of the three targets are simulated according to the BODE (Bhattacharyya Optical Decoy Evaluation) model [21] and the basic physics parameters for space objects [22]. The IR intensity has been normalized, which is shown in Figure 4.Figure 4
The intensity history for three types of exoatmospheric targets.The exoatmospheric target state can be approximately modeled as a linear Gaussian process. Therefore (3) becomes(26)xk=Fxk-1+Gnk-1,where the process noise is nk-1~N(0,0.052). F is the state transition matrix, and G denotes the noise matrix; they are defined as (27)F=1T000100001T0001,G=T220T00T220T.The sampling period T=1 s is adopted. Each target is subject to the same velocity (0.5, 0.5) pixel/frame. The initial prior probability density function (PDF) of velocity is assumed to be uniform within [-1,1] pixel/frame. It is known according to above parameters that the three exoatmospheric targets appear and move within the image region [1,100]×[1,100]. To observe the particles’ distribution conveniently, the image size for displaying is cut to 100×100 pixels. Two simulated image frames are shown in Figure 5. Some targets are surrounded by the actual dense background clutter noise.Figure 5
Two image sequences: frame 1 and frame 8.
(a)
Frame 1
(b)
Frame 8The procedure of implementing the original PHD filter based on particle filtering can refer to [20]. Here, we specify the main parameters for implementing the adaptive PHD filtering in Table 1. The number of particles representing one target Lk=2000, the number of particles representing new born targets Jk=3000, the probability of target survival ps,k=0.98, the probability of target spawning αk=0, and the probability of spontaneous target birth bk=0.01. The probability of detection PD is calculated according to (13).Table 1
The main parameters for implementing the adaptive PHD filtering.
Variable
L
k
J
k
p
s
,
k
α
k
b
k
P
D
Value
2000
3000
0.98
0
0.01
Equation (13)The raw infrared image sequence cannot be directly used to detect the targets because of background clutter. The adaptive clutter suppression filter can be used to remove the background clutter. Figure6(a) shows the target appearance curve which obviously has a nonstationary mean before the suppression of clutter. Figure 6(b) depicts an adaptive suppression curve in which the pixels distribute around zero. Comparing Figure 6(a) with Figure 6(b), the adaptive clutter suppression algorithm removed the background correlation. Thus the residual clutter/noise is approximately subject to the Gaussian distribution. It is noted that since the Gaussian model is the sufficient condition to implement the RX algorithm, the measurements produced by the clutter suppression filter can be utilized in the RX algorithm. The experimental results demonstrated that the proposed background clutter suppression algorithm performed very well.Figure 6
Adaptive background clutter suppression.
(a)
Before removing the clutter
(b)
After removing the clutterThe multispectral images are fused by using RX algorithm in each frame. In order to show the potential of multispectral fusion, the single band image and the three band images are processed, respectively. Compared with the single wavelength image, the three band images have good detection output, which was shown in Figure7(b). In contrast, the target was immersed in the single band image, which was shown in Figure 7(a). The bright spots in Figure 7(b) give a priori positions of the aiming target. The knowledge of data can be used to increase the efficiency of the particle filter. Thus the efficient proposal density function can be designed by the bright pixels. The 50 highest pixels in the image are used to design the proposal density function. The proposal density which includes the latest measurements improves sample degeneration and impoverishment. The experiment shows the effectiveness of the spectral image data fusion, which also reduced the rate of false alarm and clutter density.Figure 7
Multispectral fusion output using RX algorithm.
(a)
The multispectral fusion output (L=1)
(b)
The multispectral fusion output (L=3)The proposed APHD algorithm is performed on the spectral fusion output to detect and track multiple targets. The impact of the number of spectra on the performance of the proposed algorithm is considered to show the potential of multispectral fusion. In addition, the robustness of the proposed algorithm for low SNR targets is investigated. Two scenarios are constructed by selecting different SNR. Actually, low SNR targets are defined as the targets for which the sensor responses have a value of SNR lower than 10 dB in target detection and tracking applications [23]. Thus, two typical SNR values of 8.7 dB and 4.8 dB are considered in our work. The definition of SNR refers to [24]. The performance of this algorithm is evaluated by using the averaged optimal subpattern assignment (OSPA) [25], which is a good metric for multitarget detection and tracking. The cut-off parameter c and the order parameter p should be chosen in OSPA. To yield smooth distance curves, we have c=40 and p=2 in this paper. One hundred Monte Carlo trials are performed for each scenario. Figures 8(a) and 8(b) show the estimated number of targets and estimated error, respectively, when SNR = 8.7 dB. Figures 9(a) and 9(b) show the results distinctively when SNR = 4.8 dB.Figure 8
The estimated number of targets and OSPA (SNR = 8.7 dB).
(a)
Estimated number of targets
(b)
Curves of OSPAFigure 9
The estimated number of targets and OSPA (SNR = 4.8 dB).
(a)
Estimated number of targets
(b)
Curves of OSPAUnder the condition that SNR = 8.7 dB, the estimated number of targets is accurate and this algorithm performs very well with no delay in track termination after the targets disappeared when multiple spectrums are used. Although the output ofL=3 is not correct at frames 1, 17, and 28 in Figure 8(a), the continuing estimation of number of targets is still stable. However, the estimated number of targets of L=1 has a bias, concluding from the curve of OSPA in Figure 8(b).Under the condition that SNR = 4.8 dB, the estimated number of targets has a little bias and delay appears in track termination after the target disappeared when multiple spectrums are used, which is shown in Figure9(a). However, the estimated number of targets of L=1 has a big bias and long delay when the targets disappeared, which can be seen from the curve of OSPA in Figure 9(b). In the circumstances, it is deduced that multiple spectral fusion improved the detection and tracking performance.Figure10 shows the estimated error of the number of targets, which is calculated by T-T^, where T is the true number of targets and T^ denotes the estimated number of targets. From Figures 10(a) and 10(b) one can observe that the multispectral fusion (L=3) provides significantly higher accuracies than single waveband spectral (L=1) when estimating the number of targets.Figure 10
The estimated error of number of targets at different SNR.
(a)
The estimated error (SNR = 8.7 dB)
(b)
The estimated error (SNR = 4.8 dB)In order to show the benefit of the proposed method, we contrast the performance of the proposed approach (APHD) with the SPHD filtering in the same scenario as above (SNR = 8.7 dB,L=3). The detailed description of SPHD filtering can refer to [17]. Since the lack of memory of the SPHD filtering may result in premature death of targets, we used a spectrum with low probability of detection at the start of the iterated update rather than at the end. The probability of detection can be calculated according to the RX algorithm.We compared the performance of APHD filtering and SPHD filtering in terms of the estimated number of targets and OSPA in Figure11. It is clear in Figure 11(a) that the proposed algorithm is more stable than the SPHD filtering. The OSPA distance is adopted to compare APHD filtering with SPHD filtering in Figure 11(b), which shows that the SPHD filtering has more fluctuations due to the estimation of target number. As a matter of fact, it is difficult to decide the order of the multisensor update owing to the probability of detection. Even though we obtain the correct order of the multisensor update for the SPHD filtering in advance, the presented APHD filtering has better performances than the SPHD filtering. Thus it is concluded that the proposed APHD filtering outperforms the SPHD filtering. The corresponding results for SNR = 4.8 dB are not shown because of the poor performance of the SPHD filtering in similar cases.Figure 11
The performance comparison of APHD filtering and SPHD filtering (SNR = 8.7 dB).
(a)
Estimated number of targets
(b)
Curves of OSPA


## 6. Discussions and Future Work

The fusion of multispectral image sequence improved the detection and tracking performance in multiple targets scenario, especially when the targets had low SNR; that is, SNR < 10 dB (Figures8 and 9). We found that even though the SNR was 4.8 dB, the detection and tracking for multitarget was still robust (Figures 9(a) and 9(b)). We did observe that the higher the SNR was, the better the detection and tracking performance became (Figures 8, 9, and 10). Compared to the broadband image sequence, multispectral data fusion produced more information of the detected targets (Figures 7(a) and 7(b)).For the issues arise from multitarget tracking, previous studies have found that the PHD filtering has superiority. However, the derivation of PHD filtering was based on the single sensor measurement model. We first derived the likelihood function via multispectral data fusion (23) and used this likelihood to update the posterior density. This procedure is centralized data fusion and different from the sequential PHD filtering which is distributed data fusion. The previous studies have shown that centralized data fusion is better than distributed data fusion. Moreover, the probability of detection and the probability of false alarm are always fixed in the conventional PHD filtering. These two parameters usually have variability in the real world. We presented the calculated probability of detection and the calculated probability of false alarm via measured data ((13), (14), and (15)), which replace the two parameters in the conventional PHD filtering. Thus the posterior density in this paper is adaptive to the coming measurements. Additionally, the comparison of the proposed method with the conventional PHD filtering is given in the paper. It showed that the proposed adaptive PHD filtering outperforms the conventional sequential PHD filtering (Figures 11(a) and 11(b)).Though multispectral data fusion and PHD filtering have been successfully used to study the detection and tracking of multitargets, some limitations and potential work in future should be stated. Firstly, the RX algorithm used whitening data. So the correlation between the background pixels should be removed as much as possible. This paper proposed the linear least square algorithm. The nonlinear method removing background can be studied in the future. Secondly, when the target is far from the sensor, that is, the target has subpixel image, the multispectral fusion method based linear spectral mixing model should be considered. Finally, the particle filtering is a framework to implement the PHD filtering. But the original particle filtering has many issues, such as the degeneracy and diversity of particles. Thus other variants of particle filtering can be proposed to improve the PHD filtering in the future.


## 7. Conclusions

In this paper we have considered the multitarget detection and tracking problem arises from IR multispectral image sequences system. We have proposed a fusion tracking approach utilizing the spectral and temporal information of IR multispectral image sequences. The fusion tracking method consists of multispectral data fusion and adaptive probability hypothesis density (APHD) filtering algorithm. The influence of the number of spectrums and the SNR values on the system performance has also been studied. The experimental results have shown that the presented fusion tracking method has better detection and tracking performance than the conventional algorithm, even for targets with low SNR.



---

*Source: 179039-2015-04-20.xml*