Maneuvering Object Tracking and Movement Parameters Identification by Indirect Observations with Random Delays

Abstract

The paper presents an approach to solving the problem of unknown motion parameters Bayesian identification for the stochastic dynamic system model with randomly delayed observations. The system identification and the object tracking tasks obtain solutions in the form of recurrent Bayesian relations for a posteriori probability density. These relations are not practically applicable due to the computational challenges they present. For practical implementation, we propose a conditionally minimax nonlinear filter that implements the concept of conditionally optimal estimation. The random delays model source is the area of autonomous underwater vehicle control. The paper discusses in detail a computational experiment based on a model that is closely aligned with this practical need. The discussion includes both a description of the filter synthesis features based on the geometric interpretation of the simulated measurements and an impact analysis of the effectiveness of model special factors, such as time delays and model unknown parameters. Furthermore, the paper puts forth a novel approach to the identification problem statement, positing a random jumping change in the motion parameters values.

1. Introduction

A stochastic filtering theory provides methods for solving very different practical problems. The moving objects control tasks, including the most interesting ones such as target tracking and navigation [1], attract our attention to a greater extent. An atypical property of the models used for such tasks is the observations random time delays [2,3]. The source of these delays is the aquatic environment, and autonomous underwater vehicles (AUVs) are a practical application [4]. This rapidly developing applied area regularly generates new and original challenges for researchers. For example, sea currents are a characteristic phenomenon of the aquatic environment, and they have a specific effect on the most typical analysis algorithms [5]. The water peculiarities can affect not only the conditions but also the movement goals. At the same time, it seems that researchers pay more attention to the target tracking models and tasks than to others [6,7,8,9,10,11,12,13,14]. But first of all, the specificity of the water affects the measuring instruments. Within the cooperative scenario framework, when the AUV interacts with measuring instruments to solve the positioning problem, different measuring instruments are available for use [15,16], all based on common physical laws. These are acoustic sensors, which means that their performance is influenced by various factors, such as temperature, salinity, and water pressure [17,18,19]. As a result, the observer has a time delay in AUV’s position characteristics receiving data. This delay is random and may turn out to be quite large since the propagation velocity of the acoustic wave is not that high.

The lack of accurate measurement-performing-moment information is the studied model’s main feature. These kinds of problems are well known and not unique. Typical situations arise in the use of distributed telecommunications and global positioning systems [20,21], as well as in AUV control [22]. Generally, the time out of sync between data sources is an obstacle in the work. Among the tasks solved in such conditions, there is also a tracking problem whose solution is provided by the use of nonlinear filters [23,24]. The fundamental difference of the considered here model is that the desynchronization time is not only unknown and random, but it also changes during observation.

The idea of researching a stochastic dynamic observation system model that takes into account the acoustic signal delay factor appeared in a study described in [2]. The article in [3] presented a detailed computational experiment with the tracking problem-solving methodology based on the concept of conditionally optimal filtering [25,26] and on its evolution—conditionally minimax nonlinear filtering (CMNF) [27]. This work takes the model and experiment of [3] as a basis and adapts them to more complex conditions.

Here, we focus on the issue of AUV movement model identification. This is a question of the assumption that the moving object state model is known, which is typical for stochastic observation systems. In practice, this does not happen, but we use models with certain errors. To compensate for the impact of these errors, an error model can be included in the state equations. Researchers working in the field of robust estimation [28] have studied a large number of tasks and methods of this kind. Another approach, known as parametric systems identification [29], suggests describing model inaccuracies with parameters and evaluating them. Having obtained parameter estimates, we can consider the system identified, i.e., use parameter estimates in the motion model. The identified system can be used further to solve basic tasks, positioning, and target tracking, i.e., to use filters. Because pre-identification is not always possible, it is quite natural to combine the filtering and identification tasks. It is convenient to do this within the Bayesian approach framework, i.e., under the assumption that unknown model parameters are random. This is a widely used and well-known method that has been around for a long time [30,31]. In this article, we apply it to a model with random observation delays (Section 2). As a result, we present the optimal Bayesian filter equations (Section 3), the adapted for identification CMNF structure (Section 4), and a detailed analysis of the computational experiment, reproducing the calculations performed in [3] in the model with unknown AUV movement parameters.

2. Filtering and Model Parameters Identification Problems Statement

Three formal variants of the discrete stochastic dynamic observation system model describe motion and measuring instruments in [3]. The simple model that assumes additive disturbances in dynamics and additive noises in observations is most suitable for the practical purposes set out in this paper. It quite adequately reflects the real physics in the target application area. To the basic model, we will add the assumption of motion model incomplete information, for which we will include unknown parameters in the state equations. In addition, the observation model requires a priori knowledge of the time delay maximum possible value $T.$ It means that the reflected acoustic wave object is located at a quite large distance. Therefore, it takes some time for the acoustic signal to reach the sensor that receives it. This time is unknown and random, but, according to the proposed model, it does not exceed the value of $T.$ Since the sound velocity in the water can be considered known, this assumption boils down to estimating the maximum distance between the sensors and the observed object, which should not be difficult. Next, we assume that the object tracking start moment is $0$ and a discrete variable $t$ describes the time. This variable should take the values $t=-T,-T+1,\dots ,0,1,\dots$ Then, at the moment $t=0$, measurements corresponding to the object positions at all previous time points down to and including $-T$ will be guaranteed to be formed.

Let us consider the object position, denoted by the vector $x_t={\left({x_1}_t,\dots ,{x_p}_t\right)}^{\mathrm{\prime }}\in {\mathbb{R}}^p,$ and the unknown motion model parameters vector ${\mu =\left({\mu }_1,\dots ,{\mu }_r\right)}^{\prime }\in {\mathbb{R}}^r$. Here and hereafter, we use «$\prime \mathrm{»}$ to denote a transposed vector or matrix.

So, the position $x_t$ can consist of an object coordinate vector ${\left(x\left(t\right),y\left(t\right),z\left(t\right)\right)}^{\prime }$ and its velocity vector ${\left(v_x\left(t\right),v_y\left(t\right),v_z\left(t\right)\right)}^{\prime }$ in a geocentric inertial reference frame. The parameters $\mu$ do not change values during the observation time and may include, for example, an unknown initial position and/or average velocities. The position–parameter pair forms the dynamic system state vector ${\left(x_t^{\prime },{\mu }^{\prime }\right)}^{\prime }\in {\mathbb{R}}^{p+r}$.

Following the Bayesian approach, we assume the unknown parameter vector $\mu$ to be random and its distribution to be given. We will also consider the initial position of the object to be random, denoted by $\eta \in {\mathbb{R}}^p$. Given that the time variable $t$ starts counting time from the value $t=-T$, the initial position must be determined as $x_{-T-1}=\eta$. Next, we assume that $\mu$ and $\eta$ are independent and have probability densities $f_{\eta }$ and $f_{\mu }$, so the initial Bayesian filter condition (vector ${\left(x_{-T-1}^{\prime },{\mu }^{\prime }\right)}^{\prime }$) probability density $f_0$ takes the form

$$f_0\left(X_{-T-1},M\right)=f_{\eta }\left(X_{-T-1}\right)f_{\mu }\left(M\right),$$

where $X_{-T-1}\in {\mathbb{R}}^p$ and $M\in {\mathbb{R}}^r$ are the probability density arguments corresponding to random vectors $x_{-T-1}$ and $\mu$.

The observation vector $={\left({y_1}_t,\dots ,{y_q}_t\right)}^{\prime }\in {\mathbb{R}}^q,$ $t=0,1,\dots ,$ consists of the $q$ position $x_t$ measurements at different times, shifted from the current time by a random delay. Moreover, each element ${y_i}_t$ has its own offset, which is determined by a discrete random variable ${{\tau }_i}_t\in \left\{0,1,\dots ,T\right\}$. This means that the observation $y_t, \mathrm{r}\mathrm{e}\mathrm{c}\mathrm{e}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{d} \mathrm{a}\mathrm{t} \mathrm{t}\mathrm{h}\mathrm{e} \mathrm{m}\mathrm{o}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{t} t,$ consists of object state measurements at different time points $x_{t-{{\tau }_1}_t},\dots ,x_{t-{{\tau }_q}_t}$.

Values ${{\tau }_i}_t$ make a vector ${\tau }_t={\left({{\tau }_1}_t,\dots ,{{\tau }_q}_t\right)}^{\prime }\in {\mathbb{R}}^q$ that defines a known function ${\theta }_t\left(x_t\right)$. Thus, the dynamic system has the form:

$$\begin{array}{c}{x}_t={\phi }_t\left(x_{t-1},\mu \right)+w_t, x_{-T-1}=\eta ,\\ \begin{array}{c}{y_i}_t={\psi }_t\left(x_{t-{{\tau }_i}_t}\right)+{v_i}_t, i=1,\dots ,q,\\ {\tau }_t={\theta }_t\left(x_t\right),\end{array}\end{array}$$

(1)

where additive dynamic disturbances $w_t={({w_1}_t,\dots ,{w_p}_t)}^{\prime }\in {\mathbb{R}}^p$ and additive measurement noises $v_t={({v_1}_t,\dots ,{v_q}_t)}^{\prime }\in {\mathbb{R}}^q$ are discrete-time vector white noises having known distributions. The vectors $\eta ,\mu ,w_t,v_t$ are assumed to be independent in aggregate, and the functions ${\phi }_t$, ${\psi }_t$ and ${\theta }_t$ satisfy sufficient conditions for the filtering estimates existence from Theorems 1 and 2 in [3]. These are typical conditions for the growth rate of the functions and the disturbances distributions, which guarantee the existence of finite second moments of the processes $x_t$ and $y_t$.

Let us introduce a new composite state vector ${\left(x_t^{\prime },{\mu }_t^{\prime }\right)}^{\prime }$, combining the object position vector $x_t$ and unknown parameter vector $\mu$. To do this, we determine that under the designation $x_{t-{\tau }_t},$it refers to the vector of positions ${\left(x_{t-{{\tau }_1}_t}^{\prime },\dots ,x_{t-{{\tau }_q}_t}^{\prime }\right)}^{\prime }$ shifted in accordance with the delays ${\tau }_t={\left({{\tau }_1}_t,\dots ,{{\tau }_q}_t\right)}^{\prime }$. Now the observation system model (1) takes the following traditional dynamic system form:

$$\begin{array}{c}{x}_t={\phi }_t\left(x_{t-1},{\mu }_{t-1}\right)+w_t, x_{-T-1}=\eta ,\\ \begin{array}{c}{\mu }_t={\mu }_{t-1}, {\mu }_{-T-1}=\mu ,\\ y_t={\psi }_t\left(x_{t-{\tau }_t}\right)+v_t, {\tau }_t={\theta }_t\left(x_t\right).\end{array}\end{array}$$

(2)

For (2), we can state the filtering problem, i.e., the estimation problem for the state ${\left(x_t^{\prime },{\mu }_t^{\prime }\right)}^{\prime }$ by the observations $y_s,s=0,1,\dots ,t.$ Solving the filtering problem for this composite state vector, we will get both a tracking problem solution in the form of a current position $x_t$ filtering estimate and a solution to the movement model parameter identification in the form of a filtering estimate of the state vector second part ${\mu }_t$. This estimate is a parameter $\mu$ Bayesian estimate.

It remains to be noted that, in estimating the problem under consideration, the evaluation accuracy criterion for some estimates ${\left({\breve{x}}_t^{\prime },{\breve{\mu }}_t^{\prime }\right)}^{\prime }$ is the mean square, i.e., $\mathrm{E}\left\{{‖{\left(x_t^{\prime },{\mu }_t^{\prime }\right)}^{\prime }-{\left({\breve{x}}_t^{\prime },{\breve{\mu }}_t^{\prime }\right)}^{\prime }‖}^2\right\}$ where $\mathrm{E}\left\{x\right\}$ is the vector $x$ mathematical expectation and $‖x‖$ is the vector $x$ usual Euclidean norm. The Bayesian estimates $x_t^{\mathrm{*}}$ and ${\mu }_t^{\mathrm{*}}$ obtained in the next section represent conditional mathematical expectations

$$x_t^{\mathrm{*}}=\mathrm{E}\left\{x_t|{\mathcal{F}}_t^y\right\}, {\mu }_t^{\mathrm{*}}=\mathrm{E}\left\{\mu |{\mathcal{F}}_t^y\right\},$$

where ${\mathcal{F}}_t^y$ is the $\sigma$-algebra, generated by observations $y_s,s=0,1,\dots ,t.$

The conditional minimax filter estimates ${\left({\left({\widehat{x}}_t\right)}^{\prime },{\left({\widehat{\mu }}_t\right)}^{\prime }\right)}^{\prime }$ used further as an effective computational alternative turn out to be the closest to Bayesian in the sense of the specified quadratic criterion on a class of the given structure filters. The main attention in this filter synthesis is paid precisely to the suitable structure choice that takes into account the specific features of model (2). Section 4 of the paper contains a detailed description of this choice.

The last feature, which we will add to model (2), affects the interpretation of the problem of the tracking object movement model parameter identification. The traditional task of the system parametric identification [29] is to estimate unknown parameters. This means that it is possible to obtain an unlimited number of observations so that the current parameter estimate can be improved with each new observation, achieving estimate convergence to the true parameter value. This identifies the system. In reality, of course, the observation scope is limited. In addition, the assumption of the identified parameters constancy is also valid generally for a limited time. Thus, movement conditions change from time to time, for example, under the influence of external control factors, after which another model has to be identified. There are different ways to take such situations into account. This paper proposes the following model.

Suppose, in addition to the position $x_t$, there is a Poisson process with a known intensity. It can be enabled by extending the state as follows. For the sake of certainty, we will assume that ${x_p}_t={\lambda }_{u}P\left(u\right)$ where $P\left(u\right)$ is a standard Poisson process and the variable $u=\delta t$ denotes continuous time at the same observation interval as $t$. That is, we assume that the tracking interval $[0,T^{end}]$ is sampled in increments of $\delta$. Let ${\mu }_t^P$ be a discrete-time white noise with distribution $f_{\mu }$. Then, instead of Equation (2), we will consider a dynamic observation system model of the following form:

$$\begin{array}{c}{x}_t={\phi }_t\left(x_{t-1},{\mu }_{t-1}\right)+w_t, x_{-T-1}=\eta ,\\ \begin{array}{c}{\mu }_t=\left(1-{x_p}_t+{x_p}_{t-1}\right){\mu }_{t-1}+\left({x_p}_t-{x_p}_{t-1}\right){\mu }_t^P, {\mu }_{-T-1}={\mu }_{-T-1}^P,\\ y_t={\psi }_t\left(x_{t-{\tau }_t}\right)+v_t, {\tau }_t={\theta }_t\left(x_t\right).\end{array}\end{array}$$

(3)

Note that, formally, model (3) can be reduced to form (2) if we limit the Poisson process $P\left(u\right)$ switching jumps number. Then, it is enough to expand the vector ${\mu }_t$ putting ${\mu }_t={\left({\left({\mu }^{\left(0\right)}\right)}^{\prime },{\left({\mu }^{\left(1\right)}\right)}^{\prime },\dots ,{\left({\mu }^{\left(max\right)}\right)}^{\prime }\right)}^{\prime }$ and meaning that ${\mu }^{\left(0\right)}={\mu }_{-T-1}^P$ are initial unknown movement model parameters, ${\mu }^{\left(1\right)}$ is unknown movement model parameters generated at the first $P\left(u\right)$ jump, ${\mu }^{\left(2\right)}$ is unknown movement model parameters generated at the second $P\left(u\right)$ jump, etc. The number of jumps is limited by the value $\left(max\right).$ In this case, all ${\mu }^{\left(i\right)}$ should be independent random vectors with the distribution $f_{\mu }$. Then, from the filtering point of view, both Bayesian and conditionally minimax, the estimation statements in model (3) and in model (2) do not differ. From a practical point of view, the difference between these models is fundamental. In model (3), it is impossible to achieve the system identification goal by parameter estimation with any accuracy. First, unknown model parameters change within the same trajectory. Second, there is only a finite number of observations to parameter estimation at any interval of their constancy. Thus, when performing experimental or practical model (3) calculations, it is necessary to estimate the unknown current parameters ${\mu }^{\left(i\right)},$ estimate the $i$-th interval of constancy boundaries of the switching process $P\left(u\right),$ adapt to the current movement model within the interval, and form a accurate position $x_t$ estimate, providing a solution to the main target-tracking task.

It may be noted right away that such a movement observation model complication will not necessarily make the filtering algorithms more difficult. Indeed, a certain physical meaning is embedded in the model, namely, $x_t$ is the target position coordinates and $y_t$ is the data of this position measurements, for example, bearings and ranges. Therefore, the parameter ${\mu }_t$ importance may not be so great since these parameters have a greater impact on the object’s future positions, and for filtering, it is much more important to observe the estimating current position. Such a situation occurs in the computational experiment discussed below but only in the absence of time delays, i.e., when ${\tau }_t=0$. As soon as ${\tau }_t$ begins to differ significantly from $0,$ the prediction accuracy, for which ${\mu }_t$ estimates are needed, begins to contribute significantly to the target-tracking problem-solving.

3. Optimal Bayesian Filtering and Identification

In order to obtain the Bayesian estimate ${\left({\left(x_t^{\mathrm{*}}\right)}^{\prime },{\left({\mu }_t^{\mathrm{*}}\right)}^{\prime }\right)}^{\prime }$ for model (1) or (2), i.e., to obtain the relations for the a posteriori probability density [32], we need to repeat the proof of Theorem 1 from [3]. There are no fundamental difficulties here, so we will describe it briefly, omitting the details.

To write the filter, an extended position vector ${\mathbf{x}}_t={\left(x_{t-T}^{\prime },\dots ,x_{t-1}^{\prime },x_t^{\prime }\right)}^{\prime }$ is needed since all possible time delays must be taken into account. The desired a posteriori probability density of the vector ${\left({\mathbf{x}}_t^{\prime },{\mu }^{\prime }\right)}^{\prime }$ given the observation vector $y^t={\left(y_0^{\prime },\dots ,y_t^{\prime }\right)}^{\prime }$ denotes ${\rho }_t={\rho }_t\left({\mathbf{X}}_t,M|Y^t\right)$. For a more concise writing form, we also introduce an unnormalized a posteriori density ${\widehat{\rho }}_t={\widehat{\rho }}_t\left({\mathbf{X}}_t,M|Y^t\right)$ so that

$$\begin{array}{c}{\rho }_t\left({\mathbf{X}}_t,M|Y^t\right)={\displaystyle \frac{{\widehat{\rho }}_t\left({\mathbf{X}}_t,M|Y^t\right)}{\int {\widehat{\rho }}_t\left({\mathbf{X}}_t,M|Y^t\right)d{\mathbf{X}}_{t}dM}}.\\ \end{array}$$

(4)

Let us pay attention to the convenient notation for the a posteriori density arguments: ${\mathbf{X}}_t={\left(X_{t-T}^{\prime },\dots ,X_{t-1}^{\prime },X_t^{\prime }\right)}^{\prime }$, $M={\left(M_1,\dots ,M_r\right)}^{\prime }$ and $Y^t{=\left(Y_0^{\prime },\dots ,Y_t^{\prime }\right)}^{\prime }$. They correspond to vectors ${\mathbf{x}}_t$, $\mu$ и $y^t$.

The extended position vector ${\mathbf{x}}_t\in {\mathbb{R}}^{\left(T+1\right)p}$ and the equation for ${\mathbf{x}}_t$ has the following form:

$$\begin{array}{c}{\left({\mathbf{x}}_t\right)}_1^p={\left({\mathbf{x}}_{t-1}\right)}_{p+1}^{2p},\\ \dots \\ \begin{array}{c}{\left({\mathbf{x}}_t\right)}_{\left(T-1\right)p+1}^{Tp}={\left({\mathbf{x}}_{t-1}\right)}_{Tp+1}^{\left(T+1\right)p},\\ {\left({\mathbf{x}}_t\right)}_{Tp+1}^{\left(T+1\right)p}={\phi }_t\left({\left({\mathbf{x}}_{t-1},\mu \right)}_{Tp+1}^{\left(T+1\right)p},\right)+w_t,\end{array}\end{array}\left\{\begin{array}{c}{x}_{t-T}\\ \dots \\ \begin{array}{c}{x}_{t-1}\\ x_t={\phi }_t\left(x_{t-1},\mu \right)+w_t\end{array}\end{array}\right\}$$

where ${\left(\mathbf{x}\right)}_i^j$ denotes a sub-vector of vector $\mathbf{x}$ with elements from $i$-th to $j$-th.

Using the corresponding additive noise ${\mathbf{w}}_t={\left(0^{\prime },\dots ,0^{\prime },x_t^{\prime }\right)}^{\prime },$ we obtain the equation for ${\mathbf{x}}_t$ as:

$${\mathbf{x}}_t={\mathsf{\Phi }}_t\left({\mathbf{x}}_t,\mu \right)+{\mathbf{w}}_t.$$

(5)

Similarly, for observations $y_t$ we have the following equation:

$$y_t={\Theta }_t\left({\mathbf{x}}_t\right){\mathsf{\Psi }}_t\left({\mathbf{x}}_t\right)+v_t,$$

(6)

where the observations matrix function ${\mathsf{\Psi }}_t\left({\mathbf{x}}_t\right)\in {\mathbb{R}}^{\left(T+1\right)\times q}$ is as follows:

$${\mathsf{\Psi }}_t\left({\mathbf{x}}_t\right)=\left(\begin{array}{c}{\psi }_t^{\prime }\left({\left({\mathbf{x}}_t\right)}_1^p\right)\\ {\psi }_t^{\prime }\left({\left({\mathbf{x}}_t\right)}_{\left(T-1\right)p+1}^{Tp}\right)\\ \begin{array}{c}\begin{array}{c}\dots \\ {\psi }_t^{\prime }\left({\left({\mathbf{x}}_t\right)}_{Tp+1}^{\left(T+1\right)p}\right)\end{array}\\ \end{array}\end{array}\right)=\left(\begin{array}{c}\begin{array}{c}{\psi }_t^{\prime }\left(x_{t-T}\right)\\ \end{array}\\ \begin{array}{c}{\psi }_t^{\prime }\left(x_{t-T+1}\right)\\ \dots \end{array}\\ \begin{array}{c}{\psi }_t^{\prime }\left(x_t\right)\\ \end{array}\end{array}\right)$$

and delays matrix ${\Theta }_t\left({\mathbf{x}}_t\right)\in {\mathbb{R}}^{q\times \left(T+1\right)}$ is defined as:

$$\begin{gathered} {\left({\Theta }_t\left({\mathbf{x}}_t\right)\right)}_{i,j}=1, \mathrm{if} {\left({\theta }_t\left(x_t\right)\right)}_i=T-j+1, \\ {\left({\Theta }_t\left({\mathbf{x}}_t\right)\right)}_{i,j}=0, \mathrm{else} \\ i=1,\dots ,q, j=1,\dots ,T+1. \end{gathered}$$

It means that the element ${\left({\Theta }_t\left({\mathbf{x}}_t\right)\right)}_{i,j}$ takes the value $1$ if the observation ${y_i}_t$ corresponds to the position $x_{T-j+1}.$ Here we can enumerate the values ${\left[{\Theta }_t\right]}_k$, $k=1,\dots ,q\left(T+1\right),$ of the ${\Theta }_t\left({\mathbf{x}}_t\right)$ putting the sequence number $k$:

$$k=j_{k,1}+\left(T+1\right)j_{k,2}+\dots +{\left(T+1\right)}^{q-1}{j}_{k,q},$$

where $j_{k,i}$ is the row number of the $i$–th column of the matrix ${\left[{\Theta }_t\right]}_k$, in which there is $1$. Exactly the same, the possible values ${\left[{\theta }_t\left(x_t\right)\right]}_k={\left[{\tau }_t\right]}_k$ of the time delay vector ${\tau }_t$ are numbered as follows:

$${\left({\left[{\tau }_t\right]}_k\right)}_i=T-j_{k,i}+1, i=1,..,q.$$

For random vectors $x\in {\mathbb{R}}^p$ and $y\in {\mathbb{R}}^q$, the probability densities used below will be denoted as follows: $f_x\left(X\right)$ is the marginal density of $x$, $f_{x,y}(X,Y)$ is the density of the composite vector ${\left(x^{\prime },y^{\prime }\right)}^{\prime }$, and $f_{x|y}\left(X\right|Y)$ is the conditional density of $x$ given $y$ under the additional assumption $f_y\left(Y\right)>0$.

To find the conditional distribution of the state ${\left({\mathbf{x}}_t^{\prime },{\mu }^{\prime }\right)}^{\prime }$ defined by (5) given the observations $y_s,s=0,\dots ,t,$ defined by (6), we write down recurrent Bayesian relations for an unnormalized a posteriori probability density by analogy with the classical solution [32].

Let us constitute the desired a posteriori density in the following form:

$$\begin{array}{c}\begin{array}{c}{f}_{{\mathbf{x}}_t,\mu |y^t}\left({\mathbf{X}}_t,M|Y^t\right)={\displaystyle \frac{f_{{\mathbf{x}}_t,\mu ,y^t}\left({\mathbf{X}}_t,M,Y^t\right)}{f_{y^t}\left(Y^t\right)}}={\displaystyle \frac{f_{{\mathbf{x}}_t,\mu ,y^{t-1},y_t}\left({\mathbf{X}}_t,M,Y^{t-1},Y_t\right)}{f_{y^t}\left(Y^t\right)}}=\\ ={\displaystyle \frac{f_{{\mathbf{x}}_t,\mu ,y^{t-1}}\left({\mathbf{X}}_t,M,Y^{t-1}\right)f_{y_t|{\mathbf{x}}_t,\mu }\left(Y_t|{\mathbf{X}}_t,M\right)}{f_{y^t}\left(Y^t\right)}}=\end{array}\\ \begin{array}{c}={\displaystyle \frac{f_{{\mathbf{x}}_t,\mu ,y^{t-1}}\left({\mathbf{X}}_t,M,Y^{t-1}\right)f_{y_t|{\mathbf{x}}_t}\left(Y_t|{\mathbf{X}}_t\right)}{f_{y^t}\left(Y^t\right)}}.\end{array}\end{array}$$

In the last equation, it is taken into account that in (6), there is no explicit dependence of $y_t$ on $\mu ,$ which means equality $f_{y_t|{\mathbf{x}}_t,\mu }\left(Y_t|{\mathbf{X}}_t,M\right)=f_{y_t|{\mathbf{x}}_t}\left(Y_t|{\mathbf{X}}_t\right).$

Next, for the first multiplier $f_{{\mathbf{x}}_t,\mu ,y^{t-1}}\left({\mathbf{X}}_t,M,Y^{t-1}\right)$, we perform the following transformations:

$$\begin{gathered} \begin{array}{c}{f}_{{\mathbf{x}}_t,\mu ,y^{t-1}}\left({\mathbf{X}}_t,M,Y^{t-1}\right)=f_{x_{t-T},\dots ,x_t,\mu ,y^{t-1}}\left(X_{t-T},\dots ,X_t,M,Y^{t-1}\right)=\hfill \\ \\ =\int f_{x_{t-T-1},x_{t-T},\dots ,x_t,\mu ,y^{t-1}}\left(X_{t-T-1},X_{t-T},\dots ,X_t,M,Y^{t-1}\right)d{X}_{t-T-1}=\hfill \\ \\ =\int f_{x_t|{\mathbf{x}}_{t-1},\mu ,y^{t-1}}\left(X_t|{\mathbf{X}}_{t-1},M,Y^{t-1}\right)f_{{\mathbf{x}}_{t-1},\mu ,y^{t-1}}\left({\mathbf{X}}_{t-1},M,Y^{t-1}\right)d{X}_{t-T-1}=\hfill \\ \\ =\int f_{x_t|x_{t-1},\mu }\left(X_t|X_{t-1},M\right)f_{{\mathbf{x}}_{t-1},\mu |y^{t-1}}\left({\mathbf{X}}_{t-1},M|Y^{t-1}\right)d{X}_{t-T-1}·f_{y^{t-1}}\left(Y^{t-1}\right)=\hfill \\ \\ =\int f_{{\mathbf{x}}_{t-1},\mu |y^{t-1}}\left({\mathbf{X}}_{t-1},M|Y^{t-1}\right)d{X}_{t-T-1}·f_{x_t|x_{t-1},\mu }\left(X_t|X_{t-1},M\right)f_{y^{t-1}}\left(Y^{t-1}\right)=\hfill \\ \\ =\int f_{{\mathbf{x}}_{t-1},\mu |y^{t-1}}\left({\mathbf{X}}_{t-1},M|Y^{t-1}\right)d{X}_{t-T-1}·f_{w_t}\left(X_t-{\phi }_t\left(X_{t-1},\mu \right)\right)f_{y^{t-1}}\left(Y^{t-1}\right).\hfill \end{array} \end{gathered}$$

For the second multiplier $f_{y_t|{\mathbf{x}}_t}\left(Y_t|{\mathbf{X}}_t\right)$, it holds

$$\begin{gathered} \begin{array}{c}\hfill f_{y_t|{\mathbf{x}}_t}\left(Y_t|{\mathbf{X}}_t\right)={\displaystyle \sum _{\begin{array}{c}k=1\\ \end{array}}^{q\left(T+1\right)}}I\left({\Theta }_t\left({\mathbf{X}}_t\right)={\left[{\Theta }_t\right]}_k\right)f_{y_t|{\mathbf{x}}_t}\left(Y_t|{\mathbf{X}}_t,{\Theta }_t\left({\mathbf{X}}_t\right)={\left[{\Theta }_t\right]}_k\right)=\hfill \\ \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \sum _{\begin{array}{c}k=1\\ \end{array}}^{q\left(T+1\right)}}I\left({\theta }_t\left(X_t\right)={\left[{\tau }_t\right]}_k\right)f_{y_t|{\mathbf{x}}_t,{\tau }_t={\left[{\tau }_t\right]}_k}\left(Y_t|{\mathbf{X}}_t,{\tau }_t={\left[{\tau }_t\right]}_k\right)=\hfill \\ \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \sum _{\begin{array}{c}k=1\\ \end{array}}^{q\left(T+1\right)}}I\left({\theta }_t\left(X_t\right)={\left[{\tau }_t\right]}_k\right)f_{v_t}\left(Y_t-{\left[{\Theta }_t{\mathsf{\Psi }}_t\left({\mathbf{X}}_t\right)\right]}_k\right).\hfill \end{array} \end{gathered}$$

And finally, we get

$$f_{{\mathbf{x}}_t,\mu |y^t}={\displaystyle \frac{f_{w_t}\left(X_t-{\phi }_t\left(X_{t-1},\mu \right)\right)\sum _{k}I\left(\begin{array}{c}{\theta }_t\left(X_t\right)=\\ ={\left[{\tau }_t\right]}_k\end{array}\right)f_{v_t}\left(Y_t-{\left[{\Theta }_t{\mathsf{\Psi }}_t\left({\mathbf{X}}_t\right)\right]}_k\right)\int f_{{\mathbf{x}}_{t-1},\mu |y^{t-1}}d{X}_{t-T-1}}{f_{y^t}\left(Y^t\right)/f_{y^{t-1}}\left(Y^{t-1}\right)}}.$$

The numerator of the obtained equality is the desired non-normalized a posteriori density (4), since the denominator performs only the function of normalization, i.e.,

$$\begin{array}{c}\begin{array}{c}\begin{array}{c} \\ {\widehat{\rho }}_t\left({\mathbf{X}}_t,M|Y^t\right)=f_{w_t}\left(X_t-{\phi }_t\left(X_{t-1},\mu \right)\right)\times \end{array}\\ \times {\displaystyle \sum _{k=1}^{q\left(T+1\right)}}I\left(\begin{array}{c}{\theta }_t\left(X_t\right)=\\ ={\left[{\tau }_t\right]}_k\end{array}\right)f_{v_t}\left(Y_t-{\left[{\Theta }_t{\mathsf{\Psi }}_t\left({\mathbf{X}}_t\right)\right]}_k\right)\int f_{{\mathbf{x}}_{t-1},\mu |y^{t-1}}d{X}_{t-T-1}.\end{array}\end{array}$$

(7)

Recall that the defined above variable $k$ numbers the values of ${\left[{\tau }_t\right]}_k$ and ${\left[{\Theta }_t{\mathsf{\Psi }}_t\left({\mathbf{X}}_t\right)\right]}_k$. Calculations start with the following initial condition:

$$\begin{array}{c}\begin{array}{c}{\widehat{\rho }}_{-1}\left({\mathbf{X}}_{-1},M|Y^{-1}\right)={\widehat{\rho }}_{-1}\left({\mathbf{X}}_{-1},M\right)={\widehat{\rho }}_{-1}\left(X_{-T-1},\dots ,X_{-1},M\right)=\\ =f_0\left(X_{-T-1},M\right)f_{w_{-T}}\left(X_{-T}-{\phi }_{-T}\left(X_{-T-1},M\right)\right)·\dots ·f_{w_{-1}}\left(X_{-1}-{\phi }_{-T}\left(X_{-2},M\right)\right).\end{array}\\ \end{array}$$

Here we used the notations $f_{w_t}\left(W_t\right)$, $t=-T,-T+1,\dots ,$ for the dynamic model disturbance probability density and $f_{v_t}\left(V_t\right),$ $t=\mathrm{0,1},\dots$ for the observation error probability density.

After solving (7) and normalizing ${\widehat{\rho }}_t,$ the Bayesian estimates $x_t^{\mathrm{*}}$ and ${\mu }_t^{\mathrm{*}}$ are obtained by integrating:

$${\left({\left({\mathit{x}}_t^{\mathrm{*}}\right)}^{\prime },{\left({\mu }_t^{\mathrm{*}}\right)}^{\prime }\right)}^{\prime }=\int {\left({\mathit{X}}_t^{\prime },M^{\prime }\right)}^{\prime }{\rho }_t\left({\mathit{X}}_t,M|Y^t\right)d{\mathit{X}}_{t}dM, x_t^{\mathrm{*}}={\left({\mathbf{x}}_t^{\mathrm{*}}\right)}_{Tp+1}^{\left(T+1\right)p}.$$

4. Conditionally Minimax Nonlinear Filter

So, there is a fundamental possibility to obtain optimal estimates. This is the first conclusion from the optimal Bayesian filtering and identification relations obtained in the previous section. However, it is not possible to rely on their practical application. This is how we get the second conclusion. Indeed, even in the case of an AUV moving in a plane at a constant velocity for tracking, we would need to implement Formula (7) calculation for dimensions $p=2,$ $r=2,$ $q=4,$ and $T=50,$ which gives the dimension of the vector ${\mathbf{X}}_t$ equal to $200$ and the same number of values ${\left[{\tau }_t\right]}_k.$ It seems, at the least, rash to hope for success in calculating integrals for such dimensions with good accuracy in the real-time mode receipt of observations $y_t$.

If we are unable to implement an optimal solution, then we need to use a suboptimal estimate that has good accuracy and is practically feasible. There are many methodologies that could provide such a tool. The choice can be made from various universal solutions, such as the extended Kalman filter [33], the Gaussian sum filter [34], particle filters [35,36], and sigma point filters [37,38]. Specialized methods, such as pseudo-measurement filters [39], may be more promising in such tasks, as demonstrated in [40]. However, all these methods are intended for use in the models (5) and (6) and cannot be applied in model (3) without refinement. The huge dimensions that prevent the optimal filter use will also lead to great difficulties when using known suboptimal filters. Moreover, the characteristic and long-known tendency of such filters to diverge will exacerbate the complexity of calculations.

Only the concept of conditionally optimal filtering [25,26] allows usage to the proposed random time observations delays model without deep refinement. An important part of this concept is the initial filter structure flexibility. Anyone can take into account physical, geometric, and other laws in force in a particular observation system. The final filter structure obtained as a result of its parameter optimization guarantees the estimation of the statistical properties, namely, non-bias, non-divergence, and superiority over the a priori estimate. In the problem under consideration, the flexible filter structure allows for taking into account not only the laws governing the AUV movement and measurements but also time delays. In addition, the conditionally optimal filter optimization is grounded on the well-known Gaussian distribution minimax property [27], which provides robustness to the Kalman filter [41]. Complementing all this with the filter parameters calculation using computer simulation, we obtain effective, high-quality, and realizable solutions to both the target-tracking and movement parameter identification problems.

Based on the conditional minimax terminology [27] and the designations used in the target-tracking statement in [3], we briefly describe the CMNF synthesis stages, adapting the description to the problem under consideration for the system (1) or (3).

The desired filtering estimate ${\left({\left({\widehat{x}}_t\right)}^{\prime },{\left({\widehat{\mu }}_t\right)}^{\prime }\right)}^{\prime }$ combines the position $x_t$ estimate and the movement model random parameter $\mu$ estimate given the observations of $y^t.$ This estimate has the standard prediction–correction form:

$$\begin{array}{c}\begin{array}{c}{\widehat{x}}_t={\tilde{x}}_t+∆{\widehat{x}}_t,\end{array}\\ \begin{array}{c}{\widehat{\mu }}_t={\widehat{\mu }}_{t-1}+∆{\widehat{\mu }}_t.\end{array}\end{array}$$

Here, to compute a position prediction ${\tilde{x}}_t$, it is necessary to choose the basic prediction function ${\xi }_t\left(X,M\right),$ and to compute a correction, it is necessary to choose the basic correction function ${\zeta }_t\left(X,Y\right)$, ${X\in \mathbb{R}}^p,$ ${M\in \mathbb{R}}^r$, ${Y\in \mathbb{R}}^q$. After selecting them, the prediction ${\tilde{x}}_t$ is calculated as the ${\xi }_t={\xi }_t\left({\widehat{x}}_{t-1},{\widehat{\mu }}_{t-1}\right)$ transformation minimizing the distance to $x_t;$ the correction $∆{\widehat{x}}_t$ is calculated as the ${\zeta }_t={\zeta }_t\left({\tilde{x}}_t,y_t\right)$ transformation minimizing the distance to $x_t-{\tilde{x}}_t.$ Distance has mean square sense, as mentioned above.

Note, the correction ${\zeta }_t$ does not use the parameter estimate ${\widehat{\mu }}_{t-1}$ (there is no argument $M$), which is explained by the absence of the parameters in the observation model.

Implementing ${\tilde{x}}_t$ and $∆{\widehat{x}}_t$ functions is linear:

$$\begin{array}{c}{\tilde{x}}_t=F_t{\xi }_t\left({\widehat{x}}_{t-1},{\widehat{\mu }}_{t-1}\right)+f_t,\\ \left(\begin{array}{c}{\widehat{x}}_t\\ {\widehat{\mu }}_t\end{array}\right)=\left(\begin{array}{c}{\tilde{x}}_t\\ {\widehat{\mu }}_{t-1}\end{array}\right)+H_t{\zeta }_t\left({\tilde{x}}_t,y_t\right)+h_t,\end{array}$$

(8)

where

$$\begin{array}{c}{F}_t=\mathrm{c}\mathrm{o}\mathrm{v}\left(x_t,{\xi }_t\right){\mathrm{c}\mathrm{o}\mathrm{v}}^{+}\left({\xi }_t,{\xi }_t\right),f_t=\mathrm{E}\left\{x_t\right\}-F_t\mathrm{E}\left\{{\xi }_t\right\},\\ H_t=\mathrm{c}\mathrm{o}\mathrm{v}\left(\left(\begin{array}{c}{x}_t\\ \mu \end{array}\right)-\left(\begin{array}{c}{\tilde{x}}_t\\ {\widehat{\mu }}_{t-1}\end{array}\right),{\zeta }_t\right){\mathrm{c}\mathrm{o}\mathrm{v}}^{+}\left({\zeta }_t,{\zeta }_t\right),h_t=-H_t\mathrm{E}\left\{{\zeta }_t\right\}.\end{array}$$

(9)

In (9), in addition to the mathematical expectation notation $\mathrm{E}\left\{x\right\},$ $\mathrm{c}\mathrm{o}\mathrm{v}\left(x,y\right)$ is the covariance of $x$ and ${y,}^{\mathrm{«}+\mathrm{»}}$ is Moore–Penrose pseudoinversion [42]. At the same time, the position prediction ${\tilde{x}}_t,$ estimates ${\widehat{x}}_t$ and ${\widehat{\mu }}_t$ are unbiased with the following estimation error covariances:

$$\begin{array}{c}{\tilde{K}}_t=\mathrm{c}\mathrm{o}\mathrm{v}\left(x_t-{\tilde{x}}_t,x_t-{\tilde{x}}_t\right)=\mathrm{c}\mathrm{o}\mathrm{v}\left(x_t,x_t\right)-F_t\mathrm{c}\mathrm{o}\mathrm{v}\left({\xi }_t,x_t\right),\\ \begin{array}{c}{\widehat{K}}_t^x=\mathrm{c}\mathrm{o}\mathrm{v}\left(x_t-{\widehat{x}}_t,x_t-{\widehat{x}}_t\right)={\tilde{K}}_t-{∆\widehat{K}}_t^x, {\tilde{K}}_0=\mathrm{c}\mathrm{o}\mathrm{v}\left(\eta ,\eta \right),\\ {\widehat{K}}_t^{\mu }=\mathrm{c}\mathrm{o}\mathrm{v}\left(\mu -{\widehat{\mu }}_t,\mu -{\widehat{\mu }}_t\right)={\widehat{K}}_{t-1}^{\mu }-{∆\widehat{K}}_t^{\mu }, {\tilde{K}}_0=\mathrm{c}\mathrm{o}\mathrm{v}\left(\mu ,\mu \right),\end{array}\end{array}$$

(10)

where matrices ${∆\widehat{K}}_t^x$и ${∆\widehat{K}}_t^{\mu }$ are upper and lower diagonal blocks of the matrix $H_t\mathrm{c}\mathrm{o}\mathrm{v}\left({\zeta }_t,\left(\begin{array}{c}{x}_t\\ \mu \end{array}\right)-\left(\begin{array}{c}{\tilde{x}}_t\\ {\widehat{\mu }}_{t-1}\end{array}\right)\right)=\left(\begin{array}{cc}{∆\widehat{K}}_t^x& \dots \\ \dots & {∆\widehat{K}}_t^{\mu }\end{array}\right)$.

Basic prediction and correction function linear transformations (8) have a minimax justification as the best estimates for classes of all probability distributions with a known mean and covariance [27]. We replace the filter coefficient $F_t$, $f_t$, $H_t$, and $h_t$ analytical calculations by Monte Carlo estimates in practical calculations, i.e., instead of $\mathrm{E}\left\{x\right\}$ in (9), we substitute $\overline{\mathrm{E}}\left\{x\right\}={\displaystyle \frac{1}{N}}\sum _{i=1}^N{x}^{\left(i\right)}$ where ${\left\{x^{\left(i\right)}\right\}}_{i=1}^N$ are $x$ sample values modeled on a computer. Finally, estimate (8) existence of sufficient conditions in relation to model (1) is in [2,3].

Thus, to obtain the target-tracking problem solution, we need to specify movement and observation conditions and synthesize the CMNF by selecting the basic prediction ${\xi }_t$ and correction ${\zeta }_t$ functions.

5. Tracking and Parameter Identification of an AUV Moving in a Plane

5.1. Observation System Model

The paper [3] contains the results of the first practical study of the proposed random delayed observation model and the filtering algorithms (8) and (9). The computational experiment idea for this work is to completely reproduce the same dynamic system and measuring complex, adding only the unknown movement parameter assumption. Therefore, we are considering a simple AUV movement model in a plane with unknown velocity. AUV moves in the horizontal plane $xy$ and has custom coordinates $x\left(t\right)$ and $y\left(t\right)$ (km), i.e., the position is $x_t={\left(x\left(t\right),y\left(t\right)\right)}^{\prime }\in {\mathbb{R}}^2.$ The Gaussian vector $\eta$ gives the initial AUV’s position ${\left(x\left(-T-1\right),y\left(-T-1\right)\right)}^{\prime }$ and has the mean ${\left(6.25;12.5\right)}^{\prime }$ and the covariance $\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left\{{2.5}^2;5^2\right\}.$

The continuous movement model is sampled in increments of $\delta =0.0001$ h (hour). The movement begins at time $t=-T=-50$ and ends at time $t=1000$, i.e., the movement lasts moment $T^{end}=1000\delta =0.1$ h or 6 min.

We denote $v_x\left(t\right)$ and $v_y\left(t\right)$ AUV velocity values by directions $x$ and $y$, respectively. At first, we will assume that the velocity on each trajectory is constant, i.e., $v_x\left(t\right)=v_x\left(-T-1\right)$ and $v_y\left(t\right)=v_y\left(-T-1\right)$. Initial Gaussian vector ${\left(v_x\left(-T-1\right),v_y\left(-T-1\right)\right)}^{\prime }$ has the known mean, $\mathrm{E}\left\{v_x(-T-1)\right\}=25$ kph and $\mathrm{E}\left\{v_y(-T-1)\right\}=12.5$ kph and the covariance $\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left\{{12.5}^2;{6.25}^2\right\}$. Accordingly, we assume that the values $v_x\left(-T-1\right)$ and $v_y\left(-T-1\right)$ are unknown and subject to identification, i.e., $\mu \in {\mathbb{R}}^2$, ${\mu }_1=v_x\left(-T-1\right)$, ${\mu }_2=v_y\left(-T-1\right)$.

Then, turning to model (3), we will assume that the velocity ${\left(v_x\left(t\right),v_y\left(t\right)\right)}^{\prime }$ keeps constant values at time intervals between process ${\lambda }_{u}P\left(u\right)$ jumps where ${\lambda }_u=\lambda ={\displaystyle \frac{3}{T^{end}}},$ $u\in \left[0,T^{end}\right],$ which corresponds to three velocity changes on average over the tracking period, or an average velocity constant time of 2 min. Next jump ${\lambda }_{u}P\left(u\right)$ in interval $\left(\delta \left(t-1\right),\delta t\right)$ means that in time point $t,$ the previous constant velocity value ${\left(v_x\left(t-1\right),v_y\left(t-1\right)\right)}^{\prime }$ replaces the new value, which has the same distribution as ${\left(v_x\left(-T-1\right),v_y\left(-T-1\right)\right)}^{\prime }$.

The AUV is affected by disturbances $w_x\left(t\right)$ in the $x$ axis direction and $w_y\left(t\right)$ in the $y$ axis direction. They form a real velocity deviating from the nominal velocity ${\left(v_x\left(t\right),v_y\left(t\right)\right)}^{\prime }$ to a random value of the order of 25 kph. So, the position vector has the form:

$$\begin{array}{c}x\left(t\right)=x\left(t-1\right)+\delta \left(v_x\left(t\right)+w_x\left(t\right)\right),\\ y\left(t\right)=y\left(t-1\right)+\delta \left(v_y\left(t\right)+w_y\left(t\right)\right).\end{array}$$

(11)

Describing chaotic velocity changes, vector ${\left(w_x\left(t\right),w_y\left(t\right)\right)}^{\prime }$ is assumed to be Gaussian with the mean ${\left(0,0\right)}^{\prime }$ and the covariance $\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left\{{25}^2;{25}^2\right\}$.

The measuring complex consists of the same type of observers located at two points on the $xy$ plane: the first has coordinates $\left(0,l_y\right)$, $l_y=6.25$ km, the second has coordinates $\left(l_x,0\right)$, $l_x=12.5$ km. Assuming that acoustic sonars are used for observation, and the velocity of sound in water is equal to $v_s=5400$ kph, we can determine the maximum possible observation delay. We have already indicated above that $T=50$, i.e., the maximum value by which observations can be delayed is $18$ s. Suppose the observers are stationary acoustic beacons installed in advance. Typical measurements performed by such devices are bearing angles [43]. In order not to clutter up our model, we will measure the distance to the AUV since the receiving measurements delay depends on the distance. It is clear that with two beacons, distance determining is a simple geometric problem solution. In the experiment, we immediately assume that each observer measures the distance to the target and the bearing angle.

So, we get the measurements $d_1\left(t\right), d_2\left(t\right)$ (distances measured by the first and second observers) and $c_1\left(t\right), c_2\left(t\right)$ (direction cosines measured by the first and second observer). The Gaussian vector of observation errors ${\left(v_{d_1}\left(t\right),v_{d_2}\left(t\right),v_{c_1}\left(t\right),v_{c_2}\left(t\right)\right)}^{\prime }$ has the mean ${\left(0,0,0,0\right)}^{\prime }$ and the covariance $\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left\{{0.001}^2;{0.005}^2;{0.001}^2;{0.005}^2\right\}.$

The observers measure the AUV position characteristics at time points $t-{\left({\tau }_t\right)}_1$ and $t-{\left({\tau }_t\right)}_2$, i.e., with a delay from the current position ${\left(x\left(t\right),y\left(t\right)\right)}^{\prime }$ by the value of ${\left({\tau }_t\right)}_1$ (first) and ${\left({\tau }_t\right)}_2$ (second), that take the following values:

$$\begin{array}{c}{\left({\tau }_t\right)}_1=\mathrm{m}\mathrm{i}\mathrm{n}\left\{T,\left[\sqrt{{\left(x\left(t\right)\right)}^2+{\left(y\left(t\right)-l_y\right)}^2}/\left({\delta v}_s\right)\right]\right\},\\ {\left({\tau }_t\right)}_2=\mathrm{m}\mathrm{i}\mathrm{n}\left\{T,\left[\sqrt{{\left(x\left(t\right)-l_x\right)}^2+{\left(y\left(t\right)\right)}^2}/\left({\delta v}_s\right)\right]\right\}.\end{array}$$

(12)

In (12), the notation $\left[x\right]$ is used for the floor function of $x$ and the potential possibility of the observed object going at a distance for which the delay value becomes greater than the specified maximum $T$ is taken into account.

Thus, the observation vector $y_t={\left(d_1\left(t\right),c_1\left(t\right),d_2\left(t\right),c_2\left(t\right)\right)}^{\prime }\in {\mathbb{R}}^4$ is as follows:

$$\begin{array}{c}{d}_1\left(t\right)=\sqrt{{\left(x\left(t-{\left({\tau }_t\right)}_1\right)\right)}^2+{\left(y\left(t-{\left({\tau }_t\right)}_1\right)-l_y\right)}^2}+v_{d_1}\left(t\right),\\ \begin{array}{c}\begin{array}{c}{c}_1\left(t\right)={\displaystyle \raisebox{1ex}{$\left(y\left(t-{\left({\tau }_t\right)}_1\right)-l_y\right)$}\!\left/ \!\raisebox{-1ex}{$\sqrt{{\left(x\left(t-{\left({\tau }_t\right)}_1\right)\right)}^2+{\left(y\left(t-{\left({\tau }_t\right)}_1\right)-l_y\right)}^2}$}\right.}+v_{c_1}\left(t\right),\\ d_2\left(t\right)=\sqrt{{\left(x\left(t-{\left({\tau }_t\right)}_2\right)-l_x\right)}^2+{\left(y\left(t-{\left({\tau }_t\right)}_2\right)\right)}^2}+v_{d_2}\left(t\right),\end{array}\\ c_2\left(t\right)={\displaystyle \raisebox{1ex}{$\left(x\left(t-{\left({\tau }_t\right)}_2\right)-l_x\right)$}\!\left/ \!\raisebox{-1ex}{$\sqrt{{\left(x\left(t-{\left({\tau }_t\right)}_2\right)-l_x\right)}^2+{\left(y\left(t-{\left({\tau }_t\right)}_2\right)\right)}^2}$}\right.}+v_{c_2}\left(t\right).\end{array}\end{array}$$

(13)

To bring these designations into line with (1) the vector ${\tau }_t\in {ℝ}^4$ components must be expressed using ${\left({\tau }_t\right)}_1$ and ${\left({\tau }_t\right)}_2$ from (12) as ${{\tau }_1}_t={{\tau }_2}_t={\left({\tau }_t\right)}_1$ and ${{\tau }_3}_t={{\tau }_4}_t={\left({\tau }_t\right)}_2$.

5.2. CMNF Synthesis

The filter (8) synthesis consists of the structural functions ${\xi }_t$ and ${\zeta }_t$ and in the coefficients (9) approximate calculation by the Monte Carlo method. The proposed in [3] structure fully meets the problem under consideration. It remains only to clarify the aspects related to the identified velocity.

Thus, the basic prediction function ${\xi }_t={\xi }_t\left(X,M\right)$ has the form due to the dynamics (11), i.e.,

$${\xi }_t\left({\widehat{x}}_{t-1},{\widehat{\mu }}_{t-1}\right)={\phi }_t\left({\widehat{x}}_{t-1},{\widehat{\mu }}_{t-1},\mathrm{E}\left\{w_t\right\}\right)=\left(\begin{array}{c}\tilde{x}\left(t\right)\\ \tilde{y}\left(t\right)\end{array}\right)=\left(\begin{array}{c}\widehat{x}\left(t-1\right)+\delta {\widehat{v}}_x\left(t-1\right)\\ \widehat{y}\left(t-1\right)+\delta {\widehat{v}}_y\left(t-1\right)\end{array}\right),$$

where ${\widehat{x}}_t={\left(\widehat{x}\left(t\right),\widehat{y}\left(t\right)\right)}^{\prime }$, ${\widehat{\mu }}_t={\left({\widehat{v}}_x\left(t\right),{\widehat{v}}_y\left(t\right)\right)}^{\prime }$.

Further, to make a correction, we need to estimate (more exactly to predict) values for delays (12), which are naturally set as follows:

$$\begin{array}{c}{\left({\tilde{\tau }}_t\right)}_1=\min \left\{T,\left[\sqrt{{\left(\tilde{x}\left(t\right)\right)}^2+{\left(\tilde{y}\left(t\right)-l_y\right)}^2}/\left(\delta v_s\right)\right]\right\},\\ {\left({\tilde{\tau }}_t\right)}_2=\min \left\{T,\left[\sqrt{{\left(\tilde{x}\left(t\right)-l_x\right)}^2+{\left(\tilde{y}\left(t\right)\right)}^2}/\left(\delta v_s\right)\right]\right\},\end{array}$$

(14)

where ${\tilde{x}}_t={\left(\tilde{x}\left(t\right),\tilde{y}\left(t\right)\right)}^{\prime }$.

We define the basic correction function ${\zeta }_t={\zeta }_t\left(X,Y\right)$ as a result of recalculating the distances and cosines into Cartesian coordinates. The obtained approximate coordinates correspond to the AUV’s position with the delays ${\left({\tau }_t\right)}_1,$ ${\left({\tau }_t\right)}_2,$ so it needs to be shifted. For the shifting, we use the delay estimates ${\left({\tilde{\tau }}_t\right)}_1$ and ${\left({\tilde{\tau }}_t\right)}_2$ out of (14) and the current velocity identification estimates ${\widehat{v}}_x\left(t-1\right), {\widehat{v}}_y\left(t-1\right)$. Finally, the correction itself is the difference between the predicted position and the calculated approximate coordinates:

$${\zeta }_t\left({\tilde{x}}_t,y_t\right)=\left(\begin{array}{c}\tilde{x}\left(t\right)-\left(\sqrt{{\left(d_1\left(t\right)\right)}^2-{\left(c_1\left(t\right)d_1\left(t\right)\right)}^2}+{\left({\tilde{\tau }}_t\right)}_1{\widehat{v}}_x\left(t-1\right)\right)\\ \begin{array}{c}\tilde{y}\left(t\right)-\left(c_1\left(t\right)d_1\left(t\right)+l_y+{\left({\tilde{\tau }}_t\right)}_1{\widehat{v}}_y\left(t-1\right)\right)\\ \tilde{x}\left(t\right)-\left(c_2\left(t\right)d_2\left(t\right)+l_x+{\left({\tilde{\tau }}_t\right)}_2{\widehat{v}}_x\left(t-1\right)\right)\end{array}\\ \tilde{y}\left(t\right)-\left(\sqrt{{\left(d_2\left(t\right)\right)}^2-{\left(c_2\left(t\right)d_2\left(t\right)\right)}^2}+{\left({\tilde{\tau }}_t\right)}_2{\widehat{v}}_y\left(t-1\right)\right)\end{array}\right).$$

(15)

5.3. Computer Simulations

In all performed computational experiments, we simulated two independent samples of $N=100000$ random trajectories of the observation system (11) and (13) (in two versions (1) or (3)), CMNF (8) with basic correction (15) position estimate ${\widehat{x}}_t={\left(\widehat{x}\left(t\right),\widehat{y}\left(t\right)\right)}^{\prime }$, and velocity estimate ${\widehat{v}}_t={\left({\widehat{v}}_x\left(t\right),{\widehat{v}}_y\left(t\right)\right)}^{\prime }.$ We used the first sample to calculate the filter coefficients (9) and their theoretical accuracy (10). This sample can be called a training sample. The second sample allows us to calculate the CMNF estimate ${\widehat{x}}_t$ and ${\widehat{v}}_t$ actual quality. When the model is system (1), then ${\widehat{v}}_t={\left({\widehat{v}}_x\left(-T-1\right),{\widehat{v}}_y\left(-T-1\right)\right)}^{\prime }$ and the task is unknown constant velocity identification, more precisely its systematic component. When the model is system (3), then ${\widehat{v}}_t$ is an estimate of the current velocity, which remains constant over a certain time interval. We need to see how accurate the velocity estimate is at these intervals and how the filter responds to sudden changes in velocity.

The estimate error mean square standard deviations denoted by ${\sigma }_{\widehat{x}}\left(t\right),$ ${\sigma }_{\widehat{y}}\left(t\right),$ ${\sigma }_{{\widehat{v}}_x}\left(t\right),$ ${\sigma }_{{\widehat{v}}_x}\left(t\right),$ determine the estimation accuracy of the AUV’s position and velocity. This is practical accuracy. The values $\sqrt{{\left({\widehat{K}}_t^x\right)}_{i,i}}, i=1, 2,$ represent the theoretical accuracy of the CMNF position estimate. The values $\sqrt{{\left({\widehat{K}}_t^{\mu }\right)}_{i,i}}, i=1, 2,$ represent the theoretical accuracy of the CMNF velocity estimate. We compare these values with the actual accuracy ${\sigma }_{\widehat{x}}\left(t\right),$ ${\sigma }_{\widehat{y}}\left(t\right),$ ${\sigma }_{{\widehat{v}}_x}\left(t\right),$ ${\sigma }_{{\widehat{v}}_x}\left(t\right),$ and if they coincide (sufficiently close), we conclude that the filter synthesis is successful, i.e., there is the sufficient statistical volume of the training sample.

Our fundamental interest, of course, is models with the specified time delay value $T=50$ and two options (1) or (3) for movement parameter identification, i.e., with velocity components $v_x$ and $v_y$ standard deviations $12.5$ kph and $6.25$ kph, respectively. But in order to understand what effect the time delay and unknown movement parameter factors have on the result, we also performed calculations without time delay modeling, i.e., for $T=0,$ and without velocity identification, i.e., for zero variances ${\left(v_x\left(-T-1\right),v_y\left(-T-1\right)\right)}^{\prime }$ and assuming that they are known.

Figure 1, Figure 2 and Figure 3 begin to illustrate the obtained results. These figures represent several randomly selected trajectories from the available $100000$ trajectories of the second sample and the diversity of the AUV movement nature. The characteristics of the estimates’ accuracy are discussed below.

Figure 1 shows several typical AUV trajectories ${\left(x\left(t\right),y\left(t\right)\right)}^{\prime }$ and the corresponding CMNF position estimates ${\left(\widehat{x}\left(t\right),\widehat{y}\left(t\right)\right)}^{\prime }.$ The graphs at the top are trajectories for the switching velocity model (3), and the bottom graphs are trajectories for the constant and unknown velocity model (1). We see how, in the lower graphs, the trajectories are quite diverse and chaotic. However, they still show clearly defined directions that remain constant throughout the entire observation period. Compared with the trajectories in the upper graphs, these show only the general trend of AUV movement away from observers as provided by the model. At the same time, the nature and movement direction do not demonstrate constancy but change at different moments, in different directions, with different magnitudes.

Figure 2 complements the trajectories examples. It presents the velocity estimate ${\left(v_x\left(t\right),v_y\left(t\right)\right)}^{\prime }$ trajectories. As in the previous figure, the graphs above correspond to model (3), and below to model (1). In all the graphs presented, we can see that the filter is trying to identify an unknown velocity and does it quite successfully. But if we look at the graphs below, we can see a typical example of parametric identification where the estimates tend to converge to a specific value. On the other hand, in the graphs above, the filter does not have enough time and observations to indicate convergence. Instead, it can only demonstrate “intentions” toward convergence. It is all the more interesting to see how such a difference affects the solution of the main task of AUV tracking.

Figure 3 of this group illustrates several time delay value trajectories, ${\left({\tau }_t\right)}_1$ and ${\left({\tau }_t\right)}_2$ from (12). Here, we can assess the adequacy of the assumption $T=50$ and the variety of forms.

Figure 1, Figure 2 and Figure 3 demonstrate a qualitative picture of the estimation problem solutions. To determine the estimates’ objective quality and compare them, we need to calculate a quadratic criterion. According to (10), filter coefficients (9) determine the theoretical accuracy $\sqrt{{\left({\widehat{K}}_t^x\right)}_{i,i}}, i=1, 2,$ and $\sqrt{{\left({\widehat{K}}_t^{\mu }\right)}_{i,i}}, i=1, 2,$ so we find out these values when we synthesize the filter in the first sample.

The second sample purpose is only to verify the synthesis effectiveness, for which we need to calculate the actual accuracy ${\sigma }_{\widehat{x}}\left(t\right),$ ${\sigma }_{\widehat{y}}\left(t\right),$ ${\sigma }_{{\widehat{v}}_x}\left(t\right),$ ${\sigma }_{{\widehat{v}}_x}\left(t\right).$ For ease of interpretation, the measurement unit for the quality of the AUV’s position estimate quality is meters, and kph is left for velocity. Figure 4 shows a graph of the position estimate accuracy, and Figure 5 shows a graph of the velocity estimate accuracy. In each pair of graphs, first of all, we compare the calculation results for model (3) (left) and for model (1) (right).

All these visual and quantitative figures illustrate the real difference between models with and without velocity jumps. Of course, as seen in both Figure 4 and Figure 5, the differences between the theoretical and actual accuracy of all CMNF estimates are very small. To give a completely quality objective assessment of the performed estimates,

Table 1. Filters quality comparison in the target-tracking task.

Model	${\mathit{\sigma }}_{\widehat{\mathit{x}}}^{\left(\mathit{r}\mathit{e}\mathit{a}\mathit{l}\right)}$	${\mathit{\sigma }}_{\widehat{\mathit{x}}}^{\left(\mathit{t}\mathit{e}\mathit{o}\mathit{r}\right)}$	${\mathit{\sigma }}_{\widehat{\mathit{y}}}^{\left(\mathit{r}\mathit{e}\mathit{a}\mathit{l}\right)}$	${\mathit{\sigma }}_{\widehat{\mathit{y}}}^{\left(\mathit{t}\mathit{e}\mathit{o}\mathit{r}\right)}$
$\lambda ={\displaystyle \frac{3}{1000}}$, $T=50$	34.89	34.82	22.01	22.96
$\lambda =0$, $T=50,$ $\mu \ne \mathrm{const}$	22.74	22.69	18.20	18.17
$\lambda =0$, $T=50,$ $\mu =\mathrm{const}$ *	19.20	19.01	16.72	16.58
$\lambda ={\displaystyle \frac{3}{1000}}$, $T=0$	18.09	18.10	13.80	13.79
$\lambda =0$, $T=0,$ $\mu \ne \mathrm{const}$	14.97	14.94	12.61	12.59
$\lambda =0$, $T=0,$ $\mu =\mathrm{const}$ *	14.15	13.94	12.26	12.05

* Constant and known velocity, no identification task.

and

Table 2. Filters quality comparison in the identification task.

Model	${\mathit{\sigma }}_{{\widehat{\mathit{v}}}_{\mathit{x}}}^{\left(\mathit{r}\mathit{e}\mathit{a}\mathit{l}\right)}$	${\mathit{\sigma }}_{{\widehat{\mathit{v}}}_{\mathit{y}}}^{\left(\mathit{t}\mathit{e}\mathit{o}\mathit{r}\right)}$	${\mathit{\sigma }}_{{\widehat{\mathit{v}}}_{\mathit{x}}}^{\left(\mathit{r}\mathit{e}\mathit{a}\mathit{l}\right)}$	${\mathit{\sigma }}_{{\widehat{\mathit{v}}}_{\mathit{y}}}^{\left(\mathit{t}\mathit{e}\mathit{o}\mathit{r}\right)}$
$\lambda ={\displaystyle \frac{3}{1000}}$, $T=50$	6.79	6.78	3.79	3.76
$\lambda =0$, $T=50,$ $\mu \ne \mathrm{const}$	0.85	0.84	0.85	0.83
$\lambda ={\displaystyle \frac{3}{1000}}$, $T=0$	5.59	5.59	3.50	3.48
$\lambda =0$, $T=0,$ $\mu \ne \mathrm{const}$	0.85	0.83	0.83	0.82

below show the accuracy numerical characteristics. Namely, the quality analysis of ${\widehat{x}}_t$ gives the following values:

$$\begin{array}{c}{\sigma }_{\widehat{x}}^{\left(real\right)}={\displaystyle \frac{1}{1000}}{\displaystyle {\displaystyle \sum }_{t=1}^{1000}}{\sigma }_{\widehat{x}}\left(t\right), {\sigma }_{\widehat{y}}^{\left(real\right)}={\displaystyle \frac{1}{1000}}{\displaystyle {\displaystyle \sum }_{t=1}^{1000}}{\sigma }_{\widehat{y}}\left(t\right),\\ \end{array}$$

$$\begin{array}{c}{\sigma }_{\widehat{x}}^{\left(teor\right)}={\displaystyle \frac{1}{1000}}{\displaystyle {\displaystyle \sum }_{t=1}^{1000}}{\left({\widehat{K}}_t^x\right)}_{11}^{1/2}, {\sigma }_{\widehat{y}}^{\left(teor\right)}={\displaystyle \frac{1}{1000}}{\displaystyle {\displaystyle \sum }_{t=1}^{1000}}{\left({\widehat{K}}_t^x\right)}_{22}^{1/2}.\\ \end{array}$$

The values characterizing the quality of identification of the estimates ${\widehat{\mu }}_t$ have the following form:

$$\begin{array}{c}{\sigma }_{{\widehat{v}}_x}^{\left(real\right)}={\displaystyle \frac{1}{100}}{\displaystyle {\displaystyle \sum }_{t=901}^{1000}}{\sigma }_{{\widehat{v}}_x}\left(t\right), {\sigma }_{{\widehat{v}}_x}^{\left(real\right)}={\displaystyle \frac{1}{100}}{\displaystyle {\displaystyle \sum }_{t=901}^{1000}}{\sigma }_{{\widehat{v}}_x}\left(t\right),\\ \end{array}$$

$$\begin{array}{c}{\sigma }_{{\widehat{v}}_x}^{\left(teor\right)}={\displaystyle \frac{1}{100}}{\displaystyle {\displaystyle \sum }_{t=901}^{1000}}{\left({\widehat{K}}_t^{\mu }\right)}_{11}^{1/2}, {\sigma }_{{\widehat{v}}_y}^{\left(teor\right)}={\displaystyle \frac{1}{100}}{\displaystyle {\displaystyle \sum }_{t=901}^{1000}}{\left({\widehat{K}}_t^{\mu }\right)}_{22}^{1/2}.\\ \end{array}$$

To make the comparison as informative as possible, we have supplemented the calculations for the main models (1) and (3) with a few more simplified ones. The most complete model (3) is denoted in the tables as $\lambda ={\displaystyle \frac{3}{1000}},$ $T=50,$ and its simpler version without velocity jumps is denoted as $\lambda =0,$ $T=50,$ $\mu \ne \mathrm{const}.$ The example for the model with known velocity $v_x=12.5$ kph, $v_y=6.25$ kph complements the options for comparison, and its notation is $\lambda =0,$ $T=50,$ $\mu =\mathrm{const}.$ Table 2 does not contain this model because there is no identification in it. And for all three models, we repeated the calculation under the assumption $T=0$. Note that the calculations for the third ($\lambda =0,$ $T=50,$ $\mu =\mathrm{const}$) and sixth ($\lambda =0,$ $T=0,$ $\mu =\mathrm{const}$) models correspond to the experiment presented in [3].

Once again, note that the quantities dimensions in Table 1 are meters. It should be noted that the parameters in models (11)–(13), although they do not conflict with common sense, may have “extreme” values in some sense. This applies both to the “fast” and chaotic AUV movement, as well as to observers who are not very accurate and are located far away. This choice is meaningful because the computational experiment task was to compare models and analyze certain factors that influence the estimates’ accuracy. For this reason, such “large” parameters were chosen. Some more adequate numerical characteristics can be found in [44].

6. Conclusions

As a final discussion, we propose some conclusions highlighted by the calculation results in Table 1 and Table 2:

• With a relatively small $100000$ training sample size, the difference between the accuracy of the theoretical and actual estimates in all models is negligible (one of the challenges of using CMNF is the inability to calculate analytically the filter coefficients using formulas (9) and the computer modeling; as with any machine learning method, which, of course, is Monte Carlo method, computer simulation has applicability limitations and requires an effectiveness confirmation. The model under consideration has experimentally shown that replacing the exact coefficients of conditionally optimal filtration (9) with their computer-simulated approximations has practically no effect on either the AUV’s position estimate accuracy or movement parameter identification accuracy. It can be reasonably assumed that this is due to the exceptionally successful filter structure which, in turn, is due to the understandable and realistic nature of the problem being solved).
• The first and most significant factor affecting the accuracy of the AUV’s position estimates is the consistency of its velocity: if the velocity is constant for a long time (model (1)), then the filter has time to identify it, so there is no fundamental difference in the effectiveness of target tracking at a known velocity and at an unknown velocity; abrupt velocity jumps $\left(\lambda \ne 0\right)$ have a great impact, impairing the tracking effectiveness.
• The second major factor that affects the accuracy of AUV’s coordinate estimates is the time delay in observations, due to the deterioration in tracking accuracy in models with $T\ne 0$, which is indeed significant. It was expected that the effect of these delays would be aggravated; this is because the model becomes more complex by assuming the presence of unknown identifiable parameters.
• In general, in all models, the identification result can be considered positive, especially, in model (1), in 6 min of measuring only the bearing, the velocity can be identified with an accuracy of about 0.75–1.5 kph; in model (3), of course, it is incorrect to talk about identification because the velocity here has become one from the state variables. It is important that the obtained velocity estimates still allow us to solve the main problem of AUV tracking quite successfully.
• If it is difficult to increase the bearing measurement accuracy, and a higher tracking quality is required, then, perhaps, the only way to improve it is to measure velocity.