Adaptive highly localized waveform design for multiple target tracking

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 phase-coded (MCPC) waveforms based on Björck constant amplitude zero-autocorrelation (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 high-resolution 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 non-adaptive waveform configurations or single-carrier phase-coded CAZAC waveforms.

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 mean-squared error (MSE) or percentage of lost tracks, is poor. The location and magnitude of the AF measurement sidelobes in the delay-Doppler 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 [1-3]. 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 delay-Doppler 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 frequency-modulated 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 (IP-LPF) tracker adaptively configures multicarrier phase-coded (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 delay-Doppler 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 IP-LPF 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 delay-Doppler 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 non-zero 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 equal-length Björck constant amplitude zero-autocorrelation code (CAZAC) sequences [16-18] that are cyclic-shifted versions of one another. A Björck CAZAC sequence provides higher-resolution measurements than a linearly frequency-modulated 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 single-carrier phase-coded (SCPC) CAZAC sequences exhibit sidelobes that are spread in the delay-Doppler 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 “IP-LPF algorithm”, we provide a detailed description of the IP-LPF algorithm and its application in minimizing the predicted tracking error. In Section “Adaptive waveform selection”, we integrate the IP-LPF with a waveform configuration algorithm and demonstrate its performance in tracking multiple targets in Section “Simulation results”.

Björck CAZAC sequences

A CAZAC sequence ξ(m) with finite length M, has constant magnitude, |ξ(m)| = 1, m = 0,…,M − 1, and zero-autocorrelation, , 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 signal-to-noise ratio (SNR). They also exhibit very tight localization in the delay-Doppler plane that can enhance the range resolution and range-rate 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 delay-Doppler 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.

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 frequency-shifted CAZAC is also a CAZAC [20] and also that a sum of cyclic frequency-shifted 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 frequency-shifted 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 time-bandwidth 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 time-bandwidth 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 frequency-shift

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 non-zero 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 delay-Doppler measurement plane near the AF origin is the area that is most commonly interrogated by the IP-LPF tracker when accurately tracking a target, as we will see in Section “IP-LPF 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 non-zero 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 non-zero sidelobes exist when n is an integer multiple of Q = 130.

Positive cyclic frequency-shift

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 numerical-based 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 non-zero 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 in-between 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 non-zero 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 non-zero 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”.

Tracking model

We consider L targets moving in a two-dimensional (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 zero-mean, 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_cT_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 zero-mean, 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 zero-mean complex Gaussian with variance 2 N₀ 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 time-frequency 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, zero-mean, 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₀ (not target is assumed present) and under hypothesis H₁(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 non-zero, 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 IP-LPF, 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 delay-Doppler 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 delay-Doppler 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

where and .

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 delay-Doppler pair in the template signal representing a single partition λ. The single partition likelihood for this delay-Doppler bin is given by:

We evaluate the likelihood ratio for each delay-Doppler 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_cT_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.

Stage 2: partitions sampling

During Stage 1, we propose partitions , i = 1,…,N, from delay-Doppler 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 delay-Doppler plane derived from the proposed partitions, the measurements from sensors u = 1,…,U are given by

where is the delay-Doppler 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 delay-Doppler 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 right-hand 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 delay-Doppler 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.

IP-LPF algorithm

■ For each partition λ = 1,…,Λ and for each particle i = 1,…,N

Stage 1: Likelihood Partition Sampling

✶ Let

✶ For each sensor u = 1,…,U

✻ For and for

♦ Form using (13) and using (14)

♦ Evaluate using (15) and using (19)

♦ Sample

♦ Let

♦ Let and

✶ Calculate from and

Stage 2: Likelihood Partition Sampling

✶ For each particle ɨ =1,…,N

✻ Evaluate using (20) and using (22)

✻ Sample

✶ Let and

Particle Weighting

■ For each particle i = 1,…,N

✶ Assemble particles

✶ Evaluate using (24)

✶ For each particle ɨ =1,…,N

✻ Evaluate particle weights:

■ Estimate

■ Increment k by 1

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 mean-squared 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 delay-Doppler locations associated with it. As described in Section “IP-LPF algorithm”, in order to propose particles, we have considered a discrete finite set of delay-Doppler 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 IP-LPF proposal and subsequently weighted. Therefore, is a large set including all combinations of possible delay-Doppler locations from two of the sensors for each target and particle. The process of forming partitions from sampled delay-Doppler 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 weight-particle pairs corresponding to delay-Doppler 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 delay-Doppler 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 IP-LPF, we will follow the sampling process of the IP-LPF 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 delay-Doppler 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 delay-Doppler values in (13) and (14) with the additional constraint that for all partitions λ, particles i and sensors u. These sets of delay-Doppler locations correspond to the belief on target state as explained previously and only include delay-Doppler 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 delay-Doppler values in (13) and (14) and such that . Since the matched filter statistic variance depends on the AF of the waveform, and since the above-mentioned set of delay-Doppler 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 delay-Doppler plane that is given by the sequential tracking process of the particle filter as explained in Section “IP-LPF 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 state-of-the 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 CAZAC-based 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.

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 χ₁ = −1000 m, ψ₁ = 500 m, χ₂ = 2500 m, ψ₂ = 500 m, and χ₃ = 500 m, ψ₃ = 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).

In the first scenario, the weak target l = 2 has a cross-sectional 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₀ = 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⁸ 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₀ = 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 non-adaptive 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.

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 non-adaptive 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.

We developed the IP-LPF 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 high-resolution 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.

The authors declare that they have no competing interests.

This study was supported under MURI Grant No. AFOSR FA9550-05-1-0443.

1. N Levanon, E Mozeson, Radar Signals (Wiley, New York, 2004)

2. C Rago, P Willett, Y Bar-Shalom, Detection-tracking performance with combined waveforms. IEEE Trans. Aerosp. Electron. Syst 34(2), 612–624 (1998). Publisher Full Text

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

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

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

6. S-M 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

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

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

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

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

11. SP Sira, A Papandreou-Suppappola, D Morrell, Dynamic configuration of time-varying waveforms for agile sensing and tracking in clutter. IEEE Trans. Signal Process 7(55), 3207–3217 (2007)

12. SP Sira, A Papandreou-Suppappola, D Morrell, Time-varying 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

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

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

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

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

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

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

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

20. JJ Benedetto, JJ Donatelli, Ambiguity function and frame-theoretic properties of periodic zero-autocorrelation waveforms. IEEE J. Sel. Topics Signal Process 1(1), 6–20 (2007)

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

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

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

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

25. MI Skolnik, in Introduction to Radar Systems (McGraw-Hill, New York, 1980)

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

27. 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 EEE-5-30-2008-1, Department of Electrical Engineering, Arizona State University, 2007, http://www), . ikyriakides.net webcite

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

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