Abstract
When tracking multiple targets, radar measurements from weak targets are often masked by the ambiguity function (AF) sidelobes of the measurements from stronger targets. This results in deteriorated tracking performance and lost tracks. In this study, we consider the design of configurable waveforms whose AF sidelobes can be positioned to unmask weak targets. Specifically, we construct multicarrier phasecoded (MCPC) waveforms based on Björck constant amplitude zeroautocorrelation (CAZAC) sequences. The MCPC CAZAC waveforms exhibit wide regions in their AF surface without sidelobes and allow for selective positioning of sidelobes. We apply these waveforms in the context of a target tracker by selecting waveform parameters that minimize the expected tracking error. We show that this is accomplished by selecting the position of AF sidelobes to unmask weak targets. The target tracker is based on an independent partitions likelihood particle filter that is capable of processing the highresolution measurements resulting from the Björck CAZAC sequences and tracks a fixed and known number of targets. Using simulations, we demonstrate the improvement in tracking performance when we adaptively select the MCPC CAZAC waveforms over tracking using nonadaptive waveform configurations or singlecarrier phasecoded CAZAC waveforms.
Introduction
When tracking multiple targets using radar sensors, weak targets are often difficult to observe in the presence of strong targets. This is because the ambiguity function (AF) sidelobes of measurements from strong targets are higher than the AF mainlobe of measurements originating from weak targets. As a result, the joint tracking performance of a multitarget tracker, expressed either in terms of meansquared error (MSE) or percentage of lost tracks, is poor. The location and magnitude of the AF measurement sidelobes in the delayDoppler plane are directly related to the location and magnitude of the AF sidelobes of the transmitted waveform. The AF is in turn defined by the type of transmitted signal and its parameters [13]. Therefore, there is a need to design configurable radar waveforms and develop an adaptive radar sensor configuration technique to position sidelobes from strong target returns away from the predicted locations of weak targets.
In [2,3], the processing of the return signal was performed by partitioning the delayDoppler plane into resolution cells with fixed locations. These cells were constructed in such a way as to approximate a probability of detection contour that depends on the signal type and its parameters, which were assumed to be fixed. A detection in a resolution cell was declared based on the thresholded output of a matched filter placed on the centroid of the cell. However, the shape of the probability of detection contour is often not well approximated by a tessellating shape, resulting in measurement errors. Instead of using a fixed waveform, adaptive waveform techniques were used to minimize either the tracking error or validation gate volume in [4,5]. Moreover, in [6,7] waveform parameter adaptation was used to minimize track loss in the presence of clutter. In [8], the probability of track loss and a function of estimation error covariance was minimized by selecting both waveform parameters and detection thresholds for range and range rate tracking in clutter. In [9], the time to detect new targets was minimized by posing the problem as a partially observed Markov decision problem. In addition, in [10], the one step ahead expected information from the target kinematic model was maximized using the appropriate waveform selection.
The methods mentioned above rely on linear observation models that do not accurately represent physical systems with nonlinear characteristics. In [11,12], iterative adaptive waveform techniques were developed for nonlinear system models with a single target and using frequencymodulated waveforms. These adaptive waveform techniques assume that the measurements are processed using matched filters on a fixed grid.
In this study, we develop a tracker that selectively collects measurements based on the predicted target state (instead of a fixed grid for linear systems) to track a known and fixed number of targets. We implement the tracker using a new method that does not require the collection of measurements exhaustively on a fixed grid in the AF plane. The new independent partitions likelihood particle filter (IPLPF) tracker adaptively configures multicarrier phasecoded (MCPC) waveforms [1] so that their AF sidelobes can be positioned in such a way as to not mask weak targets in the presence of strong targets. We also develop an adaptive configuration strategy to select the MCPC waveform parameters based on the relative positioning of the targets in the delayDoppler plane. In contrast to previous methods [11,12], the proposed adaptive waveform selection method is not iterative; in contrast, waveform parameters are directly selected based on the state information on the weak target relative to the strong targets. The new IPLPF tracker uses a proposal distribution that is based on the independent partitions algorithm [13,14] and the likelihood particle filter [15]. The particle filtering framework accommodates the propagation of resolution cells to locations of interest in the delayDoppler plane based on prior information on the target state. It can also use the exact shape of the probability of detection contour as a resolution cell instead of the simplified tessellating shape. With this approach, we do not need to exhaustively collect measurements on all points of a fixed grid. Instead, the matched filter is matched to locations that are likely, given the belief about the target state; these locations are represented by each of the particles of the particle filter.
We employ MCPC waveforms in our work as their AFs can exhibit both wide regions with zero magnitude as well as nonzero sidelobes that can appropriately be positioned based on how the waveform parameter values are chosen. In order to construct these MCPC waveforms, we use multicarrier modulation and equallength Björck constant amplitude zeroautocorrelation code (CAZAC) sequences [1618] that are cyclicshifted versions of one another. A Björck CAZAC sequence provides higherresolution measurements than a linearly frequencymodulated chirp [2,3] due to its highly concentrated AF. The high concentration in the AF plane results in improved tracking performance as we demonstrated in [19] for the single target case. For the multitarget case, measurements from singlecarrier phasecoded (SCPC) CAZAC sequences exhibit sidelobes that are spread in the delayDoppler plane and could mask weak targets. Our proposed use of configurable MCPC CAZACs, on the other hand, can adaptively position the waveform sidelobes to not mask weak targets.
We configure the MCPC CAZAC parameters at each time step of the tracking scenario to minimize the predicted tracking error. The waveform parameters selected are the cyclic shift in frequency used to generate the waveform, the number of CAZAC sequences used, and the length of the CAZAC sequences. We present a computationally feasible method of selecting parameters to position sidelobes that accounts for the predicted relative positioning of the targets. Furthermore, our simulations demonstrate that the minimization of the predicted tracking error, achieved by selecting the MCPC waveform parameters, can be achieved by positioning the AF sidelobes such that a weak target is not masked in the presence of strong targets.
The rest of the article is organized as follows. In the following section, we present the MCPC CAZAC waveforms and investigate the properties of their AFs. In Section “IPLPF algorithm”, we provide a detailed description of the IPLPF algorithm and its application in minimizing the predicted tracking error. In Section “Adaptive waveform selection”, we integrate the IPLPF with a waveform configuration algorithm and demonstrate its performance in tracking multiple targets in Section “Simulation results”.
MCPC CAZAC sequences
Björck CAZAC sequences
A CAZAC sequence ξ(m) with finite length M, has constant magnitude, ξ(m) = 1, m = 0,…,M − 1, and zeroautocorrelation, , for n ≠ 0, where the addition is modulo M[17,20]. An example of a CAZAC sequence with quadratic phase is the Björck CAZAC sequence. For prime length M = 1, mod 4, it is given by [16,17]
where m, mod M (or m modulo M) is the remainder of the division m/M, and [ (m/M)] is the Legendre symbol that is given by
Björck CAZAC sequences are an attractive choice for target tracking with radar [17] as their constant amplitude allows for continuous transmission of peak power and can thus lead to increases in signaltonoise ratio (SNR). They also exhibit very tight localization in the delayDoppler plane that can enhance the range resolution and rangerate resolution of the measurements. The discrete AF of a Björck CAZAC sequence is given by [21]
where n and νare the discrete delay and Doppler parameters, respectively. Specifically, the AF exhibits a large spike at the origin (nν)=(0,0) of the discrete delayDoppler plane, with very small sidelobes. An example of the AF of a Björck CAZAC of length M = 1,741 is shown in Figure 1.
Figure 1. AF surface plot of a Björck CAZAC of length M = 1,741.
MCPC Björck CAZAC sequences
As the AF of a Björck CAZAC sequence is very highly localized, we want to exploit its properties to position AF sidelobes for multiple target tracking. In particular, we use the fact that a cyclic frequencyshifted CAZAC is also a CAZAC [20] and also that a sum of cyclic frequencyshifted sequences has an AF surface whose sidelobe locations depend on the difference in cyclic frequency shift, number of sequences, and sequence length. Note that, although cyclic permutations of CAZACs are possible in both time and frequency, we restrict our attention to frequency shifts as they result in wide zero regions in the AF plane and can better facilitate the adaptive positioning of the AF sidelobes.
The MCPC scheme combines multiple waveforms that are modulated by orthogonal carriers; the carriers are separated in frequency using orthogonal frequency division multiplexing [1]. The phase coding is required to reduce the bandwidth of the CAZAC sequence so that it meets transmission requirements. We use this scheme to form the MCPC CAZAC waveform by combining Q cyclically permuted Björck CAZAC sequences. Specifically, if we cyclic frequency shift ξ(m) in (1) using the frequency shift ζ_{q}(assumed mod M) to obtain the qth SCPC cyclic frequencyshifted CAZAC waveform
then the MCPC CAZAC waveform, modulated with carrier frequency ζ_{c}, is given by
where m = 0,…,MQ−1, and ⌊·⌋ denotes rounding down to the nearest integer. Note that we restrict ζ_{q} = qζ (mod M) in (3) as this selection of cyclic frequency shift causes the positioning of the sidelobes of the AF to depend on ζ, thus facilitating adaptive waveform configuration in our proposed algorithm. Thus, Θ = (QMζ) in (4) defines the three parameters of the MCPC CAZAC waveform.
When processing the SCPC CAZAC waveform in (3) using the AF in (2), the narrowband assumption is used that states that the transmitted waveform does not experience any time scale changes due to target motion. This assumption is valid since the timebandwidth product of the waveform can be shown to be much less than as the speed of propagation in the air c is large, where is the target range rate ( [22], p. 241). As we restrict the timebandwidth product of an SCPC sequence and an MCPC sequence to be the same, the narrowband assumption also holds for MCPC CAZACs. Note also that we double the number of possible AFs by taking the Fourier transform (FT) of each of the MCPC CAZAC waveforms that we construct. The AF of the transformed waveform is equal to the AF of the original waveform with the delay and Doppler variables interchanged. This offers a convenient method of producing additional sidelobe positioning options with little effort.
AF surface of MCPC CAZAC waveforms
The AF surface of the unmodulated MCPC CAZAC waveform in (4) is given by . Using (2), (3), and (4), the AF is given by
where is the energy of s_{Θ}(m) that is normalized to have the same energy as ξ(m). Next, we consider two separate cases of cyclic frequency shifts: ζ=0 and ζ > 0.
Zero cyclic frequencyshift
When ζ = 0, two of the exponential terms in (5) cancel out. We can also simplify the summations and where and are integers (see [23] for the derivation details). The resulting AF of the MCPC CAZAC with Θ = (QM,0) becomes
We can see that the AF in (6) is nonzero only if is a multiple of Q, thus resulting in zero AF surface regions of width Q. Although these regions can be used to reveal weak targets at selected areas in the AF plane, we also need to reduce the sidelobes near the origin of the AF. The area in the delayDoppler measurement plane near the AF origin is the area that is most commonly interrogated by the IPLPF tracker when accurately tracking a target, as we will see in Section “IPLPF algorithm”. Since we already have zero AF surface regions in the interval n = 1,…,Q − 1, we need to investigate the shape of the AF surface along the Doppler axis ν at n = 0. Setting in (6), we obtain the AF surface as
Since ξ(m) = 1 for all m, we conclude that the AF surface is nonzero only when νis an integer multiple of M. Therefore, the location of the sidelobes when n = 0 can also be chosen by adjusting the value of M. An example of this is shown in Figure 2 that depicts the AF surface of an MCPC CAZAC waveform with Θ = (130,13,0); all nonzero sidelobes exist when n is an integer multiple of Q = 130.
Figure 2. AF surface of an MCPC Björck CAZAC with parameters: (a) Θ = (Q,M, ζ) = (130,13,0) and (b) Θ = (130,13,1).
Positive cyclic frequencyshift
When ζ > 0, we obtain higher diversity in the locations of the AF sidelobes. Specifically, when ζ ≠ 0 in (5), we observe that the terms e^{j2Π(⌊(m − n) / Q)⌋)qζ/M}and repeat multiple times in the summation in Equation (5). This is due to the summation of these terms over q = 1,…,Q − 1 and the modulo M / ζ effect of the two exponential functions which is due to M / ζ < Q. Therefore, we can factor the repeating terms and rewrite (5) as
Note that q and now vary from 0 to ⌊Q/β⌋ − 1. We choose Q, M and ζ such that β = ⌊(M−1)/ζ⌋ + 1 is approximately a multiple of Qfor most choices of ζ = 1,…,M − 1. This eliminates the summation terms for q and that fall between the values of (β(⌊Q/β⌋−1) + β−1) and Q − 1. These terms were omitted in (7), which explains the use of the approximation symbol, as they only cause a negligible variation of sidelobes in the AF surface compared to the exact expression. The accuracy of the above approximation can be verified using a numericalbased analysis, i.e., the generation of the AF surface using the Matlab code used in this work which is available to the reader upon request. The summation with respect to can be simplified using , where is an integer. We then let , where and . Also, , where is an integer. We let with an integer and . This simplifies Equation (7) to
This expression shows that nonzero values of the AF exist for integer and , i.e., for integer . This provides for controlled size valleys in the AF surface.
We also examine what happens along the Doppler axis ν at zero delay and ζ > 0. If we set n = 0 or in (8), then the term reveals that the AF surface has sidelobes that periodically repeat with period βM. Evaluating the above expression at the inbetween intervals, we can obtain the AF surface sidelobe values. We then choose to use only waveforms with parameters Q, M, and ζ with relatively low sidelobe levels in their AF surface with n = 0.
When ζ = 1, β = M, larger valleys appear in the AF surface. Specifically, the AF in (8) becomes:
Using since , and, therefore, having above expression becomes
Then we note that the factor can only take the values of 0 if and −1 if since both m and take values less than β. This implies that nonzero values of the AF surface only exist at delay locations such that or are multiples of β. For which restricts to be a multiple of βwith , and since , M<Q then indices of the waveform in the AF expression summation are very limited compared to the waveforms’ length MQ (i.e., m<M). Therefore, the case where in the AF expression appears in a very small number of additive terms and is omitted in the following analysis. Letting , and using , we obtain (for details, see [23])
Since this implies that , thus with β = M, it follows that nonzero sidelobes of the AF surface appear at intervals of Q + (Q/M) in the delay. This is demonstrated in the AF surface of the MCPC CAZAC waveform with Θ = (130,13,1), as shown in Figure 2b.
In summary, the possibility of choosing the parameters Θ = (Q,M,ζ) of an MCPC CAZAC waveform, and also rotating the entire AF surface by choosing to take the FT of the waveform, enables us to position sidelobes in order to minimize the predicted MSE, as we will show in Section “Adaptive waveform selection”.
IPLPF algorithm
Tracking model
We consider L targets moving in a twodimensional (2D) plane, where the number of targets is fixed and known. The target dynamics are modeled by a linear, constant velocity model [24] given by
where is the state vector for the lth target at time k, T denotes vector transpose, x_{l,k}, y_{l,k} and , are the position and velocity in Cartesian coordinates, respectively, the matrix F is given by F =[1 Δt 0 0;0 1 0 0;0 0 1 Δt;0 0 0 1] (with each row in square brackets), Δt is the time difference between observations, and v_{l,k} is a zeromean, additive white Gaussian process with diagonal covariance matrix that models target deviations from constant velocity. The model in (9) can be used to determine the kinematic prior probability distribution function for the lth target. The multitarget state vector is expressed in terms of the state vectors of each target as . Following [13], we refer to each component x_{l,k} of X_{k} as a partition. Since we assume that the targets move independently, the multitarget kinematic prior distribution is given by .
A radar sensor collects information on the range and range rate of the targets in the scene relative to the sensor by transmitting pulses and processing the returns after they are reflected by the targets. The return waveform provides range information, in the form of time delays, and range rate information, in the form of frequency shifts of the return waveforms relative to the transmitted waveform. Assuming point targets, the range and range rate for partition l at time step k, relative to the uth sensor, u = 1,…,U are given, respectively, by [11] and , where are the Cartesian coordinates of the location of the uth sensor, and the U sensors are assumed to transmit and receive waveforms independently. The discrete time shift value and discrete Doppler shift value at the u sensor, due to the lth target, are given by [11] and , respectively, where round (·) transforms the real number to the nearest integer, c is the velocity of propagation in the medium, f_{c} is the carrier frequency, T_{s} is the sampling period, and M is the total number of waveform samples.
Matched filter statistic
At every time step k, a signal s(m), m = 0,…,M − 1, is simultaneously transmitted from each sensor in different frequency bands to avoid interference. The received signal at the uth sensor after demodulation is a linear combination of the reflections from all L targets, and it is given by
Here, ζ_{c} = f_{c}T_{s} is the discrete carrier frequency of the transmitted waveform. The sum of random complex returns, A_{l,k}, from many different target scatterers on target l are zeromean, complex Gaussian with known variance and follow the Swerling I model [25]. Each target is assumed to have a different radar cross section (RCS) [26] and thus its return signal has a different strength that is represented by the variance of A_{l,k}. It is also assumed that the return signal strength depends only on the target RCS and not on the distance between the sensor and the target; the distance is compensated for by amplifying returns that arrive later in time. The noise terms v_{u,k}(m), u = 1,…,U, are assumed to be zeromean complex Gaussian with variance 2 N_{0} and independent for each sensor. The SNR is given by [2], where l_{w} is the index of the weakest target and E_{s} is the energy of the transmitted waveform from each sensor. When no target is present, d_{u,k}(m) = v_{u,k}(m).
At the receiver, the return signal is matched filtered with a signal representing returns from Θ targets, at different time shifts and different frequency shifts , λ = 1,…,Θ. These timefrequency shifts are derived from the belief in target state using a particle filtering approach. The matched filter output is thus given by
where M_{d} > M should be large enough to accommodate a maximum delay in the signal due to a reflection from the target. Moreover, Θ in (10) equals 1 when independently proposing one partition and Θ = L if particles of L partitions are proposed. The matched filter statistic that we will use for estimation is , and it is written in terms of the AF of s(m) in Equation (2).
Measurement likelihood
The statistical properties of the matched filter statistic y_{u,k} depend on the fact that A_{l,k} and v_{u,k}(m) are independent, zeromean, and complex Gaussian. It can be shown that in (10) is also complex Gaussian with zero mean and y_{u,k} is exponentially distributed both under hypothesis H_{0} (not target is assumed present) and under hypothesis H_{1}(L targets are assumed present). The two hypothesis formulation for the measurement likelihood is thus given by
where (see [23] for derivation)
When using the MCPC CAZAC waveforms to compute the measurement likelihoods in (11), we can reduce the computational complexity by approximating the above variance expressions. In particular, since the AF sidelobes of MCPC CAZAC waveforms are zero at the locations where, according to the belief in target state, the targets are expected to be we can set AF_{s}(n,ν) to be 0 for n ≠ 0 and ν ≠ 0 in the above expressions. In addition, for the SCPC waveforms the above holds only approximately due to their nonzero, however, very low AF surface sidelobes. Also, using the fact that AF_{s}(0,0) = 1, we can let for all l (where is a nominal value that we choose, since we assume the target strength to be unknown). Based on this, we can then set and in (12).
Likelihood partition sampling
The highly concentrated AF of a Björck CAZAC sequence provides a highly concentrated likelihood proposal distribution and a high measurement accuracy. However, the proposal process needs to be modified to sample particles from the likelihood instead of the kinematic prior since the former is much more localized than the latter. To achieve this, we use a likelihood particle filter [15], where the importance density depends on the measurements rather than the kinematic prior.
We propose to integrate the use of the likelihood proposal with the independent partition (IP) particle filtering [13,14] concept to efficiently propose particles. In the IP, we propose individual partitions of the multitarget state vector, each representing the state of a single target. We then combine the more accurate partition proposals into particles. The IP algorithm is an approximation to the joint multitarget probability density particle filter [14]; the approximation is accurate when the targets are well separated in the observation space. When targets are close in measurement space, their partitions cannot be independently proposed as described above. Due to our use of the Björck CAZAC sequences that have sharply peaked AFs, the measurements can be well approximated as independent [27]. Our resulting algorithm, the IPLPF, belongs to the class of sequential partition algorithms [28]. Algorithms of this class propose partitions sequentially and then combine them into particles.
Specifically, we first independently evaluate likelihood values at discrete delayDoppler bins for each partition. Using these values, we create histograms and sample partition states. Note that we narrow our bin selection to a region of probability of almost one in the kinematic prior partition sample. This is necessary to ensure a minimum number of bins to build the histogram and to ensure that the sample from the measurements is consistent with the kinematic prior. We then evaluate partition weights by combining measurements from the different sensors using the kinematic prior. Using the normalized partition weights, we independently sample values for each partition. We combine the sampled partitions into particles, compute the weights of the particles, and estimate and resample the particles.
Note that range information from two sensors, combined with kinematic prior knowledge, can provide adequate information to estimate the position of a single target [19]. Specifically, geometry is used to find the intersection between two circles whose radii (in Cartesian coordinates) are the ranges of the sensors. The two intersections of these circles provide two coordinate locations, one of which can be selected that agrees with the kinematic prior information. For multiple targets, there are multiple of these circles for each sensor and each target, and multiple intersection points that do not correspond to true coordinate locations. Therefore, the use of three sensors can help clear the ambiguity by providing fewer intersections of three circles. In order to avoid complicated geometry, we first process the returns of two sensors and sample Cartesian coordinate target locations using the likelihood. We then weight the sampled locations with measurements from the third sensor. Our method is for a general number of sensors equal or greater than three; in this work, we kept the number of sensors to a minimum of three.
The partitions sampling based on the likelihood is performed in two stages. In Stage 1, we utilize information from only two of the sensors in order to propose a preliminary set of partitions. This avoids the complex geometry required to sample Cartesian locations from range and range rate information obtained from three or more sensors. In Stage 2, we refine our partitions selection by sampling from the preliminary set of partitions created in the first stage using information from all the sensors.
Stage 1: partitions sampling
We start by propagating each state partition without noise. We let λdenote the proposed partition at time k and l denote the partition that represents the true state of the lth target. Assuming that we use a sequential importance resampling particle filter [15], the ith state particle, i=1,…,N, is given by . Using , we obtain
We want to determine a region of delayDoppler bins that could contain observations if the true state is . This region is obtained from the spread of the kinematic prior that determines the possible states of partition λ. If we assume for simplicity that the variances and of the kinematic prior in the 2D (x,y) dimensions are equal, then with probability of almost one, the proposed particle will fall within 3σ_{x} from and within 3σ_{y} from . The maximum and minimum possible sampled x and y coordinates will then yield the maximum and minimum range. That is, if we assume that the target is at angle Π/2 with the sensor, then . The range would then increase/decrease by the amount . The delay would also assume minimum and maximum index values given by
Similarly, the minimum and maximum values can be obtained for the range rate and thus the Doppler
We form all combinations of indices for delay and Doppler that lie within the minimum and maximum delay and Doppler values using
where and . Then, we evaluate the matched filter output at each of these values and for sensor u,
Note that we have used only one delayDoppler pair in the template signal representing a single partition λ. The single partition likelihood for this delayDoppler bin is given by:
We evaluate the likelihood ratio for each delayDoppler bin as
We then obtain
and the normalized distribution
from which we sample , , j_{n} = , , for each particle i and each sensor u = 1,2. The resulting sampled range and range rate and the bias are, respectively, , (2 f_{c}T_{s}), . The values of and , in turn, yield proposed state values . This is accomplished by taking the intersection of two circles in the 2D Cartesian plane and choosing the intersection point that mostly agrees with the kinematic prior information. This process is illustrated in Figure 3. We note that the sampled partitions from Stage 1 are based on information provided from only two sensors. Therefore, some of these partitions may be incorrect, as previously explained. However, the value of can be used in Stage 2 to evaluate likelihoods for three sensors in order to sample proposal partitions and help remove partitions that have incorrectly been sampled.
Figure 3. Schematic of the likelihood proposal process. Each particle is deterministically propagated forward (left top figure, arrows 1,2), the observation points for sensors 1 and 2 are defined (right top and right bottom), one point is sampled from each observation set of each sensor (arrows 3,4), and one of the two states formed by the sampled ranges in the Cartesian coordinates that agrees more with the prior is selected (left bottom).
Stage 2: partitions sampling
During Stage 1, we propose partitions , i = 1,…,N, from delayDoppler bins associated with sensors u = 1,2. In Stage 2, we utilize the return signals transmitted by U ≥ 3 sensors to refine our choice of partitions and to more accurately represent the target state. We first describe the complexity in calculating the partition weights and the approximation we use to make the computation tractable before we provide the details on the sampling process.
After matched filtering, and using the locations in the delayDoppler plane derived from the proposed partitions, the measurements from sensors u = 1,…,U are given by
where is the delayDoppler pair that corresponds to the state and the uth sensor, and is the true target state x_{l,k}. Therefore, the single partition likelihood function for each proposed partition λ of particle ɨ =1,…,N is given by . Here, the hypothesis of particle ɨ and partition λ is that the partition state equals and not for i ≠ ɨ, while , i = 1,…,N, u = 1,…,U are the measurements obtained from matched filters at the delayDoppler location defined by the particle proposed target state vectors , i = 1,…,N. However, each likelihood for sensor u is a multivariate exponential distribution [29] that grows in dimensionality as the number of particles N increases. We approximate the likelihood for each partition to be , where denotes the likelihood that a target exists at and denotes the likelihood that a target does not exist at . In [27], we show that the covariance between measurements and , ɨ ≠ i, depends on the filter proximity (i.e., the closeness of , and , ) relative to the AF spread. Therefore, the measurement independence approximation for the Björck CAZAC is reasonable due to its concentrated AF. Using this approximation, the weights for partition λ of particle ɨ = 1,…,Nare
If we divide the righthand side by the constant , and we use (17)–(19), we obtain
where the likelihood probability functions are given in (11) for a single target. This is then normalized
where
We finally perform partition resampling, where we sample a partition index , ɨ =1,…,N, from the distribution of with replacement. The resulting selected partition has value and selection probability .
Particle weighting
After partition resampling, we assemble particles from the sampled partitions as . We weigh these particles with weights that incorporate prior and measurement information. To find the weight equation, we start by defining a measurement matrix Y_{k} that is composed of measurements from and contains measurements from each of the U sensors. Specifically,
where
The likelihood function (for a single partition case) for each proposed particle ɨ =1,…,N is given by , i = 1,…,N, u = 1,…,U. Here, the hypothesis of particle ɨ is that the state equals , while are the measurements obtained from matched filters at the delayDoppler location defined by the particle proposed target state vectors . Note that this likelihood is a multivariate exponential distribution [29]. Using similar arguments as for the single partition case, we approximate the likelihood for each particle to be , where is the likelihood given that the target state equals and is the likelihood given that no targets exist having state .
The weight of particle ɨ [15] using the aforementioned assumptions is given by
Dividing by the constant , using the likelihood in (11), and normalizing the weights by , we obtain the normalized weighs
The state estimate is thus given by . The algorithm is outlined next.
IPLPF algorithm
■ For each partition λ = 1,…,Λ and for each particle i = 1,…,N
Stage 1: Likelihood Partition Sampling
✶ For each sensor u = 1,…,U
♦ Form using (13) and using (14)
♦ Evaluate using (15) and using (19)
Stage 2: Likelihood Partition Sampling
✶ For each particle ɨ =1,…,N
✻ Evaluate using (20) and using (22)
Particle Weighting
■ For each particle i = 1,…,N
✶ For each particle ɨ =1,…,N
■ Increment k by 1
Adaptive waveform selection
In order to further improve tracking performance, we adaptively select the parameters of the MCPC CAZAC transmit waveform at each time step k so that we can minimize the predicted tracking root meansquared error (RMSE). The three MCPC CAZAC parameters we consider are in (4), where Q_{k} is the number of cyclically permuted Björck CAZAC sequences at time k, M_{k} is the length of the sequences at time k, and ζ_{k} is a parameter that controls the cyclic frequency shift at time k. The expected RMSE is given by the cost function
where the weighting matrix Cmakes the units of the cost function consistent by compensating for the differing units of the state vector. The subscript in the expectation operator E_{·}[·] shows the dependance of the expected RMSE on the random target strength vector A_{k}, the random noise matrix v_{k}, the unknown true target state X_{k}, and the estimate , given the multitarget state estimate at k−1 and the choice of Θ_{k}.
Next, we identify the set of values that the multitarget state estimate can take in terms of the delayDoppler locations associated with it. As described in Section “IPLPF algorithm”, in order to propose particles, we have considered a discrete finite set of delayDoppler locations for each partition and each particle. This set corresponds to Cartesian coordinate locations that are most likely to occur according to the kinematic prior and the set of particles , i = 1,…,N generated at the previous time step k − 1. The set is given in (13) and (14) as , and , for λ = 1,…,Λ, u = 1,…,2 and i = 1,…,N. We use index ȷ to denote a member of the set , of cardinality , consisting of the N particles that could be sampled by the IPLPF proposal and subsequently weighted. Therefore, is a large set including all combinations of possible delayDoppler locations from two of the sensors for each target and particle. The process of forming partitions from sampled delayDoppler locations is explained in Section “Stage 1: partitions sampling” and illustrated in Figure 3. Subsequently, one possible outcome of the likelihood sampling process and particle weighting is N weightparticle pairs corresponding to delayDoppler locations , λ = 1,…,Λ, u = 1,…,U, i = 1,…,N. Similarly, based on the target motion model, we can identify a discrete finite set of possible true target states X_{k}. Each possible true target state with index is a member of the set , of cardinality . is related to corresponding delayDoppler locations , l = 1,…,L, u = 1,…,U.
From the above, we may rewrite the cost function in (26) as
where the probability distributions , and p(v_{k}) are defined in the context of the motion and measurement models in Sections “Tracking model” and “Matched filter statistic”.
In order to minimize the cost function, we need to minimize the probability and the particle weights in (25) for particles (i) such that . As the ȷth set of particles results from sampling by the IPLPF, we will follow the sampling process of the IPLPF and identify the selection probability for each partition of particles . According to Section “Stage 1: partitions sampling”, we obtain values for each partition λ = 1,…,Θ and each particle i = 1,…,N by sampling delayDoppler bins from sensors u = 1,2 with probability given by (19). The values allow us to evaluate the likelihoods for the U sensors in order to sample partitions. In Section “Stage 2: partitions sampling”, we obtain partitions with selection probability given by (22). These sampled partitions are combined into particles into particles in Section “Particle weighting”. From the sampling process of each particle , we conclude that the probability of each particle being selected is . Therefore, any set of particles appears with probability . Furthermore, from (21), (22), for and from (17), (19) for we observe that the above sampling probabilities depend on the single partition likelihood ratio which using (16) is proportional to . Since , the selection probability monotonically increases with the matched filter statistic . Therefore, in order to minimize the above selection probability the matched filter statistic needs to be minimized for the delayDoppler values in (13) and (14) with the additional constraint that for all partitions λ, particles i and sensors u. These sets of delayDoppler locations correspond to the belief on target state as explained previously and only include delayDoppler locations that imply erroneous target states (i.e., AF sidelobes). Since the matched filter statistic is a random variable it is minimized by minimizing its variance, given in (12) with Λ = 1, with respect to the waveform parameters.
Next, we observe that the particle weights in (25) contain the likelihood ratio both in the numerator and denominator. This, together with the fact that the prior has a wide spread compared to the likelihood, makes the particle weights nearly constant. Therefore, particle weights cannot be significantly reduced by adjusting the waveform parameters.
Therefore, the focus is on minimizing the matched filter statistic variance in (12) with respect to the waveform parameters specifically for the delayDoppler values in (13) and (14) and such that . Since the matched filter statistic variance depends on the AF of the waveform, and since the abovementioned set of delayDoppler values refer to AF locations where the target is expected to be at time step k, excluding the AF mainlobe, the problem of minimizing the RMSE reduces to the problem of reducing AF sidelobes in the area where the target is expected to exist. This is a very well defined area in the delayDoppler plane that is given by the sequential tracking process of the particle filter as explained in Section “IPLPF algorithm”. Configuring the waveform so that zero sidelobes appear in selected areas of the AF surface was described in Section “AF surface of MCPC CAZAC waveforms”. Therefore, at each time step of the scenario, the parameters of the waveform to be transmitted at the next time step are selected such as to achieve low AF sidelobes in the areas where the weak target is expected to be found, resulting to a minimization of the expected RMSE. This is computationally efficient compared to iterative methods of waveform parameter selection [11,12]. However, the entire multitarget particle filtering method proposed is still associated with a large computational load which is not expected to reach real time operation with stateofthe art hardware which also limits the number of targets that can be tracked.
It is noted that this method works well only if the number of weak targets is low. The AF surface valleys created by these waveforms are, as expected, of limited size since the uncertainty cannot be entirely eliminated. Therefore, if multiple weak targets happen to be relatively positioned such that AF surface valleys cannot be configured to contain them then these targets will be masked. The problem of unmasking a larger number of weak targets is, therefore, an open problem and a limitation of the proposed method. Moreover, there is a prediction error associated with each target location which is estimated based on the Bayesian methodology employed. In this study, the prediction error is minimized as targets are highly localized when using the high resolution, high AF surface peaked CAZACbased waveforms. The AF surface valleys designed are then large enough to contain this uncertainty and guarantee the unmasking of a weak target.
Furthermore, it is noted that the weak targets need to have a signal strength that is well above the noise level so that they are observable. In this study, it is assumed that what keeps weak targets masked are in fact the sidelobes from stronger measurement returns and not the noise. In order to initially detect weak targets, when no prior tracking information on their state is available, a sequential selective positioning of the sidelobes over different regions of the field of view is necessary. Once a weak target is detected and the tracking process begins then the selective positioning of the sidelobes based on prior tracking information described in this study is possible to take place.
Simulation results
We consider two simulation scenarios to demonstrate the performance of tracking multiple targets. The first scenario consists of one weak target and two strong targets; the second scenario consists of three targets of equal strength. Three different types of waveforms will be used: (a) SCPC Björck CAZAC, (b) MCPC Björck CAZAC with fixed parameters, and (c) MCPC Björck CAZAC with adaptively configured parameters. Three targets move in a 2D plane. The motion is completed in 199 time steps. Three sensors located at χ_{1} = −1000 m, ψ_{1} = 500 m, χ_{2} = 2500 m, ψ_{2} = 500 m, and χ_{3} = 500 m, ψ_{3} = 0 m collect range and range rate measurements. The trajectory of the target and the location of the sensors are shown in Figure 4. The target is assumed to move according to a nearly constant velocity model with covariance matrix Q = diag (225 64 225 64).
Figure 4. Target trajectory and sensor location.
In the first scenario, the weak target l = 2 has a crosssectional area such that, for SNR varying as 5, 10, 12, 15, 17, 20 dB, . The strong targets are characterized by . The noise variance is N_{0} = 1 and the waveform energy is E_{s} = 1. In the second scenario, all three targets are observed with SNR that varies as 5,10,12,15,17,20 dB. The SCPC Björck CAZAC waveform has length M = 1,741. The choice of parameters of the MCPC waveforms was limited to combinations of values {M,Q} = {7,245},{11,154},{13,130}, and ζ = 0,1 in order to reduce computational expense in the adaptive selection process. The FT of the waveform was also used to introduce another degree of freedom by rotating the AF. All waveforms are sampled at 8 MHz and frequency modulated by f_{c} = 40 GHz. The speed of propagation of the waveforms is c = 2.997925×10^{8} m/s. For the simulations, we used N = 300 particles, initialized by drawing from a Gaussian distribution with mean the true initial target position and covariance Q_{0} = diag (1000 100 1000 100). The results were averaged over 300 Monte Carlo runs. The parameters for the adaptively configured MCPC waveform are selected at each time step as described in Section “Adaptive waveform selection”, while for the fixed MCPC the parameters were selected randomly at the beginning of the scenario.
For the first scenario, the RMSE tracking performance is shown in Figure 5a for different values of SNR and for all waveforms. The percentage of lost tracks is shown for each waveform and SNR value in Figure 5b. A lost track is declared if, for 6 consecutive steps, the tracking error exceeded 300 m. We observe that the MCPC waveform with adaptive configuration (indicated as AMCPC in all figures) clearly outperformed the SCPC Björck CAZAC and fixed MCPC waveform when considering the number of lost tracks. In terms of the RMSE, the SCPC Björck CAZAC appears to have similar performance as the adaptive MCPC case since both have the same measurement resolution. Performing well in RMSE is, however, not useful if it is accompanied with a high number of lost tracks. The nonadaptive MCPC waveform case has the lowest performance rating due to its high sidelobes that are not avoided during measurement. When there are no successful tracks, the RMSE value is shown as zero in Figure 5a,b.
Figure 5. (a) RMSE versus SNR with 95% confidence intervals and (b) percentage of lost tracks versus SNR for three waveforms when tracking one weak and two strong targets.
The corresponding results from the second scenario are demonstrated in Figures 5a and 6a. We can observe that, if the targets have equal strength, then the adaptive MCPC and SCPC CAZAC perform similarly as their AF sidelobes do not mask weak targets. The large sidelobes of the nonadaptive MCPC, however, still result in large errors. Another observation is that in the first scenario, in the adaptive MCPC case, the results are improved compared to the second scenario. This is because in the first scenario, two of the targets have higher SNR values than in the second scenario.
Figure 6. (a) RMSE versus SNR with 95% confidence intervals and (b) percentage of lost tracks for three waveforms when tracking targets of equal strengths.
Conclusions
We developed the IPLPF algorithm, a particle filtering method based on the IPs approach and the likelihood particle filter, to track a fixed and known number of targets. A particle filter selects measurements based on the belief on the target state instead of collecting measurements exhaustively on a fixed grid. Moreover, the likelihood particle filter is capable of processing measurements resulting from the use of waveforms with highresolution properties such as Björck CAZACs. In addition, we developed MCPC waveforms whose AF sidelobes can be adaptively positioned. We outlined a configuration strategy for selecting the parameter values of MCPC waveforms to position AF sidelobes such that weak targets are unmasked by minimizing the predicted MSE. We demonstrated with simulations that when tracking targets with different strengths using single Björck CAZAC or fixed parameter MCPC waveforms results in deteriorated tracking performance. On the other hand, the use of adaptively configured MCPC waveforms enables the tracking of weak targets in the presence of strong targets and offers significant tracking performance improvements.
Competing interests
The authors declare that they have no competing interests.
Acknowledgements
This study was supported under MURI Grant No. AFOSR FA95500510443.
References

C Rago, P Willett, Y BarShalom, Detectiontracking performance with combined waveforms. IEEE Trans. Aerosp. Electron. Syst 34(2), 612–624 (1998). Publisher Full Text

R Niu, P Willett, Y BarShalom, Tracking considerations in selection of radar waveform for range and rangerate measurements. IEEE Trans. Aerosp. Electron. Syst 38(2), 467–487 (2002). Publisher Full Text

DJ Kershaw, RJ Evans, Optimal waveform selection for tracking systems. IEEE Trans. Inf. Theory 40(5), 1536–1550 (1994). Publisher Full Text

DJ Kershaw, RJ Evans, Waveform selective probabilistic data association. IEEE Trans. Aerosp. Electron. Syst 33(4), 1180–1188 (1997)

SM Hong, RJ Evans, HS Shin, Control of waveforms and detection thresholds for optimal target tracking in clutter. IEEE Conference on Decision and Control, vol. 4 ((2000), pp), . 3906–3907

SM Hong, RJ Evans, HS Shin, Optimization of waveform and detection threshold for target tracking in clutter. Proceedings of the 40th SICE Annual Conference ((2005), pp), . 42–47

SM Hong, RJ Evans, HS Shin, Optimization of waveform and detection threshold for range and rangerate tracking in clutter. IEEE Trans. Aerosp. Electron. Syst 41(1), 17–33 (2005). Publisher Full Text

BF La Scala, W Moran, RJ Evans, Optimal adaptive waveform selection for target detection. International Conference on Radar ((2003), pp), . 492–496

SD Howard, S Suvorova, W Moran, Waveform libraries for radar tracking applications. International Conference on Waveform Diversity and Design (2004)

SP Sira, A PapandreouSuppappola, D Morrell, Dynamic configuration of timevarying waveforms for agile sensing and tracking in clutter. IEEE Trans. Signal Process 7(55), 3207–3217 (2007)

SP Sira, A PapandreouSuppappola, D Morrell, Timevarying waveform selection and configuration for agile sensors in tracking applications. IEEE International Conference on Acoustics, Speech, and Signal Processing, vol. 5 ((2005), pp), . 881–884

M Orton, W Fitzgerald, A Bayesian approach to tracking multiple targets using sensor arrays and particle filters. IEEE Trans. Signal Process 50(2), 216–223 (2002). Publisher Full Text

C Kreucher, K Kastella, AO Hero III, Tracking multiple targets using a particle filter representation of the joint multitarget probability density. SPIE Int. Symp. on Optical Science and Technology, vol. 5204 ((2003), pp), . 258–259

MS Arulampalam, N Gordon, S Maskell, T Clapp, A tutorial on particle filters for online nonlinear/nonGaussian Bayesian tracking. IEEE Trans. Signal Process 2(50), 174–188 (2002)

G Björck, Functions of modulus one on whose Fourier transforms have constant modulus, and cyclic nroots. Proc. NATO Adv. Study Inst. on Recent Adv. in Fourier Analysis and its Applications, ed. by JS Byrnes, JL Byrnes ((1990), pp), . 131–140

J Benedetto, I Konstantinidis, K Okoudjou, A Bourouihiya, Concatenating codes for improved ambiguity behavior. Int. Conf. on Electromagnetics in Advanced Applications ((2007), pp), . 464–467

JJ Benedetto, I Konstantinidis, J Donatelli, C Shaw, A Doppler statistic for zero autocorrelation waveforms. Conf. on Inf. Sciences and Systems ((2006), pp), . 1403–1407

I Kyriakides, D Morrell, JJ Benedetto, I Konstantinidis, A PapandreouSuppappola, Target tracking using particle filtering and CAZAC sequences. Waveform Design and Diversity Conference ((2007), pp), . 367–371

JJ Benedetto, JJ Donatelli, Ambiguity function and frametheoretic properties of periodic zeroautocorrelation waveforms. IEEE J. Sel. Topics Signal Process 1(1), 6–20 (2007)

A Kebo, JJ Benedetto, MR Dellomo, JM Sieracki, I Konstantinidis, Ambiguity and sidelobe behavior of CAZAC coded waveforms. IEEE Radar Conference ((2007), pp), . 99–103

HLV Trees, in Detection Estimation and Modulation Theory, Part III (Wiley, New York, 1971)

Kyriakides I, On the use of Monte Carlo techniques for integrated sensing and processing (PhD thesis, Arizona State University, 2008), .

GJ Foster, T Phan, JJ Petruzzo III, Track filtering of boosting targets. Proceedings of the 35th Southeastern Symposium on System Theory, vol. 35 ((2003), pp), . 450–454

MI Skolnik, in Introduction to Radar Systems (McGrawHill, New York, 1980)

EF Knott, JF Shaeffer, MT Tuley, in Radar Cross Section, ((SciTech Publishing Inc), . , Raleigh NC, 2004)

Kyriakides I, On the validity of the measurement independence approximation when using single and MCPC waveforms based on Björck CAZAC sequences in multiple target radar tracking with a particle filter (Technical Report EEE53020081, Department of Electrical Engineering, Arizona State University, 2007, http://www), . ikyriakides.net webcite

I Kyriakides, D Morrell, A PapandreouSuppappola, Sequential Monte Carlo methods for tracking multiple targets with deterministic and stochastic constraints. IEEE Trans. Signal Process 56(3), 937–948 (2008)

R Mallik, On multivariate Rayleigh and exponential distributions. IEEE Trans. Inf. Theory 49(6), 1499–1515 (2003). Publisher Full Text