# A Tightly Coupled BDS/INS Integrated Positioning Algorithm Based on Triple-Frequency Single-Epoch Observations

Authors: Fei Ye; Shuguo Pan; Wang Gao; Hao Wang; Chun Ma; Yunfeng Wang
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Vehicular dynamic positioning based on tightly coupled (TC) Global Navigation Satellite System (GNSS)/Inertial Navigation System (INS) integration in urban areas is due to either low accuracy of pseudorange or poor continuity of carrier phase, resulting in insufficient positioning performance. To enhance the stability while ensuring positioning accuracy, this paper proposed a tightly coupled Beidou Navigation Satellite System (BDS)/INS integration scheme by improving measurement modelling with triple-frequency observations: first, a stepwise single-epoch ambiguity resolution of extra-wide-lane (EWL)/wide-lane (WL) combined observations and then modelling the measurement equation with fixed WL observation instead of conventional pseudorange or carrier phase. Experiments were carried out for verification with data collected in real traffic by a measurement vehicle. The proposed method achieved single-epoch output with an RMS statistical accuracy of decimetre level of 0.152 m horizontally and 0.196 m vertically. The signal outage experiment verified that the proposed algorithm is restoring high-accuracy positioning output in single-epoch once the signal is recaptured. The proposed method obtained a positioning accuracy improvement of 43.6% horizontally and 6.2% vertically in signal outage sections compared to the conventional method. This avoids the multiepoch ambiguity searching to fix with conventional carrier-phase processing, thereby improving the positioning stability.


---

## Body

## 1. Introduction

With the rapid development of autonomous vehicles, positioning performance in the urban environment has become a hot research point. The tightly coupled (TC) Global Navigation Satellite System (GNSS)/Inertial Navigation System (INS) integration technique is used to provide great absolute positioning output, but the easily blocked satellite signals in an urban environment make the integrated system hard to perform as stably as in an open area. Urban canyons, tunnels, and boulevards covered by thick trees are typical GNSS restrict environments. In these areas, factors such as bad satellite geometric distribution, weak or blocked signal, and even multipath are the most significant error sources which must be considered [1]. Unstable GNSS observations lead to unstable error correction, which forces integrated positioning to more rely on high-end INS performance. Therefore, the current autonomous vehicle environment perception solutions prefer more sensors to perform relative positioning. The problems mentioned above directly lead to a substantial increase in hardware costs and system complexity.The principle of tightly coupled integration is using two sensors on measuring the same target and then calibrating one of the sensors with the bias measured, to obtain the final output [2]. In the case of GNSS/INS integration, both the GNSS and INS are observing the geometry distance and its changing rate between satellites and the rover station [3]. On the assumption that the GNSS observation value is stable and reliable, the GNSS observations are used to correct the INS measurements to avoid the divergence, thereby achieving calibration [4]. However, GNSS observations in an urban environment cannot guarantee its own stability. In order to improve the positioning performance of the tightly coupled integration, some studies start from the state of vehicle movement in the urban environment. For example, velocity and the small changing altitude are used as constraints for maintaining the INS error bounded in the case of GNSS outage [5]. Some scholars improved the fusion algorithm of tightly coupled integration to enhance its robustness and adaptive performance in urban environments [6, 7]. Besides, researchers also investigated innovative improvements into the environmental constraints, i.e., the context-aided framework [8]. However, the study about GNSS observations’ stability improving in complex environments is relatively small.The most common TC is using pseudorange as GNSS observation [2]. With the technique developed, the pseudorange based meter level accuracy cannot fulfil the requirement of high-accuracy applications. As a higher accuracy observation, carrier phase is first used to smooth pseudorange [9]. Though the accuracy is improved indeed, the real-time performance is reduced. Some studies directly use carrier-phase observation instead [10, 11]. However, the restriction of using carrier-phase observations is the ambiguity resolution (AR). Generally, the ambiguity resolution requires several continuous epochs to fix. And once the loose lock happened, it must be refixed again. In the case of GNSS/INS integration, during the refixing period, either use uncontrollable float carrier-phase observations or let INS reckon independently without calibration. The AR process of narrow-lane (NL) carrier phase includes Kalman filter (KF) and LAMBDA (Least-squares Ambiguity Decorrelation Adjustment) algorithm [12]. The serial use of this process and GNSS/INS integrated filtering greatly increases the complexity [13, 14].In order to find a solution that has both simple calculations as using pseudorange and high accuracy as using carrier phase, AR must be processed as simple as possible. Some use the time difference between neighbour epochs [15]. This method can receive a much higher accuracy than using pseudorange. Its error accumulates as time goes on, so this method is only suitable for short-term filtering. Some researchers use wide-lane (WL) observations to form a double-differenced (DD) measurement model and then calculate ambiguity aided by INS [16]. Considering dual frequency, stepwise WL and NL ambiguity resolutions were used for AR in precise point positioning (PPP) [17]. After BeiDou Navigation Satellite System (BDS) started to broadcast triple-frequency signal in 2012, multisystem multifrequency is feasible for a better GNSS positioning model [18] and positioning with triple frequency received wide attention [19]. Using triple-frequency carrier observations can receive a higher AR success rate than using dual frequency [20]. Because extra-wide-lane (EWL) and WL have longer wavelength, AR is realized by rounding to integer [21]. Therefore, EWL/WL technique is able to be used for single-epoch positioning in real-time kinematic (RTK) [22, 23]. Researchers keep on improving EWL/WL single-epoch AR with triple frequency based on the geometry-free and ionosphere-free (GIF) model and three-carrier ambiguity resolution (TCAR) method [24–28]. Besides, positioning with triple frequency allows not only avoiding but also detecting and repairing cycle slip [29]. Thus, EWL/WL based single-epoch positioning is ideal for improving continuity while ensuring accuracy to decimetre or subdecimetre level [22].The research we are working on aims to use the excellent performance of BDS to achieve stable high-accuracy vehicular positioning in urban environments. In this paper, we improved TC measurement modelling with three-frequency single-epoch AR technology and BDS observations. With the proposed method, on the one hand, INS is calibrated with ambiguity fixed observations once a satellite signal is observable, thereby improving the positioning accuracy and continuity. On the other hand, due to the independence of epochs, the interference will not pollute the output of neighbour epochs, thereby improving the positioning stability.The rest of the paper is organized as follows. In Section2, based on the BDS triple-frequency single-epoch observations, a tightly coupled BDS/INS integration with WL observation is proposed, and the corresponding integration architecture is given. In order to verify the effectiveness and advantages of the proposed method, a series of experiments are carried out in Section 3, together with detailed experimental procedures, and results are given. Section 4 analyses the experimental results in detail and discusses the advantages and disadvantages of the proposed method. The conclusion and future work are given in Section 5.


## 2. Materials and Methods

### 2.1. BDS Triple-Frequency Single-Epoch Observations

Triple frequency of pseudorange and carrier-phase linear combinations is formed first. Then, the EWL fixed carrier phase is obtained by the geometric-free ionosphere-free model. With help of the fixed EWL observation and TCAR algorithm, the WL ambiguity is fixed by rounding to integer; thus, the accurate WL observation is obtained. The specific progress in detail is as follows.EWL and WL are different linear combinations of observations. The combined observations have longer wavelengths, which are easier to fix the ambiguity. EWL has a longer wavelength than WL. For BDS observations, the combination coefficient of EWL is0,−1,1 and for WL is 1,−1,0 [30]. The linear combinations are formed following the equations:(1)Δ∇ϕi,j,k=if1Δ∇ϕ1+jf2Δ∇ϕ2+kf3Δ∇ϕ3if1+jf2+kf3,Δ∇Pi,j,k=if1Δ∇P1+jf2Δ∇P2+kf3Δ∇P3if1+jf2+kf3,where f1,f2,andf3 are three frequencies. i,j,andk are the combination coefficients. Δ∇ represents the double-differenced operator. ϕn and Pn are the carrier-phase and pseudorange observations of the n-th frequency in the unit of meter. The geometry-free and ionosphere-free (GIF) model is used to calculate the ambiguity parameter of EWL [31]. After double difference between stations and satellites in the short baseline, the ionospheric error was eliminated. The calculation is shown as follows:(2)Δ∇NEWL=Δ∇ϕEWL−Δ∇PEWLλEWL,where λEWL is the wavelength of EWL observation and the Δ∇NEWL is the EWL ambiguity, which is fixed by rounding to integer. The fixed EWL ambiguity will be used as a known value to participate in the calculation of WL ambiguity based on the TCAR method as follows:(3)Δ∇NWL=1λWLΔ∇ϕWL−Δ∇ϕEWL−ηEWL−ηWLΔ∇I+λEWLΔ∇NEWL,where λWL is the wavelength of WL observation and Δ∇NWL is the WL ambiguity. I is the ionospheric delay, together with the coefficient η, the influence of which is controllable. Same as Δ∇NEWL, Δ∇NWL is also rounding to fix. Once the ambiguity of WL is fixed, the carrier-phase observation is also fixed:(4)Δ∇ϕWL=Δ∇ρ+Δ∇T−ηWLΔ∇I+λWLΔ∇NWL+Δ∇εWL,where T, I, and ε are received by interpolation of area modelling of network RTK. Constants as ηEWL=−0.352 and λWL=0.8470m are given in Gao’s research [31]. The only unknown parameter is the double-differenced satellite-station geometry distance Δ∇ρ. For dynamic application, Δ∇ρ is used for modelling the measurement equation of the tightly coupled integration.


### 2.2. Tightly Coupled BDS/INS Integration

#### 2.2.1. System Model

The system model of tightly coupled integration system is established in the following form [32]:(5)x˙t=fxt+wt,where f⋅ is the nonlinear function of system dynamics with state xt and wt is the system noise vector. The discrete form of (5) is described as follows:(6)Xk=Φk,k−1Xk−1+Wk−1,where Xk is the estimated state vector in the kth epoch, while Xk−1 is the state vector of the last epoch; Φk,k−1 is the discrete form of f(∙); and Wk−1 is the discrete version of system noise vector. The system modelling is derived from the dynamic equation of INS [33, 34].Since the double-differenced observations do not need to estimate the satellite clock bias, the observable errors are only INS errors. Considering ECEF-frame, the state vector is selected as the most basic 15-state solution, which is(7)X=δr1×3δv1×3φ1×3bg,1×3ba,1×3T,where δr stands for positioning error vector in three dimensions. Similar, δv and φ stand for velocity and attitude error vectors. bg and ba are measurement drift/bias vectors of gyros and accelerometers. With the state vector determined, the system error dynamics of GNSS/INS integration is so that becomes basically an INS mechanization equation. It is written in ECEF-frame according to [34] as(8)δr˙eδv˙eφ˙e=δve+εre−2Ωiee×δve+Cbefbφe+Cbebab+εve−Ωiee×φe+Cbebgb+εφe,where Ωie is the skew-symmetric form of the Earth rotation rates and Cbe represents the rotation matrix from base-frame to ECEF, here to rotate the specific force fb to ECEF. εr, εv, and εφ indicate random walk process driving noise vectors for the position, velocity, and attitude, respectively.Initially, the GNSS/INS integrated positioning algorithm was designed for highly dynamic aircraft, and the dynamics of the vehicular platform (such as autonomous vehicles) is not as complicated as that of aircraft. Generally, in vehicular applications, the parameters such as height errors and three-dimensional acceleration do not change frequently and greatly, so the detailed modelling with those parameters can be simplified to a certain extent in practical applications. For example, it is also mentioned in [5] that the velocity of the vehicular platform in the lateral and elevation directions can be almost assumed to be zero. Besides, if low-cost MEMS IMU is used and the instrumental error of the gyro itself is far greater than the influence of other biases, then only εφe needs to be retained for the estimation of the attitude error. When the simplified system model participates in the estimation, the computation can be further reduced and lay a foundation for the feasibility of autonomous vehicle mass applications.Finally, the system model is presented based on the system error dynamics as follows:(9)δr˙δv˙φ˙b˙gb˙a=0I0000−2ΩieCbefb0Cbe00−ΩieCbe00000000000⋅δrδvφbgba+εrεvεφ00.


#### 2.2.2. Measurement Model

Similar with the system model, the measurement model of integrated positioning system describes the relationship between observation and the estimated states. The origin measurement model is normally formed as follows:(10)Zk=HkXk+Vk,whereZk represents the measurement vector at epoch k. Hk is the measurement mapping matrix, which expresses relationship between states and measurement vectors. Vk denotes measurement error. Normally, measurement vector Zk in the TC model reflects the measurement difference of two sensors on the same object. Geometry distance is the most common selection. In case of this paper, measurement vector is expressed as follows:(11)Zk=Δ∇ρG,k−Δ∇ρI,k,where Δ∇ρG,k and Δ∇ρI,k are representing DD geometry distance measured by GNSS and INS. Geometry distance between satellite and INS approximate position is used to form ρI. The BDS triple-frequency WL observation is used to form ρG instead of conventional pseudorange or carrier-phase observations. Therefore, in the rest of the section, we use ρWL to represent the observation.In the measurement equation, geometry distance received by two systems is used for modellingHk. Once Δ∇ρWL and Δ∇ρI received, the error between the two can be described by the equation as follows:(12)Δ∇ρWL−Δ∇ρI=Δ∇lG−Δ∇lIδx+Δ∇mG−Δ∇mIδy+Δ∇nG−Δ∇nIδz+ε,where Δ∇l, Δ∇m, and Δ∇n represent DD direction cosine. δx, δy, and δz mean positioning error between GNSS and INS. Subscripts G and I represent values of GNSS and INS. Let e represent the direction cosines for short, such as(13)Δ∇ek=Δ∇lG,k−Δ∇lI,kΔ∇mG,k−Δ∇mI,kΔ∇nG,k−Δ∇nI,k.Then, considering other error states, there should be ann×15-dimensional measurement matrix to describe relationship between Zk and X^k (the estimated state). n stands for satellite numbers in the epoch involved. Then, Hk is written as follows:(14)Hk=Δ∇ek,n×1201×12.The DD eliminates clock error between satellite and receiver; thus, the equivalent distance caused by clock error is no longer involved. Finally, the measurement equation is written as follows:(15)Δ∇ρWL,k1−Δ∇ρI,k1Δ∇ρWL,k2−Δ∇ρI,k2⋮Δ∇ρWL,kn−Δ∇ρI,kn=Δ∇e1,k01×12Δ∇e2,k01×12⋮Δ∇en,k01×12⋅X+ε1,kε2,k⋮εn,k.


#### 2.2.3. Integration Scheme

In GNSS/INS integration, with the usage of conventional carrier-phase observations, it is necessary to use a first KF for float ambiguity estimation and LAMBDA for fixed ambiguity searching; then, the second KF for GNSS/INS fusion can be performed. Or, in another solution, the ambiguity is appended to the state vector and estimated together. Although double KF usage is avoided, the increased ambiguity parameter greatly increases the matrix dimension and also the corresponding computation pressure. Therefore, in the case when the conventional carrier-phase observation processing is used, these two schemes have both shortcomings for an autonomous vehicle to use. With the proposed algorithm, the ambiguity is fixed in a single epoch. The benefits of high-accuracy positioning and less computation are both covered. Figure1 presents the schematic diagram. The meanings of the signs used in the diagram are the same as those in the model explanation in the rest of the sections.Figure 1
The proposed integration scheme.Cubature Kalman filter (CKF) [35] is the most popular fusion estimator for the GNSS/INS integration in recent years. Compared with the conventional KF variants, CKF has higher accuracy and less computation. As presented in the architecture, the fusion algorithm we used for BDS/INS integration is CKF. Estimation algorithms based on Gaussian filtering generally require state and measurement modelling, and the CKF is no exception. Figure 2 presents the flowchart of CKF.Figure 2
The flowchart of CKF.


## 2.1. BDS Triple-Frequency Single-Epoch Observations

Triple frequency of pseudorange and carrier-phase linear combinations is formed first. Then, the EWL fixed carrier phase is obtained by the geometric-free ionosphere-free model. With help of the fixed EWL observation and TCAR algorithm, the WL ambiguity is fixed by rounding to integer; thus, the accurate WL observation is obtained. The specific progress in detail is as follows.EWL and WL are different linear combinations of observations. The combined observations have longer wavelengths, which are easier to fix the ambiguity. EWL has a longer wavelength than WL. For BDS observations, the combination coefficient of EWL is0,−1,1 and for WL is 1,−1,0 [30]. The linear combinations are formed following the equations:(1)Δ∇ϕi,j,k=if1Δ∇ϕ1+jf2Δ∇ϕ2+kf3Δ∇ϕ3if1+jf2+kf3,Δ∇Pi,j,k=if1Δ∇P1+jf2Δ∇P2+kf3Δ∇P3if1+jf2+kf3,where f1,f2,andf3 are three frequencies. i,j,andk are the combination coefficients. Δ∇ represents the double-differenced operator. ϕn and Pn are the carrier-phase and pseudorange observations of the n-th frequency in the unit of meter. The geometry-free and ionosphere-free (GIF) model is used to calculate the ambiguity parameter of EWL [31]. After double difference between stations and satellites in the short baseline, the ionospheric error was eliminated. The calculation is shown as follows:(2)Δ∇NEWL=Δ∇ϕEWL−Δ∇PEWLλEWL,where λEWL is the wavelength of EWL observation and the Δ∇NEWL is the EWL ambiguity, which is fixed by rounding to integer. The fixed EWL ambiguity will be used as a known value to participate in the calculation of WL ambiguity based on the TCAR method as follows:(3)Δ∇NWL=1λWLΔ∇ϕWL−Δ∇ϕEWL−ηEWL−ηWLΔ∇I+λEWLΔ∇NEWL,where λWL is the wavelength of WL observation and Δ∇NWL is the WL ambiguity. I is the ionospheric delay, together with the coefficient η, the influence of which is controllable. Same as Δ∇NEWL, Δ∇NWL is also rounding to fix. Once the ambiguity of WL is fixed, the carrier-phase observation is also fixed:(4)Δ∇ϕWL=Δ∇ρ+Δ∇T−ηWLΔ∇I+λWLΔ∇NWL+Δ∇εWL,where T, I, and ε are received by interpolation of area modelling of network RTK. Constants as ηEWL=−0.352 and λWL=0.8470m are given in Gao’s research [31]. The only unknown parameter is the double-differenced satellite-station geometry distance Δ∇ρ. For dynamic application, Δ∇ρ is used for modelling the measurement equation of the tightly coupled integration.


## 2.2. Tightly Coupled BDS/INS Integration

### 2.2.1. System Model

The system model of tightly coupled integration system is established in the following form [32]:(5)x˙t=fxt+wt,where f⋅ is the nonlinear function of system dynamics with state xt and wt is the system noise vector. The discrete form of (5) is described as follows:(6)Xk=Φk,k−1Xk−1+Wk−1,where Xk is the estimated state vector in the kth epoch, while Xk−1 is the state vector of the last epoch; Φk,k−1 is the discrete form of f(∙); and Wk−1 is the discrete version of system noise vector. The system modelling is derived from the dynamic equation of INS [33, 34].Since the double-differenced observations do not need to estimate the satellite clock bias, the observable errors are only INS errors. Considering ECEF-frame, the state vector is selected as the most basic 15-state solution, which is(7)X=δr1×3δv1×3φ1×3bg,1×3ba,1×3T,where δr stands for positioning error vector in three dimensions. Similar, δv and φ stand for velocity and attitude error vectors. bg and ba are measurement drift/bias vectors of gyros and accelerometers. With the state vector determined, the system error dynamics of GNSS/INS integration is so that becomes basically an INS mechanization equation. It is written in ECEF-frame according to [34] as(8)δr˙eδv˙eφ˙e=δve+εre−2Ωiee×δve+Cbefbφe+Cbebab+εve−Ωiee×φe+Cbebgb+εφe,where Ωie is the skew-symmetric form of the Earth rotation rates and Cbe represents the rotation matrix from base-frame to ECEF, here to rotate the specific force fb to ECEF. εr, εv, and εφ indicate random walk process driving noise vectors for the position, velocity, and attitude, respectively.Initially, the GNSS/INS integrated positioning algorithm was designed for highly dynamic aircraft, and the dynamics of the vehicular platform (such as autonomous vehicles) is not as complicated as that of aircraft. Generally, in vehicular applications, the parameters such as height errors and three-dimensional acceleration do not change frequently and greatly, so the detailed modelling with those parameters can be simplified to a certain extent in practical applications. For example, it is also mentioned in [5] that the velocity of the vehicular platform in the lateral and elevation directions can be almost assumed to be zero. Besides, if low-cost MEMS IMU is used and the instrumental error of the gyro itself is far greater than the influence of other biases, then only εφe needs to be retained for the estimation of the attitude error. When the simplified system model participates in the estimation, the computation can be further reduced and lay a foundation for the feasibility of autonomous vehicle mass applications.Finally, the system model is presented based on the system error dynamics as follows:(9)δr˙δv˙φ˙b˙gb˙a=0I0000−2ΩieCbefb0Cbe00−ΩieCbe00000000000⋅δrδvφbgba+εrεvεφ00.


### 2.2.2. Measurement Model

Similar with the system model, the measurement model of integrated positioning system describes the relationship between observation and the estimated states. The origin measurement model is normally formed as follows:(10)Zk=HkXk+Vk,whereZk represents the measurement vector at epoch k. Hk is the measurement mapping matrix, which expresses relationship between states and measurement vectors. Vk denotes measurement error. Normally, measurement vector Zk in the TC model reflects the measurement difference of two sensors on the same object. Geometry distance is the most common selection. In case of this paper, measurement vector is expressed as follows:(11)Zk=Δ∇ρG,k−Δ∇ρI,k,where Δ∇ρG,k and Δ∇ρI,k are representing DD geometry distance measured by GNSS and INS. Geometry distance between satellite and INS approximate position is used to form ρI. The BDS triple-frequency WL observation is used to form ρG instead of conventional pseudorange or carrier-phase observations. Therefore, in the rest of the section, we use ρWL to represent the observation.In the measurement equation, geometry distance received by two systems is used for modellingHk. Once Δ∇ρWL and Δ∇ρI received, the error between the two can be described by the equation as follows:(12)Δ∇ρWL−Δ∇ρI=Δ∇lG−Δ∇lIδx+Δ∇mG−Δ∇mIδy+Δ∇nG−Δ∇nIδz+ε,where Δ∇l, Δ∇m, and Δ∇n represent DD direction cosine. δx, δy, and δz mean positioning error between GNSS and INS. Subscripts G and I represent values of GNSS and INS. Let e represent the direction cosines for short, such as(13)Δ∇ek=Δ∇lG,k−Δ∇lI,kΔ∇mG,k−Δ∇mI,kΔ∇nG,k−Δ∇nI,k.Then, considering other error states, there should be ann×15-dimensional measurement matrix to describe relationship between Zk and X^k (the estimated state). n stands for satellite numbers in the epoch involved. Then, Hk is written as follows:(14)Hk=Δ∇ek,n×1201×12.The DD eliminates clock error between satellite and receiver; thus, the equivalent distance caused by clock error is no longer involved. Finally, the measurement equation is written as follows:(15)Δ∇ρWL,k1−Δ∇ρI,k1Δ∇ρWL,k2−Δ∇ρI,k2⋮Δ∇ρWL,kn−Δ∇ρI,kn=Δ∇e1,k01×12Δ∇e2,k01×12⋮Δ∇en,k01×12⋅X+ε1,kε2,k⋮εn,k.


### 2.2.3. Integration Scheme

In GNSS/INS integration, with the usage of conventional carrier-phase observations, it is necessary to use a first KF for float ambiguity estimation and LAMBDA for fixed ambiguity searching; then, the second KF for GNSS/INS fusion can be performed. Or, in another solution, the ambiguity is appended to the state vector and estimated together. Although double KF usage is avoided, the increased ambiguity parameter greatly increases the matrix dimension and also the corresponding computation pressure. Therefore, in the case when the conventional carrier-phase observation processing is used, these two schemes have both shortcomings for an autonomous vehicle to use. With the proposed algorithm, the ambiguity is fixed in a single epoch. The benefits of high-accuracy positioning and less computation are both covered. Figure1 presents the schematic diagram. The meanings of the signs used in the diagram are the same as those in the model explanation in the rest of the sections.Figure 1
The proposed integration scheme.Cubature Kalman filter (CKF) [35] is the most popular fusion estimator for the GNSS/INS integration in recent years. Compared with the conventional KF variants, CKF has higher accuracy and less computation. As presented in the architecture, the fusion algorithm we used for BDS/INS integration is CKF. Estimation algorithms based on Gaussian filtering generally require state and measurement modelling, and the CKF is no exception. Figure 2 presents the flowchart of CKF.Figure 2
The flowchart of CKF.


## 2.2.1. System Model

The system model of tightly coupled integration system is established in the following form [32]:(5)x˙t=fxt+wt,where f⋅ is the nonlinear function of system dynamics with state xt and wt is the system noise vector. The discrete form of (5) is described as follows:(6)Xk=Φk,k−1Xk−1+Wk−1,where Xk is the estimated state vector in the kth epoch, while Xk−1 is the state vector of the last epoch; Φk,k−1 is the discrete form of f(∙); and Wk−1 is the discrete version of system noise vector. The system modelling is derived from the dynamic equation of INS [33, 34].Since the double-differenced observations do not need to estimate the satellite clock bias, the observable errors are only INS errors. Considering ECEF-frame, the state vector is selected as the most basic 15-state solution, which is(7)X=δr1×3δv1×3φ1×3bg,1×3ba,1×3T,where δr stands for positioning error vector in three dimensions. Similar, δv and φ stand for velocity and attitude error vectors. bg and ba are measurement drift/bias vectors of gyros and accelerometers. With the state vector determined, the system error dynamics of GNSS/INS integration is so that becomes basically an INS mechanization equation. It is written in ECEF-frame according to [34] as(8)δr˙eδv˙eφ˙e=δve+εre−2Ωiee×δve+Cbefbφe+Cbebab+εve−Ωiee×φe+Cbebgb+εφe,where Ωie is the skew-symmetric form of the Earth rotation rates and Cbe represents the rotation matrix from base-frame to ECEF, here to rotate the specific force fb to ECEF. εr, εv, and εφ indicate random walk process driving noise vectors for the position, velocity, and attitude, respectively.Initially, the GNSS/INS integrated positioning algorithm was designed for highly dynamic aircraft, and the dynamics of the vehicular platform (such as autonomous vehicles) is not as complicated as that of aircraft. Generally, in vehicular applications, the parameters such as height errors and three-dimensional acceleration do not change frequently and greatly, so the detailed modelling with those parameters can be simplified to a certain extent in practical applications. For example, it is also mentioned in [5] that the velocity of the vehicular platform in the lateral and elevation directions can be almost assumed to be zero. Besides, if low-cost MEMS IMU is used and the instrumental error of the gyro itself is far greater than the influence of other biases, then only εφe needs to be retained for the estimation of the attitude error. When the simplified system model participates in the estimation, the computation can be further reduced and lay a foundation for the feasibility of autonomous vehicle mass applications.Finally, the system model is presented based on the system error dynamics as follows:(9)δr˙δv˙φ˙b˙gb˙a=0I0000−2ΩieCbefb0Cbe00−ΩieCbe00000000000⋅δrδvφbgba+εrεvεφ00.


## 2.2.2. Measurement Model

Similar with the system model, the measurement model of integrated positioning system describes the relationship between observation and the estimated states. The origin measurement model is normally formed as follows:(10)Zk=HkXk+Vk,whereZk represents the measurement vector at epoch k. Hk is the measurement mapping matrix, which expresses relationship between states and measurement vectors. Vk denotes measurement error. Normally, measurement vector Zk in the TC model reflects the measurement difference of two sensors on the same object. Geometry distance is the most common selection. In case of this paper, measurement vector is expressed as follows:(11)Zk=Δ∇ρG,k−Δ∇ρI,k,where Δ∇ρG,k and Δ∇ρI,k are representing DD geometry distance measured by GNSS and INS. Geometry distance between satellite and INS approximate position is used to form ρI. The BDS triple-frequency WL observation is used to form ρG instead of conventional pseudorange or carrier-phase observations. Therefore, in the rest of the section, we use ρWL to represent the observation.In the measurement equation, geometry distance received by two systems is used for modellingHk. Once Δ∇ρWL and Δ∇ρI received, the error between the two can be described by the equation as follows:(12)Δ∇ρWL−Δ∇ρI=Δ∇lG−Δ∇lIδx+Δ∇mG−Δ∇mIδy+Δ∇nG−Δ∇nIδz+ε,where Δ∇l, Δ∇m, and Δ∇n represent DD direction cosine. δx, δy, and δz mean positioning error between GNSS and INS. Subscripts G and I represent values of GNSS and INS. Let e represent the direction cosines for short, such as(13)Δ∇ek=Δ∇lG,k−Δ∇lI,kΔ∇mG,k−Δ∇mI,kΔ∇nG,k−Δ∇nI,k.Then, considering other error states, there should be ann×15-dimensional measurement matrix to describe relationship between Zk and X^k (the estimated state). n stands for satellite numbers in the epoch involved. Then, Hk is written as follows:(14)Hk=Δ∇ek,n×1201×12.The DD eliminates clock error between satellite and receiver; thus, the equivalent distance caused by clock error is no longer involved. Finally, the measurement equation is written as follows:(15)Δ∇ρWL,k1−Δ∇ρI,k1Δ∇ρWL,k2−Δ∇ρI,k2⋮Δ∇ρWL,kn−Δ∇ρI,kn=Δ∇e1,k01×12Δ∇e2,k01×12⋮Δ∇en,k01×12⋅X+ε1,kε2,k⋮εn,k.


## 2.2.3. Integration Scheme

In GNSS/INS integration, with the usage of conventional carrier-phase observations, it is necessary to use a first KF for float ambiguity estimation and LAMBDA for fixed ambiguity searching; then, the second KF for GNSS/INS fusion can be performed. Or, in another solution, the ambiguity is appended to the state vector and estimated together. Although double KF usage is avoided, the increased ambiguity parameter greatly increases the matrix dimension and also the corresponding computation pressure. Therefore, in the case when the conventional carrier-phase observation processing is used, these two schemes have both shortcomings for an autonomous vehicle to use. With the proposed algorithm, the ambiguity is fixed in a single epoch. The benefits of high-accuracy positioning and less computation are both covered. Figure1 presents the schematic diagram. The meanings of the signs used in the diagram are the same as those in the model explanation in the rest of the sections.Figure 1
The proposed integration scheme.Cubature Kalman filter (CKF) [35] is the most popular fusion estimator for the GNSS/INS integration in recent years. Compared with the conventional KF variants, CKF has higher accuracy and less computation. As presented in the architecture, the fusion algorithm we used for BDS/INS integration is CKF. Estimation algorithms based on Gaussian filtering generally require state and measurement modelling, and the CKF is no exception. Figure 2 presents the flowchart of CKF.Figure 2
The flowchart of CKF.


## 3. Results

### 3.1. Experiment Setup

The vehicular test was carried out with a measurement vehicle (see Figure3). In the trunk of the vehicle, a set of reference systems was assembled. The reference system was set up with a whole NovAtel SPAN system with an ISA-100C IMU and a PP6 GNSS receiver (Figure 4). The reference system was used for evaluating the performance of the developed tightly coupled BDS/INS integrated positioning system. In addition, a GNSS receiver (Trimble BD990) and a microelectromechanical system (MEMS) INS (XW-5651) are deployed for applying the proposed algorithm (Figure 5). The receiver is specifically receiving BDS triple-frequency observations. Both reference and BDS receiver are connected with the same antenna with a signal splitter, to make sure they are processing the same signal. The specifics of IMUs used in the testing are listed in Table 1.Figure 3
The measurement vehicle for testing.Figure 4
NovAtel SPAN-ISA-100C system.Figure 5
The MEMS INS and BDS receiver for verifying the proposed algorithm. (a) The MEMS INS XW-5651. (b) The GNSS receiver Trimble BD990.Table 1
Specifications of IMUs.
ItemsISA-100CXW-5651GyroscopeAccelerometerGyroscopeAccelerometerBias0.5 deg/hr1250μg6 deg/hr0.1 mgScale factor100 ppm100 ppm1400 ppm1500 ppmRandom walk0.3 deg/√hr100μg/√hr0.3 deg/√hr0.25 mg/√hrDuring the vehicular field test, both the reference and developed systems are connected to continuous operating reference stations (CORS) network. Reference system using data was provided by CORS for high accuracy after positioning. The test system uses one of the base stations to form short-baseline double-difference observation. The offline process software Waypoint Inertial Explorer (IE) is used for processing high-accuracy positioning results. The test environment is chosen in a relatively open area with real traffic. In this test, the vehicle is driven at a normal traffic speed of 30–60 km/h.


### 3.2. Field Test

To test the positioning performance of the proposed algorithm, real observations of BDS are collected with help of the measurement vehicle and the reference stations. 2400 epochs of short-baseline differential data with triple-frequency BDS observations were involved in the calculation. The driving trajectory is shown in Figure6.Figure 6
Vehicle trajectory.Positioning error is obtained by the difference between process output and ground truth value. The corresponding visible satellite numbers are shown in Figure7. The errors of all epochs are plotted in Figure 8. The 95% empirical confidence which is shown in Figure 9 is presenting the reliability of the proposed algorithm. The skyplot for one of the epochs is presented in Figure 10. The accuracy represented with RMS statistics is given in Table 2.Figure 7
Observable satellite number.Figure 8
Positioning error of proposed TC algorithm based on triple frequency (m).Figure 9
The 95% empirical (red) confidence ellipses/intervals for (a) horizontal (N, E) position scatter and corresponding and (b) vertical (U) time series for the proposed model (green for inside and gray for outside).
(a)(b)Figure 10
Skyplot of one epoch.Table 2
Statistical positioning accuracy.
DirectionsRMSmE0.110N0.104U0.196Horizontal0.152


### 3.3. Signal Outage Test

A signal outage experiment is designed by artificially removing satellite signals in the experimental traffic section. Two lengths of common time interval used in inertial navigation experiments are designed in this experiment, where the signal outage durations of 10 s and 30 s were selected. The outages are designed as shown in Figure11.Figure 11
BDS signal outage arrangements.To verify the relationship between positioning results and visible satellite numbers, Figures12 and 13 show changes of visible satellites with triple-frequency observations and positioning error at the same time. The outage of visible satellite numbers in the figure is obtained by removing satellite observations of the corresponding epoch during data processing. The simulated blockage is easier to control the time length of outages. We focused on the positioning error curve around the outages. The positioning error and the corresponding satellite numbers are presented accordingly. The set of figures also includes the effect comparison of the proposed integration algorithm based on triple-frequency WL single-epoch observations versus the conventional algorithm. The conventional method refers to the method that the GNSS observation comes from the continuous measurement of multiple epochs. The float ambiguity is first estimated with KF and then fixed with LAMBDA.Figure 12
Positioning errors of 10 s outage compared with visible satellite numbers. (a) The corresponding visible satellite numbers. (b) The positioning error curve obtained with the proposed method. (c) The positioning error curve obtained with the conventional method.
(a)(b)(c)Figure 13
Positioning errors of 30 s outage compared with visible satellite numbers. (a) The corresponding visible satellite numbers. (b) The positioning error curve obtained with the proposed method. (c) The positioning error curve obtained with the conventional method.
(a)(b)(c)As presented in Figures12 and 13, the proposed algorithm instantly restores stable positioning output once the satellite signals are restored, while the conventional algorithm delays 4-5 epochs before returning to the stable output. Within the epochs, the divergence stopped under the effect of float ambiguity. But because the float ambiguity is uncontrollable, the positioning result is still not stable.We have counted the average statistical errors of three outages with the same length, as indicated in Table3 for 10 s outages. The 15 epochs before and after each outage were involved in statistics. In the table, the statistics of RMS, mean value, and maximum deviation range of the two methods are, respectively, carried out.Table 3
Average positioning error statistics during three 10 s signal outages.
ItemProcess methodEmNmUmRMSProp.0.2130.1720.422Conv.0.3090.2480.450MeanProp.−0.0820.0010.133Conv.−0.1640.0940.064Max. (abs)Prop.0.8830.6851.573Conv.0.9490.8651.289Prop. stands for the proposed method. Conv. stands for the conventional method.The average statistical errors of three 30 s outages are presented in Table4.Table 4
Average positioning error statistics during three 30 s signal outages.
ItemProcess methodEmNmUmRMSProp.1.7691.7784.283Conv.2.2542.1444.627MeanProp.−0.8710.7012.357Conv.−1.1300.8732.454Max. (abs)Prop.7.7087.09915.167Conv.8.8588.82016.179Prop. stands for the proposed method. Conv. stands for the conventional method.Compared with the conventional method, the positioning accuracy of the proposed method improves 43.6% horizontally and 6.2% vertically under the 10 s outages, with 27.5% in horizontal and 7.4% under 30 s outages (Table5).Table 5
The improvement of accuracy under different outage lengths.
Directions10 s (%)30 s (%)E31.121.5N30.717.1U6.27.4Horizontal43.627.5All44.128.5


## 3.1. Experiment Setup

The vehicular test was carried out with a measurement vehicle (see Figure3). In the trunk of the vehicle, a set of reference systems was assembled. The reference system was set up with a whole NovAtel SPAN system with an ISA-100C IMU and a PP6 GNSS receiver (Figure 4). The reference system was used for evaluating the performance of the developed tightly coupled BDS/INS integrated positioning system. In addition, a GNSS receiver (Trimble BD990) and a microelectromechanical system (MEMS) INS (XW-5651) are deployed for applying the proposed algorithm (Figure 5). The receiver is specifically receiving BDS triple-frequency observations. Both reference and BDS receiver are connected with the same antenna with a signal splitter, to make sure they are processing the same signal. The specifics of IMUs used in the testing are listed in Table 1.Figure 3
The measurement vehicle for testing.Figure 4
NovAtel SPAN-ISA-100C system.Figure 5
The MEMS INS and BDS receiver for verifying the proposed algorithm. (a) The MEMS INS XW-5651. (b) The GNSS receiver Trimble BD990.Table 1
Specifications of IMUs.
ItemsISA-100CXW-5651GyroscopeAccelerometerGyroscopeAccelerometerBias0.5 deg/hr1250μg6 deg/hr0.1 mgScale factor100 ppm100 ppm1400 ppm1500 ppmRandom walk0.3 deg/√hr100μg/√hr0.3 deg/√hr0.25 mg/√hrDuring the vehicular field test, both the reference and developed systems are connected to continuous operating reference stations (CORS) network. Reference system using data was provided by CORS for high accuracy after positioning. The test system uses one of the base stations to form short-baseline double-difference observation. The offline process software Waypoint Inertial Explorer (IE) is used for processing high-accuracy positioning results. The test environment is chosen in a relatively open area with real traffic. In this test, the vehicle is driven at a normal traffic speed of 30–60 km/h.


## 3.2. Field Test

To test the positioning performance of the proposed algorithm, real observations of BDS are collected with help of the measurement vehicle and the reference stations. 2400 epochs of short-baseline differential data with triple-frequency BDS observations were involved in the calculation. The driving trajectory is shown in Figure6.Figure 6
Vehicle trajectory.Positioning error is obtained by the difference between process output and ground truth value. The corresponding visible satellite numbers are shown in Figure7. The errors of all epochs are plotted in Figure 8. The 95% empirical confidence which is shown in Figure 9 is presenting the reliability of the proposed algorithm. The skyplot for one of the epochs is presented in Figure 10. The accuracy represented with RMS statistics is given in Table 2.Figure 7
Observable satellite number.Figure 8
Positioning error of proposed TC algorithm based on triple frequency (m).Figure 9
The 95% empirical (red) confidence ellipses/intervals for (a) horizontal (N, E) position scatter and corresponding and (b) vertical (U) time series for the proposed model (green for inside and gray for outside).
(a)(b)Figure 10
Skyplot of one epoch.Table 2
Statistical positioning accuracy.
DirectionsRMSmE0.110N0.104U0.196Horizontal0.152


## 3.3. Signal Outage Test

A signal outage experiment is designed by artificially removing satellite signals in the experimental traffic section. Two lengths of common time interval used in inertial navigation experiments are designed in this experiment, where the signal outage durations of 10 s and 30 s were selected. The outages are designed as shown in Figure11.Figure 11
BDS signal outage arrangements.To verify the relationship between positioning results and visible satellite numbers, Figures12 and 13 show changes of visible satellites with triple-frequency observations and positioning error at the same time. The outage of visible satellite numbers in the figure is obtained by removing satellite observations of the corresponding epoch during data processing. The simulated blockage is easier to control the time length of outages. We focused on the positioning error curve around the outages. The positioning error and the corresponding satellite numbers are presented accordingly. The set of figures also includes the effect comparison of the proposed integration algorithm based on triple-frequency WL single-epoch observations versus the conventional algorithm. The conventional method refers to the method that the GNSS observation comes from the continuous measurement of multiple epochs. The float ambiguity is first estimated with KF and then fixed with LAMBDA.Figure 12
Positioning errors of 10 s outage compared with visible satellite numbers. (a) The corresponding visible satellite numbers. (b) The positioning error curve obtained with the proposed method. (c) The positioning error curve obtained with the conventional method.
(a)(b)(c)Figure 13
Positioning errors of 30 s outage compared with visible satellite numbers. (a) The corresponding visible satellite numbers. (b) The positioning error curve obtained with the proposed method. (c) The positioning error curve obtained with the conventional method.
(a)(b)(c)As presented in Figures12 and 13, the proposed algorithm instantly restores stable positioning output once the satellite signals are restored, while the conventional algorithm delays 4-5 epochs before returning to the stable output. Within the epochs, the divergence stopped under the effect of float ambiguity. But because the float ambiguity is uncontrollable, the positioning result is still not stable.We have counted the average statistical errors of three outages with the same length, as indicated in Table3 for 10 s outages. The 15 epochs before and after each outage were involved in statistics. In the table, the statistics of RMS, mean value, and maximum deviation range of the two methods are, respectively, carried out.Table 3
Average positioning error statistics during three 10 s signal outages.
ItemProcess methodEmNmUmRMSProp.0.2130.1720.422Conv.0.3090.2480.450MeanProp.−0.0820.0010.133Conv.−0.1640.0940.064Max. (abs)Prop.0.8830.6851.573Conv.0.9490.8651.289Prop. stands for the proposed method. Conv. stands for the conventional method.The average statistical errors of three 30 s outages are presented in Table4.Table 4
Average positioning error statistics during three 30 s signal outages.
ItemProcess methodEmNmUmRMSProp.1.7691.7784.283Conv.2.2542.1444.627MeanProp.−0.8710.7012.357Conv.−1.1300.8732.454Max. (abs)Prop.7.7087.09915.167Conv.8.8588.82016.179Prop. stands for the proposed method. Conv. stands for the conventional method.Compared with the conventional method, the positioning accuracy of the proposed method improves 43.6% horizontally and 6.2% vertically under the 10 s outages, with 27.5% in horizontal and 7.4% under 30 s outages (Table5).Table 5
The improvement of accuracy under different outage lengths.
Directions10 s (%)30 s (%)E31.121.5N30.717.1U6.27.4Horizontal43.627.5All44.128.5


## 4. Discussion

The first experiment is the functional verification of the proposed method. It can be clearly seen from the satellite number that the satellite is being frequently changed or blocked. It is obvious that the frequent satellite shifting does not cause obvious interference on the positioning outputs. The statistical results in the table demonstrate that the proposed method achieves a promising sub-decimetre-level positioning accuracy in an environment with frequent switches of observable satellites, which also meets the theoretical expectations brought by the ambiguity fixed WL observations. With help of the 95% confidence statistical information, it is proved that, with the proposed tightly coupled triple-frequency BDS/INS integration, a reliable and stable subdecimetre positioning accuracy is achieved when driving in a traffic environment with vehicular dynamic.The second experiment is about the verification of positioning recovery after signal interruption. The strong competitiveness of the proposed method is when operating in urban areas where signals are easily blocked. As shown in Figures12 and 13, the conventional method performs better with sufficient satellite observable because of the NL ambiguity fixed with the LAMBDA algorithm. However, when the satellite visibility recovers from blockage, because of the multiepoch ambiguity fixing requirement, a delay of several epochs appeared before stable output was obtainable. In those epochs, although the uncontrollable float solution prevented the positioning error from continuing to diverge, the positioning result was not stable. After the ambiguity is finally fixed, the stable positioning output is recovered. In contrast, the algorithm proposed in this paper fixes the ambiguity at the moment the satellite resumes observable. The stable positioning output is restored at the first epoch that the satellite observation restored, and therefore, relatively stable positioning performance is obtained under signal outages. Such an improvement has also been proved in the statistical tables.It is worth noting that although the proposed algorithm has achieved the expected performance, an important prerequisite for achieving a single-epoch ambiguity fixing is a short baseline (less than 10 km). For a longer baseline, the impact of atmospheric errors must be considered.


## 5. Conclusions

This paper contributed improved tightly coupled BDS/INS integrated positioning utilizing BDS’s triple-frequency observation characteristics. By introducing a single-epoch ambiguity resolution technique, a measurement modelling based on BDS triple frequency is proposed, which realizes an improvement of tightly coupled integration. This method enables BDS/INS integrated positioning algorithm to have both a simple scheme and high positioning accuracy in urban areas. As a result of the vehicular test in an actual traffic environment accompanied by frequent satellite varying, the positioning accuracy is stable and achieves 0.152 m for horizontal and 0.196 m for vertical positioning. The results prove that the proposed algorithm is reliable and meets the theoretical prospect of stable decimetre level positioning accuracy. The ability to restore stable positioning output in a single epoch was verified in the signal outage experiment. Compared with the conventional method, the positioning accuracy of the proposed method improves 43.6% horizontally and 6.2% vertically under the 10s outages, with 27.5% in horizontal and 7.4% under 30s outages. By ensuring reliable performance of tightly coupled BDS/INS integration, the research is expected to achieve the goal of lowering the hardware cost but with a better performance of the whole positioning system. In future work, we will focus on improving tightly coupled integration by introducing differential intersystem biases of GNSS to make the most use of limited available satellites.



---

*Source: 6653625-2021-01-28.xml*