Pedestrian Smartphone Navigation Based on Weighted Graph Factor Optimization Utilizing GPS/BDS Multi-Constellation

Abstract

Many studies have focused on the smartphone-based global navigation satellite system (GNSS) for its portability. However, complex urban environments, such as urban canyons and tunnels, can easily interfere with GNSS signal qualities. Current smartphone-based positioning technologies using the GNSS signal still pose great challenges. Since the last satellite of the BeiDou navigation system (BDS) was successfully launched on 23 June 2020, it is possible to use a low-cost Android device to realize the localization based on the BDS signals worldwide. This research focuses on smartphone-based outdoor pedestrian navigation utilizing the GPS/BDS multi-constellation system. To improve the localization accuracy, we proposed the Weighted Factor Graph Optimization localization model (W-FGO). In this paper, firstly, we evaluate the signal qualities of the BDS via the data collected by the static experiment. Then, we structure the cost function based on the pseudo-range and the time series data for the traditional Factor Graph Optimization (FGO). Finally, we design the weight model based on the signal quality of each satellite and the time fading factor to further improve the localization accuracy of the conventional FGO method. An Android smartphone is utilized to collect the GNSS data for the evaluation and the localization. The experiment results demonstrate the superior performance of the proposed method.

1. Introduction

In recent decades, the demand for location-based services (LBS) has significantly increased, leading to the development of positioning technologies and systems [1]. The global navigation satellite system (GNSS) has been widely adopted for outdoor positioning [2,3,4]. Additionally, with the rapid advancement of mobile internet technology, smart devices have become increasingly important in the field of location-based services [5,6,7]. Since 2016, when Google announced the availability of GNSS raw data for the Android operating system starting from version 7.0, GNSS positioning using smartphones has become a popular research area [8,9].

Traditional satellite navigation systems, such as the Global Positioning System (GPS), Galileo, and GLONASS, have been operating successfully for many years, and have been the subject of numerous studies [10]. The Chinese Bei Dou navigation system (BDS) has emerged in recent years. The BDS progressed from being a demonstration navigation satellite system (BDS-1) to a regional navigation satellite system (BDS-2) by 2012. On 23 June 2020, the last satellite of the third-generation global BeiDou navigation system (BDS-3) was successfully launched. With the global deployment of the BDS constellation, scholars have begun to focus on using this system and its new signals for positioning. However, there are still relatively few studies on smartphone positioning based on the BDS [11]. This paper aims to explore smartphone-based pedestrian positioning using signals from both the BDS and GPS.

The Kalman filter (KF) is an algorithm that can improve the accuracy of observed series measurements over time, even when there are statistical inaccuracies or other sources of noise. The KF generates unknown estimated variables intending to improve the accuracy of the filter results compared with observed measurements. In recent decades, the Kalman Filter has been employed across numerous domains. Scientists have dedicated their efforts to enhance and refine the traditional KF technique. Xia et al. utilized Kalman filter to determine both the yaw misalignment and the velocity error, as well as realized the data fusion between the reduced Inertial Navigation System (R-INS) and the GNSS. These techniques are utilized for the development of intelligent and autonomous vehicles [12,13,14]. A novel approach is also proposed to estimate the Vehicle Sideslip Angle (VSA), which combines data from the GNSS and the Inertial Measurement Unit (IMU). To eliminate the gravity effects caused by vehicle roll and pitch, a square-root cubature Kalman filter (SCKF)-based vehicle attitude angle observer is designed for estimating roll and pitch. The results demonstrated that the proposed method effectively estimateed the VSA in both slalom and double-lane-change (DLC) scenarios. For GNSS positioning, the observations’ state functions and observation functions are usually nonlinear, which cannot be appropriately filtered by KF. The Extended Kalman Filter (EKF) was proposed as a mature solution for this problem. In the EKF, the first-order Taylor Expansion is simply used to linearize the nonlinear system model [15,16]. Despite its benefits, the EKF has a limitation that it only uses data from the adjacent sampling period to calculate position. Besides these data, the historical measurements are also meaningful since the time series data are correlated. To improve the positioning accuracy, we proposed a Weighted Factor Graph Optimization (W-FGO) method, which contains a weighting model and the Factor Graph Optimization (FGO) framework. This method will take advantage of the historical measurements and their correlation. In addition, the weighting model based on the signal qualities of the GNSS data will adjust the weight of different satellites in the cost function. The time fading factor will determine the importance of the data from different epochs. This research’s primary contributions are as follows:

(1)

We conducted the analysis of BDS signal qualities based on the smartphone in Nordic areas. In addition, we utilized the GPS/BDS multi-constellation data to realize pedestrian positioning based on the smartphone.

(2)

We proposed a W-FGO method for GNSS positioning. The W-FGO method consists of two parts: the FGO framework and a weighting model. By utilizing the FGO framework, we explore the influence of time-correlated measurements and states on positioning accuracy. The weighting model is designed based on signal quality and the time fading factor. By utilizing the signal quality, we can adjust the weight of different satellites’ signals, as well as the proportion between the observations and the predicted values. The time fading factor can determine the importance of the data from different epochs.

(3)

We implemented the ground tests with the Huawei Mate40 Pro and the experimental box designed by ourselves. The collected data are processed by the proposed method. The positioning results of the W-FGO method are compared with the least square method (LSM), including BDS-signal-based, GPS-signal-based, and BDS/GPS-multi-signal-based, the EKF method, and the conventional FGO method.

The rest of this paper is organized as follows. Section 2 reviews the related works on the smartphone-based pedestrian positioning utilizing the GPS and BDS systems, as well as studies related to FGO. In Section 3, we analyze the BDS signal qualities based on a static experiment. Section 4 briefly explains the traditional EKF algorithm utilized to realize the GPS/BDS multi-system positioning. In Section 5, we explain the W-FGO algorithm proposed to improve the GPS/BDS-based smartphone positioning. In Section 6, the ground tests of the pedestrian positioning are conducted to compare the performance of the LSM, the traditional EKF, the conventional FGO method, and the W-FGO. The discussion and potential future works are revealed in Section 7. Finally, the conclusion is given.

2. Related Works

2.1. The Smartphone-Based GPS/BDS Multi-System Positioning

GNSS market reports indicate that smartphones now constitute a significant portion of the installed base of GNSS devices. In response to market demands, the principles and algorithms behind smartphone positioning have emerged as common areas of research and discussion within the GNSS field. In recent years, scholars have studied the positioning performance of Single Point Positioning (SPP), Real-Time Kinematic (RTK), and Precise Point Positioning (PPP) using GNSS observations collected by a variety of smartphones and developed some new algorithms [17,18,19]. Lachapelle and Gratton explored the effectiveness of static PPP in modern smartphones, highlighting considerable progress compared to prior smart device models. They illustrated the possibility of achieving coordinate accuracy on the scale of 1 m after gathering data for a duration of 30 min [20]. As a novel satellite navigation system, the BDS positioning on a smartphone is also a research hotspot. Scientists made great efforts in BDS observation analysis. An analysis of the properties of 13 kinds of BDS DCBs and the accuracy of BDS-based GIM was conducted, using data from the International GNSS Service (IGS) and International GNSS Monitoring and Assessment System (iGMAS). The findings, covering a one-month period, reveal that the stability of BDS DCB estimations across distinct frequency bands is connected to the contributing observations. Moreover, the receiver DCB estimations exhibit higher standard deviation values compared to the satellite DCB estimations [21]. Comprehensive PPP models utilizing single-, dual-, triple-, and quad-frequency BDS observations are introduced and assessed by Jin’s group [22]. Chen et al. uses a single-frequency PPP strategy that estimates two clock biases of a smartphone (Xiaomi Mi 8) to achieve a real-time high-precision smartphone positioning [23]. The results indicate that, due to the instability of BDS-locked satellites during dynamic experiments, the positioning accuracy of BDS is inferior to that of GPS. In addition, the GPS/BDS multi-system positioning is also worthy of exploration. Sun et al. explored the relationship between the data quality of smartphones equipped with GNSS modules and the accuracy and reliability of single-frequency RTK positioning [24]. The results show that using multi-GNSS systems data, including the BDS, can effectively improve positioning performance. Additionally, The Smartphone Decimeter Challenge has already been held at the ION in 2021 and 2022. The Google Smartphone Decimeter Challenge (GSDC) is a contest focused on achieving positioning accuracy by utilizing raw GNSS data from smartphones. Smartphone GNSS data exhibit lower signal levels and increased noise in comparison to commercial GNSS receivers, making it challenging to directly apply high-precision positioning methods. The Google Smartphone Decimeter Challenge 2022 (GSDC2022) aimed to advance research in smartphone GNSS positioning accuracy. Dai’s group introduced a global optimization method using gradient descent, accounting for pseudorange, pseudorange rate, accumulated carrier phase (ADR), phone speed, and acceleration constraints at each time epoch on a track. This approach demonstrated superiority over other methods, such as precise point positioning and real-time kinematic. The solution achieved a final score of 1.499 on the private leaderboard, earning second place in GSDC2022 [25]. To sum up, the pedestrian smartphone positioning is a hot topic and worthy of investigation. In our strategy, we utilize the multi-GNSS systems, including the BDS and GPS, to realize the smartphone pedestrian positioning.

2.2. The Factor Graph Optimization

The smartphone GNSS positioning poses great challenges due to the degraded nature of the data in urban environments (e.g., multi-path, poor satellite visibility). The GNSS measurements are highly environmentally dependent and time-correlated [26,27]. However, conventional GNSS positioning methods (e.g., EKF and LSM) cannot simultaneously explore the time correlation among historical measurements and perform poorly in smartphone GNSS positioning. Recently, the FGO is developed for applications with a large number of constraints [28,29]. FGO is well known for its robustness against outliers [30]. Meanwhile, FGO constructs a global cost function to estimate a series of states, comprehensively exploring the correlation among the historical measurements and states [31]. Wen et al. tested the EKF and FGO for GNSS/INS integration in the typical urban scenario in Hong Kong. The results indicate that FGO potentially outperforms the EKF. This research concluded that the effectiveness of FGO can be attributed to two primary factors: (1) FGO employs multiple iterations in the estimation process, resulting in a more robust estimation; and (2) FGO efficiently examines the time correlation between measurements and states, relying on a batch of historical data when the measurements deviate from the Gaussian noise assumption [32]. Ng Hoi-Fung et al. integrated multi-constellation L5-band measurements into 3DMA GNSS to enhance positioning performance in urban canyons, resulting in the L1-L5 3DMA GNSS. Additionally, the study compares various approaches for estimating receiver clock biases in 3DMA GNSS. The integration of different 3DMA GNSS systems is also presented. The FGO method was introduced into 3DMA GNSS to estimate the solution iteratively rather than distributing candidates. Experiments conducted using smartphone data demonstrate that L1-L5 3DMA GNSS provides superior position solutions compared to 3DMA GNSS with only the L1-band, achieving an average positioning accuracy of within 10 m [33]. The 2021 Google Smartphone Decimeter Challenge (GSDC) took place from May to August, 2021. Suzuki’s group proposed a method for estimating a smartphone’s position by using FGO and accumulated delta range (ADR) observations from the smartphone. The incorporation of ADR allows for the estimation of highly accurate relative positions, while precise absolute positions are determined using pseudorange observations corrected with GNSS reference stations as constraints for FGO. By employing the proposed method, they assessed the smartphone’s location and participated in the competition. Their final public score was 2.86 m, securing 2nd place, while the final private score was 1.62 m, earning the 1st place [34]. Thus, it is promising to apply FGO to GNSS-based pedestrian smartphone positioning. Furthermore, the smartphone-based pedestrian positioning without outer assistance is also worthy of research. In addition, based on the traditional FGO framework, there is enormous potential for the exploration and improvement of multi-GNSS-system-based smartphone positioning.

3. The BDS Signal Quality Analysis

To capture data for signal quality analysis, we performed a static experiment using the Huawei Mate40 Pro smartphone and the GEO++RINEX Android app. The experiment involved collecting static data for BDS/GPS/Galileo/GLONASS signals at one-second intervals. The observed data were then saved in RINEX 3.03 format on the phone.

We conducted the static experiment at Aalto University near the seaside. The placement of the smartphone and surrounding environments are shown in Figure 1. During the experiment, we placed the smartphone in an open environment and collected GNSS data for three hours from 10:16:30 to 13:17:00 (UTC time). By processing the GNSS raw data, we were able to extract various features and perform an analysis of the signal quality for BDS.

The Huawei Mate40 Pro, utilized as the test smartphone, was capable of receiving satellite signals from GPS, Galileo, GLONASS, and BDS. During the static experiment, 19 BDS satellites were tracked and their details are presented in

Table 1. The BDS satellites details tracked by Huawei Mate40 Pro.

PRN	Common Name	Int. Sat. ID	Orbit
C02	BDS-2 GEO-6	2012-059A	80.3°E
C05	BDS-2 GEO-5	2010-036A	58.75°E
C08	BDS-2 IGSO-3	2011-013A	117°E
C13	BDS-2 IGSO-6	2016-021A	94°E
C14	BDS-2 MEO-5	2012-050B	between slots B-3 and B-4
C20	BDS-3 MEO-2	2017-068B	Slot B-8
C26	BDS-3 MEO-12	2018-067A	Slot C-2
C27	BDS-3 MEO-7	2018-003A	Slot A-4
C28	BDS-3 MEO-8	2018-003B	Slot A-5
C29	BDS-3 MEO-9	2018-029A	Slot A-2
C30	BDS-3 MEO-10	2018-029B	Slot A-3
C32	BDS-3 MEO-13	2018-072A	Slot B-1
C33	BDS-3 MEO-14	2018-072B	Slot B-3
C36	BDS-3 MEO-17	2018-093A	Slot C-4
C38	BDS-3 IGSO-1	2019-023A	110.5°E
C41	BDS-3 MEO-19	2019-090A	Slot B-2
C42	BDS-3 MEO-20	2019-090B	Moving to Slot B-4
C45	BDS-3 MEO-23	2019-061B	Slot-C3
C46	BDS-3 MEO-24	2019-061A	Slot C-5

. Additionally, Figure 2 intuitively illustrates the sky plot of the observed BDS satellites over the three-hour duration of the experiment, while Figure 3 depicts their visibility. It is evident that the BDS satellites exhibited excellent visibility in Nordic regions.

This section aims to evaluate the quality of static BDS data observations, along with comparative analyses of other constellations’ signals. The signal quality assessment encompasses parameters, such as Signal-to-Noise Ratio (SNR), the number of satellites tracked by the smartphone, and dilution of precision (DOP) [35].

$C/N0$ is a crucial parameter for determining the quality of global navigation satellite system signals, as it represents the normalized SNR, which is the ratio of signal power to noise power density. Vector receivers can use $C/N0$ as a priori information to verify observations and estimate observation noise. Higher $C/N0$ values indicate better GNSS signal quality.

In order to analyze the $C/N0$ of BDS, we computed their mean values and standard deviations (STD). The analysis results are presented in

Table 2. The mean values and standard deviations of $C/N0$ values of different constellations.

Constellation	Mean (dB-Hz)	STD (dB-Hz)
BDS	32.9915	5.7309
GPS	35.3457	6.0886
Galileo	30.0272	6.1583
GLONASS	35.9885	5.8280

. The average $C/N0$ value of BDS was found to be 32.9915 dB-Hz, which is lower than GPS but higher than Galileo. However, this difference can be attributed to inconsistent $C/N0$ measurements among satellites at varying elevation angles. Nonetheless, the STD values of 5.7309 dB-Hz were found to be lower than those of GPS, GLONASS, and Galileo.

However, the C/N0 value for a specific satellite is also influenced by the elevation angle and environmental factors. To further analyze the satellites signals, we also provided the details of each satellite, which is over 20°.

Table 3. The mean values and standard deviations of $C/N0$ values of different satellites.

Constellation	PRN	Elevation Angle (°)	Mean (dB-Hz)	STD (dB-Hz)
	C08	47.9∼22.2	31.5002	4.4961
	C13	54.8∼35.2	37.3910	1.8184
	C27	80.3∼20.0	34.5171	4.9459
	C28	34.8∼20.0	32.9846	6.5131
	C29	20.0∼52.3	35.9722	4.3560
	C30	34.5∼83.4	33.6805	5.2821
BDS	C32	20.0∼43.3	31.3789	5.9250
	C33	20.4∼25.6	31.3565	6.1090
	C36	45.9∼20.0	35.9820	4.0039
	C38	35.4∼20.0	29.2391	4.6902
	C41	31.9∼37.2	29.0289	3.5818
	C45	20.0∼32.5	35.2199	4.6683
	C46	39.5∼20.0	36.1796	5.0112
	G05	33.1∼20.0	36.7451	4.7027
	G08	20.0∼26.9	35.1775	3.8554
	G10	20.0∼41.0	39.3837	2.6044
	G15	20.0∼26.2	37.0806	3.6709
	G16	20.0∼56.4	39.2454	4.7941
	G18	42.1∼75.8	35.3325	3.6929
GPS	G20	20.0∼40.9	32.6049	4.6311
	G23	20.0∼59.6	34.3010	4.8063
	G25	28.5∼20.0	32.4971	4.1494
	G26	20.0∼58.8	40.8628	4.2843
	G29	74.1∼20.0	36.7590	4.3347
	G31	29.6∼20.0	38.9876	1.8416
	E01	32.1∼31.6	23.2340	4.1922
	E07	20.0∼22.9	27.3001	5.4657
	E12	20.0∼75.7	27.9136	4.5362
	E24	31.4∼52.8	34.6085	5.4665
Galileo	E25	20.0∼35.4	30.4673	6.9783
	E26	59.2∼20.1	36.2810	3.4519
	E31	67.3∼20.0	33.4713	4.4128
	E33	57.3∼79.6	29.4240	4.2281
	R01	26.3∼25.0	37.8000	1.1353
	R07	22.0∼20.0	31.9132	4.7945
	R08	21.7∼27.1	33.3159	5.7961
	R09	20.0∼36.6	32.1362	4.3685
	R14	41.3∼20.0	38.6520	1.7602
GLONASS	R15	61.8∼20.0	39.8010	3.7695
	R17	41.1∼84.1	39.1089	3.2452
	R18	20.0∼80.7	38.2912	4.3538
	R19	25.0∼34.5	25.4738	4.5936
	R23	30.5∼20.0	34.5310	2.0161
	R24	79.2∼20.0	39.5584	2.2728

illustrates the relationship between the elevation angles and the C/N0 values. Generally, the probability of high C/N0 values decreases as the elevation angle decreases. Additionally, the numbers of the BDS and GPS satellates with high C/N0 values (over 30 dB-Hz) are more than Galileo and GlONASS, which indicates that BDS and GPS has more satellites suitable for the SPP positioning.

In GNSS navigation and positioning, DOP is employed to assess the impact of the spatial distribution of observed satellites on positioning accuracy. Generally, a better distribution of satellites in the sky results in higher positioning accuracy. Therefore, lower DOP values indicate better satellite geometry and a high probability of achieving high accuracy. DOP is categorized into geometric dilution of precision (GDOP), position dilution of precision (PDOP), horizontal dilution of precision (HDOP), and vertical dilution of precision (VDOP).

Table 4. The average DOP of different constellations during the static experiment.

Constellation	GDOP	PDOP	HDOP	VDOP
GPS	2.4	2.1	1.1	1.7
BDS	1.8	1.6	0.8	1.4
Galileo	2.5	2.3	1.3	1.8
GLONASS	3.1	2.8	1.8	2.1

shows the DOP values for each constellation during the three-hour static experiment. The number of satellites being tracked with elevation angles above 10° for each constellation during the static experiment is shown in Figure 4. BDS has the highest visibility in this scenario, with more than 10 satellites visible for the majority of the time. GPS has over 8 visible satellites for most of the time as well. On the other hand, Galileo and GLONASS have inferior visibility compared to both BDS and GPS. The ascending order of the average GDOP values is as follows: BDS (1.8), GPS (2.4), Galileo (2.5), and GLONASS (3.1). These findings indicate that BDS exhibited the best distribution of satellites in the sky during the static experiment, while GPS came in second.

4. Extended Kalman Filter for GNSS Positioning

The EKF is wildly used for data fusion, which is a mature method for GNSS data filtering. To realize the data fusion of the GPS/BDS signals, the position and velocity are used to realize the one-step prediction of the state vector. The system model of the EKF is described as follows:

$$\begin{array}{ccc}\hfill {\widehat{X}}_{k,k-1}& =& f\left({\widehat{X}}_{k-1}\right)+W_k\hfill \\ \hfill Z_k& =& h\left({\widehat{X}}_{k,k-1}\right)+V_k\hfill \end{array}$$

(1)

where ${\widehat{X}}_k={[x,y,z,v_x,v_y,v_z,ctr]}^T$ is the state vector, which contains the position and the velocity in an ECEF (Earth-Centered Earth-Fixed Coordinate System) coordinate. $tr$ is the receiver clock bias. The state vector is predicted by the nonlinear function $f(·)$ during the adjacent sampling periods. The measurement vector $Z_k$ is the pseudo-range received from the GPS/BDS multi-system, which will be calculated by the nonlinear observation function $h(·)$. The process noise and measurement noise, denoted by $W_k$ and $V_k$, respectively, are modeled as zero-mean Gaussian noise and are associated with covariance matrices. The statistical properties of the noise terms are described as follows:

$$\begin{array}{c}\hfill W_k\sim N(0,Q_k)\\ \hfill V_k\sim N(0,R_k)\end{array}$$

(2)

where the process and measurement noise covariance matrices $Q_k$ and $R_k$ are positive symmetric matrices. The filtering process of EKF is shown as Equation (3):

$$\begin{array}{ccc}\hfill {\widehat{X}}_{k,k-1}& =& f({\widehat{X}}_{k-1})+W_k\hfill \\ \hfill P_{k,k-1}& =& F_k{P}_{k-1}{F}_k^T+Q_k\hfill \\ \hfill K_k& =& P_{k,k-1}{H}_k^T{(H_k{P}_{k,k-1}{H}_k^T+R_k)}^{-1}\hfill \\ \hfill {\tilde{y}}_k& =& y_k-h({\widehat{X}}_{k,k-1})\hfill \\ \hfill {\widehat{X}}_k& =& {\widehat{X}}_{k,k-1}+K_k{\tilde{y}}_k\hfill \\ \hfill P_k& =& (1-K_k{H}_k)P_{k,k-1}\hfill \end{array}$$

(3)

In the EKF method for the GNSS data fusion, $F_k=\frac{\partial f}{\partial X}{|}_{X={\widehat{X}}_{k-1}}=diag[1T;01]$ is the linearized state dynamic matrix of $f(·)$. T denotes the position updating duration. The value of T is one second.

Since we choose the pseudo-range as the observation, the observation matrix $h(·)$ is different from the position-velocity model (PVT) and needs to be linearised. The pseudo-range positioning principle is:

$$\begin{array}{c}\hfill \rho =\sqrt{{(X_s-X_0)}^2+{(Y_s-Y_0)}^2+{(Z_s-Z_0)}^2}+c{\delta }_t\end{array}$$

(4)

where $[X_s,Y_s,Z_s]$ is the position of the satellites in the ECEF coordinate, $[X_0,Y_0,Z_0]$ is the position of the GNSS receiver, and $\rho$ is the pseudo-range observation. Assume that there the number of the satellite-to-receiver pseudo-range measurement is N. The measurement matrix H is:

$$\begin{array}{c}\hfill \mathbf{H}=\left[\begin{array}{ccccccc}\frac{\partial {\rho }_1}{\partial X}& \frac{\partial {\rho }_1}{\partial Y}& \frac{\partial {\rho }_1}{\partial Z}& 0& 0& 0& 1\\ \frac{\partial {\rho }_2}{\partial X}& \frac{\partial {\rho }_2}{\partial Y}& \frac{\partial {\rho }_2}{\partial Z}& 0& 0& 0& 1\\ ⋮& ⋮& \cdots & ⋮& ⋮& ⋮& ⋮\\ \frac{\partial {\rho }_N}{\partial X}& \frac{\partial {\rho }_N}{\partial Y}& \frac{\partial {\rho }_N}{\partial Z}& 0& 0& 0& 1\end{array}\right]\end{array}$$

(5)

The partial derivatives in the above matrix are:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}& \frac{\partial {\rho }_i}{\partial X}=-{\displaystyle \frac{X_i-X_0}{\sqrt{{(X_i-X_0)}^2+{(Y_i-Y_0)}^2+{(Z_i-Z_0)}^2}}}\hfill \\ & \frac{\partial {\rho }_i}{\partial Y}=-{\displaystyle \frac{Y_i-X_0}{\sqrt{{(X_i-X_0)}^2+{(Y_i-Y_0)}^2+{(Z_i-Z_0)}^2}}}\hfill \\ & \frac{\partial {\rho }_i}{\partial Z}=-{\displaystyle \frac{Z_i-X_0}{\sqrt{{(X_i-X_0)}^2+{(Y_i-Y_0)}^2+{(Z_i-Z_0)}^2}}}\hfill \end{array}\right.\end{array}$$

(6)

where $i=1,2\cdots ,N$.

5. Weighted Factor Graph Optimization

To utilize the historical measurements and their inner correlation of GNSS data, FGO has gained much attention as an alternative method. To further improve the positioning accuracy, we improve the conventional FGO framework by introducing an adaptive weighting model.

5.1. The Factor Graph Optimization for Pedestrian Navigation

An overview of the factor graph of GPS/BDS positioning with the constraints of the pseudo-range, the velocity, and the height is shown in Figure 5. Unlike the EKF, the past states are also regarded as unknowns in the optimization method [36]. The state propagation model establishes a strong correlation between these states and measurements.

As shown in Figure 5, the FGO module consists of three kinds of factors. Firstly, the pseudo-ranges obtained from different satellites. Secondly, the location predicated by the previous state based on the velocity. Last but not least, we also use the heights between consecutive epochs as a constraint. The pedestrian navigation addressed in this paper mainly considers the daily walking patterns of pedestrians. As pedestrians’ daily walking typically occurs on flat ground, on the same floor, and in similar situations, there is little significant change in the height over a short period. Even in special situations such as slopes or stairs, the limited walking speed of pedestrians ensures that height information does not change abruptly between adjacent sampling periods. Therefore, this section uses the height information from the adjacent moments in the ENU coordinate system as a constraint to construct the error function, aiming to further enhance the optimization results.

In our strategy, we take the past five state vectors into account and construct the cost function. The cost function can be expressed as follows:

$$\begin{array}{cc}\hfill f(X,\rho ,h)& =argmin(\sum _{i=k-4}^k(X_{i+1}-F_{i,i+1}·X_i)·Q_{i+1}^{-1}·(X_{i+1}-F_{i,i+1}·X_i))\\ & +argmin(\sum _{i=k-4}^k(Z_{i+1}-H_{i+1}·X_{i+1})·R_{i+1}^{-1}·(Z_{i+1}-H_{i+1}·X_{i+1}))\\ & +argmin(\sum _{i=k-4}^k(Heigh{t}_{i+1}-Heigh{t}_i)·I_{i+1}^{-1}·(Heigh{t}_{i+1}-Heigh{t}_i)\end{array}$$

(7)

where $Z_i=({\rho }_1^i,{\rho }_2^i,\cdots ,{\rho }_m^i)$, m is the number of satellites tracked by the smartphone at the epoch i, X means the state vector including the position and velocity, $\rho$ is the pseudo-range, h indicates the height, I means the confidence coefficient matrix, and $H_i$ is the measurement matrix shown in Equation (5).

The FGO aims to find the minimal value of the cost function $f(X,\rho ,h)$. The Levenberg–Marquart (LM) algorithm is utilized to solve the optimal estimations [31,36,37]. The steps of the LM algorithm can be summarized as follows.

Step 1: Expanding the cost function

By expanding the function $f(·)$ and ignoring high order terms, we can obtain Equation (8):

$$\begin{array}{ccc}\hfill \Vert {\Phi }_{i+1}(X_i)-X_{i+1}{\Vert }_{Q_{i+1}^{-1}}^2& \widetilde{=}& (\epsilon ({\tilde{X}}_{i+1})+Jac{b}_{i+1}^{\Phi }·△X)·Q_{i+1}^{-1}·{(\epsilon ({\tilde{X}}_{i+1})+Jac{b}_{i+1}^{\Phi }·△X)}^T\hfill \\ & =& \epsilon ({\tilde{X}}_{i+1})·Q_{i+1}^{-1}·\epsilon ({\tilde{X}}_{i+1})+2·\epsilon ({\tilde{X}}_{i+1})·Jac{b}_{i+1}^{\Phi }·△X\hfill \\ & +& {(△X)}^T·{(Jac{b}_{i+1}^{\Phi })}^T·Q_{i+1}^{-1}·Jac{b}_{i+1}^{\Phi }·△X\hfill \\ \hfill \Vert h_{i+1}(X_i)-Z_{i+1}{\Vert }_{R_{i+1}^{-1}}^2& \widetilde{=}& (\sigma ({\tilde{X}}_{i+1})+Jac{b}_{i+1}^{obs}·△X)·R_{i+1}^{-1}·{(\sigma ({\tilde{X}}_{i+1})+Jac{b}_{i+1}^{obs}·△X)}^T\hfill \\ & =& \sigma ({\tilde{X}}_{i+1})·R_{i+1}^{-1}·\sigma ({\tilde{X}}_{i+1})+2·\sigma ({\tilde{X}}_{i+1})·Jac{b}_{i+1}^{obs}·△X\hfill \\ & +& {(△X)}^T·{(Jac{b}_{i+1}^{obs})}^T·R_{i+1}^{-1}·Jac{b}_{i+1}^{obs}·△X\hfill \\ \hfill \Vert heigh{t}_i-heigh{t}_{i+1}{\Vert }_{I_{i+1}^{-1}}^2& \widetilde{=}& (\varsigma ({\tilde{X}}_{i+1})+Jac{b}_{i+1}^h·△X)·I_{i+1}^{-1}·{(\varsigma ({\tilde{X}}_{i+1})+Jac{b}_{i+1}^h·△X)}^T\hfill \\ & =& \varsigma ({\tilde{X}}_{i+1})·I_{i+1}^{-1}·\varsigma ({\tilde{X}}_{i+1})+2·\varsigma ({\tilde{X}}_{i+1})·Jac{b}_{i+1}^h·△X\hfill \\ & +& {(△X)}^T·{(Jac{b}_{i+1}^h)}^T·I_{i+1}^{-1}·Jac{b}_{i+1}^h·△X\hfill \end{array}$$

(8)

where $Jac{b}_{i+1}$ denotes the Jacob matrix of the cost function. It can be calculated based on Equation (9):

$$\begin{array}{c}\hfill Jacb={\displaystyle \frac{\partial f(X,\rho ,h)}{(\partial X,\partial \rho ,\partial h)}}\end{array}$$

(9)

Equation (8) can be simplified as follows:

$$\begin{array}{c}\hfill f(X+△X)=A+2·B·△X+{(△X)}^T·C·△X\end{array}$$

(10)

Step 2: Deviating the cost function

By differentiating Equation (10) for $△X$, and assuming that the value of the deviated equation equals zero, we can obtain Equation (11):

$$\begin{array}{c}\hfill C·△X=-B\end{array}$$

(11)

Based on Equation (11), we can obtain $△X$ and update $\tilde{X}$:

$$\begin{array}{c}\hfill \tilde{X}=\tilde{X}+△X\end{array}$$

(12)

Step 3: Repeat step 1 and step 2 until the iteration count reaches the predefined threshold or the increment $△X$ reaches a predefined threshold.

5.2. The Adaptive Weighting Model

To further improve the positioning accuracy, we proposed the W-FGO algorithm by introducing a weighting model into the conventional FGO. The weighting model consists of two parts: the $C/N0$-based weighting model, and the adaptive model for the cost function. The $C/N0$-based weighting is used to evaluate the signal importance of each satellite. The adaptive model is used to choose the proportion of different components in the cost function adaptively, as well as to consider the time fading factor.

5.2.1. The $C/N0$ Weighting

Generally, the observation weight matrix $W_{obs}$ can be depicted as [38]:

$$\begin{array}{c}\hfill W_{obs}=diag({\sigma }_1^{-2},{\sigma }_2^{-2},{\sigma }_3^{-2},···,{\sigma }_m^{-2})\end{array}$$

(13)

where m is the number of observations and $\sigma$ is the observation variance.

The $C/N0$ value can be generally used as the criteria to weigh an observation. For example, Figure 6 shows the single-differenced pseudo-range residuals against the $C/N0$ for some selected GNSS satellites. It can be seen that the pseudo-range residuals become larger correspondingly with the decrease of the $C/N0$ values. It denotes that we can use $C/N0$ values to evaluate the quality of the pseudo-range observations.

The $C/N0$-based weighting model applied in this work is given as follows:

$$\begin{array}{c}\hfill {\sigma }_{C/N0}^2={\sigma }_{0,C/N0}^2\times {10}^{{\displaystyle \frac{max(C/N{0}_{max}-C/N0,0)}{10}}}\end{array}$$

(14)

where ${\sigma }_{0,C/N0}$ is the standard deviation of the pseudo-range observations, of which the value is obtained from the numerical analysis and field tests, which is 9.0 m for pseudo-range, $C/N0$ means the current $C/N0$ value of the satellite signal tracked by the smartphone, $max(·)$ is the maximum function, and $C/N{0}_{max}$ is a threshold which is set to 40 dB-Hz.

The observation weight value is the reciprocal of the observation variance. If the observation variance is determined according to Equation (14), the observation weight can be acquired based on Equation (13). Then, we utilized the observation weight $W_{obs}$ to replace $R_{i+1}^{-1}$ in Equation (7).

5.2.2. The Adaptive Weighting for Cost Function

The DOP, which represents the satellite formation, is widely utilized in satellite navigation performance prediction. DOP is derived from the error covariance, which is calculated from the LSM. According to the general definition of DOP, we can obtain Equation (16):

$$\begin{array}{c}\hfill D={(H^T·H)}^{-1}\end{array}$$

(15)

$$\begin{array}{c}\hfill DOP=\sqrt{\sum D_{ii}},i=1,2,3,4\end{array}$$

(16)

where $D_{ii}$ is the $i$ row and the column diagonal element of matrix D. WDOP, proposed by [39], considers the individual error characteristics of the measurements based on their weight, which is more reliable and accurate than DOP. The weighting of DOP can be depicted as Equation (17):

$$\begin{array}{c}\hfill W_{DOP}=diag(w_1,w_2,w_3,···,w_N),w_i={({\displaystyle \frac{{ϵ}_i}{{\sigma }_{UERE}^2}})}^{-1}\end{array}$$

(17)

where N is the index of satellites, ${ϵ}_i$ is the covariance value of pseudo-range observations, which includes both effects of range errors concerning elevation angle and inherent random noise, and ${\sigma }_{UERE}$ means the user equivalent range error (UERE), which is set to a constant value. By applying Equation (17) to Equation (15), we can obtain the weighted least squares equation and calculate the WDOP:

$$\begin{array}{c}\hfill D_{WDOP}={(H^T·W·H)}^{-1}\end{array}$$

(18)

According to Equation (16), we can utilize the diagonal components of matrix $D_{WDOP}$ to compute the WDOP.

WDOP is a performance index that indicates the quantity of positioning errors based on the geometrical deployment of all measurements. As mentioned above, WDOP is used to determine the proportion of the components in the cost function adaptively. According to Equation (7), the cost function of the FGO contains three factors. These three factors can be divided into two parts. One consists of $f(X,F)$ and $f(X,Height)$, which are based on the state differences between the current epoch and the previous epoch. That is, the difference between the predicted state based on the transition matrix F and the current state, as well as the difference between the previous height and the current height. Another one, $f(X,Z)$, is based on pseudo-range observations. That means the difference between the pseudo-range observations obtained from GNSS and the pseudo-range calculated by observation matrix $h(·)$. Introducing the WDOP into the conventional FGO function, the adaptive cost function can be derived as Equation (19). $\beta$ is an empirically determined value. By utilizing the WDOP, we can adjust the proportion between the predicted state $f(X,F)$, the height state $f(X,Height)$ and the observations $f(X,Z)$. The proportion of $f(X,F)$ will increase when the geometric formation of the satellite is poor. In contrast, $f(X,Z)$ will play a more important role in the cost function.

$$\begin{array}{ccc}\hfill f(X,\rho )& =& (1-e^{-WDOP/\beta })·(f(X,F)+f(X,Height))+e^{-WDOP/\beta }·f(X,Z)\hfill \end{array}$$

(19)

Despite the WDOP, since we take the historical data into account, the weight of the data from different epochs should also be considered. Generally, we utilized the time fading factor $\lambda$ to adjust the effect of the data from different epochs. By introduced the time fading factor, Equation (19) can be improved as:

$$\begin{array}{c}\hfill f_i(X,\rho )={\lambda }_i·((1-e^{-WDOP/\beta })·(f_i(X,F)+f_i(X,Height))+e^{-WDOP/\beta }·f_i(X,Z))\end{array}$$

(20)

$$\begin{array}{c}\hfill {\lambda }_i=(1-k_{th}/W_{length}),k_{th}=0,1,2,···,W_{length}-1\end{array}$$

(21)

where $W_{length}$ means the number of the historical epochs considered in the FGO cost function.

6. Experiments and Results

6.1. Ground Tests in Urban Areas

We implement the kinematic pedestrian experiments in urban areas with the smartphone Huawei Mate40 Pro, which is shown in Figure 7b. We utilize the GnssLogger, which is an open-source software announced by Google, to realize the satellite data collection. The experimental box designed by ourselves, which is shown in Figure 7a, contains four parts: the Novatel SPAN-CPT, the antenna, a laptop, and the battery. Novatel SPAN-CPT is a commercial GNSS/INS integration system. The laptop is employed to collect the raw SPAN-CPT measurements. By utilizing the Inertial Explorer (IE) software pronounced by the Novatel company, we can process the Novatel raw measurements and obtain the RTK trajectory as the reference, of which the accuracy can reach the centimeter level, to analyze the positioning accuracy of the proposed method. In the pedestrian experiments, we hold the Huawei Mate40 Pro and carry our experimental box, walking in different scenes in urban areas.

After data collection, we utilize the open-source code pronounced by Google, which goes with the GnssLogger, to process the GNSS raw data and calculate the position, velocity, and time (PVT). This open-source code utilizes the LSM to realize the single point positioning (SPP), based on the single frequency data (1575.42 MHz for GPS and 1561.098 MHz for BDS). In addition, we also used the EKF method to calculate the PVT results as the comparison. For the kinematic pedestrian experiments, we implemented four tests (test 1∼test 4) and collected the different datasets in complex urban areas. The trajectories of test 1∼test 4 contain different scenes, which indicate different levels of the environmental severity. Based on these datasets, we can verify the performance of the proposed method in GNSS-degrade environments. The reference trajectories of these tests are shown in Figure 8. Surrounded buildings and forests will influence the GNSS signal propagation and decrease the smartphone positioning accuracy.

As shown in Figure 8, test 1∼test 4 are conducted in complex urban areas surrounded by trees and buildings. The severity of the positioning environment varies in different tests. test 1 is conducted alongside the road, which is closer to the GNSS-blocked side. In the test 1, around half of the sky is blocked. Test 2∼test 4 is conducted in a community environment. Both sides of the trajectories are surrounded by forests and buildings, especially in the test 4. The features of these tests are also different. Test 2 and test 4 contain more complex and varied scenes. Different sections in the trajectories represent different environments. Some of them are surrounded by trees or buildings, while some of them are not. In contrast, test 1 and test 3 indicate the single scene. Test 1 represents the street scene, which has fewer blockages for the signals. Test 4 represents the alley scene, which is almost full of buildings and trees.

6.2. Results

The statistical analysis results of the horizontal position errors are presented in Table 4. To evaluate the performance of the BDS-based LSM, GPS-based LSM, GPS/BDS-based LSM, EKF, the conventional FGO method, and the proposed method W-FGO, we utilize the mean errors and the standard deviations of the horizontal positions as the evaluation indicators.

According to

Table 5. The horizontal errors for the pedestrian tests.

Dateset	Method	Mean (m)	STD (m)
	BDS-LSM	6.0199	3.5744
	GPS-LSM	5.0221	3.4315
	GPS/BDS LSM	4.3285	3.2680
test 1	EKF	2.8675	2.3140
	FGO	2.5886	1.4571
	W-FGO	1.8729	1.1016
	BDS-LSM	6.1050	3.3940
	GPS-LSM	5.2150	3.9528
	GPS/BDS LSM	3.3169	2.6891
Test 2	EKF	2.4750	2.3051
	FGO	2.0902	1.3024
	W-FGO	1.5704	0.8973
	BDS-LSM	5.7422	3.9868
	GPS-LSM	4.8661	2.4944
	GPS/BDS LSM	4.5547	2.5342
test 3	EKF	2.7858	2.0991
	FGO	2.5551	1.5375
	W-FGO	1.9408	1.5290
	BDS-LSM	6.5474	2.9093
	GPS-LSM	6.2373	2.3262
	GPS/BDS LSM	5.6119	2.0653
test 4	EKF	3.4884	1.5254
	FGO	2.7083	1.3474
	W-FGO	1.8792	0.9530

, as an emerging satellite system, the positioning performance of the BeiDou satellite system is comparable to that of GPS. By utilizing both BDS and GPS satellite signals, single-point positioning accuracy can be further improved. Compared with the BDS LSM, the mean values of the horizontal positioning error of the GPS/BDS LSM decreased by 28.09%, 45.67%, 20.68%, and 14.29% for test 1∼test 4, respectively. The STD values of the horizontal positioning of the GPS/BDS LSM compared with the BDS LSM decreased by 8.57%, 20.77%, 36.44%, and 29.01% for test 1∼test 4, respectively. On the other hand, compared with the GPS LSM, the mean values of the horizontal positioning error of the GPS/BDS LSM decreased by 13.81%, 36.40%, 7.40%, and 10.03% for test 1∼test 4, respectively. The STD values of the horizontal positioning of the GPS/BDS LSM compared with the GPS LSM decreased by 4.76%, 31.97%, almost the same, and 11.22% for test 1∼test 4, respectively. Additionally, it can be seen that the positioning accuracy of the LSM is much worse than the EKF, FGO, and W-FGO. The mean positioning errors and the STD of W-FGO are significantly lower than other algorithms, that is, the filtering performance and the positioning accuracy of W-FGO are the best among these algorithms. Figure 9 and Figure 10 can also intuitively illustrate the positioning performance of the BDS-based LSM, GPS-based LSM, GPS/BDS-based LSM, EKF, the FGO method, and W-FGO.

According to Table 5 and Figure 9 and Figure 10, it can be illustrated that the W-FGO can significantly increase positioning accuracy. Compared with the GPS/BDS LSM, the mean values of the horizontal positioning error of the W-FGO decreased by 56.73%, 52.65%, 57.39%, and 66.51% for test 1∼test 4, respectively. The STD values of the horizontal positioning of the W-FGO compared with the the GPS/BDS LSM, decreased by 68.33%, 66.63%, 39.67%, and 53.86% for test 1∼test 4, respectively. In terms of the EKF, the mean values of the horizontal positioning error of the W-FGO decreased by 34.69%, 36.55%, 30.33%, and 46.1% for test 1∼test 4, respectively. The STD values of the horizontal positioning of the W-FGO decreased by 52.39%, 74.09%, 27.16%, and 37.52% for test 1∼test 4, respectively. Regarding the FGO method, the average horizontal positioning error for W-FGO decreased by 27.65%, 24.97%, 24.04%, and 30.61% for test 1∼test 4, respectively. Additionally, the STD values of the horizontal positioning for W-FGO decreased by 24.40%, 31.10%, almost the same, and 29.27% for test 1∼test 4, respectively.

Excepting the superior positioning performance of W-FGO compared with other methods, the positioning results also demonstrate that data processed by the W-FGO are smoother and less volatile. That is, the W-FGO method enables stable positioning.

The cumulative distribution functions (CDF) of the horizontal positioning errors for test 1∼test 4 are given in Figure 11. It can be obviously demonstrated that the proportions of the W-FGO method with horizontal errors less than 3 m are over 80%, which are significant degradations compared to the FGO, EKF, GPS LSM, BDS LSM, and GPS/BDS LSM.

To further evaluate the positioning performance of the proposed method, the boxplots of horizontal positioning errors for test 1∼test 4 are given in Figure 12. The median horizontal errors of the BDS LSM are 5.7437 m, 5.9088 m, 4.7720 m, and 6.8595 m for test 1∼test 4, while that of the GPS LSM are 4.8565 m, 4.3361 m, 4.9419 m, and 5.9458 m, respectively. The median horizontal errors of the GPS/BDS LSM are 3.0627 m, 2.4849 m, 3.8506 m, and 5.1796 m for test 1∼test 4, while that of the EKF method are 2.7981 m, 1.5499 m, 1.9923 m, and 3.5386 m, respectively. The median horizontal errors of the FGO are 2.4470 m, 1.7958 m, 2.3276 m, and 2.5271 m for test 1∼test 4, respectively. As for the proposed W-FGO method, the median horizontal errors are 1.6221 m, 1.3983 m, 1.4529 m, and 1.8103 m, with improvements of (47.04%, 43.73%, 37.73%, 65.05%), (42.03%, 9.78%, 27.07%,48.84%), and (33.71%, 22.13%, 37.58%, 48.84%) without the outliers for the GPS/BDS LSM, EKF, and the FGO, respectively.

7. Discussion

Note that the iterative initial value of the first epoch is set as the reference value. Except for the first epoch, the iterative initial values are set as follows. For all of the LSM methods, the iterative initial values are set as the iterative result of the previous epoch. For other methods (EKF, FGO, W-FGO), the iterative initial values are set as the combination of the iterative result of the previous epoch and the result of LSM at the same epoch. The experimental results demonstrated the superior performance of the W-FGO method. However, there are still some limitations that need to be developed and further investigated.

(1)

For some extremely severe and varied urban scenes, such as test 4, the positioning accuracy based on GPS/BDS signals is not satisfied. By utilizing the GPS/BDS LSM method pronounced by Google, the mean error of test 4 is even over 5 m. That is, there are still various challenges to the improvement of raw data processing, which will significantly influence the optimization performance of the W-FGO method.

(2)

The LM algorithm is used to solve the nonlinear least squares problem and obtain the position. However, this method leads to a local optimum rather than a globally optimum solution. Thus, the iterative initial value has a significant influence on the accuracy of the results. In this study, we utilize the combination of the iterative result of the previous epoch and the results of the LSM method as the iterative initial value. However, due to the problem of low accuracy of the LSM we mentioned above, the positioning accuracy of the W-FGO cannot be guaranteed in some specific scenes. Obtaining the global optimum solution and enhancing the constraint of the initial value need to be further investigated.

In addition, there are still potential investigation values and prospects in the aspect of smartphone-based pedestrian positioning.

(1)

Since we focus on the investigation of pedestrian positioning, the characteristics of human beings will also make sense. The features of the specific person (e.g., height, step length, and stride frequency) are the potential constraints worthy of research.

(2)

In this study, we pay attention to single-user pedestrian positioning. Collaborative pedestrian positioning is possible due to the information exchange between our smartphones in the future. With the increase of the collaborative network, more constraints can be introduced into the W-FGO, which may be positive for the improvement of the positioning accuracy.

8. Conclusions

Pedestrian navigation based on smartphones plays an increasingly important role in modern life. The positioning accuracy is easily affected by complex urban environments (e.g., multi-path and non-line-of-sight). Since the BDS navigation system has realized worldwide deployment. It is possible to implement pedestrian navigation based on the BDS signals, which has not been thoroughly investigated. In this study, we utilize the GPS/BDS multi-constellation system to realize the smartphone-based pedestrian navigation. An adaptive W-FGO method is proposed to realize the data fusion of GNSS signals and pedestrian positioning. The W-FGO is derivate from the conventional FGO, considering the historical GNSS data and their inner correlation as well as the time fading factor. Firstly, we utilize the $C/N0$-based factor to estimate the signal quality of the observations. Then, the adaptive factor can indicate the proportion of different components in the cost function, as well as consider the influence of time fading. That is, newer data make more sense. We conducted several ground tests with the Huawei Mate40 Pro and simultaneously obtain the RTK results as the reference. The experimental results illustrate that the W-FGO method performs better than other filtering methods. This method can be further improved and extended to more complicated pedestrian navigation systems in our future research.