Introduction

The Quasi-Zenith Satellite System (QZSS) is a Japanese satellite positioning system operated by Japan Aerospace Exploration Agency (JAXA). With the orbits specially designed to increase high-elevation signals in Japan, QZSS has been developed as a GPS-complementary system to increase the availability, reliability, integrity and accuracy of positioning in the Asia-Oceania region, especially in urban canyons and mountainous areas (QZSS 2018; Wu et al. 2004a). The first QZSS satellite (QZS-1) is in orbit since 2010. With the second (QZS-2), the third (QZS-3) and the fourth QZSS satellites (QZS-4) launched, respectively, in June, August and October of 2017, QZSS has now completed a four-satellite constellation. QZSS is planned to be extended to a seven-satellite system by 2023 (Murai 2015). Except for QZS-3, which is a geostationary (GEO) satellite, the other three QZSS satellites are in quasi-zenith orbits (QZO) with an inclination of around 43 ± 4°. Having similar figure-eight shaped ground tracks, the QZO satellites spend an extended time with high elevations above Japan with a separation of 8 h passing the same region (Teunissen and Montenbruck 2017). Figure 1 [left] shows the colormap indicating the percentage of a 24-h period with at-least one QZSS satellite (out of QZS-1, -2, -3 and -4) visible above 70° of elevation on 2 February 2018. We note that this percentage indeed reaches 100% in Japan. The middle panel in Fig. 1 illustrates the ground tracks of the four QZSS satellites where J1, J2, J3 and J7 are, respectively, the PRN numbers of QZS-1, -2, -4 and -3. The colormap shown in Fig. 1 [right] indicates the percentage of a 24-h period on 2 February 2018 with all the four QZSS satellites being visible simultaneously where a cut-off elevation angle of 10° is considered. It can be seen that there is an eye-shaped area in Asia-Oceania region, covering a large part of Australia, where all the four QZSS satellites always lie above 10° of elevation.

Fig. 1
figure 1

QZSS satellites visibility. [Left] Colormap indicating the percentage of a 24-h period with at-least one QZSS satellite (out of J1, J2, J3, J7) visible above 70° of elevation. [Middle] Ground tracks of the QZSS satellites J1, J2, J3, J7. [Right] Colormap indicating the percentage of a 24-h period with four QZSS satellites visible simultaneously considering a cut-off elevation angle of 10°. All the panels are obtained based on the QZSS constellation on 2 February 2018

Full size image

The QZSS satellites send GPS-compatible signals on L1 (1575.42 MHz), L2 (1227.6 MHz) and L5 (1176.45 MHz) and augmentation signals, i.e. the L1-SAIF (Submeter-Class Augmentation with Integrity Function) signal and the LEX (L-band Experimental) signal. Signal analyses were performed by different studies after the first QZSS satellite was launched. Using the signals from QZS-1, Nie et al. (2015) observed similar multipath/noise-combined behaviors for QZSS and GPS on each of the three frequencies after forming the multipath combination (Leick et al. 2015), while Hauschild et al. (2012) reported lower standard deviations of the noise and multipath errors on QZSS L5 than those on GPS L5. In the zenith direction, both of them amounted to around 5 cm. Based on zero-baseline tests, the standard deviations of the double-differenced carrier-phase residuals were reported to be between 0.5 and 1 mm for QZSS L1, L2 and L5 (Quan et al. 2016), and Nadarajah et al. (2016) showed that the QZSS L5 code observations have smaller zenith-referenced standard deviation compared to the QZSS L1 code signals. A recent study by Steigenberger et al. (2018) investigates the QZSS signal power and clock performance.

Having reached a four-satellite constellation just recently, QZSS positioning performance has hitherto been studied only in combination with the other Global Navigation Satellite Systems (GNSSs). Comparative studies on GPS-only versus GPS + QZSS positioning can be found in Wu et al. (2004b) and Choi et al. (2015). Nadarajah et al. (2016) assessed multi-GNSS single-frequency real-time kinematic (RTK) positioning and attitude determination performance through combining QZSS L1/L5 signals with their counterparts from GPS, Galileo and IRNSS. RTK results of QZSS in combination with GPS, Galileo and Beidou were presented in, e.g. Odolinski et al. (2015), Quan et al. (2016) and Odolinski and Teunissen (2017).

With four operational QZSS satellites available in orbit, it is now possible to carry out positioning based on standalone QZSS. In this contribution, using the signals from all the four QZSS satellites on the three frequencies L1, L2 and L5, a first assessment of standalone QZSS ambiguity resolution performance and positioning accuracy is presented. After introducing the processing strategy, analysis of signal precision is performed for L1, L2 and L5. The impact of multipath on both the code and phase observations is discussed and the zenith-referenced variance matrices are presented before and after multipath mitigation. Next, considering a short-baseline RTK model, the single-epoch ambiguity resolution performance is investigated based on the QZSS triple-frequency observations, through ambiguity dilution of precision (ADOP) and integer bootstrapped success rate. It is thereby shown that instantaneous successful ambiguity resolution is feasible using both the multipath-uncorrected and -mitigated observations. Afterward, the corresponding positioning accuracies are analyzed and compared for ambiguity-float and -fixed cases under single-epoch RTK model without and with multipath mitigation. We further extend our evaluations to a longer baseline for which the differential atmospheric delays are not negligible. A summary of conclusions is provided at the end.

Processing strategy

For baselines within a few kilometers, it is assumed that the ionospheric delays are significantly reduced by forming between-receiver differences of the GNSS signals (Goad 1998). The double-differenced (DD) observed-minus-computed code and phase observations transmitted by m satellites on f frequencies to a pair of receivers at a single epoch can be formulated as (Teunissen and Montenbruck 2017)

$$E\left[ {\begin{array}{*{20}{c}} p \\ \phi \end{array}} \right]=\left[ {\begin{array}{*{20}{c}} {{e_{f}} \otimes {I_{m - 1}}}&0 \\ {{e_{f}} \otimes {I_{m - 1}}}&{\Lambda \otimes {I_{m - 1}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\widetilde {\rho }} \\ a \end{array}} \right]$$
(1)

where \({E}\left[ \cdot \right]\) represents the expectation operator, and \(\otimes\) denotes the Kronecker product (Henderson et al. 1983). The DD code and phase observations are, respectively, represented by \(p={\left[ {p_{1}^{{T}}, \ldots ,p_{{f}}^{{T}}} \right]^{T}}\) and \(\phi ={\left[ {\phi _{1}^{{T}}, \ldots ,\phi _{{f}}^{{T}}} \right]^{T}}\), where \({p_j}\) and \({\phi _j}\,\left( {j=1, \ldots ,f} \right)\) contain the \((m - 1)\) DD code and phase observations on frequency \(j\). The symbols \({e_f}\) and \({I_{m - 1}}\) stand for a vector of ones with the length of \(f\) and identity matrix with the size of \((m - 1)\), respectively. The \(f \times f\) diagonal matrix \(\Lambda\) has as its diagonal entries \({\lambda _j}\), the wavelength of frequency \(j\,(j=1, \ldots ,f)\). The \((m - 1)\)-vector \(\widetilde {\rho }\) captures the non-dispersive parameters including the DD geometric ranges \(\rho\) and DD slant tropospheric delays \(\tau\), and \(a\) denotes the \(f(m - 1)\)-vector of DD ambiguities. The dry part of the slant tropospheric delays can a priori be corrected using Saastamoinen model (Saastamoinen 1972), and the remaining residual tropospheric delays (wet part) \(\Delta \tau\) can be parametrized in the scalar Zenith Tropospheric Delay (ZTD) \(\Delta {\tau _{\text{Z}}}\). One can further linearize the link between the DD geometric ranges and the baseline components to achieve

$$\Delta \widetilde {\rho }=\Delta \rho +\Delta \tau =~D_{{m}}^{{T}}\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} { - {u^{1^{T}}}}&{{g^1}} \end{array}} \\ {\begin{array}{*{20}{c}} { - {u^{2^{T}}}}&{{g^2}} \end{array}} \\ \vdots \\ {\begin{array}{*{20}{c}} { - {u^{m^{T}}}}&{{g^m}} \end{array}} \end{array}} \right]~~\left[ {\begin{array}{*{20}{c}} {\Delta x} \\ {\Delta {\tau _{\text{Z}}}} \end{array}} \right]$$
(2)

with \(\Delta \rho\) the DD geometric range increments, \(\Delta x\) the 3-vector of baseline increments, \({u^s}\,(s=1, \ldots ,m)\) the receiver-to-satellite unit vectors, \({g^s}\,(s=1, \ldots ,m)\) the satellite elevation-dependent mapping function which maps the slant tropospheric delays from all satellites to ZTD (Ifadis 1986), and \(D_{{\text{m}}}^{{\text{T}}}=\left[ { - {e_{m - 1}},~~{I_{m - 1}}} \right]~\) is the differencing operator forming the between-satellite differences. Note in (2), the elevation angle from both receivers of the baseline to a given satellite is assumed to be the same. We also remark that the differential ZTD in (2) can be assumed negligible for short baselines, i.e. \(\Delta {\tau _{\text{Z}}} \approx 0\).

The variance matrix of the DD observations in (1) can be formed as

$$D\left[ {\begin{array}{*{20}{c}} p \\ \phi \end{array}} \right]=2 \times \left[ {\begin{array}{*{20}{c}} {\left( {{I_{f}} \otimes D_{{m}}^{{T}}} \right){Q_{{pp}}}\left( {{I_{f}} \otimes {D_{m}}} \right)}&0 \\ 0&{\left( {{I_{f}} \otimes D_{{m}}^{{T}}} \right){Q_{\phi \phi }}\left( {{I_{f}} \otimes {D_{m}}} \right)} \end{array}} \right]$$
(3)

where \({D}\left[ \cdot \right]\) denotes the dispersion operator, and the variance matrices of the undifferenced code and phase observations \({Q_{pp}}\) and \({Q_{\phi \phi }}\) read

$${Q_{{\text{pp}}}}={C_{{pp}}} \otimes {W^{ - 1}}$$
(4)
$${Q_{\phi \phi }}={C_{\phi \phi }} \otimes {W^{ - 1}}$$
(5)
$$W={\text{diag}}\left( {{{\left[ {{w^1}, \ldots ,{w^m}} \right]}^T}} \right)$$
(6)

where the weighting element \({w^s}\) to the observations of the satellite \(s\) is calculated using the elevation-dependent exponential function (Euler and Goad 1991):

$${w^s}={\left( {1+10~\exp \left( { - \frac{{{\beta ^s}}}{{10}}} \right)} \right)^{ - 2}}$$
(7)

in which \({\beta ^s}\) denotes the elevation angle in degrees from the rover to satellite \(s\). The \(f \times f\) matrices \({C_{{pp}}}\) and \({C_{\phi \phi }}\) are, respectively, variance matrices of the code and phase observable types in zenith.

Measurement setup

In this study, the processing is performed using 1 Hz code and phase QZSS data on L1, L2 and L5, collected by two baselines of different lengths in Perth, Australia. The first baseline, of around 9 m length, is formed by the stations CUAA and CUBB which are equipped with JAVAD TRE_G3TH DELTA receivers (Fig. 2). While the second baseline, of around 8 km length, is formed by the stations CUCC and UWA0 which are, respectively, equipped with a JAVAD TRE_G3TH DELTA receiver and a Septentrio PolaRx5 receiver (Fig. 3). Stations CUAA, CUBB and CUCC are connected to TRM59800.00 SCIS antennas, and station UWA0 is connected to a JAVRINGANT_DM SCIS antenna. The cut off elevation angle for all the data analysis in this contribution is set to be 10° and the broadcast ephemeris is used for the processing. Table 1 summarizes the information about the datasets in use.

Fig. 2
figure 2

Short baseline configuration. [Left] Baseline CUAA–CUBB located in Curtin University, Perth, Australia. [Middle] CUAA station. [Right] CUBB station

Full size image
Fig. 3
figure 3

Long baseline configuration. [Left] Baseline CUCC–UWA0 located in Perth, Australia. Map data @ 2018 Google (Google Earth 2018). [Middle] CUCC station. [Right] UWA0 station

Full size image
Table 1 Characteristics of the datasets used for this study
Full size table

Stochastic properties

RTK positioning based on (1)–(3) requires information about the structure of \({C_{{\text{pp}}}}\) in (4) and \({C_{\phi \phi }}\) in (5). Our aim here is to find representative values for the entries of these matrices. To do so, we first form the single-epoch, single-baseline DD code and phase residuals as follows

$$\left[ {\begin{array}{*{20}{c}} {{e_{p}}} \\ {{e_\phi }} \end{array}} \right]=\left[ {\begin{array}{*{20}{c}} p \\ \phi \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {{e_{f}} \otimes {I_{m - 1}}}&0 \\ {{e_{f}} \otimes {I_{m - 1}}}&{\Lambda \otimes {I_{m - 1}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} \rho \\ a \end{array}} \right]$$
(8)

where the \((m - 1)\)-vector of DD geometric ranges \(\rho\) is computed with the ground truth of the baseline coordinates and the broadcast satellite ephemeris. The \(f(m - 1)\)-vector \(a\) contains the reference integer ambiguities obtained from the strong baseline-known model. The residuals \({e_{p}}\) and \({e_\phi }\), in addition to the noise of the DD code and phase observations, also capture any mismodeled effects like, e.g. differential atmospheric delays and multipath. Therefore, one can write

$$\left[ {\begin{array}{*{20}{c}} {{e_{p}}} \\ {{e_\phi }} \end{array}} \right]=\left[ {\begin{array}{*{20}{c}} {{m_{p}}} \\ {{m_\phi }} \end{array}} \right]+~\left[ {\begin{array}{*{20}{c}} {{n_{p}}} \\ {{n_\phi }} \end{array}} \right]$$
(9)

in which \({m_{p}}\) and \({m_\phi }\) denote the mis-modeled effects on the DD code and phase observations respectively. Assuming that the geometry vector and ambiguity vector are deterministic, \({n_{p}}\) and \({n_\phi }\), respectively, denote the DD code and phase observations noise—a combined effect of the transmitted signals quality and the receiver architecture.

To estimate the entries of \({C_{{pp}}}\) and \({C_{\phi \phi }}\), we apply the least-squares variance component estimation (LS-VCE) method (Amiri-Simkooei et al. 2009; Teunissen and Amiri-Simkooei 2008) to the residuals \({e_{p}}\) and \({e_\phi }\) computed for the baselines CUAA–CUBB and CUCC–UWA0. As such, our estimations will capture the combined effects of the transmitted signals quality, the receiver architecture and any mismodeled effects like multipath. The multipath effect on the QZSS observations can be mitigated to a large extent through a procedure explained in the following.

Multipath

The QZSS satellite geometry repeats after one sidereal day. In addition, the surrounding environment of the static stations in use remains almost unchanged over time. Therefore, it is expected that the multipath pattern repeats after one sidereal day. As multipath corrections, one may then use the DD residuals in (8) corresponding with the same satellite geometry obtained after/before one sidereal day (Radovanovic 2000; Choi et al. 2004; Zaminpardaz et al. 2017; Zaminpardaz and Teunissen 2017). Subtracting these multipath corrections from the DD observations will largely mitigate the unwanted multipath effect, i.e. \({m_{p}} \to 0\) and \({m_\phi } \to 0\). Note, however, that subtracting the DD residuals of one single day from our dataset to diminish multipath would also increase the observations noise by a factor of \(\sqrt 2\), i.e. \({n_{p}} \to \sqrt 2 \times {n_{p}}\) and \({n_\phi } \to \sqrt 2 \times {n_\phi }\). To reduce such increase in the noise level of observations, one can make use of less noisy multipath corrections generated through averaging over the DD residuals of multiple days corresponding with the same satellite geometry. If the data of d days are used in this averaging, then the increase in the observations noise would be a factor of \(\sqrt {1+\frac{1}{d}}\). The multipath corrections then get less noisy if the data of more days become involved in its computation.

Multipath is the main source of time correlation for the short-baseline data. Therefore, to evaluate the efficiency of the multipath correction method explained above, we compute the time correlation of the QZSS signals collected by CUAA–CUBB on DOY 71 of 2018 before and after applying multipath corrections. The multipath corrections are taken from the corresponding observations on DOY 72 of 2018 (one single day). Figure 4 illustrates the time correlation of the QZSS code and phase observations on L1, L2 and L5 as functions of the time differences which are computed using the LS-VCE method. It can be seen that the significant time correlation of all the signals are dramatically reduced to negligible values upon applying the multipath corrections.

Fig. 4
figure 4

Time correlation of the QZSS code/phase observations on L1, L2 and L5 frequencies as functions of the time differences (left) before and (right) after multipath mitigation. The red dashed lines represent the 95% formal confidence interval, and a zoom-in window over the first 500 epochs is shown for each of the figures in the right panel. The data used for the plots are collected from CUAA and CUBB during UTC [02:46:39–03:46:40] on DOY 71 of 2018. The corresponding data on DOY 72 of 2018 is used for multipath mitigation

Full size image

Variance component estimation

Table 2 provides information about the zenith-referenced variance/covariance of QZSS undifferenced code and phase observations on all the three QZSS frequencies, which are the entries of \({C_{{pp}}}\) and \({C_{\phi \phi }}\). Two types of estimations are given for each parameter, specified as MP-uncorrected and MP-mitigated which are, respectively, computed based on the multipath-uncorrected and multipath-mitigated QZSS observations. The results are provided for both the CUAA–CUBB and CUCC–UWA0 baselines, by taking the average of the single-epoch LS-VCE estimations over the time intervals considered. Estimations for CUAA–CUBB are based on the data collected on DOY 71 of 2018 which are corrected for multipath using the observations on DOY 72 of 2018 (1 day), whereas estimations for CUCC–UWA0 are based on data collected on DOY 70 of 2018 which are corrected for multipath using the observations on DOY 71 of 2018 (1 day). The differences between the values corresponding with CUAA–CUBB and CUCC–UWA0 stem from several factors amongst which different receiver types and different multipath environments.

Table 2 Zenith-referenced standard deviation [m] and covariance [m2] of the QZSS undifferenced code and phase signals on L1, L2 and L5
Full size table

For both the baselines in use, our estimations confirm that the code–code and code–phase covariances corresponding with different frequency pairs are negligible, and therefore, we have skipped them in Table 2. We, therefore, consider \({C_{{pp}}}\) to be diagonal. However, as shown in Table 2, the phase–phase covariances are significant (w.r.t. phase precision) and cannot be neglected. With this in mind, in our following analyses, we make use of the fully-populated \({C_{\phi \phi }}\). We remark that for the CUCC–UWA0 baseline, as the residuals in (8) also contain non-negligible differential ZTDs, we also modeled the impact of ZTD when conducting LS-VCE estimation. Also, note that the multipath-mitigated estimations capture the noise level increase due to day-differencing. In other words, assuming that multipath would completely vanish through the mentioned day-differencing, then the multipath-mitigated results in Table 2 divided by \(\sqrt 2\) would represent the dispersion of the QZSS data in a multipath-free environment. With these in mind, we make two comparisons; once between multipath-uncorrected results and multipath-mitigated ones, and once between multipath-uncorrected results and multipath-free ones.

Comparing the results in the third (fourth) column of Table 2 with those in the fifth (sixth) column, we note that applying 1-day multipath corrections, although mitigates multipath to a large extent (see Fig. 4), increases the standard deviation of L1 and L2 code observations. Showing the best code precision, L5-code signal improves in precision upon applying multipath corrections. The degradation of L1- and L2-code precision can be circumvented if more precise multipath corrections are provided through involving the data of more than 1 day. For example, applying 6-day multipath corrections to the data of CUAA–CUBB baseline, our estimations for the L1-, L2- and L5-code standard deviation amount to 0.315, 0.242 and 0.103 m, respectively. The standard deviation of phase observations does not show significant dependency on frequency and multipath corrections.

Our second comparison is made between multipath-uncorrected standard deviations and multipath-free ones which, as previously mentioned, are obtained by dividing the results in the fifth and sixth columns of Table 2 by \(\sqrt 2\). By doing so, one can find out, when the observations are contaminated by multipath, that the code precision estimations degrade by a factor of around 1.2 for L1 and L2 signals, and by a factor of around 2.2 for the L5 signal. Multipath also deteriorates the phase precision estimations by a factor of around 1.2.

RTK performance

In this section, we analyze the capability of the current (February–March 2018) QZSS constellation, for instantaneous RTK positioning application. In doing so, we employ the QZSS observations on L1, L2 and L5 in their original form as well as their multipath-mitigated form. The data is collected at the rate of 1 Hz on DOY 69 of 2018 for CUAA–CUBB, and on DOY 72 of 2018 for CUCC–UWA0. The 1-day multipath corrections are provided as the corresponding data recorded on DOYs 70 and 73 of 2018, respectively. Figure 5 shows the 24-h skyplot view of the four QZSS satellites (J1, J2, J3 and J7) in Perth with the cut off elevation of 10°. In the following, we discuss the RTK results first for the CUAA–CUBB baseline, and then for the CUCC–UWA0 baseline.

Fig. 5
figure 5

The 24 h skyplot of the four QZSS satellites in Perth, Australia on DOY 69 of 2018. The cut off angle is set to be 10°

Full size image

CUAA–CUBB baseline

In this subsection, we present the standalone QZSS RTK results based on the data of the short-baseline CUAA–CUBB for which the differential atmospheric delays can be neglected, i.e. \(\Delta {\tau _{\text{Z}}} \approx 0\) in (2). Therefore, the only unknown parameters to be estimated are \(\Delta x\) and \(a\). The current QZSS constellation provides the users at Perth with four satellites being simultaneously visible (with the cut off angle of 10°) only for 50% of a 24-h period (see Fig. 1). However, even if four satellites are visible, they do not always form a geometry which is desirable from a positioning point of view. In this case, the position dilution of precision (PDOP) reaches extremely large values. To circumvent these extreme PDOP values, our evaluations will be carried out over the time intervals with four visible satellites considering different thresholds for PDOP computed as

$${\text{PDOP}}=\sqrt {{\text{tr}}\left\{ {{{\left( {{G^{T}}{P_{{{\text{D}}_{m}}}}WG} \right)}^{ - 1}}} \right\}}$$
(10)

with \(G= - {[{u^1},{u^2}, \ldots ,{u^m}]^{T}}\), \({\text{tr}}\{ .\} ~\) the trace operator and \({P_{{{\text{D}}_{m}}}}={D_{m}}{\left( {D_{{m}}^{{T}}{W^{ - 1}}{D_{m}}} \right)^{ - 1}}D_{{m}}^{{T}}{W^{ - 1}}\) the orthogonal projector. The single-epoch PDOP time series at Perth during DOY 69 of 2018 with the cut off angle of 10° is depicted in Fig. 6 employing the criterion of PDOP < 100.

Fig. 6
figure 6

Single-epoch PDOP time series (cf. 10) at Perth on DOY 69 of 2018. Time points satisfy the criterion of PDOP < 100

Full size image

Ambiguity resolution

Realization of RTK positioning relies on successful phase ambiguity resolution. In this subsection, we evaluate the triple-frequency QZSS ambiguity resolution performance using two measures, i.e. ambiguity dilution of precision (ADOP) and Integer Bootstrapped (IB) success rate—probability of correct integer estimation. As introduced by Teunissen (1997), ADOP is an easy-to-compute scalar diagnostic, which measures the model strength for successful ambiguity resolution. With the observational model in (1), ADOP is computed as

$${\text{ADOP}}={\sqrt {\left| {{Q_{\widehat {a}\widehat {a}}}} \right|} ^{~~\frac{1}{{f(m - 1)}}}}$$
(11)

where \({Q_{\widehat {a}\widehat {a}}}\) denotes the variance matrix of the float ambiguities and \(\left| \cdot \right|\) stands for the determinant. Figure 7 shows the single-epoch ADOP time series of the QZSS triple-frequency multipath-mitigated data collected by the baseline CUAA–CUBB on DOY 69 of 2018. The multipath corrections are provided as the corresponding data on DOY 70 of 2018. The time points satisfy the criterion of PDOP < 100. According to Odijk and Teunissen (2008), low ADOP values indicate high ambiguity success rates. As a rule of thumb, the integer least-squares (ILS) ambiguity success rate becomes larger than 99.9% when the ADOP value reduces to below 0.12 cycles. Therefore, since the ADOP value in Fig. 7 does not exceed 0.08, a high ambiguity resolution success rate is expected. Here, we remark that the ADOP time series has a different signature from that of the PDOP time series in Fig. 6. For example, just after sample 20000, while PDOP increases dramatically, ADOP shows a decreasing behavior. This demonstrates importantly that the PDOP and ADOP characteristics can be quite distinct and that one, therefore, should not confuse a poor PDOP with poor ambiguity resolution capabilities and correspondingly poor ambiguity-fixed positioning performance. Hence, a positioning-wise good geometry (i.e. small PDOP) is not necessarily needed for successful ambiguity resolution (Zaminpardarz et al. 2017).

Fig. 7
figure 7

Single-epoch ADOP time series (cf. 11) of QZSS triple-frequency multipath-mitigated data collected by the baseline CUAA–CUBB on DOY 69 of 2018. Time points satisfy the criterion of PDOP < 100

Full size image

The formal IB success rate, for the ambiguity \(f(m - 1)\)-vector in (1), is given as (Teunissen 1998)

$${\text{Formal}}~{\text{Ps}}=\mathop \prod \limits_{{i=1}}^{{f(m - 1)}} \left( {2\Phi \left( {\frac{1}{{2{\sigma _{{{\widehat {z}}_{i|I}}}}}}} \right) - 1} \right)$$
(12)

with

$$\Phi \left( x \right)= \int \limits_{{ - \infty }}^{x} \frac{1}{{\sqrt {2\pi } }}\exp \left( { - \frac{1}{2}{y^2}} \right){\text{d}}y$$
(13)

where \({\sigma _{{{\widehat {z}}_{i|I}}}}\) is the conditional standard deviation of the \(i{\text{th}}\) decorrelated ambiguity with \(~i=1, \ldots ,f\left( {m - 1} \right)\) and \(I=1, \ldots ,i - 1\). Note that IB success rate is the sharpest lower bound to the ILS success rate, which itself has the highest success rate of all admissible integer estimators (Teunissen 1999; Verhagen and Teunissen 2014).

Table 3, for the short baseline CUAA–CUBB based on the QZSS triple frequency data on DOY 69 of 2018, gives the empirical and formal IB success rates. The results are provided based on both the MP-uncorrected and the MP-mitigated data (based on 1-day multipath corrections). The formal values are computed by taking an average of the formal IB success rate (12) of all the selected processing epochs, whereas the empirical IB success rate is computed as.

Table 3 Single-epoch empirical \({\text{Ps}}\) (cf. 14) and formal \({\text{Ps}}\) (cf. 12) IB success rate for QZSS triple-frequency case
Full size table
$${\text{Empirical}}\;{\text{Ps}}=\frac{{{\text{number}}\;{\text{of}}\;{\text{the}}\;{\text{processing}}\;{\text{epochs}}\;{\text{with}}\;{\text{the}}\;{\text{correctly}}\;{\text{fixed}}\;{\text{ambiguities}}}}{{{\text{number}}\;{\text{of}}\;{\text{the}}\;{\text{processing}}\;{\text{epochs}}}}$$
(14)

To judge whether a DD ambiguity vector is correctly fixed, its corresponding IB solution is compared with the reference integer DD ambiguity vector computed based on the ILS solution of the baseline-known model. The results in Table 3 are classified based on a PDOP threshold. The reason for such classification is that for the positioning results, coming in the next subsection, we consider these thresholds for PDOP value. One should, however, bear in mind that the ambiguity resolution performance is not characterized through PDOP (Zaminpardaz et al. 2017). In Table 3, both the multipath-uncorrected and the multipath-mitigated data show a high level of success rate around 96 and 100%, respectively, revealing that instantaneous RTK is feasible using standalone QZSS if the data on all the three frequencies are employed. Also, the formal and empirical ambiguity success rates are well in agreement with each other, confirming the consistency between model and data.

Positioning performance

Having analyzed the ambiguity resolution performance of the standalone QZSS triple-frequency short-baseline model, here we investigate the corresponding ambiguity-float and -fixed positioning accuracies. Considering PDOP < 100, the deviations of the coordinate estimates from the ground truth are plotted in Fig. 8. Three types of solutions are shown in each panel; float solutions as gray dots, correctly-fixed solutions as green dots and wrongly-fixed solutions as red dots. It can be observed that the float solutions (gray), except for some intervals, match well with their corresponding 95% formal confidence intervals (blue curves). The inconsistencies between the formal and empirical float solutions can be attributed to the fact that the multipath corrections do not remove the whole multipath effect, but they capture the low-frequency multipath components to a large extent and partly the high-frequency multipath components. To elaborate more on this, we focus on the inconsistency showing up during the interval [9000, 13,000]. Our analysis shows that the presence of this inconsistency can be traced back to the high-frequency multipath residuals in the DD observations corresponding with the satellite pair J4–J3. Figure 9, illustrates the DD code residuals on L5 (cf. 8) corresponding with J4–J3 over the interval [9000, 15,000]. The top panel shows these residuals on DOYs 69 (blue) and 70 (purple) of 2018, while the bottom panel shows the difference of these residuals. It can be seen that, over the period [9000, 13,000], once the data are corrected for multipath, while the trend of the multipath is largely eliminated, the high-frequency multipath residuals are not mitigated, thus affecting the corresponding positioning results.

Fig. 8
figure 8

Deviations of the CUAA–CUBB baseline single-epoch solutions from the ground truth in (left) multipath-uncorrected and (right) multipath-mitigated cases. The gray, green and red dots represent the solutions with float, correctly-fixed and wrongly-fixed ambiguities, respectively. The 95% formal confidence intervals for the float solutions are marked by the blue curves. The QZSS triple-frequency data on DOY 69 of 2018 is used for the plots and the corresponding data on DOY 70 of 2018 is used for multipath mitigation

Full size image
Fig. 9
figure 9

Illustration of multipath residuals. [Top] DD code residuals on L5 (cf. 8) corresponding with satellite pair J4–J3 and receiver pair CUAA–CUBB over the interval [9000, 15,000] in Fig. 8. The blue graph is obtained on DOY 69 of 2018, while the purple graph is obtained on DOY 70 of 2018. [Bottom] Difference of the two graphs shown in the top panel

Full size image

Table 4 lists the single-epoch empirical and formal standard deviations of CUAA–CUBB baseline components before and after multipath mitigation for both the ambiguity-float and -fixed cases, considering different thresholds for PDOP. The ambiguity-fixed results are computed based on the correctly-fixed solutions. The formal outcomes are obtained from taking the average of all the single-epoch least-squares baseline variance matrices, whereas the empirical outcomes are obtained from the differences between the estimated baseline and the available ground truth of the CUAA–CUBB baseline. It can be observed that the meter-level coordinate errors in the ambiguity-float case are reduced to cm-level after fixing ambiguities in all three directions. In both the ambiguity-float and -fixed cases, the empirical and formal results match well with each other, validating the stochastic model used for the processing.

Table 4 The single-epoch empirical and formal standard deviations of the CUAA–CUBB baseline estimations based on QZSS triple-frequency observations collected on DOY 69 of 2018, for both the ambiguity-float and -fixed cases
Full size table

From Table 4, it can be seen that multipath mitigation has indeed improved the precision of float solutions, but marginally degraded that of fixed solutions. However, one should keep in mind that these fixed solutions would only be available in case the ambiguities can successfully be fixed, and as Table 3 shows, the probability of successful ambiguity resolution gets higher when the data is corrected for multipath. To explain the impact of multipath corrections on baseline solutions precision, we first give the single-epoch baseline variance matrix for both the ambiguity-float case \({Q_{\widehat {x}\widehat {x}}}\) as well as the ambiguity-fixed case \(Q_{\check{x}\check{x}}\)

$${Q_{\widehat {x}\widehat {x}}}=\frac{2}{{e_{{f}}^{{T}}C_{{{pp}}}^{{ - 1}}~{e_{f}}}}{\left[ {\mathop \sum \limits_{{s=1}}^{m} {w^s}\left( {{u^s} - \overline {u} } \right){{\left( {{u^s} - \overline {u} } \right)}^{T}}} \right]^{ - 1}}$$
(15)
$${Q_{\widehat {x}\widehat {x}}}=\frac{2}{{e_{{f}}^{{T}}(C_{{{pp}}}^{{ - 1}} + C_{{{\phi \phi}}}^{{ - 1}})~{e_{f}}}}{\left[ {\mathop \sum \limits_{{s=1}}^{m} {w^s}\left( {{u^s} - \overline {u} } \right){{\left( {{u^s} - \overline {u} } \right)}^{T}}} \right]^{ - 1}}$$
(16)

with \(\overline {u} =\left( {{{\mathop \sum \nolimits_{{s=1}}^{m} {w^s}{u^s}} \mathord{\left/ {\vphantom {{\mathop \sum \nolimits_{{s=1}}^{m} {w^s}{u^s}} {\mathop \sum \limits_{{s=1}}^{m} {w^s}}}} \right. \kern-0pt} {\mathop \sum \limits_{{s=1}}^{m} {w^s}}}} \right)\) being the weighted-average of the lines of sight. The difference between the baseline variance matrices corresponding with multipath-uncorrected and -mitigated data lies in the difference between the entries of the corresponding matrices \({C_{{pp}}}\) and \({C_{\phi \phi }}\) (see Table 2). Also, note that the multipath-mitigated results in Table 4 are obtained based on 1-day multipath corrections. As a result, the ratio of the term \({1 \mathord{\left/ {\vphantom {1 {e_{{f}}^{{T}}C_{{{pp}}}^{{ - 1}}~{e_{f}}}}} \right. \kern-0pt} {e_{{f}}^{{T}}C_{{{pp}}}^{{ - 1}}~{e_{f}}}}\) for multipath-uncorrected case to that of the multipath-mitigated case amounts to around 1.68, and \(\sqrt {1.68} \approx 1.3\) is then the factor by which the ambiguity-float positioning precision improves upon applying multipath corrections.

Now, we consider the ambiguity-fixed case where the ratio of the term \({1 \mathord{\left/ {\vphantom {1 {e_{{f}}^{{T}}\left( {C_{{{pp}}}^{{ - 1}}+C_{{\phi \phi }}^{{ - 1}}} \right)~{e_{f}}}}} \right. \kern-0pt} {e_{{f}}^{{T}}\left( {C_{{{pp}}}^{{ - 1}}+C_{{\phi \phi }}^{{ - 1}}} \right)~{e_{f}}}}\) for the multipath-uncorrected case to that of the multipath-mitigated case amounts to around 0.63. Therefore, it is expected that ambiguity-fixed positioning precision degrades upon applying multipath corrections by almost a factor of \(1/\sqrt {0.63} \approx 1.26\). This degradation, however, can be circumvented if one, instead of 1-day multipath corrections, makes use of less noisy multipath corrections. Based on 6-day multipath corrections, our evaluations confirm an instantaneous ambiguity resolution success rate more than 99.9%, and the consequent ambiguity-fixed results have the same precision as those based on multipath-uncorrected observations.

Among the three baseline components, the poorest precision is obtained for the East component while the best precision is achieved for the North component which can be explained by the behavior of the relative QZSS satellite geometry in Perth as follows. Here, as an example, we take the ambiguity-float case with the original data for which the single-epoch baseline variance matrix is given in (15). Figure 10 shows, over the observational period with PDOP < 100, the absolute values of the East, North and Height components of \(\sqrt {{w^s}} \left( {{u^s} - \overline {u} } \right)\). The larger such component is, the better this component can be estimated. Therefore, according to Fig. 10, one can indeed expect the precision of the baseline components can be ranked in an ascending order as East (red), Height (blue) and North (green) component.

Fig. 10
figure 10

Time series of the absolute value of \(\sqrt {{w^s}} \left( {{u^s} - \overline {u} } \right)\) (cf. 15, 16) components of the QZSS satellites w.r.t. the CUAA–CUBB baseline on DOY 69 of 2018, over the intervals satisfying PDOP < 100

Full size image

CUCC–UWA0 baseline

So far, we have investigated the QZSS performance in short-baseline RTK positioning. In this section, we extend our evaluations to a longer baseline, i.e. CUCC–UWA0 of around 8 km length (Fig. 3). For this longer baseline, using the QZSS data on three frequencies, we checked the significance of the differential ionospheric and tropospheric delays. Our evaluations showed that while the impact of atmospheric residuals on the noisy code observations can still be neglected, the differential tropospheric delays impact on the precise phase observations can reach significant levels compared to the phase precision. Therefore, when parametrizing \(\widetilde {\rho }\) in (1) for the observations of the CUCC–UWA0 baseline, in addition to \(\Delta x\), one needs to take into account \(\Delta {\tau _{\text{Z}}}\) as well.

Ambiguity resolution

Note, since we are working with the observations of four QZSS satellites, that parametrization of vector \(\widetilde {\rho }\) according to (2) would not bear additional redundancies for the single-epoch observational model in (1), implying that ADOP and ambiguity resolution performance would remain invariant w.r.t. such parametrization. Therefore, here, we assess the ambiguity resolution performance without further parametrizing \(\widetilde {\rho }\). Figure 11 shows the single-epoch ADOP time series of the QZSS triple-frequency multipath-mitigated data collected by the baseline CUCC–UWA0 on DOY 72 of 2018. The time points satisfy the criterion of four satellites being visible simultaneously. Since the ADOP value in Fig. 11 does not exceed 0.08, a high ambiguity resolution success rate is expected. Our results confirm an instantaneous ambiguity resolution success rate of 99.65%.

Fig. 11
figure 11

Single-epoch ADOP time series (cf. 11) of QZSS triple-frequency multipath-mitigated data collected by the baseline CUCC–UWA0 on DOY 72 of 2018

Full size image

Positioning performance

As was mentioned earlier, \(\widetilde {\rho }\) for the observations of the CUCC–UWA0 baseline contains significant residual tropospheric delays. However, with ZTD included as an unknown parameter, instantaneous RTK will require the observations of at least five satellites simultaneously. This in tandem with the fact that the data of only four QZSS satellites are available make single-epoch positioning infeasible for CUCC–UWA0. To tackle this issue, instead of the measurements of a single epoch, we incorporate the data of two epochs assuming that the baseline coordinates remain invariant. To realize a sufficient change of receiver–satellite geometry, these two epochs are chosen as follows.

  • Check the number of visible satellites on a given day and store the indicator of the epochs with four satellites in \(\Omega\);

  • If the size of \(\Omega\), denoted by \(k\), is odd, make it even by removing one element;

  • Define \({\Omega _1}\) containing the first \(\frac{k}{2}\) elements of \(\Omega\), and \({\Omega _2}\) containing the last \(\frac{k}{2}\) elements of \(\Omega\);

  • For two-epoch positioning, take one set of DD observations at an epoch from \({\Omega _1}\) and one set at an epoch from \({\Omega _2}\).

To assess the performance of two-epoch positioning explained above, use is made of the QZSS triple-frequency multipath-mitigated data collected by the baseline CUCC–UWA0 on DOY 72 of 2018. The multipath corrections are given by the corresponding data on DOY 73 of 2018. Through the above procedure, the maximum and minimum time differences between the epochs in the epoch pairs considered are, respectively, 09:28:31 and 14:27:19. As was discussed in the previous subsection, the corresponding ambiguities can successfully be resolved to their integer values already based on single-epoch multipath-corrected observations. Hence, based on a stronger model formed by the observations of two epochs with the assumption of baseline components being time-constant, ambiguities must successfully be fixed. In the following, we, therefore, present the results corresponding with the multipath-corrected data.

The deviations of the ambiguity-float (gray) and -fixed (green) coordinate estimates from the ground truth are plotted in Fig. 12. It can be observed that the float solutions (gray) match well with their corresponding 95% formal confidence intervals (blue curves). In this figure, we also note some sudden changes in the positioning accuracy (see blue curves). This happens as a result of sudden changes in the time difference between the two epochs in use, thereby having significant changes in satellite geometry. For example, the standard deviation drop from data number 13,570 to 13,571 is due to that fact that the time difference between the two epochs in use increases (from 9 h 34 min 25 s to 13 h 1 min 46).

Fig. 12
figure 12

Deviations of the CUCC–UWA0 baseline two-epoch solutions from the ground truth in ambiguity-float (gray) and ambiguity-fixed (green) case. The 95% formal confidence intervals for the float solutions are marked by the blue curves. The QZSS triple-frequency data on DOY 72 of 2018 is used for the plots and the corresponding data on DOY 73 of 2018 is used for multipath mitigation

Full size image

Table 5 provides the precision of baseline float and fixed solutions. To compare the results of this table with those in Table 4 (for e.g. PDOP < 100), we divide the standard deviations in Table 4 by \(\sqrt 2\) to achieve the approximate two-epoch positioning precision for the short-baseline case. One can then realize that the inclusion of ZTD as an unknown parameter will deteriorate the precision of the baseline components. The larger deterioration in the vertical direction can be explained by the strong correlation between the ZTDs and the receiver height component (Rothacher and Beutler 1998).

Table 5 The two-epoch empirical and formal standard deviations of the CUCC–UWA0 baseline estimations based on QZSS triple-frequency observations collected on DOY 72 of 2018, for the ambiguity-fixed case
Full size table

Summary and conclusions

In this contribution, for the four-satellite QZSS constellation as a stand-alone system, we provided an initial assessment of L1 + L2 + L5 RTK performance in the form of integer ambiguity resolution success rate and positioning accuracy. We analyzed the noise characteristics of the QZSS code and phase observations collected by two baselines in Perth, Australia, i.e. CUAA–CUBB with identical receiver types and a length of 9 m and CUCC–UWA0 with mixed receiver types and a length of 8 km. Our results were presented in the form of zenith-referenced variance–covariance matrices and time correlations.

It was shown that through applying multipath corrections, the significant time correlation of all the QZSS code and phase observations become negligible. Zenith-referenced standard deviations of the code observations were the largest for the L1 signal and the lowest for the L5 signal. The phase observations on different frequencies, however, had almost the same precision. Applying multipath corrections reduced the L1, L2 and L5 code standard deviation by, respectively, around 15, 20 and 50%. Our estimations of covariances between different observation types implied that the code observations on different frequencies can be considered uncorrelated with each other and also with the phase observations, both before and after multipath corrections. However, significant covariances were observed between the phase observations on different frequencies. The standard deviation of the phase observations on each frequency had very low sensitivity w.r.t multipath correction.

Having investigated the characteristics of the QZSS signals, we then presented the standalone QZSS triple-frequency RTK performance. In doing so, we first considered the short-baseline CUAA–CUBB for which the differential atmospheric delays were negligible. In Perth, Australia, the current four-satellite constellation of QZSS provides PDOP values lower than 100 for more than 10 h per day during which, ADOP values are lower than 0.08 for triple-frequency case indicating a high ambiguity success rate. It was shown that the instantaneous integer bootstrapped (IB) success rate is around 96% and 100% for the original and multipath-mitigated data, respectively, implying that the instantaneous RTK is feasible based on standalone QZSS observations.

In both the ambiguity-float and -fixed cases, the empirical and formal results matched well with each other, confirming the consistency of the model and data. Meter-level coordinate errors in the ambiguity-float case were reduced to cm-level after fixing ambiguities in all three directions. The multipath mitigation improved the float solutions precision, and only marginally degrades the fixed solution accuracy which was due to using noisy multipath corrections. It was shown that such degradation can be prevented using less noisy multipath corrections computed based on the data of multiple days. The baseline components in terms of precision can be ranked in ascending order as East, Height and North. This was explained as a result of specific QZSS satellites geometry at Perth. Our illustrations also demonstrated that PDOP and ADOP characteristics can be completely different and that one, therefore, should not confuse a poor PDOP with poor ambiguity resolution capabilities. To predict ambiguity-fixed positioning performance, one should therefore, first consult the ADOP before considering the PDOP.

Finally, the standalone QZSS positioning performance was assessed for the longer baseline CUCC–UWA0 for which the impact of differential tropospheric delays on the precise phase observations can reach significant levels compared to phase precision. As such, the observational model for this baseline included, in addition to baseline components and DD ambiguities, the zenith tropospheric delay (ZTD) as well. In order for all the unknown parameters to be estimable, we made use of the data of two epochs assuming that the baseline coordinates are time-invariant. Based on multipath-mitigated data, it was shown that the ambiguity resolution success rate is equal to 100%. Comparing the ambiguity-resolved positioning performance of CUCC–UWA0 with that of short-baseline CUAA–CUBB, the precision of the baseline components gets poorer, particularly in the vertical direction, when including ZTD as an unknown parameter.