1. Introduction {#sec1-sensors-16-02019}
===============

Precise navigation is an essential factor of modern carriers. Nowadays, most modern vehicles and aircraft depend on the Global Positioning System (GPS) for position update as they navigate. GPS utilizes the signals from navigation satellites to realize high precision locating. However, normally, GPS cannot work when the satellite signals are not available due to physical blockage, such as inside a cave or under water. The inertial navigation system utilizes the laws of Newtonian physics to realize autonomous navigation worldwide in all weather conditions. Due to INS's inherent character, it can overcome the disadvantages of GPS and widely serve in all branches of military and many civil applications \[[@B1-sensors-16-02019],[@B2-sensors-16-02019]\].

With the rapid improvement of the inertial instruments (accelerometers and gyros), gravity compensation in INS has already become an important way to further improve the navigation accuracy. Estimating the gravity field for navigation with INS is known as gravity compensation. In other words, it uses the position provided by INS to estimate the gravity disturbance vector δg \[[@B3-sensors-16-02019]\].

To achieve high accuracy for INS navigation, there are mainly three feasible ways for gravity compensation \[[@B4-sensors-16-02019]\]. First, the traditional method obtains the gravity disturbance on the trajectory using a gravitational gradiometer that senses the gradient of the potential field \[[@B5-sensors-16-02019]\]. This case depends on augmenting the INS with a sensitive instrument that has a long history of technology development, but very little operational experience. The second method is to use existing gravity field models (like DQM2000, GFZ97, EGM2008, etc.) and perform compensation on the location results directly by calculating the position errors with the gravity disturbance given \[[@B6-sensors-16-02019]\]. However, the accuracy of these models cannot satisfy the requirements of high precision INS, especially in rough topography, such as mountains, plateaus and oceanic trenches \[[@B7-sensors-16-02019]\]. The third way is to obtain the gravity disturbance using the interpolation method based on measured gravity data on the geoid, then to process the gravity disturbance to the height where INS has an upward continuation \[[@B4-sensors-16-02019],[@B8-sensors-16-02019]\]. In recent years, the most popular interpolation methods used in the geodetic and geophysical communities have been the inverse distance weighted (IDW) interpolation method and the bilinear interpolation method \[[@B9-sensors-16-02019],[@B10-sensors-16-02019]\]. A relationship exists between the position point, and the value of gravity data is different in different areas: linearity in a plain area; nonlinearity in a rough topological area. The above two interpolation methods can achieve a high accuracy when applied in a plain area, yet have bad performance for areas where fierce gravity variation exists, especially in rugged mountain regions. In this paper, a novel gravity disturbance estimation method that uses extreme learning machine (ELM) is proposed to solve this problem. ELM has extensive applications because of its simple structure and mature technology. It provides a new approach for the nonlinear approximation. This motivates this paper to apply ELM to estimate the gravity disturbance on the trajectory \[[@B11-sensors-16-02019]\]. The ELM-based gravity disturbance estimation algorithm is utilized in the training process to establish the prediction model, with the carrier position (longitude and latitude) that the INS provided as input and the gravity disturbance on the geoid as output, then processes the obtained gravity disturbance to the height (provided by the altimeter) where INS has an upward continuation \[[@B12-sensors-16-02019],[@B13-sensors-16-02019]\]. Finally, the estimated gravity disturbance on the trajectory is compensated in the INS error equations incorporated with gravity disturbance to restrain the error propagation in INS. In this paper, numerical tests and field experiments were carried out to verify the accuracy and effectiveness of the proposed gravity compensation method.

The paper is organized as follows: in [Section 2](#sec2-sensors-16-02019){ref-type="sec"}, the error analysis of the INS solution incorporated with gravity disturbance is proposed; in [Section 3](#sec3-sensors-16-02019){ref-type="sec"}, a brief review of the artificial neural network (ANN) is introduced; in [Section 4](#sec4-sensors-16-02019){ref-type="sec"}, the theory and framework of the ELM-based gravity disturbance compensation method are proposed; in [Section 5](#sec5-sensors-16-02019){ref-type="sec"}, the numerical tests are designed to prove the accuracy and superiority of ELM-based gravity disturbance estimation method; in [Section 6](#sec6-sensors-16-02019){ref-type="sec"}, field experiments in a city area and a mountain area are presented. Finally, [Section 7](#sec7-sensors-16-02019){ref-type="sec"} concludes the paper.

2. Error Analysis of INS Solution Considering Gravity Disturbance {#sec2-sensors-16-02019}
=================================================================

2.1. Definition of Gravity Disturbance Vector {#sec2dot1-sensors-16-02019}
---------------------------------------------

The gravity disturbance vector is the vector difference between the actual gravity and the normal gravity on the same point in space ([Figure 1](#sensors-16-02019-f001){ref-type="fig"}). It is divided into two parts: the tangential component (vertical deflection) and the orthogonal component (gravity anomaly) \[[@B14-sensors-16-02019]\].

In [Figure 1](#sensors-16-02019-f001){ref-type="fig"}, n is the plumb line on the geoid and perpendicular to the geoid, n' is the ellipsoidal normal on the reference ellipsoid and perpendicular to the reference ellipsoid. Suppose the gravity vector *g~p~* and the normal gravity vector *γ~p~* are at the same point *P*. The gravity disturbance vector *δg* is defined as their difference: $$\delta g = g_{P} - \gamma_{P}$$ $g_{P}$ and *γ~p~* are different in value and direction. The angle between the projection of $g_{P}$ on the north plane and *γ~p~* is defined as ζ, and the angle between the projection of $g_{P}$ on the east plane and *γ~p~* is defined as η. ζ and η together are defined as deflections of the vertical (DOVs), as shown in [Figure 2](#sensors-16-02019-f002){ref-type="fig"} \[[@B14-sensors-16-02019]\].

$$\delta g = \lbrack\Delta g_{E}\Delta g_{N}\Delta g_{U}\rbrack$$

The gravity disturbance $\delta g$ can be decomposed into Cartesian coordinate system: in the east direction, it is $\Delta g_{E}$; in the north direction, it is $\Delta g_{N}$; and in the vertical direction, it is $\Delta g_{U}$. In the gravity database, the gravity anomaly ($\Delta g_{U}$) and DOVs ($\zeta,\eta$) are given. $\Delta g_{N}$ and $\Delta g_{E}$ are calculated by the following equation: $$\left\{ \begin{array}{l}
{\Delta g_{N} = - \gamma_{0}\zeta} \\
{\Delta g_{E} = - \gamma_{0}\eta} \\
\end{array} \right.$$

In Equation (3), $\gamma_{0}$ is the value of normal gravity.

2.2. INS Error Equations Considering Gravity Disturbance {#sec2dot2-sensors-16-02019}
--------------------------------------------------------

The INS error equations when incorporated with gravity disturbance can be defined as follows \[[@B15-sensors-16-02019],[@B16-sensors-16-02019]\]: $$\begin{array}{l}
{\delta{\overset{˙}{V}}^{n} = - \phi^{n} \times f^{n} + C_{b}^{n}(\delta K_{A} + \delta A)f^{n} + \delta V^{n} \times (2\omega_{ie}^{n} + \omega_{en}^{n})} \\
{+ V^{n} \times (2\delta\omega_{ie}^{n} + \delta\omega_{en}^{n}) + \nabla^{n} + \delta g^{n}} \\
\end{array}$$ $$\left\{ \begin{array}{l}
{\delta\overset{˙}{L} = \frac{\delta V_{N}}{R_{M} + h}} \\
{\delta\overset{˙}{\lambda} = \frac{\delta V_{E}}{R_{N} + h}\sec~L + \delta L\frac{V_{E}}{R_{N} + h}\tan~L~\sec~L} \\
{\delta\overset{˙}{h} = \delta V_{U}} \\
\end{array} \right.$$ $$\overset{˙}{\phi} = \phi \times \omega_{in}^{n} + \delta\omega_{in}^{n} - C_{b}^{n}(\delta K_{G} + \delta G)\omega_{ib}^{b} + \varepsilon^{n}$$ where $V^{n}$ is the velocity in the navigation frame, $\delta V^{n}$ is the velocity error in the navigation frame, $\delta{\overset{˙}{V}}^{n}$ is the differential form of $\delta V^{n}$, $\phi^{n}$ is the attitude error, $f^{n}$ is the specific force expressed in the navigation frame, $C_{b}^{n}$ is a direction cosine matrix used to transform the body acceleration vector into the navigation frame, $\delta K_{A}$ and $\delta A$ are the scale coefficient error and the installation angle error of the accelerometer, respectively, $\nabla^{n}$ is the accelerometer bias expressed in the navigation frame, $\delta g^{n}$ is the gravity disturbance in the navigation frame, $L,\lambda,h$ are the current latitude, longitude and altitude of the body, $\delta L,\delta\lambda,\delta h$ denote the latitude error, longitude error and altitude error, respectively, $\delta\overset{˙}{L},\delta\overset{˙}{\lambda},\delta\overset{˙}{h}$ are the differential forms of $\delta L,\delta\lambda,\delta h$, $\delta V_{N},\delta V_{E},\delta V_{U}$ are the velocity errors of the north, east and vertical directions, respectively, $R_{M}$ and $R_{N}$ denote the meridian radius and prime vertical radius, respectively, $\delta K_{G}$ and $\delta G$ are the scale coefficient error and installation angle error of the gyro, respectively, and $\varepsilon^{n}$ is the gyro drift in the navigation frame; $\delta\omega_{in}^{n}$ can be expressed as follows: $$\delta\omega_{in}^{n} = \delta\omega_{ie}^{n} + \delta\omega_{en}^{n}$$ where $\omega_{ie}^{n}$ and $\omega_{en}^{n}$ are the Earth's rotation rate and the navigation frame's rotation with respect to the Earth, respectively; both are expressed in the navigation frame. Their computational formulas are defined as follows: $$\delta\omega_{en}^{n} = \begin{bmatrix}
{- \frac{\delta V_{N}}{R_{M} + h}} \\
\frac{\delta V_{E}}{R_{N} + h} \\
{\frac{\delta V_{E}\tan~L}{R_{N} + h} + \delta L\frac{V_{E}\sec^{2}~L}{R_{N} + h}} \\
\end{bmatrix}$$ $$\delta\omega_{ie}^{n} = {\lbrack 0~~ - \delta L\omega_{ie}\sin~L~~\delta L\omega_{ie}\cos~L\rbrack}^{T}$$

From Equation (4) to Equation (9), we can see that the accelerometer's output error $(\delta K_{A} + \delta A)f^{n} + \nabla^{n}$, velocity error $\delta V^{n}$ and gravity disturbance $\delta g^{n}$ are the main error sources causing INS velocity error $\delta V^{n}$. The position errors of INS are mainly caused by velocity error $\delta V^{n}$ and the coupling error between velocity and position.

According to the above analysis, gravity disturbance $\delta g^{n}$ firstly influences INS velocity accuracy through Equation (4) and then affects INS position and attitude accuracy by Equations (5) and (6). With the high-precision inertial sensors being available, errors due to gravity disturbance are equal to or even more than the errors caused by inertial sensors, so the influence of gravity disturbance on INS cannot be ignored and must be compensated to further improve the navigation accuracy of INS.

To better illustrate the effect of gravity disturbance on INS error propagation, we calculated INS north position errors caused by different values of north-south gravity vertical deflection $\zeta$ as an example using Equation (5); the results are shown in [Table 1](#sensors-16-02019-t001){ref-type="table"} and [Figure 3](#sensors-16-02019-f003){ref-type="fig"}.

From [Table 1](#sensors-16-02019-t001){ref-type="table"} and [Figure 3](#sensors-16-02019-f003){ref-type="fig"}, we can draw the conclusion that the position error of the north channel caused by gravity disturbance presents periodic variation on the Schuler cycle (about 84.4 min), and the amplitude is proportional to the value of the horizontal component of gravity disturbance. Therefore, the effect of gravity disturbance for high-precision INS cannot be ignored, and effective measures must be taken to compensate it for achieving a better navigation solution of INS.

3. Brief Review of Artificial Neural Networks {#sec3-sensors-16-02019}
=============================================

Nowadays, artificial neural networks (ANNs) have been widely used in most areas. ANNs are computational models that simulate the process of the human brain, and the structure of an artificial neuron is much simpler than a biological neuron. A neural network includes a large number of inter-connected neurons in the input and hidden layers that join to all neurons in the output layer (which is called a full synapse network). Each neuron can do some simple computation, such as summation, subtraction and multiplication. [Figure 4](#sensors-16-02019-f004){ref-type="fig"} presents the structure of an artificial neural network \[[@B17-sensors-16-02019],[@B18-sensors-16-02019],[@B19-sensors-16-02019]\].

In [Figure 4](#sensors-16-02019-f004){ref-type="fig"}, $X_{1}\cdots X_{N}$ and $O_{1}\cdots O_{L}$ are the input and output of the artificial neural network, respectively; $\alpha_{1}\cdots\alpha_{N}$ are the neuron nodes of the input layer; $\theta_{1}\cdots\theta_{q}$ are the neuron nodes of the hidden layer; $\beta_{1}\cdots\beta_{L}$ are the neuron nodes of the output layer; $\omega_{ij}$ is the connecting weight between the input layer and hidden layer; and $\omega_{ki}$ is the connecting weight between the hidden layer and output layer.

ANN is a nonlinear statistical data modeling or decision-making method, so it can be used to model complex relationships between the input and output or to find patterns in data. ANN has the abilities of being data-driven, having self-learning and self-adaption; meanwhile, it has strong capabilities of anti-jamming. Therefore, the neural network-based gravity disturbance estimation method is proposed in this paper.

4. The Theory and Framework of the ELM-Based Gravity Disturbance Compensation Method in INS {#sec4-sensors-16-02019}
===========================================================================================

In recent years, ELM has become a hot research topic for machine learning and artificial intelligence. Compared with other traditional artificial learning methods, ELM has its own advantages, which are described as follows \[[@B11-sensors-16-02019],[@B12-sensors-16-02019]\]: (1)No parameters need to be tuned except the predefined network structure;(2)ELM is capable of faster learning, and most trainings can be completed quickly;(3)ELM can achieve a high generalization performance;(4)ELM has a wide selection range of activation functions that are all piecewise continuous functions that can be used as activation functions.

Because of ELM's many advantages, in this paper, we use ELM to estimate the gravity disturbance on the trajectory based on measured gravity on the geoid and then process the gravity disturbance on the geoid to the height where INS has an upward continuation. Finally, the estimated gravity disturbance on the trajectory is compensated in the INS error equations to restrain the error propagation in INS.

4.1. Extreme Learning Machine {#sec4dot1-sensors-16-02019}
-----------------------------

The mathematical model of the ELM is described as below \[[@B20-sensors-16-02019],[@B21-sensors-16-02019]\]:

Input and output: Consider that we have a training set $\left\{ {(x_{i},t_{i})} \right\}_{i = 1}^{N}$ with N distinct examples, where $x_{i} = \left\lbrack {x_{i1},x_{i2},x_{i3},...,x_{in}} \right\rbrack^{T}$ has n inputs and $t_{i} = \left\lbrack {t_{i1},t_{i2},t_{i3},...,t_{im}} \right\rbrack^{T}$ has m outputs. Here, we define each input $x_{i}$ to be composed of longitude and latitude, described as $x_{i} = \lbrack\lambda_{i}~L_{i}\rbrack$, and the output $t_{i}$ is the gravity disturbance on the geoid, described as $t_{i} = \lbrack\Delta g_{Ei}\quad\Delta g_{Ni}\quad\Delta g_{Ui}\rbrack$.

Symbols in the network: Assume that $l$ is the number of hidden neurons, ω is the $l~ \times ~n$ input weight matrix connecting the $i - {th}$ hidden neuron and the input neurons where $\omega_{j} = {\lbrack\omega_{j1},\omega_{j2},\omega_{j3},...,\omega_{jn}\rbrack}^{T}$, $b_{j}$ is the bias for each single hidden neuron, $b$ is the $l~ \times ~1$ biases vector for hidden neurons, $\beta$ is the $l~ \times ~m$ output weight matrix connecting the $i - {th}$ hidden neuron and the output neuron where $\beta_{j} = \left\lbrack {\beta_{j1},\beta_{j2},\beta_{j3},...,\beta_{jm}} \right\rbrack^{T}$. Generally, the ELM network function is shown as \[[@B11-sensors-16-02019]\]: $$t_{i} = {\sum\limits_{j = 1}^{l}\beta_{j}}g\left( {\omega_{j} \cdot x_{i} + b_{j}} \right),~i = 1,2,...,N;$$ where $j \in \left\{ {1,2,...,l} \right\},~\omega_{j} \cdot x_{i}$ denotes the inner product of $\omega_{j}$ and $x_{i}$; the activation function *g*(*x*) is the sigmoid function in the proposed method, which is formulated as: $$g\left( x \right) = \frac{1}{1 + exp\left\lbrack {- \left( {\omega \cdot x + b} \right)} \right\rbrack}$$

Equation (10) can be simplified as the following simple form: $$H\beta = T$$ where: $$H = \begin{bmatrix}
{g(\omega_{1} \cdot x_{1} + b_{1})} & \cdots & {g(\omega_{l} \cdot x_{1} + b_{l})} \\
 \vdots & \ddots & \vdots \\
{g(\omega_{1} \cdot x_{N} + b_{1})} & \cdots & {g(\omega_{l} \cdot x_{N} + b_{l})} \\
\end{bmatrix}_{N \times l}$$

*H* is defined as the hidden layer output matrix of the network. *T* is the output matrix, and $T = {\lbrack t_{1},t_{2},...,t_{N}\rbrack}^{T}$.

In the ELM algorithm, the input weights and hidden biases are chosen at random, and the output weights are calculated through the equation below: $$\beta = H^{+}T$$ where *H*^+^ is the Moore--Penrose (MP) generalized inverse of matrix *H*.

The learning procedure of ELM can be summarized as the four steps in Algorithm 1. **Algorithm 1.** Extreme learning machine (ELM).1. Given a training set with *N* distinct examples $\mathsf{\Psi} = \left\{ {\left. \left( {x_{i},t_{i}} \right) \right|x_{i} \in R^{n},t_{i} \in R^{m}} \right\}_{i = 1}^{N}$, activation function *g*(*x*) and hidden neuron number $l$;2. Set input weights *ω* and hidden biases *b* from \[--1,1\] at random;3. Calculate the hidden layer output matrix *H* using matrix multiplication;4. Calculate the output weights $\beta = H^{+}T$ according to the Moore--Penrose generalized inverse.

4.2. The Framework of the ELM-Based Gravity Compensation Method {#sec4dot2-sensors-16-02019}
---------------------------------------------------------------

The framework of the ELM-based gravity disturbance compensation method is shown in [Figure 5](#sensors-16-02019-f005){ref-type="fig"} \[[@B22-sensors-16-02019],[@B23-sensors-16-02019]\].

The process of the ELM-based gravity disturbance compensation method is described as follows: Get the INS position. Obtain the position value of the INS, longitude (λ) and latitude (L), through the calculation of INS.Choose the gravity database. Search the adaptive gravity database (provided by Institute of Geodesy and Geophysics, Chinese Academy of Sciences) according to the position obtained by Step 1. The gravity disturbance is related to the correlation distances. This means too wide a training area will not improve the estimation result, and too small a training area will not include enough information for the estimation of the gravity disturbance. After many trainings, we found that the size of the training area set as 5′ × 5′ would have the best estimation result. Therefore, here, the gravity database is set as 5′ × 5′ and takes the position of the INS as the central point.ELM training: Set λ, L as the inputs of the ELM algorithm, and obtain the gravity disturbance on the geoid ($\Delta g_{E0},\Delta g_{N0},\Delta g_{U0}$) through the training with the gravity database obtained by Step 2.Upward continuation: Process the gravity disturbance with upward continuation to the height where the INS is. The height of the INS is obtained by the altimeter. In geographic engineering applications, the most practical upward continuation method is free air correction. The computational formula is described as follows \[[@B24-sensors-16-02019]\]: $$\Delta g = \Delta g_{0} - 0.3086H$$ where $\Delta g_{0}$ is the value of the gravity disturbance on the geoid, $\Delta g$ is the value of the gravity disturbance where the INS is and *H* is the height value between the geoid and INS.Compensate the gravity disturbance calculated by Step 4 in the INS error equations to restrain the error propagation caused by gravity disturbance.

5. Numerical Test {#sec5-sensors-16-02019}
=================

To test the accuracy of the ELM-based gravity disturbance estimation method, two areas each with 50 gravity data points derived from the Institute of Geodesy and Geophysics, Chinese Academy of Sciences, are chosen as the application regions. Additionally, the resolution of the gravity database is 1′ × 1′. Region 1 is in the Huabei plain area in China where there is a mild gravity variation; Region 2 is in the Qinling mountain area in China where there is a rough gravity variation. The gravity anomaly maps in the two application regions are shown in [Figure 6](#sensors-16-02019-f006){ref-type="fig"} and [Figure 7](#sensors-16-02019-f007){ref-type="fig"}. In each region, the gravity points are divided into two parts: 10 points are used as testing data; the other 40 points as training data.

The evaluation criterions chosen to evaluate the performance of the estimation methods are mean absolute error (MAE), mean radial error (MRE) and root mean square error (RMSE). MAE can reflect the possible estimation error range of estimation methods; MRE stands for the ensemble estimation error; and RMSE reflects the sensitivity and extreme value effect of the estimation results. The calculation formulae of these three evaluation criteria are described as follows \[[@B25-sensors-16-02019],[@B26-sensors-16-02019]\]: $$MAE = \frac{1}{n}{\sum\limits_{i = 1}^{n}{ABS(X_{m,i} - X_{e,i})}}$$ $$MRE = \frac{1}{n}{\sum\limits_{i = 1}^{n}{\lbrack ABS(X_{m,i} - X_{e,i})/X_{m,i}\rbrack}}$$ $$RMSE = \sqrt{\frac{1}{n}{\sum\limits_{i = 1}^{n}{(X_{m,i} - X_{e,i})}^{2}}}$$ where *n* is the number of the testing points and $X_{m,i}$, $X_{e,i}$ are the measured and estimated gravity disturbance values, respectively. $ABS$ is the symbol of the absolute value function.

The performances of the different estimation methods applied in two test regions are listed in [Table 2](#sensors-16-02019-t002){ref-type="table"} and [Table 3](#sensors-16-02019-t003){ref-type="table"}.

From [Table 2](#sensors-16-02019-t002){ref-type="table"} and [Table 3](#sensors-16-02019-t003){ref-type="table"}, we can see that the ELM-based estimation method had a better performance than the other two methods, especially in the mountain area with rough topography. The RMSE of the ELM estimation method improved by 23% and 44% compared with the bilinear interpolation in the plain and mountain area, respectively. Additionally, for the other two variables of gravity disturbance, we had similar estimation results as the gravity anomaly. The estimation performances of the other two variables with the ELM method were better than the two traditional methods.

6. Experiment {#sec6-sensors-16-02019}
=============

To further validate the proposed gravity compensation method, two field experiments were carried out in Shanxi province in China: one was in an urban road with a mild gravity variation, and the other one was in a mountain road with a fierce gravity variation. On the test car were carried a GPS receiver and the INS, which contained the inertial measurement unit (IMU) and the processing computer system (PSC). The PSC collects the GPS data and the IMU data for data storage, then calculates it for navigation. The altimeter and battery were also carried on the test car for measurement and electricity supply. The test car is shown in [Figure 8](#sensors-16-02019-f008){ref-type="fig"}. The performances of the inertial sensors and GPS are listed in [Table 4](#sensors-16-02019-t004){ref-type="table"}. The resolution of the gravity database is 1′ × 1′ (provided by the Institute of Geodesy and Geophysics, Chinese Academy of Sciences).

Field Test 1 was carried out in Xi'an city area in China, and field Test 2 was carried out in Qinling mountain area in China. The travel profiles, gravity anomalies and DOVs on the two trajectories are shown in [Figure 9](#sensors-16-02019-f009){ref-type="fig"} and [Figure 10](#sensors-16-02019-f010){ref-type="fig"}.

The position results without gravity compensation and results with the proposed gravity compensation method are compared with the GPS result, as shown in [Figure 11](#sensors-16-02019-f011){ref-type="fig"}. In [Figure 11](#sensors-16-02019-f011){ref-type="fig"}a--c are the position results of Test 1, [Figure 11](#sensors-16-02019-f011){ref-type="fig"}d--f are the position results of Test 2.

Compared with the position results of GPS, the maximum position errors of two tests are listed in [Table 5](#sensors-16-02019-t005){ref-type="table"}.

The results in [Figure 11](#sensors-16-02019-f011){ref-type="fig"} and [Table 5](#sensors-16-02019-t005){ref-type="table"} show that the proposed gravity compensation method can both improve the positioning accuracy in the two tests, especially in Test 2 with the rough topography, because in Test 2 along the trajectory, the value of the gravity disturbance is much bigger and has a wider variation range. According to the characteristics of the error propagation in free-INS, the error caused by the gravity disturbance increases with time. This can be verified in [Figure 11](#sensors-16-02019-f011){ref-type="fig"}. The errors shown in [Figure 11](#sensors-16-02019-f011){ref-type="fig"} are similar for both methods (without and with the gravity compensation) at the beginning of the experiment, but they diverge in the later stages. During the 2-h field experiment, the maximum value of position error improved by 13% and 29%, respectively, in Tests 1 and 2, when the navigation scheme was compensated by the proposed gravity compensation method. Therefore, the effectiveness of the proposed gravity compensation method is verified in the field experiments.

7. Conclusions {#sec7-sensors-16-02019}
==============

In this paper, a novel gravity compensation method for high-precision INS is proposed, which uses the ELM-based estimation method to obtain the optimal gravity disturbance on the track based on measured gravity data on the geoid, then processes gravity disturbance on the geoid to the height where INS has an upward continuation. Finally, the estimated gravity disturbance on the trajectory is compensated into the INS error equations to restrain the error propagation in INS. The accuracy of the proposed estimation method was evaluated using the numerical tests, which divided the test area into two parts, one was in the plain area and the other in the mountain area. The estimation results showed that the ELM-based estimation method had a better performance than the other two methods, especially in the mountain area with rough topography. The RMSE of the ELM-based estimation method improved by 23% and 44% compared with the bilinear interpolation in the plain and mountain area, respectively. To further verify the performance of the proposed gravity compensation method, two field experiments were carried out with a high-precision INS, and the experiment results showed that the positioning accuracy improved by 13% and 29%, respectively, in the plain area and mountain area with the proposed gravity compensation method. It should be noted that in this work, the experiment was conducted only in cars at a relatively low speed. In the future, a flight experiment with relative high travel speed should be carried out to further verify the efficiency of the proposed gravity compensation method.

The project is supported by the National Natural Science Foundation of China (61340044) and the Fundamental Research Funds for the Central Universities (YWF-10-01-B30). The writers would like to thank the Institute of Geodesy and Geophysics, Chinese Academy of Sciences, for providing the gravity data. Many thanks are also given to Chai at the Institute of Geodesy and Geophysics, Chinese Academy of Sciences, for his valuable comments.

In this paper, Xiao Zhou proposed the gravity compensation method for high precision INS based on extreme learning machine and did the numerical tests. Gongliu Yang analyzed the feasibility of the proposed method and modified the English. Qingzhong Cai provided the source code of inertial navigation algorithm. Jing Wang participated in the experiment and did the simulation.

The authors declare no conflict of interest.

![Description of the gravity disturbance vector.](sensors-16-02019-g001){#sensors-16-02019-f001}

![The deflection of the vertical.](sensors-16-02019-g002){#sensors-16-02019-f002}

![The north position error caused by gravity disturbance.](sensors-16-02019-g003){#sensors-16-02019-f003}

![Structure of an artificial neural network.](sensors-16-02019-g004){#sensors-16-02019-f004}

![The framework of ELM-based gravity compensation method.](sensors-16-02019-g005){#sensors-16-02019-f005}

![The gravity anomaly map in the Huabei plain area.](sensors-16-02019-g006){#sensors-16-02019-f006}

![The gravity anomaly map in the Qinling mountain area.](sensors-16-02019-g007){#sensors-16-02019-f007}

![Field test device. PSC, processing computer system.](sensors-16-02019-g008){#sensors-16-02019-f008}

![Field test trajectories on Google map.](sensors-16-02019-g009){#sensors-16-02019-f009}

![The gravity anomaly and deflections of the vertical (DOVs) of the two tests.](sensors-16-02019-g010){#sensors-16-02019-f010}

###### 

Position errors of the two tests.

![](sensors-16-02019-g011a)

![](sensors-16-02019-g011b)

sensors-16-02019-t001_Table 1

###### 

The north component of gravity and the amplitude of the north position error corresponding to the gravity vertical deflection.

                              $\mathbf{\mathbf{\zeta}}$ = 2 s   $\mathbf{\mathbf{\zeta}}$ = 8 s   $\mathbf{\mathbf{\zeta}}$ = 15 s   $\mathbf{\mathbf{\zeta}}$ = 25 s   $\mathbf{\mathbf{\zeta}}$ = 35 s
  --------------------------- --------------------------------- --------------------------------- ---------------------------------- ---------------------------------- ----------------------------------
  $\Delta g_{N}$ (m Gal \*)   9                                 38                                71                                 118                                166
  North position error (m)    117                               494                               923                                1534                               2160

\* 1 m Gal = 1 × 10^−5^ m/s^2^.

sensors-16-02019-t002_Table 2

###### 

Performance of the estimation methods in test Region 1. MRE, mean radial error; IDW, inverse distance weighted method.

  Evaluation Criterion   Estimation Methods           
  ---------------------- -------------------- ------- -------
  MAE                    0.157                0.138   0.098
  MRE                    0.059                0.058   0.032
  RMSE                   0.285                0.279   0.213

sensors-16-02019-t003_Table 3

###### 

Performance of the estimation methods in test Region 2.

  Evaluation Criterion   Estimation Methods           
  ---------------------- -------------------- ------- -------
  MAE                    0.228                0.203   0.128
  MRE                    0.076                0.056   0.041
  RMSE                   0.367                0.314   0.193

sensors-16-02019-t004_Table 4

###### 

Performance of the sensors for the experiment.

  Sensors Types            Characteristics     Magnitude (1 σ)
  ------------------------ ------------------- -----------------
  Gyroscope                Constant bias       0.003°/h
  Accelerometer            Constant bias       10 μg
  GPS velocity             Horizontal error    0.03 m/s
  Height error             0.05 m/s            
  GPS position             Horizontal error    2 m
  Height error             5 m                 
  Altimeter                Measurement error   ±5 m
  Measurement resolution   0.1 m               

sensors-16-02019-t005_Table 5

###### 

The maximum value of position error (m) compared with the GPS result.

           Without Gravity Compensation   With Gravity Compensation   Position Improvement
  -------- ------------------------------ --------------------------- ----------------------
  Test 1   1050                           913                         137 (13%)
  Test 2   1120                           790                         330 (29%)
