Tracking a Maneuvering Object by Indirect Observations with Random Delays

: A mathematical model for the target tracking problem is proposed. The model makes it possible to describe conditions when the time for an observer to receive the results of indirect observations of a moving object depends not only on the state of the observation environment but also on the state of the object itself. The source of such a model is the observation process, by stationary means, of an autonomous underwater vehicle, in which the time for obtaining up-to-date data depends on the unknown distance between the object and the observer. As part of the study of the problem, the equations of the optimal Bayesian ﬁlter are obtained. But this ﬁlter is not possible to implement. For practical purposes, it is proposed to use the conditionally minimax nonlinear ﬁlter, which has shown promising results in other complex tracking models. The conditions for the ﬁlter’s evaluation and its accuracy characteristics are given. A large-scale numerical experiment illustrating the ﬁlter’s operation and the observation system’s features with random delays are described.


Introduction
The problem of tracking a maneuvering object using a model with random factors is a typical application of stochastic filtering theory and methods [1]. In addition to navigation, these methods are widely used in other areas, sometimes entirely unexpectedly, such as medicine, industry, economics, and telecommunications [2][3][4][5][6][7]. The model and problem statement discussed in the article are inspired by the popular modern applied field of autonomous underwater vehicles (AUVs) [8]. Many works are devoted to AUV applications, most of which are related to control tasks. The special issue of [9] sufficiently illustrates the scope and magnitude of the research interest in this area. Moreover, let us pay a ention to the fact that researchers use many new methods as tools, primarily due to the development of machine learning (a typical example is work [10], and there are many such works [11]). At the same time, the conventional tasks of target tracking remain in demand; one can find many new and exciting things in them [12][13][14][15][16][17][18][19][20]. Here it is necessary to note the tasks of tracking multiple targets [21,22], the solution of which, as well as conventional statements, is provided by stochastic filtering methods. The Kalman filtering model, algorithm [23], and many different suboptimal filters have become widely used. The simplest algorithms of this kind, starting with the extended Kalman filter [24] and modern developments of the odorless or cubature Kalman filters [25][26][27], can be effective in individual applications. Such ideas have found many applications and options for development/improvement. An excellent a empt to systematize the results in this area was made in [28]. To date, many more works have appeared devoted to the research and improvement in such suboptimal filters [29][30][31][32], including those applied to target tracking tasks of AUVs [33]. However, no fundamentally new ideas have been proposed, and the initially existing shortcomings remain. The essence of these disadvantages is that none of these filters guarantees accurate characteristics for at least any broad class of models. That is, the estimates do not guarantee either the absence of bias or the limitation of the error, at least by the variance of the estimated state. This means there is always a danger that the filter will be unstable in another model and the estimate will diverge. The only guarantee is a practical experiment with a specific model and a fixed set of parameters.
It should be noted that the problem of delayed observations is well known in the field of application under discussion. Such delays are due to a prominent characteristic of acoustic meters: the dependence of the acoustic wave propagation velocity on temperature, salinity, and water pressure [34]. This a racted a ention in the above-mentioned study [20], and a solution was found by combining measurement data from acoustic sensors with information from other sensors in the onboard inertial navigation system. A more typical way to deal with delayed observations is to estimate the delay time and consider it by adjusting the usual filter [35]. Interestingly, the delay in observations is explained by nature not only for underwater observations but also in some chemical processes [36].
It should be understood that models with a deterministic delay time have limited use because, in reality, the delay time is random and significantly depends on the observation conditions, and most importantly, it changes during tracking. Moreover, these changes can be significant. In the model proposed in this article, the time delay of observations is described by a random process, which is a known function of the state of the observation system. Such a description of delays makes the model as adequate as possible when describing the results of the acoustic sensors.
At the same time, formally, the model of the observation system remains a Markov process model with discrete time. Hence, the classical Bayesian filtering procedure can be applied to it [37], i.e., obtaining recurrent relations for the a posteriori probability density is possible. The next section of the article shows these relations. However, this formal success turned out to be relevant only to justify the impossibility of practical implementation of either the optimal filter itself or its suboptimal simplifications. The reason for this is that the formal reduction of the equations of the observation system to the traditional form enormously increases the dimension.
An alternative to suboptimal filters is the conditionally minimax nonlinear filtering (CMNF) [38]. Initially and in development [39], this method focused on applications in aircraft tracking tasks, which remain relevant for modern tasks [40]. However, the same approach has shown effectiveness in specific tasks for AUVs. So, in [19], an already modified CMNF for tracking underwater vehicles was shown, but with an emphasis on the efficiency of using Doppler sensors. It turned out that the CMNF algorithm does not require refinement for the proposed model with random delay. It is only necessary to clarify the conditions of existence. Section 3 of the article contains the statement required.
It turned out to be more challenging to adapt the typical filter structure to the task of AUV tracking. A standard filter structure using a prediction via a system state model and a correction in the form of a residual did not show good results when processing the acoustic observations. Nevertheless, the flexibility of the CMNF method made it possible to adjust the filter structure, consider the physical properties of observers, and obtain a very effective estimate as a result. Within the framework of the numerical experiment presented in Section 4 of the article, we show the process of physically justified adjustments to the filter structure. Section 5 contains concluding remarks, including considerations on options for further development of the proposed method.

Observation System Model and Optimal Filtering
The model uses discrete time. It is assumed that the filtering process, i.e., the calculation of the first state estimate, begins at the moment = 0. The maximum possible delay is known and is equal to ; therefore, the state of the system begins to form at the moment = − . Thus, = − , − + 1, … ,0, 1, …, and the initial state is set at the moment − − 1. Three variants of recording the equations of the observation system are considered. The first one has the most general form: The second version of the recording is a discrete analog of the diffusion process: The la er version is the most particular, but the most frequently used. This model is used by many suboptimal filters. It is also the most convenient for optimal Bayesian filtering. The system with state-independent noise is wri en as follows: In (1)-(3) ∈ ℝ represents the system state vector; ∈ ℝ is the discrete white noise modeling perturbations of the system; ∈ ℝ is the initial state vector; ∈ ℝ is the vector of indirect observations; ∈ ℝ is the discrete white noise modeling measurement errors; vectors , , are assumed to be independent in aggregate. The model of observations delays, , is given by a vector of the same dimension as the observations themselves. The elements of the vector, ∈ ℝ , set the time delays for the corresponding components of the vector . . This notation makes accurate (1)-(3) and further relations.
The problem of estimating the state from observations , = 0, 1, … , , is considered, and the evaluation accuracy criterion is the mean square, i.e., E{‖ − ‖ }, E{ } is the mathematical expectation of the vector , and ‖ ‖ is the usual Euclidean norm of the vector .
Thus, the random observation delay ( ) is the only difference between the considered problem statement and the traditional one of filtering the state of the observation system in discrete time. This difference does not exist formally because the following shows how to bring systems (1)-(3) into the traditional recording form. Therefore, the optimal solution is a conditional mathematical expectation of relative to observations , = 0, 1, … , . So, it is enough to know the a posteriori probability density of relative to , = 0, 1, … , . Next, the a posteriori density for system (3) and = 0, 1, …, will be wri en out in the form of recurrent Bayesian relations [37] under conditions when the necessary probability densities exist, i.e., under certain restrictions on the continuity of nonlinear functions and perturbations in (3). We avoid unnecessarily cumbersome relations with no practical value by limiting the discussion of the optimal filter to the system (3). For formal conclusions, system (3) is sufficient.
The non-Markov observation system (3) can be wri en as a Markov process with an extended state vector. Formally, this extended process will have the same form (3), but with = 0. The extended state vector is further denoted by ∈ ℝ ( ) and is a composite vector, consisting of all the states of the system from time − to the current moment , i.e., = ( , … , , ) . The notation ′ is used here for the transpose operation. The equation for has the following form: where ( ) denotes a sub-vector of vector with elements from -th to -th. Denoting such a transformation of the state vector by Φ , and the corresponding additive noise by = (0 , … , 0 , ) , we obtain the equation of state in the following form: To obtain the relation for the observer, we define the matrix function Ψ ( ) ∈ ℝ ( )× as follows: Thus, the rows of this matrix contain all possible observations at time due to delays without noise. To describe the delays themselves, we also introduce the matrix function Θ ( ) ∈ ℝ ×( ) , which takes values by the following rule: Thus, the element ( , ) of the matrix Θ ( ) takes the value 1 if the observation ( ) corresponds to the delay ( ) = − + 1.
These two notations allow us to write the equation of observations in the following form: Thus, (4) and (5) are the canonical form of a Markov observation system with discrete time and additive independent noises.
Further, the values that the function Θ ( ) can take are denoted by [Θ ] , = 1, … , ( + 1). The sequence number can be set by pu ing where , is the row number of the -th column of the matrix [Θ ] , in which there is 1. The designations could be more convenient, but they are needed to represent the optimal Bayesian filter. To this end, we introduce the following rules for the designation of probabilistic characteristics. For random vectors ∈ ℝ , ∈ ℝ via ( , ) ∈ ℝ denote the composite vector using the notation ′ (transpose).
The probability densities used below will be denoted as follows: ( ) is the marginal density of , , ( , ) is the density of the composite vector ( , ) , and | ( | ) is the conditional density of x relative to y under the additional assumption ( ) > 0.
Additionally, we list all the symbols used:  Symbol is used to denote random vectors in the notations of the original observation system (3);  Symbol is used to denote arguments of probability densities corresponding to ;  Symbol is used to denote random vectors in the extended state notation in the observation system (4) and (5);  Symbol is used to denote arguments of probability densities corresponding to ;  Notations and are used for the vector of all observations up to and including the moment , i.e., = ( , … , ) , and for the corresponding probability density argument.
For the conditional distribution of the state of the observation systems (4) and (5) with respect to observations , = 0, … , , it is possible to formally write the recurrent Bayesian relations for the a posteriori probability density [37]: This expression can be used only if all probability densities exist in it. For systems (4) and (5), this is not the case for the transition density | ( | ). In addition, it is necessary to take into account the explicit form of density | ( | ). The transformations performed in the following theorem allow us to use instead of the transition density | ( | ) of the expanded state and the transition density | ( | ) of the initial state. Accordingly, we can instead compute instead ∫• , which is calculated according to the space ℝ ( ) and the integral ∫• according to the space ℝ . In addition, an explicit form of density, | ( | ), is written. Then, for the a posteriori probability density = ( | ) of the extended state = ( , … , , ) , = 0,1, … , of the system given by the relation (4), with respect to observations = ( , … , ) , = 0,1, … , and the estimates of the optimal Bayesian filter * = { | } of the system (3) state , the following recurrent equalities are performed: with an initial condition and the numbering rule for = 1, … , , defined in (6).
In the recurrent equalities of the theorem for a posteriori densities ( | ), arguments are omi ed, and the recording of the sum is simplified. Note that there are no formal reasons preventing computer implementation (7). However, there needs to be more realism in performing such a computer calculation due to the increase in dimension, the scale of which determines the value of . In the numerical experiment discussed later in the article, = 50. So, even with small = 2 and = 4, the integrals in (7) will be in the space ℝ , and the summands will be 404. With such parameters, calculations according to (7) cannot be performed in practice. This confirmed the need for realistic prospects for using classical filtering algorithms in the problem under consideration that both the extended observation system (4) and (5), and Theorem 1 were needed. is denoted as ∑ , and it is taken into account that the coefficient ( )/ ( ) is a normalizing factor. Conditions of boundedness of the second-order moments and linear growth of the functions of the system are sufficient for the existence of second-order moments of and : E{‖ ‖ } + E{‖ ‖ } < ∞, and optimal in mean square sense estimate E{ | }. □

Conditionally Minimax Nonlinear Filter
The main statements of the conditionally minimax nonlinear filtering (CMNF) theory are in the works [33,34]. Further, they are formulated briefly, considering the features of the surveillance system under consideration.
CMNF estimate of the state by the observations has the form of predictioncorrection: = + ∆ . A basic predictive function = ( ) is used to calculate the prediction . As a simple example for system (1) where ℱ denotes the distribution of the vector = ( , ) . It is assumed that ℱ ∈ ℱ( , ) is a class of all probability distributions with mean and covariance . Thus, the basic prediction and correction are refined best (in the sense of proximity in the mean square to the estimated state) under the assumption that only the second-order moment characteristics of the estimated and observed variables are known.
We emphasize that the options for choosing the structural functions and are not limited by anything. In practice, such a choice may be an independent task, the solution of which will allow taking into account the nature of a particular dynamic system. Such an example is given in the article's next section devoted to experimental research.
The worst distribution in problems (8) is the normal one, and the corresponding best estimate in the mean square sense is determined by the normal correlation theorem [41]. Thus, solutions (8), i.e., optimal in the minimax sense functions Ξ ( ) and Ζ ( ) are linear as follows: In (10), the notation cov( , ) is the covariance of and , is the Moore-Penrose pseudoinverse.
Alongside this, the prediction and the state estimate are unbiased and provide the following estimation accuracy: It is essential to understand that the relations (9) and (10) determine the conditionally optimal Pugachev filter [42][43][44], the linear structure of which is postulated initially. Accordingly, the CMNF concept complements this filter with a minimax justification of the structure. The second element of the CMNF concept is a method of practical determination of the coefficients , , and ℎ by the Monte Carlo method. The practically used filter is obtained by replacing (10) mathematical expectations and covariances with their statistical estimates from computer simulation.
The primary a ention should be paid to the fact that for the presented model of a stochastic observation system with random delays, the CMNF estimate is determined by the same fundamental relations (8)-(10), initially wri en for the canonical form of the Markov system. At the same time, delays are not ignored, and their statistical properties are considered in the filter parameters (9). A more subtle account of the phenomenon of delayed observations can be provided due to flexibility in the choice of structural functions and . This will be demonstrated in the next section of the article. Here it remains for each of the variants of the model version (1) or (2) (note that entry (3) is a simple case of (2)) to formulate the conditions for the existence of solution (10). This is conducted in the following statement. Then, the prediction and the estimate of the conditional minimax nonlinear filter (8) exist, are described by the relations (9) and (10), and ensure the quality (11).
Proof of Theorem 2. The statement is an adaptation of the existence theorems of the CMNF estimate formulated in [38]. The proof is based on two propositions. First, it is the solution of the minimax problems (8). The worst distribution of = ( , ) on the class ℱ( , ) is the normal distribution. Therefore, the minimax estimate of is the estimate of normal correlation, i.e., a linear function of observations . Secondly, it is necessary to check the sufficiency of the theorem conditions for the existence of the second-order moments of the functions used, i.e., covariances of vectors ( , ) and ( − , ) . This becomes obvious if we use the canonical Markov form for writing (1) and (2) using the extended vector defined in (4). Since ‖Θ ( )‖ = 1, the existence of the necessary moments follows from the existence theorems [38]. □

Observation System Model
For the practical study of the proposed model with random delayed observations and the filtering algorithms (9) and (10), a simple movement model in a plane with constant velocity is chosen. This example is of a model nature and does not pretend to be practical, since at this stage of research, it is important to make sure that the proposed methodology is fundamentally effective. For this reason, the model parameters are given "extreme" numerical values for practice. In particular, this applies to average speeds and maneuvering speeds. This technique puts the algorithm in much worse conditions than can be the case in practice.
It is assumed that the AUV moves at a depth in the horizontal plane . Random noise is added to the current speed at each moment, simulating chaotic maneuvering. The coordinates of the AUV trajectory are denoted, as is customary, as ( ) and ( ). Note that these designations should be distinct from the designations of the state vector and observations in the original models (1), (2), or (3). The unit of position measurements is in kilometers (km), time is in hours (h), and speeds are measured in km/h. The continuous movement model is sampled in increments of = 0.0001 h, corresponding to a frequency of about three measurements per second. For the convenience of a graphical representation of the results, it is assumed that the movement begins at time = 0 and ends at time = 1000, i.e., the movement lasts 6 min. The systematic constant AUV speed values are = 25 km/h and = 12.5. During the movement, the AUV moves an average of 3 km. The Gaussian vector gives the initial position of the AUV (0), (0) , which has the mean (6.25; 12.5) and the covariance diag{2.5 ; 5 }.
The same type of observers are located at two points on the same plane: the first has coordinates 0, , = 6,25, i.e., it is located 6,25 km from the origin on the axis, the coordinates of the second are ( , 0), = 12.5, i.e., 12.5 km from the origin on the axis. Assuming that acoustic sonars are used for observation and that the speed of sound in water is equal to = 5400 km/h, it is possible to determine the maximum possible delay of observations = 50, i.e., the maximum value by which observations can be delayed is 18 s. Accordingly, the filtering process and the first estimate appear at time , i.e., 18 s after the start of movement (we can assume after the target is detected).
The source of the observer model is typical acoustic sonars [45,46] that measure the bearing and Doppler frequency shift. So, we can assume that the azimuth angle and range are measured, i.e., the AUV acoustic signature is considered known. Each observer measures the range and the direction cosine to the target with an additive measurement error. The error vector is assumed to be the Gaussian with mean (0, 0, 0, 0) and the covariance diag{0.001 , 0.005 , 0.001 , 0.005 }. One can use the work [47] to understand the accuracy parameters of real devices. The scheme of movement and observation is illustrated in Figure 1.
Describing chaotic velocity changes vector ( ), ( ) is assumed to be Gaussian with mean (0, 0) and two variants of covariances: (1) diag{0.0025 ; 0.0025 } and (2) diag{0.01 ; 0.01 }. Accordingly, in the first case, the AUV can change the speed by about 25 km/h, and in the second, by about 100 km/h. Thus, case (2) is about a high-speed maneuvering target. For the presented calculation, it is more interesting because it gives a very diverse trajectory.
In (14), the notation [ ] is used for the floor function of , and the potential possibility of removing an observed object at a distance for which the delay value becomes greater than the specified maximum is taken into account.

CMNF Structure
To implement the CMNF algorithm, it is required to determine the structure of the filter, i.e., to select the basic prediction function = ( ) and the basic correction function = ( , ). The linear model of movement (12) does not require clarification of the typical choice of the structure due to the system, i.e., We can also try and choose a residual as the basic correction structure, i.e., In theory, even the fact that the delay is not taken into account in the residual (15) should not interfere with the operation of the filter. However, the first experiments have already shown that the structure (15), although it gives a working filter, could have a higher accuracy. The apparent idea is to supplement (15) with a delay estimate of the form (14), i.e., but two different delays did not lead to an improvement in the result. A positive result was obtained when the observations' physical nature was considered in the correction structure. Namely, if we assume that there is no noise in (13), we can determine the AUV's position by recalculating the distances and cosines into Cartesian coordinates. Such a transformation gives the following correction structure: In addition to the geometry of observations in correction (17), the delay estimate (16) is also used. Thus, (17) are two estimates of the AUV position that each observer could give independently without considering measurement errors. Combining these estimates and taking into account the impact of errors is already the filter's task, i.e., of the relations (10).
It should be noted that the demonstrated flexibility in selecting the estimate structure is the key advantage of the CMNF.

Computer Simulations
In each of the performed computer computations, we simulated two independent beams of 100,000 trajectories of the observation system (12)- (14) and CMNF estimates (9). On the first beam, the filter parameters were calculated according to the formulas in (10) and their accuracy according to Formula (11). Accordingly, the statistical average replaced mathematical expectations according to the Monte Carlo method. On the second beam, the actual quality of the CMNF estimate = ( ), ( ) was evaluated. The mean square deviations of the estimate errors determine the estimation accuracy of the coordinates of the AUV position. These values, denoted by ( ), ( ), were calculated from the second beam of trajectories, and can be called practical accuracy.

The values
, and , represent the theoretical accuracy of the CMNF estimate. We calculated these values using Formula (11) and compared them with ( ) and ( ), respectively.
Calculations were performed for two variants of maneuvering speed (let us call them, conditionally, a maneuver of 25 km/h and a maneuver of 100 km/h) and two correction options: the main (17) and the residual (15). For a qualitative assessment of the results, each of the calculations was repeated under the assumption that = 0, i.e., in a model that does not consider the observation delays. The following series of figures illustrate the results. Figure 2 shows several AUV trajectories ( ), ( ) for a maneuver of 25 km/h (left) and a maneuver of 100 km/h (right), and the corresponding CMNF estimates ( ), ( ) calculated using correction (17). This figure and the following ones show the same trajectories from the simulated beam, i.e., the same set of random variables forms the states and estimates, and the difference is only in the value of the model parameters (in Figure 2, these parameters are variances of ( ) and ( )). This figure gives a qualitative idea of the filtering results and explains why the trajectories for a maneuver of 100 km/h are more interesting to study.
A formal analysis of filtering quality is presented in Figure 3. This figure gives quantitative estimates of the effectiveness of target tracking, shows the extent of deterioration in the quality of estimates between a maneuver of 25 km/h and a maneuver of 100 km/h, and also confirms the applicability of the Monte Carlo method for the practical synthesis (training) of CMNF, since the theoretical and practical accuracy are very close.
The same illustrations, as shown in Figures 2 and 3, are repeated in the following Figures 4 and 5 with a single change: instead of correction (17), we use the correction (15). These figures support the above statement about the correction in the form of a residual that is not very suitable for the model under consideration. Not only did the filtering quality turn out to be very low, but the difference between the theoretical and practical accuracy of the CMNF is also unacceptably large. The la er circumstance can be dealt with by increasing the size of the training beam used for approximate calculations using Formula (10), but it is impossible to improve the quality compared to the theoretical ones , and , in any case. The last group of figures is given in order to illustrate the influence of the main factor of the model proposed here: the delay of observations . We should immediately note that it is not advisable to consider figures of both variants of the maneuver, so for graphical presentations, we have limited ourselves to only the maneuver of 100 km/h, as it is more beautiful. Accordingly, the same trajectories and the exact figures above are given to compare corrections (15) and (17). Figure 6 shows the trajectories, and Figure 7 illustrates the filtering accuracy.    These visual and quantitative figures illustrate the real difference between models with and without observation delays. Nevertheless, they also show that the quality of a universal filter with a simple correction that needs to consider the specifics of the movement observation model needs to be higher. Of course, as seen in Figure 7 on the right, the difference between the theoretical and practical accuracy of the CMNF estimates will disappear as the size of the training beam increases. However, the absolute advantage of the "smart" correction will remain.
Finally, the calculations' results are presented in tabular form. Namely, Table 1 shows the average values of the quality indicators for all sets of parameters. The average value is understood as averaging over time, for example E { ( )} = ∑ ( ) , etc.

Conclusions
Summing up the article, we repeat the main provisions:  In target tracking tasks where acoustic sonars are used, it is necessary to take into account the random delay of observations according to model (1) or its simplified versions, (2) and (3);  The optimal solution of the tracking problem for the simplest model, (3), described by a Bayesian filter, can be wri en out but cannot be implemented practically;  The method of conditionally minimax nonlinear filtering can be applied to a model with random observation delays;  In problems with observers of the "geometric" type, the use of simple universal filters can lead to unsatisfactory results even in simple observation conditions; taking into account the physical nature of observations can dramatically improve the quality of filtering.
The presented relatively complete picture of the model with random delays of indirect observations and the problem of filtering its state includes a theoretical description of the model and the optimal Bayesian filter, a practically consistent algorithm of conditionally nonlinear minimax filtering, an applied experiment with a model of AUV movement close to practice and observations of it with typical sonars. At the same time, the following question remains open. Many parameters used in the experimental model can indeed be known or established experimentally. But part of the parameters, and above all, the systematic component of the AUV velocity = , , cannot be known in advance. Thus, the possibility of identifying these parameters in parallel with goal tracking comes to the forefront. In addition, the constancy of speed also does not look like a realistic assumption. So, for practical purposes, we need to learn how to identify the speed in a reasonably short observation period and repeat the identification process as necessary. Further research is supposed to be devoted to this issue.