Abstract
This paper presents a noise covariance estimation method for dynamical models with rectangular noise gain matrices. A novel receding horizon least squares criterion to achieve high estimation accuracy and stability under environmental uncertainties and experimental errors is proposed. The solution to the optimization problem for the proposed criterion gives equations for a novel covariance estimator. The estimator uses a set of recent information with appropriately chosen horizon conditions. Of special interest is a constant rectangular noise gain matrices for which the key theoretical results are obtained. They include derivation of a recursive form for the receding horizon covariance estimator and iteration procedure for selection of the best horizon length. Efficiency of the covariance estimator is demonstrated through its implementation and performance on dynamical systems with an arbitrary number of process and measurement noises.
1. Introduction
Filtering and system identification are powerful techniques for building models of complex control systems in time and frequency domains [1–5]. The Kalman filtering and its variations are signal processing techniques in widely used applications, such as navigation, target tracking, vehicle state estimation, communications engineering, air traffic control, and many others [6–9]. As the well-known standard Kalman filter (KF) equations do not directly depend on the original process and measurement noise covariances and , respectively. They depend on the transformed covariances,where and are the rectangular gain matrices associated with process and measurement noises, respectively [1, 2, 6, 7]. In this paper, the original and transformed noise covariances shall be denoted by “ONC” and “TNC”, respectively.
The optimality of KF depends on the accuracy of TNC, which directly depend on the prior assumptions on the ONC. Inadequacy of the prior values for the TNC can lead to unexpected results and the KF divergence [1, 2, 7]. The adaptive filtering is one of the approaches to prevent the divergence problem when precise knowledge on the noise covariances is not available.
Adaptive Kalman filtering (AKF) is not the subject of this paper, so we present only its brief survey. The AKF represents a promising strategy for dynamical adjustment of the TNC for online estimation. Two popular types of the AKF algorithms are the innovation-based estimation approach [10–13] and the adaptive fading filtering approach [14, 15]. The latter is a type of covariance scaling method, into which suboptimal fading factors are incorporated. The innovation-based estimation approach calculates and/or , assuming that the innovation sequence of the KF represents a white noise. In this approach, one of the covariance matching, the correlation, and the maximum likelihood methods is used. In the covariance matching method at first the sample error covariance of a state is computed based on a moving window and then is estimated using the calculated sample covariance. In the same way, is computed, but using the sample covariance of the KF innovation sequence. In the correlation method, the noise covariance matrices are estimated based on the sample autocorrelations between the innovations by exploiting the relations between the estimation error covariance and the innovation covariance [11, 16, 17]. The drawback of the covariance matching and the correlation methods is that they do not guarantee the positive definiteness of the calculated covariances and require a large window of data. Lastly, implementation of the maximum likelihood method is very complicated which makes it impractical.
As noted the AKF procedures estimate not the ONC, but the TNC. Therefore the estimation of the ONC remains open. Knowledge of the ONCs is necessary to study the quality of a system model. The resulting ONCs help to understand how much environmental disturbance affects the behavior of dynamics and to evaluate the accuracy of sensor measurements. However, if the process or measurement (“noise gain matrix”) is rectangular, we then have both theoretical and practical challenges in calculating the ONC. As will be shown shortly, popular approaches including the usage of the standard least squares approach and current data may lead to unsatisfactory results. To the best of the authors’ knowledge, there are no existing results for the noise covariance estimation within the Kalman filtering framework utilizing not only current data, but finite-memory information from the most recent time interval (receding horizon).
The primary aim of this paper is to expand the previous results concerning the standard least squares covariance estimation to a new robust receding horizon covariance estimation. We propose a novel general estimation method for calculation of the ONCs with rectangular noise gain matrices and The method combines the least squares (LS) optimization and receding horizon strategy. The proposed estimator can be integrated into the AKF for simultaneous estimation of both transformed and original noise covariances.
The main contributions of the paper are listed in the following:(1)A novel receding horizon least squares optimization criterion is proposed. Using this criterion, an optimal covariance estimator for an arbitrary structure of noise gain matrices is derived.(2)Constant noise gain matrices are comprehensively investigated, including the derivation of a recursive procedure for the optimal covariance estimator with a horizon interval and a stopping condition during the iteration process.(3)Performance of the proposed estimator is illustrated on theoretical and practical examples for proving its correctness and efficiency.
This paper is organized as follows. In Section 2, two motivation examples within the Kalman filtering framework are demonstrated. Section 3 presents the standard LS noise covariance estimator and describes consequences of processing current (single) data only. In Section 4, a novel receding horizon least squares estimator (RHLSE) is proposed. It is based on the generalized LS approach and the novel receding horizon LS criterion. Here, the cardinal equations for determining the RHLSE and its equivalent modifications are derived and discussed. In Section 5, the algorithm of searching for the best horizon length improving the accuracy of the RHLSE is proposed. Section 6 demonstrates the practical usage of the RHLSE in adaptive Kalman filtering. In Section 7, performance and effectiveness of the proposed estimator are studied via a theoretical example with rectangular noise gain matrices and a motion with random velocity. Finally, we conclude the paper in Section 8. The list of main notations is given in Table 1.
2. Motivation Examples with Rectangular Noise Gain Matrices
To begin, two typical motivation examples of a dynamical system with rectangular noise gain matrices are presented.
Let us describe the noise estimation problem of present interest via two motivation examples. Suppose we have a linear discrete-time system,Then the TNCs and represent covariance matrices of the linearly transformed noises and with the rectangular gain matrices and , respectively, i.e., and .
Motivation 1 (“colored process noise”). Consider the scalar linear system with colored process noise:where is a measurement white noise.
Suppose that is a colored process noise generated by the autoregressive equation of order , ,where is a white noise.
Equation (4) may be considered as some system shaping the colored noise from a white noise . Then the following is a state-space representation of the model (4):whereWe combine (3) and (6) to obtain a system model with –dimensional augmented state and a scalar process noise ,where
Motivation 2 (“composite process white noise and colored measurement noise”). Consider the following discrete state-space model,where , , and , are two uncorrelated process white noises with covariances and , respectively; and
The state equation (9a) can be written in standard form with one composite process white noise , and rectangular gain matrix ,where , , ,
Following [2] we define an auxiliary measurement described by the new measurement equation with white noise ,where , ,
Combining (10) and (11) we obtain the standard structure of the system model (2) with rectangular noise gain matrices and
3. Least Squares (LS) Estimation of Noise Covariance Using Single Data
As noted above, the purpose of this paper is to calculate the unknown ONCs and based on the available TNCs and obtained by the AKF algorithm. In this section, we present a simple LS solution to this problem.
For a clear presentation, let us first rewrite (1) in the following form:whereIn a special case where the gain matrix is square and nonsingular, the LS solution to (12) takes the form , where the known matrix represents a current TNC at time instant But in practice, the gain matrix is rectangular, as the two motivation examples illustrated. In this rectangular case, the LS solution to (12) can be obtained using a pseudoinverse matrix.
The linear matrix equation (12) has been considered by many authors [18–21]. They have studied general solutions to (12) for special solution structures, e.g., symmetric, triangular, or diagonal using matrix decomposition techniques such as the singular value decomposition (SVD), the generalized SVD, and the canonical correlation decomposition.
Theorem 3 (see [18]). The least squares solution to matrix equation (12) takes the form
Note that the Moore–Penrose pseudoinverse always exists and is unique.
However, the known TNC in practice is usually obtained from an experiment. It is calculated by using the AKF. Since the LS solution uses only the single (approximately calculated) TNC , it is likely to be unstable or inaccurate. Therefore to overcome this disadvantage we propose a more realistic idea to use the set of data (TNCs), , for the calculation of unknown covariance . This set of data is obtained over the most recent time interval (receding horizon), , saving in a finite memory. Using these datawe can improve the LS solution (14), and, as a result, one can expect that the obtained receding horizon LS solution will be more robust against processing and measurement errors than the LS solution which utilizes only the current data
4. Receding Horizon Least Squares Estimator (RHLSE)
In this section, we propose the novel receding horizon least squares estimator referred to as RHLSE.
In literature, the estimation problem in dynamical systems having various uncertainties has been solved in different ways, including finite-memory (receding horizon or sliding window) estimations [7, 22]. As a result, the finite-memory estimators are known to be more robust against model uncertainties and numerical errors than the standard LS estimator or Kalman filter [23–25]. For this reason, this receding horizon strategy is chosen to improve the LS solution (14).
According to the strategy proposed above, the unknown solution minimizes a new receding horizon LS criterion representing the sum of squared residuals over the receding horizon , i.e.,withwhere is the receding horizon length (HL).
Theorem 4 (RHLSE for time-varying noise gain matrices). The solution of the receding horizon least squares optimization problem represented by (17) satisfies the equationwhere
The proof of Theorem 4 is given in Appendix.
Next consider the important case in which the noise gain matrix is time-invariant. Then the receding horizon LS problem (17) is rewritten aswithand we get the following result.
Theorem 5 (RHLSE for time-invariant noise gain matrices). The optimal receding horizon least squares optimization problem for constant noise gain matrices represented by (20) has the explicit solution,where is the horizon average covariance,
The proof of Theorem 5 is given in Appendix.
As we see from (21) and (22), the explicit solution depends only on the single data representing the horizon average covariance Therefore we get the following.
Corollary 6. The RHLSE (21) represents the LS estimator (14) calculated for the horizon average covariance, , i.e.,
5. Selection of Best Horizon Length for Time-Invariant Noise Gain Matrices
The optimal noise covariance (21) is sequentially calculated for each HL (iteration), . Here we propose a stopping criterion for selecting the best HL, This stopping criterion is based on the errors between iterations. Let us first represent the solutions (21) and (22) in recursive form.
The recursive formula for the horizon average covariance (22) takes the formUsing (21) and (24) the explicit solution is also rewritten in recursive form,whereFor the stopping criterion based on the norm of iteration errorwe have the following result.
Theorem 7. Let be the RHLSE (21). Then the Frobenius norm of the iteration error is calculated as
The proof of Theorem 7 is given in Appendix.
Using (28), the best HL is defined aswhere is a given tolerance.
Remark 8. As noted in Section 3, the initial covariance , calculated at , may be inaccurate. Numerical examples given below confirm this fact. Therefore to improve the covariance , we need to continue the iterative process , using the stopping condition (29).
Remark 9. In general, the stopping condition (29) does not guarantee the convergence of the iterative process (25); that is, the inequality may not hold. Besides, the number of iterations is limited, . However, (29) seems to be still useful for stopping the iterative process, as will be demonstrated by examples.
Remark 10. The stopping condition (29) can be carried out in real time. In fact, the term depends on the TNCs and for a given horizon. Therefore, the complexity of the real-time implementation of the stopping condition is not an issue.
6. Application of RHLSE
In this section, an application of the RHLSE to the estimation of process and measurement noise covariances is considered.
6.1. Kalman Filtering
The Kalman filtering framework involves estimation of the state of a discrete-time linear dynamical system with additive white Gaussian noise,where state , measurement , process noise , , and measurement noise , All matrices and are with appropriate dimensions. We assume that the initial state , and process and measurement noises are mutually uncorrelated with constant covariances (ONCs), and , respectively.
As mentioned in Introduction the KF equations depend on the known TNCs and , rather than the ONCs We have
6.2. Adaptive Kalman Filtering
In a real problem the ONCs and, as a consequence, the TNCs are unknown. In this case the AKF simultaneously estimates both the state and TNCs . The general scheme for adaptive estimation of the TNCs can be presented in Table 2.
Then the AKF equations for the system model (30) are given byTo estimate the unknown ONCs using the available TNCs we apply the RHLSE. Equation (21) for the corresponding estimates and areHere and are calculated in (32), and the best HLs for and for are determined by (29), i.e.,andrespectively.
7. Simulation Examples
In this section, we first test the efficiency of the proposed RHLSE and the stopping criterion on a theoretical example for which the true ONC is known (Section 7.1). And then, performance of the RHLSEs (33)-(35) shall be demonstrated with the example of a moving object (Section 7.2).
7.1. Theoretical Example
Let us estimate a random constant given two single sensor measurements and of corrupted by their own uncorrelated Gaussian white noises (measurement errors) and , respectively. We assume that the two sensors contain a common environmental Gaussian white noise , which is uncorrelated with and Then the dynamic equations describing this situation areHere we calculate the estimation accuracy of the proposed RHLSE in terms of the receding HL, , for the two measurement models.
Model 1 (“measurement model without environmental noise ”). or in matrix formHere the ONC and TNC take the formrespectively.
For the model (36)-(39) the measurement noise gain matrix , and and Then using (33) and (34) the RHLSE for the measurement noise covariance and the best HL take the following form,whereand is chosen empirically to give some statistical smoothing.
The simulation results are illustrated in Table 3 and Figure 1 for , , , and Our point of interest are the values of different norms of the true error , and the best HL In Table 3, a sign represents one of the five considered matrix norms, , , , , and (Note: for symmetric matrices). The columns II and III show the absolute norm (A-norm) and the relative norm (R-norm) of the iteration error , respectively.

Comments on Table 3 and Figure 1 are in order.(1a)The values of norms in Table 3 confirm the known inequalities for matrix norms of the symmetric matrix , i.e.,(1b)The best HL for each of the presented norms is the following, ; ; and So a small HL, , is suitable for Model 1 in which the sizes of the ONC and TNC are equal, (1c)The values of all norms of the true error confirm Remark 8 (see columns IV~VII). They show that the initial HL is not optimal, i.e., the single ONC is not good.(1d)The curves of all norms of the true error have a similar shape. The curves slightly fluctuate and then increase around a small HL, .(1e)The curves of the A- and R-norms are similar. Therefore, only one of them is shown in Figure 1.(1f)As noted in Remark 9, fulfillment of the condition does not guarantee the convergence. In Table 3 we can observe that the small iteration errors , correspond to the large true errors, This fact is also illustrated in Figure 2, where F-norm and scaled A-norm are plotted for a clear comparison.

Next, we consider more complicated Model 2, in which the noise gain matrix is rectangular, and the sizes of the original and transformed measurement covariances are different. In this case (12) represents an undefined system of equations.
Model 2 (“measurement model with environmental noise ”). or in matrix formHere the ONC and TNC take the formrespectively.
For the model (43)-(45),and the RHLSE (33) takes the form
Setting up , , , and , we calculate the norms of iteration error and true error The simulation results are illustrated in Table 4 and Figure 3.

In general, comments (1a)~(1f) made for Model 1 are the same for Table 4 and Figure 3, but there are some minor changes concerning (1b) and (1d).(2b)The best HLs for the norms are little different, i.e., , , , So a small HL, , is also suitable for Model 2, when the noise gain matrix is rectangular.(2d)The curves of all norms of the true error have a similar shape. The curves increase and then slightly decrease around a medium HL, .
Here is a summary of the simulation results presented for Models 1 and 2.(i)We propose to use the RHLSE with a small HL, , for both square and rectangular noise matrices Usage of short and long horizon lengths may lead to unsatisfactory results.(ii)In parallel with testing the efficiency of the stopping criterion we also study the behavior of the mean square errors (MSEs) for the covariance estimates (Model 1) and (Model 2), where is the true simulated value, is the RHLSE estimate of the element is the time instant, and is the best HL, , To calculate the MSEs the Monte-Carlo simulation with 1000 runs was used, i.e., the expectation operator in (47) is taken for the 1000 error data at each time instant . We observe that the relative errors, , for Model 1, range from 10.8% to 13.7% within the time zone , and then they decrease. At the maximum values of the and reach the values of and , respectively. Similar results, but slightly worse due to the rectangularity of the noise gain matrix , are obtained for Model 2. For this model the relative errors, range from 12% to 16.5% within the time zone . At the maximum values of the , , and amount to , and , respectively. Thus, this analysis of the MSEs shows that the RHLSE is suitable for practical applications.
7.2. Numerical Example: 1D Tracking with Random Velocity
In this section, the proposed RHLSE is used in adaptive filtering with unknown ONC ( and ).
Let the true target state at discrete time be defined as , where and denote the position and velocity in Cartesian coordinates, respectively. Assume that the constant velocity is subjected to a random disturbance The measured position can be defined as , where is a sensory error and is an environmental noise. The target state dynamics is then modeled as follows:Simulation results with usage of the AKF (32) and RHLSE (12) are illustrated in Figures 4–8 for the following model parameters:Figures 4 and 5 show the position and velocity estimates compared to the true values. The AKF appears to converge, in the sense that the estimation error tends to decrease, but slowly for velocity because the sensor measures only the position.





The key value in adaptive filtering is the TNC This covariance is plotted in Figure 6.
As we observe, the covariance is not seriously changing for . For this reason we fixed the receding horizon interval for further estimation of the ONCs and by formulas (33). We haveThe best HLs and are calculated based on the iteration errors and , respectively,andAs indicated earlier in Remark 9, the best HLs and do not guarantee a minimum of the true errors and , respectively. Therefore, in addition we calculate the norms of the true errors, and , which confirm that the proposed idea of choosing the best HL using the iterative error (53) or (54) leads to good results. Figure 7 presents the norms of the true error and the iteration error for measurement covariance, respectively. Analogously, Figure 8 presents the norms of the true error and iteration error for process covariance, respectively. We observe that as well as in the theoretical example (Section 7) the obtained best HLs and are located within a small horizon interval. In Figure 8, both curves show one value for the best HL, But in Figure 7 we observe five horizon lengths , which satisfy condition within the horizon interval , and the last four lie in the middle of the horizon interval (see Table 5). Table 5 shows that the iteration errors, , are different; however, the true errors remain almost the same. Hence, it can be concluded that the small length is acceptable.
To compare the relative MSEs,the Monte-Carlo simulation with 1000 runs was performed. The numerical values of all relative errors (55) within the time interval range from 13.8% to 16.7%. And for large time instants, , they do not change much. In this case, the maximum values of the errors , , and amount to 14.5%, 12.7%, and 13%, respectively.
8. Conclusion
In some application problems, knowledge of noise covariances is necessary to study quality of the system and measurement model. In order to estimate arbitrary process and measurement noise covariances, a receding horizon estimation criterion with a generalized least squares approach is proposed (Theorem 4).
Special attention is given for time-invariant noise gain matrices (Theorem 5). In this case the proposed RHLSE is comprehensively investigated, including the derivation of compact matrix form for recursive RHLSE (25) and an effective stopping condition for iterative estimation process (Theorem 7). The computational procedure based on the Moore–Penrose pseudoinverse for evaluation of the best noise covariances is recommended.
In view of the importance of noise covariances in practice, the proposed RHLSE is illustrated on theoretical and practical dynamical models with different types of noises. The examples show that the RHLSE with selection of the best HL within a small horizon interval yields reasonably good estimation accuracy.
Simulation analysis and comparison of the proposed estimator by several norms of estimation error show that usage of the receding horizon strategy is a very effective approach for the achievement of robustness against any uncertainties and errors.
Appendix
Proof of Theorem 4. Using the notations and , , the criterion (17) can be rewritten asDifferentiating each summand of the functional with respect to using formulas of the trace and matrix derivatives, and ,we getThen setting the result to zero, , we obtainThis completes the proof.
Proof of Theorem 5. Since the matrices and are constant, and after simple manipulations the general solution (18) takes the formNext, using the equality of pseudoinverse [18], we obtain (21).
This completes the proof.
Proof of Theorem 7. Since , and using (25) we haveorUsing the cyclic property of trace for arbitrary matrices we getNext using the cyclic property of trace of the product of three symmetric matrices , we obtainSubstituting (A.8) and (A.9) into (A.7) we getThis completes the proof.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the fund of research promotion program, Gyeongsang National University, 2017, the National Research Foundation of Korea (NRF) funded by the Ministry of Education (Grant no. NRF-2018R1D1A3A03000717), and the Ministry of Science and ICT (Grant no. NRF-2017R1A5A1015311).