An improved gravity compensation method for high-precision free-INS based on MEC–BP–AdaBoost

A gravity disturbance vector is the vector difference between the actual gravity and the normal gravity at the same point in space (figure 1). It divides into two parts: the tangential component (vertical deflection) and the orthogonal component (gravity anomaly) [19].

Zoom In Zoom Out Reset image size

Suppose the gravity vector g_P and the normal gravity vector γ_P are at the same point P and the gravity disturbance vector δg is defined as their difference:

$$\begin{eqnarray}&&\delta g={{g}_{P}}-{{\gamma}_{P}}.\end{eqnarray} \tag{ 1 }$$

The differences in direction are defined as the deflections of the vertical (DOVs) [19]. The DOVs have two components, a north–south component ζ and an east–west component η (figure 2).

$$\begin{eqnarray}&&\delta g=\left[\Delta {{g}_{\text{E}}} \Delta {{g}_{\text{N}}} \Delta {{g}_{\text{U}}}\right],\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\left\{\begin{array}{*{35}{l}} \delta {{g}_{\text{N}}}=-{{\gamma}_{0}}\zeta \\ \delta {{g}_{\text{E}}}=-{{\gamma}_{0}}\eta \end{array}.\right. \end{eqnarray} \tag{ 3 }$$

Zoom In Zoom Out Reset image size

In equations (2) and (3), Δg_N is the north component of gravity disturbance, Δg_E is the east component of gravity disturbance, Δg_U is the vertical component of gravity disturbance and γ₀ is the value of normal gravity.

The INS error equations when incorporated with gravity disturbance can be defined as follows [20]:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \delta {{\overset{\centerdot}{{V}}\,}^{\text{n}}}=-{{\phi}^{\text{n}}}\times {{f}^{\text{n}}}+C_{\text{b}}^{\text{n}}\left(\delta {{K}_{\text{A}}}+\delta A\right){{f}^{\text{n}}}+\delta {{V}^{\text{n}}}\times \left(2\omega _{\text{ie}}^{\text{n}}+\omega _{\text{en}}^{\text{n}}\right) \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+\,{{V}^{\text{n}}}\times \left(2\delta \omega _{\text{ie}}^{\text{n}}+\delta \omega _{\text{en}}^{\text{n}}\right)+{{\nabla}^{\text{n}}}+\delta {{g}^{\text{n}}}, \end{array}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\left\{\begin{array}{*{35}{l}} \delta \overset{\centerdot}{{L}}\,=\frac{\delta {{V}_{\text{N}}}}{{{R}_{\text{M}}}+h} \\ \delta \overset{\centerdot}{{\lambda}}\,=\frac{\delta {{V}_{\text{E}}}}{{{R}_{\text{N}}}+h}\sec L+\delta L\frac{{{V}_{\text{E}}}}{{{R}_{\text{N}}}+h}\tan L\sec L \\ \delta \overset{\centerdot}{{h}}\,=\delta {{V}_{\text{U}}} \end{array},\right. \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\overset{\centerdot}{{\phi}}\,=\phi \times \omega _{\text{in}}^{\text{n}}+\delta \omega _{\text{in}}^{\text{n}}-C_{\text{b}}^{\text{n}}\left(\delta {{K}_{\text{G}}}+\delta G\right)\omega _{\text{ib}}^{\text{b}}+{{\varepsilon}^{\text{n}}},\end{eqnarray} \tag{ 6 }$$

where ${{V}^{\text{n}}}$ is the velocity in the navigation frame, $\delta {{V}^{\text{n}}}$ is the velocity error in the navigation frame, $\delta {{\overset{\centerdot}{{V}}\,}^{\text{n}}}$ is the differential form of $\delta {{V}^{\text{n}}}$, ${{\phi}^{\text{n}}}$ is the attitude error, ${{f}^{\text{n}}}$ is the specific force expressed in the navigation frame, $C_{\text{b}}^{\text{n}}$ is a direction cosine matrix used to transform the body acceleration vector into the navigation frame, $\delta {{K}_{\text{A}}}$ and $\delta A$ are the scale coefficient error and installation angle error of the accelerometer respectively, ${{\nabla}^{\text{n}}}$ is the accelerometer bias expressed in the navigation frame, $\delta {{g}^{\text{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{\centerdot}{{L}}\,,\delta \overset{\centerdot}{{\lambda}}\,,\delta \overset{\centerdot}{{h}}\,$ are the differential forms of $\delta L,\delta \lambda,\delta h$, $\delta {{V}_{\text{N}}},\delta {{V}_{\text{E}}},\delta {{V}_{\text{U}}}$ are the velocity errors of the north, east, and vertical directions respectively, ${{R}_{\text{M}}}$ and ${{R}_{\text{N}}}$ denote the meridian radius and prime vertical radius respectively, $\delta {{K}_{\text{G}}}$ and $\delta G$ are the scale coefficient error and installation angle error of the gyroscope respectively, and ${{\varepsilon}^{\text{n}}}$ is the gyroscope drift in the navigation frame. $\delta \omega _{\text{in}}^{\text{n}}$ can be expressed as follows:

$$\begin{eqnarray}&&\delta \omega _{\text{in}}^{\text{n}}=\delta \omega _{\text{ie}}^{\text{n}}+\delta \omega _{\text{en}}^{\text{n}},\end{eqnarray} \tag{ 7 }$$

where $\omega _{\text{ie}}^{\text{n}}$ and $\omega _{\text{en}}^{\text{n}}$ are the Earth's rotation rate and the navigation frame's rotation with respect to Earth respectively, and both are expressed in navigation frame. Their computational formulas are defined as follows:

$$\begin{eqnarray}&&\delta \omega _{\text{en}}^{\text{n}}=\left[\begin{array}{c} -\frac{\delta {{V}_{\text{N}}}}{{{R}_{\text{M}}}+h} \\ \frac{\delta {{V}_{\text{E}}}}{{{R}_{\text{N}}}+h} \\ \frac{\delta {{V}_{\text{E}}}\tan L}{{{R}_{\text{N}}}+h}+\delta L\frac{{{V}_{\text{E}}}{{\sec}^{2}}L}{{{R}_{\text{N}}}+h} \end{array}\right],\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}\delta \omega _{\text{ie}}^{\text{n}}={{\left[\begin{array}{c} 0 & -\delta L{{\omega}_{\text{ie}}}\sin L & \delta L{{\omega}_{\text{ie}}}\cos L \end{array}\right]}^{T}}.\end{eqnarray} \tag{ 9 }$$

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

According to the above analysis, gravity disturbance $\delta {{g}^{\text{n}}}$ first causes INS velocity errors in equation (4), which leads to INS position and attitude errors by equations (5) and (6). With the rapid improvement of inertial sensors, errors caused by gravity disturbance are at least as large as errors caused by inertial sensors; therefore the influence of gravity disturbance on INS cannot be ignored, and must be compensated.

To better illustrate the effect of gravity disturbance on INS error propagation, we calculate 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 and figure 3.

Table 1. The effect of gravity vertical deflection $\zeta$ on the north component of gravity disturbance $\Delta {{g}_{\text{N}}}$ and north position error of INS.

	Gravity vertical deflection $\zeta$
	$\zeta =1\,\text{s}$	$\zeta =5\,\text{s}$	$\zeta =10\,\text{s}$	$\zeta =20\,\text{s}$	$\zeta =30\,\text{s}$
$\Delta {{g}_{\text{N}}}$ (mGala)	5	24	48	95	143
North position error (m)	62	309	617	1234	1852

^a1 mGal  =  1  ×  10⁻⁵ m s⁻²

Zoom In Zoom Out Reset image size

Table 1 and figure 3 show that the position error of the north channel caused by gravity disturbance generates periodic variation of the Schuler cycle (about 84.4 min), and the amplitude is proportional to the value of north–south gravity vertical deflection. Therefore, the effect of gravity disturbance for high- precision INS should be effectively measured and compensated in order to achieve a better navigation solution of INS.

Figure 4 illustrates the framework of the MEC–BP–AdaBoost neural-network-based gravity compensation method for INS [13–15]. Figure 4 shows that the proposed gravity compensation method for INS is composed of four modules: (1) training sample selection, (2) BP neural network optimizing selection by MEC, (3) MEC–BP–AdaBoost model training, and (4) gravity disturbance compensation.

Zoom In Zoom Out Reset image size

• (1)  
  Training sample selection. In this paper, measured marine gravity disturbance is chosen as the training data, of which roughly one-third are selected randomly as testing samples and the rest are used as training samples.
• (2)  
  Data normalization. Before the training process, the data should be normalized. Although it is important to apply neural network training on some data, biased results may be generated from such training. The biased results are inherent whenever the data trained has many outliers because neural network training ignores the impact of outliers in data. Therefore, the data must be normalized before neural network training is applied. In this paper, a normalization function 'mapminmax' is adopted. Data are normalized into [0.1,0.9] to guarantee convergence of the network.
• (3)  
  Determine the structure of the BP neural network. According to the theory of Kolmogorov [21], a three-layer network can approach any continuous function with any expected accuracy in theory. So this paper provides a three-layer BP neural network, which includes the input layer, the hidden layer and the output layer. The input layer consists of two neurons, which are the latitude (L) and longitude (λ) of the INS. The output layer consists of three neurons, which are three components of the gravity disturbance, denoted as $\Delta {{g}_{\text{E}}}$, $\Delta {{g}_{\text{N}}}$, $\Delta {{g}_{\text{U}}}$. The structure of the BP neural network is shown in figure 5. To fulfill the requirements on accuracy and training time-saving, the number of hidden layers should be selected reasonably. Each node in the network is used to calculate the inner product of the input vector and weight vector by a nonlinear transfer function to obtain a scalar result; here we choose the tansig function. And the number of single hidden-layer nodes can be calculated using the following equation:

  $$\begin{eqnarray}&&m=\sqrt{n+l}+\alpha,\end{eqnarray} \tag{ 10 }$$

  where $n$ is the number of nodes in the input layer, $l$ is the number of nodes in the output layer, and $\alpha$ is a constant value, which ranges from 1 to 10. According to equation (10), the number of neurons in the single hidden layer network is 3–13. To determine the number of nodes in the hidden layer, we test using three groups of gravity disturbance data under different numbers of hidden layer nodes, and take mean square error (MSE) as the evaluation criterion. The calculation equation of MSE is described as follows:

  $$\begin{eqnarray}&&\text{MSE}=\frac{1}{n}\underset{i=1}{\overset{n}{\sum}}\,{{\left({{y}_{i}}-{{\hat{y}}_{i}}\right)}^{2}},\end{eqnarray} \tag{ 11 }$$

  where n is the number of testing samples, ${{y}_{i}}$ is the true value of the gravity disturbance, and ${{\hat{y}}_{i}}$ is the predicted value of the gravity disturbance. The prediction results are listed in table 2. The minimum MSE of the neural network appears when m  =  10, so the number of the hidden layer nodes is set as 10.
• (4)  
  Group initialization. To obtain optimal performance of the training, MEC initializes the population randomly with five superior groups and five temporary groups. Each group includes 30 individuals bred by the uniform distribution in the solution space.
• (5)  
  Similartaxis and local competition. First, the local optimum will be sought in every group through competing with each other. Second, all individuals are dispersive around the winner, which is obtained by the above step, according to the normal distribution during the new individuals' generation. A new winner will appear by comparing the score of each individual in the new generated group. On comparing the scores of the two best winners, if the old one is higher than the new one we say that the group is mature. Otherwise, the new one will replace the old one, and the new winner's information will be recorded on the local billboard. The process is repeated until the group is stabilized.
• (6)  
  Dissimilation and global competition. The temporary group will be abandoned when its score is lower than the score of any stabilized superior one. Meanwhile, the temporary group that has a higher score than the score of any previous stabilized superior will replace the superior group obtained previously. New temporary groups will be built randomly to replace the abandoned ones. The new temporary groups will also do the similartaxis and local competition. The above process will be repeated until the global optimum is found.
• (7)  
  Interconnection weights and thresholds initialization. After decoding the global optimum, the initialization of interconnection weights and thresholds will be obtained.
• (8)  
  Weak predictor prediction. Initializing the distribution ${{D}_{t}}(i)=1/m$, $i=1,2,\cdots,m$ when the training samples have m sets of data. To determine the number of MEC–BP neural network predictors, we do a test with three groups of gravity disturbance data under different numbers of MEC–BP neural network predictors, and also compare the prediction errors of different numbers of MEC–BP neural network predictors with those of a single MEC–BP neural network predictor. The prediction performance of the test with different numbers of MEC–BP neural networks is listed in table 3. Table 3 shows that there is minimum MSE of MEC–BP–AdaBoost neural networks when there are 10 MEC–BP neural networks. So the number of MEC–BP neural networks is set as 10, and the 10 MEC–BP neural network weak predictors are trained. After the training of t weak predictors, the error sum ${{e}_{t}}$ of prediction serial $g(t)$ is obtained. The computational formula of prediction serial weight ${{a}_{t}}$ is denoted as follows:

  $$\begin{eqnarray}&&{{a}_{t}}=\frac{1}{2}\ln \left(\frac{1-{{e}_{t}}}{{{e}_{t}}}\right).\end{eqnarray} \tag{ 12 }$$
• (9)  
  Adjusting weight. The training sample weight of the next running (${{D}_{t+1}}(i)$) would be adjusted based on the weight ${{a}_{t}}$,

  $$\begin{eqnarray}&&{{D}_{t+1}}(i)=\frac{{{D}_{t}}(i)}{{{B}_{t}}}\times \exp \left[-{{a}_{t}}{{y}_{i}}{{g}_{t}}\left({{x}_{i}}\right)\right]\quad \quad i=1,2,\cdots,m.\end{eqnarray} \tag{ 13 }$$

  In the above equation, ${{B}_{t}}$ is the unitary factor, ${{g}_{t}}\left({{x}_{i}}\right)$ is the output of the network, and ${{y}_{i}}$ is the expected prediction result.
• (10)  
  Strong prediction function. The T groups of weak predictor function $f\left({{g}_{t}},{{a}_{t}}\right)$ will be obtained after T times running and finally combined into a strong prediction function $F(x)$.

  $$\begin{eqnarray}&&F(x)=\frac{{{a}_{t}}}{{\sum}_{t=1}^{T}{{a}_{t}}}f\left({{g}_{t}},{{a}_{t}}\right).\end{eqnarray} \tag{ 14 }$$
• (11)  
  Gravity disturbance estimation. The gravity disturbance $\delta g$($\Delta {{g}_{\text{E} ~ }}, \Delta {{g}_{\text{N}}}, \Delta {{g}_{\text{U}}}$) can be obtained by substituting the position coordinates into prediction function $F(x)$. And the position coordinates are obtained from the INS navigation solution of the previous time.
• (12)  
  Gravity disturbance compensation. Substitute the obtained gravity disturbance $\delta g$($\Delta {{g}_{\text{E}}}, \Delta {{g}_{\text{N}}}, \Delta {{g}_{\text{U}}}$) from step 11 into the INS solution equations to restrain the error propagation.

Zoom In Zoom Out Reset image size