Introduction

Moving measuring systems (MMSs) have been applied in surveying and city modeling (Viler et al. 2023). Initially, the core concept involves the integration of the global positioning system (GPS) and inertial navigation system (INS), providing precise centimeter-level data on position, velocity, and attitude (Ding et al. 2007). Since 2012, the China BeiDou Navigation Satellite System (BDS)-2 has further advanced the technique of integrating Global Navigation Satellite Systems (GNSS) and INS. Numerous studies have already explored the applications of BDS-2 in conjunction with INS (Han et al. 2017; Gao et al. 2018; Xiao et al. 2021). By 2020, BDS-3 began to provide global services (Zhu et al. 2023). Compared with BDS-2, BDS-3 has more advantages at the global service scale with the MEO constellation (Zhu et al. 2021) and precise point positioning (PPP)-B2b precise ephemeris augmentation (Tang et al. 2022). Therefore, BDS-3/INS has become a notable research topic (Ma et al. 2022; Xu et al. 2022).

The critical theoretical problems for BDS/INS tightly coupled (TC) include the measurement model and ambiguity resolution (AR) methods (Li et al. 2018), which determine the measurement usability, model robustness, and positioning accuracy. At some point, the measurement model and AR also undergo change. The model determines how accurately the measurements can be processed; thus, AR can be beneficial.

The BDS-3 visible satellites for standard positioning are sufficient for most districts and times (Li et al. 2020); however, reliable centimeter-level positioning puts forward more requirements on available measurements and model strength. Numerous studies have attempted to address this classical problem (Li et al. 2018), particularly in utilizing multi-frequency (Guo et al. 2016) and multi-GNSSs. Regarding multi-GNSSs, the terms of the measurement equation are constructed individually for different GNSSs, namely intra-system or loose combination (Zhang et al. 2003). If the frequencies overlap, only one reference satellite needs to be chosen for all GNSSs, and a tight combination can be carried out (Julien al. 2003). Although tight combination introduces new issues, such as intersystem bias, when compared with loose combination (Håkansson et al. 2017; Zhao et al. 2020), it is still considered robust (Tian et al. 2023).

Extensive investigations have been carried out concerning the model, methodology, and assessment of multi-GNSS. This includes studies on GPS/BDS-3 (Lv et al. 2022), BDS-3/Galileo (Guo et al. 2021), and BDS-3/GPS/GLONASS/Galileo (Ogutcu et al. 2023). However, the combination of BDS-3 and BDS-2, which share both similarities and differences in frequency signals, presents unique advantages and disadvantages and requires special attention. Further, investigation of BDS-3/BDS-2 benefits research on multi-GNSS. Several studies have made efforts to the simultaneous utilization of BDS-2 and BDS-3 in real-time kinematic (RTK) applications (Zhang et al. 2019; Zhu et al. 2021), with a primary focus on dual-frequency (DF), specifically the legacy B1I and B3I signals (Liu et al. 2021). BDS-2 and BDS-3 share the overlapping frequency of B1I, B2I(B2b), B3I, and the triple-frequency (TF) intersystem model offers significant advantage for BDS RTK in terms of measurement utilization and constellation structure; however, this aspect remains underexplored. Moreover, most studies have focused on investigating the signal characteristics; thus, an ideal environment, such as static station data (Zhang et al. 2019; Hou et al. 2019) and zero-baseline data (Deng et al. 2021), was usually used. However, the influence of the kinematic surroundings on signal quality, e.g., the pseudorange multi-path, has not been fully addressed. Further, the RTK model has been demonstrated with static station data (Shi et al. 2020), but the actual performance in city kinematic situations has not been comprehensively revealed. In addition, the BDS-3 could be enhanced with INS-providing robust position information, especially when satellite observations are disrupted, but this potential has been scarcely assessed in the context of BDS-3 (Wu et al. 2023; Ma et al. 2022). For city surveying, the mechanism of BDS/INS TC should be clarified with the advancement of BDS-3.

With the aim to promote the application of BDS/INS TC in kinematic surveying, and in light of the novel combination of BDS-3 and BDS-2, we developed the overlapping triple-frequency (TF) BDS-2/BDS-3/INS TC. This model’s performance has been thoroughly evaluated to address existing knowledge gaps and facilitate the adoption of optimal schemes for MMS.

The following sections present a theoretical exploration of ambiguity dilution of precision (ADOP) and a comprehensive demonstration of AR performance, including ADOP, Least-squares AMBiguity Decorrelation Adjustment (LAMBDA) ratio, and success rates. Furthermore, the positional and attitudinal determination performance of the proposed model was assessed.

TF BDS-3/BDS-2/INS TC model

Figure 1 presents the data flowchart of the proposed model. The state space model of the extended Kalman filter was implemented with two basic equations, namely the state equation and the measurement equation, expressed as follows:

$$\delta {\dot{\varvec{x}}} = \varvec{F}\delta \varvec{x} + \varvec{Gw}$$
(1)
$$\delta {\mathbf{L}} = {\mathbf{H}}\delta {\mathbf{x}} + {\mathbf{e}}$$
(2)

where F denotes the state transition matrix; G denotes the noise projection matrix; w denotes the system noise vector; \(\delta {\mathbf{L}}\) denotes the measurement vector, representing the bias between BDS-3 measurements and INS-predicted values; H denotes the design matrix; e denotes the noise vector; and \(\delta {\varvec{x}}\) denotes the states vector in the form of error quality.

Fig. 1
figure 1

TF BDS-3/BDS-2/INS TC model and data flowchart

Full size image

According to the error model adopted, states required, and measurements utilized, the details of these two equations can be obtained.

States equation

Equation (1) is constructed based on INS dynamic characteristics, where

$${\mathbf{F}} = \left[ {\begin{array}{*{20}c} {\mathbf{0}} & {\varvec{I}} & {\mathbf{0}} & {\mathbf{0}} & {\mathbf{0}} \\ {\mathbf{0}} & { - 2{\text{skew(}}{{\varvec{\upomega}}}_{{{\text{ie}}}}^{{\text{e}}} {)}} & {{\text{skew(}}{\varvec{R}}_{{\text{b}}}^{{\text{e}}} {\varvec{f}}_{{{\text{ib}}}}^{{\text{b}}} {)}} & { - {\varvec{R}}_{{\text{b}}}^{{\text{e}}} } & {\mathbf{0}} \\ {\mathbf{0}} & {\mathbf{0}} & { - {\text{skew(}}{{\varvec{\upomega}}}_{{{\text{ie}}}}^{{\text{e}}} {)}} & {\mathbf{0}} & {{\varvec{R}}_{{\text{b}}}^{{\text{e}}} } \\ {\mathbf{0}} & {\mathbf{0}} & {\mathbf{0}} & {\mathbf{0}} & {\mathbf{0}} \\ {\mathbf{0}} & {\mathbf{0}} & {\mathbf{0}} & {\mathbf{0}} & {\mathbf{0}} \\ \end{array} } \right],\quad {\mathbf{G}} = \left[ {\begin{array}{*{20}c} {\mathbf{0}} & {\mathbf{0}} & {\mathbf{0}} & {\mathbf{0}} \\ { - {\varvec{R}}_{{\text{b}}}^{{\text{e}}} } & {\mathbf{0}} & {\mathbf{0}} & {\mathbf{0}} \\ {\mathbf{0}} & { - {\mathbf{R}}_{{\text{b}}}^{{\text{e}}} } & {\mathbf{0}} & {\mathbf{0}} \\ {\mathbf{0}} & {\mathbf{0}} & {\varvec{I}} & {\mathbf{0}} \\ {\mathbf{0}} & {\mathbf{0}} & {\mathbf{0}} & {\varvec{I}} \\ \end{array} } \right]$$
(3)

where \({{\varvec{\upomega}}}_{{{\text{ie}}}}^{{\text{e}}}\) denotes earth’s rotation angular rate, \({\varvec{f}}_{{{\text{ib}}}}^{{\text{b}}}\) denotes the specific force output of the accelerometer, \({\varvec{R}}_{{\text{b}}}^{{\text{e}}}\) denotes the coordinate rotation matrix from the carrier body frame (b) to the earth-centered earth-fixed frame (e), and skew (V) represents the skew symmetry matrix formed by the three-dimensional vector V.

\(\delta {\mathbf{x}}_{{{\text{INS}}}} = [\delta {\mathbf{r}},\delta {\mathbf{v}},\delta {{\varvec{\uppsi}}},\delta {\mathbf{b}}_{a} ,\delta {\mathbf{b}}_{g} ]^{{\text{T}}}\), denotes the state vectors of position, velocity, and misalignment angle, accelerometer zero bias, and gyro zero bias, respectively; each vector is three-dimensional. The subscript of \(\delta {\mathbf{x}}\) denotes that these states are related to INS.

\({\mathbf{w}} = [{{\varvec{\upeta}}}_{a} ,{{\varvec{\upeta}}}_{g} ,{\mathbf{w}}_{a} ,{\mathbf{w}}_{g} ]^{{\text{T}}}\) represents the noise of the corresponding state vectors denoted by the subscript. The noise is modeled as a random walk process.

Measurement equation

The details of the BDS are summarized in Table 1. BDS-2 and BDS-3 possess overlapping frequencies, including B1I, B2b/B2I, and B3I, thereby achieving the tight combination. This allows for fixing observations to a single satellite, thus leading to an enhanced AR and improved positioning accuracy (Wu et al. 2019).

Table 1 Frequency and constellation for BDS-2 and BDS-3 (China Satellite Navigation Office 2019)
Full size table

BDS-2 and BDS-3 share the same time and coordinate frame, and the double-differenced (DD) linearized model of the observation equations at frequency k is the basic component of (2), expressed in the ECEF frame as follows:

$${\mathbf{L}}_{r,b,k}^{s,j} = \left[ {\begin{array}{*{20}c} {\Delta \varphi_{r,b,k}^{s,j} - \Delta \rho_{r,b,k}^{s,j} } \\ {\Delta P_{r,b,k}^{s,j} - \Delta \rho_{r,b,k}^{s,j} } \\ {\Delta D_{r,b,k}^{s,j} - \Delta v_{r,b}^{r,b} } \\ \end{array} } \right],\quad {\mathbf{H}}_{r,b,k}^{s,j} = \left[ {\begin{array}{*{20}c} { - {\mathbf{E}}_{r}^{s} } & {\mathbf{0}} & {\lambda_{k} } \\ { - {\mathbf{E}}_{r}^{s} } & {\mathbf{0}} & {\mathbf{0}} \\ {\mathbf{0}} & { - {\mathbf{E}}_{r}^{s} } & {\mathbf{0}} \\ \end{array} } \right],\quad {\mathbf{e}}_{r,b,k}^{s,j} = [\varepsilon_{\Delta \varphi } ,\varepsilon_{\Delta P} ,\varepsilon_{\Delta D} ]^{{\text{T}}}$$
(4)

where \(\Delta\) indicates that the following term is the between-station, between-satellite DD; the superscripts j and s denote the current and reference satellites, respectively; subscripts r and b denote the rover and base stations, respectively; subscript k denotes frequency; \(\delta {\mathbf{x}}_{{{\text{BDS}}}} = [\delta {\mathbf{r}}^{{\mathbf{e}}} ,\delta {\mathbf{v}}^{{\mathbf{e}}} ,\Delta \hat{N}_{r,b,k}^{s,j} ]^{{\text{T}}}\) denotes the state vector related to BDS measurements; \(\delta {\varvec{r}}^{e}\) and \(\delta {\mathbf{v}}^{{\mathbf{e}}}\) are as previously defined; \(N\) denotes the float ambiguity of the carrier phase;\(\lambda\) denotes the carrier phase wavelength; \(\varphi\), P, and D denote the carrier phase, pseudorange, and Doppler observation, respectively; \(\varepsilon\) denotes the measurement noise of the corresponding observations indicated by the subscript; E denotes the

line-of-sight vector of the receiver to the satellite; and \(\rho\) denotes the approximate geometric distance (GD) between the receiver and the satellite. In this model, ISBs are assumed to be absent, as demonstrated by the fact that differential-phase ISBs can be neglected in short baseline RTK (Shi et al. 2022). Additionally, if the satellite receivers on the rover and base station are the same type, the ISBs for the code and phase can be neglected (Shi et al. 2022; Hou et al. 2023).

By combining the observation equations of all satellites and frequencies, the conventional BDS positioning measurement equation is obtained. This serves as the base for augmenting the absent states related to \(\delta {\mathbf{x}}_{{{\text{INS}}}}\). Consequently, the general BDS/INS TC measurement equation is derived. In this study, INS was used to assist positioning, then \(\Delta \rho_{r,b}^{j,s}\) was replaced with \(\Delta \rho_{{{\text{INS}},b}}^{j,s}\), and \(\Delta v_{r,b}^{r,b}\) was replaced with \(\Delta v_{{{\text{INS}},b}}^{j,s}\). The subscript INS denotes the terms are INS-predicted. Additionally, the following expression should be added to (2) to obtain the INS-assisted measurement equation:

$${\mathbf{0}}_{3 \times 1} = {\mathbf{I}}_{3 \times 3} \cdot \delta {\mathbf{r}}_{{{\text{INS}}}}^{e} + \varepsilon_{{{\text{INS}}}}$$
(5)

where I denotes unit matrix.

An elevation-dependent weighting function (Herring et al. 2018) was applied to determine the noise of the pseudorange and phase as follows:

$$\sigma_{{}}^{2} {(}\theta {) = }a^{2} + b^{2} /\sin (\theta )^{2}$$
(6)

where \(\theta\) denotes the elevation angle; a and b denote empirical values, both set to 0.003 m for the phase, and multiply by 100 for the pseudorange (Takasu 2013), respectively.

The measurement noise for the \({\mathbf{r}}_{{{\text{INS}}}}^{e}\) is as follows:

$$\sigma_{{{\mathbf{r}}_{{{\text{INS}}}}^{{\mathbf{e}}} }} = br_{{{\text{prior}}}}^{{\mathbf{e}}} + \upsilon \cdot \nabla t$$
(7)

where \({\text{b}}r_{{{\text{prior}}}}^{{\mathbf{e}}}\) varies according to the position bias of the last epoch, which is set to 0.05 m if it is a fixed solution, otherwise to 0.5 m; \(\upsilon\) denotes the INS error accumulating speed which is set to 0.005 m/s; and \(\nabla t\) denotes the interval between current and last epoch.

Assessment of AR accuracy

To obtain centimeter-level accuracy, the integer-valued ambiguity should first be determined. The classical approach is LAMBDA, proposed by Teunissen (1993), which has been widely applied and developed. LAMBDA solves integer ambiguity through a search process within a constraint space, defined as

$$(\Delta \hat{N} - \Delta N){\mathbf{Q}}_{{\Delta \hat{N}}}^{ - 1} (\Delta \hat{N} - \Delta N) \le \chi^{2}$$
(8)

where \(\Delta \hat{N}\) denotes the float solutions of the ambiguities, \(\Delta N\) denotes the real value of the integer ambiguity, and \({\mathbf{Q}}_{{\Delta \hat{N}}}^{{}}\) denotes the covariance matrix of \(\Delta \hat{N}\). The search space is centered on the \(\Delta \hat{N}\), the shape and orientation are decided by \({\mathbf{Q}}_{{\Delta \hat{N}}}^{{}}\), and the space volume is limited by \(\chi^{2}\). The strength and success rate of the GNSS model are generally measured using the ADOP (Odijk and Teunissen 2008).

ADOP of various frequencies and satellites

ADOP represents the average precision of ambiguities proposed by the LAMBDA proposers and is defined as follows (Teunissen and Odijk 1997):

$${\text{ADOP}} = \left| {{\mathbf{Q}}_{{\Delta \hat{N}}} } \right|^{{\frac{1}{{{\text{2N}}_{{{\text{Amb}}}} }}}} \begin{array}{*{20}c} {} \\ \end{array} ({\text{cycles}})$$
(9)

where | | denotes the determinant, and \({\text{N}}_{{{\text{AMB}}}}\) denotes the number of ambiguities, which is equal to the degree of \({\mathbf{Q}}_{{\Delta \hat{N}}}^{{}}\).

The ADOP has been comprehensively investigated and used to assess AR (Skaloud 1998; Liu et al. 2019). Equation (9) cannot qualitatively measure the float ambiguity accuracy; hence, Teunissen (1997a) developed an ADOP closed-form expression for various GNSS models. However, a TF GNSS model has not yet been formulated. In order to thoroughly investigate how factors such as the number of satellites, observation accuracy, and frequency combinations impact our model, we have deduced the \(\left| {{\mathbf{Q}}_{\Delta N} } \right|\) expression for TF, which is provided in Online Resource 1, as follows:

$$\left| {{\mathbf{Q}}_{\Delta N} } \right| = m^{y} \left( {\prod\limits_{i = 1}^{y} {\frac{{\sigma_{{\Delta \varphi_{i} }}^{2} }}{{\lambda_{i}^{2} }}} } \right)^{m - 1} \left( {1 + \frac{{\sum\nolimits_{i = 1}^{y} {1/\sigma_{{\Delta \varphi_{i} }}^{2} } }}{{\sum\nolimits_{i = 1}^{y} {1/\sigma_{{\Delta p_{i} }}^{2} } }}} \right)^{3} \quad ({\text{cyc}}^{2(m - 1)y} )$$
(10)

where y denotes the number of frequencies, and m denotes the number of satellites, thus indicating the numbers of observations and ambiguities on one frequency are m − 1. \(\sigma_{{\Delta \varphi_{i} }}^{2}\) and \(\sigma_{{\Delta p_{i} }}^{2}\) denote the variances of the DD carrier phase and pseudorange, respectively, and \(\sigma_{{\Delta \varphi_{i} }}^{2} { = }2\sigma_{{\varphi_{i} }}^{2}\) and \(\sigma_{\Delta Pi}^{2} { = }2\sigma_{Pi}^{2}\). The observation precision at all frequencies was assumed to be identical as follows:

$$\sigma_{\Delta pi}^{2} = \sigma_{\Delta p}^{2} ,\quad i = 1,2,3$$
(11)
$$\sigma_{{\Delta \varphi_{i} }}^{2} /\lambda_{i}^{2} = \sigma_{\Delta \varphi }^{2} /\lambda^{2} ,\quad i = 1,2,3$$
(12)

Then, according to (10):

$$\left| {{\mathbf{Q}}_{\Delta N} } \right| = m^{y} \left( {\frac{{\sigma_{\Delta \varphi }^{2} }}{{\lambda^{2} }}} \right)^{(m - 1)y} \left( {1 + \frac{{\sigma_{\Delta p}^{2} }}{{\sigma_{\Delta \varphi }^{2} }}} \right)^{3} \quad ({\text{cyc}}^{2(m - 1)y} )$$
(13)

The ADOP is expressed as

$${\text{ADOP}} = \left| {{\mathbf{Q}}_{{\mathbf{N}}} } \right|^{{\frac{1}{2y(m - 1)}}} = m^{{\frac{1}{2(m - 1)}}} \left( {\frac{{\sigma_{\Delta \varphi } }}{\lambda }} \right)\left( {1 + \frac{{\sigma_{\Delta p}^{2} }}{{\sigma_{\Delta \varphi }^{2} }}} \right)^{{\frac{3}{2y(m - 1)}}} \quad ({\text{cyc}})$$
(14)

It also coincides with that of Odijk and Teunissen (2008) if the specified parameters are set according to our model.

ADOP of INS assistance

Several studies have reported ADOP improvements when using INS to assist the GNSS model (Han et al. 2017). A pseudorange with large errors is replaced with an INS-predicted high-accuracy GD; thus, the float solutions of the ambiguities can be beneficial. An expression similar to (14) can then be inferred.

$${\text{ADOP}}_{{{\text{INS}}}} = \left| {{\mathbf{Q}}_{{\mathbf{N}}} } \right|^{{\frac{1}{2y(m - 1)}}} = m^{{\frac{1}{2(m - 1)}}} \left( {\frac{{\sigma_{\Delta \varphi } }}{\lambda }} \right)\left( {1 + \frac{{\sigma_{{\Delta {\text{INS}}}}^{2} }}{{\sigma_{\Delta \varphi }^{2} }}} \right)^{{\frac{3}{2y(m - 1)}}} \quad ({\text{cyc}})$$
(15)

Then, the scale factor of improved ADOPINS versus original ADOP is

$${\text{S}}_{{{\text{ADOP}}}} = \left[ {{{\left( {1 + \frac{{\sigma_{{\Delta {\text{INS}}}}^{2} }}{{\sigma_{\Delta \varphi }^{2} }}} \right)} \mathord{\left/ {\vphantom {{\left( {1 + \frac{{\sigma_{{\Delta {\text{INS}}}}^{2} }}{{\sigma_{\Delta \varphi }^{2} }}} \right)} {\left( {1 + \frac{{\sigma_{\Delta p}^{2} }}{{\sigma_{\Delta \varphi }^{2} }}} \right)}}} \right. \kern-0pt} {\left( {1 + \frac{{\sigma_{\Delta p}^{2} }}{{\sigma_{\Delta \varphi }^{2} }}} \right)}}} \right]^{{\frac{3}{2y(m - 1)}}}$$
(16)

Assume \(\sigma_{{\Delta {\text{INS}}}}^{2} \gg \sigma_{\Delta \varphi }^{2}\), \(\sigma_{\Delta P}^{2} \gg \sigma_{\Delta \varphi }^{2}\), then \({\text{S}}_{{{\text{ADOP}}}} \approx (\frac{{\sigma_{{\Delta {\text{INS}}}}^{2} }}{{\sigma_{\Delta p}^{2} }})^{{\frac{3}{2y(m - 1)}}}\). As all the parameters \({\text{S}}_{{{\text{ADOP}}}}\) are definite, the ADOP value for each ambiguity is scaled down when introducing INS.

Assuming that \(\sigma_{\Delta \varphi }\) = 0.008 \(\sqrt 2\) m and \(\sigma_{\Delta P}^{{}}\) = 0.3 \(\sqrt 2\) m, these values were determined based on their DD residuals shown in Figs. 6 and 8. For tactical-grade IMU used in MMS, set \(\sigma_{{\Delta {\text{INS}}}}^{{}}\) to 0.05 m to draw ADOP (Fig. 2).

Fig. 2
figure 2

ADOP versus number of satellites for SF, DF, TF, and INS assistance. The lower right part of the figure is enlarged in the insert on the upper right

Full size image

The single-frequency (SF) ADOP is much larger than those of DF and TF, which is severe when the number of satellites is < 8, indicating that if SF was applied to precise positioning, a sufficient number of tracked satellites should be maintained. Combining BDS increases the number of satellites, which can be an applicable approach. Additionally, as the number of satellites increases, there is a noticeable trend of the ADOPs slowing down, suggesting limited improvements in AR with DF and TF. This phenomenon also implies that a TF-only model was chosen to construct a highly precise measurement model. Furthermore, in the case of a loose combination between BDS-3 and BDS-2, the number of measurements was one less than in the tight combination, resulting in a larger ADOP. This illustrates the advantages of opting for the tight combination.

Vehicle examination and results

Vehicle tests were conducted to evaluate the proposed model from 11:15:55 to 12:19:49 (local time) in 2022. The experiment was conducted around a large lake, strengthening the multi-path effect, which provides challenging conditions, which is plotted in Fig. 3. The BDS signals were received using SinoGNSS M66 and Novatel GNSS-850 antennas, both at the rover and base stations (presented in Fig. 4). The IMU data (sampling rate of 100 Hz) were output from a StarNeto xw-7960, with an accelerometer zero-bias instability of 20 mg and a gyroscope zero-bias instability of 0.1°/h. The random walk process noises of the accelerometer and gyro were 1.95, 1.86, and 15.2 \(ug/\sqrt {{\text{Hz}}}\) and 0.0027, 0.0025, and 0.0046 \(^\circ /\sqrt h\), respectively, as estimated using the ALLAN variance method with 6 h static data. The BDS/INS integrated processing strategies are listed in Table 2.

Fig. 3
figure 3

Vehicle trajectory illustrated as a green line in the financial center in Zhengzhou city, China. The base station and baseline are illustrated in light blue

Full size image
Fig. 4
figure 4

Base station and vehicle are both equipped with the same type SinoGNSS M66 receiver and the same antenna of Novatel GNSS-850

Full size image
Table 2 Integrated data processing strategy
Full size table

The satellite elevation angles are illustrated in Fig. 5. C38 is one of the highest satellites, and its elevation remains high during the test; thus, it was used as the reference satellite. Although the visible satellites for BDS-2 and BDS-3 can reach 10 and 11, C04, C05, C12, C24, and C45 for the entire timespan and C20, C22, C26, and C29 for most epochs were low in elevation angles, thereby lowering signal quality for the city environment. Furthermore, low-elevation satellites have higher atmospheric delay errors and poor observation quality, which affect AR. Therefore, the elevation cutoff angle was set to 30°; the available satellites are sufficient when abandoning satellites whose elevation angles are below 30°.

Fig. 5
figure 5

Satellite elevation angles observed by the vehicle receiver. Bars indicate the variation range during the test

Full size image

The number of visible BDS-2 satellites over the Asia–Pacific region is 8–12 (Li et al. 2020), and the number of BDS-3 satellites over most districts worldwide is 8–14. The joint utilization of BDS-2 and BDS-3 is beneficial for increasing the number of usable satellites and their geometric distributions.

The DD residual pseudorange errors are shown in Fig. 6. For example, for B2b, most points were within a certain range during the test time, and a few so-called gross errors exceeded the regular extent in some epochs, very likely owing to the multi-path fact (Bona 2000). This phenomenon is clearly shown for other frequencies. These statistics are presented in Table 3. Moreover, these gross errors are much more severe than those in static station cases (Zhu et al. 2023) in terms of both dispersion and extent, demonstrating that the kinematic situation should be specifically considered. In addition, the smaller bar height of BDS-3 indicates that BDS-3 is more reliable than BDS-2 in pseudorange accuracy, as previously indicated by Deng et al. (2021), owing to the effective resistance to the multi-path effect.

Fig. 6
figure 6

DD residual pseudorange errors drawn as time series curves, BDS-3 B2b, as an example (top), and the overview (bottom) of B1, B2, and B3 for BDS-2 and BDS-3 are shown. The bars denote the normal range, while the line in the center marks the mean. The endpoints of the vertical lines indicate the maximum and minimum values, respectively

Full size image
Table 3 RMSs for BDS and frequencies (m)
Full size table

The accuracy of the pseudorange is a critical factor, especially since the float solutions of ambiguity in the single-epoch AR method are determined by the pseudorange. Generally, combined ambiguities are formed to mitigate the effects of large pseudorange errors and perform sequential fixing until the narrow lane (NL) ambiguity is fixed. Nonetheless, when the INS is introduced, the disadvantages of float solution determination for original ambiguities can be avoided.

Figure 7 illustrates the ionospheric delay calculated using the B1I/B3I dual-frequency phases. The ionosphere delay introduced considerable errors in the carrier phase accuracy and should thus be reduced. This category of errors is mainly eliminated by DD processing in the model.

Fig. 7
figure 7

Original ionosphere delay at B1I of BDS-2 (top) and BDS-3 (bottom) retrieved from B1I/B3I dual-frequency measurements

Full size image

As shown in Fig. 8, the DD residual carrier phase biases observed for BDS-2 and BDS-3 do not exhibit significant differences. The average RMSs for biases across the three frequencies are tabulated in Table 3, revealing that a few epochs surpass 0.04 m, which can be regarded as the most epochs that can support effective AR. However, it is noteworthy that the biases for BDS-3 appear to be less reliable. This is evident from the greater number of data points exceeding 0.02 m when compared to BDS-2. Moreover, the RMS values for BDS-3 are slightly larger than those for BDS-2. This discrepancy may likely stem from variations in receiver performance.

Fig. 8
figure 8

DD residual carrier phase bias of B1 (left), B2 (middle), and B3 (right) of BDS-2 (top) and BDS-3 (bottom)

Full size image

The DD residual bias is much smaller than the original ionosphere delay shown in Fig. 6, implying that the ionosphere errors, phase ISBs, and unmodeled errors have been effectively corrected to a large extent; therefore, AR is free of their impacts. It also demonstrated the suitability of our model for a short-medium baseline.

Approaches have been developed for BDS to solve AR (Zhang et al. 2020). Due to the frequent interference of satellite signals in urban districts, which causes cycle slips and irregularly missing observations in frequencies and satellites, single-epoch LAMBDA could be the preferred reliable method. For conciseness, BDS-3/BDS-2 is expressed as BDS-C.

ADOP

Figure 9 illustrates the average ADOP across all epochs, with each epoch’s ADOP computed based on the actual Qn values for mutual confirmation with the theoretical characteristic revealed by (14).

Fig. 9
figure 9

Average ADOPs across all epochs for various conditions indicated by labels on the x-axis. Vertical lines indicate 3 \(\sigma\) ranges

Full size image

The ambiguity can be successfully determined if ADOP < 0.2 cycles (Teunissen 1997b); however, if ADOP is below 0.12 cycles, the possibility of achievable AR success can reach 99.9% for any number of ambiguities (Odijk and Teunissen 2008). Except for SF BDS-2 and BDS-3, whose ADOPs were 0.423 and 0.503 cycles, respectively, the ADOPs were below 0.158 cycles, indicating that for general conditions of satellites and frequencies, the effective AR could be satisfied for most epochs. Additionally, an increase in both frequencies and satellites can decrease ADOPs. However, if the achievable ADOPs are sufficiently small, further improvements are insignificant. Overall, the basic law is the same as shown in Fig. 2. In addition, BDS-2 outperforms BDS-3 because the number of available satellites is equal to or greater in most epochs, as shown in Fig. 10 (top panel). In addition, BDS-3 C60 did not broadcast B2b signals in TF situations.

Fig. 10
figure 10

Time series of LAMBDA ratios (top panel left y-axis) and number of available satellites (top panel right y-axis) for BDS-2 (left), BDS-3 (middle), and BDS-C (right) with SF, DF, and TF. The bottom histogram denotes the share (y-axis) of certain ratios (x-axis) in %

Full size image

LAMBDA ratio

The LAMBDA ratio is calculated as the square norm of the second-best candidate divided by the square norm of the best candidate in the integer least-squares search (ILS). The ambiguities are considered to be correctly resolved if the ratio is larger than the critical value in the ratio test. In general, the critical value is set to 3 in static cases (Deng et al. 2021) and 2 in RTK (Zhang et al. 2019; Liu et al. 2021). In this study, the critical value is set to 2.5.

Figure 10 illustrates the time series and distribution of LAMBDA ratios for SF, DF, TF, and BDS-2, BDS-3, and BDS-C. Note that BDS-C consistently achieves higher ratios in the majority of epochs compared to individual BDS. The instances where a single BDS exhibits a lower ratio can be offset by the contributions from other BDS.

Regarding the impact of frequencies on the ratios, DF and TF significantly outperform SF. However, the distinction between the DF and TF was not evident. In particular, the performance of TF is even weaker than that of DF for BDS-C at some epochs, possibly because DF already provides sufficient observations and adding a third frequency provides limited advantages; however, low-accuracy observations can decrease AR. Additionally, more ambiguities smooth the difference between the first two candidates of ILS (Julien et al. 2003). Since a ratio larger than critical value is sufficient for determining whether to accept the integer ambiguities, having a higher ratio primarily offers theoretical advantages with minimal practical improvement.

Success rate

The success rate is determined by dividing the number of epochs with correctly resolved ambiguities by the total number of epochs. Consequently, the success rate represents the availability of high-accuracy positioning (Liu et al 2019). Figure 11 illustrates that SF success rates were 13.65 and 12.07% for BDS-2 and BDS-3, respectively, which are far from acceptable. When INS aiding was adopted, the SF BDS could be improved by approximately 12.77 and 7.12%, respectively, indicating that very low rates limit the contribution of INS. However, SF BDS-C achieved 69.29%, demonstrating that the joint combination of BDS strongly contributes to SF AR. Furthermore, BDS-C improved the DF and TF by at least 22.75 and 20%, respectively.

Fig. 11
figure 11

Success rates for various BDS and frequencies. The total epoch was 3819

Full size image

Regarding the impact of the number of frequencies, it is noteworthy that for individual BDS systems, there is a significant improvement in the rates when transitioning from SF to DF; however, the improvement becomes negligible when the additional third frequency was introduced. Furthermore, it is important to recognize that the DF performance is closely tied to the specific frequency selection, e.g., results may vary if B1I/B2I is chosen.

Concerning INS assistance, various combinations of BDS and frequencies have the potential to be beneficial. However, the improvements observed are generally modest, typically falling within the range of 1–5%.

It can be concluded that increasing the number of satellite numbers, frequency, and INS assistance would improve AR performance to some extent. The increase in the number of satellites is remarkable; the frequency increase is only evident from SF to DF. While INS assistance positively affects all combinations, including BDS-2/INS, BDS-3/INS, and BDS-C/INS, its overall contribution is quantitatively limited.

Integrated positioning and attitude determination

The position and attitude angles are the two most significant terms for the MMS and were therefore analyzed to demonstrate the effectiveness of the BDS/INS TC. It is practically impossible to establish absolute precision in kinematic situations; therefore, the Waypoint Inertial Explorer (IE) and RTKLIB softwares were employed to demonstrate the relative data processing accuracy. IE has been extensively utilized in GNSS/INS research and can provide position, velocity, and attitude angle as reference. However, IE does not support BDS-only processing (GPS-relied), and it utilizes overall information to achieve higher AR success rates. Thus RTKLIB v2.4.3 (GNSS processing software) was introduced to provide position result for fair comparison. The results are shown in Fig. 12 and Table 4.

Fig. 12
figure 12

Time series of bias for position (top, middle) and attitude angle (bottom); the reference software is denoted in the top left corner of each subfigure. The purple points indicate the belonging epochs failed in AR. Only the epochs with successfully resolved ambiguities from the reference result are compared

Full size image
Table 4 RMS of the bias between our result and those obtained from IE and RTKLIB
Full size table

As indicated, our program achieved a comparable result to RTKLIB, with a millimeter-level position bias. When compared to the IE’s results, our positioning bias reached centimeter levels of 0.030, 0.042, and 0.065 m. However, it is worth noting that the accuracy was reduced due to significant errors in epochs where AR failed. Given the variation in the signal quality of different GNSS and the undisclosed optimization strategy employed by IE, the observed degree of bias is considered reasonable.

The heading angle is the most significant component of the attitude angle, but the most difficult term for the filter to estimate in the ground-moving situation. Thus, the heading-angle bias was larger than those of the other two angles and required a considerably longer time to converge after the initial fluctuation. However, the heading angle achieved an RMS value of 0.066°. In summary, the processing biases of position, velocity, and attitude angles are admissible, demonstrating the effectiveness of the BDS/INS TC model.

The states uncertainty, comprising the diagonal elements of the variance matrix of states, indicates the precision of estimated states (Groves 2013). Figure 13 shows the uncertainty series for position and attitude. While various jump points occur on the position curve, it is worth noting that their magnitudes are very small. This was mainly due to the failure of the AR. In contrast, the attitude uncertainty decreased rapidly during the initial period, indicating that the attitude angle converged and that the BDS/INS TC filter functioned properly.

Fig. 13
figure 13

States uncertainty for position (top) and attitude angle (bottom)

Full size image

Conclusions

To promote the BDS/INS technique with the advancement of BDS-3, we developed and assessed a TF BDS-3/BDS-2/INS TC model. The proposed model benefited from the use of the TF overlapping frequency and INS, enhancing reliability. The ADOP closed-form expression was derived to investigate various factors affecting AR, providing theoretical cross-verification of the experimental results.

An urban vehicle test was conducted. In contrast to the results obtained in static situations, the accuracy of observations and the performance of AR showed lower performance levels. This suggests the need to develop a reliable model tailored to urban environments rather than simply applying a stationary-verified RTK model.

The integration of INS can enhance the GNSS model and improve its overall reliability, although the extent of improvement depends on specific conditions. Notably, INS assistance proves most valuable when faced with challenges such as insufficient satellites and single frequency.

Furthermore, the combined use of multiple BDS systems demonstrates a substantial improvement compared to utilizing individual BDS systems. This improvement can be attributed to the increased number of tracked satellites and their structure.

The utilization of multi-frequency signals improves AR. However, it is important to note that the improvement may not counteract the negative impact because, in urban districts, as the number of signals at the frequency increases, it is more difficult to ensure every signal’s accuracy. In general, DF can be a balanced scheme.

The findings of this study demonstrate the TF BDS contribution to the future application of BDS/INS. We will extend the model presented in this study to multi-GNSS and further demonstrate the ADOP closed-form expression with multi-GNSS tests in subsequent research.