Stochastic Integration H_∞ Filter for Rapid Transfer Alignment of INS

Abstract

The performance of an inertial navigation system (INS) operated on a moving base greatly depends on the accuracy of rapid transfer alignment (RTA). However, in practice, the coexistence of large initial attitude errors and uncertain observation noise statistics poses a great challenge for the estimation accuracy of misalignment angles. This study aims to develop a novel robust nonlinear filter, namely the stochastic integration H${}_{\infty }$ filter (SIH${}_{\infty }$F) for improving both the accuracy and robustness of RTA. In this new nonlinear H${}_{\infty }$ filter, the stochastic spherical-radial integration rule is incorporated with the framework of the derivative-free H${}_{\infty }$ filter for the first time, and the resulting SIH${}_{\infty }$F simultaneously attenuates the negative effect in estimations caused by significant nonlinearity and large uncertainty. Comparisons between the SIH${}_{\infty }$F and previously well-known methodologies are carried out by means of numerical simulation and a van test. The results demonstrate that the newly-proposed method outperforms the cubature H${}_{\infty }$ filter. Moreover, the SIH${}_{\infty }$F inherits the benefit of the traditional stochastic integration filter, but with more robustness in the presence of uncertainty.

1. Introduction

The inertial navigation system (INS) is an entirely autonomous navigation system, in the sense that it precisely provides the navigation parameters for the carrier without the aid of a signal from an external device. As is well known, the performance of INS operated on a moving base is heavily influenced by the accuracy of transfer alignment, which is an important research interest of modern navigation techniques. In particular, the alignment issue between a main INS (MINS) and a slave INS (SINS) has drawn extensive attention, partially due to its widespread application in diverse carriers. For the MINS/SINS system, the unknown parameters of SINS are commonly initialized by MINS at the very beginning. The inevitable initial misalignments will lead to poor accuracy of the follow-on navigation operation, and thus, the alignment stage of SINS is necessitated. In [1], the rapid transfer alignment (RTA) on the basis of the “attitude plus velocity” matching scheme was originally proposed to estimate the attitude errors between two INSs. Then, the misalignment angles of SINS with respect to the local navigation frame can be indirectly determined. Compared with traditional transfer alignment, the RTA has superiorities in both model observability and convergence rate, and it has been widely applied for both military and civilian purpose [1,2,3].

The central concerns of RTA are the estimation accuracy and convergence rate, both of which are intimately related to the performance of the filtering algorithm. For the RTA with the linearized model and Gaussian noise, the capacity of the Kalman filter (KF) that provides optimal estimation has been authoritatively verified through numerical simulation and field tests [1]. However, the INS sometimes has to be initiated in a complicated environment, and higher demanding standards are put forward to the filtering algorithm. First, for some quick response missions, the initial misalignment angles of INS may be notably large, and significant nonlinearity will be brought into the model description of RTA [2]. That is, the nonlinear character of the filtering algorithm should be deliberately considered. In addition, the precision of observations severely influences the performance of RTA. Nevertheless, in practice, the phenomenon of uncertain observations frequently occurs due to the existence of multiple disturbances, including outliers, the time delay of data transmission, the coupling effect between the lever arm and deformation, etc. [4,5]. The suppression of negative effects induced by these disturbances is of crucial importance to the improvement of the accuracy and stability of RTA. For this purpose, accurate modeling methods have been developed as intuitive solutions for the compensation of these disturbances [6,7,8]. Yet, the complexity of the model will be severely increased, and the un-modeled errors will destroy the veracity of model. Another way is to take these disturbances as uncertain interfering input of the model, and against the negative effect by means of the robust filtering algorithm. On the premise of the bounded energy of noise, one effective method is the H${}_{\infty }$ filter, which intends to minimize the downside on estimation results in the worst case of the un-modeled disturbance and uncertainty. Contrary to KF, the H${}_{\infty }$ filter requires neither statistical noise properties, nor the exact system model and provides robust and accurate estimations, while the uncertain parameters exist in the model description [9].

Many solutions, as exemplified by the Krein space approach and the game theory approach, have been developed for the linear H${}_{\infty }$ problem [9,10,11,12]. Afterward, these solutions are naturally adopted for nonlinear systems with the utilization of the first-order linearization method; to be more precise, the extended H${}_{\infty }$ filter (EH${}_{\infty }$F) [13,14]. However, the EH${}_{\infty }$F inherits the drawbacks of the extended Kalman filter, i.e., the cumbersome computation of Jacobians and the deterioration of estimation performance in the presence of significant nonlinearity. Recent decades have witnessed the development of sample-point filters on the basis of deterministic sampling quadrature methods (as exemplified by the unscented transform, the cubature integration rule, the Gauss–Hermite quadrature rule, etc.), which either approximate the probability distribution function representing state estimates by a set of sample-points or approximate the nonlinear functions using polynomial expansions [15,16,17,18]. The emergence of these quadrature methods facilitates the arising of the derivative-free H${}_{\infty }$ filter (DFH${}_{\infty }$F), such as the unscented H${}_{\infty }$ filter, cubature H${}_{\infty }$ filter (CH${}_{\infty }$F) and sparse-grid quadrature H${}_{\infty }$ filter [19,20,21]. Compared with EH${}_{\infty }$F, the DFH${}_{\infty }$F releases the restriction of computing Jacobians and achieves high accuracy, as the posterior estimations are accurately approximated to certain order. However, the approximations given by deterministic sampling quadrature methods cannot eliminate systematic errors that exist in the solution of DFH${}_{\infty }$F [15]. That is, the estimation accuracy will be degraded while significant nonlinearity exists in the model description. Moreover, the boundedness and convergence of estimation error cannot be guaranteed due to the local validity of these quadrature methods. To overcome these disadvantages, as of late, a novel filter, named the stochastic integration filter (SIF), was proposed to provide asymptotically exact estimation for nonlinear systems [22,23,24]. The core of SIF is the stochastic spherical-radial integration rule (SSRIR), which approximates the probability distribution function of posterior estimation by a set of randomly-chosen sample-points [25]. Contrary to the deterministic sampling quadrature methods, the SSRIR can theoretically eliminate the systematic error with increasing iterations. Therefore, it can be reasonably concluded that the incorporation of SSRIR with the framework of DFH${}_{\infty }$F will result in a filter with high accuracy and strong robustness, but there were few published literature works focused on this study. This paper is devoted to developing a stochastic integration H${}_{\infty }$ filter (SIH${}_{\infty }$F) to address the RTA issue with the coexistence of significant nonlinearity and uncertain observation noise statistics. For the first time, the SSRIR is combined with the framework of DFH${}_{\infty }$F, and the resulting SIH${}_{\infty }$F achieves estimation with both accuracy and robustness. Numerical simulation and the van test are separately executed for the validation of the proposed method. The results demonstrate that the SIH${}_{\infty }$F has better performance than the previous, well-known CH${}_{\infty }$F and SIF.

The rest of this paper is organized as follows: The nonlinear model of RTA is briefly introduced in Section 2. Section 3 firstly represents the preliminaries of the H${}_{\infty }$ technique and the EH${}_{\infty }$F. Then, the derivation of SIH${}_{\infty }$F is described in detail. Numerical simulation is executed in Section 4, in order to demonstrate the superiorities of SIH${}_{\infty }$F compared with the CH${}_{\infty }$F and SIF. In Section 5, the validity of the proposed method is further verified through a van test. The conclusions are given in Section 6.

2. System Modeling of Rapid Transfer Alignment

The aim of system modeling is to describe the error propagation principle of RTA. In this section, the process model and observation model of RTA with the nonlinear characteristic are briefly introduced. The system modeling lays the foundation for the development of the filtering algorithm in the following section.

The definitions of coordinate systems used in system modeling are shown in Figure 1, where i-frame denotes the inertial frame, e-frame denotes the Earth-centered Earth-fixed frame, n-frame denotes the local navigation frame, m-frame denotes the body frame of MINS, $s_r$-frame denotes the body frame of SINS and $s_c$-frame denotes the calculated body frame of SINS. Note that all of these coordinated systems are right handed and Cartesian; the detailed specifications of these coordinate systems can be found in [2].

In this paper, the two INSs are assumed to be directly installed on the carrier, and the sketch map of misalignment angles between the MINS and SINS is represented in Figure 2, where ${\mathit{\psi }}_m$ represents the misalignment angles between the m-frame and the $s_c$-frame and ${\mathit{\psi }}_a$ represents the misalignment angles between the m-frame and the $s_r$-frame. The eventual purpose of RTA is to determine the attitude transformation matrix from $s_r$-frame to n-frame. Since the errors of MINS have been well compensated by the external aid of the device, the attitude parameters given by MINS can be viewed as reliable. In other words, the goal of RTA can be indirectly achieved by estimating the attitude errors between the two INSs, i.e., ${\mathit{\psi }}_a$. However, due to the inherent accumulative errors of sensors, only the inertial data sensed by SINS projected in the $s_c$-frame imply the misalignment information. Therefore, the description of error propagation is the prerequisite of RTA.

Generally, the process model of RTA includes the attitude error equation, velocity error equation, dynamic models of inertial sensors and actual physical misalignment angles. The attitude error equation mainly describes the function relationship between ${\mathit{\psi }}_m$ and ${\mathit{\psi }}_a$. Assuming that large misalignment angles exist in all three axes, the attitude error equation of RTA is given by [2]:

$${\dot{\mathit{\psi }}}_m=\mathbf{\Xi }\left[\left({\mathit{I}}_{3\times 3}-{\mathit{C}}_m^{s_c}{\mathit{C}}_{s_r}^m\right){\mathit{\omega }}_{n{s}_c}^{s_c}+{\mathit{C}}_m^{s_c}{\mathit{C}}_{s_r}^m{\mathit{\epsilon }}^{s_r}\right]+{\mathit{w}}_m$$

(1)

where $\mathbf{\Xi }$ is defined as:

$$\mathbf{\Xi }=\begin{array}{c}\left(\begin{array}{ccc}cos{\psi }_{my}& 0& sin{\psi }_{my}\\ sin{\psi }_{my}tan{\psi }_{mx}& 1& -cos{\psi }_{my}tan{\psi }_{mx}\\ -sin{\psi }_{my}/cos{\psi }_{mx}& 0& cos{\psi }_{my}/cos{\psi }_{mx}\end{array}\right)\end{array}$$

(2)

and ${\psi }_{m,i}$, ($i=x,y,z$) denotes the measurable misalignment angle in the corresponding axis; ${\mathit{I}}_{3\times 3}$ is the unit matrix; ${\mathit{C}}_m^{s_c}$ denotes the attitude transformation matrix from the m-frame to the $s_c$-frame; ${\mathit{C}}_{s_r}^m$ represents the attitude transformation matrix from the $s_r$-frame to the m-frame; ${\mathit{\omega }}_{n{s}_c}^{s_c}$ is the angular velocity of the $s_c$-frame relative to the n-frame projected in the $s_c$-frame; ${\mathit{\epsilon }}^{s_r}$ is the constant drift of the gyroscope with dynamic model ${\dot{\mathit{\epsilon }}}^{s_r}={\mathbf{0}}_{3\times 1}$; ${\mathit{w}}_m$ is the noise term of the attitude error equation.

On the basis of the velocity solution given by the strapdown algorithm, the velocity error equation can be described as [2]:

$$\delta {\dot{\mathit{V}}}^n={\mathit{C}}_{s_c}^n\left({\mathit{I}}_{3\times 3}-{\mathit{C}}_m^{s_c}{\mathit{C}}_{s_r}^m\right){\mathit{f}}_{i{s}_r}^{s_r}-\left(2{\mathit{\omega }}_{ie}^{n\times }+{\mathit{\omega }}_{en}^{n\times }\right)\delta {\mathit{V}}^n+{\mathit{C}}_{s_r}^n{\mathbf{\nabla }}^{s_r}+{\mathit{w}}_v$$

(3)

where $\delta {\mathit{V}}^n$ is the velocity error vector projected in the n-frame; ${\mathit{\omega }}_{ie}^n$ is the angular velocity of the e-frame relative to the i-frame projected in the n-frame; ${\mathit{f}}_{i{s}_r}^{s_r}$ is the specific force sensed by SINS projected in the $s_r$-frame; ${\mathit{\omega }}_{en}^n$ is the angular velocity of the n-frame relative to the e-frame projected in the n-frame; ${\mathit{C}}_m^n$ is the attitude transformation matrix from the m-frame to the n-frame, and it is commonly provided by MINS in real time; ${\mathbf{\nabla }}^{s_r}$ is the constant bias of the accelerometer with dynamic model ${\dot{\mathbf{\nabla }}}^{s_r}={\mathbf{0}}_{3\times 1}$; ${\mathit{w}}_v$ is the noise term of the velocity error equation. Note, the derivation of Equation (3) assumes that the acceleration induced by the lever arm effect has been well compensated.

Due to the impact of vibration and flexure, the modeling of ${\mathit{\psi }}_a$ makes a compromise between the complexity and accuracy, that is [1]:

$${\dot{\mathit{\psi }}}_a={\mathit{\eta }}_a$$

(4)

where ${\mathit{\eta }}_a$ is the noise term with covariance ${\mathit{Q}}_a$.

Choosing the state vector as $\mathit{x}={\left[\begin{array}{ccccc}{\mathit{\psi }}_m& \delta {\mathit{V}}^n& {\mathit{\epsilon }}^{s_r}& {\mathbf{\nabla }}^{s_r}& {\mathit{\psi }}_a\end{array}\right]}^T$, the process model of RTA can be represented as follows:

$$\left\{\begin{array}{cc}{\dot{\mathit{\psi }}}_m\hfill & =\mathbf{\Xi }\left[\left({\mathit{I}}_{3\times 3}-{\mathit{C}}_m^{s_c}{\mathit{C}}_{s_r}^m\right){\mathit{\omega }}_{n{s}_c}^{s_c}+{\mathit{C}}_m^{s_c}{\mathit{C}}_{s_r}^m{\mathit{\epsilon }}^{s_r}\right]+{\mathit{w}}_m\hfill \\ \delta {\dot{\mathit{V}}}^n\hfill & ={\mathit{C}}_{s_c}^n\left({\mathit{I}}_{3\times 3}-{\mathit{C}}_m^{s_c}{\mathit{C}}_{s_r}^m\right){\mathit{f}}_{i{s}_r}^{s_r}-\left(2{\mathit{\omega }}_{ie}^{n\times }+{\mathit{\omega }}_{en}^{n\times }\right)\delta {\mathit{V}}^n+{\mathit{C}}_{s_r}^n{\mathbf{\nabla }}^{s_r}+{\mathit{w}}_v\hfill \\ {\dot{\mathit{\epsilon }}}^{s_r}\hfill & ={\mathbf{0}}_{3\times 1}\hfill \\ {\dot{\mathbf{\nabla }}}^{s_r}\hfill & ={\mathbf{0}}_{3\times 1}\hfill \\ {\dot{\mathit{\psi }}}_a\hfill & ={\mathit{\eta }}_a\hfill \end{array}\right.$$

(5)

The modeling of observation relies on the matching scheme of RTA. In this paper, the measurable misalignment angles and velocity error are chosen as observations. That is:

$$\mathit{z}=\left[\begin{array}{c}\mathbf{\Theta }\left({\mathit{C}}_n^{s_c}{\mathit{C}}_m^n\right)\\ {\mathit{V}}_{se}^n-{\mathit{V}}_{me}^n\end{array}\right]$$

(6)

where $\mathbf{\Theta }(·)$ represents the function that calculates the Euler angle from the attitude transformation matrix; ${\mathit{V}}_{se}^n$ and ${\mathit{V}}_{me}^n$ denote the velocity of SINS and MINS projected in their own body frame, respectively.

The observation model is given by:

$${\mathit{z}}_k={\mathit{H}}_k{\mathit{x}}_k+{\mathit{v}}_k$$

(7)

where ${\mathit{v}}_k$ is the observation noise. The coefficient matrix ${\mathit{H}}_k$ is described as:

$${\mathit{H}}_k=\left[\begin{array}{ccccc}{\mathit{I}}_{3\times 3}& {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}\\ {\mathbf{0}}_{3\times 3}& {\mathit{I}}_{3\times 3}& {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}\end{array}\right]$$

(8)

The system modeling of RTA implies that many factors could affect the precision of observations provided by MINS; for instance, the unreliable signal from the external aid of the device, the random time-delay effect, the coupling effect of the lever arm and deformation, etc. These unpredictable disturbances may induce uncertain statistics of observation noise, which motivates the investigation of the nonlinear filtering method to improve both the accuracy and robustness of RTA.

3. Stochastic Integration H${}_{\mathbf{\infty }}$ Filter

In this section, the derivations of SIH${}_{\infty }$F are described in detail. First, the preliminaries of the H${}_{\infty }$ technique and EH${}_{\infty }$F are briefly represented. Then, the general form of DFH${}_{\infty }$F is given as the prerequisite of derivations. At last, the construction of SSRIR is introduced to evaluate the Gaussian weighted integrals in DFH${}_{\infty }$F, and SIH${}_{\infty }$F can be consequently achieved by combining the SSRIR with the framework of DFH${}_{\infty }$F.

3.1. Preliminaries of the H${}_{\infty }$ Technique

A representative nonlinear discrete-time system can be represented as follows:

$$\left\{\begin{array}{c}{\mathit{x}}_k=\mathit{f}\left({\mathit{x}}_{k-1}\right)+{\mathit{\omega }}_{k-1}\hfill \\ {\mathit{z}}_k=\mathit{h}\left({\mathit{x}}_k\right)+{\mathit{v}}_k\hfill \\ {\mathit{y}}_k={\mathit{L}}_k{\mathit{x}}_k\hfill \end{array}\right.$$

(9)

where ${\mathit{x}}_k\in {\mathbb{R}}^n$ is the state vector; ${\mathit{z}}_k\in {\mathbb{R}}^m$ is the observation vector; ${\mathit{y}}_k$ is the signal to be estimated; $\mathit{f}(·)$ and $\mathit{h}(·)$ are the process model and the observation model, respectively; ${\mathit{L}}_k$ is a known matrix, which will be replaced by the identity matrix for estimating the whole state vector in this study; the process noise ${\mathit{\omega }}_{k-1}$ and observation noise ${\mathit{v}}_k$ are assumed to be the signals with bounded energy, but unknown statistics, i.e., ${\sum }_{k=0}^{\infty }{\mathit{\omega }}_k^T{\mathit{\omega }}_k<\infty$, ${\sum }_{k=0}^{\infty }{\mathit{v}}_k^T{\mathit{v}}_k<\infty$.

Let ${\mathit{e}}_k=\parallel \widehat{{\mathit{y}}_k}-{\mathit{y}}_k\parallel$ be the estimation error; the essential purpose of the H${}_{\infty }$ technique is to estimate ${\mathit{y}}_k$ to minimize ${\mathit{e}}_k$ under the worst case of initial estimation error, process noise and observation noise. For this purpose, the game theory approach gives a celebrated cost function, so as to map the disturbances and uncertainties to the estimation error. That is:

$$J=\frac{{\displaystyle \sum _{j=0}^k}{∥\begin{array}{c}{\widehat{\mathit{y}}}_j-{\mathit{x}}_j\end{array}∥}^2}{{∥\begin{array}{c}{\mathit{x}}_0-{\widehat{\mathit{x}}}_0\end{array}∥}_{{\mathit{P}}_0^{-1}}^2+{\displaystyle \sum _{j=0}^{k-1}}{∥\begin{array}{c}{\mathit{\omega }}_j\end{array}∥}_{{\mathit{Q}}_{\omega ,j}^{-1}}^2+{\displaystyle \sum _{j=0}^k}{∥\begin{array}{c}{\mathit{v}}_j\end{array}∥}_{{\mathit{R}}_{v,j}^{-1}}^2}$$

(10)

where ${\mathit{P}}_0$ denotes the initial covariance of the state vector and ${\mathit{Q}}_{\omega ,j}$ and ${\mathit{R}}_{v,j}$ are the covariance of process noise and observation noise, respectively. The norm notation stands for ${∥\begin{array}{c}\mathit{a}\end{array}∥}_{\mathit{B}}^2={\mathit{a}}^T\mathit{B}\mathit{a}$.

The H${}_{\infty }$ filter sequentially estimates ${\widehat{\mathit{y}}}_j$ to satisfy the following inequality:

$$supJ<{\gamma }^2$$

(11)

where notation $‘sup’$ stands for the supremum; $\gamma$ represents the error attenuation parameter that is empirically set to a constant value or iteratively computed in real time [19,21].

3.2. Extended H${}_{\infty }$ Filter and Derivative-Free H${}_{\infty }$ Filter

As an intuitive extension of the linear H${}_{\infty }$ solution, the EH${}_{\infty }$F has been well studied to achieve suboptimal solutions for a nonlinear system with uncertainties [13,14]. The implementation of EH${}_{\infty }$F is represented as follows:

$${\widehat{\mathit{x}}}_{k|k-1}=\mathit{f}\left({\widehat{\mathit{x}}}_{k-1|k-1}\right)$$

(12)

$${\mathit{P}}_{k|k-1}={\mathit{F}}_{k-1}{\mathit{P}}_{k-1|k-1}{\mathit{F}}_{k-1}^T+{\mathit{Q}}_{\omega ,k-1}$$

(13)

$${\mathit{K}}_{k|k}={\mathit{P}}_{k|k-1}{\mathit{H}}_k^T{({\mathit{H}}_k{\mathit{P}}_{k|k-1}{\mathit{H}}_k^T+{\mathit{R}}_{v,k})}^{-1}$$

(14)

$${\widehat{\mathit{x}}}_{k|k}={\widehat{\mathit{x}}}_{k|k-1}+{\mathit{K}}_{k|k}\left[{\mathit{z}}_k-\mathit{h}\left({\widehat{\mathit{x}}}_{k|k-1}\right)\right]$$

(15)

$${\mathit{R}}_{e,k}^{}=\left[\begin{array}{cc}{\mathit{R}}_{v,k}+{\mathit{H}}_k^{}{\mathit{P}}_{k|k-1}{\mathit{H}}_k^T& {\mathit{H}}_k^{}{\mathit{P}}_{k|k-1}\\ {\mathit{P}}_{k|k-1}{\mathit{H}}_k^T& {\mathit{P}}_{k|k-1}-{\gamma }^2\mathit{I}\end{array}\right]$$

(16)

$${\mathit{P}}_{k|k}={\mathit{P}}_{k|k-1}-\left[\begin{array}{cc}{\mathit{P}}_{k|k-1}{\mathit{H}}_k^T& {\mathit{P}}_{k|k-1}\end{array}\right]{\mathit{R}}_{e,k}^{-1}\left[\begin{array}{c}{\mathit{H}}_k{\mathit{P}}_{k|k-1}\\ {\mathit{P}}_{k|k-1}\end{array}\right]$$

(17)

where ${\widehat{\mathit{x}}}_{k|k-1}$ is the predicted mean with covariance ${\mathit{P}}_{k|k-1}$; ${\widehat{\mathit{x}}}_{k|k}$ is the filtered mean with covariance ${\mathit{P}}_{k|k}$; ${\mathit{F}}_{k-1}$ and ${\mathit{H}}_k$ represent the Jacobian matrices of $\mathit{f}(·)$ evaluated at the filtered mean ${\widehat{\mathit{x}}}_{k-1|k-1}$ and $\mathit{h}(·)$ evaluated at the predictive mean ${\widehat{\mathit{x}}}_{k|k-1}$, respectively; ${\mathit{K}}_{k|k}$ is the gain matrix; ${\mathit{R}}_{e,k}$ is the auxiliary matrix.

Though the first-order Taylor expansion is utilized in EH${}_{\infty }$F to approximate the nonlinear model, the rough linearization method may degrade the estimation accuracy or even induce divergence in the presence of significant nonlinearity. Moreover, the computation of Jacobians is also an inevitable difficulty for the complicated model description. To overcome the limitations of EH${}_{\infty }$F, multiple DFH${}_{\infty }$Fs have been proposed on the basis of deterministic sampling quadrature methods. Without loss of generality, the general form of DFH${}_{\infty }$F is represented as follows.

Analogous to the nonlinear filter with the Gaussian assumption, the prediction step of DFH${}_{\infty }$F is given by:

$${\widehat{\mathit{x}}}_{k|k-1}=\int \mathit{f}\left({\mathit{x}}_{k-1}\right)\mathcal{N}\{{\mathit{x}}_{k-1};{\widehat{\mathit{x}}}_{k-1|k-1},{\mathit{P}}_{k-1|k-1}\}d{\mathit{x}}_{k-1}$$

(18)

$${\mathit{P}}_{k|k-1}=\int \left[\mathit{f}\left({\mathit{x}}_{k-1}\right)-{\widehat{\mathit{x}}}_{k|k-1}\right]{\left[\mathit{f}\left({\mathit{x}}_{k-1}\right)-{\widehat{\mathit{x}}}_{k|k-1}\right]}^T\mathcal{N}\{{\mathit{x}}_{k-1};{\widehat{\mathit{x}}}_{k-1|k-1},{\mathit{P}}_{k-1|k-1}\}d{\mathit{x}}_{k-1}+{\mathit{Q}}_{\omega ,k-1}$$

(19)

where $\mathcal{N}\{{\mathit{x}}_{k-1};{\widehat{\mathit{x}}}_{k-1|k-1},{\mathit{P}}_{k-1|k-1}\}$ denotes the variable ${\mathit{x}}_{k-1}$ subject to the Gaussian probability distribution function with mean ${\widehat{\mathit{x}}}_{k-1|k-1}$ and covariance ${\mathit{P}}_{k-1|k-1}$, and this notation is available for other equations.

Under the assumption that the integrated state vector needs to be estimated, the filtering step of DFH${}_{\infty }$F can be described as:

$${\widehat{\mathit{x}}}_{k|k}={\widehat{\mathit{x}}}_{k|k-1}+{\mathit{K}}_{k|k}({\mathit{z}}_k-{\widehat{\mathit{z}}}_{k|k-1})$$

(20)

$$\begin{array}{cc}\hfill {\mathit{P}}_{k|k}& ={\mathit{P}}_{k|k-1}-\left[\begin{array}{cc}{\mathit{P}}_{xz,k|k-1}& {\mathit{P}}_{k|k-1}\end{array}\right]{\mathit{R}}_{e,k}^{-1}{\left[\begin{array}{cc}{\mathit{P}}_{xz,k|k-1}& {\mathit{P}}_{k|k-1}\end{array}\right]}^T\hfill \end{array}$$

(21)

where ${\mathit{K}}_{k|k}={\mathit{P}}_{xz,k|k-1}{({\mathit{P}}_{zz,k|k-1}+{\mathit{R}}_{v,k})}^{-1}$ is the gain matrix; ${\widehat{\mathit{z}}}_{k|k-1}$ is the predicted observation with covariance ${\mathit{P}}_{zz,k|k-1}$; ${\mathit{P}}_{xz,k|k-1}$ is the cross-covariance between the state vector and observation vector. Note that these three statistics of observation involve the computation of Gaussian weighted integrals, that is:

$${\widehat{\mathit{z}}}_{k|k-1}=\int \mathit{h}\left({\mathit{x}}_k\right)\mathcal{N}\{{\mathit{x}}_k;{\widehat{\mathit{x}}}_{k|k-1},{\mathit{P}}_{k|k-1}\}d{\mathit{x}}_k$$

(22)

$${\mathit{P}}_{zz,k|k-1}=\int \left[\mathit{h}\left({\mathit{x}}_k\right)-{\widehat{\mathit{z}}}_{k|k-1}\right]{\left[\mathit{h}\left({\mathit{x}}_k\right)-{\widehat{\mathit{z}}}_{k|k-1}\right]}^T\mathcal{N}\{{\mathit{x}}_k;{\widehat{\mathit{x}}}_{k|k-1},{\mathit{P}}_{k|k-1}\}d{\mathit{x}}_k+{\mathit{R}}_{v,k}$$

(23)

$${\mathit{P}}_{xz,k|k-1}=\int \left({\mathit{x}}_k-{\widehat{\mathit{x}}}_{k|k-1}\right){\left[\mathit{h}\left({\mathit{x}}_k\right)-{\widehat{\mathit{z}}}_{k|k-1}\right]}^T\mathcal{N}\{{\mathit{x}}_k;{\widehat{\mathit{x}}}_{k|k-1},{\mathit{P}}_{k|k-1}\}d{\mathit{x}}_k$$

(24)

The auxiliary matrix ${\mathit{R}}_{e,k}$ is given by:

$$\begin{array}{c}\hfill {\mathit{R}}_{e,k}=\left[\begin{array}{cc}{\mathit{R}}_{v,k}+{\mathit{P}}_{zz,k|k-1}& {\mathit{P}}_{xz,k|k-1}^T\\ {\mathit{P}}_{xz,k|k-1}& {\mathit{P}}_{k|k-1}-{\gamma }^2\mathit{I}\end{array}\right]\end{array}$$

(25)

In contrast to the EH${}_{\infty }$F, the prediction step of DFH${}_{\infty }$F is explicitly replaced by the general Gaussian approximations, i.e., Equations (18) and (19). Furthermore, on the basis of the statistical linear error propagation method, the statistics related to Jacobian matrix ${\mathit{H}}_k$ in Equations (13), (16) and (17) are substituted by the following approximations [26]:

$${\mathit{P}}_{zz,k|k-1}=E\left[({\mathit{z}}_k-{\widehat{\mathit{z}}}_{k|k-1}){({\mathit{z}}_k-{\widehat{\mathit{z}}}_{k|k-1})}^T\right]\approx {\mathit{H}}_k{\mathit{P}}_{k|k-1}{\mathit{H}}_k^T(k=1,2,\cdots )$$

(26)

$${\mathit{P}}_{xz,k|k-1}=E\left[({\mathit{x}}_k-{\widehat{\mathit{x}}}_{k|k-1}){({\mathit{z}}_k-{\widehat{\mathit{z}}}_{k|k-1})}^T\right]\approx {\mathit{P}}_{k|k-1}{\mathit{H}}_k^T(k=1,2,\cdots )$$

(27)

3.3. Stochastic Spherical-Radial Integration Rule

It is worth noting that the general form of DFH${}_{\infty }$F requires the evaluation of Gaussian weighted integrals with the form $\mathit{I}\left(\mathit{g}\right)=\int \mathit{g}\left(\mathit{x}\right)\mathcal{N}\{\mathit{x};\widehat{\mathit{x}},\mathit{P}\}d\mathit{x}$. In general, these nonlinear integrals cannot be analytically solved, and the numerical approximation method is required. Contrary to the previous deterministic sampling methods, in this paper, the SSRIR is employed for the evaluation of $\mathit{I}\left(\mathit{g}\right)$. Note that two transformations must be implemented to convert the nonlinear integral expression to the standard form.

The first transformation relates to the integration variable, i.e., $\mathit{x}=\widehat{\mathit{x}}+{\mathit{S}}_x\mathit{u}$, where ${\mathit{S}}_x$ is the lower triangular matrix obtained from the Cholesky decomposition of covariance $\mathit{P}$. That is:

$$\begin{array}{cc}\hfill \mathit{I}\left(\mathit{g}\right)& =\int \mathit{g}\left(\mathit{x}\right)\mathcal{N}\{\mathit{x};\widehat{\mathit{x}},\mathit{P}\}d\mathit{x}\hfill \\ & =\int \mathit{g}(\widehat{\mathit{x}}+{\mathit{S}}_x\mathit{u})\mathcal{N}\{\mathit{u};{\mathbf{0}}_{n_x\times 1},{\mathit{I}}_{n_x}\}d\mathit{u}\hfill \end{array}$$

(28)

where $n_x$ is the dimension of $\mathit{x}$.

The second transformation concerns the conversion of the integral variable to the radial-spherical coordinate system. Let $\mathit{u}=r\widehat{\mathit{u}}$, where ${\widehat{\mathit{u}}}^T\widehat{\mathit{u}}=1$. The domain of definition for radius is $r\in [0,\infty )$. The integral can be rewritten as:

$$\begin{array}{cc}\hfill \mathit{I}\left(\mathit{g}\right)& =\frac{1}{2}{\int }_{-\infty }^{\infty }{\int }_{\parallel \widehat{\mathit{u}}{\parallel }_2=1}\mathit{g}(\widehat{\mathit{x}}+{\mathit{S}}_x\mathit{u}){\left|r\right|}^{n_x-1}\mathcal{N}\{r;0,1\}d\widehat{\mathit{u}}dr\hfill \end{array}$$

(29)

In this stage, $\mathit{I}\left(\mathit{g}\right)$ can be regarded as the combination of radial integration with the form ${\mathit{I}}_r\left({\mathit{g}}_r\right)={\int }_{-\infty }^{\infty }{\mathit{g}}_r\left(r\right){\left|r\right|}^{n_x-1}\mathcal{N}\{r;0,1\}dr$ and spherical integration for the unit $n_x$-sphere with the form ${\mathit{I}}_s\left({\mathit{g}}_s\right)={\int }_{{∥\widehat{\mathit{u}}∥}_2=1}{\mathit{g}}_s\left(\widehat{\mathit{u}}\right)d\widehat{\mathit{u}}$. The stochastic radial integration rule (SRIR) and stochastic spherical integration rule (SSIR) can be obtained from the following two Lemmas [25]:

Lemma 1.

Assuming a set of weights ${\left\{{\omega }_{r,j}\right\}}_{j=0}^{m_r}$ are computed by the following equation:

$${\omega }_{r,j}={\int }_{-\infty }^{\infty }{\left|r\right|}^{n_x-1}\mathcal{N}\{r;0,1\}\prod _{t=0,t\ne j}^{m_r}\frac{r^2-{\rho }_t^2}{{\rho }_j^2-{\rho }_t^2}dr,j=0,1,\cdots ,m_r$$

(30)

where ${\rho }_0=0$, and the set of sample points ${\left\{{\rho }_j\right\}}_{j=1}^{m_r}$ are extracted from the joint PDF $p({\rho }_1,\cdots ,{\rho }_{m_r})\propto {\prod }_{j=1}^{N_r}{\rho }_j^{n_x+1}\mathcal{N}\{{\rho }_j;0,1\}{\prod }_{t=1}^{j-1}{({\rho }_j-{\rho }_t)}^2({\rho }_j+{\rho }_t)$. Then, an unbiased $d_r=2m_r+1$ degree SRIR can be obtained for ${\mathit{I}}_r\left({\mathit{g}}_r\right)$, that is:

$${\overline{\mathit{I}}}_r\left({\mathit{g}}_r\right)\approx \frac{1}{2}\sum _{j=0}^{m_r}{\omega }_{r,j}\left[{\mathit{g}}_r\left({\rho }_j\right)+{\mathit{g}}_r\left(-{\rho }_j\right)\right]$$

(31)

Lemma 2.

Let ${\overline{\mathit{I}}}_s\left({\mathit{g}}_s\right)={\sum }_{i=0}^{m_s}{\omega }_{s,i}{\mathit{g}}_s\left({\widehat{\mathit{u}}}_i\right)$ be a deterministic $d_s$-degree spherical integration rule for ${\mathit{I}}_s\left({\mathit{g}}_s\right)$, and $\mathit{Q}$ is a uniformly chosen orthogonal matrix. Then, ${\overline{\mathit{I}}}_{Qs}\left({\mathit{g}}_s\right)={\sum }_{i=0}^{m_s}{\omega }_{s,i}{\mathit{g}}_s\left(\mathit{Q}{\widehat{\mathit{u}}}_i\right)$ is an unbiased $d_s$-degree SSIR for ${\mathit{I}}_s\left({\mathit{g}}_s\right)$.

Note, the selection of weights ${\left\{{\omega }_{s,i}\right\}}_{i=0}^{m_s}$ and sample points ${\left\{{\widehat{\mathit{u}}}_i\right\}}_{i=0}^{m_s}$ is not unique; various deterministic spherical integration rules are given in [25].

In this paper, a three-degree SRIR and three-degree SSIR are adopted. Especially, the three-degree spherical integration rule in Lemma 2 is given by:

$${\overline{\mathit{I}}}_{3,s}\left({\mathit{g}}_s\right)=\frac{\parallel U_{n_x}\parallel }{2n_x}\sum _{i=1}^{n_x}({\mathit{g}}_s(-{\mathit{e}}_i)+{\mathit{g}}_s\left({\mathit{e}}_i\right))$$

(32)

where $U_{n_x}$ is the surface of the unit $n_x$-sphere; $\parallel U_{n_x}\parallel =2{\pi }^{n_x/2}/\mathrm{\Gamma }(n_x/2)$ is the surface content of $U_{n_x}$; ${\mathit{e}}_i$ denotes the unit vector with “1” in the i-th element and “0” in other positions.

The SSRIR can be obtained by constructing a product of a $d_r$-degree SRIR and a $d_s$-degree SSIR, that is:

$$\begin{array}{cc}\hfill {\overline{\mathit{I}}}_{}^{\left(d_r,d_s\right)}\left(\mathit{g}\right)& =\frac{1}{4}\sum _{j=0}^{m_r}{\omega }_{r,j}\sum _{i=0}^{m_s}{\omega }_{s,i}\left[\mathit{g}\left(\widehat{\mathit{x}}\pm {\rho }_j{\mathit{S}}_x\mathit{Q}{\widehat{\mathit{u}}}_i\right)+\overline{\mathit{g}}\left(\widehat{\mathit{x}}\pm {\rho }_j{\mathit{S}}_x\mathit{Q}{\widehat{\mathit{u}}}_i\right)\right]\hfill \\ & =\sum _{k=1}^{N^{SR}}{\omega }_k^{SR}\mathit{g}\left({\mathit{\chi }}_k\right)\hfill \end{array}$$

(33)

where $N^{SR}$ denotes the total number of the sample points; ${\mathit{\chi }}_k=\widehat{\mathit{x}}\pm {\rho }_j{\mathit{S}}_x\mathit{Q}{\widehat{\mathit{u}}}_i$ is the sample point with weight ${\omega }_k^{SR}={\omega }_{s,i}^{}{\omega }_{r,j}^{}/4$.

The approximation of $\mathit{I}\left(\mathit{g}\right)$ based on SSRIR with N iterations can be obtained as:

$${\mathit{I}}_N^{\left(d_r,d_s\right)}\left(\mathit{g}\right)=\frac{1}{N}{\sum }_{l=1}^N{\overline{\mathit{I}}}_l^{\left(d_r,d_s\right)}\left(\mathit{g}\right)$$

(34)

where ${\overline{\mathit{I}}}_l^{\left(d_r,d_s\right)}\left(\mathit{g}\right)$ denotes the l-th ${\overline{\mathit{I}}}_{}^{\left(d_r,d_s\right)}\left(\mathit{g}\right)$.

In theory, the statistics of function $\mathit{g}\left(\mathit{x}\right)$, including the mean, covariance and cross-covariance, should be separately approximated by SSRIR. However, in consideration of the computational cost, these statistics are commonly approximated in a single run of SSRIR [22]. For completeness, the detailed implementation of SSRIR with constant iterations is represented as follows:

• Initialize the number of iterations $l=1$. Define the maximum number of iterations as $N_{its}$. Set the initial mean value ${\mathit{\mu }}_g={\mathbf{0}}_{n_g\times 1}$, initial covariance ${\mathit{P}}_{gg}={\mathbf{0}}_{n_g\times n_g}$ and initial cross-covariance ${\mathit{P}}_{xg}={\mathbf{0}}_{n_g\times n_g}$; where $n_g$ is the dimension of $\mathit{g}\left(\mathit{x}\right)$.
• While $1\le l\le N_{its}$, use Equation (33) to generate the SSRIR-based sample points $\mathit{\chi }={\left\{{\mathit{\chi }}_{l\left(k\right)}\right\}}_{k=1}^{N_{}^{SR}}$ and corresponding weights ${\mathit{\omega }}_{}^{SR}={\left\{{\omega }_{l\left(k\right)}^{SR}\right\}}_{k=1}^{N_{}^{SR}}$ at the l-th iteration.
• Calculate the statistics of $\mathit{g}\left(\mathit{x}\right)$ at the current iteration:

  $${\mathit{\mu }}_g=\frac{1}{N}\sum _{l=1}^N\sum _{k=1}^{N_{}^{SR}}{\omega }_{l\left(k\right)}^{SR}\mathit{g}\left({\mathit{\chi }}_{l\left(k\right)}\right)$$

  (35)

  $${\mathit{P}}_{gg}=\frac{1}{N}\sum _{l=1}^N\sum _{k=1}^{N_{}^{SR}}{\omega }_{l\left(k\right)}^{SR}\mathit{g}\left({\mathit{\chi }}_{l\left(k\right)}\right)\mathit{g}{\left({\mathit{\chi }}_{l\left(k\right)}\right)}^T-{\mathit{\mu }}_g{\left({\mathit{\mu }}_g\right)}^T$$

  (36)

  $${\mathit{P}}_{xg}=\frac{1}{N}\sum _{l=1}^N\sum _{k=1}^{N_{}^{SR}}{\omega }_{l\left(k\right)}^{SR}{\mathit{\chi }}_{l\left(k\right)}\mathit{g}{\left({\mathit{\chi }}_{l\left(k\right)}\right)}^T-{\mathit{\mu }}_x{\left({\mathit{\mu }}_g\right)}^T$$

  (37)
• Set $l=l+1$, and shift to Step 2.

For simplicity, the implementation of SSRIR can be compactly denoted as:

$$\left[{\mathit{\mu }}_g,{\mathit{P}}_{gg},{\mathit{P}}_{xg}\right]=\mathrm{SSRIR}\left[\mathit{g}\left(·\right),\widehat{\mathit{x}},\mathit{P}\right]$$

(38)

If the process model in (9) is nonlinear, Equations (18) and (19) can be substituted with:

$$\left[{\widehat{\mathit{x}}}_{k|k-1},\left({\mathit{P}}_{k|k-1}-{\mathit{Q}}_{\omega ,k-1}\right)\right]=\mathrm{SSRIR}\left[\mathit{f}\left(·\right),{\widehat{\mathit{x}}}_{k-1|k-1},{\mathit{P}}_{k-1|k-1}\right]$$

(39)

Analogously, if the observation model in (9) is nonlinear, Equations (22)–(24) can be substituted with:

$$\left[{\widehat{\mathit{z}}}_{k|k-1},\left({\mathit{P}}_{zz,k|k-1}-{\mathit{R}}_{v,k}\right),{\mathit{P}}_{xz,k|k-1}\right]=\mathrm{SSRIR}\left[\mathit{h}\left(·\right),{\widehat{\mathit{x}}}_{k|k-1}^{},{\mathit{P}}_{k|k-1}\right]$$

(40)

Eventually, Equations (20), (21), (25), (39) and (40) are composed of the integrated SIH${}_{\infty }$F method. The flowchart of the SIH${}_{\infty }$F algorithm is shown in Figure 3.

Figure 3 implies that the SSRIR is incorporated with the framework of DFH${}_{\infty }$F so as to approximate the nonlinear function in the model description. Compared with the traditional SIF, the SIH${}_{\infty }$F introduces the cost function described in Equation (10) and aims to estimate the state under the worst initial error and model uncertainty.

4. Numerical Simulation

Numerical simulation is an essential approach to analyze the performance of the proposed method. As a prerequisite, a flight trajectory is designed to generate the simulated data of MINS and SINS. Then, the simulated data of the two INSs are processed by different filtering algorithms for comparison.

4.1. Designs of the Numerical Simulation

In this subsection, a typical flight trajectory with a turning maneuver is designed, as shown in Figure 4.

Firstly, the aircraft flies 20 s with a constant velocity, then turns anticlockwise ${60}^{\circ }$ in 60 s and continually flies 20 s with a constant velocity. Note that the roll angle and yaw angle of the aircraft are changing simultaneously during the turning maneuver. The specifications of the simulated flight trajectory are given in

Table 1. The specifications of simulated flight trajectory.

Initial Attitude		Initial Velocity		Initial Position
Yaw	60${}^{\circ }$	East	−86.6 m/s	Latitude	116${}^{\circ }$ N
Pitch	0${}^{\circ }$	North	50 m/s	Longitude	${40}^{\circ }$ E
Roll	0${}^{\circ }$	Up	0 m/s	Height	1000 m

.

The initial actual physical misalignment angles are set to ${\mathit{\psi }}_a=[5^{\circ };5^{\circ };{80}^{\circ }]$ with covariance ${\mathit{Q}}_a=diag\{{({0.02}^{\circ }/h)}^2;{({0.02}^{\circ }/h)}^2;{({0.02}^{\circ }/h)}^2\}$. The errors of inertial instruments in SINS are listed as follows: the constant drift of the gyroscope is $5^{\circ }/h$, and the random walk of the gyroscope is $1^{\circ }/\sqrt{h}$; the constant bias of accelerator is 0.2 mg, and the random walk of accelerator is 0.02 m/s/$\sqrt{h}$. The attitude precision and velocity precision of MINS are ${0.1}^{\circ }$ and 0.02 m/s, respectively. The observation noise with uncertain statistics is achieved by injecting non-Gaussian noise in the simulated observations. The disturbing noise for attitude observation is assumed to satisfy the uniform distribution within $[-{0.05}^{\circ }{0.05}^{\circ }]$. Similarly, the disturbing noise for velocity observation is assumed to satisfy the uniform distribution within [−0.01 m/s 0.01 m/s].

4.2. Results and Discussion

The simulated data are processed by SIF, CH${}_{\infty }$F and SIH${}_{\infty }$F, respectively. The error attenuation parameter for CH${}_{\infty }$F and SIH${}_{\infty }$F is set to $\gamma =0.08$. The number of iterations of SIF and SIH${}_{\infty }$F is $N_{its}=20$. These methods were independently performed with 50 Monte Carlo simulations, and the root-mean-square error (RMSE) of the three misalignment angles was utilized to evaluate the performance of different filtering algorithms with statistical significance. The estimation results are shown in Figure 5, Figure 6 and Figure 7.

From Figure 5 and Figure 6, it can be seen that all three methods have a quick convergence performance for the estimation of horizontal misalignment angles, i.e., ${\psi }_{ax}$ and ${\psi }_{ay}$. Meanwhile, the estimation results of ${\psi }_{az}$ are partially estimated at the beginning of RTA and quickly converged to true values while the turning maneuver of aircraft is carried out.

Among these methods, SIF gives poor performance in terms of all three misalignment angles, as it requires accurate statistics of the noise. However, the existence of non-Gaussian noise in the observation degrades the performance of SIF. Especially, Figure 5 and Figure 6 show that SIF cannot suppress the negative effect of uncertainties at the beginning of the turning maneuver (20 s). In contrast, CH${}_{\infty }$F has better performance than SIF, since the combination of the H${}_{\infty }$ technique and cubature integration rule simultaneously suppressed the negative effect of uncertainties and nonlinearities. Even so, the estimation given by CH${}_{\infty }$F fluctuates wildly for the estimation of ${\psi }_{az}$, as is shown in Figure 7. This phenomenon demonstrates that the local validity of the cubature integration rule cannot guarantee the consistency of the estimation error with the existence of uncertainty and nonlinearity. The SIH${}_{\infty }$F provides the highest estimation accuracy compared with other methods. Besides, the estimations of SIH${}_{\infty }$F are obviously smoother than CH${}_{\infty }$F. This demonstrates that the SSRIR can approximate the nonlinear function more accurately than the cubature integration rule, and the boundedness of estimation errors can be guaranteed. As a result, the performance of SIH${}_{\infty }$F outperforms the CH${}_{\infty }$F and SIF in the statistical sense.

The statistics of RMSE in the last 20 s are used to further evaluate the accuracy and stability of different filtering algorithms. The mean values (in the time-averaged sense) and standard deviation (STD) values of 50 dependent simulations are given in

Table 2. Statistics of RMSE for 50 Monte Carlo simulations (${}^{\circ }$).

	Statistics	SIF	CH${}_{\mathit{\infty }}$F	SIH${}_{\mathit{\infty }}$F
${\psi }_{ax}$	Mean	0.2101	0.1379	0.1122
	STD	0.0481	0.0402	0.0397
${\psi }_{ay}$	Mean	0.2034	0.1489	0.1167
	STD	0.0468	0.0421	0.0411
${\psi }_{az}$	Mean	2.7150	2.1596	1.8213
	STD	0.1460	0.4266	0.4134

.

From Table 2, it is obvious that the results given by SIH${}_{\infty }$F are superior to other methods in terms of both mean and STD values, since it possesses a high ability to address the coexistence of nonlinearity and uncertainty.

The major advantage of SIH${}_{\infty }$F is the ability to address the estimation problem with significant nonlinearity. To verify the performance of the proposed method, different cases with increasing degrees of nonlinearity are designed. The initial values of ${\mathit{\psi }}_{ai},(i=x,y,z)$ are set to ${20}^{\circ }$–${50}^{\circ }$ with ${10}^{\circ }$ intervals, and the RTA based on CH${}_{\infty }$F and SIH${}_{\infty }$F were independently carried out 50 times for each case. The statistics of RMSE, including the mean and STD for the last 20 s, are utilized to evaluate the performance of the filters. The resulting statistics of the RMSE are represented in Figure 8, Figure 9 and Figure 10.

As shown in Figure 8, Figure 9 and Figure 10, the estimation errors of CH${}_{\infty }$F obviously increase with the growth of initial actual misalignment angles. The result implies that the cubature integration rule based on the deterministic sampling strategy deteriorates the approximation accuracy of Gaussian weighted nonlinear integrals. In contrast, SIH${}_{\infty }$F performs better than CH${}_{\infty }$F, since the SSRIR provides asymptotically unbiased approximations for Gaussian weighted nonlinear integrals. In other words, SIH${}_{\infty }$F alleviates the impact of significant nonlinearity, and the results demonstrate the validity of the proposed method.

5. Van Test

To further verify the performance of the investigated method, a van test was carried out. In this van test, the RTA is implemented while MINS and SINS are operated on a moving base. SINS outputs the raw data of inertial sensors, including the angular velocity and specific force, while MINS provides the observation data, i.e., the attitude and velocity. By using different filtering algorithms, the RTA procedure provides the estimation results of misalignment angles in real time. Besides, an accuracy assessment procedure is also necessitated to evaluate the performance of different filters. The structure diagram of the van test is displayed in Figure 11.

5.1. Specifications of the Van Test

The test system contains a MEMS-based INS (NAV4400), an attitude and heading reference system (XW-ADU7612) and a laptop. The NAV440 is utilized as SINS to provide raw data of inertial sensors, and its dominating errors are listed as follows: the in-run stability of the gyroscope is ${10}^{\circ }/h$; the random walk of the gyroscope is ${4.5}^{\circ }/\sqrt{h}$; the in-run stability of the accelerometer is 1 mg; the random walk of the accelerometer is 1 m/s/$\sqrt{h}$. The XW-ADU7612 is taken as MINS. It contains an inertial measurement unit, whose accumulative errors are corrected by a GPS module with two external antennas. The attitude precision and velocity precision of MINS are ${0.2}^{\circ }$ and 0.02 m/s, respectively. The two INSs were bolted on a steel plate, and the GPS antennas were installed on top of the van. The update rates of MINS and SINS are 10 Hz and 100 Hz, respectively. The laptop is used for restoring the output data of the two INSs. The configurations of the van test are shown Figure 12.

The van test was performed in the northern fifth ring road of Beijing city, China. The trajectory of the van test is shown in Figure 13, in which the blue line denotes the trajectory for the implementation of RTA. The RTA was initiated at Point A and ended at Point B, and the execution time was 120 s.

The GPS signal of MINS was intentionally interrupted two times during the test, that is the observations can be deduced to contain non-Gaussian noise. The number of captured satellites by the GPS antennas is shown in Figure 14.

5.2. Accuracy Assessment, Test Results and Discussions

For the purpose of evaluating the estimation accuracy of different methods, high precision reference misalignment angles are required. However, the benchmark of misalignment angles cannot be directly obtained in real time, and thus, the off-line accuracy assessment approach is necessitated. In this test, the relative attitude between MINS and SINS is invariable since they were rigidly bolted on the steel plate, i.e., ${\psi }_a$ can be considered as a constant vector. Therefore, the CKF-based fixed-point smoother was employed in the accuracy assessment procedure, on account of the nonlinear characteristic of the model description [27,28]. As is shown in Figure 13, the accuracy assessment procedure was initiated at Point B and ended at Point C (700 s), where the GPS functioned well and the Gaussian assumption for the observation noise is accredited. The smoothing estimations of misalignment angles are ${\overline{\psi }}_{ax}={0.18}^{\circ }$, ${\overline{\psi }}_{ay}={0.9}^{\circ }$ and ${\overline{\psi }}_{az}=-{90.5}^{\circ }$.

In this test, the comparison between SIF, CH${}_{\infty }$F and SIH${}_{\infty }$F was executed. The error attenuation parameter for CH${}_{\infty }$F and SIH${}_{\infty }$F was set to $\gamma =0.6$. The the number of iterations of SIF and SIH${}_{\infty }$F was $N_{its}=20$. The estimation errors of the three methods are represented in Figure 15, Figure 16 and Figure 17.

Figure 15, Figure 16 and Figure 17 clearly demonstrate that the estimation accuracy of SIH${}_{\infty }$F outperforms SIF and CH${}_{\infty }$F in terms of all three misalignment angles. The oscillating amplitude of estimations given by CH${}_{\infty }$F is larger than that of SIH${}_{\infty }$F, since the approximation accuracy of the cubature integration rule is lower than that of SSRIR. The SIF cannot attenuate the influence caused by uncertain statistics of observation noise and, thus, gives the lowest estimation accuracy. After 120 s, the estimate errors of the three methods are shown in

Table 3. Estimation error of the van test (${}^{\circ }$).

	SIF	CH${}_{\mathit{\infty }}$F	SIH${}_{\mathit{\infty }}$F
${\psi }_{ax}$	0.3750	0.2284	0.1145
${\psi }_{ay}$	−0.8059	−0.5242	−0.3731
${\psi }_{az}$	−0.9873	−0.6494	−0.4253

.

From Table 3, one can notice that the estimation precision of SIH${}_{\infty }$F is higher than SIF and CH${}_{\infty }$F, in which the precision of ${\psi }_{ax}$ increases by 69.5% and 49.9%, respectively, the precision of ${\psi }_{ay}$ increases by 53.7% and 28.8%, respectively, and the precision of ${\psi }_{az}$ increases by 56.9% and 34.5%, respectively.

6. Conclusions

The coexistence of large initial attitude errors and uncertain statistics of observation noise severely affects the accuracy of RTA. The major contribution of this paper is to develop a SIH${}_{\infty }$F for the purpose of improving the performance of RTA with considerations of accuracy and robustness. For the first time, the SSRIR is combined with the framework of DFH${}_{\infty }$F, and the resulting SIH${}_{\infty }$F can theoretically provide asymptotically unbiased estimation for a nonlinear system with high robustness. The numerical simulation and van test were carried out to compare the performance of the SIH${}_{\infty }$F, CH${}_{\infty }$F and traditional SIF. The results show that the newly-proposed method has the highest precision of the estimation results and robustness to the unknown statistics of observation noise. Promoted by its advantages, SIH${}_{\infty }$F can be considered as a competitive candidate method for the practical application of RTA.