State-Space Measurement Update for GNSS/INS Integrated Navigation

*is paper theoretically derives the equivalence conditions for the loosely and tightly coupled GNSS/INS integration algorithms. Firstly, the equivalence is proved when using single epoch GNSS measurements, which means the GNSS processor provides standalone solution. *en, the equivalence proof is further extended for the filtering solutions, which are usually applied for differential GNSS and precise point positioning. Based on these, different state and measurement models for GNSS/INS integration navigation are summarized, and natural differences among these models are discussed.*is indicates that once the same measurement and predict information are used, the integration would be equivalent no matter what kind of coupling schemes are used. A flight dataset with GNSS and tactical IMU data is used to evaluate the equivalence and discrepancies among four different measurement models, and the results confirm the theoretical derivations.


Introduction
Although differential global navigation satellite system (DGNSS) and precise point positioning (PPP) techniques nowadays provide viable precise positioning options, the major disadvantage of GNSS still remains: signal blockage due to obstructions in urban and built-up environments and extreme power attenuation of the signals when operated indoors [1]. Integrated navigation systems employing GNSS with other sensors, such as a self-contained inertial navigation system (INS), have the potential to provide high levels of accuracy and integrity [2,3]. e integrated systems can not only mitigate the weakness of GNSS but also bound the INS error that otherwise would grow with time when the INS operates alone [4]. e presence of multiple data sources theoretically provides functional redundancy as well as greater observability of the desired navigation states [5,6].
To achieve precise navigation, DGNSS or PPP are involved in the loosely coupled (LC) or tightly coupled (TC) GNSS/INS integration [4,[7][8][9][10][11]. In the LC approach, the GNSS solution is computed firstly using only the GNSS measurements and then combined with the IMU data in a secondary processor [2,4]. When less than 4 satellites are visible, a GNSS solution is not possible and the LC integration processor will be processed in a free-inertial model without GNSS position or velocity updates to constrain the inertial drift [12,13]. e TC approach directly combines the GNSS pseudorange, pseudorange-rate, and carrier phase measurements with the inertial measurements in an extended Kalman filter [14]. erefore, the data utilization is better than the LC approach, and carrier phase ambiguities could be reinitialized for short GPS outages much more quickly. It has been agreed that the system reliability and faulttolerance are improved with the TC approach [10,15,16]. In recent decades, tightly coupled DGNSS/INS or PPP/INS integration has become popular, because it has the advantage of providing accurate position information even when GNSS measurements are not enough for standalone processing and is theoretically optimal in a filtering sense, especially in urban navigation applications [8,16]. PPP using ionosphere-free combination or ionosphere-constraint combination has been introduced into tightly coupled integration with INS [15,17]. Constraints such as context features, GNSS baseline length, vehicle moving properties have also been considered to improve the observability, reliability, stability, and availability of the GNSS/INS integration systems [1,2,4,9,18]. However, the heavier computational burden and more data exchange among subsystems usually complicate the TC processing and reduce the efficiency.
In this paper, we proposed a unified integration mode, with the measurement update based on the state space, by which the LC and TC processors are united and both the system efficiency and precision can be guaranteed. e paper is organised as follows. e general dynamic and measurement model of GNSS/INS integration algorithms is described in Section 2. en, in Section 3, we theoretically construct the state-based measurement update model based on the equivalence derivation between the LC and TC approach involving single epoch GNSS measurements.
e state-based measurement update approach is further extended when time correlations are considered in the GNSS processor in Section 4. Based on this, the possible equivalent models for GNSS/INS integration are summarized. Experiment results are presented to demonstrate the theoretical derivation in Section 5. Finally, some concluding remarks are presented in Section 6.

Full State Model Including All Parameters.
For GNSS/ INS integration processing, the following states can be consisted in the integration state vector as given in Table 1.
For the GNSS positioning model, the general full state vector is and for the GNSS/INS integration model, the general full state vector is where variables in (1) and (2) are explained in Table 1.

Dynamic and Measurement Model for Kalman Filter.
e system error dynamics of GNSS/INS integration are based on the INS error equations [19]. e psi-angle error equations of INS are expressed in the navigation frame (nframe, local east, north, and up) for position, velocity, and attitude updating [20][21][22].
where the superscript (n) refers to the navigation frame (nframe), while the subscripts (i) and (e) denote the inertial and the earth-centred earth-fixed frame (i-frame and eframe), respectively; x n p , x n v , and ψ are, respectively, the position, velocity, and attitude error vector; ω n ie is the earth rotation angular rate; ω n in is the angular rate of the navigation coordinate system with respect to the i-frame; f n is the specific force vector; δg n is the error of the gravity vector in the n-frame; ω n en is the angular rate from the e-frame to the nframe; C n b is the coordinate transformation matrix from body frame (b-frame) to n-frame; and x ∇ and x ε are the accelerometer and gyro drift vector expressed in the b-frame, respectively.
e Kalman filter (KF) algorithm uses the measurement vector, measurement model, and dynamic model to maintain optimal estimates of the state vector. e general KF dynamic model is given as where w k is the system noise vector comprising a number of independent random noise sources and F k is the transition matrix derived from the known properties of the system. With the dynamic model describing the properties of the system's dynamic, the KF measurement model describes the functional relationship between the measurement vector and the state vector and is generally given as where subscript k is the epoch indicator; z k ∈ R m is the measurement vector; H k ∈ R m×f is the design matrix relating the measurements to the state parameters; and ε k is the measurement noise vector with a covariance matrix R k . e KF measurement update is triggered at every GNSS measurement epoch using the difference between the GNSS and the INS mechanization solution. In LC mode, the states estimated from the GNSS processor are used as measurements in KF. In TC mode, the GNSS pseudorange, pseudorange-rate, and carrier phase measurements are directly used as measurements in KF. e measurements for LC integration are the position and velocity difference between GPS and INS centre and are given as where superscript * means the state-resolved frame and variants are explained in Table 1. e solution of the above system, (4) and (5), can be expressed by prediction and measurement update (e.g., see [5]). e prediction of the state and its covariance matrix are 2 Mathematical Problems in Engineering where x k is the state prediction; x k−1 is the state estimate at the previous epoch; Φ k ≈ I + FΔt + (F(Δt) 2 /2!)+ is the state transition matrix; P k−1 is the covariance matrix of x k−1 ; Q k is the dynamic noise covariance matrix; and P k is the covariance matrix of the prediction x k . When the measurements are available, the state vector is updated. e Kalman gain K k , the updated state, and its covariance matrix are given by where x k is the KF estimate of the state and P k is the corresponding estimation covariance matrix.

State-Based Measurement Model for GNSS/ INS Integration
With the utilization of multiconstellation GNSS, integration measurements and states would increase dramatically and accordingly increase the date exchange and computation burden. To simplify the integration structure, we only consider the INS-related states in the integration model as (1)- (9) show that the common states between GNSS and INS include the position and velocity states, e relationship between x vp,G and x avp are [23,24] x where C * b is the transition matrix from b-frame to * -frame; ω b eb is the angular rate of the b-frame respect to the e-frame, resolved in the b-frame; and superscript * means the state-resolved frame and other variants are expressed in Table 1.
e matrix format is with Usually, the velocity state x v is solved in n-frame. If x p and x p,G are also given in n-frame, there is 3.1. Measurement Model at Single Epoch. As the states given in (1), the measurement model of GNSS positioning can be written as where l ∈ R m is the measurement vector; A ∈ R m×6 , B ∈ R m×n b , and C ∈ R m×n c are the design matrices related to the respective states; and x vp,G , x b , and x c are defined in Table 1. Without any prediction model and letting N IJ � I T Σ −1 J, U I � I T Σ −1 l, the least-square norm equation of (14) is Using the reduction algorithm, norm equation (15) is equal to  Other real states related to receiver, atmospheric delays, and satellite-related errors n c × 1 with Combining (11) and (16), there is With predict model (7), the integration solution is To unify the state with INS, the measurement model containing full state is Without any prediction model and letting A � AT g,i N IJ � I T Σ −1 J, U I � I T Σ −1 l, the least-square norm equation of (20) is Without enough measurements, the coefficient matrix of the norm function could be ill-conditioned and the state vector is unsolvable or can only be solved with poor precision.
In (9), only position and velocity states are considered in the integration model; the dynamic model of state x b and x c should be considered within the GNSS processor. Supposing the corresponding prediction is x b , P b and x c , P c , with P * denoting the corresponding predict variance matrix, and letting N IJ � N IJ + P U I � U I + P − 1 x, (21) is extended as Using the reduction algorithm, norm equation (22) is equal to with To investigate the solution difference between loosely and tightly coupled integration, the general mathematical model for GNSS positioning is given as where x G is the unknown state vector; l is the measurement vector corrected by an initial state va; and ε is the measurement noise vector which is usually modelled as random noise. e least-squares solution is Least-squares solution (26) will be used as the measurement model of LC integration, and mathematical model (25) will be used for the measurement model of TC integration. In the following subsection, the equivalence or discrepancy between loosely and tightly coupled algorithms will be derived.

Equivalence Proof between Loosely and Tightly Coupled
Integration for GNSS Single Epoch Solution. With subscript "LC" denoting the LC scheme, the corresponding measurement model is It is derived from the difference of the GNSS solution and INS solution. From (10), there is e corresponding vector form is erefore, it is derived as GNSS positioning processing requires for an approximate value of the state. In a GNSS/INS integration procedure, the approximate value can come from the INS time update and is given as erefore, where x G is from (26). Substituting (33) into (30), there is From (8), the LC integration solution is From (25) and (26), the corresponding TC measurement model is erefore, the design matrix is and the corrected measurement vector is with its variance being According to (38), (39), and (40), the TC integration solution is (36) and (41) reveal that the loosely and tightly coupled schemes are equivalent in terms of the solution of x k and the estimated variance P k , while the gain matrix K k satisfied that It is noted that the above equivalence is satisfied based on three implied conditions: (1) in TC processor, the INS time update is used as the initial value for the linearization of the GNSS observation function; (2) the GNSS processor provides a standalone solution, which means no dynamic model used; (3) the measurement variances of LC and TC schemes are given by (35) and (40), respectively, which build the inner connection between them.

Equivalence Proof Based on Simplified Equivalent
Equations at Single Epoch. Abovementioned general state function includes all parameters, while some of the parameters are not directly necessary for navigation, such as atmospheric parameters, GNSS receiver clock parameters, and ambiguities. Considering all parameters in the state vector would unavoidably increase the computation burden. erefore, it is reasonable to utilize a simplified state model which only contains the public state in the integration filter. Dividing the unknown states of GNSS positioning into two parts, the navigation-related and navigation-unrelated (25) is rewritten as where x vp,G denotes the GNSS position and velocity states as given in Table 1. By using the equivalent function, the above function can be rewritten as where It can be proved that the matrix R 1 is idempotent and Σ −1 R 1 is symmetric, see [25].
It was proved that the estimator of x vp,G and its covariance matrix D(x vp,G ) can be written as Mathematical Problems in Engineering which are equivalent to the solution of (26). Relying on (46), it can be proved that the LC and TC integrations are also equivalent in this case With the state of (47), the LC measurement model is erefore, the LC KF solution is Accordingly, the TC measurement model is erefore, the TC KF solution is (49) and (51) indicate that the equivalence between LC and TC is also satisfied in this case. It is noted that the equivalence between (36) and (49) and between (41) and (51) are not satisfied. In (36) and (41), the dynamic characteristic of GNSS receiver's clock bias and drift are considered in the KF predict model, while the corresponding dynamic information is ignored in (49) and (51).

Measurement Models for GNSS Measurements at Multiple Epochs
It is known that the LC integration requires less information exchange for integration but would be suspended if GNSS standalone solution failed. Comparably, the TC processor can improve the continuity and stability of the system, but it requires heavier information exchange and therefore complicates the computation. e above derivation concludes that the LC and TC algorithms are equivalent based on some implied conditions. erefore, the most efficient procedure is to trigger the LC processor when the GNSS solution is provided and return to the TC processor when the available satellites are less than 4. By this way, both the computation efficiency and the navigation reliability can be guaranteed. GNSS precise positioning using carrier phase measurements implies a certain potential dynamic model, i.e., the ambiguity is constant when there is no cycle slip, and the residual ionospheric delay is temporarily stable. In this case, the solution is not "standalone" but time correlated. erefore, the equivalence between LC and TC integrations should be rediscussed.
Firstly, the ambiguity fixing model is particularly discussed. For GNSS positioning procedure, the least-squares algorithm is applied at the 1 st epoch and the solution of position and velocity parameters is given as (46). It is noticed that (46) only makes senses when both pseudorange and carrier phase measurements are used, since the matrix M 11 would be ranked deficient if only the carrier phase measurements are considered. e single epoch least-squares solution of ambiguities is Without enough information, M 22 is singular, so usually the norm equation is recorded as At current epoch k, the normal equation is accumulated as where i and k are the epoch indicators. Accordingly, the solution of x vp,G with accumulated ambiguity float solution is and the variance matrix is Intrinsically, the observation function for solutions (54), (55), and (56) can be rewritten as where the subscript "k" is the epoch indicator. e single epoch solution of GNSS position and velocity is given by (46), and the covariance between x vp,G and x b is As presented in the previous section, it is concluded that the LC measurement mode with (46), (52), and (58) is equivalent to the TC measurement mode with (25). For these solutions, the ambiguity should be included in the integration state vector and its single epoch float solution should be used as the measurement update, so that the multiple epoch ambiguity information is reserved in the integration procedure. Otherwise, if only the position or/and velocity from GNSS standalone solution is involved in the integration procedure, the ambiguity cannot be fixed, since the ambiguity dynamic constraints are not reserved. e LC procedure with states including ambiguities and measurement update from (54), (55), and (56) is equivalent to the TC procedure of using (57) as the measurement update. is mode uses the multiple epoch ambiguity solution twice, both in the integration dynamic and measurement model. erefore, the integration solution is naturally different with that of using (46), (52), and (58). e above derivation is easily extended to the case if the GNSS procedure contains another dynamic model. Generally, the GNSS solution with a dynamic model is given as where w is the state dynamic noise, and its variance is D(w) � D(x G,0 ) � Σ x ; the measurement model is the same as (25). Supposing the solution of (59) is written as x G and D(x G ), therefore, the LC procedure of using x G and D(x G ) as the measurement update is equivalent to the TC procedure of using measurement update, as (60) e corresponding dynamic model containing IMU information is Supposing (62) Integrated estimation with (60) and (61) is also equivalent to the following estimation: (63) means that the GNSS state model (59) is finally added to the integration state model (61). is is the reason that federal KF requires information allocation for the dynamic model among each subfilter [26,27]. e same dynamic model used in each subfilter indicates that the dynamic variance in the main filter is reduced by a factor of 1/n, where n is the number of the subfilter. e models used for GNSS/INS integration are summarized in Table 2, with the equivalent LC and TC measurement models presented at the last two columns. If the GNSS solution is unavailable, the normal equation is used as the measurement model for integration. It should be noted that the equivalence is satisfied within each model and shows discrepancies among different models, i.e., models 1 and 2 are different due to the dynamic of x t . If x t is modelled as random noise and the dynamic noise tends to infinite, models 1 and 2 would be approximately equivalent. Models 3, 4, and 5 are for RTK solution. With the same dynamic model for x b , the integration solution of models 3 and 4 would be similar, while model 5 cannot fix the ambiguity.

Experiment and Analysis
e above derivation proves that the LC and TC integrations are naturally equivalent when the same information is used. Usually, the measurement stochastic model of the LC integration is empirically determined rather than derived from the GNSS procedure.
is is the reason that caused the discrepancies. A flight test dataset is used to verify the discrepancies among different models. e flight test, with trajectory shown in Figure 1, was conducted on October 17, 2011, from Bankstown Airport, Sydney to Cooma Airport, which is in the southern part of New South Wales, Australia. e GNSS data was collected by a Leica dual-frequency receiver with a sample rate of 10 Hz. e cutoff elevation angle was 10°and the number of visible GPS satellites was nine throughout the experimental period. e INS data collected the KVH's IMU embedded in the SPAN-CPT with the specification listed in Table 3, and the measurement rate was set to 100 Hz.

Comparisons of Using GNSS Pseudoranges.
Firstly, the solutions with GNSS SPP processor are compared and 15 states with position, velocity, attitude, acceleration, and gyroscope errors are selected. State model 2 in Table 2 is Mathematical Problems in Engineering select as demonstration, and four different measurement models listed in Table 4 are compared. e integration solution from the Inertial Explorer software is used as the reference solution and the position, velocity, and attitude errors of the four models are shown in Figures 2-4. Figure 2 shows that the horizontal position errors are all within 2 meters and vertical position errors are around 5 meters. e position results show some systematic biases, about 1 meter on horizontal and 2 meters on vertical direction. Generally, the solutions for models B, C, and D are quite close and much more stable than that of model A.
eoretical derivation proves that the solutions of model B and D are equivalent. e zoom-in plots in each panel of Figure 2 show that the results of models B and D are same in millimeter level and their discrepancies with model C are also small, usually around several centimeters.   Figure 1: Trajectory of the flight test.
Mathematical Problems in Engineering   Figure 4 shows tthe attitude errors of the four models, within 0.2 degrees on pitch and roll and 0.5 degrees on yaw after converging. e top panel shows that pitch errors of model A have significant biases with that of models B, C, and D. Middle and bottom panels show that biases on roll and yaw between different models are smaller.
Generally, results in Figures 2-4 show that the LC and TC integrations provide equivalent solutions once the LC measurement variances are estimated by the GNSS processor, which is a function of the TC measurement variances. Also, it shows that if empirical constant variances were used for LC measurement update, the solution accuracy and stability would degrade. Table 5 lists the mean and standard deviation (STD) values of the position, velocity, and attitude errors by using the four models. Generally, performance of model A is the worst, and performances of models B and D are almost the same and little better than model C. In terms of absolute mean values, the discrepancies among the four models are generally small, except for the yaw errors, where errors of models B, C and D are much smaller than those of model A. e differences on STD values are much more significant. For position solutions, improvements of models B, C, ND D compared with model A are around 15% on horizontal and 40% on vertical direction. e corresponding improvements on velocity solutions are much larger, around 51% on horizontal and 66% on vertical direction. e improvements on pitch and roll solutions are little, but quite significant on the yaw axis, around 25%. Results of models B and D are always the same on the third decimal place. Model C also shows little differences with models B and D, which indicates that the covariance of position and velocity can be neglected in some cases when the correlation is insignificant.

Comparison of Using GNSS Pseudoranges and Carrier
Phases. In this subsection, the integration solutions with the precise GNSS RTK processor are compared and 18 states with three lever arm error states are considered. e state model 5 in Table 2 is selected for demonstration and the four measurement models for integration in Table 4 are compared with the measurement variances of model A as R � diag 0.01 0.01 0.16 0.04 0.04 0.16 , since the GNSS processor provides RTK solutions. e corresponding position, velocity, and attitude errors of the four models are shown in Figures 5-7. Figure 5 shows that the position errors are all within 0.4 meters and the horizontal accuracies are obviously higher than the vertical accuracy. Generally, similar to the results of using pseudoranges, the errors of model A fluctuate more severely compared with other three models. e zoom-in figures show that the discrepancies among models B, C, and D are always small at millimeter level for horizontal and centimeter level for vertical position errors. Figure 6 shows the velocity errors of using the four models. e velocity error discrepancies among models are not as obvious as that of the position errors, and the results of models B and D are still close. Figure 7 shows the corresponding attitude errors. e top panel shows that the pitch errors of model A are even smaller than that of other three models, especially during the 200∼1200 second.
is may be caused by the uncertainty of the carrier phase variances. e middle panel shows that the roll errors of the four models are quite close, but the errors of model A seem slightly smaller. e bottom panel shows that the yaw solutions of models B, C, and D are obviously higher than that of model A, especially during the first 700 seconds.
Generally, Figures 2-7 reveal that models B and D can provide equivalent solutions, while their differences with model C is little but large with model A. In most cases, using the estimated variances from the GNSS processor as the integration measurements' variances can improve the accuracy and stability of the system. e corresponding absolute mean and STD values of the position, velocity, and attitude errors are listed in Table 6. Generally, results of models B and D are mostly close.
e STD values for position errors show that the estimation precision of the four models are comparable, while the absolute mean values of models B, C, and D are much better than that of model A, indicating higher estimation accuracy. e velocity estimation accuracy and   estimation accuracy and precision of models B, C, and D are improved around 25%.

Concluding Remarks
is paper proves the equivalence of the GNSS/INS loosely and tightly coupled integration based on four conditions: (1) in TC processor, the INS time update is used as the initial value for the linearization of the GNSS observation function; (2) the GNSS processor provides a standalone solution, which means no dynamic model was used for the solution; (3) if dynamic constraints are considered in the GNSS solution for the LC processor, the same dynamic information should be added to the dynamic model of the TC processor; (4) the covariance of the LC measurements should be derived from the GNSS solution instead of conventional empirical values. Upon this, different state models for integration are given to illustrate the natural differences between different filtering. is indicates that once the same measurement and predict information are used, the integration would be equivalent no matter what kind of coupling schemes are used. erefore, to reserve the estimation efficiency, accuracy, and stability, the optimal procedure is to use the LC processor when the GNSS standalone solution is available and to use the TC processor when GNSS measurements are not enough. Alternatively, one can always use the norm equation from the GNSS processor for the integration measurement update, so that a unified model can be implemented in a practical processor.

Data Availability
e GNSS and IMU raw data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare no conflicts of interest.