Riemann-Langevin Particle Filtering in Track-Before-Detect

TL;DR

This paper introduces the Riemann-Langevin particle filter, a novel approach that leverages Langevin Monte Carlo sampling to improve online state estimation in track-before-detect scenarios, especially under non-linear measurement models.

Contribution

It presents the first application of Riemann-Langevin particle filtering to track-before-detect, enhancing proposal distributions by incorporating measurement information.

Findings

Improved state estimation accuracy in low-noise scenarios.

Enhanced particle filter performance by using measurement-informed proposals.

Demonstrated effectiveness over traditional methods in challenging conditions.

Abstract

Track-before-detect (TBD) is a powerful approach that consists in providing the tracker with sensor measurements directly without pre-detection. Due to the measurement model non-linearities, online state estimation in TBD is most commonly solved via particle filtering. Existing particle filters for TBD do not incorporate measurement information in their proposal distribution. The Langevin Monte Carlo (LMC) is a sampling method whose proposal is able to exploit all available knowledge of the posterior (that is, both prior and measurement information). This letter synthesizes recent advances in LMC-based filtering to describe the Riemann-Langevin particle filter and introduces its novel application to TBD. The benefits of our approach are illustrated in a challenging low-noise scenario.

Tables

Table 1. TABLE I: Sensor parameters

Range PSF constant (R)	$1.56\cdot {10}^6$ $\left[{\text{m}}^2\right]$
Bearing PSF constant (B)	$1.88\cdot {10}^{-4}$ $\left[{\text{rad}}^2\right]$
Range resolution	$500$ $\left[\text{m}\right]$
Bearing resolution	$5\cdot {10}^{-3}$ $\left[\text{rad}\right]$
Range lower bound	$22\cdot {10}^3$ $\left[\text{m}\right]$
Range upper bound	$26\cdot {10}^3$ $\left[\text{m}\right]$
Bearing lower bound	$-\pi /6$ $\left[\text{rad}\right]$
Bearing upper bound	$\pi /6$ $\left[\text{rad}\right]$

Equations

$$\min \left\{1,\frac{\text{p}({\mathbf{z}}_{\text{k}}|{\mathbf{s}}_{\text{k}}^{\ast })\text{p}({\mathbf{s}}_{\text{k}}^{\ast }|{\mathbf{z}}_{1:\text{k}-1})}{\text{p}({\mathbf{z}}_{\text{k}}|{\mathbf{s}}_{\text{k}}^i)\text{p}({\mathbf{s}}_{\text{k}}^i|{\mathbf{z}}_{1:\text{k}-1})}\frac{\text{q}\left({\mathbf{s}}_{\text{k}}^i|{\mathbf{s}}_{\text{k}}^{\ast }\right)}{\text{q}\left({\mathbf{s}}_{\text{k}}^{\ast }|{\mathbf{s}}_{\text{k}}^i\right)}\right\}.$$

$$\text{p}({\mathbf{s}}_{\text{k}},{\mathbf{s}}_{\text{k}-1}|{\mathbf{z}}_{1:\text{k}})\propto \text{p}({\mathbf{z}}_{\text{k}}|{\mathbf{s}}_{\text{k}})\text{p}({\mathbf{s}}_{\text{k}}|{\mathbf{s}}_{\text{k}-1})\text{p}({\mathbf{s}}_{\text{k}-1}|{\mathbf{z}}_{1:\text{k}-1}),$$

$$\min \left\{1,\frac{\text{p}({\mathbf{z}}_{\text{k}}|{\mathbf{s}}_{\text{k}}^{\ast })\text{p}({\mathbf{s}}_{\text{k}}^{\ast }|{\mathbf{s}}_{\text{k}-1}^i)}{\text{p}({\mathbf{z}}_{\text{k}}|{\mathbf{s}}_{\text{k}}^i)\text{p}({\mathbf{s}}_{\text{k}}^i|{\mathbf{s}}_{\text{k}-1}^i)}\frac{\text{q}\left({\mathbf{s}}_{\text{k}}^i|{\mathbf{s}}_{\text{k}}^{\ast },{\mathbf{s}}_{\text{k}-1}^i\right)}{\text{q}\left({\mathbf{s}}_{\text{k}}^{\ast }|{\mathbf{s}}_{\text{k}}^i,{\mathbf{s}}_{\text{k}-1}^i\right)}\right\}.$$

$${\text{q}}_{\text{RL}}\left({\mathbf{s}}_{\text{k}}^{\prime }|{\mathbf{s}}_{\text{k}},{\mathbf{s}}_{\text{k}-1},{\mathbf{z}}_{\text{k}}\right)=\mathcal{N}\left({\mathbf{s}}_{\text{k}}^{\prime };{\mathbf{s}}_{\text{k}}+\frac{{ϵ}^2}{2}{\text{F}}^{-1}\left({\mathbf{s}}_{\text{k}},{\mathbf{s}}_{\text{k}-1},{\mathbf{z}}_{\text{k}}\right){\nabla }_{{\mathbf{s}}_{\text{k}}}\log \left(\text{p}({\mathbf{z}}_{\text{k}}|{\mathbf{s}}_{\text{k}})\text{p}({\mathbf{s}}_{\text{k}}|{\mathbf{s}}_{\text{k}-1})\right),{ϵ}^2{\text{F}}^{-1}\left({\mathbf{s}}_{\text{k}},{\mathbf{s}}_{\text{k}-1},{\mathbf{z}}_{\text{k}}\right)\right)$$

$$\text{F}\left({\mathbf{s}}_{\text{k}},{\mathbf{s}}_{\text{k}-1},{\mathbf{z}}_{\text{k}}\right)=-{\mathbb{E}}_{{\mathbf{z}}_{\text{k}}|{\mathbf{s}}_{\text{k}}}\left[{\Delta }_{{\mathbf{s}}_{\text{k}}}^{{\mathbf{s}}_{\text{k}}}\log \left(\text{p}({\mathbf{z}}_{\text{k}}|{\mathbf{s}}_{\text{k}})\text{p}({\mathbf{s}}_{\text{k}}|{\mathbf{s}}_{\text{k}-1})\right)\right],$$

$${\nabla }_{\mathbf{\alpha }}\log \varphi \left(\mathbf{x}\right)$$

$$\text{F}\left(\mathbf{\alpha }\right)$$

$${\mathbf{z}}_{\text{k}}=[{\text{z}}_{\text{k}}^{(1)}\cdots {\text{z}}_{\text{k}}^{(\text{J})}{]}^{\text{T}}.$$

$${\text{z}}_{\text{k}}^{(j)}={\hat{\text{z}}}_{\text{k}}^{(j)}+w_{\text{k}}^{(j)}=\text{A}{\text{h}}^{(j)}({\mathbf{s}}_{\text{k}})+w_{\text{k}}^{(j)},$$

$${\text{h}}^{(j)}({\mathbf{s}}_{\text{k}})=\exp \left(-\frac{(r_j-r({\mathbf{s}}_{\text{k}}){)}^2}{2\text{R}}-\frac{(b_j-b({\mathbf{s}}_{\text{k}}){)}^2}{2\text{B}}\right),$$

$$\text{p}({\mathbf{z}}_{\text{k}}|{\mathbf{s}}_{\text{k}})=\prod _{j=1}^{\text{J}}\text{p}({\text{z}}_{\text{k}}^{(j)}|{\mathbf{s}}_{\text{k}}).$$

$${\nabla }_{{\mathbf{s}}_{\text{k}}}\log \text{p}({\mathbf{z}}_{\text{k}}|{\mathbf{s}}_{\text{k}})=\sum _{j=1}^{\text{J}}\frac{1}{{\sigma }_w^2}{\left(\frac{\partial {\hat{\text{z}}}_{\text{k}}^{(j)}}{\partial {\mathbf{s}}_{\text{k}}}\right)}^{\text{T}}({\text{z}}_{\text{k}}^{(j)}-{\hat{\text{z}}}_{\text{k}}^{(j)}),$$

$$\text{F}\left({\mathbf{s}}_{\text{k}}\right)=\sum _{j=1}^{\text{J}}\frac{1}{{\sigma }_w^2}{\left(\frac{\partial {\hat{\text{z}}}_{\text{k}}^{(j)}}{\partial {\mathbf{s}}_{\text{k}}}\right)}^{\text{T}}\left(\frac{\partial {\hat{\text{z}}}_{\text{k}}^{(j)}}{\partial {\mathbf{s}}_{\text{k}}}\right).$$

$${\mathbf{s}}_{\text{k}}={\left[\begin{array}{cccc}{x}_{\text{k}}& {\dot{x}}_{\text{k}}& y_{\text{k}}& {\dot{y}}_{\text{k}}\end{array}\right]}^{\text{T}}.$$

$$\text{A}={\text{I}}_2\otimes \left[\begin{array}{cc}1& {\Delta }_t\\ 0& 1\end{array}\right]$$

$$\text{Q}=\left[\begin{array}{cc}{\sigma }_{a_x}^2& 0\\ 0& {\sigma }_{a_y}^2\end{array}\right]\otimes \left[\begin{array}{cc}{\Delta }_t^3/3& {\Delta }_t^2/2\\ {\Delta }_t^2/2& {\Delta }_t\end{array}\right],$$

Full text

I Introduction

Sequential state estimation in nonlinear dynamical systems is a challenging problem. Following a Bayesian approach, a closed-form expression of the conditional probability of the state (posterior) is only attainable for a restricted class of models. Therefore, methods based on numerical approximations are oftentimes employed. Among these approximations, Monte Carlo (MC) methods [1, 2, 3] are popular due to their flexibility and provable convergence guarantees. Particle filters (PFs) based on importance sampling (IS) are a straightforward implementation of MC methods in dynamical systems. However, their practical application is quickly challenged as the dimension of the state space increases. Moreover, IS-based PF usually require resampling in order to avoid sample degeneracy, restraining their parallel implementation. Due to these limitations, the integration of Markov chain Monte Carlo (MCMC) in PF has emerged as competing alternative [4, 5, 6, 7].

This paper investigates the application of a certain type of MCMC-based PF to the problem of track-before-detect (TBD) [8, 9, 10]. TBD consists in removing the detection stage typically found before the tracking module in surveillance and tracking systems. Consequently, the “unthresholded” image measurements are fed directly into the tracker. Due to the non-linearities in TBD measurement models, PFs are commonly employed. Nevertheless, the proposal distributions in PFs (based on either MCMC or IS) developed for TBD have solely exploited the prior. This is in general inefficient, but especially when most of the posterior information is contained in the measurement. In this context, the Langevin Monte Carlo (LMC), which is a MCMC whose proposal exploits measurement knowledge through the Langevin equation [11], emerges as an efficient alternative that has been recently applied in the context of classical tracking [12].

The main goal of this letter is the application of LMC-based PFs to TBD. Together with the particularization to the TBD setup, we also present several modifications to enhance the original LMC filtering in [12]. One of them is related to the fact of LMC being parameterized by a step size that must be carefully adjusted for every time step. Ideally, an LMC-based filter should systematically adapt its step size. Here, this problem is avoided employing the Riemannian MCMC [13]. Furthermore, the filter in [12] approximates the gradient of the posterior and uses an MCMC-based PF methodology that introduces a bias [14]. In this letter these limitations are addressed using the sequential MCMC approach [6, 15].

The letter is organized as follows. Section II discusses fundamentals of MCMC-based PFs and sequential MCMC, with a special focus on the acceptance probability. The Langevin proposal and its extension based on differential geometry are described in Section III. The Fisher information matrix [16], plays a key role on this extension. Hence, Section III-A derives the Fisher information matrix for non-linear Gaussian models and Section III-B particularizes it for TBD. Section IV details the Riemann-Langevin MC filter. The performance of the filter is analyzed in a original low-noise TBD application in Section V. The conclusions in Section VI close the letter.

II Preliminaries: Sequential MCMC

This section introduces notations to be used in the letter and illustrates how sequential MCMC renders an efficient MCMC-based PF. We first introduce the general expression for the probability of accepting a sample, and then show how this expression simplifies for sequential MCMC. These methods were introduced in a tracking application [6], with the goal of enabling efficient computation of the acceptance probability when performing block sampling (e.g. sampling one object at-a-time). The gains in sequential MCMC come from sampling the joint posterior at two consecutive time steps.

To be specific, let the state and measurement vectors at time step k be denoted by $\mathbf{s}_{\text{k}}$ and $\mathbf{z}_{\text{k}}$, respectively, and let $\mathbf{z}_{1:\text{k}}$ denote the sequence of measurements $\mathbf{z}_{1},\dots,\mathbf{z}_{\text{k}}$. Furthermore, the Bayesian posterior is $\text{p}({\mathbf{s}_{\text{k}}|\mathbf{z}_{1:\text{k}}})$, the prediction $\text{p}({\mathbf{s}_{\text{k}}|\mathbf{z}_{1:\text{k}-1}})$, and the measurement likelihood $\text{p}({\mathbf{z}_{\text{k}}|\mathbf{s}_{\text{k}}})$. Suppose that an MCMC algorithm is used to sample from $\text{p}({\mathbf{s}_{\text{k}}|\mathbf{z}_{1:\text{k}}})$. Let the current sample of the Markov chain be denoted by $\mathbf{s}_{\text{k}}^{i}$ and the proposed sample by $\mathbf{s}_{\text{k}}^{*}$, drawn from the proposal distribution q. A key component of a MCMC is to decide if $\mathbf{s}_{\text{k}}^{*}$ is accepted or rejected. The general expression for the acceptance probability in the Metropolis-Hastings algorithm involves the ratios of the posteriors and the proposals and is given by

[TABLE]

Clearly, unless the proposal is simply equal to the prediction, computing (1) requires evaluating the density $\text{p}({\mathbf{s}_{\text{k}}|\mathbf{z}_{1:\text{k}-1}})$.

The sequential MCMC brings flexibility in choosing the proposal while maintaining efficient computation of the acceptance probability. The key is to consider jointly the state at k and $\text{k}-1$. Within each time step of the PF, an MCMC samples from the joint density $\text{p}({\mathbf{s}_{\text{k}},\mathbf{s}_{\text{k}-1}|\mathbf{z}_{1:\text{k}}})$. Due to the factorization of the joint density,

[TABLE]

a sample from a proposal q on $\mathbf{s}_{\text{k}}$ is accepted with probability

[TABLE]

Note that the factors corresponding to $\text{p}({\mathbf{s}_{\text{k}-1}|\mathbf{z}_{1:\text{k}-1}})$ in (2) simplify. Compared to (1), the sequential MCMC’s acceptance probability (3) depends on the evaluation of the transition model $\text{p}({\mathbf{s}_{\text{k}}|\mathbf{s}_{\text{k}-1}})$ in lieu of the prediction $\text{p}({\mathbf{s}_{\text{k}}|\mathbf{z}_{1:\text{k}-1}})$.

III Riemann-Langevin proposal in Track-Before-Detect

This section introduces the application of the LMC filter to tracking applications and presents the enhancements described in Section I. But first, let us review the generalization of the Langevin proposal inspired by differential geometry (namely, the Riemann-Langevin proposal) [13].

We start with some notational conventions: $(\mathbf{s}_{\text{k}},\mathbf{s}_{\text{k}-1})$ is the current sample or state of the Markov chain, $\mathcal{N}\left(\mathbf{x};\bm{\mu},\Sigma\right)$ is a Gaussian density with mean $\bm{\mu}$ and covariance $\Sigma$, $\nabla_{\mathbf{x}}\ \text{f}(\mathbf{x})$ is the gradient of f w.r.t. $\mathbf{x}$, and F is a metric tensor representing the Fisher information matrix. With these conventions, a proposed sample $\mathbf{s}_{\text{k}}^{\prime}$ is obtained from the Riemann-Langevin proposal as given in (4) at the bottom of the page.

This proposal leverages the property that probability distributions lie in Riemann manifolds and stems from extending the Langevin equation with the metric tensor [13]. A sensible choice of the metric tensor containing curvature information about the posterior is the Fisher information matrix [16]. In sequential MCMC it is defined as [15]

[TABLE]

where $\mathbb{E}_{\mathbf{z}_{\text{k}}|\mathbf{s}_{\text{k}}}$ denotes the expected value w.r.t. $\text{p}({\mathbf{z}_{\text{k}}|\mathbf{s}_{\text{k}}})$ and $\Delta_{\mathbf{x}}^{\mathbf{x}}\ \text{f}(\mathbf{x})$ the Hessian of f w.r.t. $\mathbf{x}$.

III-A Information matrices in tracking applications

Sampling from (4) requires the information matrix of the product $\text{p}({\mathbf{z}_{\text{k}}|\mathbf{s}_{\text{k}}})\text{p}({\mathbf{s}_{\text{k}}|\mathbf{s}_{\text{k}-1}})$. If the models are Gaussian and linear, the information matrix is simply the covariance matrix. However, important tracking applications such as TBD have non-linear $\text{p}({\mathbf{z}_{\text{k}}|\mathbf{s}_{\text{k}}})$. Therefore, this section reviews the information matrix of a non-linear multi-variate Gaussian. The next section will build upon this result and provide the information matrix for TBD.

Consider a Gaussian density function $\phi$ with covariance $\Sigma$ and mean $\bm{\mu}$. Suppose that $\bm{\mu}$ is a function of a vector $\bm{\alpha}$ and that $\Sigma$ is independent of $\bm{\alpha}$. Then, the gradient and information matrix w.r.t $\bm{\alpha}$ are [17]:

[TABLE]

III-B Track-before-detect

In this section we provide the gradient and information matrix of the TBD measurement model in [9].

III-B1 Measurement model

Consider an imaging sensor (e.g. camera, radar) that measures the vector $\mathbf{z}_{\text{k}}$ composed of the scalars $\text{z}_{\text{k}}^{(j)}$ in J pixels or cells:

[TABLE]

Suppose that the strength of the measured signal under the influence of an object is constant and denoted by A. In addition, let the measurement noise $w_{\text{k}}^{(j)}$ be zero-mean Gaussian with variance $\sigma_{w}^{2}$, so that the signal-to-noise ratio (SNR) is $20\log_{10}\left(\text{A}/\sigma_{w}\right)$. Then, the measurement in a single cell is

[TABLE]

where $\text{h}^{(j)}$ is the point spread function characteristic of the sensor. While our results hold for any point spread function, a popular choice (also used for the experiments in Section V) is [9]

[TABLE]

where $r_{j}$ and $b_{j}$ denote, respectively, the range and bearing cell centroids; $r(\mathbf{s}_{\text{k}})$ and $b(\mathbf{s}_{\text{k}})$ the position in polar coordinates; and R and B are sensing parameters that depend on the resolution and the boundaries of the surveillance area [9]. Other measurement dimensions (e.g. Doppler velocity, elevation) can be included in the point spread function depending on the type of imaging sensor. For illustration, Figure 1 shows a sensor grid and measurement following this model.

According to (9), $\text{p}({\text{z}_{\text{k}}^{(j)}|\mathbf{s}_{\text{k}}})$ is Gaussian with mean $\widehat{\text{z}}^{(j)}$ and variance $\sigma_{w}^{2}$. Assuming that measurements among cells are conditionally independent, the complete likelihood is given by

[TABLE]

III-B2 Riemann-Langevin proposal

The TBD likelihood in (11) consists of a product of J Gaussians. The log-gradient and information matrix of each of these Gaussians obey (6) and (7), respectively. Since the logarithm of the product becomes the sum of logarithms, we have that the TBD gradient and information matrix are

[TABLE]

[TABLE]

IV Riemann-Langevin Monte Carlo filter

Leveraging the results in the previous sections, including the adoption of a sequential MCMC approach and the Riemann-Langevin proposal, here we provide the full description of our Riemann-Langevin MC filter. The details of all the steps required to implement our filter are given in Algorithm 1. The algorithm has an outer-most loop in k corresponding to the filtering time. The loop in $i$ represents the evolution of the Markov chain in the MCMC. Being a sequential MCMC, the MCMC contains two parts corresponding to the joint and refinement phases [15]. Regarding the rest of the symbols introduced in Algorithm 1: $\widehat{\text{p}}$ denotes the sample-based representation of a density; $\text{N}_{\text{bi}}$ the burn-in length; and $\text{N}_{\text{p}}$ the number of particles.

V Experiments

The performance of the Riemann-Langevin MC filter in Section IV leveraging second-order model knowledge is assessed in this section. In particular, we compare this filter to the traditional bootstrap PF and the standard sequential MCMC in a TBD application. We start by describing the tracking scenario (including the trajectory generation, motion and measurement models) and then the results are presented and analyzed.

The scenario is composed of a single object moving along a straight-line at a constant speed equal to 180 kilometers per hour. The object trajectory lasts 30 seconds and is shown in Figure 1. Range and bearing measurements are reported by an imaging sensor every second (i.e., sampling time $\Delta_{t}=1[s]$).

Unlike common high-dimensional experiments demonstrating the power of Hamiltonian/Langevin MC [18, 13, 15], our state space is only composed of four variables. Our motivation is to demonstrate that the Riemann-Langevin MC filter is not only useful in high-dimensions, but also in other circumstances challenging for PFs such as low-noise.

V-A Motion model

The state vector comprises two-dimensional (denoted by $x$ and $y$) position and velocity vectors:

[TABLE]

The time evolution is given by a nearly constant velocity model [19] $\mathbf{s}_{\text{k}}=\text{A}\ \mathbf{s}_{\text{k}-1}+\mathbf{v}_{\text{k}-1}$, with transition matrix

[TABLE]

where $\otimes$ denotes the Kronecker product. The noise $\mathbf{v}_{\text{k}}$ is Gaussian with zero mean and covariance

[TABLE]

where $\sigma_{a_{x}}=\sigma_{a_{y}}=0.1\ [\text{m}\ \text{s}^{-2}]$ denote scalar standard deviations of the acceleration along the $x$ and $y$ axes, respectively.

V-B Measurement model and sensor parameters

The TBD measurement model is described in Section III-B. The specification of the parameters follows. The measurement noise is $\sigma_{w}=10^{-4}$ and the $\text{SNR}=80\ [\text{dB}]$, simulating a very low-noise scenario especially challenging for PF. The sensor parameters are listed in Table I. Figure 1 depicts the sensor grid in a region of the field-of-view as well as a sample measurement.

V-C Tested methods and initialization

We compare the following methods:

The Riemann-Langevin MC filter (Section IV). 2. 2.

The bootstrap PF based on IS and resampling [1]. 3. 3.

Sequential MCMC with prior proposal [15].

The Riemann-Langevin MC filter uses 400 particles, the bootstrap PF 5000, and the sequential MCMC with prior proposal 3000. For both MCMC methods the burn-in length is 100. Initial particles are drawn from two uniform distributions: for the position the area of the distribution is $1\ [\text{km}^{2}]$; and for the velocity it is $100\ [\text{m}^{2}\ \text{s}^{-2}]$. Both distributions are centered around the ground truth.

V-D Numerical results

The root mean squared error (RMSE) of the position estimate ($x$ and $y$ axes) is shown in Figure 2. Results are obtained averaging $\text{N}_{\text{MC}}=50$ MC simulations. In all methods the estimated positions are given by the sample average over the population of particles (burn-in is discarded in the MCMC methods). The results show that the Riemann-Langevin MC filter outperforms the other methods despite its lower number of samples. Note that during the first few time steps the low performance of the Riemann-Langevin MC filter is due to its fewer number of particles. The superior overall performance stems from the Riemann-Langevin proposal, which leverages both prior and measurement information, whereas the proposals in the other methods only contain prior information.

In particular, the particle clouds in the Riemann-Langevin MC filter are diverse, whereas the particle clouds in the other methods are degenerated due to the low noise. The Riemann-Langevin proposal achieves diverse clouds via adaptation, focusing on the regions of the state space where most mass of the posterior density lies. A measure of the particles’ degeneration (dispersion) with straightforward interpretation is the number of distinct particles. Across all the MC runs, the minimum and maximum number of distinct particles at the last time step are: 363 and 384 (of 400 total) in the Riemann-Langevin MC filter; 4 and 8 (of 5000) in the bootstrap PF; 4 and 11 (of 3000) in the sequential MCMC with prior proposal.

VI Conclusion

This letter presented the application of a differential-geometric MCMC-based PF to the problem of TBD. The proposed Riemann-Langevin MC filter was, to the best of our knowledge, the first PF formulated for TBD that exploited measurement information in the proposal. The Riemann-Langevin proposal for TBD was presented along with the particular expressions for the TBD gradient and the Fisher information matrix. A low-dimensional experiment dealing with low-noise, a setup particularly challenging for PF, illustrated that the proposed filter outperformed the considered alternatives. In addition, the experiment demonstrated that the benefit of the Riemann-Langevin approach is not limited to high-dimensional filtering.

The application of the Riemann-Langevin MC filter is of course not limited to TBD. In this regard the Riemann-Langevin proposal and the associated algorithm were formulated generically, enabling their application to a large class of filtering problems, with different transition and measurement models.