1. INTRODUCTION
Network Real-Time Kinematic (RTK) positioning is a carrier phase-based technique for centimetre-level positioning accuracy of users that has been developed to extend the coverage of a single baseline RTK. Global Positioning System (GPS) errors at the user's location can be interpolated using multiple corrections from reference stations constituting a network. Those multiple corrections can be combined through various interpolation methods. The Linear Combination Model (LCM) was proposed to model the GPS spatial errors, such as satellite orbit error, ionospheric and tropospheric delay, according to the horizontal directions based on the known coordinates of reference stations and approximate user coordinates (Han and Rizos, Reference Han and Rizos1996). The Distance-based Linear Interpolation Method (DIM) (Gao et al., Reference Gao, Li and McLellan1997) calculates a weight for the correction of each reference station proportional to the inverse of the distance between a user and a reference station. The Linear Interpolation Method (LIM) has been proposed to model the ionospheric delay (Wanninger, Reference Wanninger1995). In addition, the Low-order Surface Model (LSM) is based on the partial derivative principles with respect to the coordinates of the horizontal and up directions (Schaer et al., Reference Schaer, Beutler, Rothacher, Brockmann, Wiget and Wild1999; Varner, Reference Varner2000). When interpolating corrections using the LSM, the Least Square Estimation (LSE) is used for estimating the modelling parameters. Meanwhile, the corrections of Network RTK can be interpolated by the Least Squares Collocation (LSC) based on the covariance matrix of the measurement and the cross-covariance matrix between the measurement vector and interpolated vector (Raquet, Reference Raquet1997).
Users’ positioning accuracy depends on how accurately the correction from Network RTK is modelled to the actual GPS error sources at the user location through the various interpolation methods. Among the various GPS error sources, the tropospheric delay is known to be significantly affected by height. Moreover, the heights of the reference stations can differ by hundreds of metres according to the terrain. The tropospheric delay difference due to height difference can be compensated through the tropospheric model, however, they cannot be eliminated completely since the accuracies of the tropospheric delay models are limited. Collins and Langley (Reference Collins and Langley1999) have evaluated the comparison in the performance of the UNB (University of New Brunswick) 3 tropospheric delay model to the real tropospheric delay values determined using the North American radiosonde data for ten years. These values were collected by the Forecast Systems Laboratory (FSL) of the United States National Oceanic and Atmospheric Administration (NOAA) in 1999. As a result, the mean value of the residual errors was calculated to be nearly zero. However, the standard deviation was approximately 5 cm. Moreover, it stated that there existed extreme values of residual errors, whose value is larger than ±20 cm, for 76 out of 1,011,651 cases. These large residual errors can introduce incorrect integer ambiguities, therefore, it is appropriate to select the LSM containing height-related terms to model the tropospheric delay in this case.
On the other hand, the additional consideration of height requires an increased number of reference stations since the modelling parameters are estimated by the LSE. This may not be practical since the reference stations have already been installed in a fixed location. In this paper, Minimum-Norm Estimation (MNE) is considered to estimate the unknown parameters for the height-related LSM since it provides valid solutions even if the number of unknowns is larger than the number of observations. In the meanwhile, conventional height-related interpolation models of the LSM are based on the partial derivative principles with respect to the coordinates in the horizontal and up directions. In order for the corrections of Network RTK to reflect the real tropospheric delay more closely, the relationship between the tropospheric delay and height is considered as an exponential function. Subsequently, it is linearized and further applied to the height-related interpolation models of the LSM. The resultant interpolation method has the same form as the conventional LSM that includes the first and second orders of the height difference, however it has additional nonlinear constraints on the modelling parameters that need to be estimated through the nonlinear programming.
In this research, the performance of the various height-related correction interpolation models is analysed. These models include the conventional height-related interpolation models of the LSM and LIM, those with the MNE, and those with the characteristics of the tropospheric delay over height. In the following sections, the conventional height-related interpolation models are introduced. Then, the proposed methods, which reflect the relationship between tropospheric zenith delay and height, are formulated. Two possible reference station networks in mountainous areas of North Carolina and Vermont, USA are selected for performance evaluation. The Root Mean Square (RMS) errors of the observation residuals and the position-domain errors of the users are calculated to assess the performance. As a result, the conventional methods and those with the MNE produce similar levels of performance. Furthermore, the proposed methods, which consider the relationship between the tropospheric delay and height, provide improved levels of performance for selected networks.
2. METHODOLOGY
In this section, we first formulate equations for corrections of Network RTK and conventional interpolation methods. Then we propose the interpolation method which considers the characteristics of tropospheric delay over height.
2.1. Corrections of Network RTK
Network RTK can be distinguished by the methods of Virtual Reference Station (VRS), Master-Auxiliary Concept (MAC), and Flächen Korrektur Parameter (FKP) according to how it generates corrections. Since the VRS and FKP can be realized from the MAC approach (Takac, Reference Takac2008) as seen from Figure 1, this study focuses on the MAC Network RTK and briefly formulates the corrections for L1, L2, dispersive and non-dispersive portions.
Figure 1. General process for Network RTK solutions (Takac, Reference Takac2008) and a step for applying correction interpolation method.
In the MAC Network RTK, the reference stations are distinguished as one master station and other auxiliary stations (Euler et al., Reference Euler, Keenan and Zebhause2001). All of the reference stations generate Carrier Phase Corrections (CPC) and the MAC corrections are generated by differencing CPCs between master and auxiliary stations. The multiple MAC corrections are interpolated using appropriate weighting factors and they are combined with the CPC from the master station to model the GPS errors at the approximate location of the user. The equations for the L1 CPC and MAC correction for satellite i can be calculated by subtracting estimated distance (
$\hat d$
), satellite clock bias (
$\hat b$
), integer ambiguities (
$\hat N$
) from carrier phase (Park and Kee, Reference Park and Kee2010) as
$$\eqalign{\delta \phi _{c,f,mas}^{i} &\equiv \phi _{\,f,mas}^{i} - \hat d_{mas}^{i} + \hat b_{mas}^{i} - \hat N_{\,f,mas}^{i} \lambda _f^{i} \cr &= {B_{mas}} - I_{\,f,mas}^{i} + T_{mas}^{i} + \delta R_{mas}^{i} + \delta N_{\,f,mas}^{i} {\lambda _f} \cr _{aux} {\Delta _{mas}}\delta \phi _{c,f}^{i} & \equiv{ _{aux}}{\Delta _{mas}} \left( {B - I_f^{i} + {T^{i}} + \delta {R^{i}} + \delta N_f^{ref} {\lambda _f}} \right)}$$
The detailed derivation of the MAC correction including ambiguity levelling can be found in Park (Reference Park2008). The noise error of the carrier phase is omitted for simplicity. The subscript f represents the transmission frequency and Δ indicates the single difference operator between stations. The subscripts ‘mas’, ‘aux’ and the superscript ‘ref’ represent the master, auxiliary stations and the reference satellite, respectively. The symbols of δϕ c and Δδϕ c indicate the CPC and MAC correction, respectively. In addition, B, I, T, δR and λ are the receiver clock bias, ionospheric delay, tropospheric delay, ephemeris error and wavelength. In this paper, the symbol δ means the deviation from the true value and the hat (̂) indicates the estimated value. The equation for the MAC correction contains the single-differenced estimation error of integer ambiguity of the master station. If there are three auxiliary stations, the final correction at the user location can be calculated as
$$\eqalign{\delta \hat \phi _{c,f,user}^i &= \delta \phi _{c,f,mas}^i + \sum\limits_{n = 1}^3 {_{aux\;n} {\Delta _{mas}}\delta \phi _{c,f}^i \cdot {w_n}} \cr = & \, - \hat I_{\,f,user}^i + \hat T_{user}^i + \delta \hat R_{user}^i + {B_{mas}} + \delta N_{\,f,mas}^i {\lambda _f} + {\Theta _{\,f,w}} \cr where,\ & {\Theta _{\,f,w}} \equiv \sum\limits_{n = 1}^3 {_{aux\;n} {\Delta _{mas}}\left( {\delta N_{c,f}^{ref} + B} \right) \cdot {w_n}}} $$
where the weighting factors of the auxiliary stations w can be calculated based on the configuration of the network (Dai et al., Reference Dai, Han, Wang and Rizos2003). If the user is established at a known position, the CPC of the user can be calculated as
$$\eqalign{\delta \phi _{c,f,user}^i &= \phi _{\,f,user}^i - \hat d_{user}^i + \hat b_{user}^i - \hat N_{\,f,user}^i \lambda _f^i \cr &= {B_{user}} - I_{\,f,user}^i + T_{user}^i + \delta R_{user}^i + \delta N_{\,f,user}^i {\lambda _f}.} $$
The residual error of user, r can be calculated by subtracting Equation (2) from Equation (3) as
$$\eqalign{r_{\,f,user}^i &\equiv \delta \phi _{c,f,user}^i - \delta \hat \phi _{c,f,user}^i \cr &= \left\{ { - \left( {I_{\,f,user}^i - \hat I_{\,f,user}^i} \right) + \left( {T_{user}^i - \hat T_{user}^i} \right) + \left( {\delta R_{user}^i - \delta \hat R_{user}^i} \right)} \right\} \cr &\quad+{ _{user}}{\Delta _{mas}}\delta N_f^i \cdot {\lambda _f} + \left( {_{user} {\Delta _{mas}}B + {\Theta _{\,f,w}}} \right)} $$
The single difference of Equation (4) between satellites i and j can be calculated as
$$\eqalign{^i {\nabla ^j}{r_{\,f,user}} &= {^i}{\nabla ^j}\left\{ { - \left( {{I_{\,f,user}} - {{\hat I}_{\,f,user}}} \right) + \left( {{T_{user}} - {{\hat T}_{user}}} \right) + \left( {\delta {R_{user}} - \delta {{\hat R}_{user}}} \right)} \right\}\cr& \quad+{ _{user}}{\Delta _{mas}}^i {\nabla ^j}\delta {N_f} \cdot {\lambda _f}}$$
where ∇ is the single difference operator between satellites. The last terms of user Δ mas B and Θ f,w in Equation (4) are common for all visible satellites and completely eliminated in Equation (5). Since the residual ambiguities at the master station maintain the integer nature, they are combined with the residual integer ambiguities of the user. Finally, the terms in the first brace of Equation (5) are the residuals of the GPS errors, which affect resolving integer ambiguities and position errors of the user.
The L1 and L2 CPC can be linearly combined to generate dispersive and non-dispersive components. The dispersive part does not contain the tropospheric delay and orbit errors and the non-dispersive part eliminates ionospheric delay. These two components for a satellite i can be expressed as
$$\eqalign{\delta \phi _{c,disp,mas}^i &= \displaystyle{1 \over {\gamma - 1}} \cdot \delta \phi _{c,L1,mas}^i + \displaystyle{{ - 1} \over {\gamma - 1}} \cdot \delta \phi _{c,L2,mas}^i \cr &= I_{L1,mas}^i + \Omega _{disp}^i \cr & where,\quad \Omega _{disp}^i \equiv \displaystyle{1 \over {\gamma - 1}} \cdot \delta N_{L1,mas}^i \cdot {\lambda _{L1}} + \displaystyle{{ - 1} \over {\gamma - 1}} \cdot \delta N_{L2,mas}^i \cdot {\lambda _{L2}} \cr \delta \phi _{c,nondisp,mas}^i &= \displaystyle{\gamma \over {\gamma - 1}} \cdot \delta \phi _{c,L1,mas}^i + \displaystyle{{ - 1} \over {\gamma - 1}} \cdot \delta \phi _{c,L2,mas}^i \cr &= T_{mas}^i + \delta R_{mas}^i + {B_{mas}} + \Omega _{nondisp}^i \cr & where,\quad \Omega _{nondisp}^i \equiv \displaystyle{\gamma \over {\gamma - 1}} \cdot \delta N_{L1,mas}^i \cdot {\lambda _{L1}} + \displaystyle{{ - 1} \over {\gamma - 1}} \cdot \delta N_{L2,mas}^i \cdot {\lambda _{L2}}} $$
where γ is the square of the ratio of the L1 and L2 wavelengths. The estimated dispersive and non-dispersive correction at the user location can be calculated as
$$\eqalign{\delta \hat \phi _{c,disp,user}^i &= \delta \phi _{c,disp,mas}^i + \sum\limits_{n = 1}^3 {_{aux\;n} {\Delta _{mas}}\delta \phi _{c,disp}^i \cdot {u_n}} \cr &= \tilde I_{L1,user}^i + \Omega _{disp}^i + {\Theta _{disp,u}} \cr \delta \hat \phi _{c,nondisp,user}^i &= \delta \phi _{c,nondisp,mas}^i + \sum\limits_{n = 1}^3 {_{aux\;n} {\Delta _{mas}}\delta \phi _{c,nondisp}^i \cdot {v_n}} \cr &= \tilde T_{user}^i + \delta \tilde R_{user}^i + {B_{mas}} + \Omega _{nondisp}^i + {\Theta _{nondisp,v}}} $$
where u and v represent the different weighting factors of the dispersive and non-dispersive observations since the independent interpolation methods are applied in general (Dai et al., Reference Dai, Han, Wang and Rizos2003). Therefore, the GPS errors in the dispersive and non-dispersive corrections are distinguished by tilde (~) from those in the L1 corrections. The values of Θ disp and Θ nondisp are the dispersive and non-dispersive combinations of the L1 and L2 Θ f defined in Equation (2). Using the estimated dispersive and non-dispersive correction, the estimated L1 correction at the user location can be reformed as
$$\eqalign{\delta \hat \phi _{c,L1,user}^i &= - \delta \hat \phi _{c,disp,user}^i + \delta \hat \phi _{c,nondisp,user}^i \cr &= - \tilde I_{L1,user}^i + \tilde T_{user}^i + \delta \tilde R_{user}^i + {B_{mas}} + \delta N_{L1,mas}^i \cdot {\lambda _{L1}} - {\Theta _{disp,u}} + {\Theta _{nondisp,v}}.} $$
Since the values of Θ disp and Θ nondisp are common for all visible satellites, they are eliminated when forming single differences between satellites. Finally, the single-differenced residual of the user calculated from Equation (8) has the same form as Equation (5) where the values in the first brace with hat are replaced by those with tildes. The advantage of using the dispersive and non-dispersive corrections is that the weighting factors can be determined separately according to the characteristics of the ionospheric and tropospheric delays which are the major error components of the dispersive and non-dispersive corrections.
2.2. Conventional LSM and Minimum-norm Estimation
The LSM based on the partial derivative principle was first proposed by Loomis et al. (Reference Loomis, Sheynblatt and Mueller1991) for pseudoranges. Later, Varner (Reference Varner2000) derived the multi-dimensional LSM using the partial derivative principle. Based on Varner's approach, the GPS errors, which include ionospheric delay, tropospheric delay and ephemeris error, are expressed as a function, gi
about a position vector,
$\bar p = \left( {x,y,h} \right)$
for a satellite i. The specified reference point is set at
${\bar p_0} = \left( {{x_0},{y_0},{h_0}} \right)$
, and
${g^i}(\bar p)$
represents the estimated GPS errors at the reference point. If one takes the Taylor series expansion about
${\bar p_0}$
,
${g^i}(\bar p)$
can be expressed as the following:
$$\eqalign{{g^i}(\bar p) &= {g^i}({{\bar p}_0}) + \displaystyle{{\partial {g^i}} \over {\partial x}} \cdot \left( {x - {x_0}} \right) + \displaystyle{{\partial {g^i}} \over {\partial y}} \cdot \left( {y - {y_0}} \right) + \displaystyle{{\partial {g^i}} \over {\partial h}} \cdot \left( {h - {h_0}} \right) + \displaystyle{1 \over 2}\displaystyle{{{\partial ^2}{g^i}} \over {\partial {h^2}}} \cdot {\left( {h - {h_0}} \right)^2} \cr &= {a^i} + {b^i} \cdot \Delta x + {c^i} \cdot \Delta y + {d^i} \cdot \Delta h + {e^i} \cdot \Delta {h^2} \cr where \ & \Delta x \equiv x - {x_0},\;\Delta y \equiv y - {y_0},\;\Delta h \equiv h - {h_0} \cr {a^i} &\equiv {g^i}({{\bar p}_0}),\;{b^i} \equiv \displaystyle{{\partial {g^i}} \over {\partial x}},\;{c^i} \equiv \displaystyle{{\partial {g^i}} \over {\partial y}},\;{d^i} \equiv \displaystyle{{\partial {g^i}} \over {\partial h}},\;{e^i} \equiv \displaystyle{1 \over 2} \cdot \displaystyle{{{\partial ^2}{g^i}} \over {\partial {h^2}}}}$$
where the cross correlation terms and the second order of horizontal directions are assumed to be ignored (Varner, Reference Varner2000). Equation (9) can be modified according to parameters and orders of interest as
$$\eqalign{{g^i}(\bar p) &= {a^i} + {b^i} \cdot \Delta x + {c^i} \cdot \Delta y + {d^i} \cdot \Delta h \cr {g^i}(\bar p) & = {b^i} \cdot \Delta x + {c^i} \cdot \Delta y + {d^i} \cdot \Delta h \cr {g^i}(\bar p) &= {b^i} \cdot \Delta x + {c^i} \cdot \Delta y + {d^i} \cdot \Delta h + {e^i} \cdot \Delta {h^2}.} $$
The parameters ai, bi, ci, di , and ei are unknowns to be estimated using MAC corrections by the LSE as
$$\eqalign{& \left[ {\matrix{ {{{\hat a}_i}} \cr \vdots \cr {{{\hat e}_i}} \cr}} \right] = {\left( {{W^T}W} \right)^{ - 1}}{W^T}\left[ {\matrix{ {_{aux\;1} {\Delta _{mas}}\delta \phi _{c,f}^i} \cr \vdots \cr {_{aux\;n} {\Delta _{mas}}\delta \phi _{c,f}^i} \cr}} \right] \cr & where,\;W \equiv \left[ {\matrix{ 1 \hfill & {_{aux\;1} {\Delta _{mas}}x} \hfill & {_{aux\;1} {\Delta _{mas}}y} \hfill & {_{aux\;1} {\Delta _{mas}}h} \hfill & {_{aux\;1} {\Delta _{mas}}{h^2}} \hfill \cr \vdots \hfill & \vdots \hfill & \vdots \hfill & \vdots \hfill & \vdots \hfill \cr 1 \hfill & {_{aux\;n} {\Delta _{mas}}x} \hfill & {_{aux\;n} {\Delta _{mas}}y} \hfill & {_{aux\;n} {\Delta _{mas}}h} \hfill & {_{aux\;n} {\Delta _{mas}}{h^2}} \hfill \cr}} \right]} $$
where the matrix W represents the configuration of the network. The weighting factors, which have been introduced in the previous section, can also be determined from these parameters (Dai et al., Reference Dai, Han, Wang and Rizos2003). In order to estimate parameters by the LSE, the number of stations should be greater than the number of unknowns by at least one since the single-differenced corrections between stations are used in the MAC Network RTK. Therefore, considering the additional parameter may not be practical because it requires users to receive more correction information. On the other hand, the MNE, or pseudoinverse, gives a solution that yields the minimum Euclidean norm to a system of linear equations that has multiple solutions (Strang and Borre, Reference Strang and Borre1997). Therefore, this method can be attempted when the number of unknowns is greater than the number of observations such as when estimating unknown parameters of height-related interpolation methods.
2.3. Modified LSM
The height-related interpolation methods introduced in the previous section are derived from the mathematical operation of the partial derivative principle. That is, the conventional interpolation methods do not reflect any actual behaviour of the tropospheric delay. This realization brings an idea to consider the actual characteristic of the variation of the tropospheric zenith delay over height. Bean and Dutton (Reference Bean and Dutton1966) suggested that the refractive index can be represented by an exponential function of height. Since the tropospheric zenith delay is derived from the refractive index, the tropospheric zenith delay can also be represented by an exponential function of height. In order to verify this assumption, the sum of the tropospheric hydrostatic and wet zenith delays, which are calculated by UNB (Collins, Reference Collins1999), Black (Black and Eisner, Reference Black and Eisner1984) and Saastamoinen (Saastamoinen, Reference Saastamoinen1972) tropospheric delay models are plotted in Figure 2. In addition, an exponential function of height with appropriate constants is also shown for comparison. As seen from Figure 2, the tendencies of the tropospheric zenith delays that are calculated by the tropospheric delay models and the exponential function are similar.
Figure 2. Tropospheric zenith delays over height for various tropospheric models.
Therefore, the tropospheric zenith delay is assumed to be an exponential function of height (Yin et al., Reference Yin, Huang and Xiong2008; Song et al., Reference Song, Kee, Park, Park and Seo2014) as
$${T_z}(h) = {T_0} \cdot {e^{ - kh}}$$
where T z is the tropospheric zenith delay, T 0 represents a value of the tropospheric zenith delay at sea level, h is the height, and k is an appropriate coefficient. From now on, the superscript i for satellite index is dropped for simplicity. Note that the natural constant, e in Equation (12) is different from ei in Equations (9) and (10). Equation (12) can be expressed as the height difference between the two heights of the reference stations, Δh = h 2 − h 1 as
$${T_z}({h_2}) = T({h_1}) \cdot {e^{ - k\Delta h}}$$
The slant delay can be calculated by multiplying the tropospheric zenith delay by the mapping function. The single-differenced tropospheric slant delay between the two reference stations, 1 and 2 can be expressed as
$$\eqalign{{T_s}({h_2}) - {T_s}({h_1}) & = {T_z}({h_1}) \cdot {e^{ - k\Delta h}} \cdot m({h_2}) - {T_z}({h_1}) \cdot m({h_1}) \cr & = {T_z}({h_1}) \cdot \left( {{e^{ - k\Delta h}} \cdot m({h_2}) - m({h_1})} \right)} $$
The exponential function is linearized through a Taylor series expansion considering up to the second order of the height difference as
$$\eqalign{\Delta {T_s} &\equiv {T_s}({h_2}) - {T_s}({h_1}) \approx {T_z}({h_1}) \cdot \left( {m({h_2}) - m({h_1})} \right) \cr & \quad + {T_z}({h_1}) \cdot m({h_2}) \cdot \left\{ { - k\Delta h + \displaystyle{1 \over 2}{k^2}\Delta {h^2}} \right\}} $$
In order to neglect the higher order terms in Equation (15), the maximum allowable value of height difference should be calculated. Since the non-dispersive MAC correction is used to model the tropospheric delay at user location through the use of Equation (15), we find the maximum value of the height difference which makes the absolute value of the third order of the height difference term,
$\displaystyle{1 \over {3!}}{k^3}\Delta {h^3}$
smaller than the standard deviation of the single-differenced non-dispersive correction,
${\sigma _{\Delta \delta {\phi _{nondisp}}}}$
. This condition can be expressed as
$$\displaystyle{1 \over {3!}}{T_z}({h_1}) \cdot m({h_2}).{\left( {k\Delta h} \right)^3} \le {\sigma _{\Delta \delta {\phi _{nondisp}}}}$$
where the value of
${\sigma _{\Delta \delta {\phi _{nondisp}}}}$
is approximately 8·4 mm from Equation (6) when the standard deviations of both the L1 and L2 carrier phase are set to 2 mm (Dach et al., Reference Dach, Hugentobler, Fridez and Meindl2008). The coefficient k in Equation (16) is set to approximately 0·1309 km−1 from Figure 2. The values of tropospheric zenith delay and the mapping function are set to the maximum values based on the UNB tropospheric delay model for conservative analysis. Consequently, the possible maximum value of the height difference can be calculated as approximately 1,172 m. Therefore, the height difference should not exceed 1,172 m in order to neglect the third and higher order terms of kΔh. Rearranging Equation (15) using the previously defined parameters a, d and e yields a linear equation with nonlinear constraints as
$$\eqalign{& \Delta {T_s} \equiv a + {T_z}({h_1}) \cdot m({h_2}) \cdot \left\{ {d\Delta h + e\Delta {h^2}} \right\} \cr & where,\;d = - k \lt 0,\;e = \displaystyle{1 \over 2}{d^2}} $$
where the constraint of a negative d is due to the relation between ΔT S and Δh: The larger the value of the tropospheric delay, the lower the height is. Consequently, the signs of ΔT S and Δh should be the opposite. The nonlinear constraint in Equation (17) is led from the relations between k and both d and e. Accordingly the two constraints on the parameters d and e reflect the characteristics of the tropospheric delay. Since the height-related interpolation methods were proposed to model atmosphere-related errors, the terms on the right side of Equation (17) can replace the height-related terms in Equation (9) as
$$\eqalign{& {g^i}(\bar p) \equiv {a^i} + {b^i} \cdot \Delta x + {c^i} \cdot \Delta y + T_z^i ({h_1}) \cdot {m^i}({h_2})\left\{ {{d^i} \cdot \Delta h + {e^i} \cdot \Delta {h^2}} \right\} \cr & where,\;{d^i} = - {k^i} \lt 0,\;{e^i} = \displaystyle{1 \over 2}{d^{i2}}} $$
which includes a nonlinear constraint. Therefore, the parameters ai, bi, ci, di , and ei cannot be determined using the LSE or the MNE solution. These unknowns should be solved through nonlinear programming (Byrd et al., Reference Byrd, Hribar and Nocedal1999).
2.4. Definitions of Various Interpolation Methods for Performance Evaluation
This section defines names of the previously introduced height-related interpolation methods. The LIM and LSM are used as base interpolation methods that only consider the coordinates in horizontal directions and bias. The additional terms, such as the first and second orders of height differences can be further considered in the base interpolation method. Moreover, the constraints can be given on the parameters of the first or second order of height difference. These additional considerations are indicated as ‘h’, ‘h2’, and ‘con’, respectively. Taking all the considerations described above, possible interpolation methods can be defined as
$$\matrix{{\rm LIM:} \hfill & {V = b\Delta x + c\Delta y} \hfill \cr {\rm LSM:} \hfill & {V = a + b\Delta x + c\Delta y} \hfill \cr {{\rm LIM}{\rm h:}} \hfill & {V = b\Delta x + c\Delta y + d\Delta h} \hfill \cr {{\rm LIM}\;{\rm h}\;{\rm con:}} \hfill & \matrix{V = b\Delta x + c\Delta y + {T_z} \cdot m \cdot d\Delta h \hfill \cr \qquad where,\;d \lt 0 \hfill} \hfill \cr {{\rm LSM}\;{\rm h:}} \hfill & {V = a + b\Delta x + c\Delta y + d\Delta h} \hfill \cr {{\rm LSM}\;{\rm h}\;{\rm con:}} \hfill & \matrix{V = a + b\Delta x + c\Delta y + {T_z} \cdot m \cdot d\Delta h \hfill \cr \qquad where,\;d \lt 0 \hfill} \hfill \cr {{\rm LIM}\;{\rm h2:}} \hfill & {V = b\Delta x + c\Delta y + d\Delta h + e\Delta {h^2}} \hfill \cr {{\rm LIM}\;{\rm h2}\;{\rm con:}} \hfill & \matrix{V = b\Delta x + c\Delta y + {T_z} \cdot m \cdot \left( {d\Delta h + e\Delta {h^2}} \right) \hfill \cr \qquad where,\;d \lt 0,\;e = \displaystyle{1 \over 2}{d^2} \hfill} \hfill \cr {{\rm LSM}\;{\rm h2:}} \hfill & {V = a + b\Delta x + c\Delta y + d\Delta h + e\Delta {h^2}} \hfill \cr {{\rm LSM}\;{\rm h2}\;{\rm con:}} \hfill & \matrix{V = a + b\Delta x + c\Delta y + {T_z} \cdot m \cdot \left( {d\Delta h + e\Delta {h^2}} \right) \hfill \cr \qquad where,\;d \lt 0,\;e = \displaystyle{1 \over 2}{d^2} \hfill} \hfill \cr}.$$
Table 1 shows the method for estimating parameters for each described interpolation method when the number of reference stations is set to four.
Table 1. Methods used for parameter estimation.
When interpolating dispersive and non-dispersive corrections for the user location instead of L1 and L2, the dispersive correction is modelled using the LSMs that do not have height-related terms. This is because the ionospheric delay is not significantly affected by heights, unlike tropospheric delay. This can be verified by the tropospheric delay model and ionospheric product provided by the International Global Navigation Satellite System (GNSS) Service (IGS). For users located at various heights in South Korea, the tropospheric and ionospheric delays are generated for 12 hours from 4:44:10 GMT on 20 May 2014 using the UNB tropospheric delay model and the Total Electron Contents (TEC) estimated by the IGS using GPS data collected from the reference stations distributed worldwide. The differences between the atmospheric errors of all visible satellites at each height and zero height are shown in terms of height differences in Figure 3.
Figure 3. (a) Tropospheric and (b) ionospheric delay variations for various height differences.
As seen in Figure 3 (a), the variations of the tropospheric slant delays are larger for the greater height difference and lower elevation angle. More specifically, when the difference of heights of two reference stations becomes approximately 100 m, the tropospheric slant delay difference has already reached a measure of over 10 cm for the elevation angle of 20°. On the other hand, the difference of the ionospheric delay over the height difference is relatively small as shown in Figure 3 (b): the value reaches only 2·5 cm even for a height difference of 10 km. These results support the notion that the height-related terms do not need to be considered when interpolating ionospheric delay, unlike tropospheric delay.
3. DATA AND RESULTS
Actual GPS data are collected to validate the performance of the proposed algorithms and compare the results to the conventional interpolation methods.
3.1. Data
Two reference station networks are selected in North Carolina (NC) and Vermont (VT), USA as shown in Figure 4. Among the stations of the Continuously Operating Reference Stations (CORS) that provide observation at one second intervals, these two networks are the only possible ones that can be composed of more than four stations whose baselines are 50–70 km and height differences are 200–1,000 m. The MARI station in the NC network and the VTC1 station in the VT network are chosen as the master stations since they provide best condition number for matrix W defined in Equation (11), which represents the network configuration, as shown in Table 2. The NCSW and VCAP station, which are in close proximity to the centre of each network are selected as the static users. The heights of the reference stations are shown in Figure 5.
Figure 4. The configuration of the selected networks in the (a) NC and (b) VT, USA for GPS data collection.
Figure 5. The ellipsoidal heights of the reference stations in the NC and VT networks.
Table 2. Condition numbers of the matrix W representing network configuration for different choice of master stations when considering height differences.
The GPS measurements were collected at the NC and VT networks from 13:00 to 15:14 on 23 March 2014 and from 20:00 to 23:00 on 19 April 2015 (GMT), respectively. The test periods are chosen arbitrarily to evaluate feasibility of the proposed method for independent data sets. The reference stations are equipped with Trimble NetR5 receivers. Mask angle is set to 10° and the ambiguities between master and auxiliary stations are fixed when the elevation angle is larger than 15–20° to avoid wrong fixing. The sky plot during the experiment is shown in Figure 6. The mean, maximum and minimum values of the Dilution of Precision (DOP) in horizontal and vertical directions of both networks are shown in Table 3. Corrections at the user location can be calculated using various interpolation methods. Since the users are set to reference stations with precisely known positions, the residual errors of users can be calculated after applying the corrections generated by the various interpolation methods and eliminating ambiguities which is estimated by batch process. The resultant residual errors are compared to see the advantage of the proposed method.
Figure 6. The sky plot at master station in (a) North Carolina and (b) Vermont, USA.
Table 3. Horizontal and Vertical DOP for the NC and VT networks.
3.2. Results
The L1 residuals of the user are calculated by various interpolation methods. For the conventional interpolation models without the height-related terms, the difference of the tropospheric delay between stations is calculated by the Black tropospheric delay model and applied to corrections. On the other hand, the height-related interpolation methods model the differences of tropospheric delays between two reference stations due to the height differences. The residual errors of the user for the selected interpolation methods and PRN are shown in Figure 7 and the RMS errors are shown in Table 4.
Figure 7. The residual errors of user for the network in (a) NC and (b) VT, USA.
Table 4. RMS errors of the L1 user residuals of the interpolation methods.
From Table 4 and Figure 7, the consideration of the first order height component slightly reduces the RMS error of the residual when comparing the LIM and LIM h. Comparing the LIM h and LSM h, the MNE solution provides the same performance as the LSE for both networks. The additional constraints to the first order of the height difference do not have much advantage on reducing the residual errors. Figure 8 shows the estimated parameter of the first order of height difference for LSM h and LSM h con for the NC network. It is clear that the parameter, d for the LSM h already almost satisfies the derived constraint on d even if the corresponding constraint is not given. The parameter d for LSM h and T z ·m·d for LSM h con are in the same level for comparison.
Figure 8. The estimated parameter, d for the LSM h and LSM h con for the NC network and observed satellites.
If we look at Table 4 again, the additional consideration of the second order of the height difference rather increases the RMS errors of the LIM h2 and LSM h2. In Figure 7 this is apparent for the LSM h2 for PRN26 for the NC network and PRN16 for the VT network. On the other hand, when the additional constraint is given on the second order of the height difference, the LIM h2 con and LSM h2 con further reduces the RMS errors by approximately 22% and 25% for both the NC and VT networks. In order to see the effect of the constraints on the second order of height difference, the estimated parameters of LSM h2 and LSM h2 con for the VT network are compared in Figure 9.
Figure 9. The estimated coefficients for the LSM h2 and LSM h2 con for the VT network and observed satellites.
In Figure 9, the estimated parameters of b and c for LSM h2 and LSM h2 con have a similar level of values. However, there are significant differences for other parameters a, d and e, and they are represented in different ranges of axes in Figure 9. The estimated parameter, d for LSM h2, has values of an order of 10−9 m/m, which is negligible even if it is multiplied by the height difference of an order of 102 m. Meanwhile, the parameter e dominates when considering the effect of the height difference that can be cm-level values in the correction level. The non-dispersive correction is modelled mainly by the b, c and e parameters for the LSM h2 since the corresponding coordinate differences, Δx, Δy and Δh 2 are simply larger than the value of Δh and 1 for parameters d and a. This is unnatural for our observation on the characteristics of the tropospheric delay in Section 2.3. Consequently, the RMS errors of the LIM h2 and LSM h2 are rather increased compared to the LSM h. For LSM h2 con, the parameter for the first order height difference term, T z ·m·d dominates when considering the effect of the height difference and the parameter e for LSM h2 con has smaller values of an order of 10−8 m/m2 due to the nonlinear constraints.
The user residual errors are now transformed to the position domain by the weighted LSE using an observation matrix H as
$$\eqalign{\delta {\bar x} & = \left( {{H^T}H} \right)^{ - 1}{H^T}\delta {\bar r} \ where,\;\delta {\bar r} \equiv \left[ {\matrix{ {^2 {\nabla ^1}{r_{L1,user}}{ - _{user}}{\Delta _{mas}}^2 {\nabla ^1}\delta {{\hat N}_{L1}}{\lambda _{L1}}} \cr \vdots \cr {^n {\nabla ^1}{r_{L1,user}}{ - _{user}}{\Delta _{mas}}^n {\nabla ^1}\delta {{\hat N}_{L1}}{\lambda _{L1}}} \cr}} \right] \cr & \qquad H \equiv \left[ {\matrix{ {{ - ^2}{\nabla}^1{{\bar l}_{user}}} \cr \vdots \cr {{ - ^n}{\nabla}^1{{\bar l}_{user}}} \cr}} \right].}$$
The estimated ambiguity,
$_{user} {\Delta _{mas}}^n {\nabla ^1}\delta {\hat N_{L1}}$
is eliminated from the residual error in Equation (5) to form a vector
$\delta \bar r$
. The symbol of
$\bar l$
indicates a line-of-sight vector and
$\delta \bar x$
represents the position error vector. The transformed horizontal and vertical position errors are expressed in Figure 10. Table 5 represents the horizontal 2DRMS errors and vertical 2RMS errors of user for all the interpolation methods.
Figure 10. The horizontal and vertical error for (a) the NC network and (b) the VT network.
Table 5. Horizontal (2DRMS) and vertical (2RMS) positioning error of all interpolation methods.
From Table 5, the LSM h2 con provides the most accurate vertical (2RMS) positioning result since the tropospheric delay errors are related to the vertical positioning accuracy. The LSM h2 con improves the vertical positioning accuracy by approximately 4% and 7% compared to the LSM for the NC and VT network, respectively. Compared to the LIM, it is improved by approximately 43% and 10% for both networks. The conventional LIM h and LSM h have similar or slightly improved vertical positioning accuracy compared to the LSM for both networks. Both the horizontal and vertical position results are rather degraded for the LSM h2 methods even if the MNE is used. The number of observations is less than the number of unknowns for both LSM h2 and LSM h2 con methods, however, the additional constraints based on the characteristics of the tropospheric delay variation over height can be helpful for estimating parameters.
The horizontal and vertical positioning accuracies shown above are of a statistical character. In reality, for users, the epoch-by-epoch positioning errors are important as well because they are the actual position results that users experience. Therefore, the maximum horizontal and vertical errors during the observation time span are also calculated in Table 6.
Table 6. The maximum horizontal and vertical errors of all interpolation methods.
Compared to the conventional LSM method, the results prove that the LSM h2 con is also effective in improving epoch-by-epoch positioning accuracy, especially for vertical direction, by approximately 8% and 3% for the NC and the VT network, respectively. In addition, the LSM h2 con gives much improved performance compared to the LIM method.
The different levels of performance of the two networks in the residual errors shown in Table 4 are suspected to be due to the different sizes of the network (Wang et al., Reference Wang, Feng, Higgins and Cowie2010). In addition, the overall positioning accuracies of the VT network are worse than those of the NC network because of the larger residual errors and maximum values of the DOP. Nevertheless, the proposed method, LSM h2 con, is effective for both networks in reducing residual errors and positioning accuracies especially in the vertical direction.
4. CONCLUSION
This paper proposed new approaches that consider the height terms of reference stations to interpolate corrections to the user location. The proposed approach employs the MNE to account for the height-related terms for the interpolation methods when the number of reference stations in a network is fixed. In addition, in order to reflect the actual characteristics of the tropospheric delay, derived constraints are given on the various methods to make modifications of the conventional LIM and LSM. The actual GPS data for two networks are collected to evaluate the levels of performance of the conventional and proposed methods. The results indicate that the interpolation methods using the MNE do not produce any significant advantages compared to the LSE when the results of LIM h and LSM h are compared. Even adding the constraint on the first order of the height term to the interpolation methods does not help improve the level of performance of the user in this case because the constraint is already nearly satisfied without giving it intentionally. However, if the constraints on both of the first and second order of the height terms are imposed to the LSM h2, the vertical positioning accuracy of the user is improved by maximum 43% and 7% compared to the conventional LIM and LSM respectively. This method is also expected to have an advantage in terms of positioning accuracy using the epoch-by-epoch basis, as it showed maximum 40% and 8% improvements for the maximum vertical errors in these GPS data sets compared to the LIM and LSM respectively. Moreover, when a sufficient number of reference stations is possible, the levels of performance of the LSM h, LIM h2, LSM h2, LSM 2 h con and LSM h2 con might be improved since the parameters of the methods can be calculated on the LSE basis.
The proposed interpolation method with constraints LSM h2 con will be an effective solution when precise positioning accuracy is required for Network RTK users located in mountainous areas since the heights of the reference stations and users in those areas can differ by hundreds of metres. These types of applications include monitoring the elevations of mountains, automation of agriculture in mountainous areas, and so on.
ACKNOWLEDGEMENT
This research was supported by a grant from “Development of GNSS based Transportation Infrastructure Technology (06-A03)” funded by the Ministry of Land, Infrastructure and Transport of Korean government, contracted through SNU-IAMD at Seoul National University.
