Polar Region Integrated Navigation Method Based on Covariance Transformation

Abstract

Aircraft flying the trans-arctic routes usually apply inertial navigation mechanization in two different navigation frames, e.g., the local geographic frame and the grid frame. However, this change of navigation frame will cause filter overshoot and error discontinuity. To solve this problem, taking the inertial navigation system/global navigation satellite system (INS/GNSS) integrated navigation system as an example, an integrated navigation method based on covariance transformation is proposed. The relationship of the system error state between different navigation frames is deduced as a means to accurately convert the Kalman filter’s covariance matrix. The experiment and semi-physical simulation results show that the presented covariance transformation algorithm can effectively solve the filter overshoot and error discontinuity caused by the change of navigation frame. Compared with non-covariance transformation, the system state error is thereby reduced significantly.

1. Introduction

Considering that the distance of a great circle flight route is shorter, using trans-arctic routes can accomplish great savings in flying time when aircraft make transcontinental flights. Due to the demands of flight safety, each aircraft usually uses an INS/GNSS integrated navigation system to provide high-precision navigation information. The INS/GNSS integrated navigation system has broad development prospects. Previous literature [1] proposed an integrated navigation scheme based on INS and GNSS single-frequency precision point positioning, which is expected to be an advantage for low-cost precise land vehicle navigation applications. Several researchers [2,3] have discussed the application of GNSS/INS on railways. Traditional INS/GNSS-integrated navigation algorithms are based on a north-oriented geographic frame. However, as the latitude increases, the traditional algorithms lose their efficacy in the polar region because of the meridian convergence. To solve this problem, when the aircraft is in the polar region, pilots usually plan their route based on polar-adaptable coordinate frames, such as the Earth-centered Earth-fixed frame (e-frame) [4], transversal Earth frame (t-frame) [5,6], pseudo-Earth frame [7], wander frame [8] and grid frame (G-frame) [9,10].

Although these coordinate frames are adaptable to polar regions, they cannot accomplish successful global navigation individually because some of them have specific mathematical singularities, such as the t-frame, pseudo-Earth frame, wander frame, and G-frame. These coordinate frames are usually adopted only in the polar region, and the local geographic frame (n-frame) is used as the reference navigation frame in non-polar regions. The e-frame can be used for continuous worldwide navigation. However, because the e-frame adopts Cartesian coordinates, the height channel is coupled with three rectangular coordinates but this causes position errors to diverge rapidly and brings difficulties to damping filtering. In addition, the e-frame does not have an explicit azimuth, which is inconvenient for flight route planning. Usually, the INS/GNSS integrated navigation system takes the local geographic frame as the navigation frame at low and middle latitudes and turns instead to grid frames at high latitudes. When the navigation frame is switched between different coordinate frames, such as the G-frame and n-frame, the structure of the filter changes. In this case, as another study [11] points out, if the consistency of the error state estimation cannot be guaranteed, this will cause the integrated navigation filter to overshoot and trigger error discontinuity.

However, the current research [12,13,14] on polar region navigation mainly focuses on the design of an integrated navigation algorithm within the polar region or on looking for a navigation frame to achieve worldwide navigation independently and to avoid the problem caused by switching between navigation frames. One study [15] proposed the virtual sphere n-vector algorithm and derived detailed mechanization and dynamic equations. Their virtual sphere n-vector algorithm used the surface normal vector of the ellipsoid model to represent the aircraft’s position, and did not have specific mathematical singularities. Essentially, the virtual sphere n-vector algorithm is the same as the e-frame algorithm and its azimuth definition is indistinct. The researchers of [11] and [16] proposed a hybrid polar navigation method, which accomplishes the inertial navigation mechanization in the e-frame, whereas it outputs the navigation parameters in the G-frame or t-frame. In addition, the studies of [11,16] introduce a position matrix to decouple the height channel and three rectangular coordinates, which can solve the problem of position error divergence. In this way, the continuity of global navigation is guaranteed. However, it completely changes the navigation frame of the current airborne inertial navigation system, which is not conducive to system upgrades. Papers by [17,18] both proposed indirect polar navigation methods, using a combination of the wander frame and G-frame or the t-frame to achieve smooth switching of navigation frames. However, indirect polar navigation methods did not fundamentally solve the filter consistency problem during the coordinate frames switching.

In order to solve the problem of filter discontinuity caused by the change of navigation frame, this paper proposes a polar-region airborne INS/GNSS integrated navigation method, based on covariance transformation. The transformation relationship between the system error state in the local horizontal geographic frame and that in the grid frame is deduced. Flight experiments at mid-latitudes initially proved the effectiveness of the covariance transformation method. It is difficult to conduct experiments in the polar region. A purely mathematical simulation cannot accurately reflect real aircraft situations [19]. To solve this problem, the authors of [19,20] proposed a virtual polar-region method based on the t-frame or the G-frame. In this way, the experimental data from middle and low latitude regions can be converted to the polar region. Verification by semi-physical simulations, based on the proposed method by [20], is also conducted and gives more convincing results.

This paper is organized as follows. Section 2 describes the grid-based strap-down inertial navigation system (SINS), including the mechanization and dynamic model of the grid SINS. In Section 3, the covariance transformation method is presented. In addition, Section 3 also offers a navigation frame-switching method based on the INS/GNSS integrated navigation method. Section 4 verifies the effectiveness of the proposed method through experimentation and semi-physical simulation. Finally, general conclusions are discussed in Section 5.

2. The Grid SINS

2.1. Grid Frame and Grid SINS Mechanization

The definition of the grid reference frame is shown in Figure 1. The grid plane is parallel to the Greenwich meridian, and its intersection with the tangent plane at the position of the aircraft is the grid’s north. The angle between geographic north and grid north provides the grid angle, and its clockwise direction is the positive direction. The up direction of the grid frame is the same as that of the local geographic frame and forms an orthogonal right-handed frame with the orientations at grid east and grid north.

The grid angle $\sigma$ is expressed as found in [9]:

$$\begin{array}{l}\sin \sigma =\frac{\sin L\sin \lambda }{\sqrt{1-{\cos }^{2}L{\sin }^2\lambda }}\\ \cos \sigma =\frac{\cos \lambda }{\sqrt{1-{\cos }^{2}L{\sin }^2\lambda }}\end{array}$$

(1)

The direction cosine matrix ${\mathit{C}}_e^G$ between the G-frame and the e-frame (earth frame) is as found in [9]:

$${\mathit{C}}_e^G={\mathit{C}}_n^G{\mathit{C}}_e^n$$

(2)

where n refers to the local horizontal geographic frame. ${\mathit{C}}_e^n$ and ${\mathit{C}}_n^G$ are expressed as:

$${\mathit{C}}_e^n=\left[\begin{array}{ccc}-\sin \lambda & \cos \lambda & 0\\ -\sin L\cos \lambda & -\sin L\sin \lambda & \cos L\\ \cos L\cos \lambda & \cos L\sin \lambda & \sin L\end{array}\right]$$

(3)

$${\mathit{C}}_n^G=\left[\begin{array}{ccc}\cos \sigma & -\sin \sigma & 0\\ \sin \sigma & \cos \sigma & 0\\ 0& 0& 1\end{array}\right]$$

(4)

The updated equations of the attitude, the velocity, and the position in the grid frame are expressed as:

$${\dot{\mathit{C}}}_b^G={\mathit{C}}_b^G\left[{\mathit{\omega }}_{ib}^b\times \right]-\left[{\mathit{\omega }}_{iG}^G\times \right]{\mathit{C}}_b^G$$

(5)

$${\dot{\mathit{v}}}^G={\mathit{C}}_b^G{\mathit{f}}^b-\left(2{\mathit{\omega }}_{ie}^G+{\mathit{\omega }}_{eG}^G\right)\times {\mathit{v}}^G+{\mathit{g}}^G$$

(6)

$${\dot{\mathit{R}}}^e={\mathit{C}}_G^e{\mathit{v}}^G$$

(7)

where ${\mathit{\omega }}_{iG}^G$ is the turn rate of the G-frame with respect to the i-frame.

$$\begin{gathered} {\mathit{\omega }}_{iG}^G={\mathit{\omega }}_{ie}^G+{\mathit{\omega }}_{eG}^G={\mathit{C}}_e^G{\mathit{\omega }}_{ie}^e+{\mathit{\omega }}_{eG}^G \\ {\mathit{\omega }}_{ie}^G=\left[\begin{array}{c}-{\omega }_{ie}\sin \sigma \cos L\\ {\omega }_{ie}\cos \sigma \cos L\\ {\omega }_{ie}\sin L\end{array}\right],{\mathit{\omega }}_{eG}^G=\left[\begin{array}{cc}\frac{1}{{\tau }_f}& -\frac{1}{R_y}\\ \frac{1}{R_x}& -\frac{1}{{\tau }_f}\\ \frac{\kappa }{{\tau }_f}& -\frac{\kappa }{R_y}\end{array}\right]\left[\begin{array}{c}{v}_E^G\\ v_N^G\end{array}\right] \end{gathered}$$

(8)

where $R_x$ is the radius of curvature of the grid east, $R_y$ is the radius of curvature of the grid north, and ${\tau }_f$ is the distorted radius.

Since the meridian converges rapidly in the polar region, the position of the aircraft in the polar region is usually expressed in the ECEF frame. The relationship between the coordinates $x,y,z$ and the latitude $L$ and the longitude $\lambda$ is given by:

$$\{\begin{array}{c}x=\left(R_N+h\right)\cos L\cos \lambda \\ y=\left(R_N+h\right)\cos L\sin \lambda \\ z=\left[R_N{\left(1-f\right)}^2+h\right]\sin L\end{array}$$

(9)

2.2. Dynamic Model of the Grid SINS

The mechanization of the grid SINS is accomplished in Section 2.1. Next, the Kalman filter, based on the G-frame, needs to be designed. In order to design the Kalman filter, the dynamic model of the G-frame, including three differential equations, is given below, as put forward in [10].

The attitude error is defined as:

$${\tilde{\mathit{C}}}_b^G=\left[\mathit{{\rm I}}-{\mathit{\varphi }}_{}^G\times \right]{\mathit{C}}_b^G$$

(10)

$$\left[{\mathit{\varphi }}_{}^G\times \right]=-{\tilde{\mathit{C}}}_b^G{\mathit{C}}_G^b$$

(11)

where ${\tilde{\mathit{C}}}_b^G$ is the estimated attitude, expressed in terms of the direction cosine matrix.

Differentiating Equation (11) gives:

$$\left[{\dot{\mathit{\varphi }}}_{}^G\times \right]=-{\dot{\tilde{\mathit{C}}}}_b^G{\mathit{C}}_G^b-{\tilde{\mathit{C}}}_b^G{\left({\dot{\mathit{C}}}_b^G\right)}^T$$

(12)

Substituting ${\dot{\tilde{\mathit{C}}}}_b^G$ and ${\dot{\mathit{C}}}_b^G$ from Equation (5) gives:

$$\begin{array}{cc}\hfill \left[{\dot{\mathit{\varphi }}}_{}^G\times \right]& =-{\tilde{\mathit{C}}}_b^G\left[{\tilde{\mathit{\omega }}}_{ib}^b\times \right]{\mathit{C}}_G^b+\left[{\tilde{\mathit{\omega }}}_{iG}^G\times \right]{\tilde{\mathit{C}}}_b^G{\mathit{C}}_G^b+{\tilde{\mathit{C}}}_b^G\left[{\mathit{\omega }}_{ib}^b\times \right]{\mathit{C}}_G^b-{\tilde{\mathit{C}}}_b^G{\mathit{C}}_G^b\left[{\mathit{\omega }}_{iG}^G\times \right]\hfill \\ & =-{\tilde{\mathit{C}}}_b^G\left[\delta {\mathit{\omega }}_{ib}^b\times \right]{\mathit{C}}_G^b+\left[{\tilde{\mathit{\omega }}}_{iG}^G\times \right]{\tilde{\mathit{C}}}_b^G{\mathit{C}}_G^b-{\tilde{\mathit{C}}}_b^G{\mathit{C}}_G^b\left[{\mathit{\omega }}_{iG}^G\times \right]\hfill \end{array}$$

(13)

Substituting ${\tilde{\mathit{C}}}_b^G$ from Equation (10) gives:

$$\begin{array}{cc}\hfill \left[{\dot{\mathit{\varphi }}}_{}^G\times \right]& =-\left[\mathit{{\rm I}}-{\mathit{\varphi }}_{}^G\times \right]{\mathit{C}}_b^G\left[\delta {\mathit{\omega }}_{ib}^b\times \right]{\mathit{C}}_G^b+\left[{\tilde{\mathit{\omega }}}_{iG}^G\times \right]\left[\mathit{{\rm I}}-{\mathit{\varphi }}_{}^G\times \right]-\left[\mathit{{\rm I}}-{\mathit{\varphi }}_{}^G\times \right]\left[{\mathit{\omega }}_{iG}^G\times \right]\hfill \\ & =-{\mathit{C}}_b^G\left[\delta {\mathit{\omega }}_{ib}^b\times \right]{\mathit{C}}_G^b+\left[\delta {\mathit{\omega }}_{iG}^G\times \right]-\left[\delta {\mathit{\omega }}_{iG}^G\times \right]\left[{\mathit{\varphi }}_{}^G\times \right]+\left[{\mathit{\varphi }}_{}^G\times \right]\left[\delta {\mathit{\omega }}_{iG}^G\times \right]\hfill \end{array}$$

(14)

According to Equation (12), the attitude error equation is expressed by:

$${\dot{\mathit{\varphi }}}_{}^G=-{\mathit{\omega }}_{iG}^G\times {\mathit{\varphi }}_{}^G+\delta {\mathit{\omega }}_{iG}^G-{\mathit{C}}_b^G\delta {\mathit{\omega }}_{ib}^b$$

(15)

The velocity error is defined as:

$$\delta {\mathit{v}}_{}^G={\tilde{\mathit{v}}}_{}^G-{\mathit{v}}_{}^G$$

(16)

According to Equation (6), the velocity error equation can be written as:

$$\begin{array}{cc}\hfill \delta {\dot{\mathit{v}}}_{}^G& ={\tilde{\mathit{C}}}_b^G{\tilde{\mathit{f}}}^b-\left(2{\tilde{\mathit{\omega }}}_{ie}^G+{\tilde{\mathit{\omega }}}_{eG}^G\right)\times {\tilde{\mathit{v}}}^G+{\tilde{\mathit{g}}}^G-{\mathit{C}}_b^G{\mathit{f}}^b+\left(2{\mathit{\omega }}_{ie}^G+{\mathit{\omega }}_{eG}^G\right)\times {\mathit{v}}^G-{\mathit{g}}^G\hfill \\ & =\left({\tilde{\mathit{C}}}_b^G-{\mathit{C}}_b^G\right){\mathit{f}}^b+{\mathit{C}}_b^G\delta {\mathit{f}}^b-\left(2{\tilde{\mathit{\omega }}}_{ie}^G+{\tilde{\mathit{\omega }}}_{eG}^G\right)\times \delta {\mathit{v}}^G-\left(2\delta {\mathit{\omega }}_{ie}^G+\delta {\mathit{\omega }}_{eG}^G\right)\times {\mathit{v}}^G-\delta {\mathit{g}}^G\hfill \\ & ={\mathit{f}}^G\times {\mathit{\varphi }}_{}^G+{\mathit{v}}^G\times (2\delta {\mathit{\omega }}_{ie}^G+\delta {\mathit{\omega }}_{eG}^G)-(2{\mathit{\omega }}_{ie}^G+{\mathit{\omega }}_{eG}^G)\times \delta {\mathit{v}}_{}^G+{\mathit{C}}_b^G\delta {\mathit{f}}^b\hfill \end{array}$$

(17)

Substituting ${\tilde{\mathit{C}}}_b^G$ from Equation (10) and ignoring the error of gravity vector gives:

$$\delta {\dot{\mathit{v}}}_{}^G={\mathit{f}}^G\times {\mathit{\varphi }}_{}^G+{\mathit{v}}^G\times (2\delta {\mathit{\omega }}_{ie}^G+\delta {\mathit{\omega }}_{eG}^G)-(2{\mathit{\omega }}_{ie}^G+{\mathit{\omega }}_{eG}^G)\times \delta {\mathit{v}}_{}^G+{\mathit{C}}_b^G\delta {\mathit{f}}^b$$

(18)

From Equation (7), the position error equation is as follows:

$$\delta {\dot{\mathit{R}}}^e={\mathit{C}}_G^e\delta {\mathit{v}}^G+\delta {\mathit{C}}_G^e{\mathit{v}}^G$$

(19)

where:

$$\delta {\mathit{C}}_G^e=\delta {\mathit{C}}_n^G{\mathit{C}}_e^n+{\mathit{C}}_n^G\delta {\mathit{C}}_e^n$$

(20)

According to Equation (2), $\delta {\mathit{C}}_n^G$ and $\delta {\mathit{C}}_e^n$ can be written as:

$$\delta {\mathit{C}}_e^n=\left[\begin{array}{ccc}-\cos \lambda ·\delta \lambda & -\sin \lambda ·\delta \lambda & 0\\ -\cos L\cos \lambda ·\delta L+\sin L\sin \lambda ·\delta \lambda & -\cos L\sin \lambda ·\delta L-\sin L\cos \lambda ·\delta \lambda & -\sin L·\delta L\\ -\sin L\cos \lambda ·\delta L-\cos L\sin \lambda ·\delta \lambda & -\sin L\sin \lambda ·\delta L+\cos L\cos \lambda ·\delta \lambda & \cos L·\delta \lambda \end{array}\right]$$

(21)

$$\delta {\mathit{C}}_n^G=\left[\begin{array}{ccc}-\sin \sigma & -\cos \sigma & 0\\ \cos \sigma & -\sin \sigma & 0\\ 0& 0& 1\end{array}\right]\delta \sigma$$

(22)

where $\delta \sigma$ is the grid angle error, and its dynamic equation can be obtained by differentiating Equation (1):

$$\delta \sigma =\frac{\sin \sigma \cos \sigma \cos L}{\sin L}\delta L+\frac{1-{\cos }^2\sigma {\cos }^{2}L}{\sin L}\delta \lambda$$

(23)

3. Design of an INS/GNSS Integrated Navigation Filter Model with Covariance Transformation

When an aircraft flies in the polar region, it is necessary to change navigation frames from the n-frame to G-frame, and vice versa. In addition to the transformation of navigation parameters, the integrated navigation filter also needs to transform. The Kalman filter includes the state equation and the observation equation, and its update process includes a prediction update and measure update, respectively. The Kalman filter acts to update the error state and its covariance. Different Kalman filters, designed on different navigation frames, have different filter states $\mathit{x}$ and covariance matrices $\mathit{P}$, which need to be transformed.

The filtering state at low and middle latitudes is usually expressed by:

$${\mathit{x}}^n(t)={[{\varphi }_E^n,{\varphi }_N^n,{\varphi }_U^n,\delta v_E^n,\delta v_N^n,\delta v_U^n,\delta L,\delta \lambda ,\delta h,{\epsilon }_x^b,{\epsilon }_y^b,{\epsilon }_z^b,{\nabla }_x^b,{\nabla }_y^b,{\nabla }_{\mathrm{z}}^b]}^T$$

(24)

At high latitudes, the integrated filter is designed in the grid frame. The filtering state is usually expressed by:

$${\mathit{x}}^G(t)={[{\varphi }_E^G,{\varphi }_N^G,{\varphi }_U^G,\delta v_E^G,\delta v_N^G,\delta v_U^G,\delta x,\delta y,\delta z,{\epsilon }_x^b,{\epsilon }_y^b,{\epsilon }_z^b,{\nabla }_x^b,{\nabla }_y^b,{\nabla }_{\mathrm{z}}^b]}^T$$

(25)

Then, the transformation relationship of the filtering state and the covariance matrix need to be deduced. Comparing (24) and (25), it can be seen that the states that remain unchanged before and after the navigation frame change are the gyroscope bias ${\mathsf{\epsilon }}^b$ and the accelerometer bias ${\nabla }^b$. Therefore, it is only necessary to establish a transformation relationship between the attitude error $\mathit{\varphi }$, the velocity error $\delta \mathit{v}$, and the position error $\delta \mathit{p}$.

The transformation relationship between the attitude error ${\mathit{\varphi }}_{}^n$ and ${\mathit{\varphi }}_{}^G$ is determined as follows.

According to the definition of $\delta {\mathit{C}}_b^G$:

$$\delta {\mathit{C}}_b^G=-[{\mathit{\varphi }}_{}^G\times ]{\mathit{C}}_b^G$$

(26)

From the equation, ${\mathit{C}}_b^G={\mathit{C}}_n^G{\mathit{C}}_b^n$, $\delta {\mathit{C}}_b^G$ can be expressed as:

$$\delta {\mathit{C}}_b^G=\delta {\mathit{C}}_n^G{\mathit{C}}_b^n+{\mathit{C}}_n^G\delta {\mathit{C}}_b^n=-[{\mathit{\varphi }}_{nG}^G\times ]{\mathit{C}}_n^G{\mathit{C}}_b^n-{\mathit{C}}_n^G[{\mathit{\varphi }}_{}^n\times ]{\mathit{C}}_b^n$$

(27)

Substituting $\delta {\mathit{C}}_b^G$ from Equation (26), ${\mathit{\varphi }}_{}^G$ can be described as:

$${\mathit{\varphi }}_{}^G={\mathit{C}}_n^G{\mathit{\varphi }}_{}^n+{\mathit{\varphi }}_{nG}^G$$

(28)

where ${\mathit{\varphi }}_{nG}^G$ is the error angle vector of ${\mathit{C}}_n^G$:

$$\begin{gathered} \delta {\mathit{C}}_n^G={\tilde{\mathit{C}}}_n^G-{\mathit{C}}_n^G=-\left[{\mathit{\varphi }}_{nG}^G\times \right]{\mathit{C}}_n^G \\ {\mathit{\varphi }}_{nG}^G={\left[\begin{array}{ccc}0& 0& -\delta \sigma \end{array}\right]}^T \end{gathered}$$

(29)

The transformation relationship between the velocity error $\delta {\mathit{v}}_{}^n$ and $\delta {\mathit{v}}_{}^G$ is determined as follows:

$$\delta {\mathit{v}}_{}^G={\mathit{C}}_n^G\delta {\mathit{v}}_{}^n+\delta {\mathit{C}}_n^G{\mathit{v}}^n={\mathit{C}}_n^G\delta {\mathit{v}}_{}^n-[{\mathit{\varphi }}_{nG}^G\times ]{\mathit{C}}_n^G{\mathit{v}}^n$$

(30)

From Equation (9), the position error can be written as:

$$\left[\begin{array}{c}\delta x\\ \delta y\\ \delta z\end{array}\right]=[\begin{array}{c}-\left(R_N+h\right)\sin L\cos \lambda \\ -\left(R_N+h\right)\sin L\sin \lambda \\ \left[R_N{\left(1-f\right)}^2+h\right]\cos L\end{array}\begin{array}{cc}-\left(R_N+h\right)\cos L\sin \lambda & \mathrm{c}\mathrm{os}L\cos \lambda \\ \left(R_N+h\right)\cos L\cos \lambda & \cos L\sin \lambda \\ 0& \sin L\end{array}]\left[\begin{array}{c}\delta L\\ \delta \lambda \\ \delta h\end{array}\right]$$

(31)

To sum up, the transformation relationship between the system error state ${\mathit{x}}^n(t)$ and ${\mathit{x}}^G(t)$ is as follows:

$${\mathit{x}}^G(t)=\mathit{\Phi }{\mathit{x}}^n(t)$$

(32)

where $\mathit{\Phi }$ is determined by Equations (28)–(31), and is given by:

$$\begin{gathered} \mathit{\Phi }=\left[\begin{array}{ccccc}{\mathit{C}}_n^G& {\mathit{O}}_{3\times 3}& {\mathit{\Phi }}_a& {\mathit{O}}_{3\times 3}& {\mathit{O}}_{3\times 3}\\ {\mathit{O}}_{3\times 3}& {\mathit{C}}_n^G& {\mathit{\Phi }}_b& {\mathit{O}}_{3\times 3}& {\mathit{O}}_{3\times 3}\\ {\mathit{O}}_{3\times 3}& {\mathit{O}}_{3\times 3}& {\mathit{\Phi }}_c& {\mathit{O}}_{3\times 3}& {\mathit{O}}_{3\times 3}\\ {\mathit{O}}_{3\times 3}& {\mathit{O}}_{3\times 3}& {\mathit{O}}_{3\times 3}& {\mathit{I}}_{3\times 3}& {\mathit{O}}_{3\times 3}\\ {\mathit{O}}_{3\times 3}& {\mathit{O}}_{3\times 3}& {\mathit{O}}_{3\times 3}& {\mathit{O}}_{3\times 3}& {\mathit{I}}_{3\times 3}\end{array}\right] \\ {\mathit{\Phi }}_a=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ \frac{\cos L\sin \sigma cos\sigma }{\sin L}& \frac{1-{\cos }^{2}L{\cos }^2\sigma }{\sin L}& 0\end{array}\right] \\ {\mathit{\Phi }}_b=\left[\begin{array}{ccc}0& -v_U^G& v_N^G\\ v_U^G& 0& -v_E^G\\ -v_N^G& v_E^G& 0\end{array}\right]{\mathit{\Phi }}_a \\ {\mathit{\Phi }}_c=[\begin{array}{c}-\left(R_N+h\right)\sin L\cos \lambda \\ -\left(R_N+h\right)\sin L\sin \lambda \\ \left[R_N{\left(1-f\right)}^2+h\right]\cos L\end{array}\begin{array}{cc}-\left(R_N+h\right)\cos L\sin \lambda & \mathrm{c}\mathrm{os}L\cos \lambda \\ \left(R_N+h\right)\cos L\cos \lambda & \cos L\sin \lambda \\ 0& \sin L\end{array}] \end{gathered}$$

(33)

The transformation relation of the covariance matrix is as follows:

$$\begin{array}{cc}\hfill {\mathit{P}}^G\left(t\right)& =E\left\{\left({\tilde{\mathit{x}}}^G\left(t\right)-{\mathit{x}}^G\left(t\right)\right){\left({\tilde{\mathit{x}}}^G\left(t\right)-{\mathit{x}}^G\left(t\right)\right)}^T\right\}\hfill \\ & =E\left\{\mathit{\Phi }\left({\tilde{\mathit{x}}}^n\left(t\right)-{\mathit{x}}^n\left(t\right)\right){\left({\tilde{\mathit{x}}}^n\left(t\right)-{\mathit{x}}^n\left(t\right)\right)}^T{\mathit{\Phi }}^T\right\}\hfill \\ & =\mathit{\Phi }E\left\{\left({\tilde{\mathit{x}}}^n\left(t\right)-{\mathit{x}}^n\left(t\right)\right){\left({\tilde{\mathit{x}}}^n\left(t\right)-{\mathit{x}}^n\left(t\right)\right)}^T\right\}{\mathit{\Phi }}^T\hfill \\ & =\mathit{\Phi }{\mathit{P}}^n\left(t\right){\mathit{\Phi }}^T\hfill \end{array}$$

(34)

Once the aircraft flies out of the polar region, ${\mathit{x}}^G$ and ${\mathit{P}}^G$ should be converted to ${\mathit{x}}^n$ and ${\mathit{P}}^n$, which can be described as:

$$\begin{array}{l}{\mathit{x}}^n(t)={\mathit{\Phi }}^{-1}{\mathit{x}}^G(t)\\ {\mathit{P}}^n\left(t\right)={\mathit{\Phi }}^{-1}{\mathit{P}}^G\left(t\right){\mathit{\Phi }}^{-T}\end{array}$$

(35)

The process of the covariance transformation method is shown in Figure 2. At middle and low latitudes, the system accomplishes the inertial navigation mechanization in the n-frame. When the aircraft enters the polar regions, the system accomplishes the inertial navigation mechanization in the G-frame. Correspondingly, the navigation parameters are output in the G-frame. During the navigation parameter conversion, the navigation results and Kalman filter parameter can be established according to the proposed method.

4. Experimental Results

This section presents the experiment and simulation results. A high-precision ring laser gyro inertial navigation system is used to conduct flight experiments. The ring laser gyro bias stability is less than 0.01°/h. The accelerometer bias stability is less than 20 μg. The GNSS positioning error is less than 10 m, which is used as the position reference. The update frequency of the gyro and accelerometer is 200 Hz, while the update frequency of GNSS is 1 Hz.

4.1. Flight Experiment

Flight experiments were carried out at middle latitudes. The duration was 4 h, including a half-hour’s alignment. Firstly, the experiment was carried out based on the n-frame. After flying for half an hour, the navigation frame was changed to the G-frame until the end of the flight.

In order to evaluate the effects caused by navigation frame transformation, the results that are based on a local horizontal geographic frame are used as reference. Because the flight experiments are conducted at middle latitudes, there is no algorithm error introduced by the choice of navigation frame. The reference result is reliable. The navigation results, based on the covariance transformation and non-covariance transformation, are shown in Figure 3.

As shown in Figure 3a, the change of the filter structure results in a fluctuation of the relative attitude error. Among the attitude errors, the relative yaw error reaches the maximum value of 2.2 ‘without covariance transformation. As a comparison, the covariance transformation method reduces this error to 0.3’. Figure 3b shows that the relative position error is less than 10 m, regardless of whether the covariance transformation is used. The INS/GNSS-integrated navigation filter uses the position information provided by GNSS as observations to correct the INS results, resulting in less position error. As shown in Figure 3c,d, the maximum bias error of the gyroscope with and without covariance transformation reaches 0.003°/h and 0.01°/h, respectively. The maximum bias error of the accelerometer with and without covariance transformation reaches 6 μg and 50 μg, respectively. Because of the non-zero values of the non-diagonal elements in the covariance matrix, the bias estimates of the gyroscope and accelerometer are affected by the cross-coupling of other error states. As a result, the bias estimates of the gyroscope and accelerometer also show instability.

The flight experiment was repeated six times. The results of the experiments are shown in

Table 1. The relative error, based on the covariance transformation in six experiments.

Experiment Number	Attitude Error 1/′	Position Error/m	Accelerometer Bias Estimation Error 2/μg	Gyro Bias Estimation Error 2/(°/h)
1	0.67	0.6	9.48	0.002
2	0.64	0.58	9.24	0.0023
3	0.93	0.17	1.10	0.0022
4	0.61	0.47	9.26	0.0003
5	0.66	0.44	9.36	0.0003
6	0.26	0.17	1.12	0.0003
average	0.628333	0.405	6.59	0.0012

¹ The attitude error refers to the maximum attitude error. ² The bias estimation error refers to the largest of the three gyros or accelerometers.

and

Table 2. The relative error, based on the non-covariance transformation in six experiments.

Experiment Number	Attitude Error/′	Position Error/m	Accelerometer Bias Estimation Error/μg	Gyro Bias Estimation Error/(°/h)
1	3.34	2.4	59.3	0.0091
2	2.89	2.09	62.0	0.0094
3	5.56	1.26	20.0	0.0215
4	4.45	2.47	62.5	0.0069
5	2.56	1.76	61.7	0.0038
6	5.56	0.89	22.8	0.0019
average	4.06	1.811667	48.1	0.0088

.

To sum up, when the navigation frame changes directly, the integrated navigation results show severe fluctuation, taking more than an hour to reach stability again. The lower the observability of the error state, the larger the error amplitude. The integrated navigation results, based on the covariance transformation method, do not fluctuate during the change of the navigation frame, which is consistent with the reference results. Experimental results confirm the effectiveness of the proposed algorithm.

4.2. Semi-Physical Simulation Experiment

Pure mathematical simulation is difficult to use to accurately simulate an actual situation. Thus, a virtual polar-region method is used to convert the measured aviation data to 80° latitude, to obtain semi-physical simulation data [20]. In this way, the reliability of the algorithm at high latitudes can be verified. In this simulation, the navigation result based on the G-frame is used as a reference, which will avoid the decrease of algorithm accuracy caused by the rise in latitude. The simulation results, based on the covariance transformation and non-covariance transformation, are shown in Figure 4.

As can be seen in Figure 4a, among the attitude errors, the relative yaw error is the largest. The relative yaw error reaches 5 ‘without covariance transformation. The integrated navigation result with covariance transformation has a less relative yaw error of 0.2’. As shown in Figure 4b, the relative position error is 12 m, without covariance transformation. The integrated navigation result with covariance transformation shows better stability and a smaller relative position error of 8 m. As shown in Figure 4c,d, the maximum bias error of the gyroscope with and without covariance transformation reached 0.001°/h and 0.02°/h, respectively. The maximum bias error of the accelerometer, with and without covariance transformation, reached 0.1 μg and 25 μg, respectively.

5. Conclusions

The advantage of the covariance transformation method is that it establishes the transformation relationship of the integrated navigation filter between the n-frame and G-frame. It fundamentally solves the problem of filter overshoot and error discontinuity. caused by the change in navigation frame, improving navigation accuracy when crossing the polar region. Besides this, the covariance transformation method does not change the existing navigation algorithm. The results of the flight experiment and semi-physical simulation show that the covariance transformation method is effective at any latitude. As the latitude increases, the horizontal component of the earth’s rotational angular velocity decreases, resulting in weaker observability of the yaw angle error. The error fluctuation caused by the frame switching will increase. In this case, the covariance transformation method still makes a smooth transformation via the integrated navigation filter.