1. Introduction
In many target detection and tracking applications, detection and estimation are usually carried out separately under the Bayes framework [
1]. Hypothesis tests based on the Neyman–Pearson or Bayes criterion [
2] are applied to detect the existence of the target. Sequential probability ratio test (SPRT) which makes full use of all collected data based on the continuity of the target state [
3] is also a detector. The estimators start working after the target has been declared to be detected. Examples of this type of sequential estimators include the Kalman filter (KF) [
4], the extended Kalman filter (EKF) [
5], the unscented Kalman filter (UKF) [
6], and the particle filter (PF) [
7]. Especially, the PF is a sequential Monte Carlo (SMC) method that uses a random particle system (states and weights) to approximate the relevant probability distribution. PF provides a generalized approximate solution for non-linear, non-Gaussian, and unconventional measurement problems where analytical solutions are difficult to obtain [
8]. In PF applications, importance sampling is critical. To better apply the particle filter, researchers proposed appropriate resampling methods to reduce the degeneracy of the particle system [
9,
10]. Gaussian sum particle filter is a special kind of particle filter easier to implement where the relevant probability distribution is assumed as the sum of Gaussian distributions [
11]. For low signal-to-noise (SNR) scenarios, track-before-detect (TBD) algorithms arise. In the TBD algorithm, detection is not declared for every frame. Instead, detection is declared when the estimated track is returned after processing of several data frames [
12]. In the multi-target tracking problems, multiple hypothesis on measure-to-track association can be filtered to achieve the tracking result. A comparison between multiple hypothesis Kalman filter and the multi-target particle filter can be found in Reference [
13].
In passive acoustic target surveillance applications, however, a non-cooperative target may not always radiate a signal that is observable for detection and tracking. The noisy environment increases the uncertainty of obtaining the exact target existence information from available observations. A filter that can distinguish the existence of a target by noisy observations and can be able to resume track after losing the target for a while is needed.
Recently, target tracking based on random finite set (RFS) theory provides new opportunities for state-space processing which has been widely used in many tracking examples, such as video image recognition [
14], infrared target tracking [
15], and vehicle tracking detection [
16]. Mahler and Vo systematically extended RFS theory to target tracking and made great progress in both theoretical research and applications [
17]. Bernoulli filter is one of the RFS target tracking methods where the target existence and location state is assumed to be a Bernoulli RFS. It can be used to handle a variety of problems such as multiple on/off switching systems, detection imperfections, and inaccurate estimates [
18]. Target tracking with bearings-only measurements by RFS-based filters has been studied from various aspects [
19,
20,
21,
22]. In Reference [
23], Bernoulli filter tracking along with fusion is applied in an active multi-static acoustic sensor network. Gune considered the problem of detecting and tracking a high-frequency acoustic target in the frequency-azimuth plane obtained by an acoustic vector using Bernoulli filter [
24]. However, to the best of our knowledge, modeling and filtering were conducted directly on the direction-of-arrival (DOA) estimations, and the DOA estimation process was usually neglected in the previous research.
In active radar monitoring, the measurement set of Bernoulli filter is assumed to be composed of a DOA estimation of the true target, whose error is a zero mean Gaussian variable, and several DOAs generated by clutters, which form a Poisson RFS [
25]. However, the physical basis is different in the study of passive detection and tracking of an underwater target. There are no clutters, but, in estimating target DOA by passive beamforming in a noisy environment, outputs at the sidelobes may be higher than that at the mainlobe corresponding to the target position, resulting in significant errors. To deal with this problem, we jointly use beamforming and Bernoulli filters. A more general assumption is made for the underwater acoustic signal [
26], where Gaussian noise is added directly to the received sound pressure. We conduct conventional beamforming (CBF) on the target signal band to process frequency snapshots at distributed arrays separately, and then employ the Bernoulli filter to track the target.
The main contribution of our work is that three different likelihood calculation methods based on local beamforming results are proposed for the Bernoulli filter. Firstly, we construct a DOA measurement set from the passive beamforming results (the DOA-based method). Similar to the method used in Reference [
27], to increase the probability of selecting the target DOA, we select the angles corresponding to multiple peaks of the beamforming results as the measurement set. The set contains the DOA estimation corresponding to the real target (mainlobe) and several spurious peaks caused by sidelobes. Secondly, different from the method where only DOA peaks are conserved, we propose to use the beam power as the measurement set instead, because the beam power has not only DOA peak position information but also peak intensity information. Hence, we can use beam pattern to refine the target position (or DOA, or azimuth) from the data beam power. However, the tracking results can be trapped in the end-fire direction of one array where the beam pattern is flat and ambiguous. Therefore, at last, we propose a hybrid method, aiming to combine the advantages of the DOA-based method and the beam-based method in different situations. We use the practical acoustic model to generate sound pressures in the simulations to compare the proposed methods. Filters are also used to process experimental data collected in the real world.
Throughout this paper, we use uppercase calligraphic letters to denote the RFS variable, and ∅ denotes the emptyset. For any RFS , the cardinality of is denoted by . Symbol ∖ denotes the set-difference operation. Bold lowercase letters are used to denote column vectors, and bold uppercase letters denote matrices. Spaces of variables are denoted by uppercase blackboard bold letters, where is the space of non-negative integers. We use to denote block diagonal matrix, where matrices in the brackets are diagonal blocks. is denoted as a Gaussian distribution with mean and covariance matrix , and the corresponding probability density function (PDF) with respect to is denoted as . We use f as a generic symbol for the finite set statistics (FISST) PDF of the RFS variables, while p is denoted as a generic symbol for the traditional PDF of the vector variables. For any variable or function iterated in the prediction and update process of the Bayesian filtering framework, is used to represent the prediction of from time to time k, and represents the update of with all previous states and measurements from the beginning to time k.
The rest of the paper is organized as follows. The Bernoulli filter are briefly reviewed in
Section 2, and the prediction step is provided. In
Section 3, we establish two kinds of measurement sets derived from beamforming in a multiple-array system, which leads to different update process refers to DOA-based, beam-based, and hybrid Bernoulli filter, correspondingly.
Section 4 presents the numerical implementation of the filters for practical usage. The proposed filters are implemented to process data from both simulations and real experiments. The performance results are shown in
Section 5, and the paper is finally concluded in
Section 6.
2. Bernoulli RFS Target State Model and Prediction Process
The concept of RFS was introduced by Mahler to describe the state of a randomly present target or multiple targets [
17]. RFS-based Bayes filter aims to solve problems such as detection, estimation, and data association of target tracking within a unified paradigm [
17,
28].
In the RFS-based Bayes filter, target states and sensor measurements are modeled as RFSs instead of traditional vector formulation. Suppose that targets present in the surveillance area at time k, and their state vectors are , respectively. The state of all the targets at time k is presented by an RFS on the state vector space, , where is the state vector space. All the important information about the target, including the target number () and the dynamic states of each target, has been contained. Likewise, the RFSs of sensor measurements are defined as , where is the measurement vector space.
The posterior probability density of the state variable is propagated in the predicting and updating procedure of the Bayes filter. By using the RFS state and measurement sets, we can extend the traditional Bayes filter from vector space to RFS framework. The predicting and updating process of the RFS posterior densities
can be expressed as follows [
18]:
where
is the collection of all the previous measurements and
is the likelihood function of measurement set
, given the state
. The RFS state is assumed to be Markovian with transitional density
. With the updated posterior densities
at every time step
k, RFS state
can be estimated by criteria such as maximum a posteriori probability (MAP). Although the set integrals in the RFS iteration are usually intractable, solutions can be found for joint detection and tracking of a single target.
For the case of tracking a target which may randomly be present or disappeared in a surveillance area, Bernoulli RFS can be adopted to model the target state, whose cardinality distribution is Bernoulli. Such RFS is either an empty set (with probability
) or a singleton (with probability
q), whose element
is the state vector of the target distributed following the PDF
. Thus, the FISST PDF of a Bernoulli RFS is completely determined by the target existence probability
q as
Apparently, the Bernoulli RFS can be completely determined by the corresponding existence probability q and vector PDF . Therefore, predicting and updating the RFS posterior densities are equivalent to predict and update the corresponding existence probability and vector PDF. The predicted existence probability and vector PDF are denoted by and , respectively. and are denoted as the posterior existence probability and vector PDF after update.
The kinetic characterization of the target is used in the prediction stage. This dynamic is described by the Markov transitional density (from time
to
k):
where
describes the birth probability of target with state
,
is the probability of event “target birth”, and
denotes the distribution of the state of the born target on the vector space, which is defined in advance. When prior knowledge of the target is not available, the distribution of the born target is usually assumed to be a uniform distribution over the vector space.
is the survive probability,
is the probability of “target survive”, and
is the transitional density on the vector space, which is related to the motion characteristics of the target. A typical target motion model is the constant velocity (CV) model. In a two-dimensional region, the target state vector is expressed by
, where
is the Cartesian coordinate of the target and
is the velocity component, correspondingly. The vector state transition equation of CV model is given by:
where
is the state transition matrix, and
is the zero-mean Gaussian process noise with a covariance matrix
. The diagonal blocks are given by
where
is the sampling interval from time
to time
k, and
is the potential acceleration disturbance of the moving target. Therefore, the transitional density is Gaussian,
.
With the definition of the Bernoulli RFS and its Markov transitional density, the prediction of RFS state PDF can be achieved by predicting the target existence probability and vector PDF as follows [
18]:
Equations (
7) and (
8) completely specify the prediction stage of the Bernoulli filter.
In the next section, we present the update stage for the multi-array passive beamforming system. Different kinds of measurement sets are used, which lead to different forms of the update procedure.
3. Measurement Model for Multiple Arrays and Update Process
Using a multi-array system for joint passive detection and tracking of the underwater acoustic target, we first perform beamforming on each local array to obtain the relative DOA information of the target to the local array. Subsequently, in the Bernoulli filter, with known array locations, the tracking of the target state (position and velocity in two-dimensional Cartesian coordinates) is achieved using the beamforming results of all the arrays via triangulation.
3.1. Passive Local Array Beamforming
For the convenience of practical implementation, CBF is adopted to process the collected data. We divide the signal band of interest into multiple narrow subbands and perform narrow-band CBF on every subband, and then sum up the subband results to obtain the final beam power.
By fast Fourier transform (FFT), the recorded time domain signal of a given
S-sensor array is converted to frequency-domain snapshots on
M narrow subbands. The center frequency for the
mth subband is denoted by
,
. For every subband,
L snapshots are collected to form one data frame matrix,
. The
lth column vector in matrix
is an
vector. In the plane wave assumption, for signal oriented from a target with azimuth
, elevation
, and range
r, the propagation delay difference on the array is described by the steering (or replica) vector
, where
,
, are the Cartesian positions of the array elements, and
is the wave number of the
mth frequency bin with
c being the sound speed. In the shallow water environment, depth is negligible compared to the range between target and array. That is,
, which leads to
and
. Thus, in the local array beamforming step, we only need to estimate the azimuth (DOA) of the target. Defining the data sample covariance matrix as
, where the superscript
indicates the conjugate transpose, the beam power of narrow band CBF is then
. The weight vector of CBF is exactly the driving vector,
. Summation of the beam powers of all the
M subbands leads to the basic broadband beamforming of the form
When the featured frequencies of the target are included in the summed subbands, the target DOA will be enhanced in the beam power. The estimation of DOA is then given by .
In some scenarios, for example, when SNR is low or the array sensor spacing does not meet half-wavelength requirement, strong noise at the sidelobes and grating lobes will degrade the performance of DOA estimation. The correct peak still exists in the beam power but is usually lower than the spurious peaks. In other words, the beam pattern corresponding to the correct azimuth is corrupted but remains in the presence of noises and interferences. For the purpose of making full use of the beamforming results, we propose two kinds of derivative measurements of in the following context and use them in the update procedure of Bernoulli filter.
3.2. Update Process in Multiple Array Systems
Suppose that a system equipped with
N distributed arrays is deployed in the surveillance region. Every array in the system provides a measurement set at frame
k. Denote the measurement set of the
nth array at frame
k by
,
. As the arrays are usually spatially separated, it is reasonable to make the assumption that measurements obtained at different arrays are statistically independent. Thus, the entire likelihood can be expressed as the product of local likelihood functions of individual arrays:
where
and
denote the likelihood of the
nth array.
characterizes the similarity between measurement sets to a target-oriented measurement with vector state
. On the contrary,
is the likelihood to noise. Symbol
will be replaced by specific methods in the following context. With the definition of the Bernoulli RFS and the likelihood function of the measurement RFS in Equation (
10), the update step can be expressed by updating the target existence probability
and the posterior PDF
on the vector space, respectively [
18],
where
is the measurement likelihood ratio function.
The likelihood functions and are important in the update iteration, and will be introduced in the rest part of this section.
3.3. DOA-Based Method
In this section, we discuss the construction method of a DOA-based measurement set. In the case of low SNR, the sidelobe outputs may be higher than the mainlobe in the beam power. However, our simulations showed that, if we conserve certain number of peaks, the correct DOA information of the target has a high probability to be retained, which can thus be transferred to the subsequent filter. Therefore, for each local array, we search for
largest peaks on the beam power curve
. For an individual array, the azimuths corresponding to these peaks compose the measurement set
. The measurement set can be divided into two parts
, where
is the RFS for the true target DOA, and
is the RFS for the spurious peaks. As mentioned before, true set
can be modeled as a Bernoulli RFS. When the true azimuth is not included in the selected peaks,
is an empty set. At each frame, let peak number
be randomly generated according to the Poisson distribution with the mean of
, so that the set
of spurious peaks is approximately a Poisson RFS whose PDF is defined as
where
is the density of the false DOA peak
. For the target generated true azimuth set, the conditional PDF is given by
where
is the probability of detecting the object in state
; that is, for a given state
, the probability of the corresponding azimuth peak is selected. Moreover,
is the conventional likelihood function of true azimuth measurement due to the object
. Similar to other Bernoulli filters,
is assumed to follow a Gaussian distribution
. The mean value
is the correct DOA of the target
, and the standard deviation
equals to half of the 3-dB mainlobe width of the corresponding beam pattern [
29].
In an
N-array system, let
be the azimuth peak set of the
nth array. When the target is absent,
and the measured peak set is exactly the set of spurious peaks. Therefore, the likelihood is given by
When a target is present, the likelihood is given by
The update process of the DOA-based method is completed by substituting the likelihood functions in Equations (
16) and (
17) into Equations (
11)–(
13).
3.4. Beam-Based Method
In the DOA-based method, we only retain the DOA peak positions from the beamforming outputs for filtering, but the intensity information in the beam power can also be useful. Therefore, we proposed a beam-based method in this section, where the entire beam power for of the arrays are used as the measurements, that is, , , where denotes the beam power of the signal received by the nth array at time k. According to the definition of likelihood functions, should take large value when the beam power is close to those oriented from the true direction, while gets larger when the beam power is as random as the result of the omnidirectional noise. To get a likelihood function of these features, we propose to use the concept of correlation between the data beam power and the beam pattern for some assumed signal arrival direction.
For a fixed array, the element positions in the array are known a priori. Thus, we have known beam pattern (in power) for a given DOA
of some target state
:
where
is to take the
-norm of a vector. The beam pattern of the
nth array is denoted by
. The correlation between the beam pattern and the noisy beam power can tell the likelihood of measurements to target state. In statistics, the Pearson correlation coefficient (PCC) is a measure of the linear correlation between two variables. The PCC between the measured beam power
and beam pattern
is given by
where
and
are the means of
and
with respect to
, respectively.
The PCC has a value between and , where is the perfect positive linear correlation, 0 is the no-linear correlation, and is the perfect negative linear correlation. While the likelihood functions take values in , mapping rules from PCC to likelihood functions are required. In our applications, when the local beamforming result does not match the beam pattern from the assumed direction, the PCC is negative or close to 0, indicating that there is only noise in this direction, so that the likelihood function to noise should be larger. On contrast, when the local beamforming result completely matches the beam pattern, PCC is equal to 1, indicating that there is a target in the assumed direction. The closer the PCC is to 1, the larger the target likelihood function should be set.
Note that negative PCCs are considered as mismatch here. Readers may be confused because negative correlations are usually considered as strong correlations in applications such as image recognition. That is, if the direction of a measured image are opposite to the original one, the PCC between them equals to −1, but they are still the same object (highly correlated). However, in our application, the PCC between the local beamforming result and the beam pattern being negative means that peaks appear at the DOA positions that should have been valleys, while false peaks appear at where they should not have been existed. This is not a match in our expectation. If the negative PCC was still considered as strong correlation, we would get wrong DOA estimation.
The desired mapping rules should satisfy the requirements listed above. For this purpose, we first define a non-negative PCC
, where
if
and
if
. Then, we generate likelihood functions from
as follows:
where
is an adjustable probability mapping scale parameter. The update process of the beam-based method is completed by substituting the likelihood functions in Equation (
20) into Equations (
11)–(
13).
3.5. Hybrid Method
The beam-based method is beneficial to improve the accuracy when the beamforming result highly matches the beam pattern of the target DOA. However, since there are wider peaks near the end-fire direction in the beam pattern of an individual local array, the tracking results can sometimes be trapped in the end-fire regions. If there are more than three arrays, the effect of such error can be eliminated thanks to the geometric relations between the arrays. In the case of fewer arrays, this problem leads to localization error. Although the DOA-based method may also have ambiguity problems in the end-fire direction, since only one or two DOA peaks are retained for filtering, it is less affected than the beam-based method where the operation of calculating the PCC will aggravate the performance degradation caused by the ambiguity in the end-fire direction. Therefore, the method using DOA measurements only is advantageous for protecting the true value from the effects of the end-fire direction in the filtering procedure, especially when the number of arrays is not enough. Both methods have their own advantages and disadvantages. To complement each other, we propose a hybrid approach combining these two choices. We define the new likelihood functions as the geometric mean of likelihood functions of the DOA-based and beam-based methods as follows:
where
is the portion coefficient to adjust the weights of beam-based and DOA-based likelihood in the hybrid method. When
, the function is the same as the beam-based likelihood, while
, it is exactly the same as that of the DOA-based method.
4. SMC Implementation of the Bernoulli Filter
The PF/SMC method provides a general framework for the implementation of Bernoulli filter [
18]. The PDF
is approximated by a particle system
, where
is the state of particle
i and
is the normalized particle weight, that is,
. The PDF is approximated by
where
is the Dirac function at point
. In the prediction stage, the prediction of PDF in Equation (
8) is completed by predicting the particle system
The particles are divided into two parts: newly born particles and persistent particles. The states of the birth particles are drawn from a pre-defined birth distribution, for example, uniform distribution. The birth particles are equally weighted by . The persistent particles are the extension of those particles at the previous time. Their states are drawn from the state transitional density for . The relationship between the weights of the persistent particles and the previous particle weights is given by .
Note that we use proportion rather than equality, because the weights will also be normalized in the prediction stage.
In the update equations of target existence probability, the integrals can be approximated using predicted particles as
To update PDF
, the weights of particles are updated by
At the end of the update procedure, sampling importance resampling (
particles persist) is adopted to regularize the particle system [
9]. With all these particles, the MAP estimation of the target state is approximately the weighted summation of the updated particles at time
k, that is,
.