A Strap-Down Inertial Navigation/Spectrum Red-Shift/Star Sensor (SINS/SRS/SS) Autonomous Integrated System for Spacecraft Navigation

Abstract

This paper presents a new Strap-down Inertial Navigation System/Spectrum Red-Shift/Star Sensor (SINS/SRS/SS) system integration methodology to improve the autonomy and reliability of spacecraft navigation using the spectrum red-shift information from natural celestial bodies such as the Sun, Jupiter and the Earth. The system models for SINS/SRS/SS integration are established. The information fusion of SINS/SRS/SS integration is designed as the structure of the federated Kalman filter to fuse the local estimations of SINS/SRS and SINS/SS integrated subsystems to generate the global state estimation for spacecraft navigation. A new robust adaptive unscented particle filter is also developed to obtain the local state estimations of SINS/SRS and SINS/SS integrated subsystems in a parallel manner. The simulation results demonstrate that the proposed methodology for SINS/SRS/SS integration can effectively calculate navigation solutions, leading to strong autonomy and high reliability for spacecraft navigation.

1. Introduction

Considerable research efforts have been dedicated to spacecraft navigation, resulting in various navigation techniques such as ground radio navigation, satellite navigation system, inertial navigation system (INS) and celestial navigation system [1,2]. Radio navigation is a non-autonomous navigation technique [3,4]. It is sensitive to external disturbances, since its navigation accuracy depends on the coverage area of ground stations and the propagation conditions of radio waves. The satellite navigation is an extension of the space-based radio navigation [5,6]. It is simple in implementation and highly accurate in positioning. However, its performance is vulnerable to human-induced disturbance. Consequently, the satellite navigation cannot achieve full autonomy. The Strap-down Inertial Navigation System (SINS) has a simple structure and strong autonomy, and is commonly used in vehicle navigation. Nevertheless, the SINS error accumulates over time, leading to biased or even divergent navigation solutions [7,8,9]. The celestial navigation, which uses the star sensor (SS) to derive spacecraft position and attitude, is also an autonomous navigation system. It can provide high-precision orientation information and inhibit electromagnetic interference [6,10]. However, it cannot directly measure spacecraft velocity and has a low data update rate. It also suffers from the limitation of SS, and thus cannot be used alone for positioning and navigation.

The spectrum red-shift (SRS) is a relatively new technology for autonomous navigation. Different from SS, SRS calculates the relative velocity between the spacecraft and the high-precision celestial bodies such as the Sun, Jupiter and the Earth. The velocity of the spacecraft can be obtained from the spectral information (celestial ephemeris) of the solar system, without requiring any information on ground radio and relying on spacecraft orbital dynamics equations. This method has the merits of simple implementation, high precision, strong autonomy and excellent real-time performance [11,12,13], leading to a promising solution to improve the autonomy of spacecraft navigation. However, spectrum signals can be temporarily interrupted by the occlusion of celestial bodies, and thus the navigation solution of SRS may be deteriorated or even divergent due to insufficient measurement information [14,15,16,17].

Given the complementary nature of SINS, SRS and SS, it is absolutely necessary to develop an integrated navigation system by integrating these sensors together to overcome their respective shortages, leading to an improved performance for spacecraft navigation. However, there has been very limited research on SINS/SRS/SS integrated navigation systems. Just recently, Wei et al. studied a Strapdown Inertial Navigation System/Spectral Red-Shift/Geomagnetic Navigation System (SINS/SRS/GNS) integrated navigation system [18]. However, the attitude accuracy of GNS is limited, being much lower than that of SS. This limited attitude accuracy of GNS also causes the attitude accuracy of SINS/SRS/GNS integration to be limited. GNS also involves an expensive computational process to search and match local geomagnetic maps and calculate local coordinates.

The essential concept of the SINS/SRS/SS integrated navigation system is multi-sensor data fusion. The federated Kalman filter (FKF) is a popular multi-sensor data fusion strategy in integrated navigation systems [1,19,20]. It conducts local and global filtering based on the principle of information sharing and also discards the dependence of local estimations via upper bounds [1,19,20].

However, the performance of FKF is heavily dependent on that of local filtering. At present, the extended Kalman filter (EKF), unscented Kalman filter (UKF) and particle filter (PF) are the typical filtering algorithms used for nonlinear systems. EKF is a suboptimal algorithm for nonlinear state estimation, where the system model is linearized by a Taylor expansion [21]. As the linearization process causes a significant error, the EKF solution may be biased or divergent [22]. EKF also requires the calculation of Jacobian matrix, which is difficult to achieve when measurement is strong nonlinear and systems noise is non-Gaussian. UKF is also a nonlinear filtering algorithm by approximating the probability density of state distribution based on unscented transform [21]. It reduces the linearization error of EKF and does not involve the cumbersome calculation of Jacobian matrix. However, it causes an extra computational load and is not stable in case of high-dimensional non-Gaussian systems. PF provides optimal Bayesian approximations of posterior distributions by Monte-Carlo simulation. It is suitable for strong nonlinear and non-Gaussian systems. However, it requires appropriate importance sampling, which is difficult to determine and thus may lead to degraded or divergent solutions. The central difference particle filter (CDPF) adopts the central difference Kalman filter to improve the importance sampling, leading to second-order accuracy for the mean and variance [18]. The Unscented Particle Filter (UPF) improves CDPF by implementing sampling points using unscented transformation to approximate the posterior density function for nonlinear systems. It can achieve the mean and variance in third-order accuracy. However, like PF, UPF still suffers from particle degradation [23].

This paper presents a new methodology for a SINS/SRS/SS autonomous integrated navigation system to improve the autonomy, reliability and accuracy of spacecraft navigation. The system models for SINS/SRS/SS integration are established based on the position and velocity of SRS and the attitude of SS. The information fusion for SINS/SRS/SS integration is designed as the structure of FKF, where the local state estimations of SIN/SRS and SINS/SS integrated subsystems are obtained independently by a robust adaptive unscented particle filter (RAUPF) and are further fused to generate the globally optimal state estimation for spacecraft navigation. This RAUPF adopts the concept of robust adaptive filtering in UPF to prevent particles from degeneracy. It uses the equivalent weight function and adaptive factor to improve the importance sampling resulted from unscented transformation based on the information of system state and measurement models. Simulation trails were performed to examine the efficacy of the presented methodology for SINS/SRS/SS integrated navigation system.

This paper is different from Wei’s work on SINS/SRS/GNS integration [18]. It focuses on SINS/SRS/SS integration, which is advantageous to SINS/SRS/GNS integration. Accordingly, the system models for integrated navigation in this paper are different from those in Wei’s work. Further, the filtering algorithm in this paper is also different from that in Wei’s work. This paper develops a RAUPF, while Wei’s work a robust adaptive CDPF (RACDPF), for local fusion. The RACDPF in Wei’s work reduces the computational load, which is caused by GNS for matching local geomagnetic maps and calculating local coordinates. However, since UPF has higher accuracy than CDPF, the RAUPF developed in this paper also has higher accuracy than RACDPF.

2. SRS Navigation

The spacecraft spectrum red-shift autonomous navigation is based on the characteristics of the deep space environment. It treats the optical signals of celestial bodies in the solar system as navigation information sources. According to the Doppler effect, the frequency of the spectrum emitted from a celestial body is not equal to that received by the spacecraft, which varies with the spacecraft motion with reference to the celestial body. Assume there are three celestial bodies, i.e., three light sources, in the solar system. As shown in Figure 1, the velocity of the spacecraft in the inertial coordinate system can be obtained according to the measurement of the spectrum red-shift frequency based on the vector space theory. Accordingly, the position of the spacecraft in inertial coordinate system can be obtained by integrating the velocity over time.

When the spacecraft is moving relative to a celestial body, the relationship between the light wave frequencies $f_m$ and $f_0$ from the celestial body to the spacecraft and ground station can be expressed as:

$$f_m=f_0\frac{\sqrt{1-{\left|\mathit{v}\right|}^2/c^2}}{1+\left|\mathit{v}\right|\cos \theta /c}$$

(1)

where $f_m$ and $f_0$ are the light wave frequencies from the celestial body to the spacecraft and ground station, respectively, $\mathit{v}$ is the two-dimensional velocity vector of the spacecraft relative to the celestial body, $\theta$ is the angle between velocity vector $\mathit{v}$ and the light wave vector from the celestial body to the spacecraft; $\left|\mathit{v}\right|\cos \theta$ is the radial velocity, and $c$ is the velocity of the celestial body in a vacuum.

According to Equation (1), the following relationship may be written:

$$1+\left|\mathit{v}\right|\cos \theta /c=\frac{f_0}{f_m}\sqrt{1-{\left|\mathit{v}\right|}^2/c^2}$$

(2)

Equation (2) may be further written as:

$$\left|\mathit{v}\right|\cos \theta =\frac{f_0}{f_m}\sqrt{c^2-{\left|\mathit{v}\right|}^2}-c$$

(3)

Equation (3) is actually the expression of the radial velocity:

$$\begin{array}{l}{\mathit{v}}_r=\left|\mathit{v}\right|\cos \theta \\ =\frac{f_0}{f_m}\sqrt{c^2-{\left|\mathit{v}\right|}^2}-c\end{array}$$

(4)

Thus, for the case of three celestial bodies as shown in Figure 1, the light wave frequencies $f_{m1}$, $f_{m2}$ and $f_{m3}$ received by the spacecraft can be obtained from the light wave frequencies $f_{01}$, $f_{02}$ and $f_{03}$ received and measured by the ground station.

According to Equation (4), the radial velocity of the spacecraft relative to the three celestial bodies can be expressed as:

$$\{\begin{array}{l}{v}_{r1}=\frac{f_{01}}{f_{m1}}\sqrt{c^2-{\left|{\mathit{v}}_p\right|}^2}-c\\ v_{r2}=\frac{f_{02}}{f_{m2}}\sqrt{c^2-{\left|{\mathit{v}}_p\right|}^2}-c\\ v_{r3}=\frac{f_{03}}{f_{m3}}\sqrt{c^2-{\left|{\mathit{v}}_p\right|}^2}-c\end{array}$$

(5)

where ${\mathit{v}}_p$ is the velocity vector of the spacecraft in the inertial coordinate system. By spatial geometry reasoning in relation to the three celestial bodies, the relationship between ${\mathit{v}}_p$ and $v_{r1}$, $v_{r2}$, $v_{r3}$ can be expressed as:

$$\{\begin{array}{l}{v}_{r1}=\left({\mathit{v}}_P-{\mathit{v}}_1\right)\cdot {\mathit{u}}_1\\ v_{r2}=\left({\mathit{v}}_P-{\mathit{v}}_2\right)\cdot {\mathit{u}}_2\\ v_{r3}=\left({\mathit{v}}_P-{\mathit{v}}_3\right)\cdot {\mathit{u}}_3\end{array}$$

(6)

where ${\mathit{v}}_1$, ${\mathit{v}}_2$ and ${\mathit{v}}_3$ are three celestial bodies’ velocities in the inertial frame. ${\mathit{u}}_1$, ${\mathit{u}}_2$ and ${\mathit{u}}_3$ are the unit vectors in the inertial frame, pointing to the spacecraft from each of the three celestial bodies.

The state equations in terms of the velocity vector and position vector can be described as:

$$\{\begin{array}{l}{\mathit{v}}_P=f\left({\mathit{v}}_1,{\mathit{v}}_2,{\mathit{v}}_3,u_1,u_2,u_3,v_{r1},v_{r2},v_{r3}\right)\\ {\mathit{p}}_P={\displaystyle \int {\mathit{v}}_{P}dt}\end{array}$$

(7)

where $f(\cdot )$ describes the nonlinear function. Given the initial value, the velocity vector ${\mathit{v}}_p$ of the spacecraft in the inertial frame may be acquired via Equation (7). Then, the position vector ${\mathit{p}}_P$ can be obtained by integration.

3. System Models for SINS/SRS/SS Autonomous Integrated Navigation

The SINS/SRS/SS autonomous integrated navigation system consists of integrated SINS/SRS and SINS/SS subsystems. Its navigation coordinate system is chosen as East-North-Up (E-N-U) geographic coordinate system.

3.1. System State Equation

The system state vector is defined as:

$$\mathit{X}\left(t\right)={\left[\begin{array}{ccccccccccccccc}\delta v_E& \delta v_N& \delta v_U& \delta \phi & \delta \lambda & \delta h& {\psi }_E& {\psi }_N& {\psi }_U& {\theta }_x& {\theta }_y& {\theta }_z& {\zeta }_x& {\zeta }_y& {\zeta }_z\end{array}\right]}^{\mathrm{T}}$$

(8)

where $\left(\delta v_E,\delta v_N,\delta v_U\right)$, $(\delta \phi ,\delta \lambda ,\delta h)$ and $\left({\psi }_E,{\psi }_N,{\psi }_U\right)$ are the errors in velocity, position and attitude; $\left({\theta }_x,{\theta }_y,{\theta }_z\right)$ is the constant bias of the gyros; and $\left({\zeta }_x,{\zeta }_y,{\zeta }_z\right)$ the zero bias of the accelerometers.

The system state equation of SINS/SRS/SS integration is described by:

$$\dot{\mathit{X}}\left(t\right)=f\left(\mathit{X}\left(t\right)\right)+\mathit{G}(t)\mathit{W}(t)$$

(9)

where $f(\cdot )$ is a nonlinear function given by Equation (10) [19], $\mathit{X}\left(t\right)$ is the system state vector, and $\mathit{W}\left(t\right)$ is the system noise.

$$\begin{array}{l}f\left(\mathit{X}\left(t\right)\right)=\\ \hspace{1em}\hspace{1em}\left[\begin{array}{l}{\mathit{C}}_{\omega }^{-\mathit{1}}\left[\left(\mathit{I}-{\mathit{C}}_n^c\right){\widehat{\omega }}_{in}^n+{\mathit{C}}_n^c\delta {\omega }_{in}^n-{\mathit{C}}_b^c\delta {\omega }_{ib}^b\right]\\ \left[\mathit{I}-{\left({\mathit{C}}_n^c\right)}^{\mathrm{T}}\right]{\mathit{C}}_b^c{\widehat{\mathit{f}}}_{sf}^b+{\left({\mathit{C}}_n^c\right)}^{\mathrm{T}}{\mathit{C}}_b^c\delta {\mathit{f}}_{sf}^b\\ -\left(2\delta {\omega }_{ie}^n+\delta {\omega }_{en}^n\right)\times \mathit{V}-\left(2{\widehat{\omega }}_{ie}^n+{\widehat{\omega }}_{en}^n\right)\times \delta \mathit{V}\\ +\left(2{\omega }_{ie}^n+{\omega }_{en}^n\right)\times \delta \mathit{V}+\delta g\\ \frac{v_N}{R_M+h}-\frac{\left(v_N-\delta v_N\right)}{\left(R_M-\delta R_M\right)+\left(h-\delta h\right)}\\ \frac{v_E\sec \phi }{R_N+h}-\frac{\left(v_E-\delta v_E\right)\sec \left(\phi -\delta \phi \right)}{\left(R_N-\delta R_N\right)+\left(h-\delta h\right)}\\ \delta v_U\\ {\mathbf{0}}_{1\times 7}\end{array}\right]\end{array}$$

(10)

where ${\mathit{C}}_{\omega }^{}$ is the Euler platform error angel matrix; ${\mathit{C}}_n^c$ is the transformation matrix from the navigation ($n$) to computer ($c$) frames in terms of attitude; ${\mathit{C}}_b^c$ is the transformation matrix from the body ($b$) to computer ($c$) frames in terms of attitude; ${\widehat{\omega }}_{ie}^n$ and ${\omega }_{ie}^n$ are the projections of the actual and ideal values of angular velocity into ($n$) from the Earth coordinate system ($e$) to the inertial coordinate system ($i$); $\delta {\omega }_{in}^n$ and $\delta {\omega }_{ib}^b$ are the calculation errors of ${\omega }_{in}^n$ and ${\omega }_{ib}^b$; ${\widehat{\mathit{f}}}_{sf}^b$ and $\delta {\mathit{f}}_{sf}^b$ are the real accelerometer’s force and its associated error; $\mathit{V}$ and $\delta \mathit{V}$ are the real velocity and its associated error; $\delta \mathit{g}$ is the error of gravity acceleration; $\phi$ and $h$ are the latitude and altitude values; and $R_M$ and $R_N$ are the radii of curvatures of local meridian and prime vertical.

The noise coefficient matrix is defined as:

$$\mathit{G}\left(t\right)={\left[\begin{array}{cc}\begin{array}{cc}{\mathit{C}}_{b(3\times 3)}^c& {\mathbf{0}}_{(3\times 3)}\\ {\mathbf{0}}_{(3\times 3)}& {\mathit{C}}_{b(3\times 3)}^c\end{array}& {\mathbf{0}}_{\left(9\times 7\right)}\\ {\mathbf{0}}_{(10\times 6)}& {\mathit{I}}_{(7\times 7)}\end{array}\right]}_{\left(16\times 13\right)}$$

(11)

3.2. Measurement Equation of SINS/SRS Integrated Subsystem

The architecture of SINS/SRS integrated subsystem is shown in Figure 2. In order to overcome the drawback that the velocity error of SINS is accumulated in time series, the velocity information from SRS is used to correct the velocity error of SINS. Further, a radar altimeter is used to provide the altitude information to correct the altitude channel of SINS.

The measurement of SINS/SRS integrated subsystem can be chosen as the difference between the velocities of SINS and SRS as well as the difference of altitude between the radar altimeter and SINS.

Suppose the spacecraft velocities obtained by SRS and SINS are ${\mathit{V}}_{SRS}=\left(v_{SE},v_{SN},v_{SU}\right)$ and ${\mathit{V}}_{SINS}=\left(v_E,v_N,v_U\right)$. The difference of spacecraft velocity between SRS and SINS is defined as:

$$\Delta {\mathit{Z}}_v=\left[\begin{array}{c}\Delta Z_{S1}\\ \Delta Z_{S2}\\ \Delta Z_{S3}\end{array}\right]=\left[\begin{array}{c}{v}_E-v_{SE}\\ v_N-v_{SN}\\ v_U-v_{SU}\end{array}\right]={\mathit{H}}_v\mathit{X}\left(t\right)+{\mathit{V}}_v\left(t\right)$$

(12)

where ${\mathit{V}}_v$ is the velocity measurement noise matrix, and ${\mathit{H}}_v={\left[\begin{array}{cc}{\mathit{I}}_{3\times 3}& {\mathbf{0}}_{3\times 12}\end{array}\right]}^{\mathrm{T}}$ is the measurement matrix.

Further, the difference of spacecraft altitude between the radar altimeter and SINS can be described as:

$$\Delta {\mathit{Z}}_h=\left[h_{SINS}-h_H\right]={\mathit{H}}_h\mathit{X}\left(t\right)+{\mathit{V}}_h\left(t\right)$$

(13)

where $h_{SINS}$ and $h_H$ are the spacecraft altitudes obtained by SINS and the radar altimeter; ${\mathit{V}}_h\left(t\right)$ is the measurement noise matrix of spacecraft altitude; and ${\mathit{H}}_h={\left[\begin{array}{ccc}{\mathbf{0}}_{1\times 5}& 1& {\mathbf{0}}_{1\times 9}\end{array}\right]}^{\mathrm{T}}$ is the measurement matrix.

Combining Equation (12) with Equation (13), the measurement equation of SINS/SRS integrated subsystem is established as:

$$\Delta {\mathit{Z}}_1\left(t\right)=\left[\begin{array}{c}{\mathit{H}}_v\\ {\mathit{H}}_h\end{array}\right]\mathit{X}\left(t\right)+\left[\begin{array}{c}{\mathit{V}}_v\left(t\right)\\ {\mathit{V}}_h\left(t\right)\end{array}\right]={\mathit{H}}_1\left(t\right)\mathit{X}\left(t\right)+{\mathit{V}}_1\left(t\right)$$

(14)

3.3. Measurement Equation of SINS/SS Integrated Subsystem

Figure 3 shows the framework of SINS/SS integrated subsystem, where the spacecraft attitude obtained from SS is applied to correct the attitude error of SINS. The difference of spacecraft attitude between SS and SINS is taken as the measurement of SINS/SS integrated subsystem. The measurement equation is given by:

$$\Delta {\mathit{Z}}_2\left(t\right)={\mathit{H}}_2\mathit{X}\left(t\right)+{\mathit{V}}_2\left(t\right)$$

(15)

where ${\mathit{H}}_2=\left[\begin{array}{ccc}{\mathbf{0}}_{{}_{3\times 6}}& {\mathit{I}}_{3\times 3}& {\mathbf{0}}_{{}_{3\times 6}}\end{array}\right]$, and ${\mathit{V}}_2(t)$ is the measurement noise corresponding to the measurement error of SS.

4. Information Fusion for SINS/SRS/SS Integrated Navigation

In this section, a fusion framework is designed as the structure of FKF for INS/SRS/SS integration, where the local state estimations of SINS/SRS and SINS/SS integrated subsystems are obtained by a local filter in a parallel manner and further fused to generate the global state estimation for spacecraft navigation. The local filter is implemented by the proposed RAUPF to obtain the local state estimations of SINS/SRS and SINS/SS integrated subsystems.

4.1. Fusion Framework

As a multi-sensor system, information fusion is a key element in SINS/SRS/SS integration to achieve optimal state estimation. In this paper, the fusion framework of SINS/SRS/SS integration is designed as the structure of FKF. As shown in Figure 4, the two local filters are implemented by the proposed RAUPF to calculate the locally optimal estimations ${\widehat{\mathit{X}}}_i$ and the corresponding error covariance matrices ${\Sigma }_{\widehat{X},i}(i=1,2)$ of SINS/SRS and SINS/SS integrated subsystems in a parallel manner. Subsequently, the two local optimal estimations are fused by using the federated filtering technology [6,24,25,26] to obtain the globally optimal state estimation ${\widehat{\mathit{X}}}_g$ and the associated error covariance matrix ${\Sigma }_{\widehat{X},g}$, i.e.,

$$\{\begin{array}{c}{\Sigma }_{\widehat{X},g}={({\Sigma }_{\widehat{X},1}^{-1}+{\Sigma }_{\widehat{X},2}^{-1})}^{-1}\\ {\widehat{\mathit{X}}}_g={\Sigma }_{\widehat{X},g}({\Sigma }_{\widehat{X},1}^{-1}{\widehat{\mathit{X}}}_1+{\Sigma }_{\widehat{X},2}^{-1}{\widehat{\mathit{X}}}_2)\end{array}$$

(16)

where the local state estimations of SINS/SRS and SINS/SS integrated subsystems are generated by the two local filters in a parallel manner and are subsequently fused to generate the global state estimation for INS/SRS/SS integrated navigation.

Finally, the SINS error is rectified by ${\widehat{\mathit{X}}}_g$ in real time. As mentioned previously, since the altitude channel of SINS is divergent and neither SRS nor SS can output the spacecraft altitude, the radar altimeter is used to correct the SINS altitude channel to suppress the divergence in the altitude of SINS.

4.2. Robust Adaptive Unscented Particle Filter

Consider the following nonlinear system [27,28]:

$$\{\begin{array}{l}{x}_k=f(x_{k-1},w_{k-1})\\ z_k=h(X_k,V_k)\end{array}$$

(17)

where $x_k\in R^n$ is the state vector of the system at time k, $z_k\in R^n$ is the measurement vector, $w_k\in R^n$ is the system noise with variance $R_k$, ${\mathbf{V}}_k\in R^n$ is the measurement noise with variance $Q_k$, both $f(\cdot )$ and $h(\cdot )$ are a nonlinear function, and the sampling time is $k=0,1,\cdots ,N$. The procedure of the proposed RAUPF includes the following steps:

(i)

Initialization. Draw $N$ particles according to the initial state and corresponding error covariance. For time $k=0$, let $x_{{}_0}^i~p(x_o)$, $i=1,2,\cdots ,N$, and set the initial weight as $w_{{}_0}^i=1/N$.

(ii)

For time $k=1,2,\cdots ,N$, conduct the following steps:

(a)

Calculate the equivalent weight matrix $\overline{P}$ and the adaptive factor $\alpha$. In order to enhance the robustness of measurement, the IGG algorithm [29,30] is adopted to construct equivalent weight matrix $\overline{P}$ as the following decreasing function:

$$\begin{gathered} \overline{P}=diag({\overline{P}}_1,{\overline{P}}_2,\cdots ,{\overline{P}}_k) \\ {\overline{P}}_k=(\begin{array}{ccc}{p}_k& & \left|{\epsilon }_k\right|\le k_0\\ p_k\frac{k_0}{\left|{\epsilon }_k\right|}& & k_0<\left|{\epsilon }_k\right|\le k_1\\ 0& & \left|{\epsilon }_k\right|>k_1\end{array} \end{gathered}$$

(18)

In some cases, the equivalent weight matrix can be also defined as [31]:

$${\overline{P}}_k=(\begin{array}{ccc}{p}_k& & \left|{\epsilon }_k\right|\le k_0\\ p_k\frac{k_0}{\left|{\epsilon }_k\right|}\frac{{(k_1-\left|{\epsilon }_k\right|)}^2}{{(k_1-k_0)}^2}& & k_0\le \left|{\epsilon }_k\right|<k_1\\ 0& & k_1\le {\left|{\epsilon }_k\right|}_{}\end{array}$$

(19)

where $k_0\in (1,1.5)$, $k_1\in (3,8)$, ${\epsilon }_k$ is the residual vector of measurement $z_k$, and ${\widehat{x}}_k$ is the current state estimate.

The adaptive factor is selected as:

$$\alpha =(\begin{array}{ccc}1& & \left|\Delta {\tilde{x}}_k\right|\le c_0\\ \frac{c_0}{\left|\Delta {\tilde{x}}_k\right|}\frac{{(c_1-\left|\Delta {\tilde{x}}_k\right|)}^2}{{(c_1-c_0)}^2}& & c_0\le \left|\Delta {\tilde{x}}_k\right|<c_1\\ 0& & c_1\le {\left|\Delta {\tilde{x}}_k\right|}_{}\end{array}$$

(20)

where $c_0\in (1,1.5)$, $c_1\in (3,8)$, $\Delta {\tilde{x}}_k=\frac{\Vert {\widehat{x}}_k-{\overline{x}}_k\Vert }{\sqrt{tr({\displaystyle \sum {\overline{x}}_k})}}$, $tr(\cdot )$ is the trace of the matrix, ${\overline{x}}_k$ is the state prediction, and ${\widehat{x}}_k={(A_k^T{\overline{P}}_k{A}_k+\alpha P_{{\overline{x}}_k})}^{-1}(A_k^T{\overline{P}}_k{l}_k+\alpha P_{{\overline{x}}_k}{\overline{x}}_k)$.

(b)

Calculate the Sigma points. Update the particles $\{x_{k-1}^i,p_{k-1}^i\}$ using UKF to obtain $\{{\overline{x}}_k^i,p_k^i\}$, where ${\overline{x}}_k^i$ satisfies $q(x_k^i|x_{k-1}^i,z_k)=N({\overline{x}}_k^i,p_k^i)$. Taking ${\overline{x}}_k^i$ as a new sample, $2N+1$ Sigma points are selected as [27,28]:

$${\chi }_{k-1}^i=[{\overline{x}}_{k-1}^i,{\overline{x}}_{k-1}^i+\sqrt{(N+\lambda )P_{k-1}^i},{\overline{x}}_{k-1}^i-\sqrt{(N+\lambda )P_{k-1}^i}]$$

(21)

where $\lambda ={\epsilon }^2(n+\rho )$ represents the second order scale factor, n is the dimension of the system state, $\rho$ is the adjustment coefficient, $N$ denotes the quantity of particle samples, and $\epsilon$ is the factor to measure the sample distribution in regard to the mean of the predicted state. Subsequently, UKF is used to update the particles.

(iii)

Calculate the weights $w_k^i=w_{k-1}^i\frac{p(y_k|{\overline{x}}_k^{i^{\ast }})p({\overline{x}}_k^{i^{\ast }}|{\overline{x}}_{k-1}^{i^{\ast }})}{q({\overline{x}}_k^{i^{\ast }}|{\overline{x}}_{k-1}^{i^{\ast }},z_k)}$, and normalize them as ${\tilde{w}}_k^i=w_k^i/{\displaystyle \sum _{i=1}^n{w}_k^i}$. In the above formula, ${\overline{x}}_k^i{}^{\ast }={\overline{x}}_{k|k-1}^i,+K_k^{\ast }(z_k^{}-{\overline{z}}_{k|k-1}^i)$, $K_k^{\ast }=P_{x_k{l}_k}{\overline{P}}_{l_k{l}_k}^{-1}$, and ${\overline{z}}_{k|k-1}^i={\displaystyle \sum _{j=0}^{2N}{W}_j^m}{z}_{j,k|k-1}^i$.

(iv)

Calculate the normalized estimate ${\widehat{N}}_{eff}=1/{\displaystyle \sum _{i=1}^N{({\tilde{w}}_{{}_k}^i)}^2}$. The particle degradation can be measured with the value of ${\widehat{N}}_{eff}$. If ${\widehat{N}}_{eff}$ is smaller, the particle degradation will be worse. In case particles are severely degraded, once can resample the posterior density and further allocate the same factor $\frac{1}{M}$ to every newly sampled particle.

(v)

Calculate the state estimation ${\widehat{x}}_k^{\ast }={\displaystyle \sum _{i=1}^N{\tilde{w}}_k^i}{\overline{x}}^i{{}_k}^{\ast }$. Then, repeat the above steps (ii)–(iv) for the next time point.

In the above steps, in order to achieve a more effective distribution function for the process of importance sampling, the particles resulted from unscented transformation are governed by equivalent weight matrix $\overline{P}$, $N$ and $\alpha$.

5. Simulations and Analysis

The presented methodology was assessed by conducting simulations for SINS/SRS/SS integration in terms of the flight of a spacecraft. Comparison analysis with SINS as well as SINS/SRS and SINS/SS integrated subsystems were also conducted to demonstrate the efficiency of the proposed methodology for SINS/SRS/SS integration. Further, simulations and comparison analysis with EKF, UKF, PF and UPF were also conducted to evaluate the RAUPF effectiveness for SINS/SRS/SS integration.

The navigation coordinate system is the East-North-Up geocentric coordinate system. Assume the spacecraft orbits the Earth, and the orbit parameters are described in

Table 1. The parameters of the spacecraft.

Orbit Parameters	Value
Semi-major axis	7128.042335 km
Eccentricity	0.001169
Orbit inclination	29.268°
Right ascension of ascending node	341.24°
Argument of perigee	347.128°
True anomaly	232.57°
Orbit cycle	49755 s

. A flight period of 1000 s was selected for the simulation test, where the initial position was (2,207,542 m, 3,393,318 m, −2,194,259 m), and the end position was (1,229,267 m, 3,515,079 m, −1,194,458 m). The flight trajectory is shown in Figure 5.

In the process of simulation, the SINS measurements, namely, the angular velocity increments of the gyros and the specific force of the accelerometers, are generated according to the flight trajectory of the spacecraft and the parameters of Earth rotation. The SRS measurements are generated by ASTM G173-03 Reference Spectra derived from SM-ARTS [32]. The SS measurements are obtained from the astronomical ephemeris information and the flight trajectory of the spacecraft.

The initial alignment error of SINS is zero. The initial position error, initial velocity error and initial attitude error of the spacecraft are (10 m, 10 m, 10 m), (1 m/s, 1 m/s, 1 m/s) and ($1^{\prime }$, $1^{\prime }$, $1^{\prime }$). The unscented transformation parameters are $\alpha =0.5$ and $\beta =2$. The adaptive factor calculation parameters are $c_0=1$ and $c_1=3.5$. The equivalent weight matrix calculation parameters are $k_0=1$ and $k_1=4.2$. The number of particles is $M=200$. The sensor parameters used in the simulation test are shown in

Table 2. The sensor parameters in the simulation test.

Sensor Parameters	Values
Gyro constant bias	$0.02°/\mathrm{h}$
Gyro random drift	$0.005°/\sqrt{\mathrm{h}}$
Accelerometer constant bias	0.05 mg
Accelerometer random drift	$0.005\mathrm{mg}/\sqrt{\mathrm{h}}$
Radar Altimeter measurement accuracy	5 m
SS measurement accuracy	$20^{″}$
Spectral redshift estimation accuracy	Standard deviation ≤ $1.8\times {10}^{-9}$

.

5.1. Performances of SINS/SRS and SINS/SS Integrated Subsystems

Simulation trials were conducted to evaluate the navigation performances of the SINS/SRS and SINS/SS integrated subsystems. The solutions of the SINS/SRS and SINS/SS integrated subsystems were achieved under the same conditions using RAUPF, and were further compared with the flight trajectory as reference to calculate their individual navigation errors. For comparison analysis, the SINS navigation error was also calculated in the simulation trials. Figure 6, Figure 7 and Figure 8 show the course angle errors, East velocity errors and latitude errors of SINS, and SINS/SRS and SINS/SS integrated subsystems.

Table 3. Errors of SINS/SRS and SINS/SS integrated subsystems.

Error Type		Integrated Subsystem
		SINS/SS	SINS/SRS
Attitude RMSE (′)	Pitch angle	0.886	2.7325
	Roll angle	0.86768	2.1698
	Course angle	1.0858	5.102
Velocity RMSE (m/s)	East	0.4945	0.1832
	North	0.5147	0.1987
	Up	0.5124	0.1963
Position RMSE (m)	Longitude	23.1405	2.7173
	Latitude	33.7497	2.6735
	Altitude	2.3447	2.3895

summaries the root mean square errors (RMSEs) of SINS/SRS and SINS/SS integrated subsystems.

From the above simulation results, it is obvious that the standalone SINS cannot provide a high-precision navigation solution. Its attitude error, velocity error and position error are accumulated and divergent with time. Further, for SINS/SS integrated subsystem, SS can effectively correct the attitude error of SINS with its highly accurate attitude. However, it cannot effectively correct the position error of SINS. SINS/SRS integrated subsystem is the opposite of SINS/SS integrated subsystem, where SRS can effectively correct the velocity error of SINS with its high-accurate velocity, but it cannot effectively correct the attitude error of SINS. Thus, neither SINS/SS nor SINS/SRS integrated system can provide reliable solutions for spacecraft navigation.

5.2. Performance of SINS/SRS/SS integrated navigation system

Simulation trials were further conducted under the same conditions as those in Section 5.1 to evaluate the performance of SINS/SRS/SS integration. The simulation results are shown in Figure 9, Figure 10 and Figure 11.

Table 4. Positioning errors of SINS/SRS/SS integrated navigation system.

Parameter		RMSE
Attitude (′)	Pitch	0.8399
	Roll	0.8431
	Course	0.8470
Position (m)	Longitude	2.6401
	Latitude	2.4153
	Altitude	2.5762
Velocity (m/s)	East	0.1911
	North	0.1962
	Up	0.1951

shows the SINS/SRS/SS integration errors.

It can be seen from Figure 9, Figure 10 and Figure 11 that the attitude, velocity and position errors of SINS/SRS/SS integration is within 1′, 0.6m/s and 5m. This demonstrates that SINS/SRS/SS integration overcomes the disadvantages of single navigation systems by combining the advantages of SINS, SRS and SS, leading to the improved navigation accuracy and reliability.

5.3. Performance of RAUPF

Simulation trials were also conducted under the same conditions to evaluate the performance of the proposed RAUPF in comparison with EKF, UKF, PF and UPF for SINS/SRS/SS integration. For PF, UPF and RAUPF, three different particle numbers ($M=100$, $M=200$ and $M=400$) were used in the simulation trials.

Figure 12 shows the latitude errors of EKF and UKF. It can be seen that the estimation accuracy of EKF is lower than that of UKF. This is because the linearization of the nonlinear model in EKF causes a large error to the state estimation. Although UKF improves the filtering accuracy of EKF, the improved accuracy is still limited. The reason is that UKF approximates the posterior probability distribution of the system state using the Gaussian distribution. When the posterior probability distribution of the system state is non-Gaussian, which is the case of the simulation test, the UKF performance will be significantly degraded. Therefore, both EKF and UKF have limited navigation accuracy for the spacecraft navigation.

Figure 13 shows the latitude errors of PF, UPF and RAUPF, where the particle number is $M=200$. Comparing Figure 13 with Figure 12, it is evident that all three particle filters (PF, UPF and RAUPF) have higher accuracy than both EKF and UKF. This is mainly because these three particle filters use samples to describe the a priori and a posteriori information, thus discarding the constraint that random samples must satisfy a Gaussian distribution. Further, it can also be seen that RAUPF has much higher accuracy than PF and UPF. This is because RAUPF uses the equivalent weight function and adaptive factor to control particle samples based on system state and measurement models to improve the importance sampling, leading to the enhanced accuracy.

Table 5. The simulation results of the filtering algorithms.

Filtering Algorithms	Position RMSEs/m
	$\mathit{M}=100$	$\mathit{M}=200$	$\mathit{M}=400$
EKF	-	14.7298	-
UKF	-	9.0291	-
PF	7.3213	6.4402	5.7090
UPF	5.9115	5.0808	4.3818
RAUPF	3.1862	2.4610	2.0954

summaries the RMSEs of the above filters, where the RMSEs of PF, UPF and RAUPF are subject to three different particle numbers, i.e., $M=100$, $M=200$ and $M=400$. It can be seen that even with the small number of particles ($M=100$), all three particle filters (PF, UPF and RAUPF) still have higher accuracy than both EKF and UKF.

6. Conclusions

This paper presents a new methodology of SINS/SRS/SS integration for spacecraft navigation. It establishes the system models based on the position and velocity of SRS and the attitude of SS for SINS/SRS/SS integration. It also develops an information fusion framework for SINS/SRS/SS integration, where the local state estimations of SINS/SRS and SINS/SS subsystems are obtained by RAUPF and are further fused to generate the globally optimal state estimation for spacecraft navigation. This RAUPF uses the equivalent weight function and adaptive factor to improve the importance sampling, thus leading to the enhanced accuracy of state estimation. The simulation results demonstrate that the proposed methodology for SINS/SRS/SS integration realizes the globally optimal fusion of SINS/SRS and SINS/SS subsystems based on the locally optimal fusion results of each subsystem, thus effectively improving the autonomy, reliability and accuracy for spacecraft navigation.

Future work will focus on improvement of the proposed methodology. By adopting advanced artificial intelligence such as deep learning neural networks, genetic algorithms and pattern recognitions, the proposed methodology will be improved to intelligently characterize errors and uncertainties and further automatically inhibit their disturbances on state estimation in the fusion process.