A novel GNSS precise point positioning method based on spherical harmonic expansion

ABSTRACT The Global Navigation Satellite System (GNSS) precise point positioning (PPP) from carrier-phase data is influenced by complex elements. The error sources are divided into two categories according to their relation to the signal propagation path. A new GNSS PPP algorithm based on spherical harmonic expansion (PPP-SH) is proposed, which expresses the error adjustments related to the signal propagation path through the spherical harmonic expansion, such as high-order ionosphere delay, zenith wet delay, phase wind-up delay, multipath delay and so on. Carrier-phase data of multiple continuous epochs are used to achieve the sequential sliding PPP solution, and the truncated singular value decomposition and iteration method by correcting characteristic values are used to solve the ill-conditioned equations in the PPP-SH algorithm. The reliable results can be obtained by using dual-frequency ionospheric-free carrier phase data when spherical harmonic expanding to degree 5. The standard deviations of E, N and U are 0.48, 0.43 and 0.27 cm, respectively. The addition of velocity parameters has positive influence on accuracy of positioning. The best positioning results are obtained by using the precise ephemeris and precise clock bias provided by ESA and the standard deviations of E, N and U are 0.41, 0.43 and 0.26 cm, respectively.


Introduction
Global Navigation Satellite System (GNSS) has been commonly used in navigation and positioning for its low cost, high precision, all-weather, vast coverage and great efficiency (Cheng et al., 2022;Dao et al., 2022;J. Guo et al., 2016).GNSS precise point positioning (PPP) from carrier phase data can obtain highprecision positioning results (Carlin et al., 2021;Zhang et al., 2021;Zumberge et al., 1997).
GNSS PPP has been widely used in the geodesy research.For example, an earthquake of magnitude 7.8 occurred in Turkey at 4:17 local time on 6 February 2023.Many scholars calculated the displacement of local GNSS IGS station by using PPP.A lot of scholars and scientific research institutions pay attention to higher accuracy of GNSS PPP.Many PPP software have programmed during decades years, such as RTKLIB (Takasu & Yasuda, 2009), Bernese 5.2 software (Dach et al., 2015), PRIDE PPP-AR (Geng et al., 2019) and so on.
In GNSS PPP, the main error sources are divided into two categories on the basis of whether they are related to the signal propagation path or not.Errors independent of the signal propagation path mainly include satellite clock bias, satellite orbit error, relativistic effect, receiver clock bias, earth rotation error and tidal error (Jiao & Song, 2022).Errors related to signal propagation path include ionospheric delay error, tropospheric delay error, phase wind-up and multipath effect (Hardwick & Jeffery, 1995;J. Guo et al., 2023;Ning & Yuan, 2017;Qi et al., 2021).
Ionosphere delay is important error source of PPP (J.Guo et al., 2009).A large number of electric ions distributed in the ionosphere will change the signal in the signal propagation path, resulting in errors of several metres or even tens of metres in positioning.When only single-frequency observation data can be received, using ionospheric models to eliminate ionosphere delay errors needs to be considered (Zhang et al., 2018), such as Klobuchar model (Klobuchar, 1987), but the accuracy of this model is not high.When the dual-frequency observation received by the GNSS receivers, the ionospheric delay errors are removed by using the dual-frequency ionospheric combination (Hoque & Jakowski, 2008;J. Guo et al., 2023;Qi et al., 2021), but the dual-frequency ionospheric combination is able to remove the first-order main term of the ionospheric error merely (Hoque & Jakowski, 2008).
Troposphere delay, segmented into Zenith Hydrostatic Delay (ZHD) and Zenith Wet Delay (ZWD), also has noticeable effect on results of PPP (Askne & Nordius, 1987;Zhou et al., 2022).In GNSS PPP, models are usually used to reduce the influence of tropospheric ZHD.The commonly used models are Saastamoinen model (Saastamoinen, 1972;Zhang et al., 2019), Hopfield model (Hopfield, 1969), and so on.These models are highly dependent on meteorological data and empirical parameters, and the real state of troposphere can't be represented completely by models (Penna et al., 2001), the positioning accuracy is limited.The tropospheric ZWD is usually estimated in PPP, but it is relatively difficult to estimate it accurately because of the intense vapour activity in troposphere (Bevis et al., 1996).
The effect of multipath on positioning can't be ignored (J.Y. Guo et al., 2014).After decades of development, some receivers can reduce the influence of multipath effect, but multipath effect still has a great influence on results of PPP (Dong et al., 2016), it still needs to be treated as an error related to GNSS signal propagation path (Hardwick & Jeffery, 1995).The influence of phase wind-up on results of PPP can reach sub-decimetre level at the highest (Georg, 2010), usually phase wind-up is also related to the signal propagation path, so it needs to be treated as an error related to the signal propagation path.
Spherical harmonic expansion has been a fitting means widely applied in geodesy.Considering the errors in PPP related to signal propagation path, spherical harmonic expansion can be used to infinitely approximate its true propagation path.A PPP method based on spherical harmonic expansion is proposed to correct the above errors related to signal propagation path, which uses spherical harmonic expansion to represent the error correction related to signal propagation path of GNSS signal, namely PPP-SH algorithm.Compared to traditional PPP method, PPP-SH regards the ZWD, high-order ionospheric errors, multipath errors and phase wind-up as an entirety, simplifies the equation expression and takes the errors which are ignored in the traditional PPP method into concern.The disadvantages of PPP-SH is that more parameters are led to PPP in PPP-SH, the amount of parameters is related to the degree of SH expansion.
In section 2, the PPP observation equation and error equation based on spherical harmonic expansion are derived, and the solutions to the possible problems in the PPP-SH are put forward.In section 3, the PPP-SH algorithm is verified by using GNSS observation data provided by IGS, corresponding precise ephemeris and precise clock bias.In section 4, some possible influencing factors are discussed.In section 5, research of this paper is concluded.

Methods
According to the basic principles of GNSS PPP, the observation equation of dual-frequency carrier phase is expressed as where i represents the number of the epoch, j represents the number of satellite; φ IF represents dualfrequency Ionospheric-Free (IF) carrier phase linear combination; R j i represents the geometric distance between the satellite j and the receiver at epoch i; c represents speed of light; t r represents the receiver clock bias; t s represents the satellite clock bias; ρ trop À � j i represents the tropospheric delay error of the satellite j on the signal propagation path at epoch i; ρ ion ð Þ j i represents the ionospheric delay error of satellite j on the signal propagation path at the epoch i; λ represents the wavelength; N represents the integer ambiguity; b r and b j s represent the receiver hardware delay and satellite j hardware delay, respectively, and they are processed by using IGS antenna files.ε φ represents the residual of carrier phase observation data, including phase wind-up delay, relativistic delay and so on.
Using spherical harmonic expansion to represent the error correction terms such as tropospheric delay error, ionospheric delay error and phase wind-up error, observation equation of GNSS PPP-SH can be expressed as where n represents the degree of spherical harmonic expansion, m represents the order of spherical harmonic expansion, N max is the maximum degree of the spherical harmonic expansion; P nm ðcos eÞ represents the associated Legendre function in degree n and order m, a multi-core parallel acceleration method is used to calculate the associative Legendre function; C nm and S nm respectively represent coefficients of the spherical harmonic expansion in degree n and order m; e j i and α j i respectively represent the elevation angle and the azimuth angle between the station and the satellite j at epoch i; and other parameters have the same meanings as those in Equation (1).
Assuming that the coordinate of satellite j at epoch i is and the coordinate of the station is ðX; Y; ZÞ.The motion state of the station is considered, including static and dynamic, and the geometrical dis-tanceR j i between the satellite and the station is ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where t represents time interval for two adjacent epochs, and v x , v y , v z represent the velocity components of the station in X, Y, Z three directions, respectively.If the station has no velocity component, then β ¼ 0; if the station has velocity component, then Assuming that ,Equation ( 2) can be simplified as The GNSS PPP-SH is expressed by Equation ( 4), and the coefficient of spherical harmonic expansion is directly solved as parameters to correct the error terms related to the altitude angle and azimuth angle between receiver and satellite.
The error equation of PPP-SH can be obtained from Equation (4), and expressed as where v φ IF À � j i represents the carrier phase observation data correction value of the satellite j at epoch i.
The error equation is expanded by Taylor series at the approximate coordinates ðX 0 ; Y 0 ; Z 0 Þ of the station, and the first-order term is reserved.The linearised expression of the error equation of GNSS PPP-SH is obtained and expressed as where R j i � � 0 represents geometrical distance between the approximate coordinates of the station and satellite j at epoch i; dXrepresents the correction value of the X 0 , dY represents the correction value of the Y 0 , dZ represents the correction value of the Z 0 ; and dv Xr , dv Yr , dv Zr represent the adjustment value of velocity components in three directions, respectively.
represents the initial value of integer ambiguity of satellite j at epoch i; N j i represents the integer ambiguity of satellite j at epoch i; dN represents the correction value of integer ambiguity; and d t r ð Þ i represents the receiver clock bias at epoch i.Because the clock bias of the receiver can't be well fitted by polynomial, it is taken as a parameter to be solved (Dai et al., 2001).
Assuming that Then Equation ( 6) can be simplified as The linearisation expression of error equation is expressed by matrix form as where P represents the weight matrix determined by altitude angle of satellites, P ¼ sin e σ 2 , and σ 2 represents the variance of observation.
When dual-frequency carrier phase observation data is received, IF can be applied to eliminate effects of the first-order main term of ionospheric delay (J.Guo et al., 2009;Ning & Yuan, 2017), and IF is expressed as where φ IF is the observation of the carrier phase ionosphere-free linear combination, f 1 and f 2 are frequencies corresponding to the GNSS signals, the frequency of f 1 is 1575.42MHz and the frequency of f 2 is 1227.6MHz.φ 1 and φ 2 respectively represent the carrier phase observations corresponding to two frequencies.
Assuming that coefficient matrix of error equation of PPP based on spherical harmonic expansion is B, correction vector matrix of carrier phase observation values is V, constant term matrix is L, parameters matrix is X, and they can be expressed as Because parameters solved by PPP-SH are more than PPP, it is necessary to combine multiple epochs to solve the parameters.Build a container named data pool that can hold i epoch observations based on sequential adjustment.First i epochs observations are stored into data pool and used to position when reading observations, then first epoch data in data pool is deleted when sliding window reads one epoch afterwards, and the terms calculated by first epoch data in B, V,L are deleted at the same time, other terms are reserved.
Assuming that B 0 , V 0 , L 0 are the matrices calculated by first epoch observations, firstly, B 0 is divided into three matrices, V 0 and L 0 are divided into two matrices, respectively.
0 and L 1 0 are deleted totally.Newly read epoch observations by sliding are used to calculate related terms.New coefficient matrix, correction value matrix and constant value matrix are constructed.Then splice the new matrices with the remaining matrices, respectively, and spliced matrices are used to positioning.; A slight change in the matrix results in a great change in the results of the equation, which is called illconditioned problem (Kirsche & Böckmann, 2005).In view of the ill-conditioned equation problems that may occur in the PPP-SH, the conditional-number method is used to discriminate the equation (Kirsche & Böckmann, 2005).If there is no ill-conditioned problem in the equation, least square method can be used to solve it directly.If the equation has ill-conditioned problems, truncated singular value decomposition method can be applied to obtain the initial value, and the initial value is used for the least square spectrum iteration to obtain results (Ge et al., 2016).According to the calculation process above, the Figure 1 flow chart of PPP-SH algorithm can be obtained.

Results and analysis of calculation examples
The data used in this paper are provided by Wuhan IGS Data Center (http://www.igs.gnsswhu.cn/index.php/home/data_product/igs.html),including GPS observations, IGS precise ephemeris and IGS precise clock bias.The GPS observations' sampling interval is 30 s, precise ephemeris' sampling interval is 15 min, precise clock bias' sampling interval is 5 min.The date of observations, precise ephemeris and precise clock bias is 15 January 2022.The selected IGS stations are IQAL, BJFS, BOGT and TIDB, which are located in the low, middle and high latitudes, respectively.They are all static stations.The stations are stable and wellqualified, so they can be used to verify PPP-SH algorithm.The specific information of the stations is listed in Table 1.
The PPP-SH solution strategies are listed in Table 2.The GPS 24 h observations of IQAL station on 15 January 2022, where sampling interval is 30 s, is used to positioning by PPP-SH when spherical harmonic expanded to degree 8. Considering the efficiency and accuracy of the solution, in this paper, the epoch number of the data pool is set to 15, and the results are compared with the high-precision coordinates solved by IGS on 15 ,January 2022, and reduced to ENU coordinate system.For the sake of checking accuracy of PPP-SH algorithm, results solved by PPP-SH are compared with the results of PPP solved by Bernese5.2and CSRS-PPP (https://webapp.csrs-scrs.nrcan-rncan.gc.ca/ geod/tools-outils/ppp.php),and the results are listed in Table 3.
As shown in the above table, positioning accuracy of PPP-SH in three directions is all within 4 cm, and the standard deviation of E-direction is 0.46 cm, N-direction is 0.39 cm and U-direction is 0.30 cm.The ambiguity was processed by Go-back method

Influence of velocity parameters on positioning results
To explore effects of velocity parameters on results of PPP-SH, static station is be used to positioning but added the velocity parameters.the GPS 24 h observations of IQAL station on 15 ,January 2022, where sampling interval is 30 s, is used to positioning by the PPP-SH when spherical harmonic expanded to degree 8.When other conditions are completely same, the results without velocity parameters and with velocity parameters are calculated, respectively.The results are shown in Table 4.
As listed in Table 4, after adding velocity parameters, the accuracy of E and N directions is improved, and the accuracy of U direction is decreased.Generally speaking, the addition of the velocity parameters has a positive influence on the results of PPP-SH in most cases, and its negative influence is acceptable.

Influence of spherical harmonic expansion degree on positioning results
To explore the influence of the degrees of spherical harmonic expansion on the positioning, the GPS 24 h observations of IQAL station on 15 ,January 2022, where sampling interval is 30 s, is used to positioning by the PPP-SH expanded to degree 5, 6, 7, respectively.The results of PPP-SH are shown in Table 5.
Table 5 shows that results of PPP-SH when spherical harmonic expanded to degree 5 are all in centimetre magnitude, and standard deviation of positioning in E direction is 0.48 cm, that in N direction is 0.43 cm and that in U direction is 0.27 cm.Compared with Table 4, it is obvious that the results of degree 5 are superior to that of degree 4. When the spherical harmonic expanded to the degree 6, the standard deviation in E direction, N direction and U direction is 0.46, 0.41 and 0.30 cm respectively.The results of PPP-SH when spherical harmonic expanded to the degree 6 are better in E and N directions than those of PPP-SH expanded to the degree 5, while the accuracy in U direction decreases slightly.When the spherical harmonic expanded to the degree 7, the standard deviation in E direction is 0.45 cm, that in N direction is 0.41 cm, and that in U direction is 0.32 cm.Compared with the results of PPP-SH when spherical harmonic expanded to degree 6, the accuracy of results of PPP-SH when spherical harmonic expanded to degree 7 are slightly improved in E and N directions, but the accuracy of U direction is slightly reduced.
It is acknowledged that the positioning accuracy will be improved with improvement of the degrees of spherical harmonic expansion when spherical harmonic expanding after degree 5. Therefore, if the positioning results and the calculation efficiency are considered, the better results can be obtained when spherical harmonic expanded to degree 5.It is suggested that when PPP-SH algorithm is used to positioning, if it is necessary to positioning quickly, spherical harmonic can only be expanded to degree 5 to meet the requirements.However, if the accuracy needs to be higher, it needs spherical harmonic to be expanded to a higher degree, but the efficiency of calculation will be sacrificed.

Applicability of PPP-SH for stations located in different areas
For the purpose of exploring the applicability of PPP-SH method for different latitude and longitude areas, stations located in different latitude and longitude are selected to verify the method, namely IQAL, BJFS, BOGT and TIDB.The date of the observations adopted is 15 ,January 2022, the type of observations is GPS, and the sampling interval is 30 s, spherical harmonic is expanded to degree 5, and the results of positioning are listed in Table 6.
As can be seen from Table 6, best average value can be obtained by using IQAL observations and the worst average value can be obtained by using TIDB observations in E direction.Best average value can be obtained by using IQAL observations and the worst average value can be obtained by using BOGT observations in N direction.The average value of IQAL in U direction is the best, and TIDB is the worst.The standard deviation values of IQAL are the best overall, while the standard deviation values of BOGT are the worst overall.This problem may be caused by the high content of water vapour, high air temperature, active electric ions, larger ionospheric delay and tropospheric delay in low latitude areas, resulting in poor positioning results.Therefore, when PPP-SH algorithm is used to positioning at low latitude, it is suggested to increase the degree of spherical harmonic expansion to improve the accuracy of the positioning.

Influence of precise ephemeris and precise clock bias on results of positioning
In PPP-SH, precise ephemeris and precise satellite clock bias are necessary.At present, precise ephemeris and precise clock bias are provided by many institutions all over the world (Chen et al., 2022;Jiao & Song, 2022).In this paper, precise ephemeris and precise clock bias provided by the German Research Center for Geosciences (GFZ), the Massachusetts Institute of Technology (MIT) and the European Space Agency (ESA) are used to positioning, in which the precise ephemeris and precise clock bias of GFZ and MIT are provided by Wuhan IGS Data Center.ESA precise ephemeris and precise clock bias are provided by IGS (https://cddis.nasa.gov/archive/gnss/products/2192/).The sampling interval of precise ephemerides are 15 min, and the sampling interval of precise clock biases are 5 min.The GPS 24 h observations of IQAL station on 15 ,January 2022, where sampling interval is 30 s, is used to positioning by the PPP-SH when spherical harmonic expanded to degree 5, and the initial number of epochs is 15.The results are shown in Table 7.
Table 7 shows that the best positioning results can be obtained by using the precise ephemeris and clock bias of ESA in PPP-SH algorithm; and the positioning results obtained by using the precise ephemeris and precise clock bias provided by MIT are relatively the worst.Therefore, it is suggested that the precise ephemeris and precise clock bias provided by ESA should be selected as far as possible in PPP-SH algorithm.

Conclusions
In this paper, aiming at processing of errors related to signal propagation path in GNSS PPP, a GNSS PPP method based on spherical harmonic expansion is proposed, which uses spherical harmonic expansion to express the error correction related to signal propagation path.Compared with the traditional PPP, this PPP-SH is simpler in expression and does not need to consider meteorological parameters and other factors.At present, the following conclusions can be drawn on the basis of the actual test.
(1) Using spherical harmonic function to represent the errors related to signal propagation path in GNSS PPP can obtain good positioning results.When spherical harmonic expanded to degree 8, the standard values of E, N and U directions are 0.46, 0.39 and 0.30 cm, respectively, so it can be seen that the PPP-SH is feasible.(2) With the improvement of spherical harmonic expansion degree, the accuracy of positioning results is also constantly improved.In most cases, the accuracy of position results can be cm and mm level.After the spherical harmonic expanded to degree 5, the positioning accuracy will improve with the improvement of the degree of spherical harmonic expansion.Therefore, if the positioning result and the calculation efficiency are considered, the spherical harmonic expanding to the degree 5 can obtain better results.(3) Stations located in low latitudes are more seriously affected by ionosphere and troposphere, so it is suggested to expand spherical harmonics to higher degree to obtain higher accuracy results when stations located in low latitudes positioning.(4) When other conditions are the same, the best positioning results are obtained by using precise ephemeris and precise clock bias provided by ESA, and the accuracy of results are better than that calculated by precise ephemeris and precise clock bias provided by other institutions in three directions.The positioning results obtained by using the precise ephemeris and precise clock bias provided by MIT are relatively the worst.Therefore, it is recommended to choose the precise clock bias and precise ephemeris provided by ESA as far as possible when PPP-SH is used to positioning.
This method can be applied in SPP, PPP and PPP-RTK to express the errors which are related to the signal propagation path.At present, the research of this paper still has some shortcomings and needs to be improved.
For example, the PPP-SH algorithm only expands spherical harmonic to degree 8 at the highest, and it can still try to expand spherical harmonic to a higher degree to calculate more accurate results of positioning.

Table 1 .
Information of IGS stations.

Table 2 .
The solution strategies of PPP-SH.

Table 3 .
The results of IQAL solved by PPP-SH when spherical harmonic expanded to degree 8, Bernese5.2software and CSRS-PPP online (m).

Table 5 .
The results of IQAL solved by PPP-SH when SH expanded to degree 5-7 (m).

Table 4 .
The results of IQAL solved by PPP-SH without adding velocity parameters and with adding velocity parameters (m).

Table 7 .
The results of IQAL station solved by different precise ephemerides and precise clock biases (m).

Table 6 .
The results of different areas stations solved by PPP-SH when spherical harmonic expanded to degree 5 (m).