Adaptive Stochastic Filtration Based on the Estimation of the Covariance Matrix of Measurement Noises Using Irregular Accurate Observations

: In measurement systems operating under various disturbances the probabilistic characteristics of measurement noises are usually known approximately. To improve the observation accuracy, a new approach to the Kalman’s ﬁlter adaptation is proposed. In this approach, the Covariance Matrix of Measurement Noises (CMMN) is estimated by accurate measurements detected irregularly by the mobile object observation system (from radiofrequency identiﬁers, etalon reference, ﬁxed points etc.). The problem of adaptive estimation of the observer’s noises covariance matrix in the Kalman ﬁlter is solved analytically for two cases: mutual noises correlation, and its absence. The numerical example for adaptive ﬁltration of complexing navigation system parameters of a mobile object using irregular accurate measurements is given to illustrate the effectiveness of the proposed algorithm. Coordinate estimating errors have changed in comparison with the traditional scheme from 100 m to 2 m in latitude, and from 200 m to 1.5 m in longitude.


Introduction
Currently, one of the central problems in the theory and practice of measurements is the development of algorithms for processing measurements [1], providing a given accuracy of information extracted from them under conditions of noises of various physical nature. This problem is particularly acute for stochastic dynamic processes with dynamic noisy measurements [2], moreover, in various fields of measurement technology: • measurements of state parameters of complex mechanical systems [3]; • measurements of electric motor characteristics [4,5]; • photogrammetric measurements [6]; • estimation of spacecraft orientation [7]; • measurements of motion parameters of unmanned objects [8,9]; • measurements of chemical process characteristics [10], etc.
Among the numerous methods of processing dynamic measurements under noisey conditions, focused on various specific features of certain measured processes [3,[11][12][13], the most universal and effective today are the methods of the stochastic filtration theory [4,11,[14][15][16][17][18][19], which provide an optimal evaluation of the measured process by a given criterion.
The use of filtration theory approaches including the Kalman filter [20], for the assessment of dynamic stochastic processes assumes an accurate initialization of the random noises of these processes [14,16]. At the same time, in real information and control systems exposed to various disturbing effects, the meters' stochastic noises are recognized, as a rule, approximately or fluctuate randomly [10,19,[21][22][23][24][25][26][27]. As a consequence, one of the very critical Kalman filter characteristics is the Covariance Matrix of Measurement Noises (CMMN), which straightforwardly influences the filter gain change and, consequently, the rate of convergence of the filtration process. In order to avoid a priori uncertainty, various approaches have been designed among which the following can be recognized as the main ones: 1.
introduction of empirical scale coefficients for calculating the posteriori covariance matrix and the CMMN-matrix [16,25,28,29] using the variational Bayesian approximation to construct algorithms for interacting multiple models and model-conditioned estimates [36,37] etc.
A serious drawback of the first and second approaches is the lack of strict selection criteria and scale factors, respectively, in the procedure of their calculation. The disadvantage of the third approach is the impossibility of adaptive estimation in real time, since it is necessary to pre-calculate the covariance of the updating sequence. A common disadvantage of fourth, fifith, and sixth approaches is a significant increase in computational costs.
This leads to the need to develop a fundamentally different approach that allows us to ensure the required accuracy of the Kalman filter by using an adaptive algorithm for estimating the CMMN carried out directly in the process of current filtering of the dynamic system state vector.

Task Definition
The solution of this problem is considered for a wide class of information-measuring and control systems, in which the measurements of "rough" sensors with unstable and significant measurement errors are corrected by the measurements of significantly more accurate sensors, which are considered references. The correction procedure is carried out, as a rule, at time intervals significantly exceeding the measurement step of "rough" sensors and often is irregular (random).
Examples of such systems include the following: • inertial-satellite Navigation Systems (NS), where the measurement correction of inertial NS is based on satellite NS measurements; in this case, the measurement errors of the inertial NS increase with time, and the satellite NS measurements are considered as accurate measurements of the velocity vector and coordinates of the object [38,39]; • different robots' NS, in which the correction of the navigation parameters of the robot (or person) is subject to a zero speed of his feet (or bottom of wheel) at the moment of touching the surface of the earth [40]; • information-measuring systems of different transport systems (maritime, rail, etc.), where the correction of the parameters of orientation and navigation upon the completion of a baseline (reference) point with exactly known coordinates [41,42]; • combined NS based on inertial sensors, providing navigation inside premises, and closed spaces [43,44], etc.
In connection with the above, the possibility of using the obtained accurate observations to construct an adaptive algorithm for forming the CMMN is considered. This algorithm allows us to significantly increase the accuracy and stability of the process of evaluating the state vector of the system by the Kalman filter as a whole in the time intervals between them [44].
Since accurate observations are made at discrete time points, a Kalman's filter discrete version is further considered to solve the problem. In this version, the system state vector estimation X k at a k-th time moment is determined by the equation [14,16,[28][29][30]44]: whereX k−1 is the system state vector estimation at the (k−1)-th time moment; Φ k is the transition matrix of the system state; H k is the measurement matrix that maps the space of system state vectors to the space of measurement vectors; V k is the measurement interference approximated further by a centered Gaussian sequence with an unknown covariance matrix R k , estimated from accurate observations; K k is the filter gain defined as where P k/k−1 is the extrapolated covariance matrix that characterizes the error value of the a posteriori state vector estimation; Q k is the covariance matrix of system noise that characterizes the level of its impact on the system.
Based on the filter gain Equation (2), it is necessary to formulate the objective of adaptive assessment of the desired CMMN R k from irregular accurate measurements. In this case, the matrix is determined from the condition that the vector of estimatesX k Equation (1) coincides with the exact state system vector X k (i.e., accurate observations) at the time of accurate observations receipt.
The solution of the problem is carried out in two stages. In the first stage, an adaptive algorithm is developed for the case of no cross-correlation of noise in the components of the measurement vector, that is, for the diagonal CMMN. The adaptive algorithm development for cross-correlated noise in the measurement vector for the complete (non-diagonal) CMMN is performed at the second stage.

Adaptive Filtration Algorithm for Uncorrelated Noises in the Measurement Vector
To solve this problem, we will use the full formula of the Kalman gain obtained by inserting Equation (3) into Equation (2): where for convenience of the subsequent decision we will enter notation Φ k · P k−1 · Φ T k + Q k · H T k = γ k and we will write down the coefficient K K formula as follows: In this case, Equation (1) will take the form And, with respect to the matrix R k , is a nonlinear vector equation that requires the solution of traditional numerical methods to use a very expensive multiple procedure of matrix inversion. For the possibility of Equation (5)'s analytical solution, we will carry out the following.
Since under the task conditionsX k = X k , then, entering the notation Multiplying both parts of Equation (6) by the inverse matrix, we have The standard inversion of the matrix γ k −1 is possible only in the case of a square matrix H K ; so, to preserve the generality of reasoning, we use the pseudo-inversion operation γ k + [46]. To simplify the subsequent solution, we denote the vector γ + k X * as ∆ * . Then Equation (7) easily admits an analytical solution with respect to all elements of the diagonal matrix R K ; if we consider the possibility R k ∆ * in the form ∆ * diag R kvect , where . . , n elements of the vector ∆ * and matrix R K , respectively. In this case, we have the desired equation of the vector of elements of the dispersion matrix of measurement noises in the form where ∆ −1 * diag is the inverse matrix, easily calculated due to its diagonality: Thus, the found solution of the nonlinear vector Equation (5) in the form of Equation (8) allows us to analytically solve the problem of adaptive estimation of the CMMN from accurate observations.

Adaptive Filtration Algorithm for Correlated Noises in the Measurement Vector
In the more general case of forming a measurement vector-when the interference of its components is mutually correlated, the covariance matrix R K may be non-diagonal, which somewhat complicates the procedure for its calculation. First, with the dimension of the state vector N, the dimension of the matrix R K is generally equal to N 2 , which requires additional accurate measurements at N time steps. In practice, this possibility is often quite feasible-on roads and railways with a developed system of various road landmarks, in inertial satellite navigation systems with a high frequency (up to 100 Hz) transmission of satellite messages, etc.
Secondly, in order to be able to analytically determine all the components of the matrix R K , the obtained measurements must be regrouped accordingly, followed by the solution of a new system of equations. In this case, the algorithm construction scheme takes the following form.
The exact measurements obtained at N time steps transform the estimation scheme as follows: . . .
For each i-th step of the evaluation, transformations similar to Equations (5)-(7) are performed, after which the Eq's system (9) is transformed into the system where the matrix R K is assumed to be constant at all N time steps Next, the Equation (10) is transformed as follows: from each i-th subsystem of the equations U * i = R K · ∆ * i of dimension N, the j-th equation is selected, j = 1,...,N: where U * ij is the j-th component of the vector U * i , R Kj is the j-th row of the matrix R K . After that, a system of N linear equations is formed to determine all N components of the j-th row of the matrix R K : In vector form, the Eq's system (11) can be written as . . . T |∆ * 1 ∆ * 2 ∆ * 3 · · · ∆ * N | −1 = R Kj and, accordingly, the matrix R K is Here, as before, the task of adapting the Kalman estimation algorithm based on accurate observations is solved analytically.
In conclusion, it should be noted that the property of symmetry of the matrix R K inherent in it by definition, on the one hand, reduces the number of its unknown components from N 2 to N 2 + N 2 , but, on the other hand, when solving (N − 1) systems of equations reduced with respect to Equation (11), requires the definition of (N − 1) inverse matrices, although reduced with respect to the matrix |∆ * 1 ∆ * 2 ∆ * 3 · · · ∆ * N | −1 = Λ, but, nevertheless, requiring significant computational costs compared to a single calculation Λ. This circumstance allows us to consider the algorithm described above as the most efficient in terms of computational costs.
Next, consider an example that illustrates the effectiveness of the proposed approach.

Example
The movement of any object in the geographical coordinate system (GCS) is described by the following equation [45]: where φ, λ-the object's latitude and longitude accordingly; r is the radius of the Earth; h is the object height; V y , V x are the linear velocity projections on the GCS-axes. Next, we consider that the beginning coordinates of the object's trajectory are defined as ϕ 0 = 0.8 rad, λ 0 = 0.3 rad, time interval of movement (0; 1000) s; the object moves at a constant velocity V = 20 m s −1 along a loxodromic trajectory with an azimuthal angle at a constant speed along a loxodromic trajectory with an azimuthal angle A = 0.2 rad. In the process of movement, the topography of the Earth's surface leads to a random change in the height of the object with zero expectation and standard deviation of 0.15 m. Based on the nature of the object's trajectory, the projections of its linear velocity on the GSC-axes can be determined as follows: We further assume that the navigation system (NS) of the object is built on the basis of weak integration of inertial NS (INS) and satellite NS (SNS). The measurement complex of this integrated NS receives satellite measurements at intervals of 20, 15 and 30 s, which are considered accurate. At other times, to estimate the vector of navigation parameters with a clock cycle of 0.1 s, INS measurements are used for channels λ and ϕ, for which the CMMN has the R k = 4, 2 · 10 −10 0 0 9, 5 · 10 −10 (rad) 2 . Due to the fact that the speed of movement of an object causes a slight change in its coordinates at a given time interval, it is possible to use linearized navigation equations Equation (12) [45] relative to the initial coordinates for the subsequent solution of the problem .
Since, in accordance with the example conditions, the change in the object height is a random process with zero expectation and a standard deviation of 0.15 m, then we assume h 0 = 0. Reducing the Equation (13) to a discrete form by the Eulers method with a sampling step τ: we have a filter that in the accepted notation has the form: , Since the CMMN R k is unknown under the initial conditions, the estimation of navigation variables was conducted for two cases: when the value of R k is set with an error over the entire modeling interval R k = 1, 1 · 10 −10 0 0 2, 5 · 10 −10 (rad) 2 , and when the R k value is estimated according to the proposed algorithm.
Graphs of estimation errors obtained when implementing the Equation (14) in the first version are shown in Figures 1 and 2, where a significant error in the estimation of the current coordinates of the object is obvious: the latitude error lies in the range from 100 to 350 m, the longitude from 50 to 200 m.  When designing an adaptive filter in accordance with the above algorithm, the estimation errors have the form shown in Figures 3 and 4.  It is obvious that when using the proposed algorithm, the estimation errors in comparison with the traditional scheme changed sharply to the range from 2 to 7 m in latitude, and from 1.5 to 4 m in longitude.
For a more complete assessment of the capabilities of the proposed approach, it was compared with the innovation-based adaptive estimation algorithm [47] with the same model parameters and errors in setting the CMMN. When using innovation-based adaptive estimation, the covariance matrix R k was evaluated in accordance with the equation [47] where υ k = Z k −Ẑ k/k−1 is the discrepancy between the received measurement and its predicted value; N is the number of samples in the sliding "window" of the measurement in which the covariance matrix R k is estimated; j 0 = k − N + 1 is the initial position of the sliding "window". In the case under study, N was chosen to be 25, and the sliding window was used three times for the 50th, 250th, and 800th seconds. The character of changes in estimation errors for both coordinates remained the same as in the developed algorithm, but the estimation errors themselves increased significantly in comparison to the range from 20 to 37 m in latitude and from 15 to 28 m in longitude.

Conclusions
A comparison of all the considered variants of traditional and adaptive filtration clearly shows the advantages of the proposed algorithm, despite the slight increase in computational cost (very small relative to total costs for a traditional filter; the Central Processing Unit (CPU) operation time for the combined calculation of the covariance matrix estimation and the current step of the filtration equation (14) was 1.849 s, and without the calculation of the covariance matrix estimation was 1.253 s). The simplicity, stability, convergence, and accuracy of the proposed algorithm make it possible to effectively use it for a wide range of measurements: measurements of motion parameters of unmanned objects, assessment of spacecraft orientation, measurements of chemical process characteristics, measurements of state parameters of complex mechanical systems, and others.