Passive Underwater Target Tracking: Conditionally Minimax Nonlinear Filtering with Bearing-Doppler Observations

Abstract

The paper presents an application of the Conditionally-Minimax Nonlinear Filtering (CMNF) algorithm to the online estimation of underwater vehicle movement given a combination of sonar and Doppler discrete-time noisy sensor observations. The proposed filter postulates recurrent “prediction–correction” form with some predefined basic prediction and correction terms, and then they are optimally fused. The CMNF estimates have the following advantageous features. First, the obtained estimates are unbiased. Second, the theoretical covariance matrix of CMNF errors meets the real values. Third, the CMNF algorithm gives a possibility to choose the preliminary observation transform, basic prediction, and correction functions in any specific case of the observation system to improve the estimate accuracy significantly. All the features of conditionally-minimax estimates are demonstrated by the regression example of random position estimate given the noisy bearing observations. The contribution of the paper is the numerical study of the CMNF algorithm applied to the underwater target tracking given bearing-only and bearing-Doppler observations.

1. Introduction

The class of underwater target tracking problems is rapidly widening. This process has been provoked not only by the dramatic growth of the underwater vehicle industry, resulting in a big number of multi-functional vehicles available, but also by the active development of the registration devices and improvement of their usage policies along with advanced methods of statistical data processing. The spectrum of applications is also broad and lies both in the civil [1,2,3,4,5] and military areas [6,7,8,9,10].

Without aiming to write an exhaustive review in the area of the underwater target tracking, we should mention the following works [11,12,13]. The analysis of these comprehensive contributions together with other related literature allows for determining the factors and restrictions which affect the performance in any specific applied tracking problem:

–

available sensors configuration and their characteristics,

–

the mathematical model of the couple “target dynamics–measurements”,

–

requirements to the target state estimates,

–

theoretical framework and numerical algorithms, involved into the solution to the tracking problem.

The object of this research is the long-range passive tracking problem of a maneuvering target by stationary multi-static sonar site [11]. This task looks rather challenging. First, in this case, the available observations are the bearings towards the target with or without the Doppler shift of the received signal caused by the radial target movement relative to the sonar. The target range is unobservable in the passive mode, which complicates the tracking significantly. Moreover, the bearing measurements of sonars are considerably less accurate in comparison with radars: angular errors of sonars are typically $1^{°}$–$3^{°}$ [14], while those of radars vary in the range of ${0.05}^{°}$–${0.3}^{°}$ [15]. Second, the negative feature of the sonar observations is the port-starboard ambiguity. The choice of the correct target position is nontrivial and can be done by various methods [6,16,17]. Third, at the moment of the distant target detection, the statistical uncertainty of the position is huge; this forces applying advanced mathematical methods to make the state estimate quality reasonable. Fourth, usually, the moving target is maneuvering [18], and the parameters of the maneuvers are also subject to estimation.

From the mathematical point of view, the target tracking is a filtering problem. The unobservable state is described by a stochastic differential system (SDS), and this state has to be estimated given the indirect discrete-time nonlinear noisy observations. The optimality criterion is usually the Mean Square (MS) of the state estimate error. For practical realization of the filtering process, the common practice is to apply either nonlinear modifications of the well-known Kalman filtering (KF) algorithm: linearized/extended Kalman filter (LKF/EKF) [19,20,21,22], method of pseudo measurements [23,24,25], unscented Kalman filter (UKF) [26,27], quadrature or cubature Kalman filter [28,29], and numerous versions of the particle filter [30,31]. All these methods have known issues: the nonlinear versions of the KF may diverge in specific settings [32,33,34] or provide the biased estimates [23,24], quadrature, and cubature filters are resource-demanding for high-dimensional observation systems [28], and so are the particle filters.

The aim of this paper is the presentation of the Conditionally-Minimax Nonlinear Filtering (CMNF) algorithm for underwater target tracking in the conditions stated above. This algorithm was initially proposed for the state filtering of the discrete-time stochastic finite-dimensional [35] and distributed [36] observation systems. Recently, it was successfully applied to the tracking of a maneuvering aerial target given the radar observations [37] and the navigation of the Automated Underwater Vehicle (AUV) [38].

The idea of the CMNF algorithm is quite transparent. We fix the structure of the filter in the two-step “prediction–correction” form. The transform of the estimate on the previous step (so-called basic prediction function) and the observations on the present step (observation transform and basic correction function) are predefined. The optimization procedure deals only with optimal fusion of these functions in the final filtering estimate in order to provide the best performance in the MS-sense. The minimax property of this fusion is related closely with the one of the pair “Gaussian distribution— linear estimator”. If the mutual distribution of the triplet “state–prediction term–correction term” is uncertain with known mean and covariance matrix bounded from above, then the least favorable distribution is Gaussian, and then the best (equalizing) estimate is the linear one [35,36].

The CMNF estimate has remarkable properties. First, it is unbiased, i.e., the CMNF estimate error is a centered random vector. Second, the CMNF algorithm allows for calculating the upper bound of the error covariance matrix, and this conservative choice guarantees that the algorithm is non-divergent in contrast with the EKF and UKF. Third, the CMNF algorithm gives a possibility to select the functions of basic prediction, observation transform, and correction in accordance with additional a priori information, such as the explicit form of the state dynamics, some guess concerning the form of the optimal nonlinear filter for the observation system at hand, etc. It turns out that a successful choice provides a significant gain in estimation accuracy.

Another aim of the research is to analyze the contribution of Doppler observations into the CMNF performance in comparison with the bearing-only observation version.

The paper is organized as follows. Section 2 contains a description of the target tracking problem. We present the stochastic dynamic system, which defines the target motion, and the model of the sensor measurements.

Section 3 is devoted to the theoretical background of the conditionally-minimax estimation. Traditionally, the key features and performance of the nonlinear KF algorithms like EKF, UKF, etc. are illustrated by a regression example, hence Section 3.1 contains the formulation and solution for the minimax estimation problem of the regression parameters. The mutual distribution of the estimated and observable components is uncertain: their mathematical expectation is known, and the uncertainty set of covariance matrices has a maximal element. An assertion establishing the pair “Gaussian distribution–linear estimator” as a saddle point of the optimality criterion at hand is provided. If the observable part is complemented by some additional transforms of the observations, then it is possible to outperform the accuracy of the best linear estimate. This feature is demonstrated in the regression example. In Section 3.2, we present a conditionally-minimax filtering problem statement for the stochastic continuous-discrete observation systems along with its solution.

Section 4 is the main part of the paper: it contains the comparative material of two numerical experiment groups. Section 4.1 is devoted to the estimation of the three-dimensional random vector given noisy bearing observations. We compare the performance of the best linear and conditionally-minimax estimates. The last involves the least-squares (LS) estimates of the estimated vector as additional pseudo observations. Section 4.2 contains a comparison of the tracking results of a maneuvering target by the noisy bearing-Doppler observations. The maneuver being a planar stochastic coordinate turn with unobservable random target speed and oscillating normal acceleration is described by 17 equations. The parameters of the maneuver plane are also random and unobservable. The concluding remarks are presented in Section 5.

2. Model Description and Tracking Problem Formulation

2.1. Target Dynamics

In what follows, we will operate two coordinate systems. The first one is the coordinate system in which the target motion can be defined as a planar maneuver, and the second one is the coordinate system in which the observations are collected and the estimate is calculated. We call the former “target coordinate system” $0xyz$ and the latter—“observer coordinate system” $0XYZ$. Both systems are stationary, and the translation from one system to another is given by two Euler angles $\alpha$ and $\beta$ and the shift vector ${\mathbf{R}}^0$ (see Figure 1).

The planar target motion in the target coordinate system $0xyz$ is defined by the following SDS:

$$\left\{\begin{array}{c}d{x}_t=v_{t}cos{\phi }_{t}dt,\hfill \\ d{y}_t=v_{t}sin{\phi }_{t}dt,\hfill \\ d{z}_t=0,\hfill \\ d{v}_t=0,\hfill \\ d{\phi }_t=\frac{a_t^n}{v_t}dt,\hfill \\ d{a}_t^n=(-\lambda a_t^n+\nu )dt+\mu d{W}_t,\hfill \end{array}\right.$$

(1)

where $x_t$, $y_t$, $z_t$ are the coordinates, $v_t$ is the absolute value of the target speed, which is constant, ${\phi }_t$ is the heading angle, $a_t^n$ is the normal acceleration, and $W_t$ is a standard Wiener process. The acceleration equation parameters $\lambda$ and $\nu$, as well as the disturbance intensity $\mu$, are assumed to be known. In SDS (1), the equations for $x_t$, $y_t$, $v_t$, and ${\phi }_t$ follow from the standard planar curvilinear-motion kinematic model for constant heading speed motion with the normal acceleration $a_t^n$, whose dynamics are given by the last equation. In [18], this model is described in detail (see Equations (55), (57)–(60)); the only difference is that here we use the stochastic differential notation (and hence the Wiener process instead of the white Gaussian noise).

The target starts from the origin, so the initial condition for the coordinates is ${(x_0,y_0,z_0)}^T={(0,0,0)}^T$, but other initial parameters $v_0$, ${\varphi }_0$, $a_0^n$ are assumed to be random with known distributions.

In the observer coordinate system, the target has the coordinates ${\mathbf{X}}_t={(X_t,Y_t,Z_t)}^T$ and the vector speed ${\mathbf{V}}_t={(V_t^X,V_t^Y,V_t^Z)}^T$, which are calculated from the target motion parameters ${\mathbf{x}}_t={(x_t,y_t,z_t)}^T$, $v_t$, and ${\phi }_t$ from SDS (1) by the following rules:

$$\left(\begin{array}{c}{\mathbf{X}}_t\\ {\mathbf{V}}_t\end{array}\right)=\left(\begin{array}{cc}\mathbf{A}& 0\\ 0& \mathbf{A}\end{array}\right)\left(\begin{array}{c}{\mathbf{x}}_t\\ {\mathbf{v}}_t\end{array}\right)+\left(\begin{array}{c}{\mathbf{R}}^0\\ 0\end{array}\right),$$

(2)

where ${\mathbf{v}}_t={(v_{t}cos{\phi }_t,v_{t}sin{\phi }_t,0)}^T$ and $\mathbf{A}$ is the rotation matrix

$$\mathbf{A}=\left(\begin{array}{ccc}cos\alpha & -sin\alpha cos\beta & sin\alpha sin\beta \\ sin\alpha & cos\alpha cos\beta & -cos\alpha sin\beta \\ 0& sin\beta & cos\beta \end{array}\right),$$

The parameters of the target coordinate system rotation $\alpha$, $\beta$, and shift ${\mathbf{R}}^0$ are assumed to be unknown random values with known distributions. These parameters are subject to estimation, and, to that end, we introduce the following auxiliary dynamic equations:

$$\left\{\begin{array}{c}d{\alpha }_t=0,\hfill \\ d{\beta }_t=0,\hfill \\ d{\mathbf{R}}_t^0={\left(d{R}_0^X\left(t\right),d{R}_0^Y\left(t\right),d{R}_0^Z\left(t\right)\right)}^T=0\hfill \end{array}\right.$$

(3)

with initial conditions ${\alpha }_0=\alpha$, ${\beta }_0=\beta$, ${\mathbf{R}}_0^0={\mathbf{R}}^0$.

Summing up the planar target dynamic Equation (1), the coordinate transform Equation (2), and the equations for the auxiliary variables (3), we have a 17-dimensional stochastic differential system.

2.2. Observations

The target is observed by a set of N registration devices (acoustic sensors) which measure the bearing and the Doppler frequency shift. The coordinates $\left\{{\mathbf{X}}^i\right\}=\left\{{(X^i,Y^i,Z^i)}^T\right\}$, $i=1,\dots ,N$ of these devices are assumed to be known, as well as the precision of the measurements they provide. The observations are available at known time instants $t_k$, $k=0,\dots ,K$ of the observation interval $t\in [0,T]$ with $t_0=0$ and $t_K=T$. For the sake of simplicity, we presume equal intervals between the observations: $\Delta =t_{k+1}-t_k$, $k=0,\dots ,K-1$. The observable bearing is connected with the target motion parameters by the following relations:

$$\left\{\begin{array}{c}{\displaystyle {\xi }_{t_k}^i=\frac{Z_{t_k}-Z^i}{R_{t_k}^i}+v_k^{{\xi }^i},}\hfill \\ {\displaystyle {\eta }_{t_k}^i=\frac{X_{t_k}-X^i}{r_{t_k}^i}+v_k^{{\eta }^i},}\hfill \end{array}\right.$$

(4)

where $R^i=\parallel {\mathbf{R}}^i\parallel$ is the range from the i-th observer to the target, and $r^i$ is the length of the position vector ${\mathbf{R}}^i$ projection on the plane $Z=0$:

$$\begin{array}{c}{R}_{t_k}^i=\sqrt{{(X_{t_k}-X^i)}^2+{(Y_{t_k}-Y^i)}^2+{(Z_{t_k}-Z^i)}^2},\hfill \\ r_{t_k}^i=\sqrt{{(X_{t_k}-X^i)}^2+{(Y_{t_k}-Y^i)}^2}.\hfill \end{array}$$

The observation system is presented in detail in Figure 2. As one can see, the observations ${\xi }^i=cos{\psi }^i$ and ${\eta }^i=cos{\theta }^i$ are the direction cosines of ${\mathbf{R}}^i$. Note that this structure of observations presumes the presence of the port-starboard ambiguity.

The Doppler measurements bring information concerning the radial target speed towards the sensor. The Doppler data processing makes the accuracy of the target state estimate significantly better. In the case of active tracking, the target is irradiated by an emitter with a known fixed frequency, and the receiver compares a frequency of the reflected signal with the original. In the case of passive observations, tracking becomes a challenging problem. Only the acoustic noise (so-called target acoustic signature) distorted by the target motion relative to the sensor is available for observation. The signature identification itself is an individual non-trivial problem [39,40,41,42,43,44] outside the scope of our study. To simplify the presentation, we suppose that the acoustic signature is known and determined by some known fixed frequency ${\omega }_0$, whose Doppler distortions we can register.

Thus, the Doppler shift measurements are given by the following expression:

$${\omega }_{t_k}^i=\frac{{\omega }_0}{1-V_{t_k}^i/C}+v_k^{{\omega }^i},$$

(5)

where ${\omega }_0$ is the frequency of sound emitted by the target, ${\omega }^i$ is the frequency registered by the observer i, C is the sound speed in the medium, and $V_{t_k}^i$ is the relative speed of the target with respect to the observer i:

$$V_{t_k}^i=\frac{V_{t_k}^X(X_{t_k}-X^i)+V_{t_k}^Y(Y_{t_k}-Y^i)+V_{t_k}^Z(Z_{t_k}-Z^i)}{R_{t_k}^i}.$$

The observation noise ${\left(v_k^{{\xi }^1},v_k^{{\eta }^1},v_k^{{\omega }^1},\dots ,v_k^{{\xi }^N},v_k^{{\eta }^N},v_k^{{\omega }^N}\right)}^T$ is a sequence of independent and identically distributed random vectors with uncorrelated components. The mean values and standard deviations of the noise components are assumed to be known.

2.3. Tracking Problem Formulation

Now that the target dynamics and the observation model are defined, the tracking problem can be formulated as filtering, i.e., search for the estimate ${\widehat{\mathcal{X}}}_t={\left({\widehat{\mathbf{X}}}_t,{\widehat{\mathbf{V}}}_t,{\widehat{\mathbf{x}}}_t,{\widehat{v}}_t,{\widehat{\varphi }}_t,{\widehat{a}}_t^n,{\widehat{\alpha }}_t,{\widehat{\beta }}_t,{\widehat{\mathbf{R}}}_t^0\right)}^T$ of the target motion parameters

$${\mathcal{X}}_t={\mathcal{X}}_0+{\int }_0^{t}a\left({\mathcal{X}}_s\right)ds+{\int }_0^{t}bd{\ell }_s,$$

(6)

where $a\left({\mathcal{X}}_s\right)$ and b are defined by SDS (1)–(3) given the bearing-only observations

$${\mathcal{Y}}_k^b=col\left(\left\{{({\xi }_{t_k}^i,{\eta }_{t_k}^i)}^T\right\}\right),$$

(7)

defined by Equation (4), or combined bearing-Doppler observations

$${\mathcal{Y}}_k^{bd}=col\left(\left\{{({\xi }_{t_k}^i,{\eta }_{t_k}^i,{\omega }_{t_k}^i)}^T\right\}\right),$$

(8)

defined by Equations (4) and (5), $i=1,\dots ,N$, $k=0,\dots ,K$. Here, $col(·)$ operator applied to a set of vectors $\left\{{\mathbf{a}}_i\right\}$, $i=1,\dots ,N$ results in vertically stacking of those vectors into one: $col\left(\left\{{\mathbf{a}}_i\right\}\right)={({\mathbf{a}}_1^T,\dots {\mathbf{a}}_N^T)}^T$.

3. Basics of Conditionally-Minimax Filtering

This section contains a brief description of the CM framework and the theorems necessary for the filter design [35,36].

3.1. Minimax Estimation in Regression

The CMNF algorithm is based on the solution to the following minimax estimation problem. Let us consider a compound random vector $Q\triangleq col(\mathcal{X},\mathcal{Y})\in {\mathbb{R}}^{n+m}$, where $\mathcal{X}\in {\mathbb{R}}^n$ is an unobservable estimated part, and $\mathcal{Y}\in {\mathbb{R}}^m$ is the observation. The probability distribution $\mathrm{P}$ of Q is uncertain: the mean $m_Q\triangleq {\mathbf{E}}_{\mathrm{P}}\left\{Q\right\}=col(m_{\mathcal{X}},m_{\mathcal{Y}})$ is known, but the covariance $R^Q\triangleq {cov}_{\mathrm{P}}(Q,Q)$ is unknown but bounded:

$$\underline{R}\triangleq \left(\begin{array}{cc}{\underline{R}}^{\mathcal{XX}}& {\underline{R}}^{\mathcal{XY}}\\ {\underline{R}}^{\mathcal{YX}}& {\underline{R}}^{\mathcal{YY}}\end{array}\right)⩽R^Q⩽\overline{R}\triangleq \left(\begin{array}{cc}{\overline{R}}^{\mathcal{XX}}& {\overline{R}}^{\mathcal{XY}}\\ {\overline{R}}^{\mathcal{YX}}& {\overline{R}}^{\mathcal{YY}}\end{array}\right).$$

Here, $\underline{R}, \overline{R}\in {\mathbb{R}}_{+}^{(n+m)\times (n+m)}$ are known positively semidefinite matrices and matrix inequality $A⩽B$ means that the matrix $B-A$ is positive semidefinite. The set of all feasible distributions $\mathrm{P}$ is denoted by $\mathcal{P}$. All functions $\psi : {\mathbb{R}}^m\to {\mathbb{R}}^n$ such that ${\displaystyle \underset{\mathrm{P}\in \mathcal{P}}{sup}{\mathbf{E}}_{\mathrm{P}}\left\{{\parallel \psi \left(\mathcal{Y}\right)-\mathcal{X}\parallel }^2\right\}<+\infty }$ form the class of admissible estimators $\Psi$. The performance of an estimator $\psi (·)$ given the distribution $\mathrm{P}$ is determined by the traditional MS-criterion $J(\mathrm{P},\psi )\triangleq {\mathbf{E}}_{\mathrm{P}}\left\{{\parallel \psi \left(\mathcal{Y}\right)-\mathcal{X}\parallel }^2\right\}$.

The minimax estimation problem of $\mathcal{X}$ given $\mathcal{Y}$ is to find $\widehat{\mathcal{X}}\triangleq \widehat{\psi }\left(\mathcal{Y}\right)$ such that

$$\widehat{\psi }(·)\in \underset{\psi \in \Psi }{\mathrm{Argmin}}\underset{\mathrm{P}\in \mathcal{P}}{sup}J(\mathrm{P},\psi ).$$

Theorem 1.

The minimax estimate $\widehat{\mathcal{X}}$ is linear

$$\widehat{\mathcal{X}}=\widehat{\psi }\left(\mathcal{X}\right)\triangleq m_{\mathcal{X}}+{\overline{R}}^{\mathcal{XY}}{\left({\overline{R}}^{\mathcal{YY}}\right)}^{+}(\mathcal{Y}-m_{\mathcal{Y}}),$$

(9)

where $A^{+}$ denotes matrix pseudo-inverse. On the set $\mathcal{P}\times \Psi$, the criterion J has saddle points, i.e.,

$$\underset{\psi \in \Psi }{min}\underset{\mathrm{P}\in \mathcal{P}}{max}J(\mathrm{P},\psi )=\underset{\mathrm{P}\in \mathcal{P}}{max}\underset{\psi \in \Psi }{min}J(\mathrm{P},\psi ).$$

The Gaussian distribution $\mathcal{N}(m_Q,\overline{R})$ with the expectation $m_Q$ and the maximal feasible covariance matrix $\overline{R}$ is the least favorable one, and the pair $(\mathcal{N}(m_Q,\overline{R}),\widehat{\psi })$ is one of the saddle points. Hence, the double inequality

$$J(\mathrm{P},\widehat{\psi })⩽J(\mathcal{N}(m_Q,\overline{R}),\widehat{\psi })⩽J(\mathcal{N}(m_Q,\overline{R}),\psi )$$

holds for all $(\mathrm{P},\psi )\in \mathcal{P}\times \Psi$.

The estimate $\widehat{\mathcal{X}}$ has guaranteeing accuracy:

$$\mathbf{E}\left\{(\widehat{\mathcal{X}}-\mathcal{X}){(\widehat{\mathcal{X}}-\mathcal{X})}^{\top }\right\}⩽{\overline{R}}^{\mathcal{XX}}-{\overline{R}}^{\mathcal{XY}}{\left({\overline{R}}^{\mathcal{YY}}\right)}^{+}{\overline{R}}^{\mathcal{YX}},$$

and the equality holds for any $\mathrm{P}\in \mathcal{P}$ with ${cov}_{\mathrm{P}}(Q,Q)=\overline{R}$.

The proof of Theorem 1 is given in [35,36]. The theorem plays a key role in the development of the CMNF algorithm. First, it guarantees the estimate unbiasedness and gives a possibility to calculate a conservative approximation of the accuracy characteristic, adequate to reality. Second, the minimax estimate can be obtained under a priori uncertainty of probability distribution. To calculate the minimax estimate, it is sufficient to know the mean and some upper bound of the covariance. The closer the obtained upper bound is to the real covariance, the closer the conservative accuracy approximation is to reality. Third, the extension of the admissible estimators by the nonlinear functions given initial observations only does not enhance the estimation quality. Fourth, since the least favorable distribution is Gaussian, the moments of the 1st and 2nd order form a closed set of sufficient statistics in the problem, and, at the same time, it gives a certain justification of usage of the Gaussianity assumption concerning both the prediction and correction terms in the EKF and UKF. In fact, estimate (9) is the best linear one for the least favorable distribution $\mathcal{N}(m_Q,\overline{R})$. Note that the regression estimates calculated either via statistical linearization method or UT-transform, being affine functions of the observations, can be considered as some approximations of the minimax estimate $\widehat{\mathcal{X}}$: terms $m_{\mathcal{X}}$ and $m_{\mathcal{Y}}$ along with the gain matrix ${\overline{R}}^{\mathcal{XY}}{({\overline{R}}^{\mathcal{YY}})}^{+}$ are replaced in Formula (9) by some estimates of fixed precision. The term “fixed precision” means that there is no way of asymptotic approximation of $m_{\mathcal{X}}$, $m_{\mathcal{Y}}$ and ${\overline{R}}^{\mathcal{XY}}{({\overline{R}}^{\mathcal{YY}})}^{+}$ neither by the statistical linearization nor by the UT-transform. By contrast, the CMNF algorithm uses the Monte Carlo method and gives a possibility to calculate coefficients with arbitrary precision.

The key advantage of the CM framework is the possibility to construct the estimate, which is better than the optimal linear estimate (9), amplifying the initial observation $\mathcal{Y}$ with some transform. Let us introduce a function $z=z\left(y\right): {\mathbb{R}}^m\to {\mathbb{R}}^k$, such that ${sup}_{\mathrm{P}\in \mathcal{P}}{\mathbf{E}}_{\mathrm{P}}\left\{{\parallel z\left(\mathcal{Y}\right)\parallel }^2\right\}<\infty$ and the corresponding compound vector $Q^{\prime }\triangleq col(\mathcal{X},\mathcal{Y},z\left(\mathcal{Y}\right))$. Again, the probability distribution ${\mathrm{P}}^{\prime }$ of $Q^{\prime }$ is supposed to be uncertain: the mean $m_{Q^{\prime }}\triangleq {\mathbf{E}}_{{\mathrm{P}}^{\prime }}\left\{Q^{\prime }\right\}=col(m_{\mathcal{X}},m_{\mathcal{Y}},m_Z)$ is known, the covariance ${\mathbf{R}}^{Q^{\prime }}\triangleq {cov}_{\mathcal{P}}(Q^{\prime },Q^{\prime })$ is unknown but bounded:

$$\underline{\mathbf{R}}\triangleq \left(\begin{array}{ccc}{\underline{R}}^{\mathcal{XX}}& {\underline{R}}^{\mathcal{XY}}& {\underline{\mathbf{R}}}^{\mathcal{XZ}}\\ {\underline{R}}^{\mathcal{YX}}& {\underline{R}}^{\mathcal{YY}}& {\underline{\mathbf{R}}}^{\mathcal{YZ}}\\ {\underline{\mathbf{R}}}^{\mathcal{ZX}}& {\underline{\mathbf{R}}}^{\mathcal{ZY}}& {\underline{\mathbf{R}}}^{\mathcal{ZZ}}\end{array}\right)⩽{\mathbf{R}}^{Q^{\prime }}⩽\overline{\mathbf{R}}\triangleq \left(\begin{array}{ccc}{\overline{R}}^{\mathcal{XX}}& {\overline{R}}^{\mathcal{XY}}& {\overline{\mathbf{R}}}^{\mathcal{XZ}}\\ {\overline{R}}^{\mathcal{YX}}& {\overline{R}}^{\mathcal{YY}}& {\overline{\mathbf{R}}}^{\mathcal{YZ}}\\ {\overline{\mathbf{R}}}^{\mathcal{ZX}}& {\overline{\mathbf{R}}}^{\mathcal{ZY}}& {\overline{\mathbf{R}}}^{\mathcal{ZZ}}\end{array}\right).$$

Here, $\underline{\mathbf{R}}, \overline{\mathbf{R}}\in {\mathbb{R}}_{+}^{(n+m+k)\times (n+m+k)}$ are known positively semidefinite matrices, including the lower and upper bounds $\underline{R}$ and $\overline{R}$ as the corresponding block components. The set of all feasible distributions ${\mathrm{P}}^{\prime }$ is denoted by ${\mathcal{P}}^{\prime }$. All functions ${\psi }^{\prime }: {\mathbb{R}}^{m+k}\to {\mathbb{R}}^n$ such that ${sup}_{{\mathrm{P}}^{\prime }\in {\mathcal{P}}^{\prime }}{\mathbf{E}}_{{\mathrm{P}}^{\prime }}\left\{{\parallel \psi (\mathcal{Y},z\left(\mathcal{Y}\right))-\mathcal{X}\parallel }^2\right\}<+\infty$ form the class of admissible estimators ${\Psi }^{\prime }$. The performance of estimator ${\psi }^{\prime }(·)$ given the distribution ${\mathrm{P}}^{\prime }$ is determined as $J({\mathrm{P}}^{\prime },{\psi }^{\prime })\triangleq {\mathbf{E}}_{{\mathrm{P}}^{\prime }}\left\{{\parallel \psi (\mathcal{Y},z\left(\mathcal{Y}\right))-\mathcal{X}\parallel }^2\right\}$. In view of Theorem 1, the minimax estimate ${\widehat{\mathcal{X}}}^{\prime }$ given $col(\mathcal{Y},z(\mathcal{Y}\left)\right)$ takes the form

$${\widehat{\mathcal{X}}}^{\prime }=m_{\mathcal{X}}+\left(\begin{array}{cc}{\overline{R}}^{\mathcal{XY}}& {\overline{\mathbf{R}}}^{\mathcal{XZ}}\end{array}\right){\left(\begin{array}{cc}{\overline{R}}^{\mathcal{YY}}& {\overline{\mathbf{R}}}^{\mathcal{YZ}}\\ {\overline{\mathbf{R}}}^{\mathcal{ZY}}& {\overline{\mathbf{R}}}^{\mathcal{ZZ}}\end{array}\right)}^{+}\left(\left(\begin{array}{c}\mathcal{Y}\\ z\left(\mathcal{Y}\right)\end{array}\right)-\left(\begin{array}{c}{m}_{\mathcal{Y}}\\ m_{\mathcal{Z}}\end{array}\right)\right),$$

and its guaranteeing accuracy is defined by the inequality:

$$\mathbf{E}\left\{({\widehat{\mathcal{X}}}^{\prime }-\mathcal{X}){({\widehat{\mathcal{X}}}^{\prime }-\mathcal{X})}^{\top }\right\}⩽{\overline{R}}^{\mathcal{XX}}-\left(\begin{array}{cc}{\overline{R}}^{\mathcal{XY}}& {\overline{\mathbf{R}}}^{\mathcal{XZ}}\end{array}\right){\left(\begin{array}{cc}{\overline{R}}^{\mathcal{YY}}& {\overline{\mathbf{R}}}^{\mathcal{YZ}}\\ {\overline{\mathbf{R}}}^{\mathcal{ZY}}& {\overline{\mathbf{R}}}^{\mathcal{ZZ}}\end{array}\right)}^{+}\left(\begin{array}{c}{\overline{R}}^{\mathcal{YX}}\\ {\overline{\mathbf{R}}}^{\mathcal{ZX}}\end{array}\right)\triangleq \widehat{\mathbf{R}}.$$

From the properties of nonnegatively semidefinite matrix pseudo-inverse, it follows that $\widehat{\mathbf{R}}⩽\widehat{R}$ for any admissible $z\left(y\right)$, i.e., taking into account any admissible observation transform enhances the estimation performance. The effective choice of an appropriate transform looks like some sort of art, and it should be done reasoning from the properties of a specific observation system at hand.

3.2. CMNF in Continuous-Discrete Systems

On the probability triplet with filtration $(\Omega ,\mathcal{F},\mathrm{P},\left\{{\mathcal{F}}_t\right\})$, we consider the observation system

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{X}}_t={\mathcal{X}}_0+{\int }_0^{t}a(s,{\mathcal{X}}_s)ds+{\int }_0^{t}b(s,{\mathcal{X}}_{s-})d{\ell }_s,\hspace{1em}t>0,}\hfill \\ {\displaystyle {\mathcal{Y}}_k=c(t_k,{\mathcal{X}}_{t_k},v_k),\hspace{1em}{t}_{k+1}-t_k=\Delta , k\in \mathbb{N}.}\hfill \end{array}\right.$$

Here,

–

${\mathcal{X}}_t\in {\mathbb{R}}^n$ is an unobservable state process subject to estimation with an initial condition ${\mathcal{X}}_0$,

–

${\ell }_t\in {\mathbb{R}}^n$ is an ${\mathcal{F}}_t$-adapted square-integrable Levy process,

–

${\mathcal{Y}}_k\in {\mathbb{R}}^m$ is the observation sequence, available at the equidistant instants $t_k$ with the step $\Delta$,

–

$v_k\in {\mathbb{R}}^m$ is a ${\mathcal{F}}_{t_k}$-adapted sequence of i.i.d. random vectors; $x_0$, $\left\{{\ell }_t\right\}$ and $\left\{v_k\right\}$ are supposed to be mutually independent.

We assume for the drift $a(·,·)$ and diffusion $b(·,·)$ the standard conditions guaranteeing the existence and uniqueness of a strong solution to the stochastic differential equation for ${\mathcal{X}}_t$ along with finiteness of its moments [45]. The observation function $c(·,·,·)$ is also assumed to provide the existence of the moments up to order 2 for the observations ${\mathcal{Y}}_k$.

Generally, the problem is to find an estimate ${\tilde{\mathcal{X}}}_k$ of the state ${\mathcal{X}}_k={\mathcal{X}}_{t_k}$, given the available observations up to the moment $t_k$. The estimation quality is characterized by a traditional MS-criterion.

The filter structure is partially determined by three sequences of functions $\left\{{\alpha }_k\right\}$, $\left\{{\beta }_k\right\}$ and $\left\{{\gamma }_k\right\}$:

–

${\alpha }_k\left(x\right):{\mathbb{R}}^n\to {\mathbb{R}}^n$ is basic prediction function,

–

${\beta }_k\left(y\right):{\mathbb{R}}^m\to {\mathbb{R}}^k$ is observation transform,

–

${\gamma }_k(x,z):{\mathbb{R}}^n\times {\mathbb{R}}^k\to {\mathbb{R}}^k$ is basic correction function.

The admissible filter has the following recursive form

$${\tilde{\mathcal{X}}}_k={\varphi }_k\left({\alpha }_k\left({\tilde{\mathcal{X}}}_{k-1}\right),{\gamma }_k({\tilde{\mathcal{X}}}_{k-1},{\beta }_k\left({\mathcal{Y}}_k)\right)\right),$$

(10)

where the sequence $\left\{{\varphi }_k\right\}$ of functions ${\varphi }_k(x,z):{\mathbb{R}}^n\times {\mathbb{R}}^k\to {\mathbb{R}}^n$ should be optimized. It is supposed that the functions ${\alpha }_k(·)$, ${\beta }_k(·)$ and ${\gamma }_k(·)$ of the corresponding random arguments in Formula (10) have finite moments up to order 2. The set ${\Phi }_k$ of admissible filters is also formed by all the measurable functions, such that $\mathbf{E}\left\{\parallel {\tilde{\mathcal{X}}}_k{\parallel }^2\right\}<\infty$.

The conditions above guarantee that the compound random vector

$$Q_k\triangleq col({\mathcal{X}}_k,{\alpha }_k\left({\tilde{\mathcal{X}}}_{k-1}\right),{\gamma }_k({\tilde{\mathcal{X}}}_{k-1},{\beta }_k\left({\mathcal{Y}}_k\right)))$$

has two moments $m_k^Q\triangleq \mathbf{E}\left\{Q_k\right\}$ and $R_k^Q\triangleq cov(Q_k,Q_k)$:

$$m_k^Q=\left(\begin{array}{c}{m}_k^x\\ m_k^{\alpha }\\ m_k^{\gamma }\end{array}\right),\hspace{1em}{R}_k^Q=\left(\begin{array}{ccc}{R}_k^{xx}& R_k^{x\alpha }& R_k^{x\gamma }\\ R_k^{\alpha x}& R_k^{\alpha \alpha }& R_k^{\alpha \gamma }\\ R_k^{\gamma x}& R_k^{\gamma \alpha }& R_k^{\gamma \gamma }\end{array}\right).$$

Let us denote by ${\mathcal{P}}_k$ the set of all probability distributions of $(2n+k)$-dimensional random vectors with known expectation and unknown covariance matrix bounded from above: $R_k^Q⩽{\overline{R}}_k$ with some known positive semidefinite matrix ${\overline{R}}_k$.

The conditionally-minimax filtering problem is to find

$${\widehat{\varphi }}_k\in \underset{{\varphi }_k\in {\Phi }_k}{\mathrm{Argmin}} \underset{\mathrm{P}\in {\mathcal{P}}_k}{sup} {\mathbf{E}}_{\mathrm{P}}\left\{\parallel {\varphi }_k\left({\alpha }_k\left({\tilde{\mathcal{X}}}_{k-1}\right),{\gamma }_k({\tilde{\mathcal{X}}}_{k-1},{\beta }_k\left({\mathcal{Y}}_k)\right)\right)-{\mathcal{X}}_k{\parallel }^2\right\}.$$

Theorem 2.

The CMNF has a two-step “prediction–correction” form. The prediction and the upper bound of its error covariance matrix are defined as follows:

$$\begin{array}{c}{\breve{\mathcal{X}}}_k={\overline{R}}_k^{x\alpha }{({\overline{R}}_k^{\alpha \alpha })}^{+}{\alpha }_k({\widehat{\mathcal{X}}}_{k-1})+[m_k^x-{\overline{R}}_k^{x\alpha }{({\overline{R}}_k^{\alpha \alpha })}^{+}{m}_k^{\alpha }],\hfill \\ cov({\breve{\mathcal{X}}}_k-{\mathcal{X}}_k,{\breve{\mathcal{X}}}_k-{\mathcal{X}}_k)⩽{\breve{K}}_k^{xx}\triangleq {\overline{R}}_k^{xx}-{\overline{R}}_k^{x\alpha }{({\overline{R}}_k^{\alpha \alpha })}^{+}{\overline{R}}_k^{\alpha x}.\hfill \end{array}$$

(11)

The correction and the upper bound of its error covariance matrix are defined as follows:

$$\begin{array}{cc}{\widehat{\mathcal{X}}}_k={\breve{\mathcal{X}}}_k+{\breve{K}}_k^{x\gamma }{[{\breve{K}}_k^{\gamma \gamma }]}^{+}{\breve{\gamma }}_k,\hspace{1em}\hfill & cov({\widehat{\mathcal{X}}}_k-{\mathcal{X}}_k,{\widehat{\mathcal{X}}}_k-{\mathcal{X}}_k)⩽{\widehat{K}}_k^{xx}\triangleq {\breve{K}}_k^{xx}-{\breve{K}}_k^{x\gamma }{({\breve{K}}_k^{\gamma \gamma })}^{+}{\breve{K}}_k^{\gamma x},\hfill \end{array}$$

(12)

where

$${\breve{\gamma }}_k={\gamma }_k({\widehat{\mathcal{X}}}_{k-1},{\beta }_k({\mathcal{Y}}_k))-m_k^{\gamma }-{\overline{R}}_k^{\gamma \alpha }{({\overline{R}}_k^{\alpha \alpha })}^{+}({\alpha }_k({\widehat{\mathcal{X}}}_{k-1})-m_k^{\alpha })$$

is the modified correction, and

$$\begin{array}{cc}{\breve{K}}_k^{x\gamma }\triangleq {\overline{R}}_k^{x\gamma }-{\overline{R}}_k^{x\alpha }{({\overline{R}}_k^{\alpha \alpha })}^{+}{\overline{R}}_k^{\alpha \gamma },\hspace{1em}\hfill & {\breve{K}}_k^{\gamma \gamma }\triangleq {\overline{R}}_k^{\gamma \gamma }-{\overline{R}}_k^{\gamma \alpha }{({\overline{R}}_k^{\alpha \alpha })}^{+}{\overline{R}}_k^{\alpha \gamma }.\hfill \end{array}$$

Proof of Theorem 2 follows immediately from assertions of Theorem 1 and the recursive variant of Theorem 13.1 from [46]. The proposed CMNF has the following positive properties similar to the minimax regression estimates: unbiasedness and conservative approximation of the covariance matrix. Moreover, the CMNF algorithm gives additional possibilities. First, the transformation of observations gives a chance to raise the estimation accuracy. For example, the inclusion of the powers of the original observations into the calculation process allows for computing polynomial estimates of better accuracy than provided by the best linear estimates [47]. Second, the proper choice of the basic prediction and correction functions for each specified observation system also allows for raising the estimation quality to the optimal. For instance, in the case of the linear Gaussian observation system, identical observation transform, linear prediction, and correction, the proposed CMNF turns into the classical KF [48].

There is a considerable issue of CMNF implementation. To calculate CMNF, one needs to know the expectation and covariance matrix (or at least its conservative approximation) of the compound vectors $\left\{Q_k\right\}$. The analytical calculation of these characteristics is possible only for a narrow class of observation systems, including linear Gaussian and the ones with Markov jump processes as the states [47]. In general, the calculation of these parameters is possible through numerically expensive Monte Carlo sampling. The approximations of the mean can be obtained by regular sampling; meanwhile, the covariance matrix should be approximated from above, using the asymptotic confidence regions. The desired sample size can be controlled based on the results of [49].

4. Comparative Numerical Study

4.1. Estimation of Static Cartesian Coordinates Given Bearing-Only Multisensor Observations

In this section, we present the simulation results for the static case, which can be considered as the initial step of the filtering problem for the dynamic system presented in Section 2. All the simulation in Section 4.1 and Section 4.2 was done using Python, the source code is available in public repository [50].

The position of the target is a random vector ${\mathbf{X}}_0$ with independent components. The distribution of each component is Gaussian with the following parameters (mean and standard deviation in meters):

$$X_0\sim \mathcal{N}(0,{10}^3),\hspace{1em}{Y}_0\sim \mathcal{N}(y_0,{10}^3),\hspace{1em}{Z}_0\sim \mathcal{N}(-{10}^3,{10}^2),$$

(13)

where the parameter $y_0$ takes values $\{5,10,20\}\times {10}^3$. This diversity is introduced in order to demonstrate how the estimates depend on the degree of the nonlinearity of the system at hand.

Four pairs of registration devices (acoustic sonars) are located at the points $\left\{{\mathbf{X}}^i\right\}$, $i=1,\dots ,8$:

$$\begin{array}{cccc}{(-10000,0,-25)}^T,\hfill & {(-5000,1000,-25)}^T,\hfill & {(5000,1000,-25)}^T,\hfill & {(10000,0,-25)}^T,\hfill \\ {(-10000,0,-50)}^T,\hfill & {(-5000,1000,-50)}^T,\hfill & {(5000,1000,-50)}^T,\hfill & {(10000,0,-50)}^T.\hfill \end{array}$$

(14)

The observation noise components $(v_k^{{\xi }^i},v_t^{{\eta }^i})$ are assumed to have Gaussian distribution $\mathcal{N}(0,0.02)$, which roughly approximates the sonar precision equal to $1^{°}$.

Note that, due to the port-starboard ambiguity, the bearing observations (4) do not allow an explicit inverse transformation to calculate the target coordinates, even in the absence of the noise $(v_k^{{\xi }^i},v_t^{{\eta }^i})$. Nevertheless, the CMNF algorithm allows estimating $\left\{{\mathbf{X}}_k\right\}$ with relatively high precision. First, the algorithm relates to the Bayesian estimator class, which allows for neutralizing the ambiguity not having an individual observer but the site as in the presented example. Second, the CMNF allows for fusing optimally the original observations and some of their transforms, such as least squares (LS) estimates or other regression technique results. An alternative way is not to rely entirely on the current measurements of the target position, but also to take into account some pre-collected data of such observations. This can be done by the aids of one of the machine learning (ML) regression algorithms, e.g., lasso regression [51]. This can provide not only the estimation quality gain, but also significantly reduce the amount of time necessary for the online measurements processing. While the LS method requires a nonlinear optimization problem solution for each observation round, an ML regression algorithm can provide much faster predictions with all the resource and time demanding procedures of the model fit being performed a priori. For instance, the lasso regression coefficients calculation also requires a nonlinear optimization problem solution and the whole training set processing, but, after that, the target position estimate is easily computed as a polynomial of the measurements.

In

Table 1. Coordinate-wise std of estimation error ${\sigma }^E$ for $E\in \left\{\mathrm{Linear},\mathrm{LS},\mathrm{ML},\mathrm{CM} \mathrm{LS},\mathrm{CM} \mathrm{ML}\right\}$.

	${\mathit{y}}_0=5\times {10}^3$			${\mathit{y}}_0=10\times {10}^3$			${\mathit{y}}_0=20\times {10}^3$
	${\mathit{\sigma }}_{\mathit{X}}$	${\mathit{\sigma }}_{\mathit{Y}}$	${\mathit{\sigma }}_{\mathit{Z}}$	${\mathit{\sigma }}_{\mathit{X}}$	${\mathit{\sigma }}_{\mathit{Y}}$	${\mathit{\sigma }}_{\mathit{Z}}$	${\mathit{\sigma }}_{\mathit{X}}$	${\mathit{\sigma }}_{\mathit{Y}}$	${\mathit{\sigma }}_{\mathit{Z}}$
Linear	174.75	195.99	49.57	107.18	221.11	60.02	148.89	381.16	79.23
LS	139.57	116.69	51.31	105.62	165.49	74.61	146.97	403.34	130.72
ML	146.17	125.50	46.91	106.83	166.12	60.03	147.50	373.37	79.60
CM LS	137.28	114.79	45.56	105.04	162.44	59.69	145.30	371.56	79.22
CM ML	144.33	123.98	46.71	105.98	165.50	59.75	146.53	372.62	79.22

and Figure 3, we present the standard deviations of the estimation errors for the following estimators:

• Linear: the basic minimax estimator (9), which uses only the original bearing measurements ${\mathcal{Y}}^b=col\left(\left\{({\xi }^i,{\eta }^i)\right\}\right)$, defined by relations (7). It is also the best linear estimate.
• LS: Linear squares estimate ${\widehat{\mathcal{X}}}^{LS}={\widehat{\mathcal{X}}}^{LS}\left({\mathcal{Y}}^b\right)$, calculated using the eight pairs of bearing measurements.
• ML: Lasso regression estimate ${\widehat{\mathcal{X}}}^{ML}={\widehat{\mathcal{X}}}^{ML}\left({\mathcal{Y}}^b\right)$, which uses normalized bearing measurements and their squares as features.
• CM LS: Minimax estimator (9) with original bearing measurements supplemented by the LS estimate: ${\mathcal{Y}}^{LS}=col(\left\{({\xi }^i,{\eta }^i)\right\}\cup {\widehat{\mathcal{X}}}^{LS}\left({\mathcal{Y}}^b\right))$.
• CM ML: Minimax estimator (9) with original bearing measurements supplemented by the ML estimate: ${\mathcal{Y}}^{ML}=col(\left\{({\xi }^i,{\eta }^i)\right\}\cup {\widehat{\mathcal{X}}}^{ML}\left({\mathcal{Y}}^b\right))$.

The covariances $R^{\mathcal{X}\mathcal{Y}}$, which are necessary for the minimax estimates calculation with Formula (9), for the Linear, CM LS, and CM ML estimators, were calculated on the test set of ${10}^6$ position-measurement pairs by Monte Carlo sampling. A part of this test set was also used to fit the lasso regression model for estimators ML and CM ML. The standard deviations of the estimation errors were estimated on a different test set of the same size of ${10}^6$ samples.

From the results presented in Table 1 and Figure 3, one can conclude the following:

• The CM LS estimate uniformly outperforms all others (Table 1, in bold).
• All the presented algorithms of CM-family are able to neutralize the port-starboard ambiguity.
• The error standard deviations computed for the Linear and CM LS/ML estimators confirm the conclusion of Section 3.1. Indeed, the introduction of the additional measurements’ transform increased the performance of the minimax estimate in comparison with the best linear estimate. The last is true even for the case of $y_0=20\times {10}^3$, where the $LS$ estimate demonstrated very poor quality.
• The difference in the results for various distances between the target and the observers shows that the most benefit the minimax estimate can provide when the model at hand is far from linear. Indeed, for $y_0=20\times {10}^3$, the difference between the estimate error standard deviations for linear and nonlinear minimax estimates is negligible, which is not the case for $y_0=5\times {10}^3$, where the nonlinearity of the system is more pronounced.

4.2. Tracking of Maneuvering Target Given the Bearing-Doppler Multisensor Observations

In this section, we present the results of numerical simulation of a maneuvering target motion model (1)–(3) given the bearing-only (7) and bearing-Doppler observations (8).

The planar target motion in the target coordinate system $0xyz$ is defined by the SDS (1) with the acceleration equation parameters $\lambda =0.01$, $\nu =0.0$ and $\mu =0.01$. The initial conditions are:

$$(x_0,y_0,z_0)=(0,0,0),\hspace{1em}{v}_0\sim \mathcal{R}[5,12],\hspace{1em}{\varphi }_0\sim \mathcal{N}(-\pi /2,0.1),\hspace{1em}{a}_0^n\sim \mathcal{R}[-0.2,0.2].$$

Here, the parameters to the Gaussian distribution $\mathcal{N}(m,\sigma )$ are the mean and standard deviation in meters, and $\mathcal{R}[a,b]$ is the uniform distribution in the interval $[a,b]$. The random rotation and shift parameters for the target coordinate system are as follows:

$$\alpha , \beta \sim \mathcal{N}(0,\pi /36),\hspace{1em}{R}_0^X\sim \mathcal{N}(0,{10}^3),\hspace{1em}{R}_0^Y\sim \mathcal{N}(20\times {10}^3,{10}^3),\hspace{1em}{R}_0^Z\sim \mathcal{N}(-{10}^3,{10}^2).$$

It should be noted that the introduced initial conditions for the coordinates in the target coordinate system ${\mathbf{x}}_0={(x_0,y_0,z_0)}^T$ and the shift ${\mathbf{R}}^0={(R_0^X,R_0^Y,R_0^Z)}^T$ agree with the static model (13) for $y_0=20\times {10}^3$.

In Figure 4, we demonstrate 100 sample paths of the system (1)–(3) with the parameters defined above. On the left, for clarity, we present the paths lying on the maneuver plane $0xy$ and starting from the common initial position. On the right, the same paths, being distributed according to the random initial position, are pictured in the observer coordinate system $0XYZ$. Even though the model has only one scalar disturbing process $W_t$ in the equation for the normal acceleration $a_t^n$ in SDS (1), one can notice a rather wide variety of trajectories produced by the SDS with chosen parameters. At the same time, all the presented paths demonstrate smoothness, which is characteristic of the underwater environment.

The registration devices are also located at the same points (14) as it was assumed for the static case and have the same parameters of the noises in the bearing observations. For the Doppler measurements (5), the sound speed is set to $C=1500$ m/s, and the source frequency ${\omega }_0=20$ Hz. The noise in frequency measurements is also assumed to be Gaussian: $v_k^{\omega }\sim \mathcal{N}(0,0.005)$. Notably, the parameters of both the target motion and sonar accuracy are chosen similarly to those of the classic paper [23].

In this simulation study, we have compared the standard deviation of the estimate error for the following filters:

• LS: Linear squares estimate ${\widehat{\mathcal{X}}}_k^{LS}={\widehat{\mathcal{X}}}_k^{LS}\left({\mathcal{Y}}_k^b\right)$, calculated using the eight pairs of bearing measurements ${\mathcal{Y}}_k^b=col\left(\left\{({\xi }_k^i,{\eta }_k^i)\right\}\right)$. The same as in the static case, but applied along the whole path.
• ML: Lasso regression estimate ${\widehat{\mathcal{X}}}_k^{ML}={\widehat{\mathcal{X}}}_k^{ML}\left({\mathcal{Y}}_k^b\right)$, which uses normalized bearing measurements and their squares as features. The same as in the static case, but applied along the whole path.
• Basic CMNF: the conditionally minimax nonlinear filter (11) and (12), which uses the combined bearing and Doppler measurements ${\mathcal{Y}}_k^{bd}=col(\left\{({\xi }_k^i,{\eta }_k^i,{\omega }_k^i)\right\}$, defined by (8).
• CMNF LS: the conditionally minimax nonlinear filter (11) and (12), which uses the original bearing measurements and the LS estimate ${\mathcal{Y}}_k^{LS}=col(\left\{({\xi }_k^i,{\eta }_k^i)\right\}\cup {\widehat{\mathcal{X}}}_k^{LS}\left({\mathcal{Y}}_k^b\right))$.
• CMNF LS with Doppler: the conditionally minimax nonlinear filter (11) and (12), which uses the combined measurements and the LS estimate ${\mathcal{Y}}_k^{LS}=col(\left\{({\xi }_k^i,{\eta }_k^i,{\omega }_k^i)\right\}\cup {\widehat{\mathcal{X}}}_k^{LS}\left({\mathcal{Y}}_k^b\right))$.
• CMNF ML: the conditionally minimax nonlinear filter (11) and (12), which uses the original bearing measurements and the lasso estimate ${\mathcal{Y}}_k^{ML}=col(\left\{({\xi }_k^i,{\eta }_k^i)\right\}\cup {\widehat{\mathcal{X}}}_k^{ML}\left({\mathcal{Y}}_k^b\right))$.
• CMNF ML with Doppler: the conditionally minimax nonlinear filter (11) and (12), which uses the combined measurements and the ML estimate ${\mathcal{Y}}_k^{ML}=col(\left\{({\xi }_k^i,{\eta }_k^i,{\omega }_k^i)\right\}\cup {\widehat{\mathcal{X}}}_k^{ML}\left({\mathcal{Y}}_k^b\right))$.
• A priori estimate: equal to the mean ${\mathcal{X}}_t$. The standard deviation of this “trivial” estimate can be used as a measure of other filters’ performance: if a filter has std larger then that of a “trivial” estimate, then one can conclude that this filter diverges; if a filter’s std is much less than it deserves further analysis and application.

The model is simulated on the time interval $t\in [0,100]$ s with discretization step $\Delta =1$ s. The size of the training set for the CMNF filters, as well as of the test set, is equal to ${10}^5$ sample paths.

It should be noted that the presented model turns out to be quite challenging for EKF and UKF filters, which are common to apply when dealing with such nonlinear systems. Even in milder conditions, when only the planar motion is considered, and, along with the bearing observations, the distance is available; both mentioned filters demonstrate divergence [37].

Table 2. Parameter-wise standard deviation of estimation error ${\sigma }^F\left(T\right)$ for $F\in$ {LS, ML, Basic CMNF (Doppler), CMNF LS (Doppler), CMNF ML (Doppler), a priori}.

	${\mathit{\sigma }}_{\mathit{X}}\left(\mathit{T}\right)$	${\mathit{\sigma }}_{\mathit{Y}}\left(\mathit{T}\right)$	${\mathit{\sigma }}_{\mathit{Z}}\left(\mathit{T}\right)$	${\mathit{\sigma }}_{{\mathit{V}}^{\mathit{X}}}\left(\mathit{T}\right)$	${\mathit{\sigma }}_{{\mathit{V}}^{\mathit{Y}}}\left(\mathit{T}\right)$	${\mathit{\sigma }}_{{\mathit{V}}^{\mathit{Z}}}\left(\mathit{T}\right)$	${\mathit{\sigma }}_{\mathit{v}}\left(\mathit{T}\right)$	${\mathit{\sigma }}_{\mathit{a}}\left(\mathit{T}\right)$	${\mathit{\sigma }}_{\mathit{\phi }}\left(\mathit{T}\right)$
LS	143.36	382.62	127.03	—	—	—	—	—	—
ML	148.62	362.73	89.55	—	—	—	—	—	—
Basic CMNF	19.30	39.37	40.84	0.345	0.152	0.404	0.094	0.053	0.278
CMNF LS	36.67	77.61	43.15	1.701	2.054	0.424	1.063	0.061	0.467
CMNF LS Doppler	19.30	39.37	40.83	0.345	0.152	0.404	0.094	0.053	0.278
CMNF ML	36.66	77.60	43.11	1.702	2.052	0.423	1.062	0.061	0.464
CMNF ML Doppler	19.30	39.36	40.82	0.345	0.152	0.404	0.094	0.053	0.278
a priori	—	—	—	—	—	—	2.0	0.115	0.1

presents the standard deviations of the estimation errors for all the mentioned filters on the last step, which corresponds to the end of the observation interval $T=100$ s.

In Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9, the standard deviations for the parameter estimation errors are presented with respect to the time variable $t\in [0,T]$, $T=100$ s. Table 1 shows little or no difference between the results for Basic CMNF, CMNF LS Doppler, CMNF ML Doppler filters, and between the results for CMNF LS and CMNF ML filters on the last step. The analysis of the collected data shows that all CMNF filters with Doppler measurements demonstrate almost the same behavior at all other time instants as well. Thus, in the following figures, the Basic CMNF, CMNF LS Doppler, and CMNF ML Doppler estimation error standard deviations are shown by the same line ${\sigma }_{Doppler}$. The same is true for the CMNF filters based on the bearing observations only, so the CMNF LS and CMNF ML estimation error standard deviations are also depicted by the same line ${\sigma }_{NoDoppler}$.

Figure 5 shows the superiority of the CMNF filters for the target coordinate estimation in comparison with all the static estimates presented in the Section 4.1. Together with Figure 6, it also provides evidence of a great performance boost provided by the Doppler measurements.

Figure 7 and Figure 8 demonstrate that, along with the target motion parameters, the parameters of its maneuver, namely the heading speed $v_t$ and normal acceleration $a_t$, can be estimated with a quality much better than one of the a priori estimate. The same is true for the target coordinate system shift ${\mathbf{R}}^0$ since the a priori estimates standard deviations are equal to ${\sigma }_{prior}={(1000,1000,100)}^T$.

In contrast, Figure 9 shows rather poor performance for all the filters for the estimation of the angles $\alpha$ and $\beta$: the quality is almost the same as the quality of the estimate given by the mean value $\widehat{\alpha }=\widehat{\beta }=0$. At the same time, the standard deviation of the heading angle error ${\phi }_t$ is notably less than that of the “trivial” estimate ${\widehat{\phi }}_t=\mathbf{E}\left\{{\phi }_t\right\}$.

The results presented in Table 2 and Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9 show the following:

• The conditionally minimax nonlinear filter provides much better estimates for the target coordinates than the static ones, which do not take into account the information provided by the prediction “by virtue” of the system.
• All the presented algorithms of CM-family are able to neutralize the port-starboard ambiguity.
• Doppler measurements give a significant performance boost for CMNF estimates not only for the velocity parameters, but also for the target coordinates and normal acceleration.
• In contrast with the static case, where the introduction of the additional measurements’ transform significantly increases the performance of the estimate, here, the improvement can’t be called dramatic.

5. Conclusions

In the paper, we presented the results of the CMNF algorithm applied to the solution of the maneuvering AUV tracking problem, given the combination of the bearing-Doppler measurements. The measurement structure admitted port-starboard ambiguity. The mathematical model of the target movement represented a stochastic coordinate turn with random initial parameters: maneuver plane position, target heading speed, etc. The total dimensionality of the dynamic system state was 17.

The obtained estimates of the target state, i.e., the coordinates and speed, were unbiased and had high precision. The real estimation error covariance matrix was adequate to the corresponding theoretical values. The port-starboard ambiguity was neutralized. We showed the impact of the Doppler speed measurements on the AUV state estimation performance. It should be noted that the CMNF algorithm, being applied to the bearing-Doppler measurements, allowed us to estimate the initial parameters of the maneuver mentioned above.

We consider the following directions of further research in the same area. First, the optimization of the observations transform is a possible way to the estimate enhancement. Second, the optimization problem of the registration devices configuration and positioning looks promising. Third, the application of the obtained results for the identification of the target maneuver change and successive tracking of the AUV, performing a sequence of maneuvers, is practically important. The research on the mentioned subjects is in progress.