Design and evaluation of INS/GNSS loose and tight coupling applied to launch vehicle integrated navigation

Nomenclature

1.   Introduction

Affordable modern navigation solutions integrates an Inertial Navigation System (INS, see ^[1,2]) with a receiver of a Global Navigation System by Satellites (GNSS). This integration (or coupling) improves the quality of the navigation even if the INS has a low quality Inertial Measurement Unit (IMU, with three or more gyroscopes and accelerometers), due to the fusion between the inertial navigation estimates and the independent GNSS receiver measurements. The vehicle's guidance and control (^[3,4,5]) and the flight safety ^[6] are users of these integrated navigation outputs, being required to have navigation solutions independent from the one being used for the guidance and control. In this context, the cost of the whole system is significatively reduced by allowing lower cost IMUs, which motivated the assessment of reliable INS/GNSS integration algorithms presented in this work.

The Extended Kalman Filter (EKF) has been used since its formulation for spacecraft navigation, and is typically used to combine the INS and GNSS receiver outputs. Other integration schemes have been also considered in the literature, including the Robust Unscented Kalman Filter (RUKF), the derivative UKF and the direct filtering approaches ^[7,8,9,10]. Figure 1 shows the Navigation functional block in the context of a Navigation, Guidance and Control (NGC) system of an aerospace propelled vehicle. The Guidance block compares the Trajectory block output (i.e. the desired trajectory $\vec{p}_d(t), \vec{v}_d(t)$ with initial conditoin $\vec{p}(t_0), \vec{v}(t_0)$) with the Navigation block output (i.e., the current kinematic state $\vec{p}(t), \vec{v}(t)$) and generates an output in terms of attitude reference ($\vec{\theta}_d(t), \vec{\omega}_d(t)$) to the (attitude) Control block and other output to the Propulsion/Staging block (i.e. the main engine start or cut-off events, and stages/fairing separation events). The Control block compares the guidance attitude reference with the current navigated attitude ($\vec{\theta}(t), \vec{\omega}(t)$) and computes the commands to be implemented by the actuators. The rotation and translation dynamics are determined by the Guidance, Control and derived actuation allocation functions as well as the flight environment, and generate the (unknown) physical state for the inertial sensors and GNSS receivers whose measurements are the input to the Navigation block.

A satellite launch vehicle has unique features due to its wide operational range in terms of acceleration, velocity, position, temperature, vibrations, antenna visibility, etc. For this reason, selecting constant covariance matrices for the process noises (covariance $Q$) and the external measurement noises (covariance $R$) for the EKF implementation in the Navigation block may not be realistic beyond some suborbital flights ^[11]. Due to these particular flight dynamics, the system noise statistics may be unknown and time varying, which requires specific filtering frameworks ^{[12,13,14,15]}. The measurement uncertainties are also affected by the delay of the GNSS receiver solution availability, which is more significant for higher accelerations and velocities. On ^[16] we showed convenient to use an adaptive version of the EKF to compensate these effects, which will be also used here on the GNSS receiver filtered solution. Moreover we propose a compensation for the tropospheric delay effect on the GNSS pseudoranges during the ascent through different layers of the atmosphere, including negative elevations for the line of sight which are progressively incorporate during ascent. Finally, to assess the robustness of the navigation we simulate a model of variation of IMU biases during ascent (with increasing temperature) and the possible failure on the GNSS receiver and/or antennas.

In this work we are interested in the design and evaluation of possible coupling strategies. A comparison is made between the following integrated navigation alternatives:

● Loosely coupled INS/GNSS, using as external measurements the PVT fix solution (Position, Velocity, Time) provided by the GNSS receiver (at least four satellites are necessary, see ^[17]).

● Loosely coupled INS/GNSS, using as external measurements a PVT filtered solution provided by the GNSS receiver, which assumes a model of the movement to smooth and robustify the solution. This filter on the receiver will be called EKF-GNSS. It is not necessary to have at least four satellites for a valid solution, as the EKF-GNSS propagates it, but this solution will have an error as the internal model does not know the true accelerations and angular velocities. On the other hand, there is a loosely coupled EKF filter to combine this solution and its covariance with the INS output.

● Tightly coupled INS/GNSS, using as external measurements the observables (pseudorange and delta-pseudorange) and satellite ephemeris from the GNSS receiver. It is not necessary to have four satellites and with at least one satellite the innovations will be integrated with the inertial navigation, hence this option is a priori more efficient in terms of the use of the available information.

The solutions given by all these strategies are compared using the numerical simulator framework described in ^[3], on which it is possible to simplify the evaluation of all the estimation errors and produce the failure scenarios as desired (as presented in ^[18]). This comparison is shown in Figure 2, where EKF-LC is the loosely coupled extended Kalman filter and EKF-TC is the tightly coupled extended Kalman filter. The bottom layer on this figure defines the type of GNSS receiver output, which can be a simple PVT model (white noise plus the simulated variables) or the lower level information generated by the receiver (ephemeris, pseudorange, delta-pseudoranges, etc). The middle layer is only active when the receiver delivers high level outputs as the PVT fix solution (computed with the Bancroft algorithm) or the PVT filtered solution obtained with the EKF-GNSS. Finally, the third (and higher) level contains the integrated navigation blocks (EKF-LC and EKF-TC). This comprehensive evaluation is the main achievement of the current work, including the simulation of launch vehicle's dynamics, GNSS receiver solutions, observables generation, associated ephemeris and main IMU features, which will be augmented with hardware in the loop evaluations on a future work. Other information fusion techniques have been also considered in the literature, including nonlinear filtering ^[19,20,21], which is beyond the scope of this work and will be compared in future research.

Section 2 provides the measurement model for the GNSS observables in the receiver (pseudorange, delta-pseudorange and user clock) and Section 3 provides a deterministic compensation model for the atmospheric and relativistic effects. A summary of the EKF-LC applied to a launch vehicle is given in the Appendix. It includes a description of the kinematic model of the vehicle with measurements of acceleration (from accelerometers), angular velocity (from gyroscopes), position and velocity (from the GNSS receiver solution).

On Section 4, an algorithm is proposed to compute the PVT solution using the EKF-GNSS filter, and compared with respect to the typical fix solution obtained with the well-known Bancroft algorithm (see ^[17]), which requires at least four visible satellites to obtain a result. The EKF-GNSS uses adaptive covariance matrices as a function of the last velocities to obtain a coarse estimate of the acceleration, and provides state variable covariances to describe the PVT filtered solution to be integrated as external measurements by the EKF-LC (i.e. a time varying $R$ matrix).

Section 5 develops the EKF-TC for the proposed tightly coupled INS/GNSS. The main advantage of this coupling is that the information is exploited more efficiently, as the observables are used even when the number of satellite is not enough to compute the PVT fix solution. The communication delay compensation for EKF-TC is provided in subsection 5.4.

These methods are evaluated by simulations in Section 6, showing the advantage of the PVT filtered solution over the PVT fix solution, for the loosely coupled INS/GNSS. Moreover, the tightly coupled INS/GNSS is shown better than the other options, and is finally validated with a test using hardware in the loop. The conclusions are provided in Section 7.

2.   GNSS receiver measurement model

The GNSS-based navigation requires to develop a measurement model for the pseudoranges, delta-pseudoranges and the user clock dynamics, which are used both for the analytic description and the onboard prediction of the measurements. For the onboard evaluation of these models, we consider that the receiver obtains the ephemeris of the GNSS satellites from the signal code. For each satellite we can compute the position $\underline{p}_{s, k}^e$ and velocity $\underline{v}_{s, k}^e$, both written in ${\mathcal E}$ frame at the reception time $t^R_k$. As the information in the signal code is given on ${\mathcal E}$ frame at the transmission time of each satellite, it is necessary to correct this frame using the Earth angular velocity and the propagation time. In order to build the receiver observables, we need first to define the user (i.e. receiver) clock dynamics in terms of a bias and drift.

2.1.   User clock dynamics

The clock bias $b_t$ and drift $d_t$ can be described with $\underline{x}_c: = c[b_t\, d_t]^T$ with the continuous model ^[22]:

where $\underline{\epsilon}_{c}$ is a zero-mean gaussian white noise with autocovariance $\mathbf{E}\left[ \underline{\epsilon}_c(t) \underline{\epsilon}_c^T(\tau) \right] = \delta(\tau - t)Q_c$, being:

The coefficients $h_0$ and $h_{-2}$ can be obtained as a function of the Allan variances ^[39] and can be found for typical GNSS clocks as the high quality Temperature Compensated Crystal Oscillator in ^[41]:

In the following we will construct the model for the pseudorange and delta-pseudorange using the evaluation of the user clock at the reception times $t^R_k$, hence we define $b_{t, k} = b_t(t^R_k)$ and $d_{t, k} = d_t(t^R_k)$.

2.2.   Pseudorange measurement

Let $\tilde{\rho}_{s, k}$ be the pseudorange (i.e. computed distance using GNSS measurements) between the vehicle and the GNSS satellite $s$, which is computed with the transmission time $t_{s, k}^T$, and the time of reception $t_k^R$ at the $k$-th measurement instant. Each $t_{s, k}^T$ is measured using a GNSS time system:

where $\tau_{s, k}$ is the propagation time and $\eta_{s, k}^T$ is a measurement zero-mean gaussian noise. The time $\tau_{s, k}$ is determined by the vacuum propagation time $\tau_{s, k}^G: = R_{s, k}/c$, where $R_{s, k}$ is the distance between the user location at reception time $t_k^{R}$ and the satellite $s$ location at transmission time $t_{s, k}^T$, plus the signal delay times $\tau_{s, k}^I$, $\tau_{s, k}^T$ and $\tau_{s, k}^R$, corresponding respectively to the ionospheric, tropospheric and relativistic delays (see ^[3]):

The time $t_{k}^R$ is measured with a local clock modelled with an unknown bias $b_{t, k}$ relative to the GNSS system time, which is time-varying with the unknown clock dynamic. The receiver clock bias is common for all the pseudoranges, provided that all the measurements are made at the same time, and the measured reception time is $t_k^R + b_{t, k}$. Hence the pseudorange measurement $\tilde{\rho}_{s, k}$ is defined as:

with $c$ the speed of light in vacuum. Let $\Delta\rho_{s, k}$ be the sum of atmospheric and relativistic delays, hence:

where $R_{s, k} = \|\underline{p}^e_{s, k}-\underline{p}^e_k\|$. Finally, for certain model ${\hat \Delta}\rho_{s, k}$ of $\Delta\rho_{s, k}$ it will be useful to define the compensated pseudorange measurement:

2.3.   Delta-pseudorange measurement

The Doppler effect makes possible to construct the delta-pseudorange $\tilde{\dot{\rho}}_{s, k}$, as a function of the difference $\tilde{\Delta} f_{s, k}$ between the measured carrier frequency from the $s$-th GNSS satellite $f_{R, k}^s$ and its nominal value $f^s_0$:

where $\eta_{s, k}^{f}$ is a zero-mean white gaussian measurement noise and $d_{t, k}$ is the clock drift term.

Let $\lambda_s = c/f_0^s$ be the carrier wave lenght, hence:

where we have used $\dot{R}_{s, k} = \lambda_s(f^s_0 - f_{R, k}^s)$, which is known as the Doppler effect (see ^[40]). Moreover, $\dot{R}_{s, k}$ is also given by the dot product between the relative velocity vector $\underline{v}_{rs, k}^e$ and the unit vector $\underline{\breve{e}}^e_{s, k}$ which points along the line of sight between the user and the $s$-th satellite, i.e.:

with:

and ^[18]:

3.   Delay compensation of measured observables

The compensation (2.11) requires to estimate $\Delta \rho_{s, k}$ onboard with models for the troposphere, ionosphere and relativistic effects. A simple tropospheric model will be extended for smaller and negative elevations, while for the ionosphere we will use the Klobuchar model ^[22,33]. A simple relativistic model is also provided for compensation.

3.1.   Tropospheric delay compensation

A mapping function of the tropospheric delay is given in ^[22] as a function of the user altitude $h_{v, k}$ and the line of sight elevation $\varepsilon_{s, k}$ relative to the $s$-th satellite at time $t_k$. For the elevation range $5^\circ \leq \varepsilon_{s, k} \leq 90^\circ$ the following function provides the delay in meters:

where $h_\tau : = 7518.8\; \mathrm{m}$, $\phi_\alpha: = 0.0121$ and $k_\alpha: = 2.47\; \mathrm{m}$. As for $\varepsilon_{s, k} < {\varepsilon}_\alpha = 5^\circ$ the model is not defined, lets consider the following function for all the positive elevations:

which now is valid for $\Delta \rho^{Tm}$ in the range $0 \leq \varepsilon_{s, k} \leq 90^\circ$.

For negative angles $\varepsilon_{\mathrm{min}}(h_{v, k}) \leq \varepsilon_{s, k} \leq 0$, where $\varepsilon_{\mathrm{min}} : = -\mathrm{acos}\left(R_e/(R_e+h_{v, k})\right) \leq 0$ is the minimum feasible angle not blocked by the Earth with local radius $R_e$. Here we follow the same idea taken in the simulator implementation ^[3], which composes the delay for negative angles as shown in the Figure 3. The total delay is given by the sum of delay segments '1-2'+'2-3'+'3-$s$', hence:

The point 2 was defined such that the elevation is zero which corresponds to the altitude $\bar{h}_{s, k} : = \bar{h}(\varepsilon_{s, k}, h_{v, k})$ such that:

where $R_e$ is the local Earth radius. The tropospheric delay is computed using (3.3):

On the segment between the point 3 and the satellite, the delay is computed as:

With previous definitions:

and (3.4) is computed for negative values of $\varepsilon_{s, k}$ as:

Considering (3.2), (3.3) and (3.11), for all $\varepsilon_{\mathrm{min}}(h_{v, k}) \leq \varepsilon_{s, k} \leq 90^\circ$ the proposed model for the Tropospheric Delay Compensation (TDC) becomes:

Figure 4 shows $\hat \Delta \rho_{s, k}^{T}: = \Delta \rho^{T}\left(\varepsilon_{s, k}, h_{v, k}\right)$ as a function of $\varepsilon_{s, k}$ and parameterized for different positive altitudes $h_{v}$.

Finally, a better result was obtained, by comparison with the high fidelity model given in ^[3], by replacing the gain parameter $k_\alpha$ in (3.1) and (3.3) by a gain function $k(\varepsilon_{s, k}, h_{v, k})$ as follows:

3.2.   Ionospheric delay compensation

The compensation uses the model proposed by Klobuchar ^[33] for the Ionospheric Delay Compensation (IDC), ${\hat \Delta} \rho^I_{s, k}$, as a function of coefficients broadcasted by the GPS constellation. This model does not consider the variation in altitude but only the elevation of the line of sight.

3.3.   Relativistic delay compensation

The relativistic delay affects the measurement of the pseudoranges, and is included on the estimate $\hat{\Delta} \rho_{s, k}$. It is well known that part of this delay is generated by orbit eccentricity ^[3]. This part of the effect can be modelled as a function of the ephemeris of each observed satellite $s$ at time $t_k$:

4.   GNSS receiver PVT outputs for loosely coupled integrated navigation

It is well known that to obtain the PVT fix solution it is necessary to have at least four GNSS satellites (^[17,22]), being possible to validate it with more than four GNSS satellites in view. A more accurate solution can be obtained by filtering the GNSS observables to obtain the PVT filtered solution, using certain hypothesis on the user motion dynamics as considered in ^[23,24]. To model these PVT solutions, we assume that there are $N_k$ satellites acquired by the receiver on times $t_k = t_k^R$, with the associated ephemeris, pseudorange and delta-pseudorange. At each time, the vector $\underline{\tilde{\rho}}_k \in \mathbf{R}^{N_k}$ contains the measured pseudoranges and the vector $\underline{\tilde{\dot{\rho}}}_k \in \mathbf{R}^{N_k}$ contains the measured delta-pseudoranges. Also, the vector $\underline{\tilde{E}}_k$ contains the ephemeris position and velocity coordinates of the $N_k$ satellites.

4.1.   Point solution (fix) algorithm

The equations relating the nonlinear GNSS measurements with the unknown position and clock bias can be solved iteratively in steps $i: = 1, ..., M_k$ using Newton method ^[17], and provided $N_k \geq 4$. In order to avoid this complexity, an approximate algebraic method has been typically adopted ^[22]. The solution proposed by Bancroft in ^[17] re-states the problem in terms of a second order polynomial, by using the same linear functional as defined by Lorentz on relativistic mechanics ^[25]. The PVT point solution $\underline{p}^P_k$ on frame $\mathcal E$ and the clock bias $b_{c, k}^P$, for time $t_k$, are given by the Bancroft algorithm as:

where $\underline{y} : = [y_1\, y_2\, y_3\, y_4]^T$ is given by one of the following options:

being $\lambda_{1, 2}$ the solutions at time $t_k$ of a second order polinomial. The vectors $\underline{u}, \underline{v}\in \mathbf{R}^4$ are determined by the pseudoranges and ephemeris of the GNSS satellites observed at that moment. A test of residuals selects which of the two solutions in (4.2) is the correct one and describes it with a quality number $\mathrm{q}_k^P \in \mathbf{N}$ which measures the degree of accuracy of the point solution. The point velocity solution $\underline{v}^P_k$ for time $t_k$ in frame $\mathcal E$ and the associated clock drift $d^P_{c, k}$ can be obtained in a similar way ^[26].

4.2.   Filtered solution algorithm

Nomenclature:

Here we describe a PVT filtered solution which is typically used on GNSS receivers to improve the accuracy, based on certain model of the vehicle. The particular dynamics of a satellite launcher makes convenient to use as internal states of the filter the user position, velocity and acceleration, and the user clock bias and drift ^[24]. This includes the adaptation of covariances as proposed in ^[22].

4.2.1.   Extended kalman filter for GNSS positioning

We take the following discrete model for the position, velocity and acceleration dynamics of the satellite launch vehicle:

where $\Delta_T = 1\; \mathrm{sec}$, the super-index $F$ stands for the receiver filter solution, $\underline{p}^F_k$, $\underline{v}^F_k$ and $\underline{a}^F_k$ are the position, velocity and acceleration to be estimated, while $\underline{\dot{a}}^F_k$ is a random discrete zero-mean white noise gaussian model with covariance matrix $Q_k^{\dot{a}} : = {\bf E}[\underline{\dot{a}}^F_k \underline{\dot{a}}^{F\, T}_k]$. Let $\underline{x}^{F}_k : = [\underline{p}^{F\, T}_k\, \underline{v}^{F\, T}_k\, \underline{a}^{F\, T}_k]^T$ and:

being:

with $\mathbf{U}_{9 \times 9}^m \in \mathbf{R}^{9 \times 9}$ a matrix with $\delta_{i+m, j}$ on its $i$-th row and $j$-th column, being $\delta_{i+m, j}$ the Kronecker delta, such that $\delta_{i+m, j} : = 1$ when $i+m = j$, or $\delta_{i+m, j}: = 0$ otherwise.

A simplified model of the receiver clock is defined here with a bias $b_{c, k}^F$ and drift $d_{c, k}^F$ such that:

where $\dot{d}_{c, k}^F$ is a zero-mean random process with variance $\sigma^2_{\dot{d}, k}$. Let $\underline{b}^F_{c, k} : = [{b}_{c, k}^F\, d_{c, k}^F]^T$ be the clock model state vector to be estimated, which obeys the model:

where:

The full state vector of the filter is $\underline{\chi}^F_k : = [\underline{x}_k^{F\, T}\, \underline{b}_{c, k}^{F\, T}]^T$ and satisfy:

with $\underline{\xi}_k^F : = [\underline{\dot{a}}_k^{F\, T}\, \dot{d}_{c, k}^F]^T$, where $\underline{\dot{a}}_k^F$ and $\dot{d}_{c, k}^F$ are not correlated, and:

The covariance matrix of $\underline{\xi}^F_k$ is:

and is related with the covariance of the last term in (4.13) through:

The measurement model for the pseudoranges and delta-pseudoranges of the $N_K$ satellites seen at $t_k$ is:

being $\underline{\eta}_k^F$ a zero-mean measurement noises vector and $R^F_k : = {\bf E}[\underline{\eta}_k^F\underline{\eta}_k^{F\, T}]$, while ^[22]:

where $\underline{e}_{s, k}^F$, $\underline{\breve{e}}_{s, k}^F$ and $\underline{v}_{rs, k}^F$ are defined as in (2.18) and (2.19), as a function of $\underline{\chi}^F_k$ and satellite $s$.

The EKF-GNSS filter is depicted in Figure 5, where the a priori estimate $\underline{\hat \chi}_k^{F, (k-1)}$ and the a posteriori estimate $\underline{\hat \chi}_k^{F, k}$ have respectively the error covariance matrices $P_k^{F, (k-1)}$ and $P_k^{F, k}$, while the process noises covariance ${\hat Q}_k^{\xi_F}$ and the measurement noises covariance ${\hat R}_k^{F}$ will be computed for each time $t_k$ on frame $\mathcal E$. Let $\delta \underline{y}^N_k \in \mathbf{R}^{2N_k}$ be the innovations vector at $t_k$ for the $N_k$ observed satellites. A validation logic is used for each innovation, rejecting wrong measurements as a function of the absolute value of each innovation component. Let $\delta \underline{y}_k^F \in \mathbf{R}^{2N_k^F}$ be the vector of innovations for $N_k^F \leq N_k$ satellites validated at $t_k$. The filtered solution to be used as external measurement by the EKF-LC is composed by the position $\underline{\tilde p}^e_k : = \underline{\hat p}^{F, k}_k$ and the velocity $\underline{\tilde v}^e_k : = \underline{\hat v}^{F, k}_k$. The receiver output also includes a quality index $\mathrm{q}^F_k$, the number of observed satellites, and the estimation error covariance $\hat{R}^{pv}_k$ obtained from the matrix $P_k^{F, k}$.

The output Jacobian matrix for $h_k^F(\underline{\chi}^F_k)$ becomes:

being:

with:

Notice that the matrix $B^F$ has more rows than columns, hence it is not possible to guarantee the stochastic exponential convergence using the condition in ^[29] given as $Q^{\chi_F}_k = B^FQ^{\xi_F}_kB^{FT} > \underline{q} I$ for $\underline{q} > 0$. However, following the Remarks after Theorem 3.1 in ^[29], it is not actually necessary, to obtain error bounds, that $Q^{\chi_F}_k$ were given by the covariance of the noise term, although any other positive definite matrix could be chosen as well. In ^[30] the target tracking problem is considered with an input matrix having more rows than columns in the same way.

4.2.2.   Adaptive $\hat{Q}^{\xi_F}_k$ and $\hat{R}^F_k$

During the ascent of the launch vehicle, the acceleration grows as the mass of propellant is burned, while the thrust also grows due to the less dense atmosphere. Therefore it is good to estimate the acceleration in order to bound the process noises for the EKF-GNSS. The acceleration in ${\mathcal E}$ frame will be estimated using the last two point solutions at $t_{k-1}$ and $t_{k}$:

The process noise covariance becomes:

where $\hat{\sigma}^2_{{\ddot b}_c} > 0$ is a constant parameter associated to the deviation of the clock drift derivative and:

where $\alpha_Q \in \mathbf{R}_{ > 0}$ adjusts the adaptation as a function of the acceleration and $[a, b]_m: = a$ if $a > b$, otherwise $[a, b]_m: = b$. Moreover, $\hat{\sigma}^2_{\dot{a}} > 0$ provides an estimation of the diagonal of $\hat{Q}^{\xi_F}_k$.

For significant velocity variations, the pseudoranges and delta-pseudoranges innovations are more noisy. We define the covariance matrix of measurement noises, adapted as a function of the acceleration computed in (4.23):

being $\alpha_{\rho}, \alpha_{\dot{\rho}}\in{\mathbf{R}_{ > 0}}$ adaptation sensitivity parameters, while $\hat{\sigma}_\rho^2 > 0$ and $\hat{\sigma}_{\dot{\rho}}^2 > 0$ provides an initial estimation of the diagonal components of $\hat{R}^F_k$.

5.   Tightly-coupled integrated navigation

Nomenclature:

In this section we describe the EKF-TC using pseudorange and delta-pseudorange as observables, and propose a method to compensate communication delays on the receiver output information.

5.1.   EKF-TC augmented state with receiver clock dynamics

The measurements (2.10) and (2.16) will be fused with the INS to obtain the desired kinematic state (attitude, position and velocity) and the user clock bias $b_{t}$ and drift $d_{t}$; all these variables are included on the augmented state:

The derivative of $\underline{\chi}^\rho$ is computed using (7.3) and (2.1):

with:

and the continuous-time vector of process noises $\underline{\xi}^\rho : = [\underline{\xi}^T\, \underline{\epsilon}_c^T]^T$, zero-mean gaussian white noise with autocovariance $\mathbf{E}\left[ \underline{\xi}^\rho(t) \underline{\xi}^{\rho\, T}(\tau) \right] = \delta(\tau - t)Q^\rho(t)$ where:

Let $\underline{\hat{\chi}}^\rho$ be an estimate of $\underline{\chi}^\rho$, hence we will consider the linearization:

5.2.   Observable measurement model and innovations

Lets group (2.12) and (2.16) on a single vector function $\underline{h}_k^{\rho}$ for the observables from $N_k^{\rho}$ satellites:

where $\underline{\eta}_k^{\rho}$ is a zero-mean measurement noise with covariance $R_k^{\rho} : = {\bf E}[\underline{\eta}_k^{\rho} \underline{\eta}_k^{\rho\, T}]$. Also, the Jacobian matrix $H_k^{\rho} : = \frac{\partial \underline{h}_k^{\rho}}{\partial \underline{\chi}^\rho}$ is defined with the submatrices $\frac{\partial \underline{\rho}_k^{c}}{\partial \underline{\chi}_k^\rho}$ and $\frac{\partial \underline{\dot \rho}_k}{\partial \underline{\chi}^\rho}$, with $N_k^{\rho}$ rows each:

Notice that with $N_k^{\rho}$ variables at a given sample instant $t_k$, $H_k^{\rho}$ may have a different number of rows.

Each row in $\frac{\partial \underline{\rho}_k^{c f}}{\partial \underline{\chi}^\rho}$ uses a gradient $\underline{\nabla}_{\chi} \rho_{s}^c$ defined for a given measured satellite $s$:

and each row in $\frac{\partial \underline{\dot \rho}_k^{f}}{\partial \underline{\chi}^\rho}$ uses the gradient $\underline{\nabla}_{\chi} {\dot \rho}_{s}$ of the delta-pseudorange associated to the same satellite $s$:

Hence:

where $\underline{0}_3 : = [0\, 0\, 0]^T$, and we have used the notation $\underline{\nabla}_{p}{\rho}_{s, k}^c : = \underline{\nabla}_{p}{\rho}_{s}^c\left({\underline{\chi}^\rho_k}, \underline{\tilde E}_k^\rho \right)$ for all the gradients.

5.3.   EKF-TC for INS/GNSS navigation

We consider the following definitions to implement the EKF steps as in the Appendix: $\underline{\hat \chi}_{k+1}^{\rho k}$ and $\underline{\hat \chi}_{k+1}^{\rho (k+1)}$ are the a priori and a posteriori estimates of $\underline{\chi}_{k+1}^\rho$, and for the $k+1$ step the prediction is $\underline{\hat \chi}_{k+1}^{\rho k} : = \left[ \underline{\hat \chi}_{k+1}^{k\, T}\, \, \underline{\hat x}_{c, k+1}^{k\, T} \right]^T$ where $\underline{\hat \chi}_{k+1}^{k}$ is estimated with (7.11) and $\underline{\hat x}_{c, k+1}^k$ is estimated as follows:

The covariance matrix of the a priori estimation error $P_{k+1}^{\rho k}$ is computed as in (7.12) while the a posteriori estimation error covariance matrix $P_{k+1}^{\rho (k+1)}$ is computed as in (7.15). The filter state estimation error is $\delta \underline{\chi}_{k}^{\rho}$, the process noise covariance is $Q^\rho_k$, and the augmented state transition matrix is here:

On the filtering stage at step $k+1$ we solve for $K_{k+1}^\rho$ as in (7.14) and the a posteriori estimation $\underline{\hat \chi}^{\rho (k+1)}_{k+1}$ as in (7.13). In order to compute the innovation, we also define the estimate $\underline{\hat y}_{k+1}^{\rho} : = \underline{h}_{k+1}^{\rho}\left(\underline{\hat \chi}^{\rho k}_{k+1}, \underline{\tilde E}_{k+1}^\rho \right)$.

5.4.   GNSS Receiver processing/communication delay compensation

Nomenclature:

A compensation of the delay in the communication of the GNSS receiver output is proposed here for the tightly coupled integration with the pseudoranges and delta-pseudoranges, assuming a known (or estimated) delay. The filter applies the update on time $t_k$ with measurements corresponding to a previous time $t_{pps, k}$; the difference is a delay mainly due to processing and communication tasks defined as $\Delta t_k$. Estimation ${\tilde \Delta} t_k$ of this delay is provided by the receiver output, which also contains the pseudorange and delta-pseudorange measurements and the ephemeris of the valid satellites.

Following ^[18] we consider a small enough^* $\Delta t_k = t_k-t_{pps, k}$, hence the output Jacobian matrix evaluated at $t_k$ can be approximated by the output Jacobian matrix evaluated at $t_{pps, k}$:

^*Validated for delays around $300$ ms or smaller.

and the following approximations are also valid:

where (5.15) is based on the approximation (5.13) and the smoothness of (7.14). On (5.16) an approximation of the correction at $t_k$ is obtained using the past innovation multiplied by the transition matrix $\Phi(\Delta t_k)$. As $\Delta t_k$ is small, (5.16) can be further simplified with $\Phi(\Delta t_k) \approx I$, leading to (5.17).

Using (5.13) and the smoothness of (7.14), the following approximation will be useful for the computation of the covariance $P_k^k$ in (7.15):

Finally, the remaining EKF steps are executed for prediction and filtering as defined in the Appendix by updating the state with the correction (5.17), and (5.18) for the computation of $P_k^k$. Figure 6 shows a possible algorithmic implementation:

Algorithm 5.1. Compensation of GNSS receiver output delays, for EKF-TC:

$-$ step 01: On $t_k > t_{pps, k}$ read the GNSS receiver output corresponding to the acquisitions made on time $t_{pps, k}$. We generate the following vector having only the valid observables and ephemeris:

where ${\tilde \Delta} t_k$ is the estimated delay, $\underline{\tilde{y}}^{\rho}_{pps, k} : = [\underline{\tilde{\rho}}^{cf}(t_{pps, k})^T\; \underline{\tilde{\dot{\rho}}}(t_{pps, k})^T]^T$ contains the valid observable measurements taken at $t_{pps, k}$, and $\underline{\tilde{E}}_{pps, k}$ contains the associated satellite's ephemeris.

$-$ step 02: Compute $\underline{\hat{p}}^e_{pps, k}$ by interpolation of the INS past information parameterized by the estimated delay; i.e.: $\underline{\hat{p}}^e_{pps, k} : = \underline{\hat{p}}^e_k(\tilde{\Delta} t_k, \mbox{INS}|1, ..., n)$.

$-$ step 03: Compute in the same way: $\underline{\hat{v}}^e_{pps, k} : = \underline{\hat{v}}^e_k(\tilde{\Delta} t_k, \mbox{INS}|1, ..., n)$.

$-$ step 04: Compute in the same way: $\underline{\hat{\theta}}^b_{pps, k} : = \underline{\hat{\theta}}^b_k(\tilde{\Delta} t_k, \mbox{INS}|1, ..., n)$.

$-$ step 05: Compute the estimate $\underline{\hat{y}}_{pps, k}^{\rho}$ of $\underline{y}_{pps, k}^{\rho}$ based on (2.12) and (2.16), using $\underline{\tilde{E}}_{pps, k}$, $\underline{\hat{p}}^e_{pps, k}$, $\underline{\hat{v}}^e_{pps, k}$ and $\underline{\hat{\theta}}^b_{pps, k}$.

$-$ step 06: Compute $H_{pps, k}^{\rho}$.

$-$ step 07: Assign $H_{k}^{\rho} \leftarrow H_{pps, k}^{\rho}$ (approximation (5.13)).

$-$ step 08: Execute the EKF-TC to compute the gain $K_k^{pps} = K_k^{\rho}(H_{pps, k}^{\rho})$ and approximate (5.17) and (5.18).

$-$ step 09: Assign $(\underline{\chi}^{\rho})^{k}_k \leftarrow (\underline{\chi}^{\rho})^{k-1}_k + K_k^{pps} \cdot (\underline{\tilde{y}}^{\rho}_{pps, k} - \underline{\hat{y}}_{pps, k}^{\rho})$ (approximation (5.17)).

$-$ step 10: Assign $P_k^k \leftarrow \left(I - K_k^{pps} \cdot H_{pps, k}^{\rho} \right) P_k^{k-1}$ (approximation (5.18)).

6.   Simulation results

The proposed strategies are verified by numerical simulations of a small lift two-stage satellite launch vehicle direct ascent to a Low Earth Orbit. The trajectory is shown in Figure 8, which has typical values with an acceleration profile up to $6\; {\rm g}$, and an injection altitude of $350\; {\rm km}$. For all these simulations, the vehicle is on the launch PAD during the first minute, while the injection velocity is achieved after 10 minutes of simulation.

The acceleration profile shows the take-off after the first minute. Before four minutes of simulation, the acceleration shows the first stage Main Engine Cut-Off event (which could have a cluster of thrusters ^{[34,35,36,37]}), followed by the ignition of the second stage engine. The second stage Main Engine Cut-Off event happens after 10 minutes of simulation time. The trajectory is shown in Figure 7, for a hypothetic launch scenario towards a middle inclination Low Earth Orbit, showing in white the burn of first stage engine, in cyan the burn of second stage engine and in green the injected orbit.

Table 1 shows the main performance parameters considered for the inertial measurement unit (IMU). The initial value of each bias (at $t = t_0 = -60$ seconds) is considered as a small fraction ($\approx5\%$) of the true initial value: this considers an alignment process previous to the launch. Once the vehicle is launched at $t = 0$, and during the two-minute burn period of the engine, the measurements are modified to consider the variation in temperature in the following way:

where $\delta_T\underline{b}_\omega(t)$ and $\delta_T\underline{b}_f(t)$ are ramp functions active between $t = 0$ (lift-off) and $t = 120$ seconds (main engine cut-off). Therefore these ramps are characterized by their values at $t = 0$ and $t = 120$. For this particular simulation we took:

and these values are linearly interpolated for $0\leq t\leq 120$. Notice that we have assumed zero orthogonality and scale factor errors, which are not significant relative to the errors we have modelled, as typically considered for satellite launchers (see ^[38]).

6.1.   GNSS observables simulation during ascent

On ^[3] we developed a numerical simulation tool for GNSS constellations and a receiver, including a model of each observable's acquisition as a function of the link budget and the internal behaviour of the receiver. The purpose was to have the possibility to test navigation strategies on a controlled framework, which was used to evaluate loosely coupled and tightly coupled strategies as shown in the thesis work ^[18]. The behaviour of the receiver processing was modelled with a state-machine, while the acquired satellites were modelled using the link bugdet using the antenna radiation pattern and compared with the signal to noise threshold $C/N_0$. During the ascent, new satellites appear in visibility and other satellite are lost moving towards smaller elevations. Finally, the receiver processing delay depends on the number of visible satellites, which can also be modelled. These models were later validated using a Spirent GSS8000 simulator and a multi-antenna GNSS receiver. There are four antennas assumed to be located all at the same point of the IMU location (to simplify the exposition) and pointing orthogonal to the launch vehicle longitudinal direction, equi-spaced by $90^\circ$, as shown in Figure 9. The vectors $\underline{n}_i$ are the normal of each antenna (left side in Figure 9), while a simple radiation pattern profile $G_R(\theta_{s}^a)$ (red curve shown on the right side) has been simulated as proposed in ^[18]. This simplified pattern is also applied relative to the axial direction, which means that when the vehicle is vertically launched, the signal from a satellite located exactly on the zenith will be attenuated. Finally, we have considered all the antennas located on the Launch Vehicle (LV) longitudinal axis, although they are actually located on the external surface of the fuselage: as the attitude of this vehicle is three-axes stabilized, the effect of each antenna's arm is considered negligible (see ^[18] for a more complete model).

6.2.   Loosely coupled INS/GNSS results

Here we evaluate the EKF-LC integrated navigation, using the fix obtained with the three possible inputs: ⅰ) simulated kinematics with additive noise, ⅱ) PVT solution using the Bancroft algorithm, and ⅲ) a filtered solution computed using a Kalman filter in the GNSS receiver. All these solutions are evaluated with a numerical simulator considering several effects as communication delay, GNSS data outages and antenna failures, which will be later compared with the corresponding tightly coupled solution.

6.2.1.   GNSS measurement delay compensation.

The Figure 10 shows the result with constant delay $\Delta t_k = \tilde{\Delta} t_k = 300\; \mathrm{ms}$ comparing the cases without compensation (red), with the delay compensation given in the Appendix (blue) and without delay (black, as reference). The GNSS receiver used for hardware in the loop testing ^[32] has a smaller delay than this bound, taken as a worst case.

6.2.2.   Point vs. filtered PVT solutions for the simulated GNSS receiver

For a given simulated measurements serie with pseudoranges and delta-pseudoranges, we compare the filtered solution with the point PVT solution given by the Bancroft algorithm. The deterministic atmospheric delays are not simulated or compensated, while the clock bias and drift are simulated and compensated by the filter. Moreover, we consider a zero-mean white noise with covariance $4\; {\rm m}^2$) for the pseudorange and a covariance $0.01\; {\rm m/s}^2$ for the delta-pseudorange. In Figure 11 the clock bias and drift are shown in red and black respectively, while the respective values computed by Bancroft algorithm are shown in light blue and light green, and the values computed by EKF-GNSS on the receiver are shown in blue and green. The filtered solution converges during the first seconds and is kept with a mean error close to zero. The point solution using Bancroft algorithm has more estimation error deviation. Figure 12 illustrates the autocorrelations of observables innovation sequences of some satellites during all the simulation, which shows the expected uncorrelation.

Lets compare the position and velocity of point and filtered solutions: the Figure 13 shows the position and velocity error of the point solution (light blue and light green) and the filtered solution (blue and green) respect to the simulated values; we have shown only one coordinate to simplify the graph and the upper trace shows the number of satellites involved in the solution: $N_k$ and $N_k^F$ respectively for the point and filtered ones. The dotted lines shows the estimation deviation of the EKF-GNSS filter, where it can be observed the effect of the adaptation on matrices $\hat{Q}_k^{\xi_F}$ and $\hat{R}_k^{F}$ computed with the acceleration $\underline{\hat{a}}_k^P$, as for higher acceleration there will be more velocity and position deviations.

The filtered solution also provides the covariances of the position and velocity from the EKF-GNSS covariance matrix $P_k^{F, k}$. This makes possible to use it for the external measurement noise description used by the loosely coupled integrated navigation system.

6.2.3.   Simulation of satellite tracking failures

The number of tracked satellites on a GNSS receiver is a very important indicator of the quality of the solution, and also may be the cause of the unfeasibility to compute a navigation solution at all. Here we simulate restrictions which might not be realistic as real causes of failures, but determine a behaviour with low number of satellites which may be observed during flight ^[3].

In particular, on a multi-antenna receiver with four patch antennas we consider the failure of three of them during a 30 seconds segment shown in grey in Figure 14. After antennas recovery, the visible satellite are incorporated one by one after a processing in sequence, which is modelled by a state-machine. On this simulation the point solution is not available even one minute after the failure, However, as the filter propagates there is a continuity on the availability of the filtered solution, with increasing uncertainty shown in the filter estimation covariance $P_k^{F}$, in dotted lines.

6.2.4.   Loosely coupled integrated INS/GNSS navigation with different GNSS positioning solutions

Following the statement of the problem made in Figure 2 lets compare each of the possible loosely coupled INS/GNSS implementations, where we will not consider solution availability delays to simplify the exposition. The results are summarized on Tables 2 and 3 with RMS (root-mean-square) values of the estimation errors for the position and velocity using the INS/GNSS with the different PVT sources:

PV Bancroft: Point solution using Bancrof algorithm. The matrix $\hat{R}_k$ taken by the INS/GNSS is constant and was estimated using the standad deviation of a set of outputs as shown in Figure 13.

EKF-GNSS $\mathbf{\#1}$: The Extended Kalman Filter for the receiver estimates the receiver position, velocity, clock bias and clock drift, implemented using adaptive covariances for the process noise and measurement noise. In order to initialize these matices, we use: $\hat{\sigma}_{\dot{d}_c} = 0.005\; \mathrm{m/s^2}$, $\hat{\sigma}_{\dot a} = 0.5\; \mathrm{m/s^3}$, $\hat{\sigma}_{\rho} = 2\; \mathrm{m}$ and $\hat{\sigma}_{\dot{\rho}} = 0.1\; \mathrm{m/s}$, $\alpha_Q = 0.01$, $\alpha_{\rho} = 1\; \mathrm{m s^2}$ and $\alpha_{\dot \rho} = 0.01\; \mathrm{m}$. The INS/GNSS filter uses a covariance matrix $\hat{R}_k$ obtained from the elements of the filter state covariance $P_k^{F}$ of the receiver and scaled by 10.

EKF-GNSS $\mathbf{\#2}$: Idem $\mathbf{\#1}$ but without adaptation on covariance matrices $\hat{Q}_k^{\xi_F}$ and $\hat{R}_k^F$ of the EKF-GNSS receiver filter. Lets define $\hat{\sigma}_{\dot{d}_c} = 0.005\; \mathrm{m/s^2}$, $\hat{\sigma}_{\dot a} = 1\; \mathrm{m/s^3}$, $\hat{\sigma}_{\rho} = 4\; \mathrm{m}$ and $\hat{\sigma}_{\dot{\rho}} = 0.2\; \mathrm{m/s}$. The INS/GNSS uses the same $\hat{R}_k$ matrix values defined previously for PV Bancroft.

EKF-GNSS $\mathbf{\#3}$: Idem $\mathbf{\#1}$, but the filter state vector of the receiver only contains the position and velocity; i.e., without including the acceleration^†.

^†This implies to discard the first three rows of $B^F_{pva}$ en (4.8) and using higher variances on $Q^{\xi_F}_k$ diagonal elements (e.g.: 20 times) respect to those in $\mathbf{\#1}$.

PV simulation + noise: Position and velocity receiver solutions defined as the true values plus a zero-mean white noise with small standard deviations given by $10\; \mathrm{cm}$ in position and $0.01\; \mathrm{m/s}$ in velocity. These values were also used to define $\hat{R}_k$ on the INS/GNSS algorithm.

With the exception of case EKF-GNSS $\mathbf{\#3}$, the innovations are decorrelated. Figure 15 compares the estimation errors of the Extended Kalman Filter on the INS/GNSS using the receiver measurements obtained on cases PV Bancroft (blue), EKF-GNSS $\mathbf{\#1}$ (red) and PV Simulation + noise (grey). The point and filtered solutions have similar results with different covariances $P_k$.

Figure 16 shows the same cases under a failure on three antenna patches, leaving only one antenna operative during $30$ seconds. There is an advantage on using filtered solutions which are also available during the failure, as the filter propagates previous states, although with increasing uncertainties.

6.3.   Tighty coupled INS/GNSS results

6.3.1.   Tropospheric delay compensation results

The Tropospheric Delay Compensation (TDC) proposed in subsection 3.1 is evaluated here through numerical simulation of a satellite launcher ascent. This model extends the typical range of elevations to consider even negative values, which are not available at sea level due to the need to consider a minimum positive elevation in order to avoid obscuration and multi-path effects. The minimum elevation configured at sea level was $10^\circ$ and for higher altitudes the minimum elevation is modified with the following function:

where $h_\tau : = 7518.8\; \mathrm{m}$ as before, $h_{v, k}$ is the vehicle's altitude at time $t_k$ and $R_e$ is the local Earth radius. Notice that this function is smooth and provides a value of $10^\circ$ at sea level and a minimum angle which can be negative at typical LEO altitudes.

The TDC model was applied to compensate the simulated delays obtained with a more accurate model developed in ^[3]; the results are shown in Figure 17 for a typical satellite launch ascent trajectory where the orbit injection is achieved 9 minutes after a vertical launch from sea level. The upper graph shows the evaluation of the simulated delays, and the graph below shows the error relative from the TDC model to the more accurate model in ^[3].

The error is bounded by $\pm 1\; \mathrm{m}$, which is considered a good match considering the simplification of the TDC model proposed here in comparison with the more complex model developed in ^[3] for simulation purposes.

6.3.2.   Comparison between tightly and loosely coupled integrated navigation solutions

Here we evaluate the proposed tightly coupled (TC) integrated navigation against the loosely coupled (LC) integrated navigation proposed in subsection 7 and subsection 4.1), under antenna failures as simulated previously in subsection 6.2.3.

The TC vs. LC results are shown in Figure 18 (nominal scenario) and Figure 19 (GNSS antenna failure with recovery scenario). In the nominal scenario, both estimation error profiles are good during the ascent. However, under antenna failure the small number of tracked satellites makes unfeasible to compute the Bancroft solution (i.e. the fix) and therefore the LC algorithm estimation error is increased significantly. On the other hand, the TC algorithm output continues the fusion between the INS with the GNSS observables, which is clearly seen in Figure 19.

In section 6.2.4 we have shown the results with loosely coupled navigation and different positioning sources, including the same Bancroft algorithm taken as a reference here. It can be seen that the results with tightly coupled navigation are at least similar or better compared to any of these solutions, including the receiver filtered output which includes a propagation of the solution when there is no fix. Moreover, the complexity here is focused on the integrated navigation algorithm, which has more information due to the fact that includes the INS data which is not shared with the GNSS receiver.

The Figure 21 shows the Sky Plot identifying the declination and azimuth of the satellites in view by the GNSS receiver, at the moment of lift-off, just after the failure of three antennas (only one antenna remains operative, hence the satellites in view are limited to an angular sector), and at orbit injection when all the antennas are recovered and also due to higher altitudes there are more potential satellites to track. Notice that a satellite's declination greater than $90^\circ$ appears only when the vehicle ascents, otherwise at lift-off the maximum declination is $90^\circ$ (i.e. on the horizon).

6.3.3.   GNSS observables delay compensation for tightly coupled integrated navigation

The tightly coupled (EKF-TC) integration with delay compensation (Algorithm 5.1) is evaluated with the same direct ascent trajectory considered in the Appendix. The results for a known delay $\Delta t_k = 300\; \mathrm{ms}$ are shown in Figure 22. The performance is compared with respect to the loosely coupled (EKF-LC) integrated navigation when there are at least four satellites in view, leading to similar results.

This method based on the approximation of (5.13) which was validated by simulations in ^[18].

6.3.4.   Robustness against IMU bias temperature sensitivity during ascent

During the ascent trajectory, the temperature has a significant variation on each sector, including the payload and the secondary structures where the navigation units are mounted. In particular, the gyroscopes and accelerometers have certain sensitivity with respect to the temperature time derivative ^[42]. Considering an estimate of the temperature variation, we have considered a bias ramp during two minutes after lift-off, with is added to a constant bias component considered separately since the first minute of simulation when the vehicle is on the launch platform.

With this common statement of the problem, we test two EKF-TC implementations: one with IMU bias estimation and the other without it, as considered for $\underline{\chi}_k^\rho$ respectively in (7.5) and (5.1). The Table 4 shows the position, velocity, attitude and angular velocity errors as estimated in mean square just after the upper stage main engine cut-off. The errors obtained when the IMU bias compensation is activated are also good enough to compute the first estimation of the mean orbital elements ^[43]. However, when this compensation is not active, the velocity error becomes two order of magnitude higher.

Figure 23 shows the evolution of these estimation errors during all the ascent (in black, without bias estimation, in red, green and blue the result using bias estimation). The discontinuities after the first stage burn-out, staging and second stage engine start are shown useful to improve the IMU bias estimation, increasing navigation sensibility ^[44]. On the other hand, when the IMU biases are not estimated, the innovations becomes correlated, as shown in Figure 24.

6.3.5.   Preliminary validation with hardware-in-the-loop

A first validation of the proposal using hardware was performed using a four-antenna GPS/GLONASS L1 receiver being developed ^[32], with a Spirent GSSS 8000 simulator. This allowed to assess the position, velocity and attitude integrated navigation error and the evolution of number of satellites in comparison with the numerical simulations taken as a common simulation framework for all previous simulation results shown in this work. The number of satellites used for this simulation varies between $6$ and $8$ during the ascent, which is a lower number limited by the current version of the receiver which considers an elevation mask of $10^\circ$. Figure 25 can be directly compared with the numerical results in Figure 22, showing a similar behaviour in position and velocity errors, with a higher stationary error which may be in part due to a worse visible satellite number and distribution. Notice that with the GNSS receiver hardware there is an error in position around $10\; \mathrm{m}$ in Figure 25 which does not vanish with increasing altitudes, showing that the approximation (2.12) is more uncertain respect to the models implemented on the GNSS simulator, mainly for the uncertain ephemeris of GNSS satellite orbits, the ionospheric model (which was compensated with zero alfa and beta coefficients) and relativistic model used for compensation. Regarding the attitude navigation errors, a closer analysis shows that this error is mainly explained by the attitude error in the roll axis, which is a known effect for the typical trajectory of a launch vehicle (see ^[38]). Table 5 summarizes the worst case results for each simulation type (SIL: software in the loop, i.e. numerical simulation; HIL: hardware in the loop, i.e. with GNSS RF simulator and a GNSS receiver) over all the simulation time.

7.   Conclusions

In this work we consider design trade-offs for the navigation function on a satellite launcher. A comparison has been made between different INS/GNSS couplings. The evaluation is made with numerical simulations for a small lift launch vehicle with vertical launch and a final injection into LEO orbit. The robustness of the proposal has been tested under the following scenarios:

● Low number of satellites, which may be caused by antenna failure.

● Delay in the availability of the GNSS receiver output, which may be caused by communication and/or processing delays.

● Inertial sensor bias variations, to model temperature sensitivity.

The main contributions presented in this work are:

● Compensation of the delay on the GNSS receiver output, for both EKF-LC and EKF-TC.

● Compensation of tropospheric delay in the GNSS signals, extending a simple deterministic model to allow negative elevations of the line of sight.

● Adaptation of the covariance matrices to improve the EKF-LC and EKF-GNSS results, under the wide acceleration ranges of the ascent phase.

The best solution was given by the EKF-TC, which allows to process the GNSS observables given by the receiver even when there is no PVT solution output on the receiver. It was also verified that inertial sensor bias variation can be successfully estimated. A first hardware in the loop test was made using a Spirent GSSS 8000 simulator with a multi-antenna GNSS receiver under development, and will be continued as future work.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work has been supported by the National Commission of Space Activities in Argentina (CONAE), and the Ministry of Science, Technology and Innovation in Argentina. The authors would like to thank Prof. Juan Giribet for his suggestions and criticism.

Conflict of interest

The authors declare that there are no conflict of interest.

A.   Appendix

This Appendix is based on ^[16] and is provided as a self-contained summary of the loosely coupled integrated navigation. The state vector $\underline{x}(t)$ contains as sub-vectors a parametrization of the attitude, velocity and position of a body frame ${\mathcal B}$ of the vehicle relative to the Earth Centered Earth Fixed (ECEF) frame ${\mathcal E}$:

where $\underline{\theta}_{be}^b\in \mathbf{R}^3$ represents the vector angle (i.e. the rotation axis versor times the angle of the rotation) of the attitude of the body ${\mathcal B}$ respect to ${\mathcal E}$, expressed in ${\mathcal B}$, $\underline{v}^e \in \mathbf{R}^3$ is the translational velocity and $\underline{p}^e\in \mathbf{R}^3$ the position on ${\mathcal E}$ frame.

The vector of inertial measurements is $\underline{\mu}: = [\underline{\mu}_\omega^T\ \underline{\mu}_{\rm f}^T]^T$, where $\underline{\mu}_\omega$ and $\underline{\mu}_{\rm f}$ are respectively the outputs of a set of three gyroscopes and three accelerometers mounted following the axes of a right-hand frame. Additive noises $\underline{\xi}_\omega$ and $\underline{\xi}_{\rm f}$ are defined for gyros and accelerometers respectively.

The acceleration in body frame $\underline{a}^b$ can be estimated with a gravity model and the specific force measurements: $\underline{\gamma}^b+\underline{\xi}_g$ is the true gravity seen in body frame, where $\underline{\xi}_g$ is a model uncertainty and $\underline{\gamma}^b$ is obtained from a gravity model $\underline{\gamma}^e$ in ${\mathcal E}$ frame, the vehicle's position $\underline{p}^e$ and the attitude matrix ${\bf C}_b^e = {\bf C}_b^e{\left(\underline{\theta}_{be}^b\right)}$ as $\underline{\gamma}^b = {\bf C}_b^{e\, T}\underline{\gamma}^e(\underline{p}^e)$.

Let $\underline{\xi}_\mu: = [\underline{\xi}^T_\omega\ \underline{\xi}^T_{\rm f}\ \underline{\xi}^T_g\ \underline{\xi}^T_p]^T$ where $\underline{\xi}_p$ is a random term considered for the position derivative $\underline{\dot p}^e : = \underline{v}^e + \underline{\xi}_p$ which purpose will be later explained. On the other hand, let $\underline{b} : = [\underline{b}_\omega^T\ \underline{b}_{\rm f}^T]^T\in\mathbf{R}^{6}$ be the bias of angular velocity measurements $\underline{b}_\omega$ and specific forces bias $\underline{b}_{\rm f}$ which are modelled as Markov processes, integrating random variables defined by $\underline{\xi}_b: = [\underline{\xi}^T_{b\omega}\ \underline{\xi}^T_{b{\rm f}}]^T$, where $\underline{\xi}_{b\omega}$ and $\underline{\xi}_{b{\rm f}}$ are uncorrelated gaussian processes with zero mean.

The following relations apply:

The vehicle's kinematics can be obtained as a function of the inertial measurements (see ^[11]):

where we defined the zero matrix $\mathbf{0}_{n \times m}\in \mathbf{R}^{n\times m}$, the identity matrix $\mathbf{I}_{n\times n}\in \mathbf{R}^{n\times n}$, the zero vector $\underline{\mathbf{0}}_n \in \mathbf{R}^n$, $f_\mu$ gives the kinematic coordinates derivatives as a function of the measurements and model parameters, while $B_\mu$ is a linear transformation between the measurement noises and the effects on kinematic coordinates derivatives.

The vector of the process noises is $\underline{\xi}: = [\underline{\xi}^T_\mu\ \underline{\xi}^T_b]^T$, it has zero mean and autocovariance:

where $\delta(\cdot)$ is the Dirac function ^[27] and $0 < Q^\xi(t)\in\mathbf{R}^{18\times 18}$ for all $t$. Let $\underline{\chi}: = [\underline{x}^T\ \underline{b}^T]^T\in \mathbf{R}^{15}$, be the augmented state vector, hence the differential equation in $\underline{\chi}$ can be written as in ^[11] as:

where:

with $\mathbf{0}_3 : = \mathbf{0}_{3 \times 3}$ and $\mathbf{I}_3 : = \mathbf{\bf I}_{3 \times 3}$. Notice that the random term $\underline{\xi}_p$ enforced a full rank $B(\underline{\chi})\in\mathbf{R}^{15\times 18}$.

A.1.   Extended Kalman Filter

Let $\underline{\hat \chi}(t)$ be the estimation of the augmented state. We can linearize the associated dynamic model and describe the evolution of the error around the estimation, given by $\delta \underline{\chi}: = \underline{\chi}-\underline{\hat \chi}$. The resulting system is linear time varying:

with ^[27]:

where ${\bf S}(\underline{u})$ is the matrix associated to the vector product $\underline{u}\times$, the vector $\underline{\rm {\hat f}}^b : = \underline{\mu}_{\rm f} - \hat{\underline{b}}_{\rm f}$ is an estimate of the specific force as a function of the acceleration measured in body frame and the estimate $\hat{\underline{b}}_{\rm f}$ of the measurement bias, the vector $\underline{\bf {\hat \omega}}_{bi}^b = \underline{\mu}_\omega - \hat{\underline{b}}_\omega$ is an estimate of the angular velocity in body frame as measured by the gyroscopes and compensated by the gyro bias estimation $\hat{\underline{b}}_\omega$, the matrix ${\bf {\hat C}}_b^e$ is an estimate of the attitude matrix and the vector $\underline{\bf \Omega}_{ei}^e$ is the Earth rotation velocity.

Let $t_k$, with $k\in\mathbf{N}_0$, be the nominal discrete times on which the Kalman filter processes external measurements from the receiver, and $t_k: = t_0+k\Delta_K$. The augmented state error $\delta \underline{\chi}_k: = \delta \underline{\chi}(t_k)$ verifies:

where $\underline{\xi}^\chi_{k}\in \mathbf{R}^{15}$ is a discrete-time process zero-mean gaussian white-noise ^[28]. Finally, the transition matrix $\Phi(t_{k+1}, t_{k})$ for the state $\delta \underline{\chi}_k$ can be estimated as in ^[16], using a $5th$-order polynomial.

The discrete-time vector function $\underline{\tilde{y}}_k$ represents the external measurements from the GNSS receiver, with the position and velocity vectors in ECEF frame, evaluated at time $t_k$. This can be written in terms of the augmented vector as:

where $\underline{\chi}_k: = \underline{\chi}(t_k)$, $\underline{\eta}_k$ is a vector with measurement noises, defined as a discrete-time gaussian and uncorrelated stochastic process, with zero mean and covariance matrix $R_k : = {\bf E}[\underline{\eta}_k \underline{\eta}_k^T]$, and $H$ is given as follows:

When a new measurement $\underline{\tilde{y}}_{k+1}$ arrives on time $t_{k+1}$, the prediction and filtering tasks are executed, which are detailed in the following subsections.

A.1.1.   Prediction

The prediction step on $k+1$ is defined as ^[28]:

where $B_k: = B(t_k)$, $\Phi_k : = \Phi(t_{k+1}, t_{k})$, $Q^\chi_k$ is the covariance of the discretized noise $\underline{\xi}^\chi_{k}$^‡, $\underline{x}_{\mathrm{INS}}: = [\underline{\hat \theta}^{bT}_{be} \ \underline{\hat v}^{eT} \ \underline{\hat p}^{eT}]^T$ is the output of the INS algorithm (see ^[11,27]), and $P_{k+1}^k = {\bf E}[\delta \underline{\chi}_{k+1}^k \delta \underline{\chi}_{k+1}^{k\, T}]$. The initialization $k = 0$ is made with a priori information written in terms of $\underline{\hat{\chi}}_0^0 : = {\bf E}[\underline{\chi}_0]$ and $P_0^0$.

^‡We assume that it is pre-filtered to avoid alias and preserve de-correlation (see ^[27]).

A.1.2.   Filtering

The a posteriori estimation $\underline{\hat{\chi}}_{k+1}^{k+1}$ is:

where $\hat{\underline{y}}_{k+1}^k : = H\hat{\underline{\chi}}_{k+1}^k$ is the prediction of the external measurement while:

A.2.   EKF error bound condition

In ^[29] it is shown that the discrete-time Extended Kalman Filter estimation error can be analysed in a stochastic framework with a nonlinear observability condition, which determines an exponentially bounded behaviour:

Definition 7.1. The random process $\underline{\chi}_k - \underline{\hat{\chi}}_k$ is said to obey an exponential bound in mean square if there exist positive real numbers $a, b > 0$ and $0 < c < 1$ such that for all $k\geq 0$:

Among other sufficient conditions for this behaviour (see ^[29]), the following bounds are required:

Notice that $\underline{q} > 0$ requires $B_k = B(\underline{\chi}(t_k))$ full rank, which was enforced by adding the random term $\underline{\xi}_p$. Other ways to guarantee exponentially boundedness in mean square can be found in ^[29,30].

A.3.   Delayed GNSS measurements compensation

At time $t_{pps, k}$ the internal acquisition of the GNSS signals is made, but the integrated navigation block will have this information available later, after all the processing and communication latencies, which implies a random delay $\Delta t_k: = t_k - t_{pps, k}$. Following the compensation method developed in ^[16], we assume that the GNSS receiver output contains a measure $\underline{\tilde p}^{e\, T}_{pps, k}$ of the position, a measure $\underline{\tilde v}^{e\, T}_{pps, k}$ of the velocity and an estimate $\tilde{\Delta} t_k$ of the delay. This information corresponds to the time $t_{pps, k} < t_k$, while the EKF will apply the associated correction at time $t_k$. We integrate one step by the Euler method to estimate the external measurement for $t_k$:

To this end, we consider an interpolator which estimates the velocity $\hat{\underline{v}}^e_{pps, k}$ and acceleration $\hat{\underline{a}}^e_{pps, k}$ using past INS data for $t_k - \tilde{\Delta} t_k \approx t_{pps, k}$. The external measurements covariance $R_k$ now includes receiver measurement noises taken at $t_{pps, k}$ and the effect of the propagation error up to time $t_k$ (estimated with a vector function $\underline{\varepsilon}_k^*(\underline{\hat a}^e_{pps, k}, \underline{\hat a}^e_{pps, k-1}, \tilde{\Delta} t_k)$ and a delay estimation deviation $\hat{\sigma}_{\eta_{\Delta t}}$), estimated in ^[16] as:

where $\underline{\hat{\dot{y}}}_{pps, k} : = [\underline{\hat{v}}^{e\, T}_{pps, k}\; \underline{\hat{a}}^{e\, T}_{pps, k}]^T$ is the interpolator output; see ^[16] for a complete description, on which a covariance $\hat{R}_{k}^{pv}$ is assumed as provided by the receiver, while $\hat{R}_k^{\hat{v}\hat{a}}$ is determined by the interpolator. From a practical point of view, this method does not require to read the PPS signal.