A new polar alignment algorithm based on the Huber estimation filter with the aid of BeiDou Navigation Satellite System

Abstract

For aircrafts equipped with BeiDou Navigation Satellite System/Strapdown Inertial Navigation System integrated navigation system, BeiDou Navigation Satellite System information can be used to achieve autonomous alignment. However, due to the complex polar environment and multipath effect, BeiDou Navigation Satellite System measurement noise often exhibits a non-Gaussian distribution that will severely degrade the estimation accuracy of standard Kalman filter. To address this problem, a new polar alignment algorithm based on the Huber estimation filter is proposed in this article. Considering the special geographical conditions in the polar regions, the dynamic model and the measurement model of BeiDou Navigation Satellite System/Strapdown Inertial Navigation System integrated alignment system in the grid frame are derived in this article. The BeiDou Navigation Satellite System measurement noise characteristics in the polar regions are analyzed and heavy-tailed characteristics are simulated, respectively. Since the estimation accuracy of standard Kalman filter can be severely degraded under non-Gaussian noise, a Kalman filter based on the Huber estimation is designed combining grid navigation system and generalized maximum likelihood estimation. The simulation and experiment results demonstrate that the proposed algorithm has better robustness under non-Gaussian noise, and it is effective in the polar regions. By employing the proposed algorithm, the rapidity and accuracy of the alignment process can be improved.

Introduction

The polar regions are becoming more and more important for all mankind, and the polar flight ability is of great significance.^1–3 The lack of high-precision navigation is one of the most prominent problems in polar navigation. As the most important autonomous navigation equipment, Strapdown Inertial Navigation System (SINS) can provide all kinds of high-frequency navigation information while maintaining good concealment and real-time performance. Compared with other navigation systems, SINS has better navigation performance in high latitude areas and is the optimal choice for polar navigation.^4,5 Accelerometers and gyroscopes are employed to measure the motion parameters of the aircraft. Under precise initial conditions, the navigation results can be obtained by recursive calculation.^7–9

The performance of initial alignment will affect the navigation accuracy. Two major issues need to be addressed for initial alignment in the polar regions. One is that the traditional north-oriented frame will lose effectiveness because of the rapid convergence of meridians. ¹⁰ The other is that the stationary base alignment cannot be achieved due to the decrease of the horizontal component of rotational angular velocity of the earth.

The traditional local-level north-oriented frame is widely used in the non-polar regions. Its x-axis and y-axis are both in the local level while z-axis coincides with geography up. The y-axis always points to true north as the heading reference. However, due to the rapid convergence of meridians in the polar regions, the traditional north-oriented frame cannot work properly. Although the local-level wander-azimuth frame and the local-level free-azimuth frame are proposed, there are singular values when extracting the wander-azimuth angle and longitude from the position direction cosine matrix. ¹¹ In order to solve this problem, transversal SINS mechanization and grid SINS mechanization are proposed. In the transversal SINS mechanization, there is principle error when using the earth standard sphere model. However, the calculation is more complicated when using the earth ellipsoid model. ¹⁰ Therefore, we adopt the grid SINS mechanization for polar alignment.

In the polar regions, the horizontal component of rotational angular velocity of the earth decreases as latitude increases. Therefore, stationary base alignment cannot be achieved. It is necessary to use external information to assist the moving base alignment. A new polar alignment algorithm with the aid of star sensor is proposed and the accurate information provided by star sensor is employed as the measurement to improve the performance of initial alignment. ¹² The OCTANS composed of three fiber-optic gyroscopes and three quartz accelerometers can provide the attitude information of UUV (Unmanned Underwater Vehicles) quickly and accurately, and it is used to improve the alignment accuracy of polar UUV. ¹³ An SINS/radar alignment algorithm is proposed to improve the attitude estimation accuracy. ¹ The output of DVL (Doppler Velocity Log) is employed as the measurement in the polar transversal initial alignment algorithm. ¹⁴ The transfer alignment using high-precision master inertial navigation system (INS) output to estimate the states is researched in the polar regions.^15,16 The attitude initialization with the aid of GNSS (Global Navigation Satellite System) in the earth frame could be safely performed in the polar regions. ¹⁷ However, it only studies the coarse alignment method, and it does not consider the influence of the complex polar environment on GNSS performance. Research works have shown that the multipath effect is obvious in high latitude areas and will affect the measurement noise characteristics of GNSS.^14–16,18

The global networking of BeiDou Navigation Satellite System (BDS) was completed on 23 June 2020. It is of great significance to study the global applications of BDS. BDS consists of GEO (Geostationary Earth Orbit) satellites, IGSO (Inclined Geosynchronous Orbit) satellites, and MEO (Medium Earth Orbit) satellites. Compared with MEO satellites, IGSO satellites have longer visibility in high latitude areas. BDS can provide navigation services independently in the polar regions. ¹⁹ In recent years, BDS/SINS integrated navigation system is increasingly used in many fields due to their integrated performance outweighs the drawbacks of the individual systems.^6,20,21 For aircrafts equipped with BDS/SINS integrated navigation system, BDS information can be used to initialize SINS and achieve autonomous alignment. However, BDS also has multipath effect in the polar regions. Furthermore, the performance of standard Kalman filter for BDS/SINS integrated alignment will be degraded. Particle filter and Gaussian-sum filter are introduced to overcome the limitation of standard Kalman filter on noise characteristics. However, high computational complexity and memory burden are also mentioned. ²²

A new polar alignment algorithm based on the Huber estimation filter is proposed in this article. The main contribution of this article is to analyze the BDS measurement noise characteristics in the polar regions and to propose a Kalman filter based on the Huber estimation combining grid navigation system and generalized maximum likelihood estimation. The remaining sections are arranged as follows. SINS mechanization and navigation error equation in the grid frame are represented in section “SINS mechanization and navigation error equation in the grid frame”; based on section “SINS mechanization and navigation error equation in the grid frame,” the dynamic model and the measurement model are established in section “The design of Kalman filter model for polar alignment with the aid of BDS”; section “The BDS measurement noise characteristics in the polar regions and the design of Kalman filter based on the Huber estimation” is devoted to analyze the BDS measurement noise characteristics in the polar regions. The Kalman filter based on the Huber estimation is also designed in section “The BDS measurement noise characteristics in the polar regions and the design of Kalman filter based on the Huber estimation” to implement polar alignment. The simulation and the experiment are performed in section “Simulation and experiment results.” In addition, the results are discussed in section “Simulation and experiment results.” Finally, the conclusion is given in section “Conclusion.”

SINS mechanization and navigation error equation in the grid frame

The Greenwich meridian is selected as the heading reference to define the grid frame, and the problem caused by meridians convergence in the polar regions can be avoided.

The grid frame

The definition of G frame (the grid frame) is shown in Figure 1. $L$ and $\lambda$ are the latitude and the longitude of the aircraft, respectively. Point P represents the location of the aircraft. The plane passing through point P and parallel to the Greenwich meridian is defined as grid plane. The horizontal plane of point P is defined as tangent plane. The intersecting line of grid plane and tangent plane is defined as grid north. The angle between grid north and true north is defined as grid angle $\sigma$ . Grid up coincides with geography up, and grid east makes up the right-handed Cartesian coordinate with grid north in the tangent plane. The transform relationship among e frame (the earth centered earth fixed frame—x-axis intersects the Greenwich meridian in the equatorial plane, y-axis is perpendicular to x-axis in the equatorial plane, and z-axis coincides with the earth axis), g frame (the geography frame, east-north-up), and G frame is shown in Figure 1. ${\mathit{C}}_{\mathrm{e}}^{\mathrm{g}}$ represents the direction cosine matrix from e frame to g frame, and ${\mathit{C}}_{\mathrm{g}}^{\mathrm{G}}$ represents the direction cosine matrix from g frame to G frame

$${\mathit{C}}_{\mathrm{e}}^{\mathrm{g}}=\left[\begin{array}{ccc}-\sin \lambda \hfill & \cos \lambda \hfill & 0\hfill \\ \multicolumn{1}{c}{-\sin L\cos \lambda }& -\sin L\sin \lambda \hfill & \cos L\hfill \\ \multicolumn{1}{c}{\cos L\cos \lambda }& \cos L\sin \lambda \hfill & \sin L\hfill \end{array}\right]$$

(1)

$${\mathit{C}}_{\mathrm{g}}^{\mathrm{G}}=\left[\begin{array}{ccc}\cos \sigma \hfill & -\sin \sigma \hfill & 0\hfill \\ \multicolumn{1}{c}{\sin \sigma }& \cos \sigma \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right]$$

(2)

Thus, the direction cosine matrix from e frame to G frame will be

$$\begin{array}{c}{\mathit{C}}_{\mathrm{e}}^{\mathrm{G}}={\mathit{C}}_{\mathrm{g}}^{\mathrm{G}}{\mathit{C}}_{\mathrm{e}}^{\mathrm{g}}=\hfill \\ \left[\begin{array}{ccc}-\cos \sigma \sin \lambda +\sin \sigma \cos \lambda \sin L\hfill & \cos \sigma \cos \lambda +\sin \sigma \sin \lambda \sin L\hfill & -\sin \sigma \sin L\hfill \\ \multicolumn{1}{c}{-\sin \sigma \sin \lambda -\cos \sigma \cos \lambda \sin L}& \sin \sigma \cos \lambda -\cos \sigma \sin \lambda \sin L\hfill & \cos \sigma \cos L\hfill \\ \multicolumn{1}{c}{\cos L\cos \lambda }& \cos L\sin \lambda \hfill & \sin L\hfill \end{array}\right]\hfill \end{array}$$

(3)

Figure 1.

The description of the grid frame.

The axial unit vectors for G frame are denoted as $({\mathit{e}}_{{\mathrm{G}}_{\mathrm{E}}},{\mathit{e}}_{{\mathrm{G}}_{\mathrm{N}}},{\mathit{e}}_{{\mathrm{G}}_{\mathrm{U}}})$ , $({\mathit{e}}_{\mathrm{E}},{\mathit{e}}_{\mathrm{N}},{\mathit{e}}_{\mathrm{U}})$ for g frame, and $({\mathit{e}}_{\mathrm{X}},{\mathit{e}}_{\mathrm{Y}},{\mathit{e}}_{\mathrm{Z}})$ for e frame. It can be seen from the definition of the angle $\sigma$ and Figure 1 that ${\mathit{e}}_{{\mathrm{G}}_{\mathrm{N}}}$ is orthogonal to ${\mathit{e}}_{\mathrm{Y}}$

$${\mathit{e}}_{{\mathrm{G}}_{\mathrm{N}}}·{\mathit{e}}_{\mathrm{Y}}=0$$

(4)

Projecting ${\mathit{e}}_{{\mathrm{G}}_{\mathrm{N}}}$ and ${\mathit{e}}_{\mathrm{Y}}$ to the g frame, hence

$${\mathit{e}}_{{\mathrm{G}}_{\mathrm{N}}}^{\mathrm{g}}={\mathit{C}}_{\mathrm{G}}^{\mathrm{g}}{\mathit{e}}_{{\mathrm{G}}_{\mathrm{N}}}^{\mathrm{G}}={\left[\begin{array}{ccc}\sin \sigma \hfill & \cos \sigma \hfill & 0\hfill \end{array}\right]}^{\mathrm{T}}$$

(5)

$${\mathit{e}}_{\mathrm{Y}}^{\mathrm{g}}={\mathit{C}}_{\mathrm{e}}^{\mathrm{g}}{\mathit{e}}_{\mathrm{Y}}^{\mathrm{e}}={\left[\begin{array}{ccc}\cos \lambda \hfill & -\sin L\sin \lambda \hfill & \cos L\sin \lambda \hfill \end{array}\right]}^{\mathrm{T}}$$

(6)

Substituting equations (5) and (6) into equation (4), we can obtain

$$\sin \sigma \cos \lambda -\cos \sigma \sin L\sin \lambda =0$$

(7)

Then

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

(8)

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

(9)

Substituting equations (8) and (9) into equations (2) and (3), the direction cosine matrices ${\mathit{C}}_{\mathrm{g}}^{\mathrm{G}}$ and ${\mathit{C}}_{\mathrm{e}}^{\mathrm{G}}$ can be calculated.

Grid SINS mechanization

The mechanization of grid SINS is

$$\begin{array}{c}{\stackrel{·}{\mathit{C}}}_{\mathrm{b}}^{\mathrm{G}}={\mathit{C}}_{\mathrm{b}}^{\mathrm{G}}({\mathit{\omega }}_{\mathrm{Gb}}^{\mathrm{b}}\times )\hfill \\ {\stackrel{·}{\mathit{v}}}^{\mathrm{G}}={\mathit{C}}_{\mathrm{b}}^{\mathrm{G}}{\mathit{f}}^{\mathrm{b}}-(2{\mathit{\omega }}_{\mathrm{ie}}^{\mathrm{G}}+{\mathit{\omega }}_{\mathrm{eG}}^{\mathrm{G}})\times {\mathit{v}}^{\mathrm{G}}+{\mathit{g}}^{\mathrm{G}}\hfill \\ {\stackrel{·}{\mathit{P}}}^{\mathrm{e}}={\mathit{C}}_{\mathrm{G}}^{\mathrm{e}}{\mathit{v}}^{\mathrm{G}}\hfill \end{array}$$

(10)

where ${\mathit{C}}_{\mathrm{b}}^{\mathrm{G}}$ is the direction cosine matrix from b frame (the body frame, Right-Front-Up) to G frame and ${\stackrel{·}{\mathit{C}}}_{\mathrm{b}}^{\mathrm{G}}$ represents the derivation of ${\mathit{C}}_{\mathrm{b}}^{\mathrm{G}}$ ${\mathit{\omega }}_{\mathrm{Gb}}$ is the angular velocity of b frame relative to G frame, and ${\mathit{\omega }}_{\mathrm{Gb}}^{\mathrm{b}}$ is the projection of ${\mathit{\omega }}_{\mathrm{Gb}}$ on b frame; ${\mathit{\omega }}_{\mathrm{ie}}$ is rotational angular velocity of the earth, and ${\mathit{\omega }}_{\mathrm{ie}}^{\mathrm{G}}$ is the projection of ${\mathit{\omega }}_{\mathrm{ie}}$ on G frame; ${\mathit{\omega }}_{\mathrm{eG}}$ is the angular velocity of G frame relative to e frame, and ${\mathit{\omega }}_{\mathrm{eG}}^{\mathrm{G}}$ is the projection of ${\mathit{\omega }}_{\mathrm{eG}}$ on G frame; ${\mathit{v}}^{\mathrm{G}}$ is the velocity of SINS in G frame and ${\stackrel{·}{\mathit{v}}}^{\mathrm{G}}$ represents the derivation of ${\mathit{v}}^{\mathrm{G}}$ ; ${\mathit{f}}^{\mathrm{b}}$ is the specific force in b frame measured by SINS; ${\mathit{g}}^{\mathrm{G}}$ is the projection of gravitational acceleration on G frame; ${\mathit{P}}^{\mathrm{e}}$ is the position coordinates of the aircraft in e frame and ${\stackrel{·}{\mathit{P}}}^{\mathrm{e}}$ represents the derivation of ${\mathit{P}}^{\mathrm{e}}$ ; ${\mathit{C}}_{\mathrm{G}}^{\mathrm{e}}$ is the direction cosine matrix from G frame to e frame; $({\mathit{\omega }}_{\mathrm{Gb}}^{\mathrm{b}}\times )$ represents the skew-symmetric matrix of ${\mathit{\omega }}_{\mathrm{Gb}}^{\mathrm{b}}$ , and it can be defined as:

$$({\mathit{\omega }}_{\mathrm{Gb}}^{\mathrm{b}}\times )=\left[\begin{array}{ccc}0\hfill & -{\mathit{\omega }}_{\mathrm{Gb}z}^{\mathrm{b}}\hfill & {\mathit{\omega }}_{\mathrm{Gb}y}^{\mathrm{b}}\hfill \\ \multicolumn{1}{c}{{\mathit{\omega }}_{\mathrm{Gb}z}^{\mathrm{b}}}& 0\hfill & -{\mathit{\omega }}_{\mathrm{Gb}x}^{\mathrm{b}}\hfill \\ \multicolumn{1}{c}{-{\mathit{\omega }}_{\mathrm{Gb}y}^{\mathrm{b}}}& {\mathit{\omega }}_{\mathrm{Gb}x}^{\mathrm{b}}\hfill & 0\hfill \end{array}\right]$$

(11)

${\mathit{\omega }}_{\mathrm{eG}}^{\mathrm{G}}$ can be calculated as follows

$$\begin{array}{c}{\mathit{\omega }}_{\mathrm{eG}}^{\mathrm{G}}=\left[\begin{array}{c}{\omega }_{\mathrm{eG}x}^{\mathrm{G}}\hfill \\ \multicolumn{1}{c}{{\omega }_{\mathrm{eG}y}^{\mathrm{G}}}\\ \multicolumn{1}{c}{{\omega }_{\mathrm{eG}z}^{\mathrm{G}}}\end{array}\right]=\left[\begin{array}{cc}(\frac{1}{R_{\mathrm{Mh}}}-\frac{1}{R_{\mathrm{Nh}}})\sin \sigma \cos \sigma \hfill & -\frac{{\cos }^2\sigma }{R_{\mathrm{Mh}}}-\frac{{\sin }^2\sigma }{R_{\mathrm{Nh}}}\hfill \\ \multicolumn{1}{c}{\frac{{\sin }^2\sigma }{R_{\mathrm{Mh}}}+\frac{{\cos }^2\sigma }{R_{\mathrm{Nh}}}}& -(\frac{1}{R_{\mathrm{Mh}}}-\frac{1}{R_{\mathrm{Nh}}})\sin \sigma \cos \sigma \hfill \\ \multicolumn{1}{c}{(\frac{1}{R_{\mathrm{Mh}}}-\frac{1}{R_{\mathrm{Nh}}})\frac{\sin \lambda \cos L}{\sqrt{1-{\cos }^{2}L{\sin }^2\lambda }}\sin \sigma \cos \sigma }& -(\frac{{\cos }^2\sigma }{R_{\mathrm{Mh}}}+\frac{{\sin }^2\sigma }{R_{\mathrm{Nh}}})\frac{\sin \lambda \cos L}{\sqrt{1-{\cos }^{2}L{\sin }^2\lambda }}\hfill \end{array}\right]\left[\begin{array}{c}{v}_{\mathrm{GE}}\hfill \\ \multicolumn{1}{c}{v_{\mathrm{GN}}}\end{array}\right]\hfill \end{array}$$

(12)

where $R_{\mathrm{Mh}}$ and $R_{\mathrm{Nh}}$ are the radius of curvature in meridian and prime vertical, respectively.

${\mathit{\omega }}_{\mathrm{Gb}}^{\mathrm{b}}$ can be obtained as follows

$${\mathit{\omega }}_{\mathrm{Gb}}^{\mathrm{b}}={\mathit{\omega }}_{\mathrm{ib}}^{\mathrm{b}}-{\mathit{C}}_{\mathrm{G}}^{\mathrm{b}}({\mathit{\omega }}_{\mathrm{ie}}^{\mathrm{G}}+{\mathit{\omega }}_{\mathrm{eG}}^{\mathrm{G}})$$

(13)

where ${\mathit{C}}_{\mathrm{G}}^{\mathrm{b}}$ is the direction cosine matrix from G frame to b frame, ${\mathit{\omega }}_{\mathrm{ib}}$ is the angular velocity of b frame relative to i frame (earth centered inertial), and ${\mathit{\omega }}_{\mathrm{ib}}^{\mathrm{b}}$ is the projection of ${\mathit{\omega }}_{\mathrm{ib}}$ on b frame which could be measured by gyroscopes.

Navigation error equation in the grid frame

The attitude error equation and velocity error equation in G frame can be expressed as

$$\{\begin{array}{c}{\stackrel{·}{\mathit{\varphi }}}^{\mathrm{G}}=-({\mathit{\omega }}_{\mathrm{iG}}^{\mathrm{G}}\times ){\mathit{\varphi }}^{\mathrm{G}}+{\mathit{C}}_v\times \delta {\mathit{v}}^{\mathrm{G}}-{\mathit{C}}_{\mathrm{b}}^{\mathrm{G}}\times \delta {\mathit{\omega }}_{\mathrm{ib}}^{\mathrm{b}}\hfill \\ \delta {\stackrel{·}{\mathit{v}}}^{\mathrm{G}}={\mathit{C}}_{\mathrm{b}}^{\mathrm{G}}{\mathit{f}}^{\mathrm{b}}\times {\mathit{\varphi }}^{\mathrm{G}}+({\mathit{v}}^{\mathrm{G}}\times {\mathit{C}}_v-(2{\mathit{\omega }}_{\mathrm{ie}}^{\mathrm{G}}+{\mathit{\omega }}_{\mathrm{eG}}^{\mathrm{G}})\times )\delta {\mathit{v}}^{\mathrm{G}}+{\mathit{C}}_{\mathrm{b}}^{\mathrm{G}}\times \delta {\mathit{f}}^{\mathrm{b}}\hfill \end{array}$$

(14)

where ${\mathit{\varphi }}^{\mathrm{G}}$ is the misalignment angle in G frame and ${\stackrel{·}{\mathit{\varphi }}}^{\mathrm{G}}$ represents the derivation of ${\mathit{\varphi }}^{\mathrm{G}}$ ; $\delta {\mathit{v}}^{\mathrm{G}}$ is the velocity error in G frame and $\delta {\stackrel{·}{\mathit{v}}}^{\mathrm{G}}$ represents the derivation of $\delta {\mathit{v}}^{\mathrm{G}}$ ; ${\mathit{\omega }}_{\mathrm{iG}}$ is the angular velocity of G frame relative to i frame, and ${\mathit{\omega }}_{\mathrm{iG}}^{\mathrm{G}}$ is the projection of ${\mathit{\omega }}_{\mathrm{iG}}$ on G frame; ${\mathit{C}}_v$ can be expressed as

$${\mathit{C}}_v=\frac{1}{\sqrt{R_{\mathrm{Mh}}{R}_{\mathrm{Nh}}}}\left[\begin{array}{ccc}0\hfill & -1\hfill & 0\hfill \\ \multicolumn{1}{c}{1}& 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& -\cot L\sin \sigma \hfill & 0\hfill \end{array}\right]$$

and $\delta {\mathit{\omega }}_{\mathrm{ib}}^{\mathrm{b}}$ and $\delta {\mathit{f}}^{\mathrm{b}}$ represent the measurement errors of the gyroscopes and accelerometers

$$\{\begin{array}{c}\delta {\mathit{\omega }}_{\mathrm{ib}}^{\mathrm{b}}={\mathit{\epsilon }}_{\mathrm{b}}^{\mathrm{b}}+{\mathit{\epsilon }}_{\mathrm{r}}^{\mathrm{b}}+{\mathit{\epsilon }}_{\mathrm{w}}^{\mathrm{b}}\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{\mathit{\epsilon }}}_{\mathrm{b}}^{\mathrm{b}}=0}\\ \multicolumn{1}{c}{{\stackrel{·}{\mathit{\epsilon }}}_{\mathrm{r}}^{\mathrm{b}}=-\frac{1}{T_g}{\mathit{\epsilon }}_{\mathrm{r}}^{\mathrm{b}}+{\mathit{w}}_{\mathrm{r}}^{\mathrm{b}}}\end{array},\{\begin{array}{c}\delta {\mathit{f}}^{\mathrm{b}}={\nabla }_{\mathrm{a}}^{\mathrm{b}}+{\nabla }_{\mathrm{w}}^{\mathrm{b}}\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{\nabla }}_{\mathrm{a}}^{\mathrm{b}}=-\frac{1}{T_{\mathrm{a}}}{\nabla }_{\mathrm{a}}^{\mathrm{b}}+{\mathit{w}}_{\mathrm{a}}^{\mathrm{b}}}\end{array}$$

(15)

where ${\mathit{\epsilon }}_{\mathrm{b}}^{\mathrm{b}},{\mathit{\epsilon }}_{\mathrm{r}}^{\mathrm{b}},{\mathit{\epsilon }}_{\mathrm{w}}^{\mathrm{b}}$ are the random constant, the first-order Markov drift, and the white noise of the gyroscopes, respectively; ${\nabla }_{\mathrm{a}}^{\mathrm{b}},{\nabla }_{\mathrm{w}}^{\mathrm{b}}$ are the first-order Markov drift and the white noise of the accelerometers, respectively. $T_g,{\mathit{w}}_{\mathrm{r}}^{\mathrm{b}}$ are the correlation time and the driving white noise of gyroscope first-order Markov drift, respectively. $T_a,{\mathit{w}}_{\mathrm{a}}^{\mathrm{b}}$ are the correlation time and the driving white noise of accelerometer first-order Markov drift, respectively.

The design of Kalman filter model for polar alignment with the aid of BDS

To implement SINS alignment with the aid of BDS, BDS/SINS integrated alignment filter model should be established.

BDS/SINS integrated alignment state equation

The misalignment angle ${\mathit{\varphi }}^{\mathrm{G}}$ , the velocity error $\delta {\mathit{v}}^{\mathrm{G}}$ , the random constant and first-order Markov drift of the gyroscopes ${\mathit{\epsilon }}_{\mathrm{b}}^{\mathrm{b}},{\mathit{\epsilon }}_{\mathrm{r}}^{\mathrm{b}}$ , and the first-order Markov drift of the accelerometers ${\nabla }_{\mathrm{a}}^{\mathrm{b}}$ are chosen as the states of alignment filter

$$\mathit{X}\left(t\right)={\left[\begin{array}{ccccc}{\mathit{\varphi }}^{\mathrm{G}}\hfill & \delta {\mathit{v}}^{\mathrm{G}}\hfill & {\mathit{\epsilon }}_{\mathrm{b}}^{\mathrm{b}}\hfill & {\mathit{\epsilon }}_{\mathrm{r}}^{\mathrm{b}}\hfill & {\nabla }_{\mathrm{a}}^{\mathrm{b}}\hfill \end{array}\right]}^{\mathrm{T}}$$

(16)

According to equations (14) and (15), the system state equation can be rewritten as

$$\stackrel{·}{\mathit{X}}\left(t\right)=\mathit{A}\left(t\right)\mathit{X}\left(t\right)+\mathit{G}\left(t\right)\mathit{W}\left(t\right)$$

(17)

where $\mathit{X}(t)$ is the state vector and $\stackrel{·}{\mathit{X}}(t)$ represents the derivation of $\mathit{X}(t)$ , $\mathit{A}(t)$ is the state transition matrix, $\mathit{G}(t)$ is the system noise coefficient matrix, and $\mathit{W}(t)$ is the system noise vector. $\mathit{A}(t)$ can be obtained as follows

$$\mathit{A}=\left[\begin{array}{ccccc}-({\mathit{\omega }}_{\mathrm{iG}}^{\mathrm{G}}\times )\hfill & {\mathit{C}}_v\hfill & -{\mathit{C}}_{\mathrm{b}}^{\mathrm{G}}\hfill & -{\mathit{C}}_{\mathrm{b}}^{\mathrm{G}}\hfill & {\mathit{0}}_{3\times 3}\hfill \\ \multicolumn{1}{c}{{\mathit{C}}_{\mathrm{b}}^{\mathrm{G}}({\mathit{f}}^{\mathrm{b}}\times )}& {\mathit{v}}^{\mathrm{G}}\times {\mathit{C}}_v-(2{\mathit{\omega }}_{\mathrm{ie}}^{\mathrm{G}}+{\mathit{\omega }}_{\mathrm{eG}}^{\mathrm{G}})\times \hfill & {\mathit{0}}_{3\times 3}\hfill & {\mathit{0}}_{3\times 3}\hfill & {\mathit{C}}_{\mathrm{b}}^{\mathrm{G}}\hfill \\ \multicolumn{1}{c}{{\mathit{0}}_{9\times 3}}& {\mathit{0}}_{9\times 3}\hfill & \hfill & \mathrm{diag}\left(\begin{array}{ccccccccc}0\hfill & 0\hfill & 0\hfill & -\frac{1}{T_g}\hfill & -\frac{1}{T_g}\hfill & -\frac{1}{T_g}\hfill & -\frac{1}{T_a}\hfill & -\frac{1}{T_a}\hfill & -\frac{1}{T_a}\hfill \end{array}\right)\hfill & \hfill \end{array}\right]$$

(18)

BDS/SINS integrated alignment measurement equation

The velocity difference between SINS and BDS is chosen as the measurement. The velocity of SINS and BDS in G frame can be expressed as ${\mathit{v}}^{\mathrm{G}},{\mathit{v}}_B^{\mathrm{G}}$ , and the velocity measurement equation can be defined as

$${\mathit{Y}}_v\left(t\right)={\mathit{v}}^{\mathrm{G}}-{\mathit{v}}_B^{\mathrm{G}}={\mathit{H}}_v\left(t\right)\mathit{X}\left(t\right)+{\mathit{V}}_v\left(t\right)$$

(19)

where ${\mathit{Y}}_v(t)$ is the measurement vector, ${\mathit{H}}_v(t)=[\begin{array}{ccc}{\mathit{0}}_{3\times 3}\hfill & \mathrm{diag}[\begin{array}{ccc}1\hfill & 1\hfill & 1\hfill \end{array}]\hfill & {\mathit{0}}_{3\times 9}\hfill \end{array}]$ is the measurement matrix, and ${\mathit{V}}_v(t)$ is the measurement noise vector.

The BDS measurement noise characteristics in the polar regions and the design of Kalman filter based on the Huber estimation

The standard Kalman filter is a minimum $l_2$ norm estimator. It requires accurately known system model and Gaussian noise. In fact, BDS measurement noise exhibits a non-Gaussian distribution in the complex polar environment, and the estimation accuracy of standard Kalman filter will be severely degraded. It is necessary to analyze the BDS noise characteristics in the polar regions and research the applicable filter.

Analysis of the BDS measurement noise characteristics in the polar regions

The sky plot of BDS satellites observed from the North Pole is shown in Figure 2. ²³ None of the BDS satellites passes through the polar regions. As a result, the satellite altitude angle in high latitude regions is significantly lower than in other regions. This section will provide a detailed analysis.

Figure 2.

Sky plot of BDS satellites observed from the North Pole (2020/10/30/13:00 UTC).

GEO satellites become invisible and ineffective above 65° north latitude. The maximum altitude angle of MEO and IGSO satellites decreases gradually as latitude increases. Each IGSO satellite has an 8-shaped ground track, which is distributed in the latitude range of about ±55°. Compared with MEO satellites, IGSO satellites have longer visibility and larger maximum altitude angle in high latitude areas. Therefore, IGSO satellite is taken as an example to analyze the satellite altitude angle in the polar regions.

The schematic diagram of IGSO satellite altitude angle in the polar regions is shown in Figure 3. “h,” representing the orbital altitude of IGSO satellite, is about 35,786 km. “R,” representing the radius of the earth, is about 6371 km. From the geometric relationship in Figure 3

$$\beta =\arctan \frac{d_2}{d_1}=\arctan \frac{(h+R)\sin {55}^{°}-R}{(h+R)\cos {55}^{°}}\approx {50}^{°}$$

(20)

Figure 3.

Schematic diagram of IGSO satellite altitude angle in the polar regions.

The maximum altitude angle $\beta$ of BDS satellites in the polar regions is about 50°, which is much lower than that in the middle and low latitude areas. Due to the influence of lower satellite altitude angle and variability of surface materials, the multipath effect of BDS is obvious in high latitude areas,^14–16,18 which will result in measurement noise with non-Gaussian distributions.^24,25 Heavy-tailed distribution imitated by compound Gaussian distribution is one of the typical non-Gaussian distributions, and it is widely used in research works.^25–27 Therefore, we employ the heavy-tail distribution imitated by compound Gaussian distribution to study the polar alignment. Its probability density function is

$$f\left(\omega \right)=\frac{1-\epsilon }{{\sigma }_1\sqrt{2\pi }}\exp (-\frac{{\omega }^2}{2{\sigma }_1^2})+\frac{\epsilon }{{\sigma }_2\sqrt{2\pi }}\exp (-\frac{{\omega }^2}{2{\sigma }_2^2})$$

(21)

where ${\sigma }_1,{\sigma }_2$ are the standard deviations of two Gaussian noise. The former term in equation (21) represents the contaminated Gaussian distribution, and the latter term is the contamination distribution. The parameter $\epsilon$ represents the contamination level, whose value is in [0,1]. We take ${\sigma }_1=1$ , ${\sigma }_2=3.5$ , ${\epsilon }_1=0.1$ , and ${\epsilon }_2=0.3$ as an example, and simulate with 5000 random points. The probability density curves of noises with different $\epsilon$ are depicted in Figure 4.

Figure 4.

The probability density curves of noises with different $\epsilon$ .

As shown in Figure 4, compared with standard Gaussian noise, the probability density curve of compound Gaussian noise is flatter and the probability of distribution at the tail is larger. This means that the probability of noise distribution beyond $3\sigma$ and occurrence of outliers will increase in the polar regions, which will degrade the estimation accuracy of standard Kalman filter.

The design of Kalman filter based on the Huber estimation

Huber estimation filter^28,29 is a robust filter method based on the Huber estimation, and it is a generalized maximum likelihood estimation algorithm. It is a minimum $l_1/l_2$ mixed norm estimator, while minimizing the maximum asymptotic estimation variance under compound Gaussian noise, the estimation efficiency of $l_2$ norm under Gaussian noise is also maintained. It can be used for state estimation under non-Gaussian noise. In this section, a Kalman filter based on the Huber estimation is designed for non-Gaussian measurement noise in polar BDS/SINS integrated alignment.

From the equations (17) and (19), the discrete dynamic model can be obtained

$$\{\begin{array}{c}{\mathit{X}}_{k+1}={\mathit{\Phi }}_{k+1,k}{\mathit{X}}_k+{\mathit{\Gamma }}_k{\mathit{W}}_k\hfill \\ \multicolumn{1}{c}{{\mathbf{Y}}_k={\mathit{H}}_k{\mathit{X}}_k+{\mathit{V}}_k}\end{array}$$

(22)

The generalized maximum likelihood estimation proposed by Huber requires that the residuals are independent. In order to combine the Kalman filter with Huber’s robust estimation strategy, the measurement update should be recast as a linear regression problem between the measurement and the state prediction. If the true value of the state is ${\mathit{X}}_k$ and the one-step prediction of the state is ${\mathit{X}}_{k,k-1}$ , the state prediction error can be written as ${\mathit{\delta }}_k={\mathit{X}}_k-{\mathit{X}}_{k,k-1}$ . Then, the system model can be written as

$$\left[\begin{array}{c}{\mathit{Y}}_k\hfill \\ \multicolumn{1}{c}{{\mathit{X}}_{k,k-1}}\end{array}\right]=\left[\begin{array}{c}{\mathit{H}}_k\hfill \\ \multicolumn{1}{c}{\mathit{I}}\end{array}\right]{\mathit{X}}_k+\left[\begin{array}{c}{\mathit{V}}_k\hfill \\ \multicolumn{1}{c}{-{\mathit{\delta }}_k}\end{array}\right]$$

(23)

where $\mathit{I}$ is the identity matrix.

Defining the decoupling matrix

$${\mathit{S}}_k=\left[\begin{array}{cc}{\mathit{R}}_k\hfill & \mathit{0}\hfill \\ \multicolumn{1}{c}{\mathit{0}}& {\mathit{P}}_{k,k-1}\hfill \end{array}\right]$$

(24)

where ${\mathit{R}}_k,{\mathit{P}}_{k,k-1}$ represent the measurement noise covariance matrix and the one-step prediction variance matrix, respectively.

In order to solve the linear regression problem via maximum likelihood estimation, the residuals’ vector need to be decoupled. Multiply both sides of equation (23) by ${\mathit{S}}_k^{-1/2}$ , the linear regression problem is transformed to

$${\mathit{Z}}_k={\mathit{M}}_k{\mathit{X}}_k+{\mathit{\xi }}_k$$

(25)

where

$${\mathit{Z}}_k={\mathit{S}}_k^{-1/2}\left[\begin{array}{c}{\mathit{Y}}_k\hfill \\ \multicolumn{1}{c}{{\mathit{X}}_{k,k-1}}\end{array}\right]$$

(26)

$${\mathit{M}}_k={\mathit{S}}_k^{-1/2}\left[\begin{array}{c}{\mathit{H}}_k\hfill \\ \multicolumn{1}{c}{\mathit{I}}\end{array}\right]$$

(27)

$${\mathit{\xi }}_k={\mathit{S}}_k^{-1/2}\left[\begin{array}{c}{\mathit{V}}_k\hfill \\ \multicolumn{1}{c}{-{\mathit{\delta }}_k}\end{array}\right]$$

(28)

The covariance of ${\mathit{\xi }}_k$ is the identity matrix, which means the transformed residuals are independent and consistent with the requirement of Huber’s estimation. Therefore, the generalized maximum likelihood estimation proposed by Huber can be employed to solve the problem.

The measurement update of Kalman filter based on the Huber estimation can be solved by minimizing the cost function

$$\mathit{J}\left({\mathit{X}}_k\right)=\sum _{i=1}^m\left(\rho \left({\zeta }_i\right)\right)$$

(29)

where ${\zeta }_i$ refers to the ith component of the residuals’ vector, and the score function $\rho$ is defined as

$$\rho \left({\zeta }_i\right)=\{\begin{array}{c}\frac{{\zeta }_i^2}{2},\left|{\zeta }_i\right|<\gamma \hfill \\ \gamma \left|{\zeta }_i\right|-\frac{{\gamma }^2}{2},\left|{\zeta }_i\right|\ge \gamma \hfill \end{array}$$

(30)

where $\rho (·)$ is the $l_1/l_2$ mixed norm, and $\gamma$ is a tuning parameter which was set to 1.345 in practice. ³⁰ By adopting the cost function given above, the state estimation can minimize the maximum asymptotic estimation variance for the measurement noise with non-Gaussian characteristics.

The solution of the linear regression problem is determined from the derivative of the cost function

$$\frac{\partial \mathit{J}\left({\mathit{X}}_k\right)}{\partial {\mathit{X}}_k}=\sum _{i=1}^m\frac{\rho \left({\zeta }_i\right)}{\partial {\zeta }_i}\frac{\partial {\zeta }_i}{\partial {\mathit{X}}_k}=\sum _{i=1}^m\varphi \left({\zeta }_i\right)\frac{\partial {\zeta }_i}{\partial {\mathit{X}}_k}=0$$

(31)

where

$$\varphi \left({\zeta }_i\right)=\rho \prime \left({\zeta }_i\right)=\{\begin{array}{c}\gamma ,{\zeta }_i\ge \gamma \hfill \\ {\zeta }_i,\left|{\zeta }_i\right|<\gamma \hfill \\ -\gamma ,{\zeta }_i\le -\gamma \hfill \end{array}$$

(32)

Define the weight function

$$\phi \left({\zeta }_i\right)=\frac{\varphi \left({\zeta }_i\right)}{{\zeta }_i}=\{\begin{array}{c}1,\left|{\zeta }_i\right|<\gamma \hfill \\ \multicolumn{1}{c}{\frac{\gamma }{\left|{\zeta }_i\right|},\left|{\zeta }_i\right|\ge \gamma }\end{array}$$

(33)

Define the weight matrix

$$\mathit{\psi }=\mathrm{diag}\left[\phi \left({\zeta }_i\right)\right]=\left[\begin{array}{cc}{\mathit{\psi }}_y\hfill & \mathit{0}\hfill \\ \multicolumn{1}{c}{\mathit{0}}& {\mathit{\psi }}_x\hfill \end{array}\right]$$

and substitute it into equation (31). Then, the modified one-step prediction variance matrix ${\tilde{\mathit{P}}}_{k,k-1}$ and the measurement noise covariance ${\tilde{\mathit{R}}}_k$ in Kalman filter based on the Huber estimation can be obtained

$$\{\begin{array}{c}{\tilde{\mathit{P}}}_{k,k-1}={\mathit{P}}_{k,k-1}^{1/2}{\mathit{\psi }}_x^{-1}{\left({\mathit{P}}_{k,k-1}^{1/2}\right)}^{\mathrm{T}}\hfill \\ {\tilde{\mathit{R}}}_k={\mathit{R}}_k^{1/2}{\mathit{\psi }}_y^{-1}{\left({\mathit{R}}_k^{1/2}\right)}^{\mathrm{T}}\hfill \end{array}$$

(34)

Simulation and experiment results

A new polar alignment algorithm based on the Huber estimation filter is proposed in this article. Simulation and experiment are conducted to verify the effectiveness of our approach.

Simulation results and analyses

The latitude and the longitude of the initial location are set as 87° and 126°. The simulation time and the filter period are 200 and 1 s, respectively. The aircraft flies straight in the north direction. The initial velocity is set as 100 m/s and the acceleration is set as 3 m/s² in 10–50 s. The experimental trajectory is shown in Figure 5.

Figure 5.

The experimental trajectory.

The IMU (Inertial Measurement Unit) of SINS consists of three-axis gyroscopes and accelerometers and it is fixed directly to the aircraft. The working environment is harsh. Therefore, the first-order Markov drift of the gyroscopes and accelerometers should be considered. The other parameters for simulation can be set as listed in Table 1.

Table 1.

Parameters for simulation.

Parameters	Value
Initial attitude error (°)	[0.5, 0.5, 0.5]
Gyro constant drift (°)/h	[0.005, 0.005, 0.005]
Gyro random drift (°)/h	[0.005, 0.005, 0.005]
Correlation time of gyro first-order Markov drift (s)	[3600, 3600, 3600]
White noise variance of gyro first-order Markov drift (°)/h	[0.005, 0.005, 0.005]
Accelerometer random bias (g0)	[1 × 10−5, 1 × 10−5, 1 × 10−5]
Correlation time of accelerometer first-order Markov drift (s)	[1800, 1800, 1800]
White noise variance of accelerometer first-order Markov drift (g0)	[1 × 10−5, 1 × 10−5, 1 × 10−5]

$g_0$ is the value of the gravitational acceleration.

The initial state estimation covariance matrix ${\mathit{P}}_0$ , the system noise covariance matrix $\mathit{Q}$ , and the measurement noise covariance matrix $\mathit{R}$ are set as follows

$$\begin{array}{c}{\mathit{P}}_0=\mathrm{diag}\left\{\begin{array}{c}{(150\pi /3600/180\mathrm{rad})}^2,{(150\pi /3600/180\mathrm{rad})}^2,{(300\pi /3600/180\mathrm{rad})}^2,{(0.6\mathrm{m}/\mathrm{s})}^2,{(0.6\mathrm{m}/\mathrm{s})}^2\hfill \\ \multicolumn{1}{c}{{(0.6\mathrm{m}/\mathrm{s})}^2,{(0.005\pi /3600/180\mathrm{rad}/\mathrm{s})}^2,{(0.005\pi /3600/180\mathrm{rad}/\mathrm{s})}^2,{(0.005\pi /3600/180\mathrm{rad}/\mathrm{s})}^2}\\ \multicolumn{1}{c}{{(0.005\pi /3600/180\mathrm{rad}/\mathrm{s})}^2,{(0.005\pi /3600/180\mathrm{rad}/\mathrm{s})}^2,{(0.005\pi /3600/180\mathrm{rad}/\mathrm{s})}^2}\\ \multicolumn{1}{c}{{(1\times {10}^{-5}{g}_0\mathrm{m}/{\mathrm{s}}^2)}^2,{(1\times {10}^{-5}{g}_0\mathrm{m}/{\mathrm{s}}^2)}^2,{(1\times {10}^{-5}{g}_0\mathrm{m}/{\mathrm{s}}^2)}^2}\end{array}\right\}\hfill \\ \mathit{Q}=\mathrm{diag}\left\{\begin{array}{c}{(0.005\pi /3600/180\mathrm{rad}/\mathrm{s})}^2,{(0.005\pi /3600/180\mathrm{rad}/\mathrm{s})}^2,{(0.005\pi /3600/180\mathrm{rad}/\mathrm{s})}^2\hfill \\ \multicolumn{1}{c}{{(\sqrt{\frac{2T_D}{3600}}\times 0.005\pi /3600/180\mathrm{rad}/\mathrm{s})}^2,{(\sqrt{\frac{2T_D}{3600}}\times 0.005\pi /3600/180\mathrm{rad}/\mathrm{s})}^2}\\ \multicolumn{1}{c}{{(\sqrt{\frac{2T_D}{3600}}\times 0.005\pi /3600/180\mathrm{rad}/\mathrm{s})}^2,{(\sqrt{\frac{2T_D}{1800}}\times 1\times {10}^{-5}{g}_0\mathrm{m}/{\mathrm{s}}^2)}^2}\\ \multicolumn{1}{c}{{(\sqrt{\frac{2T_D}{1800}}\times 1\times {10}^{-5}{g}_0\mathrm{m}/{\mathrm{s}}^2)}^2,{(\sqrt{\frac{2T_D}{1800}}\times 1\times {10}^{-5}{g}_0\mathrm{m}/{\mathrm{s}}^2)}^2}\end{array}\right\}\hfill \\ \mathit{R}=\mathrm{diag}\{{(0.2\mathrm{m}/\mathrm{s})}^2,{(0.2\mathrm{m}/\mathrm{s})}^2,{(0.2\mathrm{m}/\mathrm{s})}^2\}\hfill \end{array}$$

where $T_D$ is the discrete periodic.

To compare the performance of standard Kalman filter and Kalman filter based on the Huber estimation under non-Gaussian noise, assuming that BDS velocity measurement noise distribution follows the compound Gaussian distribution analyzed in section “Analysis of the BDS measurement noise characteristics in the polar regions”. Take ${\sigma }_1=0.2$ , ${\sigma }_2=20$ , and $\epsilon =0.3$

$${\mathit{V}}_{\mathrm{BDS}}:(1-\epsilon )N(0,{\sigma }_1^2)+\epsilon N(0,{\sigma }_2^2)$$

(35)

The polar alignment algorithm based on the Huber estimation filter with the aid of BDS proposed in this article is called the Proposed Algorithm. The polar alignment algorithm based on the standard Kalman filter with the aid of BDS is called the Conventional Algorithm. The attitude error between the corrected attitude and the true attitude is analyzed. The simulation results are expressed for comparison shown as Figure 6.

Figure 6.

Attitude error curves of two alignment algorithms.

The histograms of the attitude error mean and standard deviation of the two alignment methods over a period of 150–200 s are shown in Figure 7.

Figure 7.

Performance comparison histograms of two alignment algorithms: (a) mean of attitude error and (b) standard deviation of attitude error.

As shown in Figure 6, the alignment results of conventional algorithm are as follows: the roll error remains at less than 6′ after 150 s, the pitch error remains at less than 10′ after 150 s, and the yaw error remains at less than 30′ after 150 s. And the alignment results of proposed algorithm are as follows: the roll error remains at less than 6′ after 70 s, the pitch error remains at less than 10′ after 80 s, and the yaw error remains at less than 10′ after 85 s. It takes less time to achieve the same alignment accuracy using proposed algorithm than using conventional algorithm.

Combining the mean and standard deviation of the attitude error given in Figure 7, in the same period of time (150–200 s), proposed algorithm achieves better alignment accuracy. The mean of roll error, pitch error, and yaw error of proposed algorithm respectively reduces to about 79%, 15%, and 27% compared with conventional algorithm.

Based on the results shown in Figures 6 and 7, the following conclusions about the simulation results can be made. Under heavy-tailed measurement noise, the performance of proposed algorithm is significantly better than conventional algorithm. The standard Kalman filter is sensitive to heavy-tailed noise and the alignment performance is severely degraded. Proposed algorithm achieves better state estimation, illustrating its robustness to non-Gaussian noise.

Experiment results and analyses

Considering the geography restrictions, the experiment is conducted in the form of semi-physical simulation.

The IMU information include ${\hat{\mathit{\omega }}}_{\mathrm{ib}}^{\mathrm{b}}$ and ${\hat{\mathit{f}}}^{\mathrm{b}}$ . ${\hat{\mathit{\omega }}}_{\mathrm{ib}}^{\mathrm{b}}$ is the practical angular velocity measured by gyroscopes, and it is composed of the true angular velocity ${\mathit{\omega }}_{\mathrm{ib}}^{\mathrm{b}}$ and the gyro drift $\delta {\mathit{\omega }}_{\mathrm{ib}}^{\mathrm{b}}$ . ${\hat{\mathit{f}}}^{\mathrm{b}}$ is the practical specific force measured by accelerometer, and it is composed of the true specific force ${\mathit{f}}^{\mathrm{b}}$ and the accelerometer drift $\delta {\mathit{f}}^{\mathrm{b}}$

$$\{\begin{array}{c}{\hat{\mathit{\omega }}}_{\mathrm{ib}}^{\mathrm{b}}={\mathit{\omega }}_{\mathrm{ib}}^{\mathrm{b}}+\delta {\mathit{\omega }}_{\mathrm{ib}}^{\mathrm{b}}\hfill \\ \multicolumn{1}{c}{{\hat{\mathit{f}}}^{\mathrm{b}}={\mathit{f}}^{\mathrm{b}}+\delta {\mathit{f}}^{\mathrm{b}}}\end{array}$$

(36)

where the superscript b represents the b frame of SINS.

In the simulation and the experiment, ${\mathit{\omega }}_{\mathrm{ib}}^{\mathrm{b}}$ and ${\mathit{f}}^{\mathrm{b}}$ can be gained by simulation in the polar regions once the attitude variations and maneuvers of the aircraft are confirmed.^12,14,15 The gyro drift and the accelerometer drift would not change in different regions because they are inherent to SINS. Therefore, the gyro drift $\delta {\mathit{\omega }}_{\mathrm{ib}}^{\mathrm{b}}$ and the accelerometer drift $\delta {\mathit{f}}^{\mathrm{b}}$ for polar experiment can be extracted from actual measurement data in the non-polar regions. The schematic diagram of data composition for polar experiment is illustrated in Figure 8.

Figure 8.

The schematic diagram of data composition for polar experiment.

The actual measurement data can be obtained through the experiment conducted in a laboratory in the non-polar regions (31.937°N, 118.798°E). The SINS is installed on a turntable, which can provide a high-precision three-axis rotary movement. The parameters of SINS used in the experiment are listed in Table 2.

Table 2.

The parameters of SINS used in the experiment.

Parameters	Value
Gyro bias stability	0.001 (°)/h
Gyro random walk	0.00015 (°)/sqrt(h)
Accelerometer bias stability	1 × 10−5 (g0)

According to the calculation process of gyro random walk and bias stability performance, the random walk index is only related to white noise, and the bias stability index is related to white noise and first-order Markov drift. ³¹ Thus, BDS/SINS integrated alignment system state equation is the same as that in section “BDS/SINS integrated alignment state equation.” The experiment results are expressed for comparison shown as Figure 9.

Figure 9.

Attitude error curves of two alignment algorithms.

The histograms of the attitude error mean and standard deviation of the two alignment methods over a period of 150–200 s are shown in Figure 10.

Figure 10.

Performance comparison histograms of two alignment algorithms: (a) mean of attitude error and (b) standard deviation of attitude error.

As shown in Figure 9, the alignment results of conventional algorithm are as follows: the roll error remains at less than 6′ after 150 s, the pitch error remains at less than 10′ after 150 s, and the yaw error remains at less than 50′ after 150 s. And the alignment results of proposed algorithm are as follows: the roll error remains at less than 6′ after 115 s, the pitch error remains at less than 10′ after 90 s, and the yaw error remains at less than 10′ after 110 s. It takes less time to achieve the same alignment accuracy using proposed algorithm than using conventional algorithm.

Combining the mean and standard deviation of the attitude error given in Figure 10, in the same period of time (150–200 s), proposed algorithm achieves better alignment accuracy. The mean of roll error, pitch error, and yaw error of proposed algorithm respectively reduces to about 60%, 27%, and 13% compared with conventional algorithm.

As shown in Figures 9 and 10, the semi-physical simulation results further validate that proposed algorithm has better performance than conventional algorithm. The semi-physical simulation results coincide with the simulation results.

Discussion

As shown in the simulation and experiment results, the new polar alignment algorithm based on the Huber estimation filter proposed in this article is much superior to that based on the Kalman filter. The main contributions and the proposed algorithm are discussed as follows:

1. Since the traditional north-oriented frame loses its effectiveness for polar navigation, this article employs the grid SINS mechanization for aircrafts. SINS mechanization and navigation error equation in the grid frame are represented. Furthermore, the dynamic model and the measurement model are established.

2. For aircrafts equipped with BDS/SINS integrated navigation system, BDS information can be used to achieve autonomous alignment. However, due to the complex polar environment, the BDS noise characteristics are abnormal. Therefore, the BDS noise characteristics in the polar regions are analyzed and heavy-tailed characteristics are simulated, respectively.

3. Since BDS measurement noise exhibits a non-Gaussian distribution in the complex polar environment, the estimation accuracy of standard Kalman filter is severely degraded. Kalman filter based on the Huber estimation is a minimum $l_1/l_2$ mixed norm estimator, while minimizing the maximum asymptotic estimation variance under compound Gaussian distribution, the estimation efficiency of the $l_2$ norm under Gaussian noise is maintained. Therefore, the proposed algorithm suppresses the outliers significantly better than the one based on the standard Kalman filter and has better robustness under non-Gaussian noise.

The above analysis indicates that the polar alignment algorithm based on the Huber estimation filter is effective, which takes full account of the characteristics of BDS and ensures the accuracy and robustness of the polar alignment.

Conclusion

BDS measurement noise often exhibits a non-Gaussian distribution due to the complex polar environment and multipath effect. To improve the estimation accuracy, a new polar alignment algorithm based on the Huber estimation filter is proposed. The dynamic model and the measurement model of BDS/SINS integrated alignment system in the grid frame are derived. The BDS measurement noise characteristics in the polar regions are analyzed and heavy-tailed characteristics are simulated, respectively. Since the estimation accuracy of standard Kalman filter can be severely degraded under non-Gaussian noise, a Kalman filter based on the Huber estimation is designed. The simulation and experiment results demonstrate that the proposed algorithm has better robustness under non-Gaussian noise, and it is effective in the polar regions. By employing the proposed algorithm, the rapidity and accuracy of the alignment can be improved. This method can be easily extended to other GNSS.