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.1 Virtual sphere model in traditional coordinate system

According to the literature (Qin et al., Reference Qin, Chang, Tong, Wang, Huang and Zha2018), 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.