Open-source optimization method for android smartphone single point positioning

Nowadays, a chip-scale Global Navigation Satellite System (GNSS) receiver is ubiquitous in smartphones. In a smartphone GNSS receiver, the least square (LS) or Kalman Filter (KF) is implemented to estimate the position. With the aim to improve the smartphone GNSS position accuracy, we propose a position-smoothing method considering more historical information than the traditional methods, i.e., LS and KF. More past states are regarded as unknowns, and a cost function is constructed to optimize these states. An open-source smartphone dataset from Google was used for testing the proposed method. The experimental results indicate that the proposed method outperforms the other conventional methods in position errors. In addition, we open the source code. We expect that the optimization method implemented in the smartphone GNSS position smoothing application can be an illustrative example to clearly introduce such an optimization method and a reference for its implementation, which might inspire some other meaningful and exciting applications in GNSS.


Introduction
Global navigation satellite systems (GNSS) are the dominant source for positioning, velocity and time (PVT) and are widely utilized in various applications (Zhang et al. 2018). Smart devices with chip-scale GNSS receivers can generate meter-level position results in a standalone manner in a GNSS friendly environment. In GNSS, the user receivers determine the PVT solutions with a least square (LS) method or Kalman filter (KF). In the LS-based PVT determination, the extracted pseudorange and pseudorange rate measurements, together with the position and velocity of satellites in view, are utilized to compute the PVT information at each epoch independently. The inherent relationship between the velocity and the position are ignored (Jiang et al. 2020). In the KF-based PVT determination a position-velocity (PV) dynamic motion model is utilized to estimate the PVT information and the velocity is utilized to smooth the position (Xu and Hsu 2019;Jiang et al. 2021). In fact, not only are adjacent positions correlated through velocity, but also all historical position information is correlated one by one through the velocity. Utilizing the inherent relationship among these states in the time domain might further decrease the position errors. Optimizing these states estimation together at an epoch is a prospective method to obtain better positioning results.
Nonlinear optimization originates from the simultaneous localization and mapping (SLAM) algorithm and can optimize the state estimation by considering historical information (Schubert et al. 2021;Indelman et al. 2012). Unlike the KF, all the state propagation and measurement constraints are used to construct a cost function. Optimal estimations of these states are solved by minimizing the cost function through an iterative method. Such nonlinear optimization methods have been extended to various applications, and their better performance versus KF has been comprehensively assessed and analyzed (Dellaert 2012). Therefore, it is of great significance to explore what benefits the optimization method could bring to GNSS receivers. We develop a nonlinear optimization method with its application in smartphone GNSS position smoothing. Contributions are summarized as: (1) An optimization program suitable for GNSS position smoothing applications is proposed and implemented; state propagation and measurement models are employed as the constraints in optimally estimating these states. In addition, the Levenberg-Marquardt (LM) algorithm is employed to solve for the optimal states in our code.
(2) An optimization method is developed, and the resultant code will be open-source. In addition, a detailed user manual is provided together with more example experiments, and we hope it will provide a reference or an illustrative example for implementing the optimization method suitable for other GNSS applications, i.e., signal tracking or vector tracking.
The remaining sections describe the methodology, including the PVT determination equations, the KF method and the proposed optimization method. Then, the following experimental sections present a test case, the results and the quantified and in-depth analysis are revealed and presented. Finally, the conclusions are drawn, and suggestions for future work are given.

Methodology
This section introduces the model and methods utilized in the GNSS PVT determination. First, the PVT determination model with the pseudorange measurements is presented; second, the LS and KF methods are introduced in the PVT determination, and the state and measurement models are presented.

PVT determination
In a GNSS receiver, considering the clock bias, the position vector u is usually defined as where u = x u y u z u is the three-axis position in the Earth-Centered-Earth-Fixed (ECEF) coordinates. Assuming there are N in-view satellites, the relationship of the pseudorange, (1) u = x u y u z u the user position, and the position of the satellites are given by (Hofmann-Wellenhof et al. 2007) where t b denotes the user local clock bias, x (i) sv y (i) sv z (i) sv denotes the coordinates of the ith satellite in the ECEF system, the constant c denotes the light speed, (i) denotes the pseudorange of between the ith satellite and the user. The pseudorange has been corrected, which means the ionospheric, tropospheric, relativistic and group delays and satellite clock offsets have been corrected. Similarly, the velocity computation model is given by pointing from the receiver to the ith satellite, u = v x v y v z denotes the user three-axis velocity in the ECEF system, sv,z refers to the ith satellite velocity in the ECEF coordinate, ̇( i) denotes the ith satellite pseudorange rate measurement, t d denotes the user local clock drift.

Least squares method
With the PVT model in (2-3), the LS algorithm is usually utilized to compute the PVT solutions. Expanding (2) at 0 u = x 0 u y 0 u z 0 u and ignoring higher-order terms (Jiang et al. 2020), the new equation is written as Similar to (4), the relationship between the change in velocity and the change of the pseudorange rates is given by where ̇ denotes the pseudorange noise vector, and the velocity can be computed using the similar LS method. With (4-5), the PVT is determined with the LS method; more details are presented in the reference (Jiang et al. 2020).

Kalman filter
In the KF method, the Position-Velocity dynamic state transformation model and the measurement model are given as denotes the state transformation matrix, the variable T 0 denotes the position updating duration, and its value is one second. k+1 = 3×3 3×3

3×3 3×3
, and k+1 is the measurement vector which is composed of position and velocity from the LS method. k+1 denotes the process noise vector, which is assumed to be subject to Gaussian distribution with zero mean value. k+1 denotes the measurement noise vector, which is assumed to be subject to Gaussian distribution with zero mean value. With the state and measurement models, the KF estimates the position through the five basic equations, which can be found in Jiang et al. (2020).

Optimization
Unlike the KF, all the past states are regarded as unknowns in the optimization method. These states and measurements are highly correlated through the state propagation model. We construct an optimization method that considers all the past state vectors as unknowns; the model is given by where i+1 denotes covariance matrix of the process noise vector i in (6), k+1 denotes covariance matrix with covariance matrix of the measurement noise vector k+1 in (7). Optimal estimation of the state vector = 2 , ..., k+1 can be solved by finding the minimal values of (8). The key is how to solve for the optimal estimates with (8). Here, we implement the Levenberg-Marquart (LM) algorithm to solve for the optimal estimations (Dellaert 2012). Note that Gauss-Newton and Powell's dog leg methods are also usable to find similar optimal estimates. The three-step solution of the LM method is described as follows: (1) In the first step, expanding the function f ( ) at ̃ 2 , ...,̃ k+1 ̃ 2 , ...,̃ k+1 and ignoring high order terms, the following equation is obtained Equation (9) can be simplified as (10) can be re-written in a general form (2) In the second step, differentiating (14) for Δ , and assuming the differentiated equation is set to zero In addition, to avoid the problem that the might be nonsingular, a damping factor is utilized and (15) is re-written as (Dellaert 2012): An increment is computed, and the ̃ is updated as: (3) Repeat the first and second steps until the iteration count reaches a predefined threshold or the increment Δ * reaches a predefined threshold.

Test case and results
A smartphone dataset collected by Google with Xiaomi 8 is used to assess the performance of the proposed method. A highaccuracy commercial position system serves as the reference (Fu et al. 2020). To evaluate the influence of the velocity errors on the proposed method, we divide the trajectory into three parts according to the velocity errors plotted in Fig. 1 (top). Corresponding horizontal position errors from the LS, KF, and optimization methods are presented in the bottom panel.
With the position errors plotted in Fig. 1 (bottom), we observe that the KF and the optimization method effectively smooth the horizontal position errors and reduce their rooted mean square error (RMSE) values. Further statistical analysis of these horizontal position errors is presented in Table 1. In the first part, compared with the LS method, the mean values of the KF and the optimization methods decrease by 21.2% and 41.2%, the corresponding RMSE errors from KF and the optimization methods decrease by 9.9% and 31.9%, and the corresponding figures are 4.7%, 10.1%, 17.4 and 27.8%, respectively, in the third part. However, the KF and the optimization methods seem ineffective in the second part. The mean values of the horizontal position errors from the LS, KF, and optimization  methods are almost identical, and only the RMSE values show minor differences. By observing the velocity errors plotted in Fig. 1 (top), the velocity errors are much larger than that of the first part and third part, which accounts for the phenomenon that the position smoothing methods fail to work. Comparing the velocity errors from the three different parts, we can observe that the position smoothing performance relies on the velocity accuracy. Therefore, the position smoothing schemes are effective in the first and third parts. In addition, the measurement errors covariance matrix, are fixed in the experiment, which is another reason for the smoothing failure in the second part. The measurement noise covariance parameters are obviously different from those in the first and third parts.

Conclusions
We propose an optimization approach to smooth the smartphone GNSS position with its velocity information. The following conclusions can be drawn: the optimization method effectively improves the positioning accuracy. In addition, the following directions are recommended and suggested for future work: (1) As stated, the state process and measurement noise are assumed to have a Gaussian distribution. Also, the covariance parameters stay the same while performing the smoothing. Adaptively tuning the covariance matrixes under different conditions might be helpful to improve the position accuracy.
(2) We analyzed the computation load of the KF and optimization methods and found that the computation load of the optimization method is much heavier. It is of great significance to explore how to optimize and reduce the computation load, which is important to extend the optimization method to other applications, especially real-time applications. Also, a larger window size can be utilized, and more past states can be employed in the optimization, and we confirm that a larger window size is adequate to improve and smooth the position results. Here, the window size is 10 for the listed computation load. (3) The last but most important thing is that the opensource optimization method has great potential for many other GNSS signal processing applications, such as carrier tracking, vector tracking, and cooperative navigation. Once the state and measurement models are set up, the optimization method can be implemented. We aimed to present how to understand and implement the optimization method by implementing it in the simple application of smoothing the position with the velocity. Finally, we hope this study and the opensource software inspire exciting applications in GNSS. The MATLAB software is available on the GPS Tool-box website at: https:// geode sy. noaa. gov/ gps-toolb ox/. You can also contact the corresponding author Yuwei Chen (yuwei.chen@nls.fi) to access the source code.
Funding Open Access funding provided by National Land Survey of Finland.
Data availability Dataset and source codes are available from the corresponding author on reasonable request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.