Regional Ionospheric Maps with Quad-Constellation Raw Observations as Applied to Single-Frequency PPP

Abstract

Ionospheric delay is one of the most problematic errors in single-frequency (SF) global navigation satellite system (GNSS) data processing. Global/regional ionospheric maps (GIM/RIM) are thus vitally important for positioning users. Given the coexistence of multi-GNSS, the integration of quad-constellation observations is essential for improving the distribution of ionospheric penetration points (IPPs) and increasing redundant observations compared with the existing GIM products from the IGS analysis center. In this paper, quad-constellation (GPS/GLONASS/Galileo/BDS) observations are applied to set up the RIM over Australia with uncombined precise point positioning (UC-PPP) and a low-order spherical harmonic function. The generated RIMs are then introduced to ionosphere-corrected (IC) and ionosphere-weighted (IW) single-frequency PPP (IC-SFPPP and IW-SFPPP) to verify their performance in terms of positioning accuracy and convergence time. Taking the CODE GIM as a reference, the results show that the mean root mean square (RMS) of VTEC differences is 0.867 TECUs, and the quad-constellation RIM (referred as ‘RIM4′) can improve the RMS of RIMs compared to single-constellation mode at the edge of regional experiment area. The application of the RIM4 in the BDS IC-SFPPP results in a 18.38% improvement (from 100.47 cm to 82.00 cm) of 3D positioning RMS compared to the CODE-GIMs, whereas 35.36% enhancement (from 115.92 cm to 74.62 cm) of 3D positioning RMS is achievable during an active ionospheric period. Moreover, if the criterion of the convergence time is defined as when positioning errors in the horizontal and vertical directions are less than 0.3 m and 0.6 m for 20 consecutive epochs, the IW-SFPPP can significantly speed up the convergence time compared to the uncombined SFPPP; that is, the convergence time is reduced by 52.7% (from 37 min to 17.5 min), 37.2% (from 72.5 min to 45.5 min), and 37.1% (from 62.0 min to 39.0 min) in the north, east and up direction, respectively, at the 68% confidence level.

1. Introduction

The completion of the third generation of the BeiDou navigation satellite system (BDS-3) satellite on 23 June 2020 [1,2], marks the coexistence of a global satellite navigation system (GNSS) with the global positioning system (GPS) GLObalnaya NAvigatsionnaya Sputnikovaya Sistema (GLONASS) and the Galileo satellite navigation system (Galileo). The advent of four GNSS represents a vital improvement of accuracy and robustness, not only in the application of a satellite-based positioning, navigation, and timing (PNT) service, but also in the application of space atmosphere monitoring [3,4,5]. The ionosphere is a crucial component of the atmosphere and contains extensive free electrons that disrupt the propagation of satellite signals [6]. The time delay imposed by the ionosphere can reach several or even hundreds of meters, which can severely degrade the positioning performance of GNSS users [7,8,9], especially for single-frequency users. On the other hand, GNSS, with its advantages of high temporal resolution, all-weather ready, and wide-coverage, has become a powerful tool for extracting ionospheric observations and studying the spatiotemporal variability of the ionosphere [10,11]. Hence, research on ionospheric monitoring and correction has been a hot research topic. In addition to the global ionospheric map (GIM) products from the ionospheric analysis center (IAC), the Madrigal database developed by the Massachusetts Institute of Technology (MIT) Haystack Observatory provides a MIT Automated Processing of GPS (MAPGPS) TEC product with a temporal resolution of 5 min and a spatial resolution of 1° × 1° [12]. The IONOLAB software [13], the Seemala software (http://seemala.blogspot.com (accessed on 30 November 2022)), the M_DCB software [14], and the System for Ionosphere Monitoring and Research from GNSS (SIMuRG) developed by Yasyukevich et al. [15], can all help us to monitor the ionosphere with dense GNSS networks.

As the applicable field of GNSS navigation is expanding, further requirements for vertical total electron content (VTEC) modelling are put forward, one of which is the need for high precision ionosphere products [16]. The existing GIM products from the IAC are mainly generated from GPS and GLONASS data [4], with an accuracy of 2–8 total electron content units (TECUs) [17], which is equivalent to a pseudorange error of 0.3–1.3 m at the BDS B1 frequency, which struggles to meet the demand of regional users for precise positioning. As a necessary complement to the GIMs, regional ionospheric maps (RIMs) tend to be more useful in disclosing the fine-scale changes in the ionosphere over a local area. Some scholars have attempted to develop regional ionospheric models by means of polynomial, trigonometric series, spherical harmonic, and multi-surface function, and some progress has also been made [18,19,20,21]. Li et al. evaluated the capabilities of these four models in different regions and the experimental results suggest that the polynomial and spherical harmonic functions have comparable performance in terms of their performance and computation efficiency and can be a priority choice for regional ionospheric modelling [22].

With the progressive development of GNSS, there is a growing demand for high-precision positioning from the mass users, whereas the cost of the hardware is an essential factor to be considered [23,24,25]. The promotion of single-frequency (SF) receivers and low-cost antennas is an effective alternative. Therefore, with the trend of BDS globalization, how to effectively use the BDS SF observations while ensuring positioning accuracy has become a key issue. Currently, the most widely used model is the ionosphere-corrected single-frequency precise point positioning (IC-SFPPP), the accuracy of which is tightly related to the external ionospheric products. Zhang et al. [26] proposed another approach named uncombined SFPPP (UC-SFPPP), which enables joint estimation of slant ionospheric delays (SIDs) and other parameters. Accounting for the long convergence time of the UC-SFPPP owing to the weak structure of the algorithm, an ionosphere-weighted SFPPP (IW-SFPPP) model has been proposed [27,28]. In this model, the external ionospheric information is treated as a pseudo-observation to strengthen the solution, where the key is to determine the proper constraints [29].

In this research, the raw multi-constellation GNSS code and phase observations are employed by un-combined PPP to estimate the ionospheric observations which are later used for generating the RIMs over Australia, considering that more constellations will improve the distribution of the ionospheric penetration points (IPPs) and increase redundant observations in the VTEC modelling. Since the commonly-used carrier-to-code leveling (CCL) method is easily affected by code noises, multipath effects and arc length [30,31], and the un-combined dual-frequency PPP (UC-DFPPP) method, which is proven to be flexible and effective [32,33,34], are applied to extract precise ionospheric observables. For convenience of analysis, the generated RIMs are referred to as RIM4, RIM2, RIMG, and RIMC, representing the RIMs modelled from quad-constellation (GREC), dual-constellation (GC), GPS-only, and BDS-only data, respectively.

The paper is organized as follows. After the brief introduction in Section 1, the regional VTEC modeling method and two algorithmic models of the SFPPP are introduced in Section 2. In Section 3, the four RIMs are generated over Australia, and their performance is evaluated by comparing with the GIM products from CODE and by applying to the single-frequency PPP users with the IC-SFPPP and IW-SFPPP models, respectively. Finally, the conclusions are drawn in Section 4.

2. Functional Model

In this section, the functional model of two basic steps in regional VTEC modelling, namely retrieving ionospheric observables and VTEC modelling, is presented. Then, two algorithmic models of the BDS SFPPP are introduced in detail.

2.1. Multi-GNSS Raw Observation Equations

The satellite signal is delayed by the ionosphere as it passes through the atmosphere, and we can extract the SIDs by parameter estimation in the UC-PPP [32]. The rank deficient multi-GNSS raw carrier phase and pseudorange observations can be given as follows [35]:

$$\{\begin{array}{l}{P}_r^s(i)={\rho }_r^s+\left(c\cdot d{t}_r-c\cdot d{t}^s\right)+T_r+{\gamma }_i{I}_{r,1}^s+H{D}_{r,i}-H{D}_i^s+{\epsilon }_{P_i}\\ {\mathsf{\Phi }}_r^s(i)={\rho }_r^s+\left(c\cdot d{t}_r-c\cdot d{t}^s\right)+T_r-{\gamma }_i{I}_{r,1}^s+{\lambda }_i^s{N}_{r,i}^s+{\epsilon }_{{\mathsf{\Phi }}_i}\end{array}$$

(1)

where the indices $i$, $s$, and $r$ refer to the frequency band, satellite, and receiver, respectively; $d{t}_r$ and $d{t}^s$ are the clock biases of the receiver and satellite multiplied by the speed of light c, whereas $T_r$ denotes the zenith tropospheric delay; $I_{r,1}^s$ is the SIDs on L1 signal and the frequency-dependent factor ${\gamma }_i$ can be expressed as ${\gamma }_i={\left(f_1/f_i\right)}^2$, with $f_i$ the frequency of $i$; $H{D}_{r,i}$ and $H{D}_i^s$ are the frequency-related receiver and satellite code hardware delays; $N_{r,i}^s$ denotes the phase ambiguity including the satellite and receiver phase hardware delays; and ${\epsilon }_{P_i}$ and ${\epsilon }_{{\mathsf{\Phi }}_i}$ are the effects of code and phase unmodeled errors, including multipath effects and observation noise.

In the following, note that the inter-frequency biases (IFB) of GLONASS are not ignored and are modelled as a linear function of the frequency number in this paper [36,37]. In addition, the inter-system biases (ISB) are treated as time constants [38,39]. Considering the ionosphere-free (IF) combinations of satellite code hardware delays in the IGS precision satellite clock products, the observation equations can be given after parameter reformulation [28]:

$$\{\begin{array}{l}\begin{array}{l}{P}_r^{s,G}(i)={\rho }_r^{s,G}+c\cdot d{\tilde{t}}_r{}^G+T_r+{\gamma }_i{\tilde{I}}_{r,1}^{s,G}+{\epsilon }_{P_i}\\ P_r^{s,R}(i)={\rho }_r^{s,R}+c\cdot d{\tilde{t}}_r{}^G+cIS{B}_r^R+cIF{B}_r^{s,R}+T_r+{\gamma }_i{\tilde{I}}_{r,1}^{s,R}+{\epsilon }_{P_i}\\ P_r^{s,E}(i)={\rho }_r^{s,G}+c\cdot d{\tilde{t}}_r{}^G+cIS{B}_r^E+T_r+{\gamma }_i{\tilde{I}}_{r,1}^{s,G}+{\epsilon }_{P_i}\end{array}\hfill \\ P_r^{s,C}(i)={\rho }_r^{s,G}+c\cdot d{\tilde{t}}_r{}^G+cIS{B}_r^C+T_r+{\gamma }_i{\tilde{I}}_{r,1}^{s,G}+{\epsilon }_{P_i}\hfill \end{array}$$

(2)

$$\{\begin{array}{l}\begin{array}{l}{\mathsf{\Phi }}_r^{s,G}(i)={\rho }_r^{s,\mathrm{G}}+c\cdot d{\tilde{t}}_r{}^G+T_r-{\gamma }_i{\tilde{I}}_{r,1}^{s,G}+{\lambda }_i^{s,G}{\tilde{N}}_{r,i}^{s,G}+{\epsilon }_{{\mathsf{\Phi }}_i}\\ {\mathsf{\Phi }}_r^{s,R}(i)={\rho }_r^{s,R}+c\cdot d{\tilde{t}}_r{}^G+cIS{B}_r^R+T_r-{\gamma }_i{\tilde{I}}_{r,1}^{s,R}+{\lambda }_i^{s,R}{\tilde{N}}_{r,i}^{s,R}+{\epsilon }_{{\mathsf{\Phi }}_i}\\ {\mathsf{\Phi }}_r^{s,E}(i)={\rho }_r^{s,E}+c\cdot d{\tilde{t}}_r{}^G+cIS{B}_r^E+T_r-{\gamma }_i{\tilde{I}}_{r,1}^{s,E}+{\lambda }_i^{s,E}{\tilde{N}}_{r,i}^{s,E}+{\epsilon }_{{\mathsf{\Phi }}_i}\end{array}\hfill \\ {\mathsf{\Phi }}_r^{s,C}(i)={\rho }_r^{s,C}+c\cdot d{\tilde{t}}_r{}^G+cIS{B}_r^C+T_r-{\gamma }_i{\tilde{I}}_{r,1}^{s,C}+{\lambda }_i^{s,C}{\tilde{N}}_{r,i}^{s,C}+{\epsilon }_{{\mathsf{\Phi }}_i}\hfill \end{array}$$

(3)

where symbols with identifiers indicate the reparametrized estimable parameters (seen in

Table 1. The expressions for estimable parameters in Equations (2) and (3).

Parameter	Form
$d{\tilde{t}}_r{}^G$	$d{\tilde{t}}_r{}^G=d{t}_r^G+d_{r,I{F}_{12}}^G$
$IS{B}_r^Q$	$IS{B}_r^Q=(d_{r,I{F}_{12}}^{\mathrm{Q}}-d_{r,I{F}_{12}}^G)+(d{\mathrm{T}}^Q-d{\mathrm{T}}^G) (Q\ne G)$
${\tilde{I}}_{r,1}^{s,Q}$	${\tilde{I}}_{r,1}^{s,Q}=I_{r,1}^{s,Q}-\frac{c}{{\gamma }_2-1}\left(DC{B}_r^Q-DC{B}^{s,Q}\right)$
${\lambda }_i^{s,Q}{\tilde{N}}_{r,i}^{s,Q}$	$\begin{array}{c}{\lambda }_i^{s,Q}{\tilde{N}}_{r,i}^{s,Q}={\lambda }_i^{s,Q}{N}_{r,i}^{s,Q}+c\left(d_{I{F}_{12}}^{s,Q}-d_{r, I{F}_{12}}^Q\right)\\ -\frac{c{\gamma }_i}{{\gamma }_{{}_2}-1}\left({\mathrm{DCB}}_r^Q-{\mathrm{DCB}}^{s,Q}\right)\end{array}$

); ${\rho }_r^{s,*}$ is the geometric distance between the satellite $s,*$ and receiver $r$, where $G$, $R$, $E$ and $C$ represent GPS, GLONASS, Galileo, and the BDS satellite, respectively; $T_r$ denotes the zenith tropospheric delay; $IS{B}_r^Q$ contains the corresponding differences in receiver hardware delays ($d_{r,I{F}_{12}}^{\mathrm{Q}}-d_{r,I{F}_{12}}^G$) and time datums ($d{T}^Q-d{T}^G$); and ${\tilde{I}}_{r,1}^s$ contains the satellite and receiver differential code bias (DCB) [40] ($DC{B}_r=H{D}_{r,1}-H{D}_{r,2}$, $DC{B}^s=H{D}_1^s-H{D}_2^s$).

2.2. Regional Ionospheric VTEC Modeling

With the SIDs retrieved from the general multi-GNSS PPP model, the slant total electron contents (STEC) can be calculated as [41]:

$$STEC=I_{r,1}^s\cdot \frac{f_1^2}{40.28}$$

(4)

where $f_1$ stands for the frequency of the carrier. To model the ionosphere over a region, the STEC is mapped into VTECs using the single-layer model (SLM) assumption [42], which can be expressed as:

$$\{\begin{array}{l}VTEC=Mf\cdot \mathrm{STEC}\\ Mf=\cos \left(\arcsin \left(\frac{R}{R+H}\sin (Z)\right)\right)\end{array}$$

(5)

where $Mf$ is the mapping function and $Z$ denotes the zenith angle of a satellite at the receiver, and $H$ refers to the height of the assumed single-layer ionosphere ($H=450 \mathrm{km}$), whereas $R$ is the earth’s radius ($R=6378137 \mathrm{m}$).

In Refs. [22,43], the performance of the polynomial and spherical harmonic models for VTEC modeling is comprehensively assessed within several typical regional areas. The results show that, in small regions, the performances of all models are insensitive to the model orders but sensitive to the ionosphere activity. On the whole, the polynomial and spherical harmonic function are comparable in terms of their performance and computation efficiency. For the regional VTECs retrieved from the original observations above, we use a low-order spherical harmonic function to describe its spatiotemporal variation, which can be presented as follows [33]:

$$VTEC={\displaystyle \sum _{n=0}^{n_{\max }}{\displaystyle \sum _{m=0}^n{\overline{P}}_{nm}}}(\sin {\phi }^{\prime })\left(a_{nm}\cos (m{\lambda }^{\prime })+b_{nm}\sin (m{\lambda }^{\prime })\right)$$

(6)

where ${\overline{P}}_{nm}$ denotes the normalized associated Legendre function; ${\phi }^{\prime }$ and ${\lambda }^{\prime }$ are the geomagnetic latitude and longitude of the IPPs in the solar fixed reference frame, respectively; and $a_{nm}$ and $b_{nm}$ represent the normalized coefficients. Moreover, a set of RIM coefficients is estimated hourly and the segmented linear method is used to connect adjacent ionospheric TEC models.

Another issue that needs to be addressed is that SIDs extracted from the UC-PPP model still have an inseparable receiver DCB (RDCB) and satellite DCB (SDCB), of which the columns in the design matrix have linear dependencies on each other. To eliminate this rank-deficiency, the condition that the sum of the DCBs for all satellites is zero is imposed as a constraint [44,45,46] and can be expressed as:

$${\displaystyle \sum _{i=1}^{Nu{m}^Q}SD}C{B}^Q=0$$

(7)

where $Nu{m}^Q$ represents the number of visible satellites in the selected region for constellation Q and the remaining RDCB are estimated as constants. With the above theories, the DCB parameters and ionospheric coefficients can be derived by the least-squares method.

2.3. Functional Model of SFPPP

According to Equations (1)–(3), the observations with the first frequency of BDS can be expressed as [26]:

$$\{\begin{array}{l}{P}_{r,1}^s={\rho }_r^s+c\cdot d{\overline{t}}_r+{\overline{I}}_{r,1}^s+T_r+{\epsilon }_{P_1}\hfill \\ {\mathsf{\Phi }}_{r,1}^s={\rho }_r^s+c\cdot d{\overline{t}}_r-{\overline{I}}_{r,1}^s+T_r+{\lambda }_1^s\cdot {\overline{N}}_{r,1}^s+{\epsilon }_{{\mathsf{\Phi }}_1}\hfill \end{array}$$

(8)

where symbols with identifiers indicate the reparametrized estimable parameters in

Table 2. The expressions for estimable parameters in Equation (8).

Parameter	Form
$d{\overline{t}}_r$	$c\cdot d{\overline{t}}_r=c\cdot d{t}_r+H{D}_{r,1}$
${\overline{I}}_{r,1}^s$	${\overline{I}}_{r,1}^s=I_{r,1}^s+\frac{c}{{\gamma }_2-1}\cdot DC{B}^s$
${\lambda }_i^s{\overline{N}}_{r,i}^s$	${\lambda }_1^s{\overline{N}}_{r,1}^s={\lambda }_1^s{N}_{r,1}^s-H{D}_{r,1}+\frac{c}{{\gamma }_{{}_2}-1}{\mathrm{DCB}}^s$

. The biased SIDs ${\overline{I}}_{r,1}^s$ contains the SDCBs after the introduction of the external precise satellite orbit and clock products, and the estimable receiver clock $d{\overline{t}}_r^{}$ absorbs the receiver hardware delay at the first frequency $H{D}_{r,1}^{}$.

The SFPPP algorithm can be divided into IC-SFPPP, UC-SFPPP, and IW-SFPPP, depending on the different treatment of the ionospheric delays ($I_{r,1}^s$). First, after correcting the SIDs directly with the ionospheric VTEC products and fixing the SDCBs to the day-mean values from the Chinese Academy of Sciences (CAS), the functional model of IC-SFPPP can be expressed as [47]:

$$\{\begin{array}{l}{P}_{r,1}^s={\rho }_r^s+c\cdot d{\overline{t}}_r+T_r+{\epsilon }_{P_1}\hfill \\ {\mathsf{\Phi }}_{r,1}^s={\rho }_r^s+c\cdot d{\overline{t}}_r+T_r+{\lambda }_1^s\cdot {\overline{N}}_{r,1}^s+{\epsilon }_{{\mathsf{\Phi }}_1}\hfill \end{array}$$

(9)

Unlike IC-SFPPP, the external ionospheric information is used as the pseudo-measurement in IW-SFPPP and can be given as [28]:

$$\{\begin{array}{l}{P}_{r,1}^s={\rho }_r^s+c\cdot d{\overline{t}}_r+T_r+I_{r,1}^s+\frac{f_2^2}{f_1^2-f_2^2}DC{B}^s+{\epsilon }_{P_1}\hfill \\ \begin{array}{l}{\mathsf{\Phi }}_{r,1}^s={\rho }_r^s+c\cdot d{\overline{t}}_r+T_r-I_{r,1}^s+{\lambda }_1^s\cdot {\overline{N}}_{r,1}^s+{\epsilon }_{{\mathsf{\Phi }}_1}\\ {\tilde{l}}_{r,1}^s=I_{r,1}^s+{\epsilon }_{ion}\end{array}\hfill \end{array}$$

(10)

where ${\tilde{l}}_{r,1}^s$ denotes the VTEC information obtained by external ionospheric products, whereas ${\epsilon }_{ion}$ is the noise of the pseudo-observations.

Accounting for the accuracy limitations of VTEC products, virtual ionospheric observations assume a more significant role at the initiation of filtering to accelerate convergence and gradually reduce the weights to obtain better localization. This time-varying weighting scheme can be expressed as [48]:

$${\sigma }^2={\sigma }_0^2+\alpha \cdot \mathsf{\Delta }t$$

(11)

where ${\sigma }_0^2$ denotes the initial variance; $\alpha$ is variance change rate over time; and $\mathsf{\Delta }t$ refers the time interval from the start moment of the Kalman filtering.

3. Datasets and Experiments

In this section, the experiment datasets and processing strategies are introduced firstly. Then, the reliability of the established RIMs are validated by comparing with the CODE GIM products, and their application performance in the SFPPP is analyzed using the IC-SFPPP and IW-SFPPP models.

3.1. Data Collection and Processing Strategy

In this study, we collected quad-constellation (GREC) dual-frequency (DF) observations of 30 reference stations from Geoscience Australia (GA), during DOY 183–189 in 2021 for regional ionospheric VTEC modelling, and the distribution of reference stations is shown in Figure 1. Moreover, the BDS B1 observations of 10 additional monitoring stations (blue triangle) have been selected to verify the SFPPP positioning performance with the introduction of the RIMs.

Two conditions need to be considered when selecting these reference stations; one is that they should be well distributed across the region to ensure good coverage of the IPPs, and the other is that observations should be available for all four constellations. The latter requirement is to ensure fair comparisons for both RIM evaluation as well as for single-frequency PPP validation. More specifically, in the experiment analysis, the performance of four different RIMs over Australia, which are referred to as RIM4, RIM2, RIMG, and RIMC, are evaluated comparatively by comparing with the GIM products from CODE and by applying the single-frequency PPP users with the IC-SFPPP and IW-SFPPP models. To ensure fair comparison, the processing strategies of ionosphere modelling as well as the SF-PPP users are all the same, the only difference is the input of GNSS observations for ionosphere modelling, which means that the four RIMs are modelled from quad-constellation (GREC), dual-constellation (GC), GPS-only, and BDS-only data, such that the inclusion of this condition avoids possible performance differences caused by the different distribution of reference stations.

The distribution of IPPs in the ionospheric VTEC modelling for the first hour of DOY 183 in 2021 is shown in Figure 2. It can be seen that IPPs are well-distributed in Australia, whereas the IPP coverage is poor at the boundaries of the survey area, especially over the ocean. Although the inclusion of multi-GNSS data has improved this trend, it still has a negative impact on ionospheric modelling. In addition, with the completion of the BDS-3 network, there are more visible BDS satellites in the Australian one, which provides an opportunity to monitor ionospheric changes with BDS-only data.

The data processing can be divided into three steps: (1) precise ionospheric information inversion by the multi-GNSS UC-PPP, (2) regional VTEC modelling, and (3) the performance verification of RIMs in BDS SFPPP. In the UC-PPP procedure, the final GFZ products are used with a time interval of 5 min and 30 s for precise orbit and clock products, respectively. The time lag of the final GFZ products is 12–18 days. To eliminate the convergence time, the PPP results are smoothed with forward and backward Kalman filtering in ionospheric information retrieval. The experiments integrate observations from GPS C1W-C2W, GLONASS C1P-C2P, BDS C2I-C6I, and Galileo C1X-C5X. It should be noticed that GLONASS requires additional consideration of the IFB parameters, which are modelled as a linear function of frequency in this paper [49]. Meanwhile, the ISB parameters are treated as time-constant [28]. In the regional VTEC modelling, a spherical harmonic function of four orders and four degrees is chosen with an interval of one hour [33]. The data processing is carried out in the post-processing mode. For selected monitoring stations, the average coordinates of the static DFPPP results over 1 day have been used as the reference position.

Table 3. Data processing strategies in the experiments.

Items	Processing Strategies
1. PPP processing
Systems	Quad-system (G + R + E + C) joint solution
Weighting strategies	A priori accuracy of 0.003 and 0.3 m for zenith raw phase and code observations is assumed; GLONASS code IFB: linear function of frequency number;
Clock and orbits	GFZ products
Troposphere modeling	Saastamoinen model and random walk process
Elevation cutoff	10°
Sampling interval	30 s
Slant ionospheric delays	Estimated as random-walk noises (0.04 m2/s)
Kalman filter	Smoothed with forward and backward filtering
Phase ambiguities	Estimated as float constants
2. Regional ionospheric modeling
Slant ionospheric delays	Retrieved from quad-system UC-DFPPP
Function for VTEC modeling	Spherical harmonic function (4 orders and 4 degrees, 1-h interval)
Weighting strategies	Elevation-dependent weighting
Separation of SDCB and RDCB	Zero-mean condition
Estimated DCB types	GPS: C1WC2W GLO: C1PC2P GAL: C1XC5X BDS: C2IC6I
3. SFPPP processing
Position	Kinematic
Signal selection	BDS2/BDS3: B1
Clock and orbits	GFZ products
Troposphere modeling	Saastamoinen model and random walk process
Elevation cutoff	15°
Sampling interval	30 s
Ionosphere delay	IC-SFPPP: Corrected IW-SFPPP: Slant Delay + Constraint UC-SFPPP: Estimated
Satellite DCB	Corrected

summarizes the specific processing strategies of the three steps. Since PPP users are more likely to be in realistic environments where a low-elevation multipath is present, we use an elevation cut-off angle of 15 degrees on the user side.

3.2. Validation of Regional VTEC Modeling

3.2.1. Regional VTEC Comparison with CODE-GIM

Based on the estimated ionospheric coefficients, the gridded RIMs over Australia can be generated hourly. First, the reliability of the RIMs was validated and assessed by comparing with the CODE-GIMs. Six of these comparison maps was selected for presentation during the experiment, which covered the entire day. Figure 3 shows the difference in VTEC between RIM4 and CODE-GIM with the RMS statistics shown at the top of each panel. As illustrated, the TECs calculated by RIM4 conforms well to the CODE-GIMs, with most biases within two TECUs, except for the marginal areas where IPPs are insufficient.

Next, we present the RMS time series of VTEC differences between all RIMs and the CODE-GIM in Figure 4. The RMS statistics of RIMs from BDS3-only (RIMC), GPS-only (RIMG), GPS/BDS-3 combined (RIM2), and four-system combined solutions (RIM4) are indicated by purple, yellow, red, and blue dotted lines, respectively. It can be seen that the four RIMs are all in good agreement with the CODE-GIMs with the RMS of VTEC differences within two TECUs. The RIMC has the largest deviation with an average RMS of 1.081 TECUs, whereas the corresponding RMS of RIMG is 0.932 TECUs. As the CODE-GIM is generated from a combination of GPS and GLONASS observations, the RMS of RIM2 (0.862 TECUs) and RIM4 (0.867 TECUs) are comparable and better than the RIMC and RIMG. Considering that the final products from CODE have a nominal accuracy of two–eight TECUs, it can be concluded that the method of the RIMs with the quad-constellation raw observations performs well and can enhance the reliability of the VTEC model in comparison to single-constellation RIMs.

To further compare the accuracy of the four RIMs, the VTEC results for each ionospheric model at local time 12:00 (hh:mm) (a) and 22:00 (hh:mm) (b) for the central longitude of Australia (135°E) are shown in Figure 5. We can clearly see that the VTEC values decrease with the increase of latitude, and that the VTEC values at noon are greater than that at night. Since the ionosphere is more active at noon during the day, it can be seen that the difference between RIMs is more pronounced at 12:00 (hh:mm) than 22:00 (hh:mm), especially at the edges of experiment area. For comparison purposes, the locations of the reference stations are also depicted in the diagram (grey triangles). It is clear that in the area of dense station distribution as shown between two vertical black dotted lines (−12°S to −41°S), the VTEC values are essentially the same for the different models. Combining Figure 4 and Figure 5, it is apparent that the accuracy of RIMG and RIMC is already relatively high with more stations, and that the improvement in performance of the RIM4 is mainly at the edges of the survey area and over the ocean. In addition, due to its high volume of visible satellites in the Asia-Pacific region, the BDS has ionospheric monitoring capabilities equivalent to the GPS over Australia. Therefore, monitoring the ionospheric variations with BDS satellites is also an alternative when the number and signal of GPS satellites is poor.

3.2.2. Assessment of Correction Accuracy and SDCB Stability for RIMs

The ionospheric correction percentages of RIMs at each monitoring stations are calculated relative to the ‘true’ TEC. The correction percentage ‘PER’ in this paper can be defined as [50]:

$$\mathrm{PER}=\left(1-\frac{1}{N}{\displaystyle \sum _{i=1}^N\frac{\mathrm{abs}\left({\mathrm{TEC}}_{\mathrm{RIM} }^i-{\mathrm{TEC}}_{\mathrm{true}}^i\right)}{{\mathrm{TEC}}_{\mathrm{true}}}}\right)\cdot 100\%$$

(12)

where ${\mathrm{TEC}}_{\mathrm{RIM} }^i$ denotes the calculated TEC values from RIMs, whereas ${\mathrm{TEC}}_{\mathrm{true}}^i$ denotes the reference ‘true’ TEC values from CODE-GIM. It should be noted that when the absolute value of the deviation of the calculated values from the reference values exceeds the ${\mathrm{TEC}}_{\mathrm{true}}^i$, the $\mathrm{PER}$ is defined as zero. The statistics of $\mathrm{PER}$ and RMS for four RIMs results are listed in

Table 4. PER and RMS of RIMs compared to the CODE-GIM at monitoring stations.

VTEC Models	PER (%)		RMS (TECU)
	Mean	Min	Mean	Max
RIM4	91.19	85.15	0.685	1.078
RIM2	91.07	84.56	0.694	1.082
RIMG	89.00	81.42	0.817	1.487
RIMC	86.58	78.51	1.058	1.832

. We can find that more constellation observations would contribute to the ionospheric modelling, that is, RIM4 offers optimal percentage (91.19%) and correction accuracy (0.685 TECUs), whereas RIMC offers worst amendment percentage (86.58%) and correction accuracy (1.058 TECUs).

As the by-products of ionospheric VTEC modelling, the stability of SDCBs can be considered to reflect the reliability of the RIMs and can be defined as:

$$S^i=\sqrt{\frac{{\displaystyle \sum _{d=1}^D{\left(SDC{B}^{d,i}-\overline{SDC{B}^i}\right)}^2}}{D-1}}$$

(13)

where $S^i$ denotes the stability of the SDCB for satellite $i$, and $SDC{B}^{d,i}$ represents the estimated SDCB value of satellite $i$ at day $d$, whereas $\overline{SDC{B}^i}$ denotes the mean value of the SDCB of satellite $i$ over the experiment period.

We counted the stability of SDCBs from four processing modes, which are listed in

Table 5. The stability of SDCBs in different VTEC models (ns).

SDCB	Stability of SDCBs
	Single	GC	GREC
G	0.1391	0.1348	0.1321
C	0.1629	0.1553	0.1540

. It can be noted that the stability of the GPS SDCBs is better than the BDS SDCBs in the single-constellation mode. With the additional observations, the stability of SDCBs are improved. Among them, RIM4 provides the best stability of SDCBs, both for BDS (0.1540 ns) and GPS (0.1321 ns) satellites.

3.3. Validation of BDS SFPPP with RIMs

3.3.1. Assessment of Positioning Performance in BDS IC-SFPPP

In this section, we apply the generated four RIMs, as well as CODE-GIM, to the BDS kinematic IC-SFPPP and comparatively evaluate the positioning performance of different IC-SFPPP results. The datasets on DOY 184, 2021, from two stations, BULA (−22.91°S, 139.90°E) and JAB2 (12.66°S, 132.89°E), are selected as the user stations. The former station is in the central region of Australia, whereas the latter one is at the edge. The true coordinates of two stations are given by the average coordinates of the static DFPPP results over 1 day. The time series of positioning errors in five different SFPPP modes are shown in Figure 6, with statistics information shown at the top. Firstly, it is shown that the positioning accuracy of SFPPP with CODE-GIM correction are worse than the other four cases, and the RMS at BULA is 0.46 m, 0.33 m, and 0.72 m in the north, east, and upwards directions, respectively. Overall, the time series of positioning errors while introducing the RIMs into the SFPPP model are more precise and the horizontal performance is better than the vertical one. The statistics information shows that the different RIM corrections in the IC-SFPPP model have little effect on the positioning performance at BULA, whereas the 3D RMS of RIM4, RIM2, RIMG, and RIMC is 0.471 m, 0.461 m, 0.527 m, and 0.518 m, respectively. That is to say that the RIMG and RIMC can already provide ionospheric information with high spatial resolution to serve SFPPP users with decimeter accuracy, whereas the improvement of the positioning accuracy by using RIM4 and RIM2 is only several centimeters in each direction. In contrast at JAB2, the 3D RMS of the RIM4 (0.574 m) is significantly better than the RIMG (1.087 m) and the RIMC (1.007 m), especially in the upwards direction. Additionally, the errors in RIMG and RIMC increase abruptly at some moments (marked with red arrows in the figure), which are due to fewer IPPs at the edge in the regional ionospheric modelling. Therefore, it is concluded that RIM with multi-constellation observations can effectively increase the number of IPPs and compensate for the lack of accuracy for the RIMs at the edges.

The positioning errors at all monitoring stations are then calculated while introducing CODE-GIM and RIM4 into the SFPPP model. The distributions of positioning errors in the three directions are shown in Figure 7 with statistics information shown at the top. It is shown that the STD of IC-SFPPP with CODE-GIM correction are 0.21 m, 0.26 m, and 0.70 m in north, east, and upwards direction, respectively, whereas IC-SFPPP with RIM4 correction performs better than that with CODE-GIM, with STD of 0.17 m, 0.21 m, and 0.54 m in the three directions. In the upwards direction, the mean errors and STDs are improved from 0.43 m and 0.70 m to 0.22 m and 0.54 m, respectively. In general, the IC-SFPPP with the RIM4 has a better precision and better positioning performance.

Then, according to the CODE GIM products, the mean VTEC values of all grid points over Australia (0°S–50°S, 110°E–160°E) are calculated and shown in Figure 8. As shown, the VTEC values fluctuate smoothly throughout 2021, and increase sharply in 2022, indicating that the ionosphere is more active in 2022. To further evaluate the application effectiveness of the RIMs in different ionospheric periods, observations from DOY 183–189, 2021 (Period I), and DOY 77–83, 2022 (Period II), are collected and utilized in the IC-SFPPP model in the following.

Table 6. RMS statistics of the BDS IC-SFPPP (cm) and percentage increase (%) compared to the CODE-GIM for the Period I and Period II.

	CODE	GREC		GC		G		C
		RMS	Imp.	RMS	Imp.	RMS	Imp.	RMS	Imp.
Period I:
Horizontal	55.92	46.84	16.24%	46.10	17.56%	50.34	9.98%	48.44	13.38%
Vertical	79.52	66.41	16.49%	67.53	15.08%	70.41	11.46%	68.89	13.37%
3D	100.47	82.00	18.38%	82.61	17.78%	87.90	12.51%	85.15	15.25%
Period II:
Horizontal	57.12	46.79	18.08%	46.38	18.80%	48.01	15.95%	50.23	12.06%
Vertical	97.63	56.85	41.77%	58.89	39.68%	67.65	30.71%	63.22	35.25%
3D	115.92	74.62	35.36%	76.02	34.42%	84.27	27.30%	83.17	28.25%

presents the RMS and percentage improvement compared to the CODE-GIM in the horizontal, vertical, and 3D components for both periods.

It is shown in Table 6 that during Period I, the RIMC achieved a positioning accuracy of 48.44 cm and 68.89 cm in the horizontal and vertical directions respectively, while improving by 13.38% and 13.37% compared to the CODE-GIMs. With its coupling to the BDS observations, the RIMC is better than the RIMG in terms of positioning performance. Considering RIM4 and RIM2 is generated with more redundant observations, the positioning performance of BDS IC-SFPPP with the RIM4 and RIM2 is better than the RIMG and RIMC. As for the RIM2, its 3D RMS has improved by 17.78%, which is comparable to the RIM4 (18.38%). This demonstrates that the redundant observations from the dual-constellation are already adequate for high precision regional ionospheric modelling applications. In addition, we find that in the Period I, when the ionosphere is less active, the improvement of less than 20% in the RIMs is limited, whereas during Period II, the contribution of RIMs to SFPPP increases substantially. The 3D RMS for each mode is 74.62 cm, 76.02 cm, 84.27 cm, and 83.17 cm, with an improvement of 35.36%, 34.42%, 27.30%, and 28.25% respectively.

3.3.2. Assessment of Convergence Performance in BDS IW-SFPPP

In this section, BDS IW-SFPPP is utilized to comparatively evaluate the RIMs applications in the terms of positioning accuracy and convergence time. The observations from all 30 reference stations and 10 monitoring stations are collected in the analysis. Furthermore, the convergence time is defined as when the positioning errors in the horizontal and vertical directions are less than 0.3 m and 0.6 m, respectively [51], for 20 consecutive epochs. Figure 9 shows the positioning error of the horizontal (top) and upwards (bottom) directions at station STHG on DOY 184, 2021, whereas the UC-SFPPP results and corresponding statistics are also shown for comparison. It is shown that the convergence time of UC-SFPPP in the horizontal and vertical direction is 43.0 min and 44.5 min, respectively, which is reduced to 23.5 min and 14 min by using the RIM4 in the IW-SFPPP model. Undoubtedly, RIM4 play a crucial role in the early filter processing with the adoption of additional ionospheric pseudo-observations, while its impact on the positioning accuracy is weakened slightly after convergence.

To further compare the application performance of four different RIMs, we divide each RINEX observation file over 24 h into six periods and evaluate the convergence performance every 4 h to ensure that the filtering has enough time to converge to the threshold value. Then, the IW-SFPPP performance in terms of convergence time is evaluated at the median and 68% quantiles in kinematic mode [3]. Figure 10 and

Table 7. Convergence times (minimum) of kinematic SFPPP positioning errors (BDS2/BDS3).

Ionospheric Correction	Median Convergence Time			68% Quantiles Convergence Time
	N	E	U	N	E	U
No Correction (UC)	22.5	52.0	41.5	37.0	72.5	62.0
RIM4	9.5	25.5	22.5	17.5	45.5	39.0
RIM2	9.5	27	22.5	18.0	44.5	41.5
RIMG	8.5	29.0	23.0	18.0	47.5	41.0
RIMC	9.0	28.5	24.5	18.5	47.0	41.5

presents the statistics of 1680 samples (7 days × 40 stations × 6 initializations) in three directions. Obviously, both the median and 68% quantile results indicate that the convergence time of the BDS UC-SFPPP can be accelerated using the RIMs in the IW-SFPPP model, and RIM4 shows the fastest convergence overall. Among the median results, the time reduction by using the RIM4 is 13.0 min (57.8%), 26.5 min (51.0%), and 19 min (45.8%) in the north, east, and upwards direction, respectively, whereas at the 68% confidence level, the convergence time was reduced by 52.7% (from 37 min to 17.5 min), 37.2% (from 72.5 min to 45.5 min), and 37.1% (from 62.0 min to 39.0 min) compared with the UC-SFPPP performance. Additionally, the RIMC performs at the same level as the RIMG. This suggests that the accuracy of the RIM generated by a single system is sufficient for accelerating the convergence of the IW-SFPPP compared with the UC-SFPPP, whereas the more constellation observations used for RIM products, the faster the convergence time of IW-SFPPP.

4. Discussion

In this paper, GNSS observations from 30 reference stations and 10 monitoring stations in the GA network were collected. A regional ionospheric map over Australia can be generated hourly based on the estimated VTEC coefficients. The reliability of the RIMs produced in four different processing modes has been investigated by comparison with the final CODE GIM products. Firstly, in terms of VTEC differences, the four RIMs show good agreement with CODE-GIMs, with RMS values of less than 2 TECUs during the experiment. The mean RMS of VTEC differences for RIM4, RIM2, RIMG, and RIMC are 0.867, 0.862, 0.932, and 1.081 TECUs, respectively. The VTEC values for the different modes are generally consistent in the center of experiment, whereas they vary considerably at the edge of the experiment area, especially between local time 12:00 (hh:mm) and 14:00 (hh:mm) when the ionosphere is more active. Then, the ionospheric correction percentages of RIMs are calculated relative to the ‘true’ TEC, which is defined as CODE-GIM products. The ionospheric correction percentages of RIMs at each monitoring station are verified, and an average amendment percentage of 91.19% could be achieved in RIM4. In addition to the coefficients of the VTEC model, the SDCBs are also estimated in this study, and the experiments show that RIM4 can improve the SDCBs, with STDs of 0.1321 ns and 0.1540 ns for BDS and GPS satellites, respectively.

We applied the generated RIMs to the IC-SFPPP and IW-SFPPP, to validate the positioning performance of SFPPP when introducing external RIM information. The BDS IC-SFPPP experiments have demonstrated that integrating multi-GNSS data in the ionosphere modelling can improve the distribution of IPPs, especially for the stations located at the edges (e.g., JAB2). The RIM4 brings more redundant observations and its positioning performance in BDS IC-SFPPP is better than the corresponding results for a single-constellation solution. Taking the CODE-GIMs as a reference, the 3D positioning RMS with application of RIM4, RIM2, RIMG, and RIMC is 82.00 cm, 82.61 cm, 87.90 cm, and 85.15 cm, respectively, resulting in an improvement of 18.38%, 17.78%, 12.51%, and 15.25% during the less active ionospheric period, whereas during an active ionospheric period, the 3D RMS for each mode is 74.62 cm, 76.02 cm, 84.27 cm, and 83.17 cm, with an improvement of 35.36%, 34.42%, 27.30%, and 28.25%, respectively. Then, we apply the RIMs to BDS IW-SFPPP, and if the criterion of the convergence time is defined as when positioning errors in the horizontal and vertical directions are less than 0.3 m and 0.6 m for 20 consecutive epochs, the IW-SFPPP can significantly speed up the convergence time compared to the uncombined SFPPP; that is, at the 68% confidence level, the convergence time of the RIM4 is reduced by 52.7% (from 37 min to 17.5 min), 37.2% (from 72.5 min to 45.5 min), and 37.1% (from 62.0 min to 39.0 min) in the north, east, and upwards directions, respectively.

5. Conclusions

With the increasing abundance of GNSS data, further application requirements for ionospheric VTEC modeling and delay correction services are put forward. Considering the current status of multi-GNSS development, we have integrated GPS, GLONASS, BDS, and Galileo observations in regional ionospheric modelling. Additionally, a general DFPPP method is employed to extract multi-GNSS SIDs, replacing the traditional CCL model. The method improves the robustness and reliability of the parameter estimates by means of least squares, allowing for better integration of the quad-constellation raw observations. In this paper, the derived ionospheric information is used to describe the spatial and temporal variability of VTEC over Australia based on the low-order spherical harmonic function, along with DCBs that are estimated as by-products.

The experimental results show that by fusing multi-GNSS observations, users can obtain better maps of the regional ionosphere, especially for marginal regions. In addition, we applied the generated RIMs to the BDS IC-SFPPP and IW-SFPPP. Taking the CODE-GIMs as a reference, the positioning performance of IC-SFPPP with RIM4, RIM2, RIMG, and RIMC improved by 18.38%, 17.78%, 12.51%, and 15.25%, respectively. This enhancement is even more pronounced during an active ionospheric period, with an increase of 35.36%, 34.42%, 27.30%, and 28.25%, respectively. For IW-SFPPP, using the ionospheric information calculated by RIM as pseudo-observations with additional time-varying weights can effectively reduce the convergence time. Accordingly, users can obtain the regional ionospheric products based on quad-constellation raw observations, which better reflects the fine-scale variability of the regional ionosphere and achieves a better positioning performance in the SFPPP.