Measurement-converted Kalman filter tracking with Gaussian intensity attenuation signal in wireless sensor networks

Abstract

In this article, the target tracking problem in a wireless sensor network with nonlinear Gaussian signal intensity attenuation model is considered. A Bayesian filter tracking algorithm is presented to estimate the locations of moving source that has unknown central signal intensity. This approach adopts a measurement conversion method to remove the measurement nonlinearity by the maximum likelihood estimator, and a linear estimate of the target position and its associated noise statistics obtained by the Newton–Raphson iterative optimization steps are applied into the standard Kalman filter. The Monte Carlo simulations have been conducted in comparison with the commonly used extended Kalman filter with an augmented state that consists of both the original target state and the augmentative central signal intensity. It is observed that the proposed measurement-converted Kalman filter can yield higher accurate estimate and nicer convergence performance over existing methods.

Introduction

Recent advances in micro-electro-mechanical systems, networking systems, and embedded microprocessor technologies have drawn tremendous interests and applications of wireless sensor networks (WSNs). A common issue in WSNs is the source localization, which is inherent to many monitoring applications, such as target tracking, infrastructure monitoring, habitat sensing, rescue, and emergency response.^1–3 The objective of source localization is to sequentially estimate the positions of a moving target by sensing the signal emitted from the target at a subset of distributed network sensors. These sensors collaborate by reporting sensing measurements to a signal processing center which subsequently estimates the source locations according to the received measurement data.

There exist many methods for source localization that are based on measurement models such as time delay of arrival (TDOA), ⁴ direction of arrival (DOA), ⁵ and received signal strength (RSS). ⁶ Generally, the signal intensity measurements are very convenient and economical to locate a moving target ⁷ since no additional sensor functionalities and measurement features are required compared to DOA or TDOA modalities. ⁸ It is known that the intensity of signal attenuates as a function of distance from the source. In the energy-based acoustic source localization applications, the acoustic energy decays exponentially with respect to the distance from source in a free space. ⁹ In the medical ultrasound applications, the echo signal reflects attenuation and scattering properties of tissue that may correlate with various disease states, and the power spectra of received signals are approximately Gaussian. ¹⁰ In the laser manufacturing applications, the nonlinear optical interactions with matter are available from laser beams because of their high power and spectral brightness, and the transverse profile of the optical intensity of the beam can be approximated by a Gaussian function. ¹¹ In the monitoring applications of forest fire, a time-varying elliptical forest fire spread model is adopted to fit to empirical data. This model considers a decreasing intensity Gaussian distribution from the head of fire to the fire tail. ¹² In this article, a type of Gaussian attenuation model is presented to facilitate more accurate characterization of signal intensity. Along this direction, we will consider the source localization problem of a moving target with a Gaussian attenuation model of signal intensity.

Previously, several least-squares methods for energy-based source localization using an isotropic signal energy attenuation model are proposed by Li and Hu ⁹ and Deng and Liu. ¹³ A more generalized statistical model of energy observation is derived by Meesookho et al., ¹⁴ in which an arbitrary acoustic source is characterized by an unknown correlation function, and a weighted direct/one-step least-squares-based algorithm incorporating the dependence of unknown parameters is investigated to reduce the computational complexity. The energy-based sensor network source localization problem is solved by the maximum likelihood (ML) estimation coupled with incremental optimization method in the work by Rabbat and Nowak; ¹⁵ however, this approach cannot converge globally and even diverge due to the problem being highly nonconvex. A new normalized incremental optimization algorithm is developed to avoid trapping into the local minima using a diminishing stepsize in the work by Shi and He. ¹⁶ In the work by Blatt and Hero, ¹⁷ the single-source localization problem is formulated as a convex feasibility problem instead of nonlinear least squares and giving a fast converging localization solution by finding a point at the intersection of some convex sets using projection onto convex sets method. A rough location estimation method is designed to localize an isotropic energy source using measurements from distributed sensors based on kernel averaging techniques. ¹⁸ In the works by Niu and Varshney ¹⁹ and Ozdemir et al., ²⁰ ML-based target localization approaches using the quantized data of signal intensity measurements are proposed to reduce the communication requirements in WSNs and employing a sequential quadratic programming–based grid search within the surveillance region and a reasonable interval of source signal power to find an approximate maximum location. In the work by Masazade et al., ²¹ an iterative source localization based on the approximation of the posterior probability distribution function of source location using a Monte Carlo method is proposed for WSNs with quantized source energy measurements. However, due to the lack of knowledge of the target’s motion, ²² these methods hardly obtain nice localization performance without using prediction and filtering algorithms. In the work by Sheng and Hu, ²³ a sequential source localization method using particle filter is put forward to estimate and track multiple target locations. The large number of particles to represent the non-Gaussian probability density function (pdf) of the target locations involve expensive calculations. A weighted extended Kalman filter (EKF) algorithm based on recursive weighted least-squares optimization for WSNs is presented by Wang and Wu, ²⁴ and the source location is calculated iteratively by taking a weighted average of the local estimates based on the tasking sensor nodes’ reliability. In the work by Mysorewala et al., ²⁵ an EKF-based localization method using an elliptical forest fire spread model is developed, but the linearization procedure of nonlinear measurement function in the EKF algorithm may incur larger errors in the true posterior mean and covariance of the target state.

To overcome the drawback of EKF, many measurement conversion methods ²⁶ have been proposed to transform the nonlinear measurement models into linear ones and estimate the covariances of the converted measurement noises and then use the Kalman filter, which clearly has better accuracy and consistency than the EKF. ²⁷ These conversion methods always require both the range and bearing measurements (and possibly other parameters) of the target, which cannot be applicable to the case where only the signal intensity measurements are available for target tracking. In this article, we propose a new measurement conversion method using ML estimation with Gaussian intensity attenuation model to compute the converted measurement noise. A linear estimate of the target position can be obtained by an iteration that amounts to re-linearization of the measurement equation to reduce the effects of linearization errors, and an approximate covariance matrix of the converted measurement noise is achieved based on Gaussian approximation of the localization error.

We begin by introducing our signal models and formulating the target tracking problem in WSNs in section “Models and problem formulation,” where the signal intensity attenuates Gaussianly with an unknown central signal intensity. We delineate the EKF benchmark that will come handy in the subsequent derivations, and an EKF with augmented state is designed to estimate the augmented state vector that consists of both the original target state and the augmentative central signal intensity parameter in section “Augmented-state EKF tracking algorithm.” In section “Measurement-converted Kalman filtering algorithm,” we lay out the measurement conversion method using ML estimation, and the new converted target position will be followed by a Kalman filter for recursive estimation of the target state. The performance evaluations of the proposed algorithms are conducted by simulation comparisons in section “Numerical simulations.” Finally, some conclusions are given in section “Conclusion.”

Models and problem formulation

Consider an ad hoc WSN that tracks a target in the two-dimensional plane. The state evolution of target is described by a multivariate discrete-time random process^22,28

$${\mathit{X}}_k={\mathit{F}}_k{\mathit{X}}_{k-1}+{\mathit{G}}_k{\mathit{w}}_{k-1}$$

(1)

where ${\mathit{X}}_k\in {\mathbb{R}}^{\aleph }$ is the state vector of the target at the $k\mathrm{th}$ time step, consisting of the position vector ${\mathit{x}}_k=[x(k),y(k){]}^{\mathrm{T}}$ and the velocity vector ${\stackrel{·}{\mathit{x}}}_k=[\stackrel{·}{x}(k),\stackrel{·}{y}(k){]}^{\mathrm{T}}$ . The input noise ${\mathit{w}}_k$ is a Gaussian random vector with zero mean and a covariance matrix ${\mathit{Q}}_k$ ; ${\mathit{F}}_k$ and ${\mathit{G}}_k$ are the state transition matrix and input matrix, respectively.

Signal model

We assume that N identical sensor nodes with the same noise statistics are densely deployed over a sensing field. The position of the $n\mathrm{th}$ ( $1\le n\le N$ ) sensor node, denoted by ${\mathbf{\rho }}_n=[x_n,y_n{]}^{\mathrm{T}}$ , is assumed to be known. The signal emitted from a point source propagates omnidirectionally through the medium, and the intensity of the source will attenuate along the propagation path like an elliptical shape as shown in Figure 1, which can be represented by a Gaussian attenuation model

$$h_n({\mathit{x}}_k)=a\exp \left\{-\frac{1}{2}{({\mathit{x}}_k-{\rho }_n)}^{\mathrm{T}}{\Sigma }^{-1}({\mathit{x}}_k-{\rho }_n)\right\}$$

(2)

where a is the signal intensity at the center of point source and $\mathbf{\Sigma }$ is a variance that determines the signal propagation range and orientation.

Figure 1.

Signal intensity attenuation model represented by a Gaussian function.

Upon sensing a signal intensity greater than a trigger threshold value r, the sensor acts as tasking node to perform tracking task, and the set of indices of tasking sensors at time k can be denoted by

$$L_k=\left\{n|h_n({\mathit{x}}_k)\ge r,1\le n\le N\right\}$$

and then the signal intensity received at the $n\mathrm{th}$ sensor during time k can be expressed as

$$z_n(k)=h_n({\mathit{x}}_k)+v_n(k)n\in L_k$$

(3)

where $v_n(k)$ is the observation noise which is modeled as zero-mean additive white Gaussian noise random variable with variance ${\sigma }_n^2$ and mutually independent with ${\mathit{w}}_k$ and ${\mathit{X}}_k$ . As the source moves in the sensing field, $L_k$ will change with respect to k. Denote ${\ell }_k$ to be the cardinal number of $L_k$ (i.e. the number of elements in $L_k$ ); the sensor measurements may be represented in a matrix form as follows

$$\begin{array}{l}{\mathit{z}}_k={\mathit{h}}_k({\mathit{x}}_k)+{\mathit{v}}_k=\left[\begin{array}{c}{h}_1({\mathit{x}}_k)\\ h_2({\mathit{x}}_k)\\ ⋮\\ h_{{\ell }_k}({\mathit{x}}_k)\end{array}\right]+\left[\begin{array}{c}{v}_1(k)\\ v_2(k)\\ ⋮\\ v_{{\ell }_k}(k)\end{array}\right]\\ =a{\mathit{e}}_k({\mathit{x}}_k)+{\mathit{v}}_k=a\left[\begin{array}{c}{e}_1({\mathit{x}}_k)\\ e_2({\mathit{x}}_k)\\ ⋮\\ e_{{\ell }_k}({\mathit{x}}_k)\end{array}\right]+\left[\begin{array}{c}{v}_1(k)\\ v_2(k)\\ ⋮\\ v_{{\ell }_k}(k)\end{array}\right]\end{array}$$

(4)

where $e_n({\mathit{x}}_k)=\exp \{-1/2({\mathit{x}}_k-{\mathbf{\rho }}_{\mathbf{n}}{)}^{\mathrm{T}}{\mathbf{\Sigma }}^{-1}({\mathit{x}}_k-{\mathbf{\rho }}_n)\},1\le n\le {\ell }_k$ is the Gaussian attenuation function of signal intensity with respect to the $n\mathrm{th}$ sensor, and ${\mathit{e}}_k({\mathit{x}}_k)=[e_1({\mathit{x}}_k),\dots ,e_{{\ell }_k}({\mathit{x}}_k){]}^{\mathrm{T}}$ . The subscripts of h and e in equation (4) refer to the indices within $L_k$ rather than node indices among all N sensor nodes, and then, the covariance matrix of ${\mathit{v}}_k$ is

$${\mathit{R}}_k=\mathrm{diag}\left({\sigma }_1^2,\mathrm{\dots },{\sigma }_{{\ell }_k}^2\right)$$

(5)

The practical signal sources with such attenuation model can be found in ultrasounds, ¹⁰ lasers, ¹¹ and fire sources. ¹²

Problem formulation

For simplicity, we only consider the problem of tracking a single target, but nevertheless our proposed tracking algorithm is also applicable to multi-target tracking as the targets can be classified based on the target signatures. ²⁹ When the target moves through the monitored area, these tasking nodes that have detected the target dynamically form a cluster. One of the nodes in the cluster will be selected as the cluster head, and each cluster member transmits its signal observation to the cluster head. The cluster head serves as the fusion center of the signal and information processing. It selects a set of tasking nodes from the cluster for tracking computation and then broadcasts the estimated target information to its neighbor nodes for the subsequent estimation. ³⁰ More elaborate clustering routing protocols^31–33 can be considered but are beyond the scope of this article.

Let ${\mathit{Z}}_{1:k}$ denote the measurements from initial time up to and including time k, that is, ${\mathit{Z}}_{1:k}=\{{\mathit{z}}_1,\mathrm{\dots },{\mathit{z}}_k\}$ . The problem is for the cluster head to estimate the target state ${\mathit{Z}}_{1:k}=\{{\mathit{z}}_1,\dots ,{\mathit{z}}_k\}$ , denoted by ${\hat{\mathit{X}}}_{k|k}$ , given the measurements ${\mathit{Z}}_{1:k}$ with the unknown central signal intensity parameter a by way of a recursive Bayesian filtering formulation.

Augmented-state EKF tracking algorithm

The basic filtering solution to the state estimation problem can be obtained by a two-stage recursive process of prediction and update. Since the unknown parameter a is not included in the ordinary state variables ${\mathit{X}}_k$ , the clairvoyant EKF cannot be directly applied to the additional parameter estimation case. The augmented model provides a simple method to adjust the state model without altering the filtering framework. In this section, we delineate an augmented-state extended Kalman filter (AEKF) algorithm with the Gaussian signal attenuation model, in which the unknown central signal intensity parameter will be to act on the augmented state, and in turn act on the original model states; then correspondingly the underlying EKF algorithm is not changed and retains most or all of its properties. ³⁴

Let ${\mathit{\theta }}_k=[x(k),\stackrel{·}{x}(k),y(k),\stackrel{·}{y}(k),a{]}^{\mathrm{T}}$ be a new augmented state vector. As a constant, the parameter a can be assumed to vary slowly, and then, the corresponding augmentative evolution model of ${\mathit{\theta }}_{\mathbf{k}}$ is given by

$${\mathit{\theta }}_k={\tilde{\mathit{F}}}_k{\mathit{\theta }}_{k-1}+{\tilde{\mathit{G}}}_k{\tilde{\mathit{w}}}_{k-1}$$

(6)

where ${\tilde{\mathit{w}}}_k~\mathcal{N}(0,{\tilde{\mathit{Q}}}_k)$ , ${\tilde{\mathit{Q}}}_k=\mathrm{diag}({\mathit{Q}}_k,{\sigma }_a^2)$ in which the variance ${\sigma }_a^2$ can be used as a tuning parameter, and

$${\tilde{\mathit{F}}}_k=\left[\begin{array}{ccccc}1& {\mathrm{\Delta }}_k& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1& {\mathrm{\Delta }}_k& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right],{\tilde{\mathit{G}}}_k=\left[\begin{array}{ccc}\frac{{\mathrm{\Delta }}_k^2}{2}& 0& 0\\ {\mathrm{\Delta }}_k& 0& 0\\ 0& \frac{{\mathrm{\Delta }}_k^2}{2}& 0\\ 0& {\mathrm{\Delta }}_k& 0\\ 0& 0& {\mathrm{\Delta }}_k\end{array}\right]$$

In the above equation, ${\mathrm{\Delta }}_k$ is the sampling time interval between two successive time steps k and $k-1$ . Using Bayes’ rule, one has

$$p\left({\mathit{\theta }}_k|{\mathit{Z}}_{1:k}\right)={\beta }_{k}p\left({\mathit{\theta }}_k|{\mathit{z}}_{1:k-1}\right)p\left({\mathit{Z}}_k|{\mathit{\theta }}_k\right)$$

(7)

where the proportionality factor ${\beta }_k^{-1}=p({\mathit{Z}}_k|{\mathit{Z}}_{1:k-1})$ is independent on ${\mathit{\theta }}_k$ .

Assume that the previous estimate of target state ${\hat{\mathit{X}}}_{k-1|k-1}$ and its error covariance matrix ${\mathit{P}}_{k-1|k-1}$ are known at step k. Correspondingly, the augmented state vector estimate at the time $k-1$ can be denoted as ${\hat{\mathit{\theta }}}_{k-1|k-1}=[{\hat{\mathit{X}}}_{k-1|k-1}^{\mathrm{T}},{\hat{a}}_{k-1}{]}^{\mathrm{T}}$ , and its associated error covariance ${\tilde{\mathit{P}}}_{k-1|k-1}=\mathrm{diag}({\mathit{P}}_{k-1|k-1},{\varrho }_{k-1}^2)$ , where ${\hat{a}}_{k-1}$ and ${\varrho }_{k-1}^2$ are, respectively, the estimate and error covariance of unknown parameter a. Now, given the distribution $p({\mathit{\theta }}_{k-1}|{\mathit{Z}}_{1:k-1})$ $=\mathcal{N}({\hat{\mathit{\theta }}}_{k-1|k-1},{\tilde{\mathit{P}}}_{k-1|k-1})$ , the joint distribution of ${\mathit{\theta }}_k$ and ${\mathit{\theta }}_{k-1}$ conditioned on the previous measurements ${\mathit{Z}}_{1:k-1}$ is

$$\begin{array}{l}p({\mathit{\theta }}_{k-1},{\mathit{\theta }}_k|{\mathit{Z}}_{1:k-1})=p({\mathit{\theta }}_k|{\mathit{\theta }}_{k-1})p({\mathit{\theta }}_{k-1}|{\mathit{Z}}_{1:k-1})\\ =\mathcal{N}(\left[\begin{array}{c}{\widehat{\mathit{\theta }}}_{k-1|k-1}\\ {\tilde{\mathit{F}}}_k{\widehat{\mathit{\theta }}}_{k-1|k-1}\end{array}\right],\\ 0.9\left[\begin{array}{cc}{\tilde{\mathit{P}}}_{k-1|k-1}& {\tilde{\mathit{P}}}_{k-1|k-1}{\tilde{\mathit{F}}}_k^{\mathrm{T}}\\ {\tilde{\mathit{F}}}_k{\tilde{\mathit{P}}}_{k-1|k-1}& {\tilde{\mathit{F}}}_k{\tilde{\mathit{P}}}_{k-1|k-1}{\tilde{\mathit{F}}}_k^{\mathrm{T}}+{\tilde{\mathit{G}}}_k{\tilde{\mathit{Q}}}_{k-1}{\tilde{\mathit{G}}}_k^{\mathrm{T}}\end{array}\right])\end{array}$$

(8)

and thus, the predictive distribution of ${\mathit{\theta }}_k$ , given the measurement history up to time step $k-1$ , can be calculated as

$$\begin{array}{l}p\left({\mathit{\theta }}_k|{\mathit{Z}}_{1:k-1}\right)=\mathcal{N}\left({\widehat{\mathit{\theta }}}_{k|k-1},{\tilde{\mathit{P}}}_{k|k-1}\right)\\ =\mathcal{N}\left({\tilde{\mathit{F}}}_k{\widehat{\mathit{\theta }}}_{k-1|k-1},{\tilde{\mathit{F}}}_k{\tilde{\mathit{P}}}_{k-1|k-1}{\tilde{\mathit{F}}}_k^{\mathrm{T}}+{\tilde{\mathit{G}}}_k{\tilde{\mathit{Q}}}_{k-1}{\tilde{\mathit{G}}}_k^{\mathrm{T}}\right)\end{array}$$

(9)

since it is the marginal distribution of ${\mathit{\theta }}_k$ (see “Lemma A2” in Appendix 2). Consider a first-order Taylor series expansion approximation of the measurement function ${\mathit{h}}_k$ around the predictive augmented state ${\hat{\mathit{\theta }}}_{k|k-1}$ in equation (4)

$${\mathit{z}}_k\simeq {\mathit{h}}_k\left({\hat{\mathit{\theta }}}_{k|k-1}\right)+\tilde{\mathit{H}}\left({\hat{\mathit{\theta }}}_{k|k-1}\right)[{\mathit{\theta }}_k-{\hat{\mathit{\theta }}}_{k|k-1}]+{\mathit{v}}_k$$

(10)

where $\tilde{\mathit{H}}({\hat{\mathit{\theta }}}_{k|k-1})=(\partial \mathit{h}({\mathit{x}}_k)/\partial {\mathit{\theta }}_k^{\mathrm{T}}){|}_{{\mathit{\theta }}_k={\hat{\mathit{\theta }}}_{k|k-1}}$ is a ${\ell }_k\times 5$ matrix, and its specific outcome can be found in Appendix 2. Note that ${\mathit{\theta }}_k$ and ${\mathit{v}}_k$ conditioned on the previous measurements ${\mathit{Z}}_{1:k-1}$ are Gaussian and independent, and using equation (10), the joint distribution of ${\mathit{\theta }}_k$ and ${\mathit{z}}_k$ conditioned on ${\mathit{Z}}_{1:k-1}$ can be approximated as

$$\begin{array}{c}p({\mathit{\theta }}_k,{\mathit{z}}_k|{\mathit{Z}}_{1:k-1})=p\left({\mathit{z}}_k|{\mathit{\theta }}_k\right)p\left({\mathit{\theta }}_k|{\mathit{Z}}_{1:k-1}\right)\hfill \\ \multicolumn{1}{c}{\simeq \mathcal{N}(\left[\begin{array}{c}{\hat{\mathit{\theta }}}_{k|k-1}\hfill \\ \multicolumn{1}{c}{{\mathit{h}}_k({\hat{\mathit{\theta }}}_{k|k-1})}\end{array}\right],[\begin{array}{c}{\tilde{\mathit{P}}}_{k|k-1}\hfill \\ \multicolumn{1}{c}{\tilde{\mathit{H}}({\hat{\mathit{\theta }}}_{k|k-1}){\tilde{\mathit{P}}}_{k|k-1}}\end{array}}\\ \multicolumn{1}{c}{\begin{array}{c}{\tilde{\mathit{P}}}_{k|k-1}{\tilde{\mathit{H}}}_k^{\mathrm{T}}({\hat{\mathit{\theta }}}_{k|k-1})\hfill \\ \multicolumn{1}{c}{\tilde{\mathit{H}}({\hat{\mathit{\theta }}}_{k|k-1}){\tilde{\mathit{P}}}_{k|k-1}{\tilde{\mathit{H}}}^{\mathrm{T}}({\hat{\mathit{\theta }}}_{k|k-1})+{\mathit{R}}_k}\end{array}])}\end{array}$$

(11)

From equation (11), it is clear that the posterior distribution of ${\theta }_k$ conditioned on ${\mathit{Z}}_{1:k}$ can be found by the conditional distribution of jointly Gaussian random variables as (see “Lemma A2” in Appendix 2)

$$\begin{array}{c}p\left({\mathit{\theta }}_k|{\mathit{Z}}_{1:k}\right)\simeq \mathcal{N}({\hat{\mathit{\theta }}}_{k|k},{\tilde{\mathit{P}}}_{k|k})\hfill \end{array}$$

(12)

where

$${\hat{\mathit{\theta }}}_{k|k}={\hat{\mathit{\theta }}}_{k|k-1}+{\tilde{\mathit{K}}}_k[{\mathit{z}}_k-{\mathit{h}}_k\left({\hat{\mathit{\theta }}}_{k|k-1}\right)]$$

(13)

$${\tilde{\mathit{P}}}_{k|k}={\tilde{\mathit{P}}}_{k|k-1}-{\tilde{\mathit{K}}}_k\tilde{\mathit{H}}\left({\hat{\mathit{\theta }}}_{k|k-1}\right){\tilde{\mathit{P}}}_{k|k-1}$$

(14)

$$\begin{array}{l}{\tilde{\mathit{K}}}_k={\tilde{\mathit{P}}}_{k|k-1}{\tilde{\mathit{H}}}_k^{\mathrm{T}}\left({\widehat{\mathit{\theta }}}_{k|k-1}\right)\\ {\left[\tilde{\mathit{H}}\left({\widehat{\mathit{\theta }}}_{k|k-1}\right){\tilde{\mathit{P}}}_{k|k-1}{\tilde{\mathit{H}}}^{\mathrm{T}}\left({\widehat{\mathit{\theta }}}_{k|k-1}\right)+{\mathit{R}}_k\right]}^{-1}\end{array}$$

(15)

The above equations lead to the state update ${\hat{\mathit{X}}}_{k|k}$ of the target state ${\mathit{X}}_k$ with the current measurement ${\mathit{z}}_k$ at time k from the augmented state estimate ${\hat{\mathit{\theta }}}_{k|k}$ and its corresponding error covariance ${\mathit{P}}_{k|k}$ from ${\tilde{\mathit{P}}}_{k|k}$ . The recursive process of computing ${\hat{\mathit{X}}}_{k|k}$ is described by the following AEKF in Algorithm 1, where $\mathit{A}[·]$ takes entries of the matrix $\mathit{A}$ .

Algorithm 1: Augmented-state EKF algorithm.
Require: ${\hat{\mathit{X}}}_{0|0}$ , ${\mathit{P}}_{0|0}$ , ${\hat{a}}_0$ and ${\varrho }_0^2$
for: $k=1,2,\cdots$ do
1: Construct augmented state ${\hat{\mathit{\theta }}}_{k-1|k-1}$
2. Construct augmentative error covariance ${\tilde{\mathit{P}}}_{k-1|k-1}$
prediction step: ${\hat{\mathit{\theta }}}_{k|k-1}$ and ${\tilde{\mathit{P}}}_{k|k-1}$
3: Predict augmented state estimate ${\hat{\mathit{\theta }}}_{k|k-1}$ using equation (9)
4: Predict augmentative error covariance ${\tilde{\mathit{P}}}_{k|k-1}$ using equation (9)
update step: ${\hat{\mathit{\theta }}}_{k|k}$ and ${\tilde{\mathit{P}}}_{k|k}$
5: Linearize measurement function ${\mathit{h}}_k({\mathit{x}}_k)$ using equation (10)
6: Update augmented state estimate ${\hat{\mathit{\theta }}}_{k|k}$ by the measurement ${\mathit{z}}_k$ using equation (13)
7: Update augmentative error covariance ${\tilde{\mathit{P}}}_{k|k}$ using equation (14)
Ensure: ${\hat{\mathit{X}}}_{k|k}$ and ${\mathit{P}}_{k|k}$
8: Obtain the target state estimate ${\hat{\mathit{X}}}_{k|k}={\hat{\mathit{\theta }}}_{k|k}[1:4]$
9: Obtain the error covariance of target state ${\mathit{P}}_{k|k}={\tilde{\mathit{P}}}_{k|k}[1:4][1:4]$
end for

Measurement-converted Kalman filtering algorithm

The linearization errors of measurement function approximated by the Taylor series expansion in the EKF algorithm may build up and result in the filtering divergence. ³⁵ In this section, we propose a measurement conversion method to transform the nonlinear measurement model into linear one by a ML estimation method, in which the iterated re-linearization of measurement equation is used to reduce the effects of linearization errors. With the converted linear measurement model, a standard Kalman filter is applied for recursive estimation of the target state.

The conditional pdf of the measurement ${\mathit{z}}_k$ , given ${\mathit{x}}_k$ from equations (4) and (5), can be written as follows

$$\begin{array}{l}p({\mathit{z}}_k|{\mathit{x}}_k)=\mathcal{N}\left({\mathit{h}}_k({\mathit{x}}_k),{\mathit{R}}_k\right)\\ =\frac{1}{{(2\pi )}^{{\ell }_k/2}{(\det \{{\mathit{R}}_k\})}^{1/2}}\\ ·\exp \{-\frac{1}{2}{({\mathit{z}}_k-{\mathit{h}}_k({\mathit{x}}_k))}^{\mathrm{T}}{\mathit{R}}_k^{-1}({\mathit{z}}_k-{\mathit{h}}_k({\mathit{x}}_k))\}\end{array}$$

(16)

The negative log-likelihood function of $p({\mathit{z}}_k|{\mathit{x}}_k)$ is proportional to a quadratic form

$$\begin{array}{cc}L({\mathit{x}}_k)\hfill & =\frac{1}{2}{({\mathit{z}}_k-{\mathit{h}}_k({\mathit{x}}_k))}^{\mathrm{T}}{\mathit{R}}_k^{-1}({\mathit{z}}_k-{\mathit{h}}_k({\mathit{x}}_k))\hfill \end{array}$$

(17)

the above equation ignores the constant terms not dependent on ${\mathit{x}}_k$ , and then, the ML estimate of ${\mathit{x}}_k$ can be obtained by minimizing $L({\mathit{x}}_k)$ . Note that the measurement function ${\mathit{h}}_k({\mathit{x}}_k)$ has an unknown parameter a, and the solution of the optimization problem must lie on a stationary point where

$$\frac{\partial L({\mathit{x}}_k)}{\partial {\mathit{x}}_k}=0,\mathrm{and}\frac{\partial L({\mathit{x}}_k)}{\partial a}=0$$

(18)

The second condition will lead to the following relation of ${\mathit{x}}_k$ and a

$$\begin{array}{c}a=[{\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k){\mathit{R}}_k^{-1}{\mathit{e}}_k({\mathit{x}}_k{)]}^{-1}{\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k){\mathit{R}}_k^{-1}{\mathit{z}}_k\hfill \end{array}$$

(19)

Substituting equation (19) into equation (17), the unknown parameter a can be eliminated, and giving a modified negative log-likelihood function

$$\begin{array}{c}\multicolumn{1}{c}{L\prime ({\mathit{x}}_k)=\frac{1}{2}{\mathit{z}}_k^{\mathrm{T}}(\mathit{I}-{\mathit{e}}_k({\mathit{x}}_k){[{\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k){\mathit{R}}_k^{-1}{\mathit{e}}_k({\mathit{x}}_k)]}^{-1}{\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k){\mathit{R}}_k^{-1}{)}^{\mathrm{T}}{\mathit{R}}_k^{-1}·(\mathit{I}-{\mathit{e}}_k({\mathit{x}}_k){[{\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k){\mathit{R}}_k^{-1}{\mathit{e}}_k({\mathit{x}}_k)]}^{-1}{\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k){\mathit{R}}_k^{-1}){\mathit{z}}_k}\\ =\frac{1}{2}{\mathit{z}}_k^{\mathrm{T}}{\mathit{R}}_k^{-1}{\mathit{z}}_k-\frac{1}{2}{\mathit{z}}_k^{\mathrm{T}}{\mathit{R}}_k^{-1}{\mathit{e}}_k({\mathit{x}}_k){[{\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k){\mathit{R}}_k^{-1}{\mathit{e}}_k({\mathit{x}}_k)]}^{-1}{\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k){\mathit{R}}_k^{-1}{\mathit{z}}_k\hfill \end{array}$$

(20)

The minimization problem of equation (20) can be equivalently stated as a nonlinear unconstrained optimization problem, and the common solution to the optimization problem is found through the following Newton–Raphson iterative method ³⁶

$${\widehat{\mathit{x}}}_{k|k}^t={\widehat{\mathit{x}}}_{k|k}^{t-1}-{\left[\frac{{\partial }^{2}L\prime ({\mathit{x}}_k)}{\partial {\mathit{x}}_k^{\mathrm{T}}\partial {\mathit{x}}_k}{|}_{{\mathit{x}}_k={\widehat{\mathit{x}}}_{k|k}^{t-1}}\right]}^{-1}\frac{\partial L\prime ({\mathit{x}}_k)}{\partial {\mathit{x}}_k}{|}_{{\mathit{x}}_k={\widehat{\mathit{x}}}_{k|k}^{t-1}}$$

(21)

where t is the iteration step, and the initial value ${\hat{\mathit{x}}}_{k|k}^0={\hat{\mathit{x}}}_{k|k-1}$ . The ML estimation problem equivalently is a nonlinear least-squares problem $L\prime ({\mathit{x}}_k)=\frac{1}{2}{\gamma }^{\mathrm{T}}({\mathit{x}}_k)\gamma ({\mathit{x}}_k)$ , where $\gamma ({\mathit{x}}_k)={\mathit{R}}_k^{-1/2}(\mathit{I}-{\mathit{e}}_k({\mathit{x}}_k)[{\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k){\mathit{R}}_k^{-1}{\mathit{e}}_k({\mathit{x}}_k){]}^{-1}{\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k){\mathit{R}}_k^{-1}){\mathit{z}}_k$ . The first-order Taylor series expansion of measurement function in the standard EKF is consistent with the Gauss–Newton method to solve the nonlinear least-squares problem. ³⁷ Unlike Newton–Raphson method, the Gauss–Newton algorithm approximates the Hessian matrix of $L\prime ({\mathit{x}}_k)$ by ignoring the second-order derivative terms of the residuals ${\gamma }^{\mathrm{T}}({\partial }^2\gamma /\partial {\mathit{x}}_k[j]\partial {\mathit{x}}_k[i])$ . As a consequence, the convergence of the Gauss–Newton method depends on whether the omitted second-order derivative terms of the residuals are large parts of the Hessian. When

$$\left|\frac{\partial {\gamma }^{\mathrm{T}}}{\partial {\mathit{x}}_k[j]}\frac{\partial {\gamma }^{\mathrm{T}}}{\partial {\mathit{x}}_k[i]}\right|>>\left|{\gamma }^T\frac{{\partial }^2\gamma }{\partial {\mathit{x}}_k[j]\partial {\mathit{x}}_k[i]}\right|$$

cannot be satisfied, it is shown that the Gauss–Newton method may not be locally convergent at all. ³⁸

In equation (21), the $i\mathrm{th}$ component of the gradient of the logarithm function $L\prime ({\mathit{x}}_k)$ with respect to ${\mathit{x}}_k[i]$ is easily given as

$$\begin{array}{l}\frac{\partial L\prime ({\mathit{x}}_k)}{\partial {\mathit{x}}_k[i]}={\left[{\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k){\mathit{R}}_k^{-1}{\mathit{e}}_k({\mathit{x}}_k)\right]}^{-2}{\left[{\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k){\mathit{R}}_k^{-1}{\mathit{z}}_k\right]}^2\\ \frac{\partial {\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k)}{\partial {\mathit{x}}_k[i]}{\mathit{R}}_k^{-1}{\mathit{e}}_k({\mathit{x}}_k)-{\left[e_k^{\mathrm{T}}({\mathit{x}}_k){\mathit{R}}_k^{-1}{\mathit{e}}_k({\mathit{x}}_k)\right]}^{-1}\\ \left[{\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k){\mathit{R}}_k^{-1}{\mathit{z}}_k\right]\frac{\partial {\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k)}{\partial {\mathit{x}}_k[i]}{\mathit{R}}_k^{-1}{z}_k\end{array}$$

(22)

Then, the element at the $i\mathrm{th}$ row and $j\mathrm{th}$ column of the Hessian matrix of the function $L\prime ({\mathit{x}}_k)$ can be computed as

$$\begin{array}{c}\frac{{\partial }^{2}L\prime ({\mathit{x}}_k)}{\partial {\mathit{x}}_k[j]\partial {\mathit{x}}_k[i]}={[{\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k){\mathit{R}}_k^{-1}{\mathit{e}}_k({\mathit{x}}_k)]}^{-2}{[{\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k){\mathit{R}}_k^{-1}{\mathit{z}}_k]}^2\hfill \\ \frac{{\partial }^2{\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k)}{\partial {\mathit{x}}_k[j]\partial {\mathit{x}}_k[i]}{\mathit{R}}_k^{-1}{\mathit{e}}_k({\mathit{x}}_k)\hfill \\ +{[{\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k){\mathit{R}}_k^{-1}{\mathit{e}}_k({\mathit{x}}_k)]}^{-2}{[{\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k){\mathit{R}}_k^{-1}{\mathit{z}}_k]}^2\hfill \\ \frac{\partial {\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k)}{\partial {\mathit{x}}_k[i]}{\mathit{R}}_k^{-1}\frac{\partial {\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k)}{\partial {\mathit{x}}_k[j]}\hfill \\ -{[{\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k){\mathit{R}}_k^{-1}{\mathit{e}}_k({\mathit{x}}_k)]}^{-1}[{\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k){\mathit{R}}_k^{-1}{\mathit{z}}_k]\frac{{\partial }^2{\mathit{e}}_k^{\mathrm{T}}({\mathit{x}}_k)}{\partial {\mathit{x}}_k[j]\partial {\mathit{x}}_k[i]}{\mathit{R}}_k^{-1}{\mathit{z}}_k\hfill \end{array}$$

(23)

The detailed derivation of equation (23) is put in the Appendix 2.

Considering a second-order Taylor series expansion approximation of $L({\mathit{x}}_k)$ around ${\hat{\mathit{x}}}_{k|k}^t$ in equation (17) and combining with equation (18) yield

$$\begin{array}{l}L({\mathit{x}}_k)\simeq L\prime ({\widehat{\mathit{x}}}_{k|k}^t)+\frac{1}{2}{\left({\mathit{x}}_k-{\widehat{\mathit{x}}}_{k|k}^t\right)}^{\mathrm{T}}\\ {\left[\frac{{\partial }^{2}L\prime ({\mathit{x}}_k)}{\partial {\mathit{x}}_k^{\mathrm{T}}\partial {\mathit{x}}_k}{|}_{{\mathit{x}}_k={\widehat{\mathit{x}}}_{k|k}^t}\right]}^{-1}\left({\mathit{x}}_k-{\widehat{\mathit{x}}}_{k|k}^t\right)\end{array}$$

(24)

Using the above approximation, it can follow from equations (16) and (17) that

$$\begin{array}{l}p\left({\widehat{\mathit{x}}}_{k|k}^t|{\mathit{x}}_k\right)\simeq \alpha \exp \\ \left\{-\frac{1}{2}{\left({\mathit{x}}_k-{\widehat{\mathit{x}}}_{k|k}^t\right)}^{\mathrm{T}}{\left[\frac{{\partial }^{2}L\prime ({\mathit{x}}_k)}{\partial {\mathit{x}}_k^{\mathrm{T}}\partial {\mathit{x}}_k}{|}_{{\mathit{x}}_k={\widehat{\mathit{x}}}_{k|k}^t}\right]}^{-1}\left({\mathit{x}}_k-{\widehat{\mathit{x}}}_{k|k}^t\right)\right\}\end{array}$$

(25)

for some constant $\alpha$ . That is, the estimate error ${\mathit{x}}_k-{\hat{\mathit{x}}}_{k|k}^t$ of the position state ${\mathit{x}}_k$ can be approximated as a Gaussian distribution with zero mean and covariance matrix

$${\mathbf{\sigma }}_{\epsilon }^2(k)={\frac{{\partial }^{2}L\prime ({\mathit{x}}_k)}{\partial {\mathit{x}}_k^{\mathrm{T}}\partial {\mathit{x}}_k}|}_{{\mathit{x}}_k={\hat{\mathit{x}}}_{k|k}^t}$$

(26)

which is described by the following linear representation in the target state ${\mathit{X}}_k$

$${\hat{\mathit{x}}}_{k|k}^t=\left[\begin{array}{cccc}1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill \end{array}\right]{\mathit{X}}_k+{\epsilon }_k=\mathbf{C}{\mathit{X}}_k+{\epsilon }_k$$

(27)

where ${\epsilon }_k~\mathcal{N}(0,{\mathbf{\sigma }}_{\mathbf{\epsilon }}^2(k))$ . We consider the estimate of the position state ${\hat{\mathit{x}}}_{k|k}^t$ as a new converted measurement and then it can be found that the converted measurement function $\mathbf{C}$ is linear, correspondingly ${\epsilon }_k$ is the converted measurement noise. The converted measurement model will be used in Kalman filtering for the recursive update of the target state ${\mathit{X}}_k$ .

From equation (9), the predictive distribution of ${\mathit{X}}_k$ , given the measurement history up to time step $k-1$ , is

$$\begin{array}{c}p\left({\mathit{X}}_k|{\mathit{Z}}_{1:k-1}\right)=\mathcal{N}({\hat{\mathit{X}}}_{k|k-1},{\mathit{P}}_{k|k-1})\hfill \\ \multicolumn{1}{c}{=\mathcal{N}({\mathit{F}}_k{\hat{\mathit{X}}}_{k-1|k-1},{\mathit{F}}_k{\mathit{P}}_{k-1|k-1}{\mathit{F}}_k^{\mathrm{T}}+{\mathit{G}}_k{\mathit{Q}}_{k-1}{\mathit{G}}_k^{\mathrm{T}})}\end{array}$$

(28)

By Bayes’ rule, and using the fact ${\hat{\mathit{x}}}_{k|k}^t$ is independent with ${\mathit{Z}}_{1:k-1}$ given ${\mathit{X}}_k$ , since ${\hat{\mathit{x}}}_{k|k}^t$ is generated from ${\mathit{z}}_k$ , then we have the following joint probability distribution of $({\mathit{X}}_k,{\hat{\mathit{x}}}_{k|k}^t)$ , and

$$\begin{array}{c}p({\mathit{X}}_k,{\hat{\mathit{x}}}_{k|k}^t|{\mathit{Z}}_{1:k-1})=p\left({\hat{\mathit{x}}}_{k|k}^t|{\mathit{X}}_k\right)p\left({\mathit{X}}_k|{\mathit{Z}}_{1:k-1}\right)\hfill \\ \multicolumn{1}{c}{\simeq \mathcal{N}(\left[\begin{array}{c}{\hat{\mathit{X}}}_{k|k-1}\hfill \\ \multicolumn{1}{c}{\mathbf{C}{\hat{\mathit{X}}}_{k|k-1}}\end{array}\right],\left[\begin{array}{cc}{\mathit{P}}_{k|k-1}\hfill & {\mathit{P}}_{k|k-1}{\mathbf{C}}^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{\mathbf{C}{\mathit{P}}_{k|k-1}}& \mathbf{C}{\mathit{P}}_{k|k-1}{\mathbf{C}}^{\mathrm{T}}+{\mathbf{\sigma }}_{\mathbf{\epsilon }}^2(k)\hfill \end{array}\right])}\end{array}$$

(29)

thus, the posterior distribution of ${\mathit{X}}_k$ conditioned on the history measurements ${\mathit{Z}}_{1:k-1}$ and the current converted measurement ${\hat{\mathit{x}}}_{k|k}^t$ can be obtained in the same way as for the AEKF in equation (12)

$$\begin{array}{c}p({\mathit{X}}_k|{\mathit{Z}}_{1:k-1},{\hat{\mathit{x}}}_{k|k}^t)\simeq \mathcal{N}({\hat{\mathit{X}}}_{k|k},{\mathit{P}}_{k|k})\hfill \end{array}$$

(30)

where

$${\hat{\mathit{X}}}_{k|k}={\hat{\mathit{X}}}_{k|k-1}+{\mathit{K}}_k[{\hat{\mathit{x}}}_{k|k}^t-\mathbf{C}{\hat{\mathit{X}}}_{k|k-1}]$$

(31)

$${\mathit{P}}_{k|k}={\mathit{P}}_{k|k-1}-{\mathit{K}}_k\mathbf{C}{\mathit{P}}_{k|k-1}$$

(32)

$${\mathit{K}}_k={\mathit{P}}_{k|k-1}{\mathbf{C}}^{\mathrm{T}}{[\mathbf{C}{\mathit{P}}_{k|k-1}{\mathbf{C}}^{\mathrm{T}}+{\mathbf{\sigma }}_{\mathbf{\epsilon }}^2(k)]}^{-1}$$

(33)

The recursive computing process of the measurement-converted Kalman filter (MCKF) is summarized in the following Algorithm 2.

Algorithm 2: Measurement-converted Kalman filter.
Require: ${\hat{\mathit{X}}}_{0|0}$ and ${\mathit{P}}_{0|0}$
for: $k=1,2,\cdots$ do
1: Compute the converted measurement ${\hat{\mathit{x}}}_{k|k}^t$ using equation (21)
2: Compute the converted noise covariance ${\sigma }_{\epsilon }^2(k)$ using equation (26)
prediction step: ${\hat{\mathit{X}}}_{k|k-1}$ and ${\mathit{P}}_{k|k-1}$
3: Predict state estimate ${\hat{\mathit{X}}}_{k|k-1}$ using equation (28)
4: Predict error covariance ${\mathit{P}}_{k|k-1}$ using equation (28)
update step: ${\hat{\mathit{X}}}_{k|k}$ and ${\mathit{P}}_{k|k}$
5: Update state estimate ${\hat{\mathit{X}}}_{k|k}$ by the converted measurement ${\hat{\mathit{x}}}_{k|k}^t$ using equation (31)
6: Update error covariance ${\mathit{P}}_{k|k}$ using equation (32)
end for

We remark that the AEKF is initialized by additionally guessing the initial estimate ${\hat{a}}_0$ and error covariance ${\varrho }_0^2$ of unknown parameter a, so that its convergence performance is extremely sensitive to the initialization accuracy. In contrast the MCKF not only has no additional assumption on the initial state of unknown parameter a, and the nice convergence can be partly assisted by the fact that the initial estimate of the iterative process of measurement conversion comes from the Kalman predictor which is typically good, as such the sampling time interval is not too long for convergence. Typically, only a few iterations are sufficient to get the extremum ³⁹ because the longer iteration steps cannot promise more accurate estimate.

Numerical simulations

To show the efficiency of the proposed MCKF, it is applied to a target tracking application in comparison with the AEKF, the traditional extended Kalman filter (TEKF) without measurement conversion, and the pure ML estimator without Kalman filtering. Error performance of the filters is evaluated with Monte Carlo simulations in this section.

The monitored area is 50 m × 50 m and covered by $N=50$ sensor nodes randomly, the trigger threshold of each sensor $r=0.1$ , and the measurement noise covariance of each sensor ${\sigma }_n^2=5$ . A moving target travels along a circular arc, and the sampling interval ${\mathrm{\Delta }}_k=0.1\mathrm{s}$ (10 Hz). The process noise ${\mathit{w}}_k$ corresponds to the variable acceleration of the target and can be approximated by a white Gaussian sequence with zero mean and covariance matrix ${\mathit{Q}}_k=\mathrm{diag}(10,10)$ . The unknown parameter a is assumed to vary slowly in the AEKF, and its tuning covariance is set to ${\sigma }_a^2=0.00001$ ; in the TEKF, a will be approximately computed by the ML estimator before the state update step.

In the experiments, the source central signal intensity is set at $a=25$ , and $\mathbf{\Sigma }=\left[\begin{array}{cc}40\hfill & 10\hfill \\ \multicolumn{1}{c}{5}& 30\hfill \end{array}\right]$ . The initial state estimate and the corresponding covariance matrix for the MCKF and TEKF, respectively, are chosen to be

$$\begin{array}{c}{\hat{\mathit{X}}}_{0|0}=[\begin{array}{cccc}0\mathrm{m}\hfill & 5\mathrm{m}/\mathrm{s}\hfill & 25\mathrm{m}\hfill & 5\mathrm{m}/\mathrm{s}\hfill \end{array}{]}^{\mathrm{T}}\hfill \\ {\mathit{P}}_{0|0}=0.01\times \mathrm{diag}(\begin{array}{cccc}1\hfill & 1\hfill & 1\hfill & 1\hfill \end{array})\hfill \end{array}$$

Additionally, the initial estimate of the variable a in the AEKF is roughly assumed to the trigger threshold of the sensor ${\hat{a}}_0=r=0.1$ , and the associated variance is set to a large value ${\varrho }_0^2=2$ so that the estimate converges quickly and the influence of the initial guess ${\hat{a}}_0$ soon will be negligible.

In total, 100 Monte Carlo runs with the above noise condition are performed. The results are summarized in the following. Figure 2 shows the tracking trajectories of the four tracking algorithms over a typical simulation run. Figure 3 describes the corresponding mean value of the population of target estimated location biases $\parallel {\mathit{x}}_k-{\hat{\mathit{x}}}_{k|k}\parallel$ for each tracking algorithm, where ${\hat{\mathit{x}}}_{k|k}$ comes from the estimated state ${\hat{\mathit{X}}}_{k|k}$ of the target in the MCKF and TEKF; in the AEKF, ${\hat{\mathit{x}}}_{k|k}\in {\hat{\mathit{\theta }}}_{k|k}$ , and in the ML and MCKF, $t=2$ . It is quite clear that the MCKF has the best tracking performance for most time steps, while the ML performs worse in comparison to the filters, and this is because the ML estimation method only utilizes the measurement information but the other filtering methods additionally utilize the target moving information. Figure 4 further depicts the average estimated location bias of the whole tracking trajectory in each realization, and the average estimated location biases of the ML are also more than the MCKF, the AEKF as well as the TEKF; the average estimated location biases of the TEKF are always greater than the MCKF and AEKF. These results indicate that the ML hardly obtains nice estimation accuracy without using prediction and filtering algorithms.

Figure 2.

Tracking trajectories of the target under different tracking algorithms: the target travels along a circular arc in a WSN of 50 nodes.

Figure 3.

The average estimated location biases of 100 realizations under different tracking algorithms: the solid line with marker represents the mean value of population of target estimated location biases for all the time steps.

Figure 4.

The average estimated location bias of the whole tracking trajectory in each realization: the solid line with marker represents the average estimated location bias for all the time steps.

Figure 5 gives the mean of quadratic sums of the position state estimation errors in each time step for the four methods, and the distributions of tracking errors of the position states along x- and y-directions together with the corresponding error covariance ellipses of position state are plotted in Figure 6. The error covariance of the position state is defined as

$${\mathit{P}}_{k|k,\mathrm{xy}}=\left[\begin{array}{cc}{\mathit{P}}_{k|k}[1,1]\hfill & {\mathit{P}}_{k|k}[1,3]\hfill \\ \multicolumn{1}{c}{{\mathit{P}}_{k|k}[3,1]}& {\mathit{P}}_{k|k}[3,3]\hfill \end{array}\right]$$

where ${\mathit{P}}_{k|k}[i,i]$ denotes the diagonal element of the error covariance matrix ${\mathit{P}}_{k|k}$ . The major and minor axes of each ellipse correspond to the square roots of eigenvalues of the error covariance matrix, and the angles between these axes and x-axis are determined by the orientation of corresponding eigenvectors. We can see that the mean-squared error (MSE) of the position state estimation in the ML is the worst among all the algorithms. The MSE of the MCKF is lower than the AEKF and TEKF, and this is because the MCKF can achieve smaller error covariance among them, as such the tracking errors of the position states in the MCKF for most realizations lie in the ellipse. In contrast, the MSE of the TEKF is larger than the AEKF and MCKF, and this is because the estimation error statistics of unknown central signal intensity parameter are not considered in the TEKF, which would reduce the filtering performance. This result indicates that the MCKF based on the measurement conversion method yields significant performance over the TEKF and AEKF based on the Taylor series linearization of measurement function.

Figure 5.

The mean-squared errors of the position state estimation under different tracking algorithms: the solid line with marker describes the mean value of the MSE population of 100 realizations.

Figure 6.

The tracking errors along x- and y-directions and the error covariance ellipses of the position state: the sample data points are the population of 100 realizations of $x(k)-{\hat{x}}_{k|k}$ and $y(k)-{\hat{y}}_{k|k}$ at some random time steps.

Figure 7 shows the average estimate of the source central signal intensity a of 100 realizations at each time step for each tracking algorithm. Specially, the average estimates of a in the AEKF are given for three different initial state cases of a. As shown in Figure 7, the convergence rate of the estimate $\hat{a}$ in the AEKF depends on the initial state of unknown parameter a, and when the error in the initial state estimate is consistent with the initial state covariance in the case that ${\hat{a}}_0=0.1,{\varrho }_0^2=2$ , the AEKF will have a good convergence performance. In contrast, the convergence rate of the estimate $\hat{a}$ in the AEKF is lower than the MCKF, the TEKF, and the ML, this is because that the ML, TEKF, and the MCKF directly obtain better estimate of $\mathit{a}$ from the current measurements using the ML estimation method.

Figure 7.

The average estimate of a of 100 realizations: the black solid line without marker represents the real source signal intensity a, and the average estimates of a in the AEKF are given for three different initial state cases of a: (a) ${\hat{a}}_0=0.1,{\varrho }_0^2=2$ , (b) ${\hat{a}}_0=100,{\varrho }_0^2=2$ , and (c) ${\hat{a}}_0=0.1,{\varrho }_0^2=0.02$ .

To investigate the stability of the proposed MCKF, another simulation is conducted that uses the same setting as above except for different measurement noise and monitoring range conditions as well as different initial state cases of a. Table 1 lists the average root mean square of the target estimated position error (RMSE) $\sqrt{\mathbb{E}\{\parallel {\mathit{x}}_k-{\hat{\mathit{x}}}_{k|k}{\parallel }^2\}}$ in all the time steps under different scenarios. As r increases, the monitoring range of each sensor will diminish, and then, the number of sensors that detect the target decreases and the RMSEs of these tracking algorithms become higher. In the low measurement noise condition, the impact of the noise has negligible magnitude and then the observations are very accurate; thus, all the tracking algorithms attain the same RMSE. However, under the high measurement noise conditions, the measurement errors become bigger, the RMSE of the ML becomes the worst, the RMSE of the AEKF increases as the initial covariance of a decreases, but the RMSE of the proposed MCKF is still lower than the AEKF and TEKF. These results indicate that the MCKF can yield higher accurate estimate and nicer convergence performance than the AEKF, TEKF, and the ML.

Table 1.

Average RMSEs under different measurement noise and monitoring range conditions with different initial state cases of a.

Measurement noise	Average RMSE
	ML	MCKF	TEKF	AEKF ${\varrho }_0^2=0.02$	AEKF ${\varrho }_0^2=0.2$	AEKF ${\varrho }_0^2=2$
${\sigma }_n^2=0.1,r=0.1$	0.0936	0.0737	0.0760	0.0799	0.0754	0.0749
${\sigma }_n^2=0.1,r=5$	0.1011	0.0786	0.0804	0.0853	0.0799	0.0789
${\sigma }_n^2=1,r=0.1$	0.2923	0.1866	0.2183	0.2557	0.1972	0.1889
${\sigma }_n^2=1,r=5$	0.3165	0.1986	0.2276	0.2718	0.2098	0.2010
${\sigma }_n^2=5,r=0.1$	0.6560	0.3476	0.4293	0.7432	0.4101	0.3572
${\sigma }_n^2=5,r=5$	0.7095	0.3694	0.4408	0.8014	0.4357	0.3791

RMSE: root mean square of the target estimated position error; ML: maximum likelihood; MCKF: measurement-converted Kalman filter; TEKF: traditional extended Kalman filter; AEKF: augmented-state extended Kalman filter.

Conclusion

In this article, we consider a practical Gaussian attenuation model to facilitate more accurate characterization of source signal intensity. With the nonlinear Gaussian intensity attenuation model, two Bayesian filtering tracking algorithms (AEKF and MCKF) are developed. The Taylor series linearization error of nonlinear measurement function in the AEKF would lead to unstable tracking performance, and its convergence performance is extremely sensitive to the initialization accuracy. The measurement conversion method in the MCKF is used to remove the measurement nonlinearity using the ML estimation method, and a linear estimate of the target position is obtained by an iteration that amounts to re-linearization of the measurement equation to reduce the effects of linearization errors. The converted measurement and its associated noise statistics are then used in a standard Kalman filter for recursive update of the target state. Simulation results have shown that the MCKF can yield higher accurate estimate and nicer convergence performance than the AEKF and the TEKF as well as the commonly used ML estimation approach. The discussion in this work is limited to the case of single target. Multiple target scenarios will be addressed in future works.