Research on Accelerating Single-Frequency Precise Point Positioning Convergence with Atmospheric Constraint

Abstract

An increasing number of researchers have conducted in-depth research on the advantages of low-cost single-frequency (SF) receivers, which can effectively use ionospheric information when compared to dual-frequency ionospheric-free combination. However, SF observations are bound to increase the unknown parameters and prolong the convergence time. It is desirable if the convergence time can be reduced by external information constraints, for example atmospheric constraints, which include ionosphere- or troposphere constraints. In this study, ionospheric delay constraints, tropospheric delay constraints, and their dual constraints were considered. Additionally, a total of 18,720 test experiments were performed. First, the nearest-neighbor extrapolation (NENE), bilinear- (BILI), bicubic- (BICU), and Junkins weighted-interpolation (JUNK) method of Global Ionospheric Map (GIM) grid products were analyzed. The statistically verified BILI in the percentage of convergence time, average convergence time, and computation time consumption of them shows a good advantage. Next, the influences of global troposphere- and ionosphere-constrained on the convergence time of SF Precise Point Positioning (PPP) were analyzed. It is verified that the ionosphere-constrained (TIC2) has significant influence on the convergence time in the horizontal and vertical components, while the troposphere-constrained (TIC1) has better effect on the convergence time in the vertical components within some thresholds. Of course, the dual constraint (TIC3) has the shortest average convergence time, which is at least 46.5% shorter in static mode and 5.4% in kinematic mode than standard SF PPP (TIC0).

1. Introduction

The concept and technology of Precise Point Positioning (PPP) were first proposed and implemented for the Global Positioning System (GPS) by the American Jet Propulsion Laboratory (JPL) in the late 1990s [1]. PPP has attracted significant interest over the intervening years due to its high accuracy without needing a specific reference station, providing correctional information, simple operations, and cost effectiveness due to reductions in labor and equipment costs. Therefore, it has been extensively used in the areas, for instance, determines the precise orbit, surface ice flow speed, as well as positioning, navigation and timing (PNT) applications [2,3,4]. The traditional PPP models are based on Ionosphere-Free (IF) combinations with dual-frequency raw phase and code observations for the removal of the first-order effect of ionospheric refraction [5]. However, the second- and third-order ionospheric effects still exist and they may cause measurements errors of sub-centimeters in GPS [6,7]. Therefore, traditional PPP technology based on IF combinations cannot obtain ionospheric information, and the ionospheric error is not completely eliminated. In the last few years, PPP used raw observations has received an increasing amount of research attention because of its several advantages as compared to traditional IF PPP, particularly the development of new multi-frequency GNSS [8].

The ionosphere refers to the atmospheric space starting at 60 km above the ground and extending to the magnetosphere. The Global Navigation Satellite System (GNSS) signal will generate ionospheric delay after passing through the ionosphere. A large number of models are currently being used to describe the delay that the ionosphere produces for electromagnetic signals propagating from satellites to receivers. If the delay is not corrected, it can have an important impact on the positioning accuracy of GNSS [9], especially for single-frequency PPP (SF PPP). Shi et al. proposed an improved method in which the deterministic representation is further refined by a stochastic process for each satellite with an empirical model for its power density [10]. The results of this method show that the single- and dual-frequency PPP exhibited enhanced convergence time, and the positioning accuracy of SF data is only improved by 25% [11]. Abd Rabbou et al. [12] developed an SF PPP model, and the improved model uses between-satellite-single-difference quasi-phase constrained GNSS observations. The GRoup and PHase Ionospheric Correction (GRAPHIC) method uses an SF code and carrier phase data to form an IF combined observation [13]. The positioning accuracy can reach several centimeters based on the GRAPHIC method, but the method requires two hours of convergence time [14].

The International GNSS Service (IGS) Ionosphere Working Group (IWG) was created in 1998. Several analysis centers (CODE (Center for Orbit Determination in Europe), ESA (European Space Agency), JPL (Jet Propulsion Laboratory), UPC (Universitat Politecnica de Catalunya), CAS (Chinese Academy of Sciences), and WHU (Wuhan University)) produce and release the post-processed Global Ionospheric Map (GIM) products, where the format is IONosphere EXchange (IONEX) [15,16]. The data portion of the GIM products utilized contain a total of 25 maps, with the latitude variation ranging from 87.5° to −87.5° with an interval −2.5° and the longitude ranging from −180° to 180° at 5° intervals. The accuracy of the GIM is improved based on the Spherical Harmonic (SH) expansion model [17]. The VTEC of each ionospheric grid point is obtained while using different interpolation methods, such as bilinear interpolation [18], inverse distance weighted function [19], and Kriging interpolation [20]. Lanyi and Roth [21] proposed a polynomial model for single-station TEC derivation, and a single-station receiver bias can be estimated while using this model [22]. This method was used by Lu et al. [23] to study the effects of ionospheric shell height on GPS-based TEC derivation by a single station, while Kao et al. [24] applied this method to multi-stations. Chen et al. [25] analyzed the applicability of the sophisticated Klobuchar model for VTEC in China. Liu et al. [26] had based on GPS to observe and analyze the fluctuation characteristics of TEC over China.

In GNSS data processing, the slant tropospheric delay on the signal propagation path between the satellite and the receiver is usually mapped to the zenith direction via a mapping function. The Zenith Tropospheric Delay (ZTD) consists of the Zenith Hydrostatic Delay (ZHD) and the Zenith Wet Delay (ZWD, i.e., the zenith non-hydrostatic delay). The models of tropospheric delay estimation usually adopt the Hopfield model and Saastamoinen model. There are roughly two types of mapping functions. The first type are empirical models, which merely require the epoch time and the approximate coordinates of the receiver station, such as the Niell Mapping Function (NMF) [27] and the Global Mapping Function (GMF) [28]. The second type of mapping function is based on a large number of weather model analyses at the epoch of the observations, such as Isobaric Mapping Functions (IMF) [29] and Vienna Mapping Functions (VMF1) [28]. Huang et al. [30] proposes a new Asian single site tropospheric correction model, called the Single Site Improved European Geostationary Navigation Overlay Service model (SSIEGNOS).

SF PPP cannot be combined or differenced to be eliminated or attenuated by a part of the error, like dual-frequency PPP. If there are some information related to the parameters, and this information is accurate enough, the precision of traditional SF PPP can be largely improved and the convergence time can be shortened. Zhang et al. [31] studied real-time GIM and its application in SF positioning, Aggrey and Bisnath [32] studied the effect of atmospheric-constrained on the convergence time of dual- and triple-frequency PPP, and Gao et al. [33] applied the Inertial Navigation System (INS) to the ionosphere-constrained PPP to overcome the drawbacks that accompany unexpected and unavoidable substandard observation environments.

This study uses the GIM products and the tropospheric zenith path delays from the IGS as the constrained information for SF PPP. The organization of this study is as follows: the next step details the methods for mathematical models of standard SF GPS PPP, troposphere-constrained SF GPS PPP, and ionosphere-constrained SF GPS PPP, and the details of the four interpolation methods of GIM products. Afterwards, in Section 3 we compare and analysis the influence of interpolation methods and constraint methods on SF GPS PPP convergence time, and finally in Section 4 we draw conclusions. Later, we will study ionospheric delay and tropospheric delay prediction models to provide virtual atmospheric delay observations for real-time PPP and also provide a priori information for the constraint processing.

2. Methods

The GPS PPP observation models are first derived. Afterwards, different interpolation methods for GIM products are introduced.

2.1. GPS PPP Observation Models

This study uses an undifferenced and uncombined GPS PPP model as compared with the traditionally used IF-PPP model. These models include a standard SF GPS PPP model, troposphere-constrained SF GPS PPP models, ionosphere-constrained SF GPS PPP models, and SF GPS PPP models with troposphere- and ionosphere-constrained.

2.1.1. Standard Single-frequency GPS PPP Observation Model

The distance from the satellite to the receiver can be measured while using pseudorange and carrier phase observations, with the following expression [34]:

$$\{\begin{array}{ll}{P}_r^s\hfill & ={\rho }_r^s+\mathrm{c}\cdot (d{t}_r-d{t}^s)+I_r^s+M_w\cdot Z_w+M_h\cdot Z_h+d_r-d^s+{\epsilon }_r^s\hfill \\ L_r^s& ={\lambda }_1\cdot {\Phi }_r^s\\ & ={\rho }_r^s+\mathrm{c}\cdot (d{t}_r-d{t}^s)-I_r^s+M_w\cdot Z_w+M_h\cdot Z_h+{\lambda }_1\cdot (N_r^s+b_r-b^s)+{\xi }_r^s.\end{array}$$

(1)

where, $P_r^s$ and $L_r^s$ are the original pseudorange and the carrier phase ${\Phi }_r^s$ is multiplied by wavelength ${\lambda }_1$ for the specific receiver $r$ and satellite $s$; ${\rho }_r^s$ is the geometric distance from the satellite to the receiver; $\mathrm{c}$ is the speed of light in vacuum; $d{t}_r$ and $d{t}^s$ are the clock error of the receiver and satellite; $I_r^s$ is the slant ionospheric delay on the first frequency ${\mathrm{f}}_1$; $M_w$ is the wet mapping function; $Z_w$ is the zenith wet delay; $M_h$ is the hydrostatic mapping function; $Z_h$ is the zenith hydrostatic delay; $d_r$ is the frequency-dependent receiver uncalibrated code delay (UCD) with respect to satellite $s$; $d_{}^s$ is the frequency-dependent satellite UCD; $N_r^s$ is the integer phase ambiguity; $b_r$ and $b_{}^s$ are the uncalibrated phase delays (UPDs) for the receiver and satellite, which is frequency-dependent; and, ${\epsilon }_r^s$ and ${\xi }_r^s$ are the sum of measurement noise caused by the pseudorange and carrier phase observations and the error caused by the multipath effect. Other errors have been modeled in advance.

The original pseudorange and carrier phase observations equations are linearized, and the receiver clock offset only absorbs common part of frequency dependent receiver UCDs. Supposing that $m$ satellites are tracked simultaneously at a certain epoch by the receiver $r$, the standard SF PPP model can be written, as follows [35]:

$$\left[\begin{array}{l}{p}_r^1\hfill \\ l_r^1\hfill \\ ⋮\hfill \\ p_r^m\hfill \\ l_r^m\hfill \end{array}\right]=\left[\begin{array}{lllll}-{\mathbf{\mu }}_r\hfill & {\mathit{I}}_t\hfill & M_w\hfill & \mathit{K}\hfill & \mathit{R}\hfill \end{array}\right]\left[\begin{array}{l}\mathit{x}\hfill \\ d{\overline{t}}_r\hfill \\ Z_w\hfill \\ {\mathit{I}}_r^s\hfill \\ {\overline{\mathit{N}}}_r^s\hfill \end{array}\right]+\left[\begin{array}{c}{\epsilon }_r\\ {\xi }_r\end{array}\right],\begin{array}{cc}& {\mathit{Q}}_L\end{array},$$

(2)

with

$$\{\begin{array}{l}d{\overline{t}}_r=d{t}_r+d_r\hfill \\ {\overline{N}}_r^s={\lambda }_1\cdot (N_r^s+b_r-b^s)+d_{IF}^s-d_r,\hfill \end{array}$$

(3)

where ${\mathbf{\mu }}_r$ is the unit vector of the coordinate component between the receiver and the satellite; $\mathit{x}$ is the vector of the receiver position increments relative to the priori position; and, ${\mathit{I}}_t$ is a unit vector of $2\cdot m$ rows and one column, corresponding to the receiver clock parameter $d{\overline{t}}_r$. In matrix $\mathit{K}$, the elements for the corresponding $p_r^s$ and $l_r^s$ are 1 and −1, respectively, corresponding to the ionospheric parameter ${\mathit{I}}_r^s$. $\mathit{R}$ is the matrix corresponding to the ambiguity parameters ${\overline{\mathit{N}}}_r^s$, and the elements for the corresponding $p_r^s$ and $l_r^s$ are 0 and 1, respectively. ${\mathit{Q}}_L$ is the stochastic model of the observed minus computed (OMC) observables. IGS precise satellite clock was calculated by the IF combination observables and the IF combination of satellite UCDs was absorbed by the satellite clock offsets, as follows [2]:

$$d{t}_{I{F}_{12}}^s=d{t}^s+d_{I{F}_{12}}^s,$$

(4)

$$d_{I{F}_{12}}^s={\mathsf{\alpha }}_{12}\cdot d_1^s+{\mathsf{\beta }}_{12}\cdot d_2^s,$$

(5)

$${\mathsf{\alpha }}_{12}=\frac{{({\mathrm{f}}_1^{\mathrm{s}})}^2}{{({\mathrm{f}}_1^{\mathrm{s}})}^2-{({\mathrm{f}}_2^{\mathrm{s}})}^2},$$

(6)

$${\mathsf{\beta }}_{12}=\frac{{({\mathrm{f}}_2^{\mathrm{s}})}^2}{{({\mathrm{f}}_1^{\mathrm{s}})}^2-{({\mathrm{f}}_2^{\mathrm{s}})}^2}.$$

(7)

2.1.2. Troposphere-Constrained Single-Frequency GPS PPP Observation Model

According to the tropospheric delay product published by IGS, a virtual observation is used in the observation equation, and a troposphere-constrained SF GPS PPP is then constructed. The observation equation is expressed, as follows:

$$\left[\begin{array}{l}{p}_r^1\hfill \\ l_r^1\hfill \\ ⋮\hfill \\ p_r^m\hfill \\ l_r^m\hfill \\ {\stackrel{\wedge }{T}}_r\hfill \end{array}\right]=\left[\begin{array}{lllll}-{\mathbf{\mu }}_r\hfill & {\mathit{I}}_t\hfill & M_w\hfill & \mathit{K}\hfill & \mathit{R}\hfill \\ \mathbf{0}\hfill & \mathbf{0}\hfill & {\mathit{I}}_T\hfill & \mathbf{0}\hfill & \mathbf{0}\hfill \end{array}\right]\left[\begin{array}{l}\mathit{x}\hfill \\ d{\overline{t}}_r\hfill \\ Z_w\hfill \\ {\mathit{I}}_r^1\hfill \\ {\overline{\mathit{N}}}_r^s\hfill \end{array}\right]+\left[\begin{array}{c}{\epsilon }_r\\ {\xi }_r\\ {\epsilon }_{r,trop}\end{array}\right],\begin{array}{ccc}& {\mathit{Q}}_L,& {\mathit{Q}}_T\end{array}.$$

(8)

where ${\mathit{I}}_T$ is a unit vector of $\mathrm{m}$ rows and one column, corresponding to the zenith wet delay parameter $Z_w$; ${\stackrel{\wedge }{T}}_r$ is derived from external tropospheric products, the products for each IGS station at 5 min intervals provided by the IGS Analysis Center; ${\epsilon }_{r,trop}$ is the corresponding noise; ${\mathit{Q}}_T$ denotes the stochastic model of virtual troposphere observables; and, $\mathbf{0}$ is a zero matrix.

2.1.3. Ionosphere-Constrained Single-Frequency GPS PPP Observation Model

The ionosphere-constrained SF GPS PPP can add the GIM product as a virtual observation to constrain the ionospheric parameters as compared with the standard SF GPS PPP. The constraint equation for the observation is as follows [35]:

$$\left[\begin{array}{l}{p}_r^1\hfill \\ l_r^1\hfill \\ ⋮\hfill \\ p_r^m\hfill \\ l_r^m\hfill \\ {\stackrel{\wedge }{I}}_r^1\hfill \\ ⋮\hfill \\ {\stackrel{\wedge }{I}}_r^m\hfill \end{array}\right]=\left[\begin{array}{lllll}-{\mathbf{\mu }}_r\hfill & {\mathit{I}}_t\hfill & M_w\hfill & \mathit{K}\hfill & \mathit{R}\hfill \\ \mathbf{0}\hfill & \mathbf{0}\hfill & \mathbf{0}\hfill & {\mathit{I}}_I\hfill & \mathbf{0}\hfill \end{array}\right]\left[\begin{array}{l}\mathit{x}\hfill \\ d{\overline{t}}_r\hfill \\ Z_w\hfill \\ {\mathit{I}}_r^1\hfill \\ {\overline{\mathit{N}}}_r^s\hfill \end{array}\right]+\left[\begin{array}{c}{\epsilon }_r\\ {\xi }_r\\ {\epsilon }_{r,ion}\end{array}\right],\begin{array}{ccc}& {\mathit{Q}}_L,& {\mathit{Q}}_I\end{array}.$$

(9)

where ${\stackrel{\wedge }{I}}_r^s$ is derived from external ionospheric products with the corresponding noise ${\epsilon }_{r,ion}$; ${\mathit{I}}_I$ is a unit vector of $m$ rows and one column; and, ${\mathit{Q}}_I$ denotes the stochastic model of virtual ionospheric observables.

2.1.4. Troposphere- and Ionosphere-Constrained Single-Frequency GPS PPP Observation Model

If the external troposphere and ionospheric products are both used as virtual observations to constrain the observation equation, the following observation equation is obtained:

$$\left[\begin{array}{l}{p}_r^1\hfill \\ l_r^1\hfill \\ ⋮\hfill \\ p_r^m\hfill \\ l_r^m\hfill \\ {\stackrel{\wedge }{T}}_r\hfill \\ {\stackrel{\wedge }{I}}_r^1\hfill \\ ⋮\hfill \\ {\stackrel{\wedge }{I}}_r^m\hfill \end{array}\right]=\left[\begin{array}{lllll}-{\mathbf{\mu }}_r\hfill & {\mathit{I}}_t\hfill & M_w\hfill & \mathit{K}\hfill & \mathit{R}\hfill \\ \mathbf{0}\hfill & \mathbf{0}\hfill & {\mathit{I}}_T\hfill & \mathbf{0}\hfill & \mathbf{0}\hfill \\ \mathbf{0}\hfill & \mathbf{0}\hfill & \mathbf{0}\hfill & {\mathit{I}}_I\hfill & \mathbf{0}\hfill \end{array}\right]\left[\begin{array}{l}\mathit{x}\hfill \\ d{\overline{t}}_r\hfill \\ Z_w\hfill \\ {\mathit{I}}_r^1\hfill \\ {\overline{\mathit{N}}}_r^s\hfill \end{array}\right]+\left[\begin{array}{c}{\epsilon }_r\\ {\xi }_r\\ {\epsilon }_{r,trop}\\ {\epsilon }_{r,ion}\end{array}\right],\begin{array}{cccc}& {\mathit{Q}}_L,& {\mathit{Q}}_T,& {\mathit{Q}}_I\end{array}.$$

(10)

2.2. Different Interpolation Methods for GIM Products

It is necessary to interpolate the ionospheric grid point from GIM products to obtain the observation station VTEC before performing the ionosphere-constrained modeling. In this study, a comparative analysis of four interpolation methods is performed. These methods include nearest-neighbor extrapolation [36], bilinear interpolation [18], Junkins weighted interpolation, and bicubic interpolation. The principles of the Junkins weighted and bicubic interpolation methods are described below.

2.2.1. Junkins Weighted Interpolation

Figure 1 shows the schematic of Junkins interpolation. In Figure 1, the latitude and longitude of the interpolation point are $b$ and $l$, that is, $p(b,l)$, and the values of the four GIM grid points adjacent to it are $p_{11}(b_1,l_1),p_{12}(b_1,l_2),p_{21}(b_2,l_1)$, and $p_{22}(b_2,l_2)$. It is assumed that the TEC values at $p_{11},p_{12},p_{21}$, and $p_{22}$ are $TE{C}_{11},TE{C}_{12},TE{C}_{21}$, and $TE{C}_{22}$, respectively. The VTEC of point $p$ can be interpolated according to the VTEC of the surrounding four grid points. The interpolation formula used is as follows:

$$f(b,l)=w_1(x,y)\cdot TE{C}_{11}+w_2(x,y)\cdot TE{C}_{12}+w_3(x,y)\cdot TE{C}_{21}+w_4(x,y)\cdot TE{C}_{22},$$

(11)

The weighted function $w_n(x,y),\begin{array}{cc}& n=1,2,3,4\end{array}$ is as follows:

$$w_1(x,y)=\omega (1-x_p,1-y_p),$$

(12)

$$w_2(x,y)=\omega (x_p,1-y_p),$$

(13)

$$w_3(x,y)=\omega (x_p,y_p),$$

(14)

$$w_4(x,y)=\omega (1-x_p,y_p),$$

(15)

where

$$\omega (A,B)=A^2{B}^2(9-6A-6B+4AB),$$

(16)

$$x_p=\frac{l-l_1}{l_2-l_1},\begin{array}{cc}& y_p=\frac{b-b_1}{b_2-b_1}\end{array}.$$

(17)

We can use nearest-neighbor extrapolation when the $p$ point only has two proximate GIM grid points.

2.2.2. Bicubic Interpolation

The bicubic interpolation [37] method uses 16 adjacent points for interpolation (Figure 2). Bicubic interpolation is an extension of cubic interpolation in a two-dimensional space. The cubic interpolation kernel is an approximation of the convolution interpolation of the ideal sampling sinc function at [−2,2] that is based on the cubic polynomial [38]. Typically, cubic interpolation produces interpolation coefficients that are based on a third-order polynomial:

$$h(q)=\{\begin{array}{ll}1-(w+3){\left|q\right|}^2+(w+2){\left|q\right|}^3,& 0\le \left|q\right|\le 1\hfill \\ -4w+8w\left|q\right|-5w{\left|q\right|}^2+w{\left|q\right|}^3,& 1\le \left|q\right|\le 2\hfill \\ 0,& 2\le \left|q\right|\hfill \end{array},$$

(18)

where $q$ represents the distance between the interpolation point and the reference point, $w$ is a tunable parameter, and the best result is obtained by verifying $w=-0.5$ [38] with a large amount of data. Under these conditions, Equation (18) is simplified to:

$$h(q)=\{\begin{array}{ll}1-2.5{\left|q\right|}^2+1.5{\left|q\right|}^3,\hfill & 0\le \left|q\right|\le 1\\ 2-4\left|q\right|+2.5{\left|q\right|}^2-0.5{\left|q\right|}^3,\hfill & 1\le \left|q\right|\le 2\\ 0,\hfill & 2\le \left|q\right|\end{array},$$

(19)

Bicubic interpolation is similar to bilinear interpolation. It is also decomposed into two one-dimensional interpolations: horizontal interpolation and vertical interpolation. Figure 2 shows the interpolation process. First, vertical interpolation is performed to obtain four virtual values, and then the values of the $p$ point are obtained via interpolating based on the four virtual values. The calculation process is as follows:

$$\{\begin{array}{c}{P}_1=C_1{P}_{11}+C_2{P}_{21}+C_3{P}_{31}+C_4{P}_{41}\\ P_2=C_1{P}_{12}+C_2{P}_{22}+C_3{P}_{32}+C_4{P}_{42}\\ P_3=C_1{P}_{13}+C_2{P}_{23}+C_3{P}_{33}+C_4{P}_{43}\\ P_4=C_1{P}_{14}+C_2{P}_{24}+C_3{P}_{34}+C_4{P}_{44}\end{array},$$

(20)

namely

$$\left[\begin{array}{cccc}{P}_1& P_2& P_3& P_4\end{array}\right]=\left[\begin{array}{cccc}{C}_1& C_2& C_3& C_4\end{array}\right]\left[\begin{array}{cccc}{P}_{11}& P_{12}& P_{13}& P_{14}\\ P_{21}& P_{22}& P_{23}& P_{24}\\ P_{31}& P_{32}& P_{33}& P_{34}\\ P_{41}& P_{42}& P_{43}& P_{44}\end{array}\right],$$

(21)

$$\begin{array}{ll}F& =L_1{P}_1+L_2{P}_2+L_3{P}_3+L_4{P}_4\\ & =\left[\begin{array}{cccc}{P}_1& P_2& P_3& P_4\end{array}\right]\left[\begin{array}{c}{L}_1\\ L_2\\ L_3\\ L_4\end{array}\right]\\ & =\left[\begin{array}{cccc}{C}_1& C_2& C_3& C_4\end{array}\right]\left[\begin{array}{cccc}{P}_{11}& P_{12}& P_{13}& P_{14}\\ P_{21}& P_{22}& P_{23}& P_{24}\\ P_{31}& P_{32}& P_{33}& P_{34}\\ P_{41}& P_{42}& P_{43}& P_{44}\end{array}\right]\left[\begin{array}{c}{L}_1\\ L_2\\ L_3\\ L_4\end{array}\right],\end{array}$$

(22)

where

$$\begin{gathered} C_1=h(b-b_1),C_2=h(b-b_2),C_3=h(b_3-b),C_4=h(b_4-b) \\ {L}_1=h(l-l_1),L_2=h(l-l_2),L_3=h(l_3-l),L_4=h(l_4-l) \\ \left[\begin{array}{cccc}{P}_{11}& P_{12}& P_{13}& P_{14}\\ P_{21}& P_{22}& P_{23}& P_{24}\\ P_{31}& P_{32}& P_{33}& P_{34}\\ P_{41}& P_{42}& P_{43}& P_{44}\end{array}\right]=\left[\begin{array}{cccc}TE{C}_{11}& TE{C}_{12}& TE{C}_{13}& TE{C}_{14}\\ TE{C}_{21}& TE{C}_{22}& TE{C}_{23}& TE{C}_{24}\\ TE{C}_{31}& TE{C}_{32}& TE{C}_{33}& TE{C}_{34}\\ TE{C}_{41}& TE{C}_{42}& TE{C}_{43}& TE{C}_{44}\end{array}\right] \end{gathered}$$

then

$$f(b,l)=\left[\begin{array}{cccc}h(b-b_1)& h(b-b_2)& h(b_3-b)& h(b_4-b)\end{array}\right]\left[\begin{array}{cccc}TE{C}_{11}& TE{C}_{12}& TE{C}_{13}& TE{C}_{14}\\ TE{C}_{21}& TE{C}_{22}& TE{C}_{23}& TE{C}_{24}\\ TE{C}_{31}& TE{C}_{32}& TE{C}_{33}& TE{C}_{34}\\ TE{C}_{41}& TE{C}_{42}& TE{C}_{43}& TE{C}_{44}\end{array}\right]\left[\begin{array}{c}h(l-l_1)\\ h(l-l_2)\\ h(l_3-l)\\ h(l_4-l)\end{array}\right].$$

(23)

We can use bilinear interpolation when the number of GIM grid points around point $p$ is less than 16.

Tropospheric delay interpolation is similar to inverse distance-weighted interpolation [19], and it is performed according to the relationship between the tropospheric delay value and time in the tropospheric products.

3. Experimental Data Set and Analysis

The PPP results based on raw observations are evaluated with 78 global IGS stations for the whole month of September 2017 to verify the improved troposphere- and ionosphere-constrained for GPS SF PPP performance. For every station, the 24-h observations were divided into eight 3-h sessions to evaluate the performance of PPP. Therefore, a total of 18720 experiments were tested and compared. Figure 3 shows the distribution of these stations in the world. Figure 4 show a flowchart of our study procedure, while

Table 1. Models and strategies used in single-frequency Precise Point Positioning (SF PPP).

Item.	Models/Strategies
Constellation	GPS
Data span	1–30 Sep 2017
Observations	GPS raw code and phase observations
Frequency selection	L1
Sampling rate	30 s
Satellite cut-off elevation	7°
Positioning mode	Static/kinematic (simulated kinematic)
Estimator	Kalman filter
Satellite orbit and clock	Final by IGS orbit and clock products
Weighing strategy	A priori precision of 0.003 and 0.3 m for raw phase and code [35]
Satellite DCB	Correct using CODE Differential Code Bias (DCB) products
Ionospheric delay	Klobuchar model [39] (priori) or GIM products from IGS (constrains)
Ionospheric noise	Random walk [10]
Tropospheric delay	Saastamoinen model [40] (priori) or products from IGS (constrains)
Tide displacement	Solid Earth, pole and ocean tide were corrected [41]
Sagnac effect	Corrected model [42]
Relativistic effect	Corrected model [43]
Satellite and receiver antenna	Phase Center Offset (PCO) and Variations (PCVs) values from igs14.atx
Reference coordinates	IGS SINEX files
Phase ambiguities	Float values

lists the models and strategies that were used in this study.

The PPP performance in terms of convergence time in the horizontal and vertical components is evaluated at a different confidence level (the 68% and 95%) in two modes (static and kinematic). The criterion for judging convergence is when the positioning error in the horizontal and vertical components is less than 0.5 m (95% level) and 0.3 m (68% level), respectively [11].

3.1. The Effect of Interpolation Method on Ionosphere-constrained GPS Single-frequency PPP

For the convenience of description,

Table 2. The four interpolation schemes.

Abbreviation	Item	Abbreviation	Item
NENE	Nearest-neighbor extrapolation	BILI	Bilinear interpolation
BICU	Bicubic interpolation	JUNK	Junkins weighted interpolation

shows the abbreviations, according to the four interpolation methods that were introduced above. For GIM grid products, this study discuss the convergence time of the four interpolation methods under ionosphere-constrained and the average time of calculation of 78 stations.

Figure 5 and Figure 6 show the percentage of convergence time. From these two pictures, whether in the 95% or 68% level, and whether it is static mode or kinematic mode, the convergence time is less than 10 min, the percentage of the BILI method calculation result is slightly lower than the other three methods, and in the subsequent convergence time less than 20 min, 30 min, etc., the percentage of BILI method calculation results has been higher than the other three methods, as shown in the green histogram. That is to say, the convergence time of the BILI interpolation method is the shortest when compared to the other three methods. Figure 7 shows the statistical results of the root mean square (RMS) error in the horizontal and vertical components in the static and kinematic mode, respectively. The RMS errors are based on the statistics of the last 15 min of the position solution error [44]. In the static mode, the percentage of BILI method calculations exceeds 40% in the horizontal and vertical components when the RMS is less than 5 cm, and the other three methods are lower than 40% and, in the subsequent statistical stages, the percentage of the BILI method is higher than the other three methods to reach 100%. In the kinematic mode, the percentage of BILI method calculation results in the horizontal component is 60% when the RMS is less than 15 cm, and the other three methods are lower than this value; while, in the vertical component, the BILI method exceeds 40%, and the other three methods not reaching 40%, and the percentage of BILI methods in other statistical stages is higher than the other three methods. In summary, the BILI method has higher final positioning accuracy than the other three methods in both the static and kinematic modes in the horizontal or vertical component.

The average convergence time (removed the data that the convergence time is equal to 180 min, that is, the data that does not convergence within three hours) is shown in

Table 3. Average convergence time of GPS static and kinematic SF PPP in different interpolation schemes.

Scheme	Static SF PPP (min)				Kinematic SF PPP (min)
	H(95%)	H(68%)	V(95%)	V(68%)	H(95%)	H(68%)	V(95%)	V(68%)
NENE	3.21	15.10	17.20	40.60	18.45	77.34	78.46	134.10
BILI	3.14	14.40	15.18	32.81	16.10	73.00	68.95	131.50
JUNK	3.87	17.55	19.44	42.92	19.00	78.30	79.44	134.60
BICU	3.42	15.70	17.68	41.30	19.30	78.30	79.01	134.30

, as compared to the JUNK method with the longest average convergence time, the average convergence time of BILI method at the 95% level is improved from 3.87 to 3.14 min (18.9%) in horizontal (H(95%)) and by from 19.44 to 15.18 min (21.9%) in the vertical components (V(95%)) in the static mode, while, in the kinematic mode, it increased by 15.3% in the horizontal component and 13.2% in the vertical component. At the 68% level, the average convergence time in the horizontal component (H(68%)) significantly reduced by 17.9% from 17.55 to 14.40 min and by 23.6% from 42.92 to 32.81 min in the vertical component (V(68%)) in the static mode, and in the kinematic mode, lessened by 6.8% and 2.3% in the horizontal component and vertical component, respectively. The calculation time of these 18720 tests was normalized to FLRS station (

Table 4. Calculation time of FLRS station in SF PPP in different interpolation schemes.

Interpolation Method	Static (sec)	Kinematic (sec)
NENE	6.905	6.821
BILI	6.550	6.512
BICU	6.775	6.797
JUNK	6.429	6.474

). We used a Lenovo computer, which is configured with Windows 8.1 professional computer with an i5-3470 Central Processing Unit (CPU) and 4GB of installed Random Access Memory (RAM). From Table 4, it can be concluded that the BILI interpolation method calculation time is almost the same as the Junkins weighted interpolation method, which takes the shortest time. According to all of the above analysis, the interpolation method that is adopted for GIM products in the following ionosphere-constrained model is bilinear interpolation.

3.2. Single-Frequency PPP with the Constraints

In this section, we analyse the effects of single-frequency GPS on unconstrained, troposphere-constrained, ionosphere-constrained, and their dual constrains on convergence time.

Table 5. The four constraint schemes.

Abbreviation	Item
TIC0	Standard single-frequency PPP model
TIC1	Troposphere-constrained single-frequency PPP model
TIC2	Ionosphere-constrained single-frequency PPP model
TIC3	Troposphere- and ionosphere-constrained single-frequency PPP model

summarizes the four constraint schemes.

Figure 8 and Figure 9 show the percentage statistics of the effects of different constraint methods on the convergence time of GPS. It can be concluded from these two figures that in the horizontal component, the convergence time of the TIC3 has the largest percentage in the same range, and the TIC1 has less influence on the convergence time. In the vertical component, the percentage of convergence time of the TIC3 in the same range is still the largest, but the influence of TIC1 on the convergence time is higher than the percentage of TIC2 in a certain range. On the whole, the convergence performance of TIC0 is worse in both components when not using GIM products or tropospheric products as constraints. The convergence performance is improved as compared to TIC0 after considering TIC1 or (and) TIC2 in GPS SF PPP processing.

At the same time,

Table 6. Average convergence time of GPS static and kinematic SF PPP in different constrained schemes.

Scheme	Static SF PPP (min)				Kinematic SF PPP (min)
	H(95%)	H(68%)	V(95%)	V(68%)	H(95%)	H(68%)	V(95%)	V(68%)
TIC0	27.17	42.11	23.21	39.34	61.63	103.07	81.84	134.81
TIC1	26.07	40.29	18.17	27.98	58.99	100.90	72.28	131.58
TIC2	3.15	14.39	15.18	32.82	16.06	73.29	69.01	131.40
TIC3	3.04	14.08	9.57	21.05	15.70	71.53	57.41	127.59

shows that TIC1 performs a notably better than TIC0 in the vertical component and marginally better than TIC0 in the horizontal component. TIC2 performs notably better than TIC0 and TIC1 in the horizontal component, while, in the vertical component, TIC2 is better than TIC0 and worse than TIC1 only in 68% level static mode. TIC3 performs best among the handling schemes. Compared with TIC0, the average convergence time of TIC1 at the 95% level is improved from 27.17 to 26.07 min (4.0%) in horizontal and by from 23.21 to 18.17 min (21.7%) in the vertical components. At the 68% level, the convergence time in the horizontal component reduced by 4.3% from 42.11 to 40.29 min and by 28.8% from 39.34 to 27.98 min in the vertical component. However, in kinematic mode, the average convergence time of TIC1 is not much higher than that of TIC0, with a maximum increase of 11.7% and a minimum increase of 2.1%. TIC2 has a larger boost when compared with TIC0. In static mode, it is increased by 16.6% in the vertical component by 68% level, and the maximum increase is 88.4% in the horizontal component of 95% level; in kinematic mode, they are 2.5% and 73.9%, respectively. TIC3 has the largest boost as compared with the other two methods. In static mode, it reaches 88.8% from 27.17 to 3.04 min and the lowest is 46.5% from 39.34 to 21.05 min In the kinematic mode, the minimum increase is 5.4%, and the maximum increase is 74.5%.

Three stations (FLRS, TRO1, and ZECK station) were randomly selected from 78 IGS stations to measure RMS error in the horizontal component, vertical component, and three-dimensional. It can be seen from Figure 10 and Figure 11 that the final positioning error reaches the sub-decimetre level in the kinematic mode, and it reaches the centimetre level in the static mode. However, the final positioning errors are substantially close in the four schemes.

4. Conclusions

In this study, we discussed the effects of four methods, including NENE, BILI, BICU, and JUNK, for the percentage of convergence time, the average convergence time and the calculation time of single station of SF TIC2 PPP. The numerical results show that: (1) the percentage of the convergence time of the estimated parameters through the BILI method is the shortest in every condition; and, (2) in the static mode, when the RMS is less than 5 cm, only the percentage of the BILI method calculation exceeds 40%. While in the kinematic mode, when the RMS is less than 15 cm, only the BILI method is more than 60% and 40% in the horizontal and vertical component, respectively. Under the other identical thresholds of the RMS, the percentage of the convergence time of the estimated parameters through BILI method is higher than the other three methods; (3) in the static mode, the average convergence time of the estimated parameters through BILI method is reduced by at least 2.2% when compared to the other three methods and the maximum is shortened by 23.6%; while, in the kinematic mode, they are 1.9% and 16.6%, respectively; and, (4) for the average computation time of a station, the BILI method has 0.121 s more than the shortest calculation time in the static mode and 0.038 s in the kinematic mode. Therefore, we chose the BILI method for calculating the TEC at the required time when using the TIC2 or TIC3 calculation based on the GIM grid products.

In order to verify the different constraint methods, a total of 18,720 tests were performed with 78 stations. The experimental results revealed the following findings: (1) the TIC1 convergence time percentage is slightly higher than the TIC0 in each convergence time threshold, while the TIC2 is larger and the TIC3 is the largest. That is to say, the convergence can be accelerated under the constraint conditions; (2) as compared with TIC0 method, the average convergence time of the TIC1 method is shortened by at least 4.0% and, at most, 28.9% in static mode. The TIC2 method has a larger percentage of shortening, at least up to 16.6%, the maximum up to 88.4%. However, in the vertical component, the TIC1 method is larger than the TIC2 method at the 68% level. Compared with TIC0, the TIC3 method is shortened by 46.5% at least and 88.8% at most. Additionally TIC3 method is better than TIC1 method and TIC2 method. In kinematic mode, the percentage of time that is shortened after the constraint is less than in static mode. The average convergence time of the TIC1 method is at least 2.1% shorter than the TIC0 method, and the maximum is shortened by 11.7%. The TIC2 method has a larger percentage of shortening compared to the TIC1 method, with a minimum of 2.5% and a maximum of 73.9%. The TIC3 method is superior to the TIC1 method and the TIC2 method, with a minimum of 5.4% reduction and a maximum reduction of 74.5% as compared to TIC0; and, (3) through comparative analysis, it is concluded that, the final positioning error of TIC0, TIC1, TIC2 and TIC3 methods are basically the same, no method is better. After the above analysis and conclusions, under the constraints of atmospheric error products, convergence can be accelerated, and there is no impact on the final accuracy improvement.