Improvement Analysis of a Height‐Deviation Compensation‐Based Linear Interpolation Method for Multi‐Station Regional Troposphere

In network real‐time kinematic positioning of multi‐reference station, the spatial and temporal distribution of tropospheric delay is affected by both horizontal and elevation. The traditional modeling strategy of regional troposphere takes more consideration of the horizontal factor, and the incomplete consideration of the elevation factor will lead to the problem of reduced modeling accuracy, especially in the face of the scene with large regional height deviation. Based on the traditional linear interpolation method (LIM), a simple and effective height‐deviation compensation‐based linear interpolation method (HCLIM) for regional tropospheric is proposed. The modeling accuracy of troposphere and the positioning accuracy of user RTK in large height deviation region are significantly improved. The method was verified based on six experimental subnets with large height deviations from a provincial continuously operating GNSS reference stations network in central China. The results showed that: For GPS satellite modeling, compared with the traditional LIM method, the average modeling accuracy improvement rate of HCLIM method is (84.5%, 75.5%, 59.3%, 26.7%) in the elevation angle range of (10–30°/30–40°/40–50°/50–90°). For BDS satellite, the average modeling accuracy improvement rate of HCLIM method in the above four elevation angles is (83.3%, 70%, 50%, 23.5%). For the positioning performance of user RTK, The horizontal positioning accuracy and RTK fixing rate were similar under the two methods, while HCLIM method showed only slight improvement. However, in the U direction, LIM method showed obvious systematic bias, while HCLIM method showed consistent positioning accuracy, which was improved to 82.8% compared with LIM method.

characteristics of simple technology implementation, stable service performance and short initialization time, so it has been widely used in engineering applications.Network RTK technology models are numerous, the basic principle is similar.As a branch of network RTK technology, virtual reference station (VRS) technology has been widely promoted because of its simple implementation principle and good compatibility (Hu et al., 2003;P. Wang et al., 2022).Figure 1 shows the basic service flow of VRS mode: the user first uploads its general location to the server.Then the server modeled and interpolated the atmospheric information at the user's location according to the atmospheric information fixed under the extracted reference station baseline.Finally, a virtual observation value is constructed near the user and sent to the user.The user can construct an ultra-short baseline RTK by using his own observation value and the virtual observation value received by the server side, so as to obtain the fixed solution position with centimeter-level accuracy.
As can be seen from Figure 1, selecting an appropriate method to model and interpolate accurate regional atmospheric delays is an important step, which directly affects the quality of virtual observations broadcast to the user.Therefore, scholars have proposed a variety of modeling methods linear combination model (Han & Rizos, 1996) is used to model satellite orbit errors, ionospheric and tropospheric delays and other GPS spatial errors.The distance-based linear interpolation model (DIM) (Gao et al., 1997) interpolates the ionospheric delay at the user's location by inverse distance weighting between the user and the reference station.Wanninger also first proposed another location-based linear interpolation model (LIM) (Wanninger, 1995) to model the regional ionospheric delay, which requires at least three reference stations in the region of the rover station to be modeled.The low-order surface model (LSM) (Fotopoulos & Cannon, 2001) built on the principle of partial derivatives of horizontal and elevation directions can be used to model distance-related errors in the reference station network, including ionospheric and tropospheric errors.In addition, the LIM model is essentially equivalent to the two-dimensional interpolation form of the LSM model.The Least squares collocation model (LSC) (Odijk et al., 2000) based on the observed quantity covariance matrix and the mutual covariance matrix between the observed vector and interpolation vector can also model the distance-dependent error within the network (Al-Shaery et al., 2011;Dai et al., 2003) comprehensively evaluated the accuracy of these interpolation algorithms, and the results showed that, in addition to the DIM model, the interpolation effect was slightly worse, the performance of various interpolation algorithms was close, and the LIM model was slightly better in all aspects.
Compared with ionospheric delay, tropospheric delay is highly dependent on elevation in GNSS precision positioning.In other words, tropospheric delay, as a spatial correlation error of GNSS baseline, is not only affected by baseline length (with the increase of baseline length, the correlation of tropospheric delay will decrease and cause large residual error, thus affecting the solution of ambiguity parameters and user positioning accuracy), but also affected by the elevation difference at both ends of the baseline (Wielgosz et al., 2011).Some scholars have pointed out that the correlation coefficient between baseline tropospheric delay and elevation component can reach about 0.9, and the correlation is especially obvious in areas with large height difference due to severe topographic relief (Wielgosz et al., 2011;Yin et al., 2008).In addition, the magnitude of tropospheric delay itself is also highly correlated with the elevation angle of the satellite.Specifically, the smaller the elevation angle of the satellite, the longer the path through the troposphere of the satellite ranging signal, and the larger the magnitude of the corresponding delay (Wielgosz et al., 2011).However, the above mentioned modeling methods are basically plane modeling based on plane coordinates and only consider the horizontal distribution of tropospheric delay.In fact, when several reference stations build baselines for tropospheric region modeling, if the height difference between each station is too large, it will obviously reduce the accuracy of elevation direction, and even cause systematic errors.In view of this problem, relevant scholars have also conducted research and given some solutions.LSM method can improve the modeling accuracy of tropospheric delay in elevation direction by adding first-order or second-order elevation factors into the modeling equation (Fotopoulos & Cannon, 2001).However, since the original troposphere is not adjusted directly, the effect is not good in the scenario of large height difference.Moreover, due to the increase of equation coefficients, at least four stations or more reference stations are required to participate in the modeling.It does not meet the requirement of at least three sites in the delaunay triangulation network in VRS mode.To solve the problem of LSM (Song et al., 2016), proposed a modified LSM strategy, that is, the zenith tropospheric delay (ZTD) was constructed as an exponential function related to height difference, and then the exponential function was used to replace the first-order and second-order elevation factors in the traditional LSM.The modeling accuracy was improved by about 10% in the elevation direction, and the modeling could be completed with at least three reference stations.But this approach is a little more complicated to implement; In literature (Zhao et al., 2021), an ECDIM model is proposed on the basis of the traditional DIM model, that is, the double difference (DD) tropospheric delay correction between the rover station and the master reference station is added to the traditional modeling strategy, and good results are achieved, but the improvement in the elevation direction of the RTK is not obvious.Based on the traditional LIM model, a MLIM modeling method was proposed in literature (Pu et al., 2021;Y. Wang et al., 2023), that is, the DD tropospheric delay of the baseline was compensated by the model first, and then the traditional LIM modeling was carried out, achieving a good effect.However, the modeling residuals and user positioning accuracy were not discussed in detail (Zheng & Feng, 2005) uses reference network data to interpolate residual ZTDs at a user's location.Then, by using ZTD empirical tropospheric model correction and ZTD residual correction, the influence of inaccurate tropospheric empirical model can be greatly reduced.However, this method mainly targets the Regional Area Differential GPS Positioning (RADGPS) scenario.
To sum up, in view of the problem of reduced modeling accuracy of tropospheric elevation component caused by height-deviation between reference stations, the above mentioned modeling methods either do not consider the influence of height-deviation, or the correction effect is not good, or the implementation is too complex.In this paper, based on the studies of existing scholars and LIM modeling method, we propose a height-deviation compensation-based linear interpolation method (HCLIM) method that takes into account height-deviation compensation.The method is simple and effective, and the specific steps are as follows: First, the tropospheric delay extracted from the baseline was compensated for the part related to the height difference with the global pressure and temperature model (GPT2) (Lagler et al., 2013) and the Saastamoinen model (Saastamoinen, 1972).That is, the tropospheric delay of each participating modeling baseline is compensated to a uniform elevation plane (the elevation plane of the master reference station can be uniformly adopted); Then, the traditional LIM model is used to model and interpolate the DD troposphere delay between the rover station and the master reference station.Finally, the interpolated tropospheric delay is compensated by a prior high-precision tropospheric model to the elevation plane of the rover station, which is the final interpolated tropospheric delay at the position of the rover station (VRS position).The internal interpolation of tropospheric delay can be used for the construction of virtual observations in Figure 1 and subsequent user RTK positioning.The structure of this paper is as follows: In Section 2, first, the baseline solution strategy of reference station network and the high-precision DD tropospheric delay extraction method under fixed solution are presented, and then the implementation details of traditional LIM method and improved LIM (HCLIM) method are presented respectively.In Section 3, the experimental data collected to verify the algorithm are given, and the parameter estimation strategies of the server and the user and the corresponding random model settings are given.In Section 4, three sets of experiments are carried out according to the collected data, and all experiments compare the modeling effects of HCLIM method and LIM method.Finally, the conclusion and suggestions and the prospect of the next step are given.

The Fixing of the Baseline Ambiguity of the Reference Station Network and Tropospheric Extraction
For the network RTK, the solution of medium and long baselines between reference stations is the foundation for the follow-up work.The key of baseline calculation is to solve the ambiguity under the constraint of precise coordinates.For a long time, many scholars have conducted a large number of studies on this, among which the combination of ionosphere-free (IF) observations and wide lane (WL) observations is the most commonly used solution strategy (Liu et al., 2021).That is, the WL combination has a long wavelength and is less affected by observation noise, so it is easy to be fixed.The WL ambiguity is fixed first, and then the narrow lane (NL) ambiguity can be determined by combining with the IF combination ambiguity (Xu et al., 2018).However, the IF combination not only eliminates the ionospheric first-order term, but also amplifies the noise.In practical applications, it takes a long time for initialization to complete the separation of NL ambiguity from various noises.Therefore, the convergence rate of ambiguity is still slow.Moreover, prior ionospheric corrections or external constraints cannot be added to the functional model after linear combination of observations.In this paper, the DD non-combination function model is adopted, and the ambiguity of L2 frequency is transformed into the difference between L1 frequency ambiguity and WL ambiguity.Then L1 and WL ambiguity are directly taken as parameters to be estimated in the equation (M.Zhang et al., 2017).The non-combination strategy has two advantages: (a) L1 and WL ambiguities are estimated at the same time, and it is very convenient to continue to search for other ambiguities and update the normal equation after partial ambiguity is fixed.Compared with IF combination strategy, the complexity of the model is greatly reduced; (b) This non-combinative model retains ionospheric parameters, which facilitates the addition of the external ionospheric weighted constraint strategy (Mi et al., 2019;Odolinski & Teunissen, 2017;Wielgosz, 2011).Previous studies have shown that ionospheric weighted constraint strategy is one of the most effective methods for fast relative positioning of GNSS (Odijk, 2000b;Odolinski & Teunissen, 2017;Paziewski, 2016;Paziewski & Wielgosz, 2014;Wielgosz, 2011;R. Zhang et al., 2022), that is, the prior ionospheric correction information and its variance are introduced into the function model as virtual observations, which can greatly improve the fixed effect of ambiguity.The uncombined DD observation equation based on dual-frequency observation data is as follows: (1) where P and L represents pseudo-range and carrier phase observations, L = λφ; p and q represent the reference satellite and non-reference satellite identifier, r and b represent the rover station and reference station identifier, respectively; ρ is the geometric distance between the receiver and the satellite (for the network RTK, since the coordinates are known precisely, this parameter is a known value); I 1 represents the slant path ionospheric delay at f 1 ; w and m represent the zenith tropospheric wet delay of the corresponding stations and wet delay mapping function of the single-difference between the satellites; λ is the carrier wavelength; B 1 is the ambiguity of L1 frequency (NL ambiguity), B W is the ambiguity of WL, B W = B 1 − B 2 ; ε P and ε L represent the noise of pseudo-range and carrier observation, respectively.In addition, in order to further improve the fixed effect of ambiguity parameters, virtual atmospheric weighted constraint equations were added for ionospheric and tropospheric parameters to assist the rapid convergence of atmospheric parameters (Odijk & Teunissen, 2010).For the specific estimation strategies of each parameter and the setting of the random model, please see the Net-RTK item in Table 1 of Section 3.2 of this paper.The baseline solution error observation equation after unified naturalization is as follows: According to Equation 2, the equation can be solved by Kalman filter.It should be noted that the solution order of ambiguity still adopts the strategy of WL first and NL later.After the WL ambiguity is fixed, Equation 3 is used to update the normal equation (Teunissen, 2000;L. Wang et al., 2017), and then the LAMBDA method (Chang et al., 2005;Teunissen, 1995) is used to fix the NL ambiguity.Because the float NL ambiguity updated with fixed WL ambiguity has higher precision and is easier to fix. where , and  B represent the float WL ambiguity, fixed WL ambiguity, float NL ambiguity, and fixed NL ambiguity, respectively; Q refers to the corresponding covariance matrix.
After the NL ambiguity parameters are fixed, the high-precision DD tropospheric delay on the slant path of each fixed satellite can be directly extracted according to Equation 4 (Liu et al., 2021):

Traditional Tropospheric Delay Modeling Strategy (LIM)
According to Formula 4, the DD tropospheric delay of fixed satellites under each baseline can be extracted and then the regional atmosphere modeling can be carried out.Since LIM method uses the least squares estimation method (LSQ) to obtain the modeling coefficients, the solution of the two plane modeling coefficients usually requires the participation of at least three reference stations.Three reference stations around the location of the rover station can construct three baselines, among which the reference station closest to the rover station can be selected as the master reference station to participate in the generation of the final virtual observation value.
Based on the traditional LIM model, the formatted formula for calculating the DD tropospheric delay of a non-reference satellite q under a certain baseline is as follows (Cui et al., 2018): where (x r , y r ) and (x b , y b ) are respectively the plane positions of the rover station and reference station of a baseline involved in modeling under Gaussian projection; x rb = x r − x b and y rb = y r − y b are the Gaussian plane coordinate difference of the baseline; a T and b T are the plane atmosphere modeling coefficient of a non-reference satellite q; The meanings of other upper and lower indexes are consistent with those in Formula 1. Suppose that satellite q has selected n baselines to participate in atmospheric modeling within the regional scope, then the regional tropospheric modeling formula of satellite q based on the LIM model is shown as Formula 6: According to Equation 7, the DD tropospheric plane modeling coefficient under LIM model can be solved.By putting the plane modeling coefficients obtained from Formula 7 and the coordinates of the rover station and the master reference station into Formula 5, the DD tropospheric delay between the rover station and the master reference station can be interpolated: where,  Δ∇   represents the interpolated DD tropospheric delay between the rover station and the master reference station.M and V All represent the identification of the master reference station and the rover station.Δx MV and Δy MV represent the difference of plane coordinates between the master reference station and the rover station.

Improved Tropospheric Delay Modeling Strategy (HCLIM)
In GNSS data processing, Compared with the ionosphere, the temporal and spatial characteristics of tropospheric delay are significantly different, which is influenced by both horizontal and elevation.As can be seen from Section 2.2, the traditional LIM modeling strategy only uses plane coordinates for interpolation, ignoring the influence of elevation factors.When the height difference between the reference station and the rover station is large, the modeling accuracy will obviously be reduced.The large modeling residual may bring systematic deviation to the user coordinates, especially the elevation direction component, and even affect the fixed effect of the ambiguity in the user RTK solution.Therefore, we propose to add a tropospheric height-deviation compensation mechanism on the basis of LIM modeling, that is, when the height-deviation between the rover station and the reference station is large, a prior high-precision tropospheric model should be used in the modeling and interpolation process to compensate the deviation.
Assume that three reference stations M, A, and B construct a minimum delaunay triangle participating in the tropospheric modeling of the region, where site M is the master reference station of the rover station V (usually the nearest station is the master reference station).Assuming that the height difference between the rover station V and the three reference stations is large, the modeling and interpolation process of troposphere needs to compensate the height-deviation.As shown in Figure 2, the flow chart of HCLIM method is given with reference to VRS mode, which can be roughly divided into the following three steps: Step 1: The tropospheric delay of each baseline extracted from Equation 4was compensated for the part related to the height difference by using the global temperature and pressure model GPT2 (Lagler et al., 2013) and the Saastamoinen model (Saastamoinen, 1972).That is, the tropospheric delay of each participating modeling baseline is compensated to a unified elevation plane (the elevation plane of the master reference station is uniformly adopted here).Take the baselines established by sites A and B as an example, as shown in Equation 9: where  Δ∇   (Δℎ) represents the DD tropospheric delay of satellite q directly extracted from Equation 4 after the ambiguity of baseline AB is fixed (p is the reference satellite).
where, h M and h B denote the elevation of the master reference station M and reference station B; ZTD B (h M ) and ZTD B (h B ) represent the ZTD of station B at elevations h M and h B calculated by the global pressure and temperature model (GPT2) and the Saastamoinen model, respectively;     (ℎ ) and     (ℎ ) respectively represent the tropospheric projection function of satellite q and p at the plane coordinates of site B with the altitude h M .
(ℎ) and     (ℎ) respectively represent the tropospheric projection function of satellite q and p at the plane coordinates of site B with the altitude h B .The solution of tropospheric projection function is based on GMF projection function (Böhm et al., 2006); The meaning of symbols in Equation 13 can be referred to Equation 12.It should be pointed out that although the ZTD calculated by the prior model itself has model error, the difference of ZTD at different heights of the same plane position obtained by the prior model can effectively eliminate the tropospheric systematic error caused by height-deviation (Song et al., 2016;Zhao et al., 2021).Equations 12 and 13 are substituted into Equations 10 and 11, and then Equations 10 and 11 are substituted into Equation 9,  Δ∇ T   (ℎ ) is the DD tropospheric delay of baseline AB eventually compensated to the elevation plane h M of the master reference station, as shown in Equation 14: With reference to baseline AB, the same tropospheric height-deviation compensation operation is performed on the other two baselines, that is, the alternative group of modeling baselines that are all naturalized as elevation plane h M is obtained.
Step 2: Based on the alternative group of modeling baselines obtained in step 1, the traditional LIM model given in Section 2.2 is adopted for tropospheric modeling.Then the DD tropospheric delay between the rover station V and the master reference station M is interpolated.At this time, the interpolated tropospheric delay is the interpolation uniformly normalized to the elevation plane h M where the master reference station is located, as shown in Equation 15: where h V and h M are the elevation values of the rover station V and the master reference station M respectively, and Δh MV is the height difference between the two stations.
Step 3: The DD tropospheric delay interpolated in step 2 based on the elevation plane of the master reference station is compensated again by the prior high-precision tropospheric model to the elevation plane h V of the rover station V.With reference to Equation 14, the compensated DD tropospheric delay between the rover station and the master reference station is shown in Equations 16 and 17: where  Δ∇ T   (ℎ  ) is interpolated DD tropospheric delay after two elevation compensation (one compensation before modeling and one compensation after interpolation).This tropospheric value can be used for the construction of virtual observations in Figure 1 and subsequent user RTK positioning.It has effectively eliminated the systematic error of regional tropospheric modeling caused by the height-deviation of the reference station, especially for the low elevation angle satellite, the tropospheric modeling accuracy has been greatly improved.Experimental demonstration will be carried out in the following chapters.

Data Sets
To evaluate the service effect of the HCLIM tropospheric modeling strategy in the large height difference region.This paper collected the data of a provincial CORS network in central China.Figure 3 shows the general distribution diagram and elevation statistics of all reference sites (different colors represent different elevations).The whole region presents a significant terrain distribution of high in the west and low in the east.In this paper, several sites with high elevation difference in the west are selected to construct six experimental subnets.Figures 4 and 5 respectively show the detailed geographical distribution and elevation statistics of the six experimental subnets.It can be seen from Figure 5 that the height difference between the six rover stations and the reference stations within each subnet varies greatly.The height difference between the six rover stations and the reference stations within each subnet ranges from 1,060 m (YC06 station) to 220 m (SY12 station).
The data collection period in this paper is 21 January 2022/00:00:00-24:00 (GPST).The cutoff elevation angle is set at 10° and the sampling interval is 10 s.In addition, all stations can receive GPS/BDS-2/BDS-3 satellites.In this paper, GPS + BDS dual system was used for algorithm processing.Figure 7 shows the skyplot and PDOP value sequence in the data collection period by taking YC06 site as an example.
In addition, the ionospheric extension is another key atmospheric correction information that affects the differential positioning of the network RTK region.In order to eliminate the interference of ionospheric delay on the experimental analysis results of tropospheric modeling, the ionosphere is relatively quiet during the period of experimental data selected in this paper.Figure 6 shows the I95 ionospheric index (Wanninger, 2004) of the covered area calculated based on the joint networking of 18 reference stations and six rover stations in Figure 4, which can directly reflect the ionospheric activity degree of the reference station area during the whole data acquisition period: If the I95 index is  less than 1, it means that the regional ionosphere is relatively calm in this time period; if the I95 index is between 1 and 2, it means that there is a disturbance of small magnitude in the regional ionosphere; if the I95 index is between 2 and 4 or greater than 4, it indicates that there is a strong or very strong disturbance in the regional ionosphere.As can be seen from Figure 6, the ionosphere is relatively stable throughout the data testing period, and high precision ionospheric modeling results can be obtained, which will not affect the analysis of experimental results of tropospheric modeling.

Processing Strategies
Table 1 shows the specific parameter estimation and random model setting strategies of the server (Net-RTK) and user (VRS-RTK) in the baseline calculation process.As can be seen, the main difference between the server and the user is that the server needs to estimate ionospheric and tropospheric parameters due to the long baseline.The client does not need to consider the ionosphere and troposphere parameters, because it has been modified when the server generates the virtual reference value, but it needs to estimate the coordinate parameters.

Analysis of Modeling Effect of Tropospheric Atmospheric Delay
In order to evaluate the modeling effect of the proposed HCLIM modeling method in this paper, different modeling results were compared based on the data of six experimental subnets in Figure 4. Figures 8-13 shows the modeling results of six rover stations (YC06, ES01, XF12, YC21, XF06, SY12).For multi-angle comparative analysis, for each rover station, different data solving periods, different elevation angle types (from left to right, high elevation angle satellites, medium elevation angle satellites low elevation angle satellites), and satellites of different systems (the top three GPS satellites, the bottom three BDS-2/BDS-3 satellites) are selected in this paper.In addition, the six simulated rover stations selected in Figure 4 are actually reference stations of CORS network whose coordinates are precisely known in Figure 3.The server can directly extract high-precision tropospheric results with fixed ambiguity by using baseline solution mode (according to Formula 4), and the tropospheric results directly extracted can be used as the tropospheric reference value in Figures 8-13.
From Figures 8-13, the modeling results of different stations in different time periods, different elevation angle types, and different satellite systems are compared.The following rules can be obtained: (a) The size of tropospheric delay is negatively correlated with the elevation angle, that is, with the increase of the elevation angle, the tropospheric delay tends to 0, and vice versa.This is understandable, because the smaller the elevation angle of the satellite, the longer the path through the troposphere of the satellite ranging signal, the greater the delay of the ranging signal, the greater the calculated tropospheric delay magnitude; (b) The results of the two modeling methods differ greatly under different types of elevation angle.Specifically, the results of the two modeling methods are basically in good agreement with the true values under high elevation angle satellites.In the medium elevation angle satellite, the deviation between LIM modeling result and true value will increase with the elevation angle decreasing.Under low elevation angle satellite, there will be a large deviation between LIM modeling results and true values.However, the modeling accuracy of HCLIM is basically balanced under any elevation angle type, and can maintain a good agreement with the true value.Take the YC06 site shown in Figure 8 as an example.The maximum deviation between LIM modeling results and true values of G32 satellite and C20 satellite (high elevation angle) is about 5 and 3 cm respectively.The maximum deviation between LIM modeling results and true values of G25 and C32 satellites (medium elevation angle) is about 10 and 7 cm respectively.The maximum deviation between LIM modeling results and true values of G12 satellite and C38 satellite (low elevation angle) is about 16 and 25 cm respectively.However, the modeling results of HCLIM of several types of satellite in YC06 site are basically consistent with the true values, maintaining a maximum deviation of about 3 cm.
The modeling error RMS of each satellite under each station of Figures 8-13 under the two modeling methods are statistically analyzed (The error here is the result of the two methods compared to the reference value), and the corresponding results are shown in Figure 14.The following rules can be observed: (a) First, whether for GPS system or BDS system, the modeling residual RMS under LIM method shows an increasing trend with the lower satellite elevation angle.However, there is no obvious correspondence between the modeling residual RMS and the elevation angle under HCLIM method, which also indicates the consistency of the modeling accuracy of HCLIM.(b) Second, for satellites of any elevation angle type under 6 rover stations, the residual RMS of modeling under HCLIM method are all smaller than LIM.Especially for satellites with low elevation angle, the modeling accuracy is improved significantly.This is consistent with the rule of (a), indicating the superiority of HCLIM method in improving the tropospheric modeling accuracy of low elevation angle satellites.(c) Finally, compared with the other five rover stations, the height-deviation between YC06 station and other reference stations in its subnet is the largest (up to 1,060 m).The modeling error RMS of each type of satellite under LIM method is significantly larger than that of the same type satellite of other rover stations, but the modeling accuracy of HCLIM method is basically consistent with that of other stations.The superiority of HCLIM method in improving the tropospheric modeling accuracy in the scenario of large height difference of reference station is also demonstrated.

Analysis of Modeling Residuals of Tropospheric Atmospheric Delay
Based on the analysis in Section 4.1, it can be seen that tropospheric delay magnitude and modeling accuracy are significantly corresponding to elevation angle.In order to further compare and analyze the improvement of modeling effect of HCLIM method compared with LIM method.Figure 15 shows the tropospheric modeling errors of all GPS satellites at each rover station in the whole data test period (24 hr) under the two modeling methods according to the elevation angle.The RMS statistical diagram of modeling errors of all satellites in each elevation angle range (10-30°/30-40°/40-50°/50-90°) is also compared on the right side of Figure 15.lists the RMS statistical results of GPS satellite modeling errors under the two modeling strategies on the right of Figure 15.After averaging the mean values of all sites, the modeling accuracy improvement of HCLIM method compared with LIM method was given under different elevation angle ranges.Figure 16 and Table 3 show the comparison results of BDS2/BDS-3 satellites under different modeling strategies.
By combining Figures 15 and 16, Tables 2 and 3, we can see the following rules: (a) No matter the GPS system or BDS system, the modeling accuracy of LIM method is strongly correlated with the elevation angle, and the modeling residuals increase with the decrease of elevation angle.The modeling residuals of satellites with low elevation angle are especially large.Taking YC06 station as an example, the maximum modeling residuals of satellites below 20° can reach about 0.5 m.Compared with LIM method, the correlation between the modeling accuracy of HCLIM method and the elevation angle is not obvious.Only the satellite in the range of low elevation angle has a slightly larger modeling residual value, but the overall satellite modeling residual value fluctuates within the maximum 0.05 m.This is consistent with the analysis in Section 4.1; (b) It can be seen from the statistical results of satellite modeling error RMS of the two modeling strategies in different elevation angle ranges that, since the residual modeling error of LIM method for satellites with low elevation angle is large, HCLIM method can improve the modeling accuracy more obviously in the elevation angle range of 10°-30°.In the range of 50°-90° elevation angle, the modeling accuracy of the two methods is basically similar, but HCLIM method still has the advantage of accuracy.For GPS systems, The average residual RMS of LIM method and HCLIM method in the range of elevation angles (10-30°/30-40°/40-50°/50-90°) are (0.103 , 0.049, 0.027, 0.015 m) and (0.016 , 0.012, 0.011, 0.011 m), respectively.The accuracy improvement rates were 84.5%, 75.5%, 59.3%, and 26.7%, respectively.For the BDS system, The average residual RMS of LIM method and HCLIM method in the range of elevation angles (10-30°/30-40°/40-50°/50-90°) are (0.114 , 0.05, 0.028, 0.017 m) and (0.019, 0.015, 0.014, 0.013 m), respectively.The accuracy improvement rates were 83.3%, 70%, 50%, and 23.5%, respectively.(c) Compared with other rover stations, the height deviation between YC06 station and other reference stations in its subnet is the largest (1,060 m).The modeling error RMS of LIM method at each elevation angle range is significantly larger than that of other rover stations, but the modeling error RMS of HCLIM method at YC06 site is not significantly different from the results of other stations, which is consistent with the argument in Section 4.1.The effectiveness of HCLIM method in compensating the height-deviation of tropospheric delay is also demonstrated.

Analysis of User RTK Positioning Performance
The final service effect of network RTK needs to be reflected by user RTK positioning.The LIM and HCLIM methods can be used to interpolate the tropospheric atmospheric delay of the rover station under different modeling strategies, and then the virtual observation values under different modeling strategies can be constructed and sent to the rover station for positioning comparison of user RTK.Since the six simulated rover stations selected in Figure 4 are actually reference stations of CORS network with known exact coordinates in Figure 3, RTK results of each rover station can be directly compared with external coincidence accuracy.Figure 17 compares the time series of RTK external coincidence positioning errors of fixed solutions of each rover station in LIM mode and HCLIM mode.Table 4 gives statistics on positioning errors and RTK fixed rate in Figure 17.
As can be seen from the results in Figure 17 and Table 4: (a) The horizontal positioning accuracy and RTK fixing rate of each rover station under the LIM and HCLIM modeling methods are relatively high, and the results under the two methods do not show much difference.The average RTK fixing rate of the rover station under LIM and HCLIM methods and the average positioning error RMS in the E and N directions were (98.6%, 1.3 cm, 1.4 cm) and (99.4%, 1.0 cm, 1.1 cm), respectively.The positioning results of HCLIM showed slight improvement.(b) Compared with the results of horizontal positioning, the LIM method presents obvious systematic deviation in the orientation of elevation.The positioning RMS of 6 rover stations (YC06, ES01, XF12, YC21, XF06, SY12) in the U direction in LIM mode were (33.2 , 12.6, 8.1, 12.6, 10.1, 14.2 cm), and the average RMS was 15.1 cm.This is obviously because the LIM method does not consider the effect of height deviation when modeling the tropospheric atmosphere at the server end, thus introducing a large modeling residual for the user RTK solution.The positioning accuracy of HCLIM method in U direction is consistent.The positioning RMS of six rover stations in U direction are (2.0, 2.4, 3.2, 3.0, 2.4, 2.7 cm), and the average RMS is 2.6 cm.Compared with LIM, the accuracy is improved to 82.8%.(c) Compared with other rover stations, the height deviation between YC06 station and other reference stations in its subnet is the largest (1,060 m).The systematic deviation of LIM method in U direction is significantly greater than that of other stations (33.2 cm).However, there is no significant difference between the results of U orientation of YC06 site and other sites by HCLIM method.This is consistent with the demonstration in Sections 4.1 and 4.2, and also demonstrates the effectiveness of HCLIM method in compensating the height-deviation of tropospheric delay.

Conclusions
Aiming at the problem of the reduction of regional tropospheric modeling accuracy caused by height deviation between reference stations, this paper proposes a simple and effective regional tropospheric linear interpolation method (HCLIM) based on the traditional LIM modeling method, which takes height deviation compensation into account: First, prior to modeling, a prior high-precision tropospheric model is used to uniformly compensate the tropospheric delay of the baseline of the modeling alternative group to the plane where the master reference station is located.Then the traditional LIM model is used to model and interpolate the tropospheric delay at the rover station location.Finally, the interpolated tropospheric delay is compensated to the plane where the rover station is located.In this paper, data of six experimental subnets with large height deviation in a provincial CORS network in central China were collected for method verification.The modeling effects of traditional LIM and The accuracy improvement rates were 84.5%, 75.5%, 59.3%, and 26.7%, respectively.For the BDS system, The average residual RMS of LIM method and HCLIM method in the above four height angles are (0.114, 0.05, 0.028, 0.017 m) and (0.019, 0.015, 0.014, 0.013 m), and the accuracy improvement rates are (83.3%,70%, 50%, 23%), respectively 0.5%.2. For user RTK positioning, the horizontal positioning accuracy and RTK fixing rate of LIM and HCLIM methods are relatively high, and the results of the two methods do not show much difference.The average RTK fixing rate of the rover station under LIM and HCLIM methods and the average positioning error RMS in the E and N directions were (98.6%, 1.3 cm, 1.4 cm) and (99.4%, 1.0 cm, 1.1 cm), respectively.The positioning results of HCLIM showed slight improvement; However, the results of LIM method in elevation direction showed obvious systematic deviation.The average positioning error RMS in U direction of the six rover stations in LIM mode is 15.1 cm, which is obviously caused by the fact that LIM method does not consider the influence of height deviation when modeling the DD tropospheric atmosphere delay at the service end, thus introducing a large modeling residual for the user RTK solution.However, the positioning accuracy of HCLIM method in the U direction is consistent, with an average RMS of 2.6 cm, which is improved to 82.8% compared with LIM method.3.In addition, for the station with the largest height deviation (YC06), both the atmospheric model residual built on the server side and the positioning error in the U direction of the user are significantly higher than those of other stations under LIM method.However, the atmospheric modeling and positioning results of the HCLIM method at YC06 station are not significantly different from those at other stations, which demonstrates the effectiveness of the HCLIM method in compensating the tropospheric delay height deviation.
It should be noted that the HCLIM method presented in this paper only studies the tropospheric modeling under the maximum height difference between the rover station and the reference station network of about 1,000 m, and has obtained a good modeling effect.In the case of larger height deviation, the improvement of the method recommended in this paper needs further research.In addition, this paper only collected sites in a mid-latitude region in central China for research and analysis.As for low latitude and high latitude regions, the improvement of the recommended method in this paper also needs further research.

Figure 1 .
Figure 1.Flowchart of network RTK services mode (virtual reference station).
Δ   (ℎ) and  Δ   (ℎ) represents the decomposed single difference (SD) tropospheric delay between the satellites of stations B and A; h B and h A are the elevation values of the two sites respectively.Δh AB is the height difference between the two stations.Add tropospheric height-deviation compensation for the SD tropospheric delay of the two stations respectively, as shown in Equations 10 and 11: ∇ T   (ℎ) = ∇   (ℎ) − Δ∇  adj, (10) ∇ T   (ℎ) = ∇   (ℎ) − Δ∇  adj, (11) where,  Δ T   (ℎ) and  Δ T   (ℎ) denote the inter-satellite SD tropospheric delay of the two stations with height-deviation compensation, respectively. Δ∇  adj, and  Δ∇  adj, represent two compensation quantities respectively.The specific derivation process is shown in Equations 12 and 13:

Figure 2 .
Figure 2. Flowchart of virtual reference station model based on height-deviation compensation-based linear interpolation method method.

Figure 3 .
Figure 3. Site geographic distribution and elevation statistics of a CORS network (six experimental data networks selected on the left).

Figure 4 .
Figure 4. Detailed geographic distribution maps of six experimental data networks.

Figure 5 .
Figure 5. Detailed elevation statistics of six experimental data networks.

Figure 6 .
Figure 6.I95 ionospheric index of the reference station area during the data test period.

Figure 8 .
Figure 8.Comparison of two tropospheric modeling results of some satellites at YC06 station under different elevation angle types during a certain period.

Figure 9 .
Figure 9.Comparison of two tropospheric modeling results of some satellites at ES01 station under different elevation angle types during a certain period.

Figure 10 .
Figure 10.Comparison of two tropospheric modeling results of some satellites at XF12 station under different elevation angle types during a certain period.

Figure 11 .
Figure 11.Comparison of two tropospheric modeling results of some satellites at YC21 station under different elevation angle types during a certain period.

Figure 12 .
Figure 12.Comparison of two tropospheric modeling results of some satellites at XF06 station under different altitude Angle types during a certain period.

Figure 13 .
Figure 13.Comparison of two tropospheric modeling results of some satellites at SY12 station under different elevation angle types during a certain period.

Figure 14 .
Figure 14.Tropospheric modeling error RMS statistics for some satellites of all rover stations (different satellite systems and different elevation angle types) under the two modeling strategies.The percentage above the bar graph represents the accuracy improvement of the height-deviation compensation-based linear interpolation method method compared to the LIM method.

Table 2
RMS Statistics of GPS Satellite Modeling Errors of Rover Stations in Different Elevation Angles Under LIM (M1) and Height-Deviation Compensation-Based Linear Interpolation Method (M2) Modeling Strategies

Table 3
RMS Statistics of BDS Satellite Modeling Errors of Rover Stations in Different Elevation Angles Under LIM (M1) and Height-Deviation Compensation-Based Linear Interpolation Method (M2) Modeling Strategies

Table 4
Statistical Results of External Coincidence Fixed RMS and Fixed Rate of RTK Solutions Under LIM (M1) and Height-Deviation