Abstract
Network real-time kinematic (NRTK) coverage is defined as the area inside a station network. In conventional NRTK, the distance between stations is limited to 100 km, thus restricting the coverage of NRTK. In this study, we propose the utilization of an ionospheric-free combination and the application of a Kriging weighting model to mitigate tropospheric delay to extend the coverage of NRTK through network expansion. A network with station distances exceeding 100 km was constructed, and the residual errors, along with the success-fix rate of integer ambiguities, were analyzed on both sunny and rainy days to confirm the potential for network expansion using the proposed method. The results confirm that the success-fix rate increased by up to 44.3% on rainy days, compared with that of the traditional interpolation method. Furthermore, a high level of performance in integer ambiguity resolution can be maintained within the expanded network, regardless of the weather conditions.
- ambiguity resolution
- ionospheric-free combination
- kriging-based weighting model
- network real-time kinematics (RTK)
1 INTRODUCTION
Network real-time kinematics (NRTK), a centimeter-level precise positioning technique based on multiple reference stations, has been proposed to expand the coverage of RTK (Wubbena & Bagge, 2002; Rizos & Han, 2003; El-Mowafy, 2012). NRTK approaches are classified into virtual reference station, Flachen Korrektur parameter, and master-auxiliary concept (MAC) approaches, according to the method for generating corrections (Park & Kee, 2010). Each NRTK method provides an NRTK correction that can be used to further remove the remaining measurement errors after the RTK correction from the master station has been applied. Therefore, NRTK users can achieve precise centimeter-level positioning at a greater distance from the master station compared with RTK, for which the distance is limited to 20 km. However, the positioning performance of NRTK deteriorates as the distance between reference stations increases. Therefore, in a network, it is recommended that reference stations be located at a distance of less than 100 km. Because the coverage of NRTK is defined inside the network, the NRTK coverage is limited.
Positioning performance deteriorates as the distance between the NRTK reference stations increases because of the increase in residual error. Residual errors are primarily caused by ionospheric delay, tropospheric delay, and satellite orbit errors. These errors are spatial decorrelation errors, and their magnitude varies depending on the locations at which the global navigation satellite system (GNSS) signal is collected. Among these errors, ionospheric delay exhibits the largest variation (Park, 2008) and is thus the most significant limiting factor in network expansion (Grejner-Brzezinska et al., 2005; Wielgosz et al., 2005). However, with the realization of dual-frequency measurements, users can eliminate ionospheric delays using an ionospheric-free (IF) combination. Accordingly, an NRTK using an IF combination has been proposed (Zhang et al., 2009; Song, 2016; Song et al., 2016). After ionospheric delay, the next factor limiting network expansion is tropospheric delay (Feng & Li, 2008).
In NRTK based on the MAC approach, a MAC correction is generated for each auxiliary station to remove the remaining measurement error after the carrier-phase correction (CPC) generated at the master station has been applied. With the use of dual-frequency measurements, MACs can be separated into dispersive MACs (which include ionospheric delay) and nondispersive MACs (which consist of tropospheric delay and satellite orbit error). Unlike ionospheric delay, which can be removed through an IF combination, one must compensate for the tropospheric delay by applying nondispersive MACs. Users combine and apply nondispersive MACs from auxiliary stations through an interpolation method, and the choice of interpolation method determines the positioning performance (Euler et al., 2001, 2004; Zebhauser et al., 2002). Traditional interpolation methods include the distance-based linear interpolation method (DIM), low-order surface model (LSM), linear interpolation method (LIM), and linear combination model (LCM) (Gao et al., 1997; Wanninger, 1995; Han & Rizos, 1996; Wubbenna, 1996). Reportedly, the performances of the LSM, LIM, and LCM are similar to each other and superior to that of the DIM (Dai et al, 2003). When the network is expanded, simply extending the coverage of NRTK through an IF combination is feasible if traditional interpolation methods can adequately compensate for the tropospheric delay of users. However, even within a network of conventional size, the tropospheric delay may exhibit irregular variations if changes in regional weather conditions occur. These variations lead to a decrease in positioning performance when traditional interpolation methods are used, as these methods assume that tropospheric delay within the network is planar (Shin et al., 2014). As the network size increases, the likelihood of including regional weather condition changes increases. Therefore, it is necessary to use interpolation methods that can account for irregular variations in tropospheric delay. When irregular variations of errors arise within the network, interpolation methods based on the Kriging method can improve accuracy compared with traditional methods (Geisler, 2006). The Kriging method has been used to interpolate ionospheric and tropospheric delays within GNSS reference station networks, particularly when there are significant variations in the ionosphere and troposphere within the network (Blanch et al., 2004; Zheng & Feng, 2005; Ma et al., 2020; Darugna et al., 2021). Additionally, the Kriging method has been utilized to calculate weightings for MAC application in NRTK, resulting in the proposal of the Kriging weighting model. This model has shown improved performance compared with traditional interpolation methods, particularly in environments with bias in tropospheric delay (Al-Shaery et al., 2010, 2011). Furthermore, it has been found that using the Kriging weighting model can enhance user correction accuracy, especially at network boundaries and for low-elevation satellites (Kim et al., 2017).
With the assumption of a constant altitude within the network, it is expected that the Kriging weighting model will outperform traditional interpolation methods when regional weather condition changes cause irregular changes in the troposphere. Furthermore, as the network size increases, there is a higher likelihood of including regional weather condition changes, making the Kriging weighting model an appropriate interpolation method for expanding the coverage of NRTK. Therefore, in this study, the positioning performance of NRTK in an expanded network was analyzed for an IF combination and the Kriging weighting model proposed by Kim et al. (2017), confirming the potential for expanding coverage. For the Kriging weighting model, variograms of tropospheric delay are required. A variogram represents the degree of spatial dependence of a random spatial field and is defined as the variance of difference between field values at two locations in the random spatial field. Herein, first, a method for modeling variograms using meteorological data from weather monitoring stations around the analysis area is described, and modeling results are presented. Next, a network is constructed such that the distance between reference stations is more than 100 km. Subsequently, the user residual error is analyzed, and the performance of integer ambiguity resolution is evaluated with a comparison to results obtained via the LSM.
The remainder of this paper is organized as follows. Section 2 introduces the NRTK correction generation algorithm, the user algorithm with an IF combination, and the Kriging weighting model used in this study. Section 3 introduces the method for modeling variograms using the Kriging weighting model and provides modeling results. Section 4 presents feasibility test results for the proposed method, obtained via actual measurements. Finally, Section 5 presents the conclusions of this study.
2 NRTK ALGORITHMS
In this section, the NRTK algorithm for generating corrections is introduced. Next, the user algorithms are explained. Finally, an interpolation method for combining the MACs is introduced.
2.1 Correction Generation Algorithm
The NRTK correction in the MAC method consists of the CPC of the master station and the MAC of the auxiliary stations. The CPC is generated from carrier-phase measurements at the master station. The CPC is calculated according to Equation (1):
1
where the superscript j refers to the j-th satellite; ϕ refers to the carrier-phase measurement; 
represents the distance between the broadcast satellite orbit and reference station; 
is a clock correction in the ephemeris; δb is the residual clock error after the clock correction has been applied; I, T, δR, B, and N refer to ionospheric delay, tropospheric delay, satellite orbit error, receiver clock error, and integer ambiguity, respectively; λ refers to the wavelength; and εϕ refers to the noise of the carrier-phase measurement. The MAC is created as the difference in the CPC between the master and auxiliary stations. The MAC for the k-th auxiliary station is calculated as follows:
2
where 
is the MAC correction between the k-th auxiliary station and the master station for the j-th satellite; 
indicates the difference between stations for the j-th satellite; and 
indicates the double difference between satellites and stations. The double-differenced integer ambiguities between the reference stations are determined and removed, in addition to the CPC difference. If the MAC for the L1/L2 dual frequency is generated, the MAC of the two frequencies can be combined to separate the MAC into dispersive and nondispersive MACs, which are expressed in Equations (3) and (4), respectively:
3
4
where 
represents the ratio of the squares of the frequencies of L1 (f1) and L2 (f2).
2.2 User Algorithm with an IF Combination
To avoid issues such as float ambiguity and wavelength reductions, the user algorithm is structured in two steps rather than a simple application of the IF combination of carrier-phase measurements. Figure 1 shows a flow chart of the user algorithm.
Flow chart of the user algorithm
First, the Melbourne–Wubbena (MW) combination is used to determine wide-lane (WL) integer ambiguities. Next, the integer ambiguities of the IF combination are converted into single-frequency integer ambiguities via predetermined WL integer ambiguities (Song, 2016). Users generate double-differenced pseudorange measurements and carrier-phase measurements by applying a pseudorange correction (PRC) generated in the same manner as the CPC. Equation (5) shows the result obtained by applying the PRC and CPC:
5
where 
indicates the difference between satellites for the user. MW combinations are generated from the WL combinations of the carrier-phase measurements and narrow-lane (NL) combinations of the pseudorange measurements calculated using Equation (6):
6
where c is the speed of light. If the multipath error is ignored, the MW combinations include only WL integer ambiguities and measurement noise. After the WL integer ambiguities are determined, IF combinations are generated with GPS L1/L2 carrier-phase measurements, using the following equation:
7
As shown in Equation (7), the integer ambiguities of the IF combinations are converted into those of the L2 frequency measurements based on the predetermined WL integer ambiguity. Next, the tropospheric delay is mitigated through the application of a nondispersive MAC. Users apply the interpolation method to combine the nondispersive MACs of the auxiliary stations. Before the nondispersive MACs are combined, given that altitude affects the tropospheric delay, an altitude correction of the tropospheric delay is applied (Pu et al., 2021):
8
In Equation (8), 
denotes the tropospheric delay calculated by the University of New Brunswick UNB3 model at the location of the auxiliary station, and 
denotes the tropospheric delay calculated for the location of the auxiliary station but with the altitude of the master station (hM) (Leandro et al., 2006, 2008):
9
Equation (9) presents the expression for the combined nondispersive MAC at the user location, where w refers to the weighting value. Because the traditional interpolation method calculates weightings based on the locations of the user and the reference station, all satellites have the same weightings. However, because the Kriging weighting model, introduced in the next section, has different weightings for each satellite, these weightings are expressed separately in Equation (9). Finally, the correction from Equation (9) is applied to the IF combination measurements of Equation (7) to determine the L2 integer ambiguities and calculate the position.
2.3 Interpolation Methods
2.3.1 Traditional Interpolation Methods
Traditional interpolation methods for NRTK include the DIM, LIM, LSM, and LCM. The DIM uses the reciprocal of the distances between the reference stations and user as weightings:
10
Equation (10) is utilized for the DIM, where l indicates the distance between users at the reference stations. Because the performances of the LIM, LSM, and LCM are known to be similar, only the LSM is introduced in this study. The LSM calculates weightings by assuming that the network’s internal error is a low-order surface. In this study, the first-order LSM was used because altitude correction of the tropospheric delay was performed using the user algorithm:
11
Equation (11) is used to compute the weightings based on the LSM; here, X indicates the eastward distance from the master station, and Y indicates the northward distance. As shown in Equations (10) and (11), the weightings of the traditional interpolation methods are determined by the location of the reference stations and user. Therefore, the weightings are the same for all satellites and do not change unless the reference stations and user locations change. Because the performance of the LSM is known to be better than that of the DIM, the LSM was used for comparison herein (Dai et al, 2003).
2.3.2 Kriging Weighting Model
The Kriging weighting model used in this study calculates weightings by applying ordinary Kriging to the vertical tropospheric delay (Kim et al., 2017). When nondispersive MACs are applied, the residual error of the tropospheric delay must be mitigated. Therefore, the weightings must satisfy Equation (12):
12
To apply ordinary Kriging to the vertical tropospheric delay, the slant delay must be converted to a vertical delay via a mapping function, as follows:
13
where TV refers to the vertical delay and F refers to the mapping function. In this study, the Niell mapping function was employed. Herein, a transformation matrix H was constructed with a mapping function using the relationship between the vertical delay and slant delay, as follows:
14
When ordinary Kriging is applied, the vertical delay is separated into a constant trend component 
and a spatial residual component 
, as follows (Webster & Oliver, 2007):
15
where 
denotes the measurement noise. The residual component and measurement noise are uncorrelated Gaussian random vectors. By substituting Equation (13) into Equation (12) and using the relationship in Equation (15), one can calculate the constraint and covariance of the weighting vector as follows:
16
17
where 〈·〉 indicates the covariance operator; C, 
, 
, c0, and ch refer to the covariances of the Kriging residual calculated from the variogram; and V represents the covariance of the measurement noise. Kriging is a type of best linear unbiased predictor, an estimator that satisfies an unbiased solution while minimizing variance. Therefore, the equation for weighting that minimizes covariance while satisfying the constraints is as follows:
18
In contrast to traditional interpolation methods, the weightings of each satellite are calculated. Even when the reference stations and user locations do not change, the weights change when the variogram changes. Thus, the weightings can be calculated by considering weather changes within the network. A variogram of the vertical tropospheric delay is required to apply the Kriging weighting model. Standard models exist for the ionosphere, such as the variogram used in the wide-area augmentation system (Blanch, 2003). However, regional changes are significant in the troposphere. Therefore, an appropriate variogram suitable for the data analysis environment is required.
3 VARIOGRAM MODELING FOR THE KRIGING WEIGHTING MODEL
This section introduces the modeling process for the variogram of vertical tropospheric delay. The variogram in this study represents the degree of spatial dependence of the vertical tropospheric delay. This section also presents variogram modeling results.
3.1 Variogram Modeling Using Meteorological Data
To generate a variogram, dense vertical tropospheric delay data for the analysis area are required. Herein, meteorological data were collected from weather stations, and vertical tropospheric delay data were obtained via the Saastamoinen model (Saastamoinen, 1972). Next, an empirical variogram was modeled.
The Automated Surface Observing Systems (ASOS) program, operated by the National Oceanic and Atmospheric Administration (NOAA)’s national weather service, obtains meteorological data from weather stations across the United States. In this study, because the NRTK network was constructed from reference stations located in North Carolina and Tennessee, data were collected from the ASOS weather stations in those and surrounding areas. Figure 2 shows the distribution of the weather stations used in the analysis.
Station map for variogram modeling
The weather stations provided meteorological data at 5-min intervals. The meteorological data included sea level pressure, temperature, and relative humidity. To calculate the vertical tropospheric delay using the Saastamoinen model, the meteorological data must be processed. For this, the temperature in Celsius must be converted to absolute temperature, as follows:
19
where τK denotes the absolute temperature at the station location and τc denotes the temperature in Celsius. To calculate the wet delay, the dew point temperature was calculated using the Magnus equation, as follows:
20
where Rh refers to the relative humidity and τdew denotes the dew-point temperature at the station location (Bolton, 1980). Next, we consider the following equation:
21
In Equation (21), e refers to the water vapor pressure at the station location (Younes, 2016). The tropospheric delay can be calculated at each weather station location from the meteorological data calculated in this manner and the Saastamoinen model. However, because the tropospheric delay depends on the station altitude, unification to a single altitude is necessary for variogram modeling to represent the statistical characteristics of horizontal changes. Herein, the altitude was unified to 0 m above sea level to allow us to use the sea level pressure as is. The altitude conversion formula for the temperature and water vapor pressure is expressed in Equation (22) (Mendes, 1999; Askne & Nordius, 1987):
22
In Equation (22), τ0,k denotes the absolute temperature at sea level, e0 denotes the water vapor pressure at sea level, R denotes a dry air constant (R = 257.05J/kg/K), Gn denotes the gravitational acceleration (Gn = 9.806m/s), hsea denotes a station altitude above sea level, and β, α refer to the lapse rate of temperature and water vapor pressure, respectively, which can be obtained from the UNB3 model. After altitude correction, the vertical tropospheric delay was calculated via the Saastamoinen model:
23
where ZHD represents the zenith hydrostatic delay, ZWD is the zenith wet delay, ZTD is the zenith tropospheric delay, P0 denotes the sea level pressure, θ denotes the station latitude, and hell represents the ellipsoidal altitude of the station. Because the ellipsoidal altitude differs from the altitude above sea level, a conversion is required. A conversion was performed by calculating the geoid height based on a geoid model. Finally, a variogram of the vertical tropospheric delay calculated at sea level was modeled. An empirical variogram of the vertical tropospheric delay was used in this study. First, a variogram value was calculated according to the distance between two points in space, with the variogram value being equal to half the square of the difference between the vertical tropospheric delays of the two points:
24
where η refers to the variogram value, x refers to the location of one point, and lag is the distance between the two points. Using the variogram value obtained in this manner, the average 
was calculated based on the number of data points (L) in a certain lag section, and the empirical variogram was modeled by fitting the average of each section. Various fitting models can be used for the empirical variograms. Exponential fitting was used by referencing previous studies (Al-Shaery et al., 2010, 2011), as follows:
25
3.2 Variogram Modeling Results
Because the troposphere has strong regional characteristics that depend on weather conditions, variograms may vary depending on weather conditions. Therefore, in this study, variograms for 4 days with different weather conditions were modeled. Table 1 summarizes the information for these data sets.
Data Set Information
Figure 3 shows a precipitation radar image provided by NOAA. The area marked in the figure indicates the location of the NRTK reference station network.
Precipitation radar image for each data set
Figures 4 and 5 show the variogram modeling process for data sets 1 and 2. According to the procedure described above, the vertical tropospheric delay was calculated using meteorological data for the weather station locations, as shown in Figures 4(a) and 5(a). Afterward, the average variogram value was calculated according to Equation (24). The variogram values were categorized into bins with a width of 20 km over the lag range, with the maximum lag range set to 300 km to ensure sufficient data for each bin. Thus, average variogram values were computed for each bin. The calculated average of variogram values underwent exponential fitting via Equation (25) to obtain the nugget, sill, and range parameters. Figures 4(b) and 5(b) illustrate both the average variogram values and the results of exponential fitting.
Empirical variogram modeling of data set 1: (a) station map with ZTD; (b) modeling result
Empirical variogram modeling of data set 2: (a) station map with ZTD; (b) modeling result
Figure 6 shows the variograms for each data set. Table 2 summarizes the coefficients for each variogram. Data sets 1 and 2 were acquired only 2 days apart, and the analysis time was set such that the satellite environments were the same. These two data sets are different in terms of their weather conditions, and the obtained variograms for the two data sets behaved differently depending on the weather conditions. By contrast, the variograms of data sets 2, 3, and 4 are similar, even though they present a time difference of approximately 1 year. Thus, the variograms of tropospheric delay can vary, depending on weather conditions. Consequently, to apply the Kriging weighting model, a variogram appropriate for the analysis date is required. However, given that similar variograms can be obtained for similar weather conditions, if a representative variogram for each weather condition is provided, users can apply the variogram directly according to their weather conditions when using the Kriging weighting model.
Comparison of empirical variograms for each data set
Coefficients of Empirical Variograms
4 NRTK TEST RESULTS
To analyze the NRTK positioning performance when the IF combination and Kriging weighting model are applied, the residual error of the double-differenced carrier-phase measurements was analyzed. Subsequently, the performance of the integer ambiguity resolution was analyzed via the least-squares ambiguity decorrelation adjustment (LAMBDA) (Teunissen, 1993, 1995; Teunissen et al., 1997). The results were compared with those obtained by the LSM.
4.1 Network Setting
A reference station network was constructed such that the distance between stations was more than 100 km. In Figure 7, the reference stations depicted by black triangles constitute the network, whereas the reference stations depicted by red triangles are regarded as users inside the network. Table 3 summarizes the information for each reference station and its distance from the master reference station.
NRTK station map
Reference Station and User Information
To minimize the effect of equipment differences, reference stations with the same Trimble receiver and antenna were used. The residual errors of NCSW and GAST, which are the closest to and farthest from the master station, respectively, are presented, along with the statistical values for all other users.
4.2 Comparison Between Sunny and Rainy Days
In this study, because the ionospheric delay error was removed via the IF combination, the tropospheric delay accounted for most of the residual error. As mentioned earlier, because tropospheric delay is affected by weather conditions, data sets 1 and 2, which have weather condition changes on nearby days, were compared and analyzed. Because the satellite geometry was set to be the same by adjusting the time, only the effects of changes in the tropospheric delay error could be analyzed. The WL integer ambiguities were determined through the MW combination. By employing the MW combination, both ionospheric and tropospheric delays are removed, thereby maintaining the performance of WL integer ambiguity resolution even with network expansion. Consequently, with the assumption that WL integer ambiguities are resolved, the residuals of the IF combination and the performance of the N2 integer ambiguity resolution were analyzed. In addition, we compare the residual error of single-baseline RTK, NRTK with a conventional weighting model such as the LSM, and NRTK with Kriging weighting applied. Because the reference stations located within the network were set as users, the true positions of the users are known. Therefore, the residual error for single-baseline RTK was calculated from the true position and integer ambiguities predetermined through the batch process:
26
We note that the term 
is computed based on the true position of the user, ensuring that the residual error computed by Equation (26) is primarily attributed to tropospheric residual errors after tropospheric correction has been applied. Equation (26) presents the residual error of the IF combination when only the CPC of the master station is applied. The distance term and integer ambiguity are removed. To assess variations attributable to weather conditions, altitude correction was employed to mitigate the impact of tropospheric delay discrepancies resulting from differences in altitude. However, because this correction addresses the same residual error as the identical correction in single-baseline RTK, 
is defined as the residual error of single-baseline RTK. The residual errors of single-baseline RTK are shown in Figures 8 and 9, where the mask angle was set to 20°.
Residual error of single-baseline RTK (5/28/2022, sunny): (a) NCSW; (b) GAST
Residual error of single-baseline RTK (5/26/2022, rainy): (a) NCSW; (b) GAST
In Figures 8 and 9, the red line represents half of the NL wavelength (5.35 cm). A comparison of the residual errors of NCSW and GAST, which are the closest and farthest stations from the master station, respectively, revealed that the residual error increases with increasing distance. In addition, the residual error is greater on rainy days. One must compensate for these residual errors by applying a nondispersive MAC:
27
Equation (27) is the residual error of NRTK when the correction defined in Equation (9) is applied to the residual error of single-baseline RTK. This correction varies depending on the interpolation method used. Therefore, the residual error was analyzed according to the interpolation method. The residual errors of NRTK for a sunny day are shown in Figures 10 and 11.
Residual error of NRTK for NCSW (5/28/2022, sunny): (a) LSM; (b) Kriging
Residual error of NRTK for GAST (5/28/2022, sunny): (a) LSM; (b) Kriging
The residual errors decreased when the Kriging weighting model was used for NCSW. Moreover, there was a decrease in the incidence of residual errors that exceeded half of the NL wavelength. The fact that the residual error exceeds half a wavelength implies that the integer ambiguities may be incorrectly determined. Therefore, employing the Kriging weighting model can mitigate the occurrence of such incorrect integer ambiguities. In the case of GAST, the residual errors of the two interpolation methods were similar and were found to be smaller overall than those for NCSW. This trend can be viewed as an aspect that is distinct from the residual error of single-baseline RTK, as the choice of interpolation method directly impacts the residual error in NRTK. Further analysis of this aspect will be provided in Section 4.4. Table 4 presents the root mean square (RMS) of residual errors for all users within the network, organized by their distance from the master station.
Comparison of RMS of Residual Error for All Users Using NRTK with LSM and Kriging (5/28/2022, Sunny)
Now we assess the impact of the weighting method on ambiguity resolution. We have observed that the residual error can be reduced by applying the Kriging method instead of the LSM. To isolate the impact of the residual error on the ambiguity determination, distance terms in the IF combination were removed, and float solutions were calculated. Subsequently, LAMBDA was applied epoch by epoch:
28
Considering the noise of IF combination measurements, float solutions were computed by using weighted least squares (WLS):
29
When LAMBDA is applied to the float solutions, candidate integer ambiguities can be derived. At each epoch, the ratio test was conducted. For epochs that passed the test, the integer ambiguities were determined. Among the commonly used thresholds for the ratio test, a threshold of 2/3 was used (Teunissen, 2013). Two states were defined for epochs in which integer ambiguities were determined. An epoch was defined as being in the “Correct” state when the integer ambiguities determined through the batch process and those determined through LAMBDA match, and an epoch was defined as being in the “Wrong” state if the integer ambiguity differs for any satellite. Figures 12 and 13 present the ambiguity resolution results for NCSW and GAST, respectively.
NRTK ambiguity resolution results for NCSW (5/28/2022, sunny): (a) LSM; (b) Kriging
NRTK ambiguity resolution results for GAST (5/28/2022, sunny): (a) LSM; (b) Kriging
In Figures 12 and 13, the value is 0 at each epoch when in the correct state and 1 when in the wrong state. As predicted from residual error analysis, the use of the Kriging weighting model for NCSW significantly reduced the incidence of the wrong state compared with the LSM. With the success-fix rate defined as the ratio of epochs for which integer ambiguities are correctly determined to those determined overall, it was confirmed that the success-fix rate is 95.1% for the LSM and 99.9% for the Kriging weighting model. For GAST, residual errors did not exceed half of the NL wavelength for either interpolation method. As a result, integer ambiguities were correctly determined in all epochs, and a success-fix rate of 100% was confirmed. Table 5 presents the success-fix rates of all users within the network, organized in order of distance from the master station.
Comparison of Success-Fix Rate of All Users Using NRTK with LSM and Kriging (5/28/2022, Sunny)
When the LSM is used, there is a slight decrease in the success-fix rate; however, regardless of the interpolation method used, a success-fix rate of over 95% was confirmed for all users. In Figure 12, it can be observed that epochs in which the wrong state persists do not occur continuously. Therefore, it is anticipated that applying sequential filters such as the Kalman filter in user algorithms can further increase the success-fix rate. Based on these results, it is expected that even when the IF combination and traditional interpolation methods are used on sunny days, there is potential for network expansion. However, by using the Kriging weighting model, the success-fix rate can be improved by up to 4.8%. Furthermore, because this approach achieves a success-fix rate of nearly 100% for all users, using the Kriging weighting model is advantageous for network expansion.
The fixed user position error was computed based on fixed integer ambiguities. Equation (30) represents the measurement equation, and Equation (31) indicates how the fixed user position is computed based on WLS:
30
31
During the process of ambiguity resolution, the distance terms were removed by using the true position of the user, ensuring that the integer ambiguities were determined independently of the user’s position. This approach was applied with the intention of isolating the impact of residual error on the ambiguity resolution. With this method of determining integer ambiguities, residual errors are absorbed into the errors of the integer ambiguities. Therefore, even when incorrectly fixed integer ambiguities are used, the position error may not be significantly large. For this reason, the analysis was conducted only for epochs with correctly fixed integer ambiguities.
Figures 14 and 15 present the position errors for NCSW and GAST, respectively, in the east–north–up (ENU) frame. In Figure 14, for NCSW, there exist differences in position accuracy depending on the interpolation method used, which is attributed to the residual error affecting the position accuracy after the correction has been applied. Nevertheless, when the integer ambiguities are correctly determined, centimeter-level position accuracy can be achieved. Therefore, the success-fix rate was utilized as a metric to analyze the positioning performance of NRTK users. Next, the results for a rainy day are presented in the same order as those for a sunny day.
NRTK position error for NCSW in the ENU frame (5/28/2022, sunny): (a) LSM; (b) Kriging
NRTK position error for GAST in the ENU frame (5/28/2022, sunny): (a) LSM; (b) Kriging
Figures 16 and 17 show the residual error of NRTK for a rainy day. Evidently, the residual error increased compared with that for a sunny day. Specifically, the residual error for NCSW increases when the LSM is used. Therefore, determining integer ambiguities when using the LSM was expected to be challenging. However, the residual error can be reduced by using the Kriging weighting model for NCSW. Although the Kriging weighting model was used, some epochs have a residual error that exceeds half of the NL wavelength; however, the incidence of such epochs is lower compared with that obtained when the LSM is used. Similarly, for GAST, although the residual error increases compared with the sunny day results, some epochs have a residual error that exceeds half of the NL wavelength; however, the application of both interpolation methods results in similar levels of residual error.
Residual error of NRTK for NCSW (5/26/2022, rainy): (a) LSM; (b) Kriging
Residual error of NRTK for GAST (5/26/2022, rainy): (a) LSM; (b) Kriging
Table 6 summarizes the RMS of the residual error for all users within the network for a rainy day. Overall, the residual error increased for the rainy day compared with that for a sunny day. In particular, when the LSM was used, the residual error increased by a factor of two for all reference stations. However, when the Kriging weighting model was used, the residual error increased by 20%; thus, the performance degradation was less significant compared with that of the LSM. Figures 18 and 19 present ambiguity resolution results for NCSW and GAST, respectively, on a rainy day, and Table 7 presents the success-fix rates of all users within the network, organized in order of distance from the master station.
NRTK ambiguity resolution results for NCSW (5/26/2022, rainy): (a) LSM; (b) Kriging
NRTK ambiguity resolution results for GAST (5/26/2022, rainy): (a) LSM; (b) Kriging
Comparison of RMS of the Residual Error of All Users Using NRTK with LSM and Kriging (5/26/2022, Rainy)
Comparison of Success-Fix Rate of All Users Using NRTK with LSM and Kriging (5/26/2022, Rainy)
The success-fix rate decreased significantly when the LSM was used for the rainy day. In particular, the rate varied significantly depending on the user’s location. Even when the Kriging weighting model was used, the success-fix rate decreased. Nevertheless, the success-fix rate for all users remained higher than 94% and increased by up to 44.3% compared with the LSM results. Earlier, it was confirmed that on sunny days, there is potential for network expansion when using the IF combination and traditional interpolation methods. However, during turbulent weather conditions, such as those encountered on a rainy day, it is difficult to maintain the performance of integer ambiguity resolution when using traditional interpolation methods, indicating difficulty in network expansion. Nevertheless, it was confirmed that using the Kriging weighting model with weather-appropriate variograms allows for high performance of integer ambiguity resolution, even on rainy days. These results suggest the possibility of expanding the NRTK network regardless of weather conditions by using the IF combination and Kriging weighting model.
Figures 20 and 21 display the position error for NCSW and GAST, respectively, in the ENU frame. Even on rainy days, achieving centimeter-level position accuracy is possible if only the integer ambiguities are correctly determined. However, as evident in Figure 20, when the LSM is used for NCSW, the success-fix rate is low, making it difficult to consistently achieve centimeter-level position accuracy. Ultimately, expanding the NRTK coverage requires stable centimeter-level position accuracy within the expanded network. Therefore, as proposed in this study, using the IF combination and Kriging weighting model would enable NRTK coverage expansion regardless of weather conditions.
NRTK position error for NCSW in the ENU frame (5/26/2022, rainy): (a) LSM; (b) Kriging
NRTK position error for GAST in the ENU frame (5/26/2022, rainy): (a) LSM; (b) Kriging
4.3 Analyses of Other Rainy Days
In the previous section, results for data sets 1 and 2 were compared. Because the two data sets were collected 2 days apart and maintained the same satellite geometry, the environmental factors other than weather conditions were assumed to be the same. In this section, data sets 3 and 4, which were collected on rainy days, were analyzed to assess whether the results of data set 2 in the previous section are reproducible. The analysis results of data set 3 are presented in Figures 22 and 23 and Tables 8 and 9.
Residual error for NCSW (6/12/2023, rainy): (a) single-baseline RTK; (b) LSM; (c) Kriging
Residual error for GAST (6/12/2023, rainy): (a) single-baseline RTK; (b) LSM; (c) Kriging
Comparison of RMS of Residual Error of All Users Using NRTK with LSM and Kriging (6/12/2023, Rainy)
Comparison of Success-Fix Rate of All Users Using NRTK with LSM and Kriging (6/12/2023, Rainy)
The analysis results for data set 4 are presented in Figures 24 and 25 and Tables 10 and 11.
Residual error for NCSW (6/19/2023, rainy): (a) single-baseline RTK; (b) LSM; (c) Kriging
Residual error for GAST (6/19/2023, rainy): (a) single-baseline RTK; (b) LSM; (c) Kriging
Comparison of RMS of Residual Error of All Users Using NRTK with LSM and Kriging (6/19/2023, Rainy)
Comparison of Success-Fix Rate of All Users Using NRTK with LSM and Kriging (6/19/2023, Rainy)
The results of data sets 3 and 4 were evidently similar to those of data set 2, despite the presence of different environmental factors other than weather conditions. Compared with the LSM, using the Kriging weighting model resulted in reduced residual errors and improved performance for integer ambiguity resolution. In particular, when the Kriging weighting model was used, a high success-fix rate was maintained even on rainy days within the expanded network. However, when the LSM was used, the success-fix rate decreased, with significant performance degradation observed for NCSW, the station closest to the master station. Unlike single-baseline RTK, in NRTK, the choice of interpolation method, instead of the distance to the master station, influences the residual error. Therefore, it is predicted that the performance degradation, particularly when the LSM is used and the station is close to the master station, is due to the characteristics of the LSM. This aspect will be discussed further in the next subsection.
4.4 Discussion
Table 12 summarizes the weightings of all satellites at the initial epoch for data sets 1 and 2 for NCSW.
Weightings of the Initial Epoch Calculated with the LSM and Kriging Weighting Model (Data Sets 1 and 2)
Evidently, the weightings for each satellite were different when the Kriging weighting model was used. As suggested in previous research, the Kriging weighting model considers the distribution of the vertical tropospheric delay using a mapping function, i.e., the optimal weightings are calculated for each satellite (Kim et al., 2017). Therefore, the residual error can be reduced compared with that of the LSM. Moreover, because traditional interpolation methods utilize only the locations of the reference stations and the user, changes in weather conditions cannot be included, as the weightings are fixed. By contrast, the Kriging weighting model can utilize different variograms depending on the weather conditions to calculate weightings that consider weather condition changes. Consequently, the residual error on rainy days is reduced, thus improving the integer ambiguity resolution performance. Moreover, as the network size increases, it becomes more likely that changes in internal weather conditions will be captured within the network. Thus, an interpolation method that can account for these changes is required to extend the coverage of NRTK.
To extend NRTK coverage, the error must be uniform, rather than small at only one point within the network. A uniform residual error implies a uniform performance of integer ambiguity resolution. When the LSM is used, the magnitude of the uneven residual error within the network was confirmed; in particular, the residual error increased as the location approached the master station because of the LSM constraint. As shown in Equation (11), the sum of the weightings calculated by the LSM must be unity, which is confirmed in Table 11. When the sum of the weights is unity, the influence of the master station error resulting from the MAC combination is removed. When the number of auxiliary stations is three and the sum of the weights is unity, the residual error can be calculated by expanding Equation (12), as follows:
32
As indicated in Equation (32), the error of the master station is removed, and the remaining error is expressed as the sum of the errors of the user and auxiliary stations. Therefore, when LSM weightings are applied, the residual error is interpreted as the difference between the interpolated tropospheric delay of the auxiliary stations and the user error. Therefore, when the auxiliary stations are closer, the residual error is smaller; additionally, the distribution shape of the residual error depends on the arrangement of the auxiliary stations.
Herein, as only a few reference stations were available within the network, the distribution of residual errors was analyzed by using simulation data. The distribution of the residual error of the tropospheric delay inside the network was confirmed by the satellite orbit in the initial epoch of data set 1. For this purpose, the location of the user was set on a grid. The UNB3 model was used for tropospheric delay. Because the residual error is the largest for low-elevation-angle satellites, the results of the LSM and Kriging weighting models were compared for pseudorandom number (PRN) 29 (elevation angle: 20°), the lowest-elevation satellite, assuming the worst case. Figure 26 shows the simulation settings, and Figure 27 shows the analysis results.
Simulation setting: (a) grid map; (b) tropospheric delay (PRN 29)
Residual error map: (a) LSM; (b) Kriging
The simulation analysis revealed that even when the UNB3 model, which has no regional changes in tropospheric delay, was used, the residual error distribution within the network was nonuniform when the LSM was applied. As expected, the error increased as the distance from the auxiliary stations increased. This phenomenon can also be expected when using the LIM and LCM, which are traditional interpolation methods that share the same constraints. However, when the Kriging weighting model was used, the residual error distribution was uniform. Because the actual troposphere has large regional variations, the nonuniformity of the residual error distribution is expected to increase when the LSM is used. Regardless of the weather conditions inside the NRTK network, the residual error due to tropospheric delay must be small and must satisfy the uniform performance of integer ambiguity resolution within the network. The analysis confirmed that the previous two conditions were satisfied when the Kriging weighting model was used. Therefore, the NRTK coverage can be expanded by using the Kriging weighting model.
5 CONCLUSIONS
In general, the coverage of NRTK is equal to the size of the reference station network. Therefore, to expand coverage, the network must be expanded by increasing the distance between reference stations. Conventional NRTK approaches recommend that the distance between reference stations be less than 100 km. The main cause of this distance limitation is the ionospheric delay; however, with the availability of dual-frequency measurements, this limitation can be eliminated by combining measurements. In addition to ionospheric delay, tropospheric delay also limits the expansion of network size. Unlike the ionospheric delay, the tropospheric delay cannot be removed through a combination; rather, one must compensate for the tropospheric delay by applying a correction. However, as the network expands, the likelihood of an irregular distribution of tropospheric delays occurring within the network increases because of changes in regional weather conditions. In this case, error compensation based on traditional NRTK interpolation methods is limited because these methods assume a linear variation in the tropospheric delay over the area within the network. However, because the Kriging weighting model utilizes the statistical properties of tropospheric delay, the weightings can be calculated by considering the uncertainty due to the irregular distribution within the network. For this reason, this study proposed the use of an IF combination and the Kriging weighting model to expand the NRTK coverage.
To use the Kriging weighting model, variograms of vertical tropospheric delay were modeled. Because the variogram of tropospheric delay depends on weather conditions, meteorological data were collected from weather stations in the eastern United States on sunny and rainy days. Subsequently, the ZTD at each station location was calculated based on the collected meteorological data and the Saastamoinen model. Finally, the empirical variogram was modeled via an exponential model with the ZTD of all stations when the weather conditions were maintained. The modeling revealed that the variograms of sunny and rainy days were different, whereas the variograms were similar for all three rainy days considered. Thus, when the Kriging weighting model is used, the variogram for each analysis day need not be modeled; however, a representative variogram can be used according to weather conditions.
After variogram modeling, a network with a distance of more than 100 km between reference stations was constructed, and NRTK was applied by selecting five internal reference stations as users. Analysis of the residual error of the IF combination revealed that, regardless of the weather conditions, using the Kriging weighting model resulted in smaller residual errors compared with that for the LSM, with fewer occurrences exceeding half of the NL wavelength. After the residual error had been assessed, the performance of integer ambiguity resolution was analyzed. Because centimeter-level position accuracy can be achieved when integer ambiguities are correctly determined, the success-fix rate was used as a performance metric. With LAMBDA applied epoch by epoch, it was observed that for a sunny day, both the LSM and the Kriging weighting model maintained a high success-fix rate of over 90% within the expanded network; moreover, using the Kriging weighting model resulted in an improvement of up to 4.8%. This finding indicates the potential for expanding NRTK coverage even when using the IF combination and traditional interpolation methods on a sunny day. However, on rainy days, there was a significant decrease in the success-fix rate when the LSM was used, with some users experiencing a decrease to approximately 50%. In contrast, when the Kriging weighting model was used, a high level of over 90% was maintained even on rainy days, and the success-fix rate was improved by up to 44.3%. The high success-fix rate with the Kriging weighting model confirms the possibility of achieving stable centimeter-level position accuracy. Expanding NTRK coverage requires maintaining stable centimeter-level position accuracy within the network. Therefore, it is concluded that using the IF combination and Kriging weighting model allows for expanded NRTK coverage, regardless of the weather conditions.
In this study, the use of an IF combination and Kriging weighting model was proposed to expand the coverage of NRTK. The IF combination used in this study employed GPS L1/L2 dual frequencies and determined WL integer ambiguities in advance through the MW combination, allowing for use of the NL wavelength. However, the short wavelength posed constraints on integer ambiguity resolution. By utilizing measurement combinations that maintain a longer wavelength while removing ionospheric delay errors, it is anticipated that not only can the positioning performance within the network be improved, as proposed in this study, but further coverage expansions may be possible as well. Furthermore, in this study, the reference station was selected as the user to investigate the potential for coverage expansion regardless of user algorithms. Residual error analysis and epoch-by-epoch integer ambiguity resolution were conducted using true position values. Hence, there were limitations in confirming the position error that could actually occur when integer ambiguities were incorrectly fixed. In future work, the performance of user positioning will be analyzed for internal users of the expanded network when sequential filters such as the Kalman filter are applied, to examine whether centimeter-level position accuracy can be achieved using the method proposed in this study.
HOW TO CITE THIS ARTICLE
Kim, B. -G., Kim, D., Song, J., & Kee, C. (2024). Expanding network RTK coverage using an ionospheric-free combination and kriging for tropospheric delay. NAVIGATION, 71(3). https://doi.org/10.33012/navi.662
ACKNOWLEDGMENTS
This research was supported in part by the Unmanned Vehicles Core Technology Research and Development Program through the National Research Foundation of Korea and the Unmanned Vehicle Advanced Research Center funded by the Ministry of Science and ICT, the Republic of Korea, contracted through the Seoul National University Future Innovation Institute (no. 2020M3C1C1A01086407). This research was also supported in part by the Institute of Advanced Aerospace Technology at Seoul National University. The Institute of Engineering Research at Seoul National University provided research facilities for this work.
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Jump to section
Related Articles
Cited By...
- No citing articles found.