A Combined Gravity Compensation Method for INS Using the Simplified Gravity Model and Gravity Database

In recent decades, gravity compensation has become an important way to reduce the position error of an inertial navigation system (INS), especially for a high-precision INS, because of the extensive application of high precision inertial sensors (accelerometers and gyros). This paper first deducts the INS’s solution error considering gravity disturbance and simulates the results. Meanwhile, this paper proposes a combined gravity compensation method using a simplified gravity model and gravity database. This new combined method consists of two steps all together. Step 1 subtracts the normal gravity using a simplified gravity model. Step 2 first obtains the gravity disturbance on the trajectory of the carrier with the help of ELM training based on the measured gravity data (provided by Institute of Geodesy and Geophysics; Chinese Academy of sciences), and then compensates it into the error equations of the INS, considering the gravity disturbance, to further improve the navigation accuracy. The effectiveness and feasibility of this new gravity compensation method for the INS are verified through vehicle tests in two different regions; one is in flat terrain with mild gravity variation and the other is in complex terrain with fierce gravity variation. During 2 h vehicle tests, the positioning accuracy of two tests can improve by 20% and 38% respectively, after the gravity is compensated by the proposed method.


Introduction
Achieving accurate navigation is crucial for the modern motion carrier. Nowadays the most common way to navigate for most modern vehicles and aircrafts is utilizing the global positioning system (GPS). The significant advantage of GPS is its real-time high precision positioning with the help of satellites' signals. However, GPS cannot receive the satellites' signals in some circumstances due to the physical blockages, such as in a tunnel or underwater. The inertial navigation system (INS) utilizes Newton's law to realize autonomous navigation worldwide, and the initialization errors propagate throughout the trajectory with the increase in time. Although the long-term navigation accuracy of free-INS cannot compare with that of GPS, it is still essential during times of loss of GPS. Because INS can work in all weather conditions, has high anti-jamming capacity, and radiates no signal to the outside, INS has been used for many military and civil applications [1]. In recent years, with the emergence of cold atom interferometry, developing high precision INS has become possible. It is predicted that the accelerometer and gyro of this INS can have an accuracy of 10 −8 m/s 2 and 10 −6 deg/h, respectively, which is at least 100 times better than the current INS. The significant improvement of INS, especially on inertial sensors, leaves the gravity compensation as the most important part effecting the positioning accuracy of INS, particularly for rough topological areas [2].
Normally, there are mainly three viable means to realize gravity compensation for INS [3]. Firstly there is the conventional method that uses gravitational gradiometers to gain the gravity disturbance on the trajectory [4,5]. This case is hard to apply to engineering applications because of the high cost of devices. The second method is using gravity models (like WGS84, DQM2000, EGM2008, etc.) to calculate the gravity disturbance and compensate it into the error equations of INS [6,7]. This method is easy to operate but the complexity of these models leads to a long time needed for a solution, which will hinder the real-time navigation of INS. The third method is to obtain the gravity disturbance with the help of the interpolation method based on a gravity database and compensate it into error equations of INS incorporated with gravity disturbance, which has a relatively ideal compensation result when gravity data is available [8][9][10][11][12]. Considering the characteristics of the second and the third method, in this paper, a combined gravity compensation method for INS using the simplified gravity model and gravity data is proposed to restrain the error propagation of INS. At last, two vehicle tests were applied to prove the effectiveness and precision of the proposed gravity compensation method.
The rest of this paper is organized as follows: Section 2 states the error analysis of INS solution, considering gravity disturbance. Section 3 presents the principle of the simplified gravity model. Section 4 introduces the theory and frame of the gravity data-based compensation method using ELM. Section 5 exhibits the framework of this combined gravity compensation method in INS. Section 6 performs vehicle tests in two different regions. Section 7 concludes the whole paper.

Definition of Gravity Disturbance Vector
As Figure 1 shows, the gravity disturbance vector (GDV) is the vector difference between the normal gravity and the actual gravity on the same point in space. GDV is composed of two parts: the orthogonal component is denoted as gravity anomaly and the tangential component is denoted as vertical deflection [13]. It is predicted that the accelerometer and gyro of this INS can have an accuracy of 10 −8 m/s 2 and 10 −6 deg/h, respectively, which is at least 100 times better than the current INS. The significant improvement of INS, especially on inertial sensors, leaves the gravity compensation as the most important part effecting the positioning accuracy of INS, particularly for rough topological areas [2]. Normally, there are mainly three viable means to realize gravity compensation for INS [3]. Firstly there is the conventional method that uses gravitational gradiometers to gain the gravity disturbance on the trajectory [4,5]. This case is hard to apply to engineering applications because of the high cost of devices. The second method is using gravity models (like WGS84, DQM2000, EGM2008, etc.) to calculate the gravity disturbance and compensate it into the error equations of INS [6,7]. This method is easy to operate but the complexity of these models leads to a long time needed for a solution, which will hinder the real-time navigation of INS. The third method is to obtain the gravity disturbance with the help of the interpolation method based on a gravity database and compensate it into error equations of INS incorporated with gravity disturbance, which has a relatively ideal compensation result when gravity data is available [8][9][10][11][12]. Considering the characteristics of the second and the third method, in this paper, a combined gravity compensation method for INS using the simplified gravity model and gravity data is proposed to restrain the error propagation of INS. At last, two vehicle tests were applied to prove the effectiveness and precision of the proposed gravity compensation method.
The rest of this paper is organized as follows: Section 2 states the error analysis of INS solution, considering gravity disturbance. Section 3 presents the principle of the simplified gravity model. Section 4 introduces the theory and frame of the gravity data-based compensation method using ELM. Section 5 exhibits the framework of this combined gravity compensation method in INS. Section 6 performs vehicle tests in two different regions. Section 7 concludes the whole paper.

Definition of Gravity Disturbance Vector
As Figure 1 shows, the gravity disturbance vector (GDV) is the vector difference between the normal gravity and the actual gravity on the same point in space. GDV is composed of two parts: the orthogonal component is denoted as gravity anomaly and the tangential component is denoted as vertical deflection [13]. In Figure 1, gp and γp are defined as the gravity vector and the normal gravity vector at the point P, respectively. The difference between gp and γp is defined as gravity disturbance vector δg: In Figure 1, g p and γ p are defined as the gravity vector and the normal gravity vector at the point P, respectively. The difference between g p and γ p is defined as gravity disturbance vector δg: As Figure 2 shows, the difference in direction is denoted as the deflections of the vertical (DOVs) [14,15]. The DOVs are composed of two components, ζ and η are represented as the north-south component and the east-west component, respectively.
where, ∆g N , ∆g E , and ∆g U are presented as the north, east, and vertical components of gravity disturbance, respectively. γ 0 is the value of normal gravity. As Figure 2 shows, the difference in direction is denoted as the deflections of the vertical (DOVs) [14,15]. The DOVs are composed of two components, ζ and η are represented as the north-south component and the east-west component, respectively.
where, ΔgN, ΔgE, and ΔgU are presented as the north, east, and vertical components of gravity disturbance, respectively. γ0 is the value of normal gravity.

INS Error Equations Incorporated with Gravity Disturbance
The following INS error equations are considered with gravity disturbance [16][17][18][19]: (2 ) g n n n n b n n n b A ie en n n n n ie en in the above formulas,

INS Error Equations Incorporated with Gravity Disturbance
The following INS error equations are considered with gravity disturbance [16][17][18][19]: in the above formulas, δV n is defined as the velocity error in navigation frame, φ n is defined as attitude error of the motion carrier, f n is defined as the specific force in navigation frame, C n b is a direction cosine matrix which can be used to transform the body acceleration vector into the navigation frame, δK A and δA are defined as scale coefficient error and installation angle error of accelerometer respectively, δG and δK G are defined as installation angle error of gyros and scale coefficient error respectively, ∇ n is defined as the bias of accelerometers in navigation frame, δg n is defined as the gravity disturbance in navigation frame and ε n is the drift of gyros in navigation frame, δω n in can be calculated from the following equation: where ω n ie and ω n en are presented as earth's rotation rate and the navigation frame's rotation with respect to earth respectively; both are expressed in the navigation frame. They can be calculated from the following formula: where L, λ, and h are denoted as the latitude, longitude, and altitude of the motion carrier respectively, R M , R N are denoted as meridian radius and prime vertical radius, respectively. From Equations (4)-(9), we could draw the conclusion that INS velocity errors δV n are mainly caused by the accelerometer's output error (δK A + δA) f b + ∇ n , velocity error δV n , and gravity disturbance δg n . And the velocity error δV n is the main error source that causes the position errors of INS.
Based on the above analysis, we can see that gravity disturbance δg n first leads to INS velocity errors in Equation (4), then further affects the position and attitude accuracy of INS through Equations (5) and (6). With the improvement of the accuracy of inertial sensors, INS solution errors due to gravity disturbance are becoming more and more significant and thus cannot be neglected; some measures must be taken to compensate for it.
Here we calculate the north position error caused by three typical values of north-south gravity vertical deflection ζ as an example to better represent the influence of gravity disturbance on INS. And the simulation results are shown in Table 1 and Figure 3. Table 1 and Figure 3 show that the north position error of INS due to gravity vertical deflection propagates according to the Schuler cycle (about 84.4 min), and the amplitude of the position error is in direct proportion to the value of gravity vertical deflection. Thus, the influence of gravity disturbance in a high-precision INS cannot be neglected; it must be compensated to realize a more accurate navigation result.

The Principle of the Simplified Gravity Model
For a simplified reference, the Earth can be assumed to take the shape of an ellipsoid of revolution with a uniform mass distribution. There are some Earth models that can describe the Earth's characteristics, such as WGS84, DQM2000, EGM2008, etc., which are normally used to describe the Earth's gravitational field and are expressed with high degree spherical harmonic polynomials. Normally, in most inertial navigators, the WGS84 reference ellipsoid is used to calculate

The Principle of the Simplified Gravity Model
For a simplified reference, the Earth can be assumed to take the shape of an ellipsoid of revolution with a uniform mass distribution. There are some Earth models that can describe the Earth's characteristics, such as WGS84, DQM2000, EGM2008, etc., which are normally used to describe the Earth's gravitational field and are expressed with high degree spherical harmonic polynomials. Normally, in most inertial navigators, the WGS84 reference ellipsoid is used to calculate the normal gravity and the order of it only expands to 2 degrees. However, with the development of the gravity models, WGS84 cannot perfectly describe the gravitational field to meet the demands of high accurate inertial navigators. To further restrain the errors caused by the gravity model, the more accurate gravity model (EGM2008) is used to replace the WGS84 model to calculate the normal gravity. Although the EGM2008 model can well describe the characteristics of the Earth, it is hard to apply this gravity model with a high number of degrees to the solution of the INS because of the time and space complexity [7]. Thus, here a simplified gravity model is proposed to solve this problem.

Expression of Spherical Harmonic Gravity Model (SHM)
The gravitational potential of the earth to external point is a harmonic function, it can be represented by the infinite order number of a spherical harmonic function [20,21]: where (ρ, θ, λ) is the spherical coordinate with the center of the earth as the origin of the coordinate; f is denoted as the Newton's gravitational constant; a e is the major radius of reference ellipsoid; M is the mass of the Earth; P nm (cos θ) is associated Legendre polynomials, C nm , S nm are denoted as constant coefficients of the spherical harmonic of degree n and order m. Therefore, if these coefficient values are known, the value of the gravitational potentiation can be obtained. In order to get the value of the gravitational potentiation, the number of coefficients is limited. Thus, Equation (10) can be rewritten as follows: (C nm cos(mλ) + S nm sin(mλ)P nm (cos θ))], (11) where N is the highest order of the known coefficients.

The Selection of the Degree for the Simplified Gravity Model
The covariance between the gravity potential and harmonic coefficients considers the degree of their correlation. The effect of certain coefficients on the gravity potential is reflected by a high value of the covariance. Therefore, we could theoretically use it to determine a suitable degree of SHM. In the gravity potential theory, Equation (11) can be changed into the following equation by multiplying both sides with a specific P nm (cos θ) cos(mλ) or P nm (cos θ) sin(mλ) and integrates over the unit sphere [22,23]: where dσ = sin θdθdλ is defined as the surface element of the unit sphere. Defining L C and L S as a linear functional, so the Equation (12) can be simplified as follows: The linear functional is defined as the continuous analogue to the usual concept of a linear function in n-dimensional vector space. Therefore, the covariance between C nm , S nm and T satisfy the following condition: where K(P) is the covariance of the gravity potential on point P, described as where R is the Earth's average radius, ρ is the geocentric radius distance of P, and K n is the degree variance of coefficients: Consequently, the covariance between the gravity potential and the harmonic coefficients can be described as The maximum number of degrees of EGM2008 is 2160 degrees, however, the high order terms have only little influence on the results of calculation. Thus, to better illustrate the number of orders of EGM2008 affecting on the results of calculation, here we set the maximum order (N max = 70) to calculate the absolute magnitudes of cov(C nm , V) P and cov(S nm , V) P . Since the covariance of the points around the world are similar, the point P(θ = 45 • , λ= 45 • ) on the Earth's surface is used for discussion. When N max = 70, the absolute magnitudes of cov(C nm , V) P and cov(S nm , V) P range from 10 −3 to 10 −15 and 10 −3 to 10 −13 , respectively. To better visualize the results, Figures 4 and 5 are plotted with common logarithms of |cov(C nm , V) P | and |cov(S nm , V) P | respectively, shown as follows:

The Principle of ELM.
In recent years, ELM has been becoming a popular research method to solve many engineering problems. ELM has its own distinguishing features, which are depicted as follows [24]: (1) Only the predefined network structure needs to be modulated; (2) ELM has the ability to do fast learning; (3) High generalization performance can be achieved through ELM; (4) A wide selection range of activation functions can be used in ELM. In Figures 4 and 5, horizontal ordinate (m) and longitudinal coordinate (n) represent the number of degrees of the EGM2008. In these two figures, the darker place means the greater value of logarithms of |cov(C nm , V) P | and |cov(S nm , V) P |. From these two figures, we could see that the values of the covariance are mainly focused below 12 degrees, which indicates the value of gravity mainly focuses on the SHM's lower degree. Thus, in this paper, we set the maximum order of the simplified SHM as N max = 12 to calculate the normal gravity.

The Principle of ELM
In recent years, ELM has been becoming a popular research method to solve many engineering problems. ELM has its own distinguishing features, which are depicted as follows [24]: (1) Only the predefined network structure needs to be modulated; (2) ELM has the ability to do fast learning; (3) High generalization performance can be achieved through ELM; (4) A wide selection range of activation functions can be used in ELM.
Due to its many advantages, in this paper, an estimation method based on ELM is used to estimate the gravity disturbance on the trajectory of the motion carrier, then compensate it into the INS solution equations to improve the navigation accuracy.
The principle of the ELM is described as follows [25]: . , x in ] T has n inputs and t i = [t i1 , t i2 , t i3 , . . . , t im ] T has m outputs. Here we define each x i has longitude and latitude two values, described as x i = [λ i L i ] and the output t i is the gravity disturbance on geoid, described as t i = [∆g Ei ∆g Ni ∆g Ui ]. Here l is denoted as the number of hidden neurons, ω is denoted as the l × n input weight matrix where ω j = [ω j1 , ω j2 , ω j3 , . . . , ω jn ] T , b is the l × 1 biases vector and β is the l × m output weight matrix where β j = β j1 , β j2 , β j3 , . . . , β jm T . Generally, the ELM network function is shown as where j ∈ {1, 2, . . . , l}, ω j · x i is defined as the inner product of ω j and x i , the activation function g(x) chooses sigmoid function in the proposed method, which is formulated as Equation (18) can be simplified as the following simple form: where H is denoted as the hidden layer output matrix of the network. T is the output matrix and According to ELM theory, the input weights and hidden biases can be selected at random, and the output weights can be obtained by calculating the following formula: where H + is the Moore-Penrose (MP) generalized inverse of matrix H.

The Procedure of the Data-Based Gravity Disturbance Compensation Method in INS Using ELM
The process of the data-based gravity disturbance compensation method for INS using ELM is described as follows: 1.
Obtain the motion carrier's position. Get the position information (λ, L) of the motion carrier through the INS calculation.

2.
Choose the gravity data base. Find the suitable gravity data base taking the position obtained by step 1 as the center. Normally, 5 × 5 gravity grid data base is chosen as the training database. 3.
ELM training. Set the motion carrier's position information (λ, L) as the inputs of the network and acquire the gravity disturbance on geoid (∆g E0 , ∆g N0 , ∆g U0 ) through training the gravity data base acquired by step 2. 4.
Upward continuation. Calculate the gravity disturbance with upward continuation to the point where INS is. The height of INS is acquired through altimeter. In the geographic engineering application, the most practical upward continuation method is free air correction. The computational equation is described as follows [26]: where ∆g is the gravity disturbance where the motion carrier is and ∆g 0 is the gravity disturbance on the geoid. The units of ∆g and ∆g 0 are milligal. H is the height between the geoid and the motion carrier. The unit of H is meter.

The Framework of the Combined Gravity Compensation Method for INS
The flow chart of the combined gravity compensation method for INS is illustrated in Figure 6: disturbance on the geoid. The units of g and  0 g are milligal. H is the height between the geoid and the motion carrier. The unit of H is meter.

The Framework of the Combined Gravity Compensation Method for INS
The flow chart of the combined gravity compensation method for INS is illustrated in Figure 6:

Position calculation
Attitude array calculation Attitude calculation data ELM Simplified gravity model g Figure 6. The flow chart of the combined gravity compensation method in the inertial navigation system (INS).
The process of the combined gravity compensation method is described as follows: 1. Obtain the motion carrier's position. Get the position information of the motion carrier through the INS calculation. 2. Calculate the gravity and the gravity disturbance. Based on the position information from step 1, the gravity g and the gravity disturbance g  are obtained respectively through the simplified gravity model and ELM training. 3. Compensate the gravity and the gravity disturbance from step 2 into INS equations, considering the gravity disturbance to restrain the error propagation.

Experiment
To prove the effectiveness and feasibility of the combined gravity compensation method in INS, two vehicle tests were applied in the city of Xi'an in China and denoted as test 1 and test 2. Test 1 and test 2 were applied in the areas with flat terrain and complex terrain respectively. On the test vehicle, an altimeter, a GPS receiver, a set of batteries, a high-precision INS, and a data processing computer (PC) were carried. The high-precision INS and the batteries were fixed on the overplate at the trunk. The computer was used to collect the data of GPS and INS in real time, then the data was resolved The process of the combined gravity compensation method is described as follows: 1.
Obtain the motion carrier's position. Get the position information of the motion carrier through the INS calculation.

2.
Calculate the gravity and the gravity disturbance. Based on the position information from step 1, the gravity g and the gravity disturbance δg are obtained respectively through the simplified gravity model and ELM training.

3.
Compensate the gravity and the gravity disturbance from step 2 into INS equations, considering the gravity disturbance to restrain the error propagation.

Experiment
To prove the effectiveness and feasibility of the combined gravity compensation method in INS, two vehicle tests were applied in the city of Xi'an in China and denoted as test 1 and test 2. Test 1 and test 2 were applied in the areas with flat terrain and complex terrain respectively. On the test vehicle, an altimeter, a GPS receiver, a set of batteries, a high-precision INS, and a data processing computer (PC) were carried. The high-precision INS and the batteries were fixed on the overplate at the trunk. The computer was used to collect the data of GPS and INS in real time, then the data was resolved for navigation. The test vehicle is shown in Figure 7. The accuracy of the inertial sensors and GPS are listed in Table 2.
The travel profiles, gravity anomalies and DOVs of the two tests are illustrated in Figures 8 and 9. In Figure 9, the gravity anomalies and DOVs of the two tests were obtained from the estimation based on the measured gravity data using ELM training. The gravity data is provided by the Institute of Geodesy and Geophysics; Chinese Academy of sciences. for navigation. The test vehicle is shown in Figure 7. The accuracy of the inertial sensors and GPS are listed in Table 2. The travel profiles, gravity anomalies and DOVs of the two tests are illustrated in Figures 8 and  9. In Figure 9, the gravity anomalies and DOVs of the two tests were obtained from the estimation based on the measured gravity data using ELM training. The gravity data is provided by the Institute of Geodesy and Geophysics; Chinese Academy of sciences.   for navigation. The test vehicle is shown in Figure 7. The accuracy of the inertial sensors and GPS are listed in Table 2. The travel profiles, gravity anomalies and DOVs of the two tests are illustrated in Figures 8 and  9. In Figure 9, the gravity anomalies and DOVs of the two tests were obtained from the estimation based on the measured gravity data using ELM training. The gravity data is provided by the Institute of Geodesy and Geophysics; Chinese Academy of sciences.   From Figure 9, we can see that the values of DOVs and gravity anomaly in the area of test 1 are both larger than that of test 2. To better prove the effectiveness and accuracy of the combined gravity compensation method, the maximum position errors were used to prove this. Here we regarded the position result of GPS as the true value, and respectively compared it with the position results compensated with three different gravity compensation methods: 1 only with the reference ellipse (shown as blue colour); 2 with the reference ellipse and gravity disturbance (shown as black colour); 3 with the proposed compensation method (shown as red colour).The position error results are shown in Figure 10 According to Figure 10, the maximum position errors of INS compensated with three different methods in two tests are listed in Table 3.
From Figure 10 and Table 3, we could draw the conclusion that the proposed gravity compensation method performed well in two tests, especially in test 2 with a wider variation range of gravity disturbance. During the 2 h vehicle tests, compared with the solution results compensated From Figure 9, we can see that the values of DOVs and gravity anomaly in the area of test 1 are both larger than that of test 2. To better prove the effectiveness and accuracy of the combined gravity compensation method, the maximum position errors were used to prove this. Here we regarded the position result of GPS as the true value, and respectively compared it with the position results compensated with three different gravity compensation methods: 1 only with the reference ellipse (shown as blue colour); 2 with the reference ellipse and gravity disturbance (shown as black colour); 3 with the proposed compensation method (shown as red colour).The position error results are shown in Figure 10 According to Figure 10, the maximum position errors of INS compensated with three different methods in two tests are listed in Table 3.
From Figure 10 and Table 3, we could draw the conclusion that the proposed gravity compensation method performed well in two tests, especially in test 2 with a wider variation range of gravity disturbance. During the 2 h vehicle tests, compared with the solution results compensated with the reference ellipse only, the maximum value of position errors can reduce by 20% and 38% respectively, in test 1 and test 2, with the proposed gravity compensation method. Thus, the effectiveness of this combined gravity compensation method for INS is proved.

Conclusions
This paper proposed a combined gravity compensation method for high-precision INS. The method utilizes the simplified gravity model to subtract the normal gravity from the navigation equations, and meanwhile compensates the gravity disturbance through the ELM training based on measured gravity data. The effectiveness and accuracy of the combined gravity compensation method were verified by vehicle tests, which were applied with a high-precision INS, and the solution results show that the maximum position errors can reduce by 20% and 38%, respectively, in flat terrain and complex terrain. It should be noted that in this work the tests were carried on an experimental vehicle running on the surface of the Earth with relatively low speed. Therefore, in the follow-up study, a flight test with a relatively high travel speed should applied to further test the applicability of this gravity compensation method.