Analysis and Discussion on Ellipsoidal and virtual sphere models used for transverse arrangement INS

Abstract

Transverse arrangement is one of the main methods used in the polar inertial navigation system (INS). In the traditional algorithm, the calculation of using the earth ellipsoid model is complex, while using the earth sphere model cannot satisfy a high-accuracy application. Therefore, an approach based on the virtual sphere model is proposed, which has been proved in simulation experiments to reduce the computational complexity and maintain the same accuracy as the ellipsoid algorithm, but its accuracy has not yet been proved in theory. Starting from the basic principles of the ellipsoid and virtual sphere model algorithm, this paper compares the key formulations of the two. Finally, it is proved that the two arrangements are actually the same.

1. Introduction

Polar navigation has been a hot topic in recent years, owing to the Arctic routes development, the exploitation of resources and the military forces in the polar region. The inertial navigation system (INS) is a passive autonomous navigation system that is not disturbed by the external environment and has become the most important navigation technology in the polar region because of the complex electromagnetic environment (Bian et al., 2020). Due to special factors, such as the convergence of longitude and the close collinearity of the vector between the earth rotation angular rate and gravity, the INS under the traditional arrangement faces problems such as the sharp increase of heading error and calculation overflow. Therefore, polar navigation methods such as grid navigation and transverse navigation need to be used.

The application of the transverse method has a long history. In 1958, the US ‘Nautilus’ nuclear submarine equipped with the N6A INS successfully crossed the Arctic (Curtis and Slater, 1959), which used the transverse navigation method. Initially, the earth sphere model was used in the INS arrangement, which led to an approximation error, so it was considered not suitable for long period navigation (Lin et al., 2019). After several decades of technology development, the transverse navigation method based on a proposed ellipsoid model solved this problem (Xu and Dou, 2014; Yao et al., 2016). It is worth noting that the ellipsoid transverse algorithm is much more complex than the traditional algorithm, but the complexity of these algorithms in the spherical model is basically the same, which will be discussed in detail in Section 2.

To optimise the algorithm and reduce the complexity of the solution, Qin et al. (2018) proposed the transverse algorithm with the virtual sphere model. In simulations, the solving accuracy is commensurate with the transverse ellipsoid arrangement, and the computational complexity is much lower than the latter. Many scholars have also carried out subsequent research, such as constructing integrated navigation with the algorithm (Ge et al., 2021) or proposing new navigation methods inspired by the algorithm (Liu et al., 2020).

It is considered that the simulation results are inadequate to support the confirmation and promotion of the virtual sphere model algorithm. Therefore, after analysing the essence of the ellipsoid algorithm and the virtual sphere model algorithm, this paper compares the calculation formulae of key variables in the two arrangements and proves that the two arrangements are the same. This finding is expected to promote the application and development of inertial navigation methods in the future.

2. Analysis of ellipsodial transverse navigation algorithms

2.1 Transverse latitude and transverse longitude

The latitude and longitude of the Earth used in navigation are measured from geodesy. The latitude of any point on the ellipsoid surface is defined as the angle between the normal vector on the ellipsoid and the equatorial plane, and the longitude is defined as the angle between the projection of the normal vector on the ellipsoid in the equatorial plane and the axis, which connects the geocentre and the prime meridian in the equatorial plane.

Similarly, the definition of transverse longitude and latitude is shown in Figure 1.

Figure 1. Definition of transverse latitude and transverse longitude.

The Earth ellipsoid model equation can be expressed as

(1)\begin{equation}F(x,y,z) = \frac{{{x^2}}}{{{a^2}}} + \frac{{{y^2}}}{{{a^2}}} + \frac{{{z^2}}}{{{b^2}}} - 1 = 0\end{equation}

and the normal vector on M (x ₀, y ₀, z ₀), i.e. O'M, is

(2)\begin{equation} \boldsymbol{n} = ({F_x},{F_y},{F_z}){|_M} = \left( {{x_0},{y_0},\frac{{{z_0}}}{{1 - {e^2}}}} \right)\end{equation}

where e denotes the eccentricity of the ellipsoid, $e = {{\sqrt {{a^2} - {b^2}} } / a}$.

The transverse latitude at M is the angle between the normal vector n and the X-O-Z plane, so the point with the same transverse latitude ${\varphi ^t}$ satisfies Equation (3), which is an elliptical cone with an opening to the Y-axis:

(3)\begin{equation}{F_1}(x,y,z) = \sin {\varphi ^t} - \frac{y}{{\sqrt {{x^2} + {y^2} + {{\left( {\frac{z}{{1 - {e^2}}}} \right)}^2}} }} = 0\end{equation}

Another way to define the transverse latitude is called auxiliary transverse latitude ${\bar{\varphi }^t}$:

(4)\begin{equation}F^\prime_{1} (x,y,z) = \sin {\bar{\varphi }^t} - \frac{y}{a} = 0\end{equation}

It is shown that the auxiliary transverse latitudes are the elliptical intersection formed by planes parallel to X-O-Z and the ellipsoid, which refers to the definition method of the grid heading, and has the advantage of simple calculation and representation, but it is not the latitude measured in geodesy.

The transverse longitude at M is the included angle between O'M’ and the Z-axis, so the point with the same transverse longitude ${\lambda ^t}$ satisfies Equation (5), which is the plane passing through the Y-axis:

(5)\begin{equation}{F_2}(x,y,z) = \tan {\lambda ^t} - \frac{{x(1 - {e^2})}}{z} = 0\end{equation}

The transverse meridian and parallel, that is, the curve formed by the points of the same transverse longitude (latitude) on the ellipsoid surface, are derived by simultaneously solving F and F ₂ (or F ₁). Figure 2 shows some transverse meridians, parallels and auxiliary transverse parallels on the surface of the ellipsoid, where the blue line represents the transverse meridians, the yellow line represents the transverse parallels and the green line represents the auxiliary transverse parallels. The complexity and irregularity can be seen intuitively.

Figure 2. Schematic diagram of ellipsoid transverse meridians and parallels.

2.2 Difficulties in constructing transverse geographic coordinate system

The tangent direction of transverse meridians, parallels and auxiliary transverse parallels passing through M(${\varphi ^t}$,${\lambda ^t}$) are to be calculated first.

The relationship between the Earth Centred Earth Fixed and the transverse longitude and latitude is

(6)\begin{equation}\left\{ \begin{array}{@{}l@{}} x = {R_N}\cos {\varphi^t}\sin {\lambda^t}\\ y = {R_N}\sin {\varphi^t}\\ z = {R_N}(1 - {e^2})\cos {\varphi^t}\cos {\lambda^t} \end{array} \right.\end{equation}

Suppose the surface normal vector determined by Equation (3) is n₁, and the direction (tangent) vector of the transverse parallel is ${\boldsymbol{T}_1} = \boldsymbol{n} \times {\boldsymbol{n}_1}$; similarly, suppose the surface normal vector determined by Equation (5) is n₂, then the direction (tangent) vector of the transverse meridian is ${\boldsymbol{T}_2} = \boldsymbol{n} \times {\boldsymbol{n}_2}$. In combination with Equation (6), the transverse longitude and latitude can be used to represent the direction (tangent) vectors T₁, T₂ and T₃, which respectively represent the direction (tangent) vector of transverse parallels, transverse meridians and auxiliary transverse parallels:

(7)\begin{equation}\left\{ \begin{array}{@{}l@{}} {\boldsymbol{T}_1} = [\cos {\lambda^t}(1 - {e^2}{\cos^2}{\varphi^t}), - {e^2}\sin {\varphi^t}\cos {\varphi^t}\sin {\lambda^t}\cos {\lambda^t}, - (1 - {e^2})\sin {\lambda^t}]\\ {\boldsymbol{T}_2} = [ - \sin {\varphi^t}\sin {\lambda^t},\cos {\varphi^t}(1 - {e^2}{\cos^2}{\lambda^t}), - (1 - {e^2})\sin {\varphi^t}\cos {\lambda^t}]\\ {\boldsymbol{T}_3} = [\cos {\lambda^t},0, - \sin {\lambda^t}] \end{array} \right.\end{equation}

It is difficult to define an orthogonal transverse geographic coordinate system with reference to traditional methods because T₁ and T₂ are not perpendicular. In application, the tangent direction T₃ of the auxiliary transverse parallels is generally defined as the east direction of the transverse geographic system. The north direction of the transverse geographic system, which is not the direction determined by T₂, is determined by the right-hand rule with the up direction (normal vector on the ellipsoid).

Therefore, the transverse eastward velocity not only leads to the change of transverse longitude, but also changes the transverse latitude; the transverse north velocity not only leads to the change of transverse latitude, but also changes the transverse longitude, which is the direct reason for the complexity of the arrangement in a transverse ellipsoid INS. The deeper reason concluded is that the earth ellipsoid model is a rotating ellipsoid around the Z-axis, not around the X-axis or Y-axis, so the transverse meridians and parallels are not as regular as the traditional ones.

2.3 Key algorithms of ellipsoid transverse arrangement

Given the attitude and velocity update of an ellipsoid transverse arrangement are mostly similar to the traditional algorithms, this section mainly focuses on the position update algorithm and the calculation of the angular velocity of the transverse geographic system (t system) relative to the earth system (ECEF system) caused by the carrier motion, $\boldsymbol{\omega }_{et}^t$. Since predecessors have already studied the algorithm, most of the derivation details will be omitted here.

According to the definitions of ellipsoidal transverse longitude, latitude and the ellipsoidal transverse geographic coordinate system (t system), $\boldsymbol{C}_t^e$ is calculated as

(8)\begin{equation}\boldsymbol{C}_t^e = \left( {\begin{array}{@{}ccc@{}} {\cos {\lambda^t}}& { - \sin {\varphi^t}\sin {\lambda^t}}& {\cos {\varphi^t}\sin {\lambda^t}}\\ 0& {\cos {\varphi^t}}& {\sin {\varphi^t}}\\ { - \sin {\lambda^t}}& { - \sin {\varphi^t}\cos {\lambda^t}}& {\cos {\varphi^t}\cos {\lambda^t}} \end{array}} \right)\end{equation}

According to the matrix differential formula,

(9)\begin{equation}\dot{\boldsymbol{C}}_t^e = \boldsymbol{C}_t^e(\boldsymbol{\omega }_{et}^t \times )\end{equation}

then

(10)\begin{equation}\boldsymbol{\omega }_{et}^t = {( - {\dot{\varphi }^t},{\dot{\lambda }^t}\cos {\varphi ^t},{\dot{\lambda }^t}\sin {\varphi ^t})^\textrm{T}}\end{equation}

For the calculation of change rate of transverse longitude ${\dot{\lambda }^\textrm{t}}$ and latitude ${\dot{\varphi }^\textrm{t}}$, a decomposition formula with Earth Cartesian position vector r^e is to be used:

(11)\begin{equation}{\dot{\boldsymbol{r}}^e} = {\boldsymbol{v}^e} = \boldsymbol{C}_t^e{\boldsymbol{v}^t} = \frac{{\textrm{d}{\boldsymbol{r}^e}}}{{\textrm{d}{\varphi ^\textrm{t}}}}{\dot{\varphi }^t} + \frac{{\textrm{d}{\boldsymbol{r}^e}}}{{\textrm{d}{\lambda ^\textrm{t}}}}{\dot{\lambda }^t}\end{equation}

Equation (11) can be divided into two parts for the convenience of calculation:

(12)\begin{equation}\left\{ \begin{array}{@{}c@{}} \boldsymbol{C}_t^e{(v_{tE}^t,0,0)^{\rm T}} = \dfrac{{d{\boldsymbol{r}^e}}}{{d{\varphi^t}}}{{\dot{\varphi }}^t}_{tE}+ \dfrac{{d{\boldsymbol{r}^e}}}{{d{\lambda^t}}}{{\dot{\lambda }}^t}_{tE}\\ \boldsymbol{C}_t^e{(0,v_{tN}^t,0)^{\rm T}} = \dfrac{{d{\boldsymbol{r}^e}}}{{d{\varphi^t}}}{{\dot{\varphi }}^t}_{tN}+ \dfrac{{d{\boldsymbol{r}^e}}}{{d{\lambda^t}}}{{\dot{\lambda }}^t}_{tN} \end{array} \right.\end{equation}

where ${\dot{\varphi }^t}_{tE}$, ${\dot{\lambda }^t}_{tE}$ denotes the change rate of transverse latitude and longitude caused by eastward transverse velocity, respectively, and ${\dot{\varphi }^t}_{tN}$, ${\dot{\lambda }^\textrm{t}}_{\textrm{tN}}$ denotes the change rate of transverse latitude and longitude caused by northward transverse velocity, respectively. Substituting in Equations (6) and (8) and combining the expression of prime vertical circle R_N denoted by transverse longitude and latitude:

(13)\begin{equation}{R_N} = \frac{a}{{\sqrt {1 - {e^2}{{\cos }^2}{\varphi ^t}{{\cos }^2}{\lambda ^t}} }}\end{equation}

The position updating formula of ellipsoid transverse arrangement can be derived, which can also be used in the calculation of $\boldsymbol{\omega }_{et}^t$:

(14)\begin{equation}\left\{ \begin{array}{@{}l@{}} {{\dot{\varphi }}^t} = \dot{\varphi }_{tE}^t + \dot{\varphi }_{tN}^t = \dfrac{{V_{tN}^t(1 - {e^2}{{\sin }^2}{\lambda^t} - {e^2}{{\cos }^2}{\varphi^t}{{\cos }^2}{\lambda^t}) + {e^2}V_{tE}^t\sin {\varphi^t}\sin {\lambda^t}\cos {\lambda^t}}}{{{R_N}(1 - {e^2})}}\\ {{\dot{\lambda }}^t} = \dot{\lambda }_{tE}^t + \dot{\lambda }_{tN}^t = \dfrac{{V_{tE}^t\sec {\varphi^t}(1 - {e^2}{{\cos }^2}{\lambda^t}) + {e^2}V_{tN}^t\tan {\varphi^t}\sin {\lambda^t}\cos {\lambda^t}}}{{{R_N}(1 - {e^2})}} \end{array} \right.\end{equation}

3. Construction of virtual sphere model and its application in transverse navigation

3.1 Virtual sphere model in traditional coordinate system

According to the literature (Qin et al., 2018), the idea of virtual sphere model design can be extracted: (a) based on the fact that the radius of the earth is closely related to its velocity in inertial navigation calculation; (b) a ‘virtual sphere’ model can be regarded as an artificial ball with radius of curvature in Prime Vertical. It can be concluded that the transition of the velocity and the unified construction of radius are of key importance.

This can be further explained in the traditional coordinate system. In the movement of the surface carrier, the curvature radius related to the eastward velocity in the INS position update is the Prime Vertical radius R_N, and the curvature radius R_M of the Prime Meridian circle is related to the northward velocity. If the northward velocity is multiplied by the ratio of R_N to R_M, which is defined as the virtual northward velocity, and the virtual northward velocity in the position update equation is to be divided equally with the original eastward velocity by R_N, then the constructed virtual resultant velocity is equal to the combination of the virtual northward velocity and the original eastward velocity, as if the carrier moves on the sphere with radius R_N at the virtual resultant velocity. See Figure 3.

Figure 3. Transformation from ellipsoid to virtual sphere.

Figure 3(a) shows the motion in the real situation, with the geographical eastward speed v_E and northward speed v_N.

Figure 3(b) shows the virtual sphere model corresponding to Figure 3(a), with the position (i.e. the longitude and latitude) of the carrier unchanged. The radius of the virtual sphere is R_N, and the carrier speed is replaced by the virtual speeds v_E′ and v_N′:

(15)\begin{equation}\left\{ \begin{array}{@{}l@{}} {v_E}^\prime = {v_E}\\ {v_N}^\prime = \dfrac{{{R_N}}}{{{R_M}}}{v_N} \end{array} \right.\end{equation}

The position update formula is

(16)\begin{equation}(\dot{\varphi },\dot{\lambda }) = \left( {\frac{{{v_N}}}{{{R_M}}},\frac{{{v_E}}}{{{R_N}\cos \varphi }}} \right) = \frac{1}{{{R_N}}}({v_N}^\prime , {v_E}^\prime \sec \varphi )\end{equation}

In the traditional arrangement, there is obviously no difference between the virtual sphere method and the general algorithm, which can only be regarded as a new way to understand the problem. However, in the transverse arrangement of the ellipsoid, the advantages of applying the virtual sphere model are revealed.

3.2 Key algorithms of INS transverse arrangement based on virtual sphere model

The core idea of the arrangement is to transfer the transverse velocity to the geographical system for changing the northward velocity v_N, and then the processed resultant velocity is transferred back to the transverse system (t system) but substitutes into the virtual sphere with R_N as the radius for transverse position update. The process avoids complex calculation of many related curvature radii.

Therefore, the relationship between the virtual transverse velocity ${\boldsymbol{v}^{t\prime}} = {(v{_{tE}^{t\prime}},v{_{tN}^{t\prime}},0)^T}$ and the real transverse velocity ${\boldsymbol{v}^t} = {(v_{tE}^t,v_{tN}^t,0)^T}$ is

(17)\begin{equation}{\boldsymbol{v}^{t\prime} } = \boldsymbol{C}_g^t\left( {\begin{array}{*{20}{c}} 1& {}& {}\\ {}& {\dfrac{{{R_N}}}{{{R_M}}}}& {}\\ {}& {}& 1 \end{array}} \right)\boldsymbol{C}_t^g{\boldsymbol{v}^t}\end{equation}

where $\boldsymbol{C}_g^t$ is the transformation matrix from transverse geographic coordinate system to geographic coordinate system (g system):

(18)\begin{equation}\boldsymbol{C}_g^t = \left( {\begin{array}{*{20}{c}} {\cos \sigma }& { - \sin \sigma }& 0\\ {\sin \sigma }& {\cos \sigma }& 0\\ 0& 0& 1 \end{array}} \right)\end{equation}

and

(19)\begin{equation}\left\{ \begin{array}{@{}c@{}} \cos \sigma = \dfrac{{ - \sin {\varphi^t}\cos {\lambda^t}}}{{\sqrt {1 - {{\cos }^2}{\varphi^t}{{\cos }^2}{\lambda^t}} }}\\ \sin \sigma = \dfrac{{\sin {\lambda^t}}}{{\sqrt {1 - {{\cos }^2}{\varphi^t}{{\cos }^2}{\lambda^t}} }} \end{array} \right.\end{equation}

Using the virtual transverse velocity and the radius R_N, the transverse position updating equation and the calculation formula of $\boldsymbol{\omega }_{et}^t$ can be written as

(20)\begin{equation}({\dot{\varphi }^t},{\dot{\lambda }^t}) = \frac{1}{{{R_N}}}(v{_{tN}^{t\prime}},v{_{tE}^{t\prime}}\sec {\varphi ^t})\end{equation}

(21)\begin{equation}\boldsymbol{\omega }_{et}^t = \frac{1}{{{R_N}}}{( - v{_{tN}^{t\prime}},v{_{tE}^{t\prime}},v{_{tE}^{t\prime}}\tan {\varphi ^t})^T}\end{equation}

3.3 Equivalence proof of virtual sphere model algorithm and ellipsoid algorithm on the ellipsoid surface

It is enough to prove that the calculation formulae for position updating between two arrangements are the same.

For the algorithm based on the virtual sphere model, substituting Equations (13), (18) and (19) into Equation (17) and combining Equation (22) yields

(22)\begin{equation}{R_M} = \frac{{{R_N}(1 - {e^2})}}{{1 - {e^2}{{\cos }^2}{\varphi ^t}{{\cos }^2}{\lambda ^t}}}\end{equation}

The expression of virtual transverse velocity is calculated by

(23)\begin{equation}\left\{ \begin{array}{@{}l@{}} v{_{tE}^{t\prime}} = \dfrac{{V_{tE}^t[ - 2 + {e^2} + {e^2}\cos (2{\lambda^t})] - {e^2}V_{tN}^t\sin (2{\lambda^t})\sin {\varphi^t}}}{{2( - 1 + {e^2})}}\\ v{_{tN}^{t\prime} } = \dfrac{{( - 4 + 3{e^2})V_{tN}^t - {e^2}[V_{tN}^t\cos (2{\lambda^t}) - 2V_{tN}^t{{\cos }^2}{\lambda^t}\cos (2{\varphi^t}) + 2V_{tE}^t\sin (2{\lambda^t})\sin {\varphi^t}]}}{{4( - 1 + {e^2})}} \end{array} \right.\end{equation}

Substituting Equation (23) into Equation (20) and making a simplification, the change rate of transverse longitude and latitude is shown to be completely equal to Equation (14).

3.4 Equivalence proof of virtual sphere model algorithm and ellipsoid algorithm with height change

The cases discussed above are all on the ellipsoidal surface. To expand the application scope of the algorithm, it is necessary to deduce a formula for when the height changes.

With altitude h, Equation (6) should be changed to

(24)\begin{equation}\left\{ \begin{array}{@{}l@{}} x = ({R_N} + h)\cos {\varphi^t}\sin {\lambda^t}\\ y = ({R_N} + h)\sin {\varphi^t}\\ z = [{R_N}(1 - {e^2}) + h]\cos {\varphi^t}\cos {\lambda^t} \end{array} \right.\end{equation}

and Equation (17) should be changed to

(25)\begin{equation}{\boldsymbol{v}^{t\prime}} = \boldsymbol{C}_g^t\left( {\begin{array}{*{20}{c}} 1& {}& {}\\ {}& {\dfrac{{{R_N} + h}}{{{R_M} + h}}}& {}\\ {}& {}& 1 \end{array}} \right)\boldsymbol{C}_t^g{\boldsymbol{v}^t}\end{equation}

In addition, calculation formulae such as Equation (16) should change R_N into R_N + h, while other calculation formulae will not be changed. Through software derivation, the transverse longitude and latitude change rate of the virtual sphere model algorithm and the ellipsoid model algorithm are derived to be equal.

So far, it is proved that the virtual sphere model algorithm is equivalent to the ellipsoid algorithm in transverse arrangement INS.

4. Conclusion and further discussion

Through analysing the principle of the INS ellipsoid transverse arrangement, this paper points out the root causes of the complexity of the ellipsoid transverse arrangement. Aiming at the newly proposed virtual sphere model algorithm, this paper briefly expounds the core idea and reveals the principle that it can avoid complex operation. The difference between the two algorithms is only the calculation involved in a position update, so the equivalence of the two algorithms is proved, which provides a sufficient theoretical foundation for the application of the virtual sphere model algorithm.

The following are points for further research.

1. 1. The proof of the equivalence may broaden the method of solving the problems on INS update with the unconventional coordinate system, which is the path from a special geographic coordinate system (e.g. oblique coordinate systems and other pseudo coordinate systems) to a conventional geographic coordinate system, to a virtual sphere system, and finally back to the special coordinate system itself, and this is a certain contribution to the development of the inertial navigation algorithm if proved effective.

2. 2. Some implicit relations between ellipsoid and sphere in mathematical analysis may be revealed by the proof, and potential practical value in other fields may be found.