*Article* **Modified Interpolation Method of Multi-Reference Station Tropospheric Delay Considering the Influence of Height Difference**

**Yakun Pu 1,2, Min Song 1,\* and Yunbin Yuan <sup>1</sup>**


**Abstract:** In network real-time kinematic (NRTK) positioning, atmospheric delay information is critical for generating virtual observations at a virtual reference station (VRS). The traditional linear interpolation method (LIM) is widely used to obtain the atmospheric delay correction. However, even though the conventional LIM is robust in the horizontal direction of the atmospheric error, it ignores the influence of the vertical direction, especially for the tropospheric error. If the height difference between the reference stations and the rover is large and, subsequently, tropospheric error and height are strongly correlated, the performance of the traditional method is degraded for tropospheric delay interpolation at the VRS. Therefore, considering the height difference between the reference stations and the rover, a modified linear interpolation method (MLIM) is proposed to be applied to a conventional single Delaunay triangulated network (DTN). The systematic error of the double-differenced (DD) tropospheric delay in the vertical direction is corrected first. The LIM method is then applied to interpolate the DD tropospheric delay at the VRS. In order to verify the performance of the proposed method, we used two datasets from the American NOAA continuously operating reference stations (CORS) network with significant height differences for experiments and analysis. Results show that the DD tropospheric delay interpolation accuracy obtained by the modified method is improved by 84.1% and 69.6% on average in the two experiments compared to the conventional method. This improvement is significant, especially for low elevation satellites. In rover positioning analysis, the traditional LIM has a noticeable systematic deviation in the up component. Compared to the conventional method, the positioning accuracy of the MLIM method is improved in the horizontal and vertical directions, especially in the up component. The accuracy of the up component is reduced from tens of centimeters to a few centimeters and demonstrates better positioning stability.

**Keywords:** GPS; NRTK; VRS; tropospheric delay; interpolation

#### **1. Introduction**

The traditional single-baseline real-time kinematic (RTK) positioning technology has been developed over the last several decades. However, due to the influence of atmospheric errors, orbital errors, and other distance-related errors, the distance between the rover and the reference station can only be about 10 km for centimeter-level positioning. The emergence of network real-time kinematic (NRTK), which is based on multiple reference stations, has expanded the scope of services for precise real-time positioning [1–3]. Virtual reference station (VRS) technology has been widely applied in NRTK data processing. When the ambiguities between reference stations are fixed in the network, distance-dependent errors, such as the double-differenced (DD) ionospheric and tropospheric delay of baselines, can easily be calculated. Using the estimated error among the reference stations, we can obtain

**Citation:** Pu, Y.; Song, M.; Yuan, Y. Modified Interpolation Method of Multi-Reference Station Tropospheric Delay Considering the Influence of Height Difference. *Remote Sens.* **2021**, *13*, 2994. https://doi.org/10.3390/ rs13152994

Academic Editor: Stefania Bonafoni

Received: 9 June 2021 Accepted: 26 July 2021 Published: 29 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the corresponding error between the reference stations and the VRS near the rover through the interpolation method. Finally, the generated virtual observations and the rover's actual observations can be used for short baseline relative positioning, and the centimeter-level coordinates of the rover can be quickly obtained [4,5].

After extracting the atmospheric error from the baseline of the reference stations correctly, selecting the appropriate algorithm to generate the virtual observations at the VRS becomes a crucial issue. Therefore, researchers have proposed several interpolation models. Wanninger originally proposed a linear interpolation method (LIM) that requires at least three reference stations around the rover. The dual-frequency phase observations and the known coordinates of the reference station are applied to derive the regional DD ionospheric model [6]. The LIM can also build the model of other distance-dependent errors, such as tropospheric delay [7]. Han and Rizos introduced a linear combination model (LCM) of single-differenced (SD) observations, mainly used to eliminate orbital errors and that can be applied to model atmospheric errors [8]. The distance-based linear interpolation method (DIM), proposed by Gao et al., interpolates the rover's ionospheric delay based on the distance between the reference stations and the rover [9]. Wübbena applied a suitable trend surface, called the low-order surface model (LSM), to simulate the distancedependent error trend in the network. The coefficients of the LSM model can be obtained from the observation data of the reference stations by the least-squares adjustment [10]. The least-squares collocation method can also be applied to interpolate the atmospheric delay [11–13]. A comprehensive comparison of the advantages and disadvantages of these interpolation methods was undertaken by Dai et al. [14], Fotopoulos et al. [15], Wu et al. [16], and Al-Shaery et al. [17]. All results showed that the effects of these interpolation methods are similar. It is difficult to say which one performs the best. However, The LSM requires at least four stations to resolve the model coefficients. If the network is formed by Delaunay triangulation in VRS mode, it is difficult to use LSM in the smallest triangular unit. The rover station is generally surrounded by three reference stations in the traditional Delaunay triangulated network (DTN), so the LIM method is easier to implement for a triangulation model.

With the increase of the baseline length, the correlation of tropospheric errors among the stations decreases, which leads to large residual errors in DD observations. Significant residual errors will affect the ambiguity resolution (AR) and positioning quality. Since the tropospheric delay is highly correlated with elevation, its correlation coefficient exceeds 0.9, which will seriously affect the estimation of the station positions [18–20]. Landau used the experimental data of the German SAPOS network and found that when an obvious height difference exists between the rover and reference stations, the tropospheric delay systematic error caused by the height deviation can reach 6.8 cm [21]. If the interpolation algorithms mentioned above are applied in the smallest triangular unit directly, they only consider the horizontal distribution of the tropospheric errors. The tropospheric delay interpolated at the VRS is strongly constrained in the horizontal plane that is formed by the reference stations. The tropospheric delay in the vertical direction may be significantly inaccurate. Wu et al. proposed that the multi-baseline tropospheric delay interpolation method has better accuracy than the traditional LIM in the star network [22]. In addition, based on the BP artificial neural network method, Qiu et al. used a neural network to better correct the tropospheric delay and build the spatial tropospheric delay error model [23].

The above-mentioned interpolation methods either did not consider the systematic error of the tropospheric delay caused by the height difference between the user station and the reference station or required using multiple baselines and complex modeling methods, which are not convenient in practice. In our study, we proposed an improved linear interpolation method that was simple to implement and that was effective. Taking the height difference between the reference stations and the rover into account, the prior model of the troposphere was used to correct the estimated DD tropospheric delay on the network baseline first. After that, the LIM was used to interpolate the tropospheric delay between the VRS and the master reference station. Finally, the interpolated value was applied to generate virtual observations and positioning.

The main objective of this research is to improve the accuracy of the tropospheric delay correction using the interpolation algorithm in the NRTK of the Delaunay triangulated network (DTN) model in order to obtain more stable and more accurate positioning results through the modified method, especially in the vertical direction. The structure of this paper is as follows: In Section 2, the calculation method of the DD tropospheric delay of the network baseline is introduced, and the traditional LIM and MLIM methods are presented in detail. In Section 3, using the experimental data from the triangular units, the accuracy of the tropospheric delay correction when applying both the conventional and modified methods is analyzed with different satellite elevation angles. The difference in the positioning results from the two methods is compared. The results of the analyses are then discussed in Section 4. Finally, the conclusions are drawn, and further research directions are suggested in Section 4.

#### **2. Materials and Methods**

Here, we first describe the VRS's tropospheric delay acquisition process, which includes the calculation of the DD tropospheric delay among the reference stations and the traditional modeling method of the tropospheric delay at the VRS. The modified method is then developed in Section 2.2, using the prior tropospheric model to correct the tropospheric elevation system error due to the height difference between the reference and the rover stations and combining the LIM interpolation algorithm with the calculated DD tropospheric delay to obtain the DD tropospheric delay between the VRS and the master station.

#### *2.1. Calculation and Modeling Process of the Tropospheric Delay*

In order to accurately obtain the tropospheric delay of the network baselines, we established the DD observation equations to obtain the wide-lane (WL) ambiguities, which can be expressed as follows:

$$\begin{cases} L\_{rb,WL}^{ij} = \rho\_{rb}^{ij} + T\_{rb}^{ij} + \frac{f\_1}{f\_2} (I\_{rb,1}^j - I\_{rb,1}^i) + \lambda\_{WL} (N\_{rb,WL}^j - N\_{rb,WL}^i) + \varepsilon\_{rb,WL}^{ij} \\ L\_{rb,1}^{ij} = \rho\_{rb}^{ij} + T\_{rb}^{ij} - (I\_{rb,1}^j - I\_{rb,1}^i) + \lambda\_1 (N\_{rb,1}^j - N\_{rb,1}^i) + \varepsilon\_{rb,1}^{ij} \\ P\_{rb,1}^{ij} = \rho\_{rb}^{ij} + T\_{rb}^{ij} + (I\_{rb,1}^j - I\_{rb,1}^i) + \epsilon\_{rb,1}^{ij} \\ P\_{rb,2}^{ij} = \rho\_{rb}^{ij} + T\_{rb}^{ij} + \frac{f\_2^i}{f\_2^2} (I\_{rb,1}^j - I\_{rb,1}^i) + \epsilon\_{rb,2}^{ij} \\ \overline{I\_{rb,1}^j} = I\_{rb,1}^j - I\_{rb,1}^i \end{cases} \tag{1}$$

where the subscripts *r* and *b* identify the reference stations, and the superscripts *i* and *j* denote the reference and rover satellites, respectively. *Lij rb*,*WL*, *<sup>L</sup>ij rb*,1, *<sup>P</sup>ij rb*,1, and *<sup>P</sup>ij rb*,2 represent the DD carrier and the pseudorange observations of the *WL*, and *L*1, *P*1, and *P*2 are in meters. *ρ ij rb* is the DD geometric distance from the satellite to receiver. *<sup>T</sup>ij rb* is the DD slant tropospheric delay of each satellite pair. *I<sup>i</sup> rb*,1 and *I j rb*,1 are the SD ionospheric delays on *L*1. *N<sup>i</sup> rb*,*WL*, *<sup>N</sup><sup>j</sup> rb*,*WL*, *<sup>N</sup><sup>i</sup> rb*,1, and *<sup>N</sup><sup>j</sup> rb*,1 are the SD *WL* and *L*1 ambiguities. *f*<sup>1</sup> and *f*<sup>2</sup> are the frequencies of *L*1 and *L*2, and *λWL* and *λ*<sup>1</sup> are the wavelengths of *WL* and *L*1. *ε ij rb*,*WL*, *ε ij rb*,1, *e ij rb*,1, and *e ij rb*,2 denote all other errors, including the noise of the DD carrier and the pseudorange. *I ij rb*,1 is DD ionospheric pseudo-observation.

Using *WL* ambiguity obtained by Equation (1), the *L*1 ambiguity was then estimated using the following observation equation:

$$\begin{cases} \begin{array}{c} \mathcal{L}^{ij}\_{rb,IF} = \rho^{ij}\_{rb} + T^{ij}\_{rb} + \frac{\mathbb{C}}{f\_1 + f\_2} (\mathcal{N}^{j}\_{rb,1} - \mathcal{N}^{i}\_{rb,1}) + \frac{\mathbb{C}f\_2}{f\_1^2 + f\_2^2} \mathcal{N}^{ij}\_{rb,NL} + \mathcal{e}^{ij}\_{rb,IF} \\\ P^{ij}\_{rb,IF} = \rho^{ij}\_{rb} + T^{ij}\_{rb} + \mathcal{e}^{ij}\_{rb,IF} \end{array} \tag{2}$$

where *IF* denotes the ionosphere-free combination. *Nij rb*,*WL* is DD WL ambiguities. The other symbols are the same as in Equation (1).

The prerequisite for obtaining the DD tropospheric delay is to resolve the ambiguity correctly. Generally, the baseline difference of the reference stations is tens of kilometers or even more than 100 km, and the AR is difficult. Therefore, we used Equation (1) to estimate the SD *WL* ambiguity, and the SD *L*1 ambiguity was then estimated by applying Equation (2). The ionospheric delay is one of the major factors restricting the AR among the reference stations. This effect is handled with a temporal and spatial correlation weighting scheme for the ionospheric pseudo-observations in Equation (1) [13,24]. The DD ionospheric pseudo-observations can be obtained from the broadcast Klobuchar ionospheric model [25] or the IGS global ionosphere models (GIM). We selected the GIM model because it has better accuracy than the broadcast Klobuchar model [26]. On the other hand, the UNB3m model combined with Niell mapping functions (NMF) was used to compute the hydrostatic tropospheric corrections [27,28]. The wet tropospheric delay was estimated as the relative tropospheric zenith delay (RZTD) parameter [29].

The Kalman filter was adopted to estimate the float SD *WL*, *L*1 ambiguity, the SD ionospheric delay, and the RZTD. The DD ambiguity was obtained from the SD ambiguity through the transition matrix, and the LAMBDA method was applied to obtain the integer of the DD *WL* and *L*1 ambiguities [30]. The *L*2 DD ambiguity could then be directly derived from the linear relationship between the *WL* and *L*1 DD ambiguities.

After we obtained the *L*1 and *L*2 DD integer ambiguities, the tropospheric slant delay corresponding to each satellite pair between *r* and *b* was calculated as follows:

$$\hat{T}\_{rb}^{ij} = \frac{f\_1^2}{f\_1^2 - f\_2^2} (L\_{rb,2}^{ij} - \lambda\_1 \tilde{N}\_{rb,1}^{ij}) - \frac{f\_2^2}{f\_1^2 - f\_2^2} (L\_{rb,2}^{ij} - \lambda\_2 \tilde{N}\_{rb,2}^{ij}) - \rho\_{rb}^{ij} \tag{3}$$

where *T*ˆ*ij rb* represents the calculated DD tropospheric delay. *<sup>N</sup>ij rb*,1 and *<sup>N</sup>ij rb*,2 denote fixed DD ambiguities on *L*1 and *L*2 in units of cycles. Observation noise, multipath effects, and other errors are ignored here.

Generally, at least three reference stations around the rover are used to build a DD tropospheric interpolation model, and the reference station closest to the rover is selected as the master reference station, while the other two are auxiliary reference stations. The interpolation model for applying the LIM can be expressed as

$$
\hat{T} = B X\_{\text{ab}} \tag{4}
$$

$$
\hat{T} = \begin{bmatrix}
\hat{T}\_{M,1}^{ij} \\
\hat{T}\_{M,2}^{ij} \\
\vdots \\
\hat{T}\_{M,n-1}^{ij}
\end{bmatrix} B = \begin{bmatrix}
\Delta X\_{M,1} & \Delta Y\_{M,1} \\
\Delta X\_{M,2} & \Delta Y\_{M,2} \\
\vdots & \vdots \\
\Delta X\_{M,n-1} & \Delta Y\_{M,n-1}
\end{bmatrix} \chi\_{ab} = \begin{bmatrix}
a \\
b
\end{bmatrix} \tag{5}
$$

where *i*, *j* denote the reference and rover satellite, subscript *M* is the master reference station, and 1, 2, ... , *n* − 1 represent the auxiliary reference stations. Δ*X* and Δ*Y* are the difference of the horizontal coordinates between the master and auxiliary reference stations. *T*ˆ is the calculated DD tropospheric delay, and *B* is the coefficient matrix. When the DD tropospheric delay is determined by (3), the coefficients *a* and *b* in (5) can be estimated by least-squares adjustment:

$$\mathcal{X}\_{ab} = \left(B^T B\right)^{-1} B^T \mathcal{T} \tag{6}$$

The DD tropospheric delay between the VRS and the master reference station can then be obtained as follows:

$$\begin{bmatrix} T\_{M,V}^{\vec{\text{j}}} \end{bmatrix} = \begin{bmatrix} \ \Delta X\_{M,V} & \Delta Y\_{M,V} \end{bmatrix} \cdot X\_{\mathsf{a}\mathsf{b}} = \begin{bmatrix} \ \Delta X\_{M,V} & \Delta Y\_{M,V} \end{bmatrix} \cdot \begin{bmatrix} B^T B \end{bmatrix}^{-1} B^T \hat{T} \tag{7}$$

where *Tij <sup>M</sup>*,*<sup>V</sup>* is the interpolated DD tropospheric delay between the master reference station and the VRS for the satellite pair *i*, *j*. Δ*XM*,*<sup>V</sup>* and Δ*YM*,*<sup>V</sup>* are the difference in the horizontal coordinates between the VRS and the master reference station.

#### *2.2. The Modified Tropospheric Interpolation Algorithm*

In GNSS data processing, the ionosphere is assumed to be a single-layer model above the earth, and its height is generally 350 km. For a small area, the point where the satellite signal passes through the single layer will form an approximate plane. Therefore, when ionospheric changes are not active, using conventional LIM to interpolate the ionospheric delay will achieve a sufficient accuracy in the NRTK. Compared to the ionospheric delay, the tropospheric temporal and spatial distribution characteristics are significantly different. The tropospheric delay is affected by both horizontal and height factors. It can be seen from Section 2.1 that the traditional interpolation method only uses horizontal coordinates to interpolate tropospheric errors, but ignores the vertical element. This method constrains the interpolation of the tropospheric delay to the height level constituted by reference stations. When the height of the rover is much higher or lower than the selected reference stations, it causes the positioning performance to decrease, especially in the up component. Therefore, we proposed that the tropospheric system error caused by the height difference between the reference stations and rovers should be corrected first. The LIM was then used for the DD tropospheric delay modeling.

Assuming that the reference stations *M*, *A,* and *B* form a triangular unit and assuming that station *M* is selected as the master reference station, *A* and *B* are the two auxiliary reference stations. The coordinates of the VRS are the result of the standard point positioning (SPP) of the rover. Suppose the height of the VRS is much higher or lower than the reference stations. In that case, the tropospheric delay between the master station and the auxiliary reference stations calculated by Equation (3) needs to be corrected in advance before applying the LIM. Take the baseline between *A* and *M* as an example. The calculated DD tropospheric delay is

$$T\_{MA}^{ij}(\Delta h\_{MA}) = T\_A^{ij}(h\_A) - T\_M^{ij}(h\_M) \tag{8}$$

where *Tij MA*(Δ*hMA*) is the calculated DD tropospheric delay of the reference and rover satellite, *<sup>i</sup>* and *<sup>j</sup>* are on the baseline *MA* with the height difference <sup>Δ</sup>*hMA*. *<sup>T</sup>ij <sup>A</sup>*(*hA*), and *Tij <sup>M</sup>*(*hM*) are the SD tropospheric delays of the reference and rover satellites *i*, *j* at stations *A* and *M*, with heights of *hA* and *hM*, respectively.

The SD tropospheric delay correction equations for station *A* and *M* are as follows:

$$
\overline{T}\_A^{ij}(h\_A) = T\_A^{ij}(h\_A) + T\_{\text{Cor}A} \tag{9}
$$

$$
\overline{T}\_M^{ij}(h\_M) = T\_M^{ij}(h\_M) + T\_{CorM} \tag{10}
$$

where *Tij <sup>A</sup>*(*hA*) and *<sup>T</sup>ij <sup>M</sup>*(*hM*) are the corrected SD tropospheric delays for stations *A* and *M*. The specific derivations of the correction terms *TCorA* and *TCorM* are as follows:

$$\begin{split} T\_{\text{CorA}} &= T\_A^{ij}(\mathbf{h}\_V) - T\_A^{ij}(\mathbf{h}\_A) \\ &= ZTD\_A(\mathbf{h}\_V) \times \left( MF\_A^\dagger(\mathbf{h}\_V) - MF\_A^i(\mathbf{h}\_V) \right) - ZTD\_A(\mathbf{h}\_A) \times \left( MF\_A^\dagger(\mathbf{h}\_A) - MF\_A^i(\mathbf{h}\_A) \right) \end{split} \tag{11}$$

$$\begin{split} T\_{\text{CevM}} &= T\_M^{lj}(\mathbf{h}\_V) - T\_M^{lj}(\mathbf{h}\_M) \\ &= ZTD\_M(\mathbf{h}\_V) \times \left( MF\_M^{i}(\mathbf{h}\_V) - MF\_M^{i}(\mathbf{h}\_V) \right) - ZTD\_M(\mathbf{h}\_M) \times \left( MF\_M^{j}(\mathbf{h}\_M) - MF\_M^{i}(\mathbf{h}\_M) \right) \end{split} \tag{12}$$

where *TCorA* and *TCorM* are the correction items, and *hA* and *hV* are the height of station *A* and the VRS, respectively. *ZTDA*(*hA*) and *ZTDA*(*hV*) represent the *ZTD* at station *A* with height *hA* and *hV*, respectively, which can be calculated using the prior model.

*MF<sup>j</sup> <sup>A</sup>*(*hA*), *MF<sup>i</sup> <sup>A</sup>*(*hA*), *MF<sup>j</sup> <sup>A</sup>*(*hV*) and *MF<sup>j</sup> <sup>A</sup>*(*hV*) denote the tropospheric mapping functions corresponding to the satellites *i*, *j* at station *A* with the height *hA* and *hV*, respectively.

Although the *ZTD* calculated by the prior model has model errors, the difference of the *ZTD* obtained by the prior model at the same horizontal position but at different heights can effectively eliminate the tropospheric system error caused by the height difference. Therefore, the *TCorA* and *TCorM* obtained by the prior model still have a high correction accuracy. We substituted Equations (11) and (12) into Equations (9) and (10), respectively, and replaced *Tij <sup>A</sup>*(*hA*) and *<sup>T</sup>ij <sup>M</sup>*(*hM*) in Equation (8) with *<sup>T</sup>ij <sup>A</sup>*(*hA*) and *<sup>T</sup>ij <sup>M</sup>*(*hM*) in Equations (9) and (10). Equation (8) then became

$$
\overline{T}\_{MA}^{ij}(\Delta h\_{VV}) = \overline{T}\_A^{ij}(h\_A) - \overline{T}\_M^{ij}(h\_M) = T\_A^{ij}(h\_V) - T\_M^{ij}(h\_V) \tag{13}
$$

where *<sup>T</sup>ij MA*(Δ*hVV*) is the corrected DD tropospheric delay of the reference and rover satellites *i*, *j* on the baseline *MA*.

The mean error of the *ZTD* calculated by the UNB series model in North America is about 2 cm [31], which represents good accuracy. The NMF and other mapping functions have similar accuracy. Taking the experimental data from North America into account, the tropospheric systematic deviation correction was calculated using the UNB3m model with the NMF in this paper. After correcting the systematic error of the DD tropospheric delay in advance, and then using the LIM for interpolation, a DD tropospheric delay with no height difference deviation on the VRS baseline and the master station could be obtained.

#### **3. Results**

In this chapter, we first describe the two experimental datasets used to verify the modified interpolation algorithm. After that, in the first experiment, we chose different periods and satellites in various triangular units in Dataset One to analyze the tropospheric delay interpolation effect of the traditional and modified methods. In the second experiment, several satellites in the same period in the same triangle unit in Dataset Two were selected to analyze the interpolation effects. Finally, the positioning performance after adopting the traditional and proposed methods was analyzed.

#### *3.1. Experimental Data*

The experimental data were from the US CORS network. We chose 21 stations as reference stations and another ten stations as simulated rovers for experimental analysis. These data were divided into two experimental datasets.

Dataset One consisted of six independent triangular units, and each unit contained a simulated rover. Dataset Two only had one triangle unit, but it included four simulated rover stations. The plane distribution and baseline length of Dataset One and Dataset Two are shown in Figures 1 and 2. The reference station closest to each rover was selected as the master reference station; they are connected by the green dashed line in the figure. Figures 3 and 4 show the station height distribution of the two datasets. The observation data on day of year (DOY) 180 in 2020 were selected. The sampling interval was set to 30 s, and the cut-off angle was set to 15◦.

**Figure 1.** Geographical distribution of the selected stations in Dataset One.

**Figure 2.** Geographical distribution of the selected stations in Dataset Two.

**Figure 3.** The elevation distribution of the stations in Dataset One.

**Figure 4.** The elevation distribution of the stations in Dataset Two.

The coordinates of the reference stations and the simulated rovers were obtained from the US CORS website and were used as known coordinates in the experiments. The Kp index, proposed by Bartels, is an indicator of geomagnetic activity and is updated every three hours [32]. We checked the Kp value of DOY 180 in 2020; the maximum was 2 and the minimum was 0, which meant that the ionosphere was quiet on that day. Therefore, the DD ionospheric delay obtained by the LIM had high accuracy, and it did not affect the

tropospheric delay interpolation. Multipath effects, antenna phase center correction, and other errors were ignored.

The slant DD tropospheric delay between the auxiliary and master stations is calculated using Equation (3). The traditional LIM and modified LIM are then applied to obtain the slant DD tropospheric delay between the rover and master stations. The true value of *ZTD* between the rover and the master reference stations could be obtained using the latest Canadian Spatial Reference System precise point positioning (CSRS-PPP) service, provided by Natural Resources Canada [33]. Since the estimation accuracy could reach 0–2 cm [34,35], the slant DD tropospheric delay of each satellite pair calculated in combination with the mapping function could be used as the reference true value. The whole experiment was conducted in post-processing mode.

#### *3.2. Analysis of Tropospheric Delay Interpolation Results with Dataset One*

The triangular units of Dataset One were distributed in different areas, and the height distribution of the reference stations and rovers also had apparent differences. They were used to verify the effectiveness of the proposed method compared to the traditional LIM. The DD tropospheric delay was obtained using the CSRS-PPP as the true value. The values of the root mean square (RMS) of the tropospheric delay obtained by the two methods were analyzed and compared.

In each triangle unit, we selected different periods, different elevation angle changes, and three continuously observed satellites to analyze the effect of tropospheric interpolation. The results were obtained using both the traditional LIM and the modified LIM, which are denoted here as TLIM and MLIM. A comparison with the true value is shown in Figures 5–10.

**Figure 5.** The DD tropospheric delay obtained by different methods (**a**) and elevation angles (**b**) for satellites G10, G18, and G32 at the rover station P478.

**Figure 6.** The DD tropospheric delay obtained by different methods (**a**) and elevation angles (**b**) for satellites G15, G18, and G21 at the rover station EWPP.

**Figure 7.** The DD tropospheric delay obtained by different methods (**a**) and elevation angles (**b**) for satellites G08, G14, and G20 at the rover station P603.

**Figure 8.** The DD tropospheric delay obtained by different methods (**a**) and elevation angles (**b**) for satellites G01, G14, and G31 at the rover station P225.

**Figure 9.** The DD tropospheric delay obtained by different methods (**a**) and elevation angles (**b**) for satellites G16, G22, and G26 at the rover station WINT.

**Figure 10.** The DD tropospheric delay obtained by different methods (**a**) and elevation angles (**b**) for satellites G07, G11, and G13 at the rover station P121.

We selected three satellites with different elevation angle range changes on each rover for the analysis. The trend of elevation angle changes for these selected satellites was also different. Overall, from Figures 5–10, we can see that as the elevation angles of each satellite increased, the DD tropospheric delay became smaller, and vice versa. These selected satellites can be divided into three types: low-elevation, medium-elevation, and high-elevation angle satellites.

For high elevation angle satellites, we took G10, G18, and G14 in Figures 5–7 as examples. The DD tropospheric delay obtained using TLIM at the VRS had a slight offset from the true value. However, the DD tropospheric delay corrected using MLIM almost coincided with the true value. For G14, which had the largest elevation angle increase, the results of the proposed method were the most consistent with the true value.

Similarly, we took G18, G20, and G14 in Figures 5, 7 and 8 as medium elevation angle satellites for analysis. These three satellites showed a downward trend in the range of medium elevation angles. The offset between the DD tropospheric delay obtained by the TLIM and true value was obvious. As the elevation angles of these satellites decreased, these offsets became larger. The changing trend of the TLIM value and the true value was not the same. When MLIM was adopted, the offsets between the interpolation results of the MLIM and true value were mostly eliminated, and a better correction effect was obtained.

From Figures 7, 9 and 10, the low elevation angle satellites G08, G16, and G11 were analyzed. Their elevation angles were relatively low, and they showed almost no change trend. Therefore, their tropospheric delay change trends according to the TLIM, MLIM, and true values were not noticeable. However, their tropospheric delay values were considerable. This is because the lower the elevation, the longer the satellite signal propagation path becomes, and the delay of signal propagation will increase the range of the troposphere. The most significant offsets between TLIM and the true value occurred with this type of elevation satellite. Although the DD tropospheric delay obtained by the MLIM still had a small difference from the true value, it effectively eliminated this offset compared to the TLIM.

It can be seen that whether the DD tropospheric delay interpolated by the proposed method is from high elevation or low elevation satellites, the offset with the true value is smaller. For all of the satellites, the trend of the MLIM results was consistent with the true value.

The statistical analysis of the results of the TLIM and MLIM for each rover is shown in Figure 11. Overall, the RMS errors of the MLIM are smaller than those of the TLIM. This indicates that the accuracy of the interpolated tropospheric delay improved after applying the proposed method. The RMS errors of the MLIM on each continuously observed satellite showed a significant improvement. The maximum increases were achieved on the G32, G29, G08, G11, G16, and G11 satellites at stations P478, EWPP, P603, P225, WINT, and P121, respectively, with an improvement of 90.1%, 92.5%, 90.0%, 93.2%, 90.6%, and 78.5%.

**Figure 11.** RMS errors of the DD tropospheric delay using modified and traditional methods for all non-reference satellites at all rover stations. The percentages on each bar are the improvements of the MILM with respect to the TLIM.

To analyze all tropospheric interpolation accuracies for all of the satellites with different elevation angles during the selected period for each rover, Figure 12 shows the TLIM and MLIM modeling error curves for all of the satellites with different elevation angles. Overall, from Figure 12, we can see the accuracy of the tropospheric modeling using TLIM increased with increasing satellite elevation. However, the accuracy of the MLIM was always better than that of the TLIM at any elevation interval. For low elevation satellites, the error of the TLIM was obvious. For example, the error forP603 reached 0.2 m, but the error from the MLIM was still minimal. The mean error values of the TLIM were 0.017 m, −0.027 m, 0.068 m, 0.029 m, 0.045 m, and 0.010 m, and the corresponding standard deviations (STD) were 0.014 m, 0.023 m, 0.048 m, 0.023 m, 0.031 m, and 0.006 m, respectively. The accuracy of the MLIM was significantly improved, especially for low elevation satellites. The mean errors of each rover were 0.002 m, 1.65 × <sup>10</sup>−<sup>5</sup> m, 0.007 m, 0.005 m, 0.005 m, and 0.003 m, and the corresponding STDs were 0.002 m, 0.004 m, 0.006 m, 0.004 m, 0.004 m, and 0.002 m, respectively. The discontinuous part in Figure 12 indicates that there were no visible or ambiguity fixed satellites within that elevation angle range.

**Figure 12.** The DD tropospheric delay interpolation errors obtained by the TLIM and MLIM (left and right) with different satellite elevation angles at different rover stations: (**a**) P478, (**b**) EWPP, (**c**) P603, (**d**) P225, (**e**) WINT, and (**f**) P121.

The mean RMS errors of the DD tropospheric delay obtained by the two methods for all satellites are listed in Table 1. The statistical analysis showed that the RMS of the TLIM reached tens of centimeters, while the RMS of the interpolation using the proposed method was only several millimeters. The improvement reached 87.5%, 87.9%, 89.5%, 85.4%, 87.9%, and 66.7%, respectively. This indicated that the modified algorithm performed better than the traditional algorithm, especially for satellites with lower elevation angles.

**Table 1.** The mean RMS and improvement percentage of the interpolated DD tropospheric delay using the TLIM and MLIM at each rover station in Dataset One.


*3.3. Analysis of the Tropospheric Delay Interpolation Results with Dataset Two*

Dataset Two only had one triangle unit, but it contained four rovers with different elevations. We also compared and analyzed the performance of the TLIM and MLIM at these stations.

We selected three identical satellites, G10, G18, and G27, at these four rovers to analyze the effect tropospheric interpolation during the same period. The results are illustrated in Figure 13.

**Figure 13.** Elevation angles (upper row) for the reference satellites G10, G18, and G27. The DD tropospheric delay obtained using the TLIM and the PPP method (three lower rows) at different rover stations: (**a**) P532, (**b**) P539, (**c**) P602, and (**d**) MASW.

As shown in Figure 13, the interpolated DD tropospheric delay by the MLIM was highly consistent with the true value for G10 at the highest elevation. For G18, which had an elevation angle on a downward trend, although there were still some deviations compared to the true value, the MLIM results were also closer to the true value. Analyzing the satellite G27, although its elevation angle was increasing, it had the lowest elevation angle among the three satellites. Its TLIM results showed that the effect of directly interpolating the tropospheric delay was not satisfactory. However, after adopting the MLIM, the value of the DD tropospheric delay interpolation was consistent with the true value. Therefore, the modified interpolation method had a better correction effect, even for low-elevation satellites.

The RMS of the DD tropospheric delays obtained using the TLIM and MLIM for each rover during the selected period are shown in Figure 14. The maximum improvement values of the MLIM at the four rovers compared to the TLIM were 93.0%, 83.4%, 78.3%, and 79.3%, respectively.

Further, troposphere modeling errors for all of the satellites over the selected observation period are shown in Figure 15. The mean errors of the TLIM were 0.011 m, 0.027 m, 0.029 m, and 0.021 m, and the corresponding STDs were 0.010 m, 0.027 m, 0.029 m, and 0.020 m, respectively. The accuracy of the MLIM was dramatically improved, with mean errors of −0.002 m, −0.006 m, 0.009 m, and 0.006 m, and the corresponding STDs were 0.007 m, 0.007 m, 0.008 m, and 0.006 m, respectively.

The mean RMS of the DD tropospheric delay using the traditional and the modified methods for all of satellites in Dataset Two is shown in Table 2. It can be seen that the mean RMS error of the TLIM of each station, P532, P539, P602, and MASW, were 0.012 m, 0.033 m, 0.035 m, and 0.026 m, respectively. However, the mean RMS error of the MLIM were 0.005 m, 0.008 m, 0.010 m, and 0.007 m. These values represent improvements of 58.3%, 75.6%, 71.4%, and 73.1%, respectively, which are consistent with the requirements of the accuracy of DD tropospheric delay interpolation that need to be met in order to generate VRS observations in the NRTK.

**Figure 14.** RMS errors of the DD tropospheric delay using modified and traditional methods for all non-reference satellites at all rover stations. The percentages on each bar are the improvements of the MILM with respect to the TLIM.

**Figure 15.** DD tropospheric delay interpolation errors obtained by the TLIM and MLIM (left and right) with different satellite elevation angles at different rover stations (from top to bottom): P532, P539, P602, and MASW.


**Table 2.** The mean RMS errors and improvement percentage of the interpolated DD tropospheric delay using the TLIM and MLIM at each rover station in Dataset Two.

#### *3.4. Comparison of the Positioning Results of Dataset One and Dataset Two*

According to the proposed and traditional tropospheric interpolation delay results, the positioning performance of these rovers using the VRS was investigated in this section. The DD tropospheric delay was obtained through both the TLIM and MLIM. We then used these interpolated atmospheric delays to generate virtual carrier phase and pseudorange observations, which were applied to the final rover positioning. The RMS errors of the positioning results were analyzed in this paper. We expected that the proposed method would show a better positioning performance, especially in the up component and under significant height differences.

#### 3.4.1. Dataset One

The positioning error of the TLIM and MLIM is shown in Figure 16, using station P225 as an example. It can be seen that the positioning error of the MLIM showed a significant improvement in the up component. However, there was little difference between the TLIM and MLIM in the horizontal direction. Only the results of the fixed solutions were counted here, and the MLIM had more fixed solutions than the TLIM. The positioning results of each rover in Dataset One are given in Table 3. In terms of the RMS, the accuracy in the three directions for all of the rovers was improved after using the MLIM. There were slight improvements in the horizontal direction, but remarkable improvements were achieved in the up component. All the RMS errors of the up component obtained using the MLIM were within 5 cm. The most obvious improvement occurred at P603, in which the accuracy of the RMS, changing from 24.95 cm to 4.24 cm, improved by 83%. It can be seen from Figure 16 that the positioning result became more stable after the proposed method was applied.

**Figure 16.** The TLIM (**a**) and MLIM (**b**) positioning errors with respect to the epoch during the selected data period at rover P225.


**Table 3.** Positioning RMS statistical results of each rover station in Dataset One.

#### 3.4.2. Dataset Two

Similarly, taking MASW as an example, it can be seen from Figure 17 that the accuracy improvement in the up direction was obvious. Table 4 shows the statistical positioning results of all of the rovers in Dataset Two. Analyzing the RMS results of the four rover stations, the accuracy of the TLIM was within 2 cm in the east and north directions. The accuracy of the up component was about 17 cm. After applying the MLIM, the positioning accuracy of all rovers was within 1 cm in the east and north directions. The accuracy in the up direction was within 5 cm, reaching the general accuracy level of the NRTK positioning. Figure 17 and Table 4 showed that in the case of the height distribution of Dataset Two, the positioning results using the modified interpolation methods were reliable and stable. %clearpage

**Figure 17.** TLIM (**a**) and MLIM (**b**) positioning errors with respect to the epoch during the selected data period at rover MASW.


**Table 4.** Positioning RMS statistical results of each rover station in Dataset Two.

The results of these experiments demonstrate that the accuracy of the traditional LIM DD tropospheric delay at the VRS deteriorates when a relatively significant height difference exists between the rover and the reference stations. The results also confirm that as the satellite elevation angle drops, the delay of the satellite signal passing through the troposphere increases, and the effect of interpolating the tropospheric delay will worsen.

It can be seen from these results that the more significant the height difference between the rover and the reference station, the worse the tropospheric accuracy obtained by using TLIM is, and the better the correction effect by using MLIM is. The average RMS error of the conventional method was 0.035 m. The MLIM's average RMS was only 0.006 m.

According to the results when using VRS observations for rover positioning, the modified method showed a slight improvement in the horizontal direction compared to the traditional method. However, the accuracy of the proposed method was significantly improved in the vertical direction. Therefore, when building or selecting reference station sites, one should thoroughly consider whether the height distribution of the reference station is approximately consistent. In addition, the height difference between the potential users and the reference stations should be as small as possible. Otherwise, a method similar to the MLIM should be used to solve the impact of the height difference in the NRTK.

It should be noted that this experiment only verified that a considerable height difference between the reference station and the rover would cause systematic errors when modeling the DD tropospheric delay in the NRTK. The heights of the reference stations and the rovers were quite different: Dataset One was from −29 m to about 1900 m, and Dataset Two was from 34 m to 713 m. The maximum height difference between the rover and reference stations of each triangle unit is from 245 m to 1158 m in Dataset One, and the maximum height difference is from 499 m to 679 m in Dataset Two. In addition, when the LIM is used to interpolate the atmospheric delay, the spatial distance and height difference between the rover and each reference station should be treated as the weighting factor for interpolating the DD tropospheric delay in each triangle unit. This perspective should be studied through further experiments in the future.

#### **4. Conclusions**

This study confirmed that a height difference between reference stations and rovers leads to a decrease of interpolation DD tropospheric delay accuracy. Therefore, based on the tropospheric delay correction model and traditional LIM, we proposed the MLIM correction method for the triangular VRS mode. This method corrects the systematic error of the tropospheric delay caused by height differences in advance. It then applies the LIM to interpolate the corrected value to obtain the DD tropospheric delay between the rover and the master stations. We compared the accuracy of TLIM and MLIM by modeling the DD tropospheric delay by using two datasets of triangular units in VRS mode. Finally, the positioning performance of the two methods was evaluated.

In the first experiment, we selected six independent triangular units in different areas. The accuracy of the DD tropospheric delay modeling, which used traditional LIM, was

inferior. This phenomenon is evident for satellites with low elevation angles. Compared to traditional LIM, the accuracy of interpolated DD tropospheric delay was significantly improved after using the modified method, especially for the satellites with lower elevation angles. The second experiment considered only one triangular unit, but it included rover stations with different elevation distributions. From the perspective of the tropospheric interpolation results, as the height of each rover increased, the accuracy of the tropospheric interpolation using the conventional method deteriorated. After adopting the MLIM, the RMS error was significantly improved.

In terms of positioning performance, the up component accuracy obtained by the TLIM showed a noticeable systematic deviation, and the worst results in the two datasets could still reach the decimeter level. However, this situation improved significantly after using the MLIM, and the positioning results of all of stations in the up direction reached just a few centimeters. In the horizontal direction, although the accuracy of the proposed algorithm was similar to that of the traditional algorithm, there was still a slight improvement. In addition, the positioning result of the proposed method is more stable than that of the conventional method.

In summary, the proposed method significantly improved the tropospheric delay interpolation and the positioning performance in the NRTK triangulation network. It should be noted that this paper only analyzed the influence of the DD tropospheric delay interpolation in the NRTK where the height distribution of the reference stations and the rover stations were quite different, and further research should be devoted to investigating how the height difference of the reference stations and the spatial distance between the rover and reference stations can affect the DD tropospheric interpolation accuracy in the NRTK.

**Author Contributions:** Y.P. provided the initial idea and wrote the manuscript; Y.P. and M.S. designed and performed the research; M.S. helped with the writing and revising, and Y.Y. and M.S. partially financed the research. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Key Research Program (No. 2016YFB0501900) and China Natural Science Funds (No. 41604033 and 41874093). The third author is supported by the Wang Kuancheng Pioneer Talents Project Lu Jiaxi International Team Program.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Publicly available datasets were analyzed in this study. The GPS observation data from the NOAA CORS network are available at https://www.ngs.noaa.gov/ UFCORS/ (accessed on 26 April 2021). The GPS broadcast ephemeris data are available at https: //cddis.nasa.gov/archive/gnss/data/daily/ (accessed on 26 April 2021). The GIM products from IGS can be obtained at https://cddis.nasa.gov/archive/gnss/products/ionex/ (accessed on 26 April 2021).

**Acknowledgments:** The authors gratefully acknowledged the Canadian Spatial Reference System Precise Point Positioning (CSRS-PPP) service, the NOAA CORS network for providing the GNSS data, and the IGS for providing the GIM products.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

