1. Introduction

Celestial navigation is an autonomous navigation technology based on the inertial frame of indestructible natural celestial bodies. It has the advantages of no electromagnetic interference, no dependence on external facilities, and no accumulation of navigation errors over time (Wang, Reference Wang2007). Although global navigation satellite systems (GNSS) are improving, they are vulnerable to electromagnetic interference, so celestial navigation is still regarded as an important maritime navigation technology and is widely used (Vulfovich and Fogilev, Reference Vulfovich and Fogilev2010).

Celestial navigation originated from ship navigation. In 1875, the altitude difference method was proposed by Saint-Hilaire, which marked the maturity of modern celestial navigation theory. As a main navigate mode in navigation for hundreds of years, the altitude difference method observes celestial bodies by using a sextant and calculates astronomical vessel position (AVP) according to the difference between measured altitude and the calculated altitude of celestial bodies (Pepperday, Reference Pepperday1992). In practice, the altitude difference method has been gradually improved. In the traditional altitude difference method based on position lines, it is necessary to have an ‘estimated position’ to start the calculation, which sometimes causes difficulties (Spencer, Reference Spencer1990). To solve this problem, a mathematical method was presented to obtain an AVP directly as the intersection point of the position circles of the celestial bodies observed with a sextant (Chiesa and Chiesa, Reference Chiesa and Chiesa1990). On this basis, an astronomical running fix algorithm was introduced to solve the problem that the second observation was taken after the observer has traveled an appreciable distance along a rhumb line (Williams, Reference Williams1990;  Matti, Reference Matti1990). In 1994, the experimental position accuracy of celestial fixes were given by comparing with global positioning system (GPS) results, and the conclusion was that the standard deviation of celestial AVP observed by sextant was about 1525 m (Bovens, Reference Bovens1994). To improve computing speed and convenience, the simultaneous equal-altitude equations method (SEEM), based on the trigonometric equation solving technique, was developed to directly calculate the AVP without an additional graphical procedure (Hsu et al., Reference Hsu, Chen and Chang2003, Reference Hsu, Chen and Chang2005), and the corresponding computation programs was proposed (Chen and Hsieh, Reference Chen and Hsieh2011). However, these algorithms mostly emphasise AVP determination only, and could not realise the real-time vehicle position, speed and course update simultaneously. An alternative celestial navigation method was invented, referring to the principle of satellite orbit determination, which provided the means for determining the values of all four parameters, describing a vessel's rhumb-line track across the Earth: latitude and longitude at a specific time, course and speed (Kaplan, Reference Kaplan1995). Later, the approach was extended to cover the case where observations were taken over several legs of a voyage, and a single least squares (LS) solution for all sailing parameters was developed (Kaplan, Reference Kaplan1996). This algorithm was then added to the Navigation Sensor System Interface (NAVSSI) for Navy shipboard use (Kaplan, Reference Kaplan1999).

With advances is photoelectric technology and digital image processing technology, a star sensor which can shoot multiple stars simultaneously to realise celestial navigation was developed. Accordingly, the process of celestial navigation also developed from artificial observation, mapping and calculating to the mode with completely automatic data acquiring, processing and data-out. From this brief history of marine celestial navigation system development, it can be seen that the equipment has experienced several stages of development from sextant, star tracker, to star sensor, and data processing translated from artificially calculating to automatically calculating. To improve the visibility of stars at sea, wide-angle lens camera which could simultaneously image all the stars above the horizon has been used in maritime celestial navigation (Li et al., Reference Li, Zheng, Zhang, Yuan, Lian and Zhou2014). Meanwhile, a wide-angle lens camera could also be used to image the moon to determine vessel heading when stars are invisible (Pu et al., Reference Pu, Li, Zheng and Zhan2018), where the corresponding algorithm used to determine the astronomical vessel position and heading are both based on the LS algorithm.

Although several real-time algorithms for maritime celestial navigation have been presented (Kaplan, Reference Kaplan1999), they usually have two problems. Firstly, all the celestial observations are processed with equal weight directedly without considering their different precisions, which limits the improvement of celestial navigation accuracy. Secondly, the ship is often in a nonstationary running state because of speed and course change, and other situations. To solve these problems, this paper proposes a robust adaptive maritime celestial navigation algorithm. First, the algorithm weights all the observed values according to the observation accuracy, which take full advantage of high-precision observations, limit the use of low-precision ones and exclude outliers, then more reliable, effective and meaningful parameter estimations can be obtained. Secondly, by introducing the adaptive factor, the automatic adjustment of the dynamic balance between astronomical observation values and the ship motion prediction values is realised, as well as the navigation error increase problem of the ship's nonstationary running state is solved.

2. Corresponding coordinate system

2.1 Equatorial coordinate system

The celestial coordinate system is a spherical coordinate system introduced to accurately represent the projection position of celestial bodies on the celestial sphere. The intersection line obtained by magnifying the equator of the earth to the celestial sphere is called the ‘celestial equator.’ The equatorial coordinate system is based on the celestial equator. If the center of the earth is taken as its coordinate origin, it is called the geocentric equatorial coordinate system. The $O{X_E}$-axis in the geocentric equatorial coordinate system $O - {X_E}{Y_E}{Z_E}$ points to the spring equinox $\gamma$, the $O{Z_E}$-axis points to the north pole P, and the $O{Y_E}$-axis is in the celestial equator plane and forms a right-handed coordinate system with the $O{X_E}$-axis and $O{Z_E}$-axis, as shown in Figure 1.

Figure 1. Equatorial coordinate system

2.2 Image coordinate system

The image coordinate system $O - {X_I}{Y_I}$ is based on the image plane. The origin O locates in the lower-left corner of the image. The $O{X_I}$-axis points to the right side of the image. The $O{Y_I}$-axis points to the top of the image, as shown in Figure 2.

Figure 2. Image coordinate system

2.3 Ship body-fixed frame

The ship body–fixed frame $O - {X_B}{Y_B}{Z_B}$ is based on the deck plane of a ship. The $O{X_B}$-axis points to the starboard of the hull. The $O{Y_B}$-axis points to the bow. The $O{Z_B}$-axis is perpendicular to the deck plane of a ship and points upwards, as shown in Figure 3.

Figure 3. Ship body–fixed frame

3. Basic maritime navigation models

3.1 Ship motion equation

The ship motion model could be generally divided into a dynamic model and a kinematic model. The dynamic model is derived based on the force suffered on ships, and the kinematics model describes the trajectory of the ship by using the dynamic position or velocity vector. In the course of ship sailing, the forces exerted on the ship are extremely complicated because both stormy waves and the power provided by the ship itself will affect the ship's motion. So the ship dynamic model is hard to be established accurately while the kinematic model could be established accurately due to the speed and course change is subtle in a short period. For satellites, aircraft, vehicles and other carriers, their state is represented by a three-dimensional (3D) position vector and a 3D velocity vector, but when the ship sailing on the sea surface without calculating the elevation, only two-dimensional (2D) position and 2D speed can express its motion state. The ship position vector r represented by the geographical coordinate system is

(3.1)\begin{equation}r = {({\varphi \textrm{,}\lambda } )^T}\end{equation}

The velocity vector $\dot{r}$ is

(3.2)\begin{equation}\dot{r} = {({{V_N}\textrm{,}{V_E}} )^T}\end{equation}

where $\varphi$ = astronomical latitude, $\lambda$ = astronomical longitude, ${V_N}$ = northward velocity component, and ${V_E}$ = eastward velocity component of the ship. Then the uniform motion model transfer from moment $t - 1$ to moment t could be represented as

(3.3)\begin{equation}\left[ {\begin{matrix} {r\textrm{(}t\textrm{)}}\\ {\dot{r}\textrm{(}t\textrm{)}} \end{matrix}} \right] = \left[ {\begin{matrix} I& {D\varDelta t}\\ 0& I \end{matrix}} \right]\left[ {\begin{matrix} {r\textrm{(}t - 1\textrm{)}}\\ {\dot{r}\textrm{(}t - 1\textrm{)}} \end{matrix}} \right] + W\textrm{(}t\textrm{)}\end{equation}

where, I = unit matrix, $\varDelta t$ = time interval from moment $t - 1$ to moment t, $W\textrm{(}t\textrm{)}$ = noise vector, and D is

(3.4)\begin{equation}D = \left[ {\begin{matrix} {\dfrac{1}{R}}& 0\\ 0& {\dfrac{1}{{R\cos \varphi }}} \end{matrix}} \right]\end{equation}

where R = earth radius. If

(3.5)\begin{align} {X_t} & = \left[ {\begin{matrix} {r\textrm{(}t\textrm{)}}\\ {\dot{r}\textrm{(}t\textrm{)}} \end{matrix}} \right]\end{align}

(3.6)\begin{align}{X_{t - 1}} & = \left[ {\begin{matrix} {r\textrm{(}t - 1\textrm{)}}\\ {\dot{r}\textrm{(}t - 1\textrm{)}} \end{matrix}} \right]\end{align}

(3.7)\begin{align}{\Phi _{t\textrm{,}t - 1}} & = \left[ {\begin{matrix} I& {D\varDelta t}\\ 0& I \end{matrix}} \right]\end{align}

here I = unit matrix, then the state equation of ship motion could be simplified as

(3.8)\begin{equation}{X_t} = {\Phi _{t\textrm{,}t - 1}}{X_{t - 1}} + {W_t}\end{equation}

where ${X_t}$ = ship's state parameter vector at moment t, ${X_{t - 1}}$ = ship's state parameter vector at moment $t - 1$, ${\Phi _{t\textrm{,}t - 1}}$ = state-transition matrix, and ${W_t}$ = Gaussian white-noise vector.

3.2 Celestial observation equation

In addition to the ship motion equation, it is necessary to establish an astronomical observation model for filtering processing. In the course of ship navigation, the celestial zenith distance calculated from the star map is regarded as the observation. The observation vector at moment t is

(3.9)\begin{equation}{L_t} = {[{{z_1}\textrm{,}{z_2} \ldots {z_i} \ldots {z_{{n_t}}}} ]^T}\end{equation}

where ${z_i}$ = celestial zenith distance of the $i$th star in the star map generated at moment t, and ${n_t}$ = number of stars recognised in the star map. Then the observation model of ship celestial navigation can be expressed as (Li et al., Reference Li, Zheng, Zhang, Yuan, Lian and Zhou2014)

(3.10)\begin{equation}{L_t} = {A_t}{X_t} + {e_t}\end{equation}

where ${A_t}$ is a ${n_t} \times 4$ design matrix, and ${e_t}$ is a ${n_t} \times 4$ n observed noise matrix. It is assumed that the ship moves in a straight line at a constant speed in a short time. ${A_t}$ could be expressed as

(3.11)\begin{equation}{A_t} = \left[ {\begin{matrix} {{a_1}}& {{b_1}}& 0& 0\\ {{a_2}}& {{b_2}}& 0& 0\\ \vdots & \vdots & \vdots & \vdots \\ {{a_n}}& {{b_n}}& 0& 0 \end{matrix}} \right]\end{equation}

where

(3.12)\begin{align}{a_i} & ={-} \frac{{\cos \varphi _{0} \sin {\delta _i} - \sin \varphi _{0}\cos\delta_{i}\cos H_{i}}}{{\sqrt {1 - P_{i0}^2} }}\end{align}

(3.13)\begin{align}{b_i} & = \frac{{\cos \varphi _{0}\cos \delta _{i}\sin H_{i}}}{{\sqrt {1 - P_{i0}^2} }}\end{align}

(3.14)\begin{align}{P_{i0}} & = \sin \varphi _{0}\sin \delta _{i} + \cos \varphi _{0}\cos \delta _{i}\cos H_{i}\end{align}

where ${\varphi _0}$ = initial value of ship's latitude, it could be extrapolated from the previous position; ${\delta _i}$ = stars’ declination; ${H_i}$ = stars’ hour angle at the observation moment.

In Equation (3.10), the expectation of ${e_t}$ is zero, the covariance matrix is

(3.15)\begin{equation}{\Sigma _{{e_t}{e_k}}} = \left\{ {\begin{matrix} {{\Sigma _t}}& {t = k}\\ 0& {t \ne k} \end{matrix}} \right.\end{equation}

There are some assumptions: ${e_t}$ = Gaussian white noise sequence, initial value of state ${X_0}$; observation noise ${e_t}$ and model noise ${W_t}$ are uncorrelated; and ${e_t}$ and ${W_t}$ are statistic independent. Therefore

(3.16)\begin{equation}{\Sigma _{{e_i}{W_j}}} = E\textrm{(}{e_i}{W_j}^T\textrm{)} = 0\end{equation}

The optimal value of the ship state parameter could be calculated by using the kinematic model information and observation information.

4. Adaptively robust celestial navigation

Several algorithms can utilize the kinematic model and observation information to estimate the ship state parameters, in which the Kalman filter is the most basic method. To calculate the optimal value of the ship state parameters, the kinematic equation (8) and the observation equation (10) of the ship could be rewritten as the error equation (Li et al., Reference Li, Zheng, Li, Yu and Wang2013)

(4.1)\begin{align}{V_{{{\bar{X}}_t}}} & = {\hat{X}_t} - {\bar{X}_t}\end{align}

(4.2)\begin{align}{V_t} & = {A_t}{\hat{X}_t} - {L_t}\end{align}

where ${\hat{X}_t}$ = ship motion state vector to be evaluated, ${\bar{X}_t}$ = predicted state vector at moment t, ${V_{{{\bar{X}}_t}}}$ = residual vector of ${\bar{X}_t}$, and ${V_t}$ = residual vector of ${L_t}$. The objective function could be constructed by the LS principle as

(4.3)\begin{equation}\Omega \textrm{(}t\textrm{)} = {V_t}^T{P_t}{V_t} + {V_{{{\bar{X}}_t}}}^T{P_{{{\bar{X}}_t}}}{V_{{{\bar{X}}_t}}} = \textrm{minimum}\end{equation}

where ${P_t}$ and ${P_{{{\bar{X}}_t}}}$ = weight matrix of ${L_t}$ and ${\bar{X}_t}$, respectively, and

(4.4)\begin{equation}{P_t} = {\Sigma _t}^{ - 1}\quad {P_{{{\bar{X}}_t}}} = {\Sigma _{{{\bar{X}}_t}}}^{ - 1}\end{equation}

The LS estimation of ship motion state vector ${\hat{X}_t}$ could be obtained according to the objective function as

(4.5)\begin{equation}{\hat{X}_t} = {({{A_t}^T{P_t}{A_t} + {P_{{{\bar{X}}_t}}}} )^{ - 1}}({{A_t}^T{P_t}{L_t} + {P_{{{\bar{X}}_t}}}{{\bar{X}}_t}} )\end{equation}

Therefore, the position and velocity estimations of the ship could be calculated in real time by using the stated formula. However, the standard Kalman filter solution is the LS solution in essence and has two implicit assumptions, one of which is the observed value vector and the state prediction vector are all subject to the normal distribution, and the other one is the ship movement is relatively smooth. In the actual navigation process, both the observation vector and the ship's motion model information may be abnormal, and the ship may also have a large range of nonstationary states, such as variable speed and steering. All of these are difficult to described by a unified model. In this case, it is hard to get high-precision of ship state parameter estimation with a standard Kalman filter and even the filtering divergence problem may arise. Because of these problems, the adaptively robust filtering algorithm could be adopted for processing. The principle of robust estimation is applied to the observed value when the observation is abnormal, and the motion model information is adopted as a unified adaptive factor to adjust its overall contribution to the state parameters when the ship motion model has abnormal errors (Yang et al., Reference Yang, He and Xu2001).

The formula (Equation 4.3) is transformed into an extremum function that can be used to control observation anomalies and abnormal state information as

(4.6)\begin{equation}\Omega \textrm{(}t\textrm{)} = {V_t}^T{\bar{P}_t}{V_t} + {\alpha _t}{V_{{{\bar{X}}_t}}}^T{P_{{{\bar{X}}_t}}}{V_{{{\bar{X}}_t}}} = \textrm{minimum}\end{equation}

In Equation (4.6), $\bar{P}$ is the robust equivalent weight matrix for ${L_t}$, and it is the adaptive estimation of the weight matrix of the celestial zenith observation vector. $\bar{P}$ is a diagonal matrix for celestial navigation observations from a star sensor and its diagonal elements can be expressed as (Li et al., Reference Li, Pu, Zhan, Zheng, Zhang, Wang, Luo and Wang2018)

(4.7)\begin{equation}\bar{p}({{{\tilde{v}}_i}} )= \begin{cases} 1 & {|{{{\tilde{v}}_i}} |\le {k_0}}\\ {\dfrac{{{k_0}}}{{|{{{\tilde{v}}_i}} |}}{{\left( {\dfrac{{{k_1} - |{{{\tilde{v}}_i}} |}}{{{k_1} - {k_0}}}} \right)}^2}} & { {k_0} \le |{{{\tilde{v}}_i}} |\le {k_1}}\\ 0 & {{k_1} \le |{{{\tilde{v}}_i}} |} \end{cases} \end{equation}

where ${k_0}$ and ${k_1}$ = 1⋅5 and 3⋅0, generally; ${\tilde{v}_i}$ = standardised residual

(4.8)\begin{equation}{\tilde{v}_i} = \frac{{{v_i}}}{{{m_{{v_i}}}}}\end{equation}

In Equation (4.8), ${v_i}$ = LS residual of observation, ${m_{{v_i}}}$ = mean square error of ${v_i}$, and is obtained by

(4.9)\begin{equation}{m_{{v_i}}} = {\sigma _0}\sqrt {\frac{1}{{{p_i}}} - {A_i}{{({{A^T}PA} )}^{ - 1}}{A_i}^T}\end{equation}

In Equation (4.9), ${\sigma _0}$ = mean square error of unit weight, which can use theoretical or empirical values, ${p_i}$ = initial right of the observed value, and ${A_i}$ = $i$th line of the matrix A as well as the coefficient vector of the $i$th error equation.

In Equation (4.6), ${\alpha _t}$ is the adaptive factor, and its value range is $0 \le {\alpha _t} \le 1$, which could be expressed as

(4.10)\begin{equation}{\alpha _t} = \begin{cases} 1 & {{m_0} \ge |{\varDelta {{\tilde{X}}_t}} |}\\ {\dfrac{{{m_0}}}{{|{\varDelta {{\tilde{X}}_t}} |}}{{\left( {\dfrac{{{m_1} - |{\varDelta {{\tilde{X}}_t}} |}}{{{m_1} - {m_0}}}} \right)}^2}} & {{m_0} \le |{\varDelta {{\tilde{X}}_t}} |\le {m_1}}\\ 0 & {{m_1} \le |{\varDelta {{\tilde{X}}_t}} |} \end{cases} \end{equation}

where the value of ${m_0}$ can be 1⋅0~1⋅5 and ${m_1}$ can be 3⋅0~8⋅0 according to actual situations. In Equation (4.10), $\varDelta {\tilde{X}_t}$ is

(4.11)\begin{equation}\varDelta {\tilde{X}_t} = \frac{{||{{{\tilde{X}}_t} - {{\bar{X}}_t}} ||}}{{\sqrt {tr{\Sigma _{{{\bar{X}}_t}}}} }}\end{equation}

where ${\Sigma _{{{\bar{X}}_t}}}$ = covariance matrix of ${\bar{X}_t}$, and ${\tilde{X}_t}$ = robust estimation solution of current state parameters of the ship motion state that based on the observation, namely

(4.12)\begin{equation}{\tilde{X}_t} = {({{A_t}^T{{\bar{P}}_t}{A_t}} )^{ - 1}}{A_t}^T{\bar{P}_t}{L_t}\end{equation}

where the equivalent weight matrix ${\bar{P}_t}$ of observations could be determined by Equation (4.7). Substituting the error equation Equation (4.1) and Equation (4.2) into Equation (4.6), the adaptively robust filtering solution of the ship motion state parameter vector is obtained as

(4.13)\begin{equation}{\hat{X}_t} = {({{A_t}^T{{\bar{P}}_t}{A_t} + {\alpha_t}{P_{{{\bar{X}}_t}}}} )^{ - 1}}({{A_t}^T{{\bar{P}}_t}{L_t} + {\alpha_t}{P_{{{\bar{X}}_t}}}{{\bar{X}}_t}} )\end{equation}

It can be seen from the calculation process of the adaptively robust filtering that if the celestial zenith distance observation of a celestial body contains outliers, then the corresponding equivalent weight elements will decrease, thus the influence of the observed outliers on the estimation of ship state parameters could be controlled. If the kinematics model of the ship is disturbed or the ship is undergoing nonstationary processes, such as speed change and steering, then the corresponding adaptive factor will decrease, then the influence of the state model information anomalies on the valuation of ship state parameters can be controlled. Therefore, the robust adaptive filtering algorithm can adjust itself automatically to obtain the high-precision ship motion-state parameters. The recursive solution of the state vector by the identity transformation of a matrix is

(4.14)\begin{equation}{\hat{X}_t} = ({I - {{\bar{K}}_t}{A_t}} ){\bar{X}_t}\textrm{ + }{\bar{K}_t}{L_t}\end{equation}

where I = unit matrix, and ${\bar{K}_t}$ = filter gain matrix, which could be expressed as

(4.15)\begin{equation}{\bar{K}_t} = \frac{1}{{{\alpha _t}}}{\Sigma _{{{\bar{X}}_t}}}{A_t}^T{\left( {\frac{1}{{{\alpha_t}}}{A_t}{\Sigma _{{{\bar{X}}_t}}}{A_t}^T + {{\bar{\Sigma }}_t}} \right)^{ - 1}}\end{equation}

In Equation (4.15), ${\Sigma _{{{\bar{X}}_t}}}$ = covariance matrix of ${\bar{X}_t}$, ${\bar{\Sigma }_t}$ = covariance matrix of ${\tilde{v}_i}$, then the estimated value ${\hat{X}_t}$ of the ship motion state parameters could be obtained through the stated process including two position parameters ($\varphi_{t},\lambda_{t}$) and speed parameters ($V_{Nt},V_{Et}$). The speed of the ship can be calculated accordingly as

(4.16)\begin{equation}{V_t} = \sqrt {{V_{Nt}}^2 + {V_{Et}}^2}\end{equation}

Based on the true north direction, the ship heading can be calculated as

(4.17)\begin{equation}H = \textrm{arctan}\frac{{{V_{Nt}}}}{{{V_{Et}}}}\end{equation}

In Equation (4.17), the signs of the northward velocity ${V_{Nt}}$ and the eastward velocity ${V_{Et}}$ are respected to determine the quadrant of the heading angle H. So far, the navigation parameters, such as the position, course and speed of the ship, have been obtained and the navigation trajectory of the ship could be obtained if continuous navigation is carried out for a period of time.

5. Experimental analysis

To verify the correctness and effectiveness of the adaptively robust filter algorithm proposed in this paper, ship-based celestial navigation experiment data is used for analysis. On 13 October 2017, a celestial navigation experiment was performed at the Bohai Sea by utilising a super-wide-angle-lens camera (Li et al., Reference Li, Zheng, Zhang, Yuan, Lian and Zhou2014). In the experiment, the celestial navigation system was anchored at the top deck edge on the Yongxingdao passenger liner. The size of the vessel is 167⋅5 meters in length and 25⋅2 meters in width, and the displacement of the vessel is about 23,000 tons; therefore, the ship's swing was so small that it was almost negligible. The vessel is shown in Figure 4.

Figure 4. Yongxingdao passenger liner

We set the exposure time of the fish camera as 0⋅3 s and continuously shooting 500 images. Because the transit time for each image by the USB 3⋅0 interface is approximately 5⋅5 s, it took a total of about 50 min. A typical fisheye star image is shown in Figure 5.

Figure 5. Typical fisheye star image

Except for stars, the moon, cloud, and light were all photographed in the same image, which results in a very complex background for star extraction and recognition. In Figure 5, the bright area in the upper-right corner of the image are the pixels projected by the moon and clouds, the bright part of the image in the surrounding area is projected by some lights on the ship, the bright pixels in the sky area are projected by the stars. Although dozens of stars can be extracted from the image, there is inevitably a lot of noise and errors in the coordinates of these stars. Therefore, these errors eventually lead to a large number of random errors, systematic errors and gross errors in star-zenith distance observations.

After star point extraction, identification and zenith distance calculation of the above 500 images, the LS, extended Kalman filtering (EKF) and robust adaptive extended Kalman filtering (RAEKF) algorithms proposed in this paper were used to calculate, respectively. In addition, GNSS observations were used as references to evaluate the accuracy and effectiveness of the three algorithms. After polynomial fitting analysis, GNSS track residuals were generally within the range of ±30 m, and the root mean square error (RMSE) of 9⋅89 m.

Figure 6 shows the position relationship between the AVP calculated by the three algorithms and the GNSS reference track, and Figure 7 shows the variation of the position errors calculated by the three algorithms with coordinated universal time (UTC). As can be seen from the two figures, the LS results have the highest dispersion, and the RAEKF results are closest to GNSS reference values; this is mainly due to the robustness of the equivalent weight and adaptive factor.

Figure 6. Vessel track and position results

Figure 7. Vessel position errors

Because the LS algorithm uses a single image to solve the problem, only the AVP can be obtained, but will not include the speed and course results (Figures 8, 9). The GNSS speed (Figure 8) and course (Figure 9) were obtained by a single antenna receiver based on pseudo-range navigation solutions, which have a large vibration amplitude. Therefore, the 8th polynomial fitting value of them can be used as the reference. The difference between the speed and course results and their GNSS reference value is shown in Figure 10 and Figure 11, respectively.

Figure 8. Vessel speed results

Figure 9. Vessel course results

Figure 10. Vessel speed errors

Figure 11. Vessel course errors

After statistics, the RMSE of the three algorithms are shown in Table 1.

Table 1. ME and RMSE of three algorithms

After statistics, the mean error (ME) and mean squared error (MSE) of the three algorithms are shown in Table 1.

It can be seen from these figures and Table 1 that when compared with the commonly used algorithm of LS and EKF, the RAEKF algorithm adopts an equivalent weight model and adaptive factors, which can effectively reduce the impact of random error, systematic error and gross error, and finally obtain high-precision navigation results.