Sensitivity of Algorithms for Estimating the Gravity Disturbance Vector to Its Model Uncertainty

The work considers the results of filtering and smoothing of the gravity disturbance vector horizontal components and focuses on the sensitivity of these results to the model parameters in the case when the inertial-geodesic method is applied in the framework of a marine survey on a sea vessel.


Introduction
The complexity of the Earth's surface and its interior is the reason why the direction of the gravity vector ! g (vertical line) does not coincide with the direction of the normal gravity vector ! n at points on the Earth's physical surface. This difference is characterized by the gravity disturbance vector (GDV) (Peshekhonov et al. 2017;Torge 2001). The problem of determining all three components of this vector is a problem of vector gravimetry. Note that the vertical component of the GDV coincides with the free-air gravity anomaly up to the accuracy of correction for the height anomaly (Jekeli 1999(Jekeli , 2016. The horizontal components of the GDV are deviations of the vertical, which are determined by two angles of deviation: in the planes of the Meridian Ÿ and the first vertical˜ (Peshekhonov et al. 2017;Torge 2001). Among the methods used to determine the GDV horizontal components on a moving base, including marine vessels, the gravimetric method, based on the measurement of gravity anomalies, and the astronomical-geodetic method, based on the comparison of astronomical and geodetic coordinates, and its modifications are particularly widespread (Koneshov et al. 2016;Peshekhonov et al. 1995;Hirt et al. 2010). Also well-known are the method O. A. Stepanov · D. A. Koshaev · O. M. Yashnikova ( ) · A. V. Motorin · L.P. Staroseltsev CSRI Elektropribor, JSC, ITMO University, St. Petersburg, Russia of gravitational gradiometry, based on the measurement of the second derivatives of the geopotential, and the method of satellite or aircraft altimetry, which makes use of measurements of trajectory altitude (Koneshov et al. 2016). We cannot but mention methods based on the construction of global models of the Earth's gravity field using the data on perturbations of satellite orbits (satellite-to-satellite tracking and other satellite missions). Combinations of these methods (for example, the astronomical-gravimetric method) are used (Peshekhonov et al. 2017;Koneshov et al. 2016). The inertial-geodetic method, which is one of the modifications of the astronomical-geodetic method, is used in the framework of the vector gravimetry problem solution. It is based on a comparison of astronomical latitude and longitude ®, oe, generated by the inertial navigation system (INS), and geodetic latitude and longitude B, L, obtained from the GNSS data (Koneshov et al. 2016;Emel'yantsev et al. 2015). Thus, the horizontal components of the GDV are defined as: In contrast to the conventional astronomical-geodetic method, in the inertial-geodetic INS method, the INS generates both astronomical coordinates and their derivatives, which makes it possible to use complementary velocity measurements (Koneshov et al. 2016;Dmitriev 1991).
The idea of using precision INS to determine horizontal components is based on the use of the relation between the errors of the INS output navigation parameters and the GDV components (Dmitriev 1997;Schultz and Winokur 1969;Li and Jekeli 2008). The analysis of this error, which is actually a methodical one, opens up the possibility, in principle, of determining the GDV horizontal components directly relative to the reference value at a reference point by solving the estimation problem using differences in the measurements of INS and GNSS.
The effectiveness of solving such a problem in real time (in filtering mode) was analyzed earlier in (Anuchin 1992;Nesenyuk et al. 1980;Staroseltsev and Yashnikova 2016). At the same time, experience shows that application of the smoothing mode, which does not require obtaining estimates of the desired parameters in real time, significantly increases the accuracy of the problem solution. Note that the feature of the smoothing problem is that when obtaining estimates at a current time moment, we can use both past and future measurements (Stepanov and Koshaev 2010;Loparev and Yashnikova 2012). The specifics of the inertial-geodetic method make it possible to implement optimal Kalman filtering and smoothing algorithms. However, such algorithms are based on the error models of the GNSS, INS and its sensors, as well as GDV horizontal components themselves in the form of random processes. To derive an adequate model of GDV horizontal components, it is necessary to have a priori information about the nature of the field in the survey area, for example, such parameters as the GDV variance and the correlation radius. Incorrect assignment of parameters can significantly reduce the accuracy of the resulting estimate. The analysis of the algorithm sensitivity is aimed to study the decrease in the estimation accuracy in the case that stochastic models in an algorithm are defined incorrectly Koshaev 2003, 2011). Sensitivity of algorithms for determining the vertical component of the GDV was discussed in (Stepanov and Motorin 2019;Stepanov et al. 2015a, b).
The paper considers the results of filtering and smoothing of the GDV horizontal components and focuses on the sensitivity of these results to the model parameters in the case when the inertial-geodetic method is applied in the framework of a marine survey on a sea vessel. The studies involve direct calculation of the covariance matrices of estimation errors (without using the Monte Carlo method).

Description of the Error Model
To solve the problems of filtering and smoothing in this study, we used the model of INS errors described in (Gusinsky et al. 1996(Gusinsky et al. , 1997. The state vector of this model includes errors in constructing the local vertical, errors in generation of the East and North velocity components, errors of the sensing elements, and the model of GDV horizontal components. Loosely coupled INS-GNSS system is assumed with coordinate, velocity, orientation angle, and sensor errors feedback for the filtering mode during survey. The fixed interval smoothing for every tack was fully done after survey. The approach to creation of an integrated system to measure GDV horizontal components implies very high requirements imposed on the instrumental errors of the INS. Gyroscope drift and measurement errors of linear acceleration must not exceed 0.0001 deg./h and 1 arcsec, respectively (Emel'yantsev et al. 2015;Nesenyuk et al. 1980). Therefore, the coefficients of the gyroscope error model were described by first-order Markov processes with the standard deviation (SD) at a level of 3 10 5ı /h and a correlation interval of 40 h; nonwhite-noise errors of the accelerometers were set by the stationary first-order Markov processes with an SD at a level of 1 arcsec and a correlation interval of 10 h; the additive white noise of the accelerometers was also set at a level of 1 arcsec/ p Hz. These processes define stability of inertial sensors. The initial root mean square errors (RMSE) of the sensors were set equal to the steady-state values of the corresponding Markov process SD. The RMSE of the INS initial alignment were set at a level of 1 arcsec for tilt angles and at a level of 0.1 m/s for the horizontal components of the velocity.
As noted above, for the solution of the estimation problem by the inertial-geodetic method in the formulation considered here, models of the GDV horizontal components must be given in the form of random processes describing the variability of the components along the vehicle trajectory (Staroseltsev and Yashnikova 2016;Nash and Jordan 1978). The expressions of spectral densities for the longitudinal S "L (!) and transverse S "T (!) horizontal components of the GDV in the rectilinear motion of the vehicle with a constant velocity V were assumed to be similar to those given in (Nesenyuk et al. 1980): where " ,˛", d are SD, damping coefficient, and inverse correlation radius, respectively. For the solution of the estimation problem, the model included measurements from the INS and the GNSS receiver, which are parts of the integrated system. It was assumed that the position measurements of the GNSS had white noise at a level of 3 cm/ p Hz and a component represented by a stationary first-order Markov process with a SD of 3 cm and a correlation interval of 1 h. Velocity measurements of the GNSS had only a white-noise error at a level of 3 cm/s/ p Hz.

Estimating the Accuracy of Determining the GDV Horizontal Components: The Results
To simulate the solution of the smoothing problem, we modified software in the Matlab package Koshaev 2003, 2011), which allowed us to calculate the RMSE of filtering and smoothing estimates and display their plots. Error-tolerant algorithms were used to calculate error covariance matrices based on UD decomposition. Owing to this, all the necessary calculations were performed without additional computational errors.
In the simulation, it was assumed that the vessel was moving North along the meridian at a speed of 10 knots; the initial inertial longitude was 30 ı E; the initial latitude was 60 ı N. The solution interval was 24 h, and the discreteness -10 s.
We set different characteristics of the gravitational field variability: from weakly disturbed with a of 5 arcsec and a correlation radius of 20 nmi to the highly disturbed with an RMSE of 20 arcsec and a correlation radius of 5 nmi. The steady-state values of the RMSE in the estimation of the GDV horizontal components obtained in the simulation for filtering and smoothing algorithms are shown in Table 1.
Using the smoothing mode to calculate the GDV horizontal components provides a significant gain in accuracy compared to the filtering mode. As may be seen from the table, for a weakly disturbed field (SD of 5 arcsec) only filtering algorithms can be used because there is practically no gain from smoothing. This fact is easy to explain: it is well known that the estimation accuracy of constant values (random bias) in the filtering and smoothing modes are the same. As for a strongly disturbed field (SD of 20 arcsec), the use of smoothing algorithms can give almost a two-fold increase in the accuracy of the estimation problem solution.

Estimating the Sensitivity of Algorithms for Determining the GDV Horizontal Components to the Model Parameters: The Results
For sensitivity analysis, we used the same parameters of the GDV model as in Table 1. When calculating the sensitivity for each combination of real (true) and estimated values of a model parameter, we determined three different RMSE in the estimation of the GDV horizontal components. The first one is the RMSE of the error in estimating the optimal algorithm, which corresponds to the estimate in the case that the model specified in the algorithm corresponds to the real one. This RMSE characterizes the maximum achievable estimation accuracy of the problem under given conditions.
The second one is the real RMSE that correspond to the real estimation accuracy of the algorithm tuned to an incorrect model of a disturbed field. Real RMSE calculated using method from Koshaev 2003, 2011). The method is based on the solution of the covariance equation for the augmented vector, which includes both optimal and suboptimal estimates of state vectors corresponding to correct and incorrect models. This RMSE will always be higher than the RMSE for the optimal algorithm.
The third RMSE are those calculated using the data from the covariance equation of the algorithms that stand for the characteristics of the accuracy of the obtained solutions. If for different values of the parameters of the real models in a certain range the calculated RMSE is not lower than the real one, it may be stated that the algorithm has a guaranteeing property. We assume suitable to provide such RMSE along with others because it is one of the ways to judge about accuracy in real survey.  An example of simulation for optimistic tuning of the algorithm, that is, such that the algorithm is tuned to a less variable field is shown in Fig. 1. It is a case of motion in a strongly disturbed field (SD of the field is 20 arcsec, correlation radius -5 nmi) with a setting to a weakly disturbed field (SD of the field -5 arcsec, correlation radius -20 nmi).
Analyzing the curves in Fig. 1, we can draw the following conclusions: for the motion in a highly disturbed field with the filter, which use incorrect weakly disturbed model. The estimation accuracy will be approximately 1.5 times worse than the optimal one and 2-3 times worse than the calculated one, which will be shown by the covariance equation of the suboptimal filter. These conclusions are also valid for the case of smoothing algorithms. However, due to the fact that the smoothing mode is more accurate than the filtering mode, the difference in the absolute values of the RMSE in estimating the GDV orizontal components between the suboptimal and optimal algorithms is not so significant here. As can be seen from Fig. 1, for the preset parameters, the guaranteeing property of the filter does not hold either.
The situation that arises is not satisfactory for estimation of the GDV horizontal components, which is why an attempt was made to find such a filter setting that would ensure min- with the filter and smoother tuned to a field with SD of 15 arcsec and correlation radius of 10 nmi imal loss to the optimal algorithm both in weakly disturbed and strongly disturbed fields. Figure 2 presents similar plots for the motion in a strongly disturbed field with a SD of 20 arcsec and a correlation radius of 5 nmi with the filter setting for the field SD of 15 arcsec and a correlation radius of 10 nmi. With this setting, the suboptimal filter insignificantly loses to the optimal one; however, the guaranteeing property of estimation does not hold either. The results of the simulation of motion in a weakly disturbed field (the field SD -5 arcsec, the correlation radius -20 nmi) with the same setting are presented in Fig. 3.
It is obvious that in this case, too, the real RMSE of the error in estimating the GDV horizontal components is close to the optimal one. However, the calculated RMSE significantly differs from both the optimal and real RMSE.  Fig. 3 RMSE of GDV horizontal components for the case of motion in a weakly disturbed field (field SD -5 arcsec, correlation radius -20 nmi) with the filter and smoother tuned to a field with SD of 15 arcsec and the correlation radius of 10 nmi

Conclusion
The accuracy of the problem solution in determining GDV horizontal components by the inertial-geodetic method on a marine moving vessel in the filtering and smoothing modes has been studied. It is shown that the use of the smoothing mode for estimation of GDV horizontal components leads to a significant gain in accuracy in a strongly disturbed field as compared to the filtering mode. However, for a weakly disturbed field, there is practically no gain from smoothing. Sensitivity of the geodetic method algorithms for estimating GDV horizontal components to the accuracy of setting the model of these components was analyzed. The analysis has shown that with optimistic filter tuning (tuning to a field that is less variable than in reality), suboptimal algorithms can significantly lose in accuracy to the optimal algorithm, and they do not provide guaranteeing estimation. At the same time, within the known limits of changes in the correlation radiuses and the RMSE of the field in the survey area, it is possible to choose the model parameters so that the suboptimal algorithm will not differ noticeably in accuracy from the optimal one. However, it is difficult to obtain an adequate calculated characteristic of accuracy. The solution to this problem might be in the development of adaptive algorithms.