Joint Detection and DOA Tracking with a Bernoulli Filter Based on Information Theoretic Criteria.

In this paper, we study the problem of the joint detection and direction-of-arrival (DOA) tracking of a single moving source which can randomly appear or disappear from the surveillance volume. Firstly, the Bernoulli random finite set (RFS) is employed to characterize the randomness of the state process, i.e., the dynamics of the source motion and the source appearance. To increase the performance of the detection and DOA tracking in low signal-to-noise ratio (SNR) scenarios, the measurements are obtained directly from an array of sensors and allow multiple snapshots. A track-before-detect (TBD) Bernoulli filter is proposed for tracking a randomly on/off switching single dynamic system. Secondly, since the variances of the stochastic signal and measurement noise are unknown in practical applications, these nuisance parameters are marginalized by defining an uninformative prior, and the likelihood function is compensated by using the information theoretic criteria. The simulation results demonstrate the performance of the filter.

1. Introduction

Detection and direction-of-arrival (DOA) estimation using an array of sensors are important topics in signal processing and have many applications, such as in radar, sonar, communications and acoustics [1,2,3,4,5,6,7]. Considering moving sources in a noisy environment [1,2], it is desirable to use both the spatial and temporal information for better performance, since the DOAs between consecutive time steps are highly correlated. In [7,8,9,10], recursive Bayesian approaches are proposed for tracking the DOA of moving sources. However, most existing DOA tracking schemes assume that the number of sources is known and fixed. This assumption is often violated in many practical situations, since a source can appear or disappear from the surveillance volume during the observation interval. Considering an unknown and time-varying number of sources, two frameworks have emerged for tracking: track-after-detect (TAD) and track-before-detect (TBD).

In a TAD system, this task is always decomposed into separated modules. A pre-processing module is applied to obtain a set of detected points. Afterwards, a multitarget tracking (MTT) procedure is conducted on these detected points. These detected points are generated by the thresholding stage, which could result in a loss of information. The performance degrades dramatically under a low signal-to-noise ratio (SNR) or a small number of snapshots [11]. For the pre-processing module, there have been several methods developed for detecting and localizing sources using an array of sensors. The conventional beamforming method and the minimum-variance-distortionless-response (MVDR) method [12], which rely on picking peaks above a selected threshold, are the most widely used. Maximum likelihood based DOA estimation algorithms are investigated in [13] coincide with Akaike information criterion (AIC) and the minimum description length (MDL) for detection. The Bayesian predictive densities (BPD) approach was developed in [14] using a prior probability density function. Further, high-resolution subspace based methods have been studied in [15,16] with an unknown number of sources.

A TBD system tracks targets directly using the raw data, without a pre-processing module. Indeed, any detection is conducted after tracking. More information is reserved, which allows better performance in the challenging environment. A TBD particle filter is firstly proposed in [10] for DOA tracking, coupled with a reversible jump Monte Carlo Markov chain (RJMCMC) step to handle fluctuations in the source number. This method can be only carried out using a single snapshot. Hence, it is sensitive to incorrect model order initialization and degrades heavily when it comes to noise [17].

Recently, Mahler introduced the concept of random finite sets (RFS) to handle uncertainty in the number of targets for tracking [18]. Since the implementation of the RFS formulation of the optimal Bayes filter for multiple dynamic systems is computationally very demanding, some moment-based approximations have emerged recently for multisource DOA tracking: the probability hypothesis density (PHD) filter and cardinalized PHD (CPHD) filter. In [17], a PHD filter is investigated using the separable observations. A CPHD filter is investigated in [19] with known variances of signal and noise and only allows a single snapshot during recursion. In [20], a multi-Bernoulli filter for DOA tracking is proposed, by using the multiple signal classification (MUSIC) pseudo-spectrum as the likelihood function.

Usually, it is not known whether the target exists or not in a particular surveillance volume that is of interest. Our aim is to determine, from the measurement, the existence of the target and its state. The Bernoulli filter provides the optimal Bayes filter for a single dynamic system which randomly switches on and off [21]. The main feature of the Bernoulli filter is that the underlying state is treated as a set (which can be empty or singleton) instead of a vector augmented with the binary existence variable. Since it has no analytic solution, particle filter implementation provides a solution to Bernoulli recursions under non-linear/non-Gaussian cases, also known as the Bernoulli particle filter (BPF) [22]. The use of Bernoulli filter for TBD systems can be found in host of applications, such as multiple input multiple output (MIMO) radar [23], acoustic or speech sources [24,25], and similar.

In this paper, we propose a track-before-detect Bernoulli filter (TBD-Ber) for single source detecting and DOA tracking based on information theoretic criteria. The novelty of this work is twofold. Firstly, we consider an array of sensors measurement model, by using the measurements directly obtained from the array elements’ raw data. Since more information is reserved, the performance of detection and DOA tracking is improved in a noisy environment. Furthermore, it has no need of the probability of miss-detect and false alarms, which are hard to obtain in passive radar/sonar applications. Secondly, the variances of the stochastic signal and measurement noise are unknown in practical applications, known as nuisance parameters. Both the signals and noise variances can be marginalized by defining an uninformative prior. Since this is improper and unreliable for detection, we use information theoretic criteria (AIC and MDL) to compensate for the likelihood function, resulting in a penalty term. Simulation experiments are carried out and showcase the performance in challenging environments where the SNR is low and the number of snapshots is small.

The paper is organized as follows. Section 2 introduces the measurement model using an array and the formulation of Bayesian recursion based on Bernoulli RFS. Section 3 presents the derivation of the update equations and the likelihood function based on information theoretic criteria. Particle implementation is also presented. The performance metric and discussions are presented in Section 4. Simulated experiments are organized in Section 5. Finally, conclusions and future works are discussed in Section 6.

2. Background

2.1. Measurement Model Using an Array of Sensors

Consider an array composed of sensors with arbitrary locations and arbitrary directional characteristics. Assume that a narrowband signal s(k), with a known center frequency ω0, is impinging on the array from DOA , at discrete time k. For simplicity, assume that source and sensors are located in the same plane and that the source is in the far-field of the array. The received signal model for an array at time k can be written as where z(k) is a M × 1 received signal vector, n(k) is a M × 1 vector of additive noise, s(k) is a signal for a single source at the reference point, is the M × 1 manifold vector, is the propagation delay between the reference point and mth sensor to a wave front impinging from direction θ. The superscript T represents transposition.

For a source located at the far-field that moves relatively slowly, the DOA parameter is approximately stable during a small number of snapshots at each time step. At time k, N snapshots are taken into account and the corresponding DOA is θ. The received data at time k with previous N snapshots can be written by the following matrix notation:where the received data Z is a M × N matrix. Then, Equation (1) can be written as:where the source signal data is a 1 × N vector and additive noise data is M × N matrix, separately written as:

Assume both the signal S and the noise N are independent, identically distributed (i.i.d.), with a zero-mean and complex Gaussian process distributed. The source signal and the noise are independent, where I denotes a M × M identity matrix, and are the real variances of signal and noise, respectively. The distribution of the received signal also follows a zero-mean, complex Gaussian process Z~CN(0,R), where the covariance matrix is . The superscript H represents conjugate transposition. Since the true covariance matrix R is not available in practice, the measurement model is the sample-covariance matrix of received signal, given as:

2.2. Bernoulli RFS Formulation

In practice, the source of interest can enter and exit the surveillance region at random instances. Therefore, it is important to consider techniques that can jointly detect and track a source. Recently, Mahler introduces finite set statistics (FISST) [18] that provides the tool for mathematical representation and a convenient probabilistic model for the representation. The state of a time-varying number of sources can be simply represented by a finite-set valued random variable, also known as a random finite set (RFS).

The probability density function (PDF) of an RFS can be specified by a cardinality distribution and a family of symmetric joint distribution , where and denotes the distribution of its elements over the state space , conditioned on cardinality n. The PDF of a RFS X is defined as and its set integral is defined as

It is straightforward to show that f(X) integrates to one (as it should, being a PDF).

In this paper, the Bernoulli RFS X is employed to model the state of a source using a probability q to be a singleton distributed according to the ‘spatial’ PDF s(x) over the state space , given as X = {q, s(x)}. Thus, the probability of being empty is equal to 1 − q. Consider that a source with DOA θ is moving with velocity (in degree/s), the state vector is constructed by . The posterior PDF, at time step k, can be represented by a Bernoulli RFS with PDF

Considering the dynamics of source presence and absence, the dynamic model is modeled as a Bernoulli Markov process. Conditional upon , the target can re-enter the scene with probability pR, and occupy a kinematic state x with PDF , or remain absent from the scene with probability 1 − pR,. Conditional upon , the source can survive to the next time step with probability and transition to x with PDF , or disappear with probability . The dynamic model can be expressed as where the is the traditional transition density when the source survives. In this paper, the constant velocity (CV) model is employed to model the source motion when the source survives and given as where the coefficient matrix F and G are defined by where ΔT represents the time period in seconds between the previous and current time step, and is a zero-mean real Gaussian process used to model the turbulence on the source velocity. Such a constant velocity model has been widely used for DOA tracking problems [7,8,26]. For more complicated trajectories and faster moving sources, a constant acceleration model is used to model the source dynamics [27].

It is worthy of mention that the measurements of TBD filters are different from the standard multi-target tracking algorithms (i.e., TAD filters) which are the standard (points) measurements. The non-standard (intensity) measurement Z is a matrix that always exists within a fixed dimension given by (2). Hence, the likelihood function is defined as where denotes that the received signal is pure noise. Assuming the posterior PDF of the source state at the time step k − 1 is known, this is given as . The predict and update equations based on Bernoulli RFS modeling become the following

Prediction

Update

In the predict equation, is the transition density defined by (10), and is the prior distribution for the current time step. Consequently, the current distribution of the state can be obtained recursively by using this Bayesian recursion.

3. Bernoulli Filtering

3.1. Bernoulli Filter

Following from (10), the predicted PDF can also be written in the form of Bernoulli RFS, given as where the predicted probability of the existence and predicted ‘spatial’ PDF are:

Following from (15), the posterior probability of the existence and posterior ‘spatial’ PDF of the posterior PDF are given as

3.2. Likelihood Function Based on Information Theoretic Criteria

Related to the measurement model, the measurement noise process is assumed to be Gaussian. Conditional upon the likelihood function can be written as

If we have the priors of the variances of the source signal and noise, the likelihood function is straightforward. However, these parameters are always unknown in many radar/sonar applications. In order to marginalize these nuisance parameters, the measurement is split into the two complementary subspaces [13,28]. The signal subspace is spanned by the columns of matrix a(θ), and the orthogonal space is referred to as the noise subspace. According to this decomposition, the measurement vector z is the split into two subspace vectors by where z is the source signal subspace scalar and z is the (M − 1) × 1 noise subspace vector. Note that denotes an M × M unitary coordinate transformation matrix. According to these definitions we have andwhere is the M × 1 vector and is the M × (M − 1) matrix, denoting orthogonal vectors that span the source signal and noise subspaces, respectively. and denote the orthogonal projection on the source signal subspace and noise subspace, respectively. Since the transformation (22) is linear, the total likelihood function can be modified using:where and , is the Jacobian of the transformation Equation (25) and is equal to 1. Since source signal and noise are independent zero mean complex Gaussian processes, their likelihood functions are given as and where denotes a trace of a matrix, superscript represents a conjugate, and .

Intuitively, the maximum likelihood estimators of the variances of source signal and noise subspace components are and , respectively. If we do not have strong prior beliefs about nuisance parameter, noninformative prior proportional to arbitrary constants is a convenient way to reflect our ignorance for the nuisance parameter. However, such arbitrary constants are improper and unreliable for detection [29]. The information theoretic criteria have been introduced by Akaike [30], Schwartz and Rissanen [31] to compensate the arbitrary resulting penalty function of the criterion. Akaike proposed Akaike information theoretic criteria (AIC), which give the minimum AIC. Schwartz and Rissanen proposed a minimum description length (MDL), which yields the minimum code length. The idea of using the information theoretic criteria with Bernoulli Filter was first presented in [32] for sensor control.

Given a set of N observations and a family of the model defined by PDF , the generalized likelihood function can be written in the form of a penalized logarithmic, given as where is the maximum likelihood estimate of the nuisance parameter vector , and k is the number of free adjusted parameters in . Substituting the MDL criterion, the log-likelihood functions of the source signal and noise components are and

Summing up (30) and (31), the total generalized likelihood function based on MDL is given as

Following from the transformation (22), we have where is the sample-covariance matrix given by (6), and also

Taking the trace of both sides, we have

Summing up (33) and (34), taking the determinant of both sides, we have

Substituting into (32) and ignoring the constant MN, the log-likelihood function based on MDL can thus be written as where the log-likelihood function conditional upon can also be obtained using MDL criterion.

Substituting AIC criterion and separating by a factor of 2, the log-likelihood function based AIC is straightforward, given as

By using the log-likelihood function based on either criterion, the variances of source signal and noise can consequently be neglected.

3.3. Particle Implementation

The particle filter provides a general framework for the implementation for the Bernoulli filter [21,22]. At time step k − 1, the posterior PDF can be approximated by the q − 1 and L number of weighted particles , where is the state vector of particle j and is its weight. Thus, the posterior PDF is given as:where is the Dirac delta function concentrated at a. Then, the predicted probability of existence q| can be computed using Equation (17). The predicted ‘spatial’ PDF involves two terms. Assuming the probability of survival, is a constant, and the surviving source can be represented by drawing particles j = 1, …, J, given as

Assuming the probability of re-entering, is a constant, and the re-enter source can be represented by drawing particles j = 1, …, B, given as

Since there is no prior knowledge of the re-enter source, the predicted DOA θ is a uniform distribution in the state space , represents a source that can appear anywhere in the surveillance area.

After the union of two terms, the predicted ‘spatial’ PDF can be approximated by

Subsequently, the posterior probability of existence q| and ‘spatial’ PDF s|(x) can be obtained using (19) and (20), given as andwhere is the integral, given as where and are computed by the log-likelihood functions related to Equations (37) or (38) based on either MDL (MDL-TBD-Ber) or AIC (AIC-TBD-Ber) criterion.

Considering a low SNR and small number of snapshots scenarios, the mainlobe of the distribution of likelihood function is usually flat and spread. In order to enhance the high-likelihood area, we exponentially weight the penalised log-likelihood function using a constant value r to enhance performance, given as for all and . Thus, the likelihood function becomes more peaky and amenable to our problem. The choice of weighting factor can be determined based on simulations in Section 4. Finally, the particles are resampled J times and weights are equals to 1/J. If the posterior probability of existence follows q > 0.5, a source exists, and the corresponding posterior DOA is equal to

The pseudo-code of the proposed BPF is presented in Algorithm 1.

4. Discussions

4.1. Performance Metric

In order to evaluate the detection and DOA tracking performance simultaneously, a metric should be chosen that not only computes the error between the estimated and true angles, but also can quantify penalties when a source that is present is missed or a source that does not exist is falsely detected. In this section, we use the optimal sub-pattern assignment (OSPA) error metric [33] to evaluate performances. A penalty value is employed in OSPA to transfer the cardinality error into the state error and then OSPA is able to present the performance on source number estimation as well as source DOA estimation.

For joint detection and DOA tracking, the construction of the OSPA distance between two finite sets X and Y with a cardinality of at most one is as follows

The cut-off parameter determines the relative weighting of the penalties assigned to cardinality and localization errors and a moderate cut-off value will be employed in our paper.

4.2. Computational Complexity

The complexity of TBD-Ber filters are generally similar to that of the standard particle filter which can be found in [34,35]. Each particle is drawn, updated and resampled, in a similar manner. The main difference lies in the computation of the likelihood for each particle and the existence probability. From (37) and (38), the computation of the likelihood requires evaluations of the particles.

For comparison, we choose the following classical array signal processing methods. Firstly, the MVDR method is selected due to its popularity and simplicity. According to the MVDR criterion, a spatial spectrum is generated and a set of detectors is obtained by picking peaks above a selected threshold. Since it relies on estimating and inverting a covariance matrix, it requires . Intuitively, the computation load severely increases when the number of elements increases. Comparatively, grid searching is avoided for the TBD-Ber filters by drawing state samples randomly over the state-space.

Secondly, a maximum likelihood (ML) estimator coinciding with AIC or MDL criteria is selected. More specifically, it requires computing ML estimates for a series of nested models and selecting the one that best fits the underlying criteria. Usually, an alternating maximization (AM) technique is employed to transform such multi-dimensional optimization problems into a sequence of much simpler one-dimensional optimization problems that are iterated until convergence. However, it still has high computational complexity since it requires ML estimates for all candidate models.

4.3. Weighting Factor r

In our algorithm, the likelihood function is enhanced by exponentially weighting via a factor r in (26). Figure 1 shows the average OSPA error for TBD-Ber approaches (AIC-TBD-Ber and MDL-TBD-Ber) versus different weighting factors r = [1, 3, 4, 5, 7, 9, 15] over 100 Monte Carlo runs. The detailed simulation setup is provided in Section 5. As we can see, both TBD-Ber approaches have better performances when r is set around 3. It is worth mentioning that slightly larger or smaller than 3 will not lead to a significant difference.

Figure 1

Optimal sub-pattern assignment (OSPA) error averaged over time steps (100 Monte Carlo runs) versus different weighting factors r = [1, 3, 4, 5, 7, 9, 15]. Akaike information criterion track-before-detect Bernoulli filter (AIC-TBD-Ber); minimum description length track-before-detect Bernoulli filter (MDL-TBD-Ber); signal-to-noise ratio (SNR).

5. Simulations

In this section, several simulations are organized to investigate the performance of the proposed algorithm developed in this paper. Considering an unknown and time-varying number of sources, the performance is demonstrated by a single run and a Monte Carlo simulation in different experimental environments. Note that the RJMCMC method [10] can only be carried out using a single snapshot and completely fails in low SNR. In order to make a fair comparison with the proposed method which tracks the DOAs in time rather than just detect in each time step, we post-process the detections of MVDR and AIC/MDL based ML estimator using a standard Bernoulli-particle filter [21], resulting in MVDR-TAD-Ber, AIC-TAD-Ber, and MDL-TAD-Ber.

The array is a standard uniformly linear array with a spacing of using sensors. A single narrowband source signal is generated with frequency . Noise variance is assumed to be unknown and fixed. SNR is computed as . The exponentially weighting factor is which has been studied in Section 4.3. The following simulations assume there is no source present from time step 1 to time step 15, one source appears from time step 16 to time step 40 and disappears from time step 41 to time step 50. The initial DOA of the source is . The DOA velocity is around with variance .

The parameters for the TBD-Ber approaches are as follows: number of particles , and , . Since we have no prior of the source, the initial DOA is assumed to be uniformly distributed over , and the initial probability of existence is .

The parameters for the TAD-Ber approaches are as follows: number of particles , and , the probability of detection is due to low SNR simulations and the variance of the measurement noise equal to 1. The birth model is also a uniform distribution from to .

The MVDR-TAD-Ber detection threshold is a selected factor times the median MVDR power over DOA. The searching step in the pre-processing module of TAD-Ber approaches is . It is worth mentioning that it is not easy to determine the clutter rate (number of false detectors) in practical applications. We simulated several times to find a better choice of clutter rate: MVDR-TAD-Ber equals , AIC-TAD-Ber equals , and MDL-TAD-Ber equals .

5.1. Single Runs

This section gives the detection and tracking results by one trail comparing with different methods. Figure 2 presents the DOA tracking results with and snapshots . Red dots are the true source localizations and the blue circles are the tracking results. Black ‘x’s are the detections of the TAD-Ber filters which correspond to the input of a standard BPF. Intuitively, without tracking procedure, the black ‘x’s are unsatisfactory by using classical array signal processing algorithms. Although the clutter rate is carefully selected, the performance of TAD-Ber filters are still affected by the false detectors as a consequence. For the TBD-Ber filters, they not only reduce the false detectors, but also are able to detect and localize the source correctly when a source exists. It is more clear in Figure 3d,e that the probability of source existence remains low when the source does not exist. After 16th step, rapidly rises close to 1.0 when the source appears and remains close to 1.0 to time step 40.

Figure 2

Direction-of-arrival (DOA) tracking results with signal to noise ratio and snapshots . (a) minimum-variance-distortionless-response (MVDR)-TAD-Ber; (b) AIC-track-after-detect (TAD)-Ber; (c) MDL-TAD-Ber; (d) AIC-TBD-Ber; (e) MDL-TBD-Ber.

Figure 3

Detection results with and snapshots . (a) MVDR-TAD-Ber; (b) AIC-TAD-Ber; (c) MDL-TAD-Ber; (d) AIC-TBD-Ber; (e) MDL-TBD-Ber.

Figure 4 and Figure 5 consider a more severe environment where and snapshots . Accompanied with the decrease of SNR, the measurements are seriously distorted by the noise. Intuitively, the performance of TAD-Ber filters is seriously affected by the noise environment. MDL-TBD-Ber suffers from underestimation, and hence miss detects many points at the start of the trajectory. However, AIC-TBD-Ber is able to localize the DOA accurately and consistently lock on the trajectory when source exist. The better performance of can be seen in Figure 5.

Figure 4

DOA tracking results with and snapshots . (a) MVDR-TAD-Ber; (b) AIC-TAD-Ber; (c) MDL-TAD-Ber; (d) AIC-TBD-Ber; (e) MDL-TBD-Ber.

Figure 5

Detection results with and snapshots . (a) MVDR-TAD-Ber; (b) AIC-TAD-Ber; (c) MDL-TAD-Ber; (d) AIC-TBD-Ber; (e) MDL-TBD-Ber.

5.2. Monte Carlo Runs

This section presents the average performance via Monte Carlo simulations. Figure 6 and Figure 7 show the detection results for different methods under different environments. Note that the y-axis denotes the probability of detection among 100 trails. Since there is no source from time step 1 to 15 and time step 41 to 50, the y-axis denotes the probability of false alarm . Since the source exists from time step 16 to 40, the y-axis denotes the probability of correct detection . TBD-Ber filters reduce the probability of false alarm compared to TAD-Ber filters. MDL-TBD-Ber performance is affected by underestimation errors under low SNR and small number of snapshots, and thus results in an attendant loss in detection sensitivity ( is very low). AIC-TBD-Ber demonstrates better performance in all simulated scenarios.

Figure 6

Detection performance for different methods (100 Monte Carlo runs) when and number of snapshots .

Figure 7

Detection performance for different methods (100 Monte Carlo runs) when and number of snapshots .

Figure 8 shows the average OSPA error versus different SNR via 100 trails. OSPA errors are averaged over time steps and 100 trails. It presents the superior performance of TBD-Ber filters compared to TAD-Ber filters. In the high SNR region, MVDR-TAD-Ber show a comparable performance when SNR is 8 dB since the detections become so accurate that the accuracy could catch up with TBD-Ber filters. MDL-TBD-Ber gives better performance when . The performance of AIC-TBD-Ber is affected by overestimation errors and slightly worse than MDL-TBD-Ber. In the low SNR region, the performance of MDL-TBD-Ber degrades rapidly due to low SNR and small snapshots, while AIC-TBD-Ber retains better performance compared to other filters. Simulations verify the proposed track-before-detect Bernoulli filter is a more robust detection and tracking of a single source in comparison if the TAD approaches especially in noisy condition.

Figure 8

OPSA error averaged over time steps (100 Monte Carlo runs) versus SNR. Vertical panels (top to bottom) vary in number of snapshots .

The comparison of computational time for 50 time steps over 100 Monte Carlo simulations is given in Table 1. The simulations are carried out on a computer with an i5 processor and 8 GB of RAM, and all programs were coded and run in MATLAB.

Table 1

Computation time of different algorithms per iteration.

Algorithm	Mean Duration ± Standard Deviation
MVDR-TAD-Ber	0.078±0.0002
AIC-TBD-Ber	0.037±0.0010
MDL-TBD-Ber	0.042±0.0014
AIC-TAD-Ber	0.102±0.0023
MDL-TAD-Ber	0.094±0.0018

6. Conclusions

In this paper, the authors focus on solving the problem of joint detection and DOA tracking using an array of sensors. For this problem, we propose the use of measurements obtained from array elements’ raw data. We implemented the Bernoulli filter, which emerged from the RFS framework. Since more information is reserved, the performance of detection and DOA tracking is improved by using these TBD measurements. The performances of TBD-Ber and TAD-Ber approaches are evaluated based on simulations which verify the robustness of the proposed TBD-Ber filters compared to TAD-Ber filters.

Our future work will focus on more complicated situations, such as propagation path loss, or a moving source that might move to near-field from the array, etc.