Improving the Accuracy of Regional Ionospheric Mapping with Double-Di ﬀ erence Carrier Phase Measurement

: Regional augmentation systems for a global navigation satellite system (GNSS) provide an ionospheric map correction to the user in order to remove the ionospheric delay error. Measurements are collected from multiple reference stations to estimate the ionospheric map. During this process, the pseudorange measurement error of a reference station causes an error in the correction, which is more evident at edge areas and causes a large error for low-elevation satellites. In this study, an ionospheric modeling algorithm was developed that uses the carrier phase with the pseudorange to greatly reduce the error. The integer-resolved double-di ﬀ erence carrier phase can be obtained through ambiguity resolution method, and the measurement is directly utilized in ionospheric modeling. The performance of the developed method was tested in simulations and with real data for validation. The results of users at various locations showed that the method e ﬀ ectively improved the accuracy of the correction.


Introduction
Ionospheric delay is a major source of error for localization systems based on a global navigation satellite system (GNSS). Multi-frequency receivers can remove ionospheric delay by using the dispersive properties of the ionosphere. Recently, a low-cost, dual-frequency receiver has been developed. However, most low-cost GNSS receivers can receive only a single frequency signal and cannot remove the ionospheric error by themselves. GPS satellites broadcast the parameters of the Klobuchar ionospheric model [1]. Users can compensate for the ionospheric error by using the broadcasted parameter; however, this usually removes only 50% of the total delay. For more precise navigation, an augmentation system is needed.
For a code-based augmentation system, undifferenced ionospheric delay correction is necessary to remove the ionospheric delay error. Local area differential GNSS (LADGNSS) provides a pseudorange correction that is the sum of the signal-in-space (SIS) error and propagation error [2]. It is effective when the user is near a reference station. However, if the user is several hundred kilometers away from the reference station, the residual error after the correction increases because of the spatial decorrelation property of the error sources. To overcome this problem, wide-area differential GNSSs such as satellite-based augmentation systems (SBASs) correct each error in vector quantities [3]. Because the correction is provided in vectors, users can apply the appropriate pseudorange correction according to their position. For the ionospheric error, a 2D ionospheric map is provided to the user. Ionospheric correction is also important for carrier phase-based augmentation systems, such as real-time kinematics (RTK). Like in code-based augmentation systems, standard RTK performance is degraded as the baseline length becomes larger. Network RTK was developed to overcome this limitation. After ambiguities among reference stations are fixed, error corrections are estimated and transmitted to users. For ionospheric delay error, one of the most effective methods is providing double-difference ionospheric corrections by interpolating those obtained from reference stations [4][5][6]. A predictive regional ionosphere model was also proposed and investigated for ionospheric correction [7].
In this study, we focus on generating accurate ionospheric correction in the form of a 2D ionosphere map. To generate the correction data, multiple reference stations were distributed over the service area. Dual-frequency measurements were collected from multiple reference stations and GNSS satellites to obtain the ionospheric delay at various locations. If the ionosphere is assumed to be a thin shell, the measured slant ionospheric delay is converted to a vertical delay [8]. Then, a 2D ionospheric map can be constructed with a specific model.
Several reasons limit the accuracy of the ionospheric map correction. First, the thin-shell approach cannot perfectly emulate the real world [9]. In addition, the modeling resolution can affect the performance. When an ionospheric storm occurs, the ionosphere shows very high dynamic temporal and spatial variations, which aggravates the above limitations [10]. Another limitation is the measurement error of the reference station, which was the main focus of this study.
When the ionospheric map is estimated for an augmentation system, the ionospheric delay is usually measured according to the pseudorange. The absolute ionospheric delay of each signal path can be acquired with the pseudorange. While this measurement provides the advantage of being easy to use for ionospheric modeling, the error can be a problem. Because the pseudorange has a relatively high measurement error owing to the noise and multipath effect, an estimation error occurs in the ionospheric map. Reference stations try to mitigate this error by using a multipath-mitigation antenna or tracking loop [11]. Despite this effort, the pseudorange measurement error can be on the order of meters. There can be severe measurement errors, especially when a satellite has low elevation angle [12]. Because of the large measurement error at low elevations, the edges of the estimated ionospheric map have degraded results compared with the central area.
To overcome the correction error caused by the pseudorange measurement, an ionospheric mapping algorithm was developed in this study that uses both the pseudorange and double-difference carrier phase for the regional augmentation system. Because the carrier phase has integer ambiguity, it is difficult to use directly. The ambiguity resolution technique can be used to measure the integer-resolved double difference for ionospheric modeling [13,14]. The much smaller error of the carrier phase measurement is exploited to combine the pseudorange and double-difference carrier phase. By properly weighting each measurement, the ionospheric map can be estimated with improved modeling results at the edge areas. The performance of the proposed algorithm was verified through simulations and tests with real data.
The remainder of this paper is organized as follows. Section 2 introduces the main concepts of this study through simple diagrams. Section 3 describes the details of the ionospheric modeling procedure. Section 4 presents the results of the simulation and tests with real data. Section 5 concludes the paper.

Benefits of the Double-Difference Carrier Phase for Ionospheric Modeling
This study focused on regional GNSS ionospheric correction. If the ionosphere is assumed to be a thin shell, a vertical ionospheric delay map can be estimated over the service area. Measurements at various ionospheric pierce points (IPPs) are collected from multiple reference stations and satellites. The main idea is to directly utilize both the pseudorange and precise difference measurements in ionospheric modeling. This section briefly introduces the difference in the results when different measurements are used. Figures 1 and 2 show conceptual diagrams of the ionospheric modeling when the pseudorange and carrier phase are used. The x-axis is the IPP of each measurement. For simplicity, the ionosphere is represented in the shape of a trapezoid. The yellow and orange colors indicate the true ionosphere and modeling results, respectively.  Figure 1a shows an example of ionospheric modeling using only the pseudorange. The vertical ionospheric delay obtained from the pseudorange is denoted with blue arrows. The pseudorange has a measurement error of meters that depends on the satellite elevation angle; the measurements collected from a high-elevation satellite have a small error, while the measurements from a low-elevation satellite have a large error. When the regional ionosphere is modeled with GNSS reference stations, most measurements are obtained from high-elevation satellites at the center of the coverage. This improves the modeling accuracy at the center. At the edge of the coverage, the measurement quality is not good because the data are collected from low elevation angles. This makes the modeling results at the edges less accurate than those inside.
The carrier phase has a much smaller measurement error than the pseudorange, but it has integer ambiguity that makes it difficult to use directly. A simple way to use carrier phase is a smoothing filter [15]. Even though it has an unknown bias because of the integer ambiguity, the bias is constant over time. The rate of change in the range can be obtained precisely with the carrier phase, which can be utilized in the smoothing filter. The pseudorange measurement error can be mitigated with a smoothing filter; however, there is still a large error at the beginning of filtering. Thus, the filtered measurement also has large errors when the satellite has a low elevation angle. Because the user is located in an area with reference stations, the geometric relation between the user and each satellite is similar to those of the reference stations. The satellite elevation angle is similar for the reference stations and user. Thus, if a correction message is generated with only the pseudorange, the accuracy of the ionospheric correction for low-elevation satellites is degraded.
Unlike with the pseudorange, it is difficult to obtain the undifferenced ionospheric delay with the carrier phase. However, the integer ambiguity can be determined in the double-difference domain with ambiguity resolution methods to obtain the ambiguity-resolved double-difference ionospheric delay [13]. The characteristics of each measurement are clearly different. The pseudorange can be used to measure the undifferenced ionospheric delay, but the error is large. The double-difference carrier phase has a much smaller error but has only differenced information. In this study, the pseudorange and double-difference carrier phase were combined to take advantage of each measurement. Figure 1b shows the ionospheric modeling result with the pseudorange and double-difference carrier phase. The double-difference carrier phase-based ionospheric delay is denoted with the red arrow. Even though the pseudorange-based measurement has a large error, its characteristics are well-known. This allows each measurement to be weighted appropriately. The double-difference carrier phase can be used to precisely measure the difference in ionospheric delay at different IPPs. Based on these properties, by giving a higher weight to measurements at higher elevation angles and applying precise difference information, the accuracy of the correction at the edge areas can be improved. This alleviates the performance degradation of the correction for satellites with low elevation angles.  Figure 1a shows an example of ionospheric modeling using only the pseudorange. The vertical ionospheric delay obtained from the pseudorange is denoted with blue arrows. The pseudorange has a measurement error of meters that depends on the satellite elevation angle; the measurements collected from a high-elevation satellite have a small error, while the measurements from a low-elevation satellite have a large error. When the regional ionosphere is modeled with GNSS reference stations, most measurements are obtained from high-elevation satellites at the center of the coverage. This improves the modeling accuracy at the center. At the edge of the coverage, the measurement quality is not good because the data are collected from low elevation angles. This makes the modeling results at the edges less accurate than those inside.
The carrier phase has a much smaller measurement error than the pseudorange, but it has integer ambiguity that makes it difficult to use directly. A simple way to use carrier phase is a smoothing filter [15]. Even though it has an unknown bias because of the integer ambiguity, the bias is constant over time. The rate of change in the range can be obtained precisely with the carrier phase, which can be utilized in the smoothing filter. The pseudorange measurement error can be mitigated with a smoothing filter; however, there is still a large error at the beginning of filtering. Thus, the filtered measurement also has large errors when the satellite has a low elevation angle. Because the user is located in an area with reference stations, the geometric relation between the user and each satellite is similar to those of the reference stations. The satellite elevation angle is similar for the reference stations and user. Thus, if a correction message is generated with only the pseudorange, the accuracy of the ionospheric correction for low-elevation satellites is degraded.
Unlike with the pseudorange, it is difficult to obtain the undifferenced ionospheric delay with the carrier phase. However, the integer ambiguity can be determined in the double-difference domain with ambiguity resolution methods to obtain the ambiguity-resolved double-difference ionospheric delay [13]. The characteristics of each measurement are clearly different. The pseudorange can be used to measure the undifferenced ionospheric delay, but the error is large. The double-difference carrier phase has a much smaller error but has only differenced information. In this study, the pseudorange and double-difference carrier phase were combined to take advantage of each measurement. Figure 1b shows the ionospheric modeling result with the pseudorange and double-difference carrier phase. The double-difference carrier phase-based ionospheric delay is denoted with the red arrow. Even though the pseudorange-based measurement has a large error, its characteristics are well-known. This allows each measurement to be weighted appropriately. The double-difference carrier phase can be used to precisely measure the difference in ionospheric delay at different IPPs. Based on these properties, by giving a higher weight to measurements at higher elevation angles and applying precise difference information, the accuracy of the correction at the edge areas can be improved. This alleviates the performance degradation of the correction for satellites with low elevation angles. Figure 2 shows the overall procedure of the proposed method. Raw GNSS measurements are collected from multiple reference stations. During the preprocessing stage, cycle slip detection is carried out to prevent an unexpected jump in the carrier phase measurement. The cycle slip can be detected by monitoring the acceleration of the ionospheric delay because it is less than 10 −4 m/s 2 [16][17][18]. Unexpected outliers are detected by using the second derivative of the geometry-free combination of the carrier phase as a monitoring value. A smoothing filter and differential code bias (DCB) calibration are needed to mitigate the pseudorange measurement error. A divergence-free hatch filter is implemented to address noise and the multipath effect [16,19]. It is a type of smoothing filter that utilizes dual-frequency carrier phase measurement to remove the divergence between the pseudorange and the carrier. When the ionospheric delay is measured with the pseudorange, the DCB can be a large error source [20]. The hardware delay of the GNSS satellite and receiver have different values depending on the signal frequency. The difference in hardware delays of the pseudorange measurement obtained from two different frequencies is called the DCB; it can be several tens of nanoseconds, and should be removed. The DCBs can be estimated using a geometry-free combination of the pseudorange measurement. The geometry combination of the pseudorange comprises the ionospheric delay, satellite DCB, and receiver DCB. Because the variation of DCBs is usually very small [21], it can be assumed to be constant for a day. With a specific model for vertical ionosphere, vertical ionosphere and DCBs can be estimated simultaneously. In this study, a spherical harmonics model was adopted for the ionospheric model [21,22]. Using the measurements collected in one day, DCBs were estimated with the batch process.

Ionosphere Modeling Algorithm
The double-difference carrier phase is measured by pairing two reference stations. Among the stations, one station was selected as the master and the others were selected as auxiliaries. By pairing the master and each auxiliary, integer ambiguities were estimated. Because the surveyed position of reference stations was known, only integer ambiguities were estimated. To prevent the distance between stations from becoming too large, the reference station located at the center was selected as the master. Wide lane ambiguity was resolved by the Melbourne-Wubbena combination. Then, both ambiguities of L1 and L2 frequency were estimated using the iono-free combination [23]. Based on these, the double-difference ionospheric delay of each pair could be obtained.

Measurement
In this study, the L1 and L2 frequencies of the Global Positioning System were measured. Equations (1)-(4) represent the raw pseudorange and carrier phase measurement [24]: The pseudorange-based measurement and ambiguity-resolved carrier phase measurement were input into the ionospheric model. The final goal of this process is to generate a vertical ionospheric delay map. Each measurement was represented by a vertical ionosphere model, which is a spherical harmonics model. Because the pseudorange measurement has a slant delay, it can be represented by multiplying the mapping function to the vertical model. The double-difference carrier phase is a linear combination of four different slant delays. It can be represented by converting the vertical model to a slant model and taking the double-difference combination of the slant model. The vertical ionospheric map is generated using both the pseudorange and double-difference carrier phase measurements from all reference stations. Sections 3.1 and 3.2 describe the basic GPS measurement and ionospheric model, respectively. Section 3.3 presents the observation equation for estimating the modeling parameters.

Measurement
In this study, the L1 and L2 frequencies of the Global Positioning System were measured. Equations (1)-(4) represent the raw pseudorange and carrier phase measurement [24]: where ρ is the pseudorange, φ is the carrier phase, d is the geometric distance, I is the ionospheric delay for the L1 frequency, γ is the ratio between the squares of the signal frequencies, T is the tropospheric delay, b is the satellite clock bias, B is the receiver clock bias, b Tx is the satellite hardware delay, b Rx is the receiver hardware delay, ε is the noise and multipath error, N is the integer ambiguity, and λ is the wavelength. Subscripts 1 and 2 indicate the L1 and L2 frequencies, respectively. Because the ionosphere has a dispersive property, the L1 and L2 frequencies have different ionospheric delays. Based on this property, the ionospheric delay can be obtained by geometry-free combination. The ionospheric delay of the L1 frequency is calculated as follows: where I ρ,DCB is the pseudorange-based ionospheric delay measurement with DCBs, I is the true ionospheric delay, DCB Tx ρ1/ρ2 and DCB Rx ρ1/ρ2 are the satellite DCB and receiver DCB, respectively, and ε I ρ is the measurement noise. After the DCBs are estimated and removed, the measurement used for ionospheric modeling is represented as follows: where ε b Tx and ε b Rx are the estimation errors of the satellite and receiver DCBs, respectively. The carrier phase-based ionospheric measurement is represented as follows: This is similar to Equation (5) except for the integer ambiguity term. The double difference of this measurement is expressed by Equation (8). During the double differencing process, both the satellite and receiver hardware delays are removed.
where ∇∆I φ,N is the double-difference ionospheric delay measurement with integer ambiguity. If ambiguities are fixed correctly, the ambiguity-resolved measurement is expressed as follows: where ∇∆I φ is the ambiguity-resolved double-difference ionospheric delay measurement and ε ∇∆I is the measurement error.

Ionospheric Model
The thin-shell assumption was adopted for the ionospheric modeling. The slant ionospheric delay can be converted to the vertical ionospheric delay by using a simple function called the obliquity factor [8]: where I S is the slant ionospheric delay, I V is the vertical ionospheric delay, and Q is the obliquity factor. Various methods have been developed for ionospheric modeling. Grid-based interpolation algorithms are a well-known example, which include the inverse distance weight (IDW) method, planar fit method, and kriging method [25][26][27][28]. These methods have different ways of calculating weights for each measurement. The IDW calculates simple weight depending on the distance. Kriging uses the statistical model to determine the weight of each measurement. Although they have differences, the basic principle is the same: they estimate the ionospheric delay at each grid point by interpolating the measurements within a certain distance limit from the grid point [26,27]. Including the double-difference measurement in grid-based algorithms is difficult because the double-difference measurement comprises four different measurements at four different locations, and there are many cases where all four different measurements are not within a certain distance limit.
Spherical harmonics was selected as the modeling method in this study [28,29]. This model represents spatial data by the linear combination of the basis functions. Distributed measurements are used to estimate the whole region at once. Spherical harmonics is suitable for double-difference measurements.
First, the spherical harmonics for the vertical ionospheric delay is expressed as follows: where I V is the vertical ionospheric delay, θ is the latitude of the IPP, λ is the longitude of the IPP, P nm is the associated Legendre function, C nm and S nm are spherical harmonics coefficients, and N is the degree. Measurements collected from various locations are used to estimate the coefficients. With the spherical harmonics model, the pseudorange-based slant delay measurement is represented as follows: The double-difference measurement is a linear combination of four different slant delays:

Observation Equation
If multiple reference stations and satellites are assumed, an observation matrix can be constructed that uses all pseudorange and carrier phase measurements. Equations (14)-(16) are the measurement vectors: where z is the measurement vector for ionospheric modeling, z ρ is the pseudorange-based measurement vector, and z ∆∇φ is the double-difference carrier-based measurement vector. The superscript indicates the satellite, and the subscript indicates the reference station. M is the master station, and i is an auxiliary station. The master station is paired with each auxiliary station for the double-difference measurement; r indicates the reference satellite. Because two types of measurement with significantly different error levels are used for the estimation, it is important to construct a covariance matrix to properly weight each measurement according to the error level. The covariance matrix of the measurements is expressed as follows: where V is the covariance matrix of all measurements, V ρ is the covariance matrix of the pseudorange-based measurement, and V ∆∇φ is the covariance matrix of the double-difference carrier phase-based measurement. If the pseudorange and carrier phase are assumed to have no correlation, the off-diagonal terms of V are zeros. The detailed elements of V ρ are expressed below: where σ l,l 2 is the diagonal element of V ρ , σ l,k 2 is an off-diagonal element of V ρ , l is a row number, k is a column number, σ ρ1 2 and σ ρ2 2 are the variances of the L1 and L2 pseudorange measurement errors, respectively, and σ b Tx 2 and σ b Rx 2 are the variances of the estimation errors for the satellite and receiver DCBs, respectively. The noise and multipath effect of the L1 and L2 frequency pseudorange are not very different, and an empirical code noise and multipath (CNMP) model has been developed that depends on the smoothing time [30,31]. With this empirical model, σ l,l 2 can be expressed as follows: where σ CNMP 2 is the empirical CNMP variance. σ CNMP was determined so that all the empirical errors can be bounded. Figure 3 shows σ CNMP against the smoothing time. Initially, it was approximately 3 m, and continuously decreased. After 12,000 s, it plateaued at a value of 0.12 m.
Remote Sens. 2019, 11, x FOR PEER REVIEW 9 of 21 The elements of the covariance matrix of the double-difference measurement are expressed in Equations (21)- (24). Because the double-difference measurements were made with one master station, all measurements were related with the reference satellite of the master station. In addition, there are many correlations based on the terms shared by the measurements.  I  I  I  I   I  I  I  I   I  I  I  I   I  I  I  I is the variance of the carrier phase-based ionospheric delay collected from the i -th The elements of the covariance matrix of the double-difference measurement are expressed in Equations (21)- (24). Because the double-difference measurements were made with one master station, all measurements were related with the reference satellite of the master station. In addition, there are many correlations based on the terms shared by the measurements.
where σ (I φ ) j i 2 is the variance of the carrier phase-based ionospheric delay collected from the i-th station and j-th satellite, σ i ∆ r M I φ 2 and σj ∇ r M I φ 2 are the variances of the single-difference measurements for the receiver and satellite, respectively, and σj i ∆∇ r M I φ 2 is the variance of the double-difference measurement from the i-th auxiliary station and j-th satellite. According to Equation (9), σ (I φ ) j i 2 can be expressed as follows: where σ 2 are the variances of the L1 and L2 carrier phase measurements, respectively. The carrier phase had a much smaller measurement error than the pseudorange with a noise error of millimeters and multipath error of centimeters. The maximum size of the multipath error can reach a quarter of its wavelength [32]. To conservatively bound the error, the 3-sigma of the L1 and L2 carrier phase measurement errors was set to a quarter of each wavelength. That is, one sigma was set to 1.6 cm for the L1 signal and 2 cm for the L2 signal.
The collected measurements were used to estimate the spherical harmonics coefficients. Equations (26) and (27) show the coefficient vector and observation equation, respectively: where x is the coefficient vector, H is the full observation matrix, H ρ is the observation matrix for the pseudorange-based ionospheric delay measurement, and H ∆∇ is the observation matrix for the double-difference carrier phase-based ionospheric delay measurement. The detailed elements of the observation matrix for the pseudorange and double-difference carrier phase measurement are given below: With the above measurement equation and covariance matrix, the weighted least-square estimation equation is expressed as follows: wherex is the estimated coefficient vector. If a user has coefficients, they can use them directly to reconstruct the estimated ionospheric delay map. The estimated coefficients were assumed to be delivered to the user as correction data. Then, the user can calculate the ionospheric correction by substituting the locations of IPPs into Equation (12).

Test Results
Tests were carried out with simulation and real data to verify the improvement with the proposed method. Reference stations and users were selected for an augmentation system in South Korea. An ionospheric correction was generated with five reference stations and the accuracy was analyzed. Only GPS satellites were used in this test. Methods using the pseudorange only were also implemented for comparison with the proposed method. The main comparison target was a spherical harmonics model that used only the pseudorange. Two grid-based algorithms were also implemented: IDW and Kriging. Under the assumption of a real-time process, correction data were estimated every epoch from measurements taken at the corresponding time. Figure 4 shows the reference stations and users for the simulation test. The reference stations were set equal to those used for the test with actual data. To improve the reliability of results, multiple users across all service areas were selected. The users were uniformly distributed over the reference station network, and their numbers are shown in the figure. Simulation data were generated using IONosphere map EXchange (IONEX) format data [33]. On 1 September 2014, the simulation data were generated for 12 h in the day time (9 AM to 9 PM). The satellite positions during the test period were used to calculate the IPPs of the reference stations and users, and ionospheric delays were generated by interpolation of the IONEX data.  Figure 4 shows the reference stations and users for the simulation test. The reference stations were set equal to those used for the test with actual data. To improve the reliability of results, multiple users across all service areas were selected. The users were uniformly distributed over the reference station network, and their numbers are shown in the figure. Simulation data were generated using IONosphere map EXchange (IONEX) format data [33]. On 1 September 2014, the simulation data were generated for 12 h in the day time (9 AM to 9 PM). The satellite positions during the test period were used to calculate the IPPs of the reference stations and users, and ionospheric delays were generated by interpolation of the IONEX data. The vertical ionospheric map results were analyzed in order to demonstrate the effect of the double-difference measurement on the ionospheric modeling. Figure 5 shows the true ionospheric map, and Figures 6-8 show the modeling results. The values are the amount of delay appearing at the L1 frequency. The IPPs of all measurement are marked with dark circles. For a given ionospheric map and measurements, spherical harmonics-based modeling was carried out with only the pseudorange, only the double-difference carrier phase, and both the pseudorange and double-difference carrier phase.  Figure 6 shows the modeling results with the pseudorange only and the error. As expected, the accuracy inside the modeling area was good, and significant degradation occurred at the edge. The vertical ionospheric map results were analyzed in order to demonstrate the effect of the double-difference measurement on the ionospheric modeling. Figure 5 shows the true ionospheric map, and  Figure 4 shows the reference stations and users for the simulation test. The reference stations were set equal to those used for the test with actual data. To improve the reliability of results, multiple users across all service areas were selected. The users were uniformly distributed over the reference station network, and their numbers are shown in the figure. Simulation data were generated using IONosphere map EXchange (IONEX) format data [33]. On 1 September 2014, the simulation data were generated for 12 h in the day time (9 AM to 9 PM). The satellite positions during the test period were used to calculate the IPPs of the reference stations and users, and ionospheric delays were generated by interpolation of the IONEX data. The vertical ionospheric map results were analyzed in order to demonstrate the effect of the double-difference measurement on the ionospheric modeling. Figure 5 shows the true ionospheric map, and Figures 6-8 show the modeling results. The values are the amount of delay appearing at the L1 frequency. The IPPs of all measurement are marked with dark circles. For a given ionospheric map and measurements, spherical harmonics-based modeling was carried out with only the pseudorange, only the double-difference carrier phase, and both the pseudorange and double-difference carrier phase.  Figure 6 shows the modeling results with the pseudorange only and the error. As expected, the accuracy inside the modeling area was good, and significant degradation occurred at the edge.  Figure 6 shows the modeling results with the pseudorange only and the error. As expected, the accuracy inside the modeling area was good, and significant degradation occurred at the edge. Figure 7 shows the modeling results with the double difference only and the error. Unlike the pseudorange-only case, the errors inside and at the edges were not significantly different. The error was smaller at the edges; however, the overall results were biased, and the error inside was greater than that of the pseudorange-only case. Because this case used only differenced data, bias could occur. Figure 8 shows the modeling results using both measurements and the error. This case showed better results for all areas compared with the previous cases. The accuracy in the center was similar to that of the pseudorange-only case, and the error at the edge areas was mitigated.

Simulation Test
Remote Sens. 2019, 11, x FOR PEER REVIEW 12 of 21 than that of the pseudorange-only case. Because this case used only differenced data, bias could occur. Figure 8 shows the modeling results using both measurements and the error. This case showed better results for all areas compared with the previous cases. The accuracy in the center was similar to that of the pseudorange-only case, and the error at the edge areas was mitigated.
(a) (b)  than that of the pseudorange-only case. Because this case used only differenced data, bias could occur. Figure 8 shows the modeling results using both measurements and the error. This case showed better results for all areas compared with the previous cases. The accuracy in the center was similar to that of the pseudorange-only case, and the error at the edge areas was mitigated.
(a) (b)  than that of the pseudorange-only case. Because this case used only differenced data, bias could occur. Figure 8 shows the modeling results using both measurements and the error. This case showed better results for all areas compared with the previous cases. The accuracy in the center was similar to that of the pseudorange-only case, and the error at the edge areas was mitigated.
(a) (b)  The results in the user range domain were analyzed. Four different algorithms were used in this test: (1) IDW, (2) kriging, (3) spherical harmonics with the pseudorange only, and (4) spherical harmonics with both the pseudorange and double-difference carrier. The user ionospheric correction was calculated for the modeling results of each algorithm. The residual error in the user range domain was computed by comparing the true ionospheric delay with the computed user correction. For each user, the residual errors of the four different algorithms were calculated, and their RMS values are presented in Table 1. All user data were collected to calculate the overall results. Algorithm 4 (i.e., the proposed method) showed a smaller error than the other methods for all users. The RMS error of all user data for Algorithms 1-4 were 0.37, 0.41, 0.23, and 0.16, respectively. Figure 9 shows the error distribution of each algorithm. The errors were assumed to have a Gaussian distribution; the mean and standard deviation were calculated, and the probability distribution function was plotted. The green, blue, cyan, and red colors indicate the results of Algorithms 1-4, respectively.  Algorithm 4 (i.e., the proposed method) showed a smaller error than the other methods for all users. The RMS error of all user data for Algorithms 1-4 were 0.37, 0.41, 0.23, and 0.16, respectively. Figure 9 shows the error distribution of each algorithm. The errors were assumed to have a Gaussian distribution; the mean and standard deviation were calculated, and the probability distribution function was plotted. The green, blue, cyan, and red colors indicate the results of Algorithms 1-4, respectively. User 5 (at the center of the region) showed the smallest difference between Algorithms 3 and 4. Because the IPP distribution of User 5 was located inwards compared with those of the other users, the pseudorange error had a relatively small effect. To clarify where improvement occurred, users were grouped, and the residual error of each group was separated by the elevation angle. Figure 10 shows the RMS error of each elevation angle interval. Figure 10a shows the results of User 5 at the center, and Figure 10b shows the results of the other users at the edges of the service area. The errors of all algorithms tended to increase as the elevation angle decreased. For User 5, Algorithm 4 provided similar results as Algorithm 3. For Users 1-4 and 6-9, the increased error at low elevation angles was more evident for Algorithms 1-3. While the other algorithms performed worse at the edges, Algorithm 4 performed relatively consistently. These results imply that the proposed algorithm can improve the accuracy of low-elevation satellites for users at edge areas. User 5 (at the center of the region) showed the smallest difference between Algorithms 3 and 4. Because the IPP distribution of User 5 was located inwards compared with those of the other users, the pseudorange error had a relatively small effect. To clarify where improvement occurred, users were grouped, and the residual error of each group was separated by the elevation angle. Figure 10 shows the RMS error of each elevation angle interval. Figure 10a shows the results of User 5 at the center, and Figure 10b shows the results of the other users at the edges of the service area. The errors of all algorithms tended to increase as the elevation angle decreased. For User 5, Algorithm 4 provided similar results as Algorithm 3. For Users 1-4 and 6-9, the increased error at low elevation angles was more evident for Algorithms 1-3. While the other algorithms performed worse at the edges, Algorithm 4 performed relatively consistently. These results imply that the proposed algorithm can improve the accuracy of low-elevation satellites for users at edge areas.

Real Data Test
While IONEX data is useful for checking the performance in both the vertical and slant domains, it has limitations. The resolution is not sufficient for expressing both spatial and temporal details. To show the effectiveness of the proposed method in the real world, a test was carried out with an actual dataset, with five reference stations and five users. Receiver Independent Exchange Format (RINEX) data were obtained from the archives of the National Geographic Information Institute in South Korea [34]. Figure 11 shows the reference stations and users. The users were numbered as shown in the figure. All the selected reference stations were equipped with a Trimble NETR9 receiver and geodetic antenna because the users also needed high-quality measurements for accurate analysis. The position, receiver type, and antenna type of the reference stations and users are listed in Table 2. For auxiliary stations, the baseline length is also provided.
Considering the ionospheric activity, two days were selected. First, data were collected on 1 September 2014, under quiet conditions. Next, data were collected on 12 September 2014, under more active conditions. The daily average Kp index was 2.5 for the first data and 4 for the second data. For each day, ionospheric modeling was conducted with the measurements from the reference stations. The results were applied to the users, and the residual errors were analyzed. Figure 11. Reference station and user (real data). Table 2. Reference station and user information. Baseline length is listed for auxiliary stations.

Real Data Test
While IONEX data is useful for checking the performance in both the vertical and slant domains, it has limitations. The resolution is not sufficient for expressing both spatial and temporal details. To show the effectiveness of the proposed method in the real world, a test was carried out with an actual dataset, with five reference stations and five users. Receiver Independent Exchange Format (RINEX) data were obtained from the archives of the National Geographic Information Institute in South Korea [34]. Figure 11 shows the reference stations and users. The users were numbered as shown in the figure. All the selected reference stations were equipped with a Trimble NETR9 receiver and geodetic antenna because the users also needed high-quality measurements for accurate analysis. The position, receiver type, and antenna type of the reference stations and users are listed in Table 2. For auxiliary stations, the baseline length is also provided.

Real Data Test
While IONEX data is useful for checking the performance in both the vertical and slant domains, it has limitations. The resolution is not sufficient for expressing both spatial and temporal details. To show the effectiveness of the proposed method in the real world, a test was carried out with an actual dataset, with five reference stations and five users. Receiver Independent Exchange Format (RINEX) data were obtained from the archives of the National Geographic Information Institute in South Korea [34]. Figure 11 shows the reference stations and users. The users were numbered as shown in the figure. All the selected reference stations were equipped with a Trimble NETR9 receiver and geodetic antenna because the users also needed high-quality measurements for accurate analysis. The position, receiver type, and antenna type of the reference stations and users are listed in Table 2. For auxiliary stations, the baseline length is also provided.
Considering the ionospheric activity, two days were selected. First, data were collected on 1 September 2014, under quiet conditions. Next, data were collected on 12 September 2014, under more active conditions. The daily average Kp index was 2.5 for the first data and 4 for the second data. For each day, ionospheric modeling was conducted with the measurements from the reference stations. The results were applied to the users, and the residual errors were analyzed.   Error (m) Figure 11. Reference station and user (real data). Considering the ionospheric activity, two days were selected. First, data were collected on 1 September 2014, under quiet conditions. Next, data were collected on 12 September 2014, under more active conditions. The daily average Kp index was 2.5 for the first data and 4 for the second data. For each day, ionospheric modeling was conducted with the measurements from the reference stations. The results were applied to the users, and the residual errors were analyzed.
This study was focused on the effect of the double-difference carrier phase on ionospheric modeling. Thus, the integer ambiguity needed to be correctly resolved. To avoid unexpected results caused by wrong integers, all integers were estimated in batches using all measurements collected throughout the test. The estimated integers were validated against the double-difference residuals. DCBs were also estimated in batches, and the errors were neglected.
Unlike the simulation, showing the error in the vertical ionospheric map for the real environment was difficult because of the challenge in making a true ionospheric map. Thus, only a range domain analysis was carried out. The precise ionospheric delay of each user was required to analyze the results in the range domain. The leveled carrier phase measurement was utilized as the ground truth data when the residual errors were calculated [35,36]. For a continuous arc with no cycle slip, the integer ambiguity is constant. In addition, the hardware biases are nearly constant over time. Therefore, there is a constant bias between the pseudorange-based measurement and carrier phase measurement, which is expressed in Equations (6) and (7). These equations can be used to compute a constant bias by averaging the difference between the two measurements. Then, the leveled carrier phase measurement can be expressed as follows: where I leveled is the leveled carrier phase ionospheric delay measurement and ε leveled is the error necessary for averaging over the entire arc and can only be calculated by a post process. Because the error is much smaller than that of the pseudorange measurement, it can be used as the ground truth. First, the test was conducted under quiet conditions. Similar to the previous simulation results, Table 3 presents the RMS of the residual error in the range domain for the real data test. Figure 12 presents the error distribution of each algorithm.
The proposed method (i.e., Algorithm 4) showed the best results for all users. The RMS errors of all user data for Algorithms 1-4 were 0.57, 0.58, 0.32, and 0.24 m, respectively. Users 1 and 5 were located relatively inward compared with the other users. With Algorithm 3, the RMS error of Users 1 and 5 was 0.20 m. On the other hand, the errors increased for users that were located outward. The RMS error of Users 2-4 was 0.37 m. With Algorithm 4, the RMS error of Users 1 and 5 was 0.19 m, which is not significantly different from that of Algorithm 3. However, the algorithms showed differences for Users 2-4. The proposed algorithm mitigated the increased error, and the RMS error of Users 2-4 was 0.27 m. Table 3. RMS of the residual error (real data, 1 September 2014). Users are grouped into two groups depending on their position. Users 1 and 5 were located relatively inward, and Users 2-4 were located outward. The proposed method (i.e., Algorithm 4) showed the best results for all users. The RMS errors of all user data for Algorithms 1-4 were 0.57, 0.58, 0.32, and 0.24 m, respectively. Users 1 and 5 were located relatively inward compared with the other users. With Algorithm 3, the RMS error of Users 1 and 5 was 0.20 m. On the other hand, the errors increased for users that were located outward. The RMS error of Users 2-4 was 0.37 m. With Algorithm 4, the RMS error of Users 1 and 5 was 0.19 m, which is not significantly different from that of Algorithm 3. However, the algorithms showed differences for Users 2-4. The proposed algorithm mitigated the increased error, and the RMS error of Users 2-4 was 0.27 m. Figure 13 shows the relation between the error and elevation angle. Similar to the simulation results, the error tended to increase as the elevation angle decreased. While Algorithms 3 and 4 had similar results for Users 1 and 5, there were clear differences for Users 2-4. Above 75°, the RMS errors were 0.19 m for Algorithm 3 and 0.15 m for Algorithm 4. Around 10°, the error of Algorithm 3 increased drastically, while the degradation was mitigated with Algorithm 4 with RMS errors of 0.79 and 0.45 m, respectively.  The second test was carried out under active conditions. The test results are represented in a manner similar to the previous results. Table 4 lists the RMS of the residual error in the range domain for the real data test. Figure 14 presents the error distribution of each algorithm. The second test was carried out under active conditions. The test results are represented in a manner similar to the previous results. Table 4 lists the RMS of the residual error in the range domain for the real data test. Figure 14 presents the error distribution of each algorithm.   Because the ionosphere was more active, the overall error increased compared with the first test. Similar to the previous results, Algorithm 4 exhibited the best results. The RMS errors of all user data for Algorithms 1-4 were 1.07, 1.01, 0.41, and 0.31 m, respectively. For Users 1 and 5, located relatively inward, Algorithm 3 and 4 showed similar results. However, Algorithm 4 showed improved results for Users 2-4, located relatively outward. The RMS errors of Algorithm 3 and 4 were 0.49 m and 0.37 m, respectively. Figure 15 shows the relation between the error and the elevation angle. Compared with the quiet day, the error further increased as the elevation angle decreased for all algorithms. However, Algorithm 4 still showed a decrease in the error for satellites with a low elevation angle. In particular, it was more evident for Uses 2-4. Around 10°, the RMS error was 1.1 m for Algorithm 3 and 0.8 m for Algorithm 4. These results clearly indicate that the proposed algorithm improved the ionospheric correction accuracy for satellites with a low elevation angle when users were at the edges of the service area. Because the ionosphere was more active, the overall error increased compared with the first test. Similar to the previous results, Algorithm 4 exhibited the best results. The RMS errors of all user data for Algorithms 1-4 were 1.07, 1.01, 0.41, and 0.31 m, respectively. For Users 1 and 5, located relatively inward, Algorithm 3 and 4 showed similar results. However, Algorithm 4 showed improved results for Users 2-4, located relatively outward. The RMS errors of Algorithm 3 and 4 were 0.49 m and 0.37 m, respectively. Figure 15 shows the relation between the error and the elevation angle. Compared with the quiet day, the error further increased as the elevation angle decreased for all algorithms. However, Algorithm 4 still showed a decrease in the error for satellites with a low elevation angle. In particular, it was more evident for Uses 2-4. Around 10 • , the RMS error was 1.1 m for Algorithm 3 and 0.8 m for Algorithm 4. These results clearly indicate that the proposed algorithm improved the ionospheric correction accuracy for satellites with a low elevation angle when users were at the edges of the service area.

Discussion
The main area of interest of this study was the edge region of the ionospheric correction map. A relatively small number of measurements are distributed across the edge region, such that the level of performance is degraded [37,38]. The small number of measurements has two adverse effects on the estimation performance. First, it increases the effect of any measurement error. When a large number of measurements are acquired, the error can be averaged, such that the influence of any measurement error will be smaller. However, when only a small number of measurements are acquired, this averaging effect is less likely to occur, such that measurement errors may have a large effect on the estimation results. In addition, the majority of the measurements made in the edge area are obtained from low-elevation satellites, which incur a relatively large error. These large errors associated with low-elevation satellites can have a considerable impact on the estimation results under this vulnerable condition. Second, an undersampling problem occurs [25,39]. To accurately estimate the shape of the ionosphere, it is desirable that measurements be collected from various locations. Measurements collected at the edge region alone are not sufficient, in that a large degree of error may arise in the correction. This problem may be more serious in situations in which the ionosphere is highly active, such as in the event of ionospheric storms. The proposed algorithm can lessen the performance degradation that results from these measurement errors. However, the problem of undersampling cannot be alleviated with the proposed method, and additional research is needed.
There are several issues that need to be addressed to further realize this research. The first is to improve the reliability and availability of ambiguity resolution for long baselines. Because there is a considerable distance between the reference stations, the spatial decorrelation of any error affects the performance of real-time ambiguity resolution [40]. If the ambiguity is incorrectly determined, the proposed algorithm can give rise to a large estimation error because it places a considerable weighting on the double-differenced carrier phase. To focus on the effectiveness of the proposed ionospheric modeling algorithms, we used pre-estimated results with a batch process preventing any unexpected errors caused by the ambiguity being incorrect. To enable the application of the new algorithm to a real-time augmentation system, this problem must be overcome. Second, a comparison with existing algorithms considering message capacity is needed. The results presented herein were obtained by directly applying the estimated ionospheric map to the users. Depending on the application, the data rate may be limited. For example, the data rate of SBAS is only 250 bps [41]. As such, there is a need to apply a correction message according to the message standard and to analyze the result that appears when it is transmitted to the user at the actual transmission rate.

Discussion
The main area of interest of this study was the edge region of the ionospheric correction map. A relatively small number of measurements are distributed across the edge region, such that the level of performance is degraded [37,38]. The small number of measurements has two adverse effects on the estimation performance. First, it increases the effect of any measurement error. When a large number of measurements are acquired, the error can be averaged, such that the influence of any measurement error will be smaller. However, when only a small number of measurements are acquired, this averaging effect is less likely to occur, such that measurement errors may have a large effect on the estimation results. In addition, the majority of the measurements made in the edge area are obtained from low-elevation satellites, which incur a relatively large error. These large errors associated with low-elevation satellites can have a considerable impact on the estimation results under this vulnerable condition. Second, an undersampling problem occurs [25,39]. To accurately estimate the shape of the ionosphere, it is desirable that measurements be collected from various locations. Measurements collected at the edge region alone are not sufficient, in that a large degree of error may arise in the correction. This problem may be more serious in situations in which the ionosphere is highly active, such as in the event of ionospheric storms. The proposed algorithm can lessen the performance degradation that results from these measurement errors. However, the problem of undersampling cannot be alleviated with the proposed method, and additional research is needed.
There are several issues that need to be addressed to further realize this research. The first is to improve the reliability and availability of ambiguity resolution for long baselines. Because there is a considerable distance between the reference stations, the spatial decorrelation of any error affects the performance of real-time ambiguity resolution [40]. If the ambiguity is incorrectly determined, the proposed algorithm can give rise to a large estimation error because it places a considerable weighting on the double-differenced carrier phase. To focus on the effectiveness of the proposed ionospheric modeling algorithms, we used pre-estimated results with a batch process preventing any unexpected errors caused by the ambiguity being incorrect. To enable the application of the new algorithm to a real-time augmentation system, this problem must be overcome. Second, a comparison with existing algorithms considering message capacity is needed. The results presented herein were obtained by directly applying the estimated ionospheric map to the users. Depending on the application, the data rate may be limited. For example, the data rate of SBAS is only 250 bps [41]. As such, there is a need to apply a correction message according to the message standard and to analyze the result that appears when it is transmitted to the user at the actual transmission rate. Third, it is necessary to find the optimal weighting for each measurement. In this study, the variance of each measurement was assigned conservatively. It is expected that accuracy could be further improved if a measurement model that reflects the actual error level of each measurement was applied. Finally, in addition to accuracy, integrity performance is also an important factor affecting the augmentation system. The integrity problem of the algorithm should be investigated by utilizing as much long-term data as possible.

Conclusions
An ionospheric modeling algorithm is proposed that uses both the pseudorange and double-difference carrier phase measurement. When the ionospheric correction is generated with only the pseudorange, the performance is degraded for users at the edges of the service area, particularly for satellites with a low elevation angle. The proposed algorithm mitigated this degradation and provided relatively consistent results. The improvement was verified through tests using simulation and real data. The proposed algorithm showed the smallest error in all tests. Especially, it was more effective for satellites with a low elevation angle. In the real data test, the error decreased by about 50% for a quiet day and 30% for an active day.
While the proposed algorithm could improve the ionospheric correction accuracy, there are several issues with regard to its application in an augmentation system. The reference stations were located several hundred kilometers apart, and the wrong integer could be estimated. Because a wrong integer can generate a large error, robust real-time estimation and validation of the integer ambiguity should be implemented. Further analysis should also be carried out regarding the message bandwidth.