## Abstract

A smartphone with a highly sensitive antenna receiving numerous unhealthy measurements suffers from non-line-of-sight (NLOS) reception and multipath effects. 3D mapping-aided (3DMA) GNSS has been proven to be effective in urban environments. However, the multipath effect remains challenging for urban positioning. In nature, the new GNSS civilian L5-band signal with a shorter chip length shows a much better resistibility to multipath than the conventional L1-band signal. Therefore, this study integrated the multi-constellation L5-band measurements into 3DMA GNSS to improve the positioning performance in urban canyons, namely the L1-L5 3DMA GNSS. Furthermore, this study compares different approaches on the receiver clock biases estimation for 3DMA GNSS. Finally, the integration of different 3DMA GNSSs is presented. The experiments conducted using smartphone data show that the L1-L5 3DMA GNSS is available for a better position solution than the 3DMA GNSS with L1-band only, thereby achieving a positioning accuracy within 10 m on average.

## 1 INTRODUCTION

Positioning has become essential for our daily routine, and we rely heavily on navigation aid to figure out our current position or the correct direction in which to go. The Global Navigation Satellite System (GNSS) is one of the most popular tools of location-based services (LBS), like smartphone applications.

However, it is a challenge for positioning with the GNSS in urbanized cities. High-rise buildings are often constructed by materials with strong reflectivity (glass, concrete, etc.) that can easily reflect or block GNSS signals.

On the other hand, non-line-of-sight (NLOS) reception and multipath effects are commonly observed in deeper urban canyons, which are more challenging for GNSS positioning. The NLOS reception means that buildings or other obstacles block the line-of-sight (LOS) signals, and only the reflected signals can be received.

In addition, the reflected signal(s) result in longer traveling distance or pseudorange delay. Consequently, there is a positive pseudorange error when NLOS reception occurs. In comparison, the multipath effect occurs at the receipt of both LOS and reflected signals by the receiver (Groves, 2013).

When both direct and reflected signals are fed into the correlator, constructive or destructive interference may have an impact on the measurement output. Therefore, the multipath error can be either positive or negative for the pseudorange measurements, which can deteriorate the positioning performance and result in an error of more than 50 m in an urban environment (Hsu, 2018).

Furthermore, the GNSS receivers equipped for smartphones are also available for measurements with very low carrier-to-noise ratio densities due to the low-cost hardware and baseband processing algorithm designs. These low-quality measurements may not be excluded from their positioning algorithms, which consequently greatly deteriorate the positioning performance as unhealthy measurements.

The multipath effect is identifiable using advanced antenna techniques, such as dual-polarization antennas (Groves et al., 2010; Jiang & Groves, 2014; Sun et al., 2021) and antenna arrays (Kubo et al., 2017; Seco-Granados et al., 2005), though they are not particularly feasible for smartphone applications.

Different types of strategies have been developed to aid positioning in the urban environment, for example, building geometry can be used to identify unhealthy measurements, which is one of the popular approaches.

Various technologies are available for the easy retrieval of the 3D building model by combining the satellite images and airborne LiDAR which provide 2D building contour and building height, respectively (Sohn & Dowman, 2007).

A complete review of the establishment of large-scale 3D building models can be found in the previous reference (Wang, 2013), which makes us quickly generate 3D building models with large coverage areas, and more open-sourced building models are available on the market for free.

The available 3D building model makes it possible to aid urban positioning, namely 3D mapping-aided (3DMA) GNSS (Groves, 2016). The 3DMA GNSS is capable of predicting satellite visibility and transmission path propagation, which be divided into two main categories: *shadow matching* (SDM) 3DMA GNSS (Wang et al., 2013, 2015) and *ranging-based* 3DMA GNSS, such as likelihood-based ranging (LBR; Groves et al., 2020), ray-tracing GNSS (Hsu et al., 2016; Miura et al., 2015), and Skymask 3DMA (SKY; Ng et al., 2020a). These algorithms will be discussed in Section 2.

The nature of different 3DMA GNSS algorithms varies. Specifically, shadow matching can provide a better performance in the across-street direction, while the ranging-based 3DMA GNSS can offer a better performance in the along-street direction as the LOS satellites distribute.

The complementary nature of effective performance in a different direction (Groves et al., 2020) proposed the integration of the shadow matching and ranging-based 3DMA GNSS. As a matter of fact, it can provide a 10 m accuracy on average in the urban environments.

One of the key points in achieving robust ranging is that 3DMA GNSS lies in the estimation of the receiver clock biases. There are two major approaches: 1) the single difference between two satellites from the same constellation, and 2) the popular weight least squares.

Taking advantage of 3D models, in fact, the LOS-visible satellite measurements can be identified and used for estimation. The pros and cons of the two approaches are not compared and discussed in the literature. In view of this, it is taken as one of the objectives of this paper.

It can be seen from the 3DMA GNSS that it is available for the mitigation of the impact of NLOS reception, but not for the multipath effect. Based on the fact that 3DMA GNSS can propagate the reflected path independently, and the assumption that the receiver can only receive this reflected signal for the multipath effect, the receiver receives both the LOS and reflected signals simultaneously. The reflected signals distort the code correlation peak, failing to equalize the power of the early and late channels (Groves, 2013). Both constructive and destructive interference occur in an in-phase reflection with positive ranging error and an out-of-phase reflection with negative ranging errors.

Therefore, an approach is required to mitigate the multipath interference. GNSS L5-band measurements are a practical solution to mitigate the multipath effect compared to the commonly used L1 C/A signal. The design of the L5-band signal leads to a higher chipping rate and shorter chip length compared to the open-service civil L1-band signal (Leclère et al., 2018).

In this case, if the multipath error exceeds one chip length, the amplitude of the autocorrelation between the incoming multipath code signal and local replica would be greatly reduced based on the property of the spreading code (Hegarty & Kaplan, 2005). Thus, the accuracy of the code tracking error discriminated by a commonly used three-channel correlator in terms of the LOS signal could be efficiently improved as its autocorrelation function is unlikely to be distorted by the multipath effect.

As the one chip length of the L5 signal is less than that of the L1 C/A signal, the L5 signal should mitigate a smaller multipath error compared to the latter; the code Doppler frequency can be measured to be closer to the truth. Then, according to the theory for the design of the receiver, both the pseudorange and the carrier-phase measurements can be improved (Dierendonck, 1996).

Taking advantage of the evolution of chip design and manufacturing technologies, the L5-band measurement is already available on low-cost consumer-level receivers. With the development of GNSSs, besides the L1-band signals, smartphones can also receive L5-band signals from different GNSS constellations, such as the L5 signal of GPS/QZSS (ARINC, 2021; Quasi-Zenith, 2021), the E5a signal of Galileo (Galileo, 2021), and the B2a signal of BeiDou-3 (BeiDou, 2017).

This paper employs the L5-band measurements for ranging-based 3DMA GNSS which includes Skymask 3DMA and likelihood-based ranging GNSS. After that, the estimation of the receiver clock bias inside the pseudorange is discussed for ranging-based 3DMA GNSS. Finally, this study proposes the integrated solution of ranging-based 3DMA GNSS and shadow matching.

Dual-frequency GNSS (GPS, Beidou, and Galileo) smartphone-level measurements recorded in Hong Kong urban canyons were utilized for evaluation. Based on the experiment results, we conclude that the L5 measurements can mitigate multipath effects and further improve 3DMA GNSS.

In the meantime, it can be seen from our experiment results that the L5 measurements have an impact on the positioning performance by about 60% epochs with an improvement of more than5monaverage.

To conclude, this paper makes three contributions:

Introduces the multi-constellation L5-band measurements into 3DMA GNSS using a smartphone

Compares the different approaches of the receiver clock estimation for 3DMA GNSS

Proposes the 3DMA GNSS integration method on shadow matching and ranging-based 3DMA, including Skymask 3DMA and likelihood-based ranging 3DMA GNSS

The remainder of this paper is organized as follows: Section 2 details the L1-L5 3DMA GNSS which consists of the shadow matching and the ranging-based 3DMA GNSS, as well as the usage of L5-band measurements; the experiment results and discussion are given in Section 3; and finally, concluding remarks and future work are summarized.

## 2 L1-L5 3DMA GNSS OVERVIEW

The flowchart for the development of 3DMA GNSS is shown in Figure 1.

In most of the 3DMA GNSS algorithms, position hypothesis candidates are distributed around the initial position. In this paper, the weighted least squares (WLS) is selected, and candidates 𝐩𝑗=1…𝐽 are sampled with a 40 m radius and 2 m separations based on pre-generated candidates considering the Skymask.

Consequently, the sampled candidates come together in grid form, which only includes the buildings’ outside. The sampling radius is determined heuristically. LOS/NLOS classification and range correction are performed on each position candidate. The position solution is then finally estimated using weighted arithmetic means and the candidates with their score (determined by means of an individual 3DMA GNSS or their combination).

### 2.1 L1- and L5-band measurements

Due to the limited number of L5-band measurements at the current stage, we only replace the L1-band measurements (pseudorange and carrier-to-noise ratio 𝐶∕𝑁_{0}) with the L5-band once it is available. The pseudorange 𝜌 for satellite 𝑖 (i.e., 𝜌^{𝑖}) is given by:

1

Similar to the carrier-to-noise ratio for satellite 𝑖, (𝐶∕𝑁_{0})^{𝑖} can be expressed as follows:

2

Furthermore, since the chip length of the L5-band signal is ten times shorter than that of the L1-band signal, the multipath errors on the L5-band measurement are also 10 times smaller than that of the L1-band measurement (Ng et al., 2020b) to the maximum extent. We also modify the weighting model (Realini & Reguzzoni, 2013) at a tuning factor of 10 for L5-band measurement (Ng et al., 2020b) which will be further discussed in Section 3.

### 2.2 Shadow matching

*Shadow matching* (Groves, 2011; Wang et al., 2015) determines position by comparing the visibility consistency between the received satellites and the prediction building model. At the same time, shadow matching compares the received satellite and visibility prediction on each positioning hypothesis candidate to locate the position with the highest similarity as the solution (Wang et al., 2012).

In theory, the received satellites are assumed to be LOS satellites, while the non-received ones that appear on ephemeris would have been assumed to be one of the NLOS satellites. Some researchers classify the measurements that are NLOS reception via reflection or diffraction from the received satellite with an intelligent classifier (Sun et al., 2018; Yozevitch & Moshe, 2015).

In Groves et al. (2020), 𝐶∕𝑁_{0} of the received satellite is used to determine the probability to be LOS, thereby mitigating the influence of NLOS reception. The implementation of the shadow matching used is based on the work of Groves et al. (2020) in this study, which is also presented in Appendix A.

Shadow matching outputs the score of 𝑆𝑗,𝑆𝐷𝑀 at candidate *j*, which is used to calculate the positioning solution with Equation (17), or integrate with other 3DMA GNSS as described in Section 3.

### 2.3 Ranging-based 3DMA GNSS

For ranging-based 3DMA GNSS methods, the distribution of the positioning hypothesis candidates is the same as shadow matching. Instead of satellite visibility, the pseudoranges of all the received satellites are simulated for each candidate. Each simulated pseudorange is then compared with the observed pseudorange. Higher scores were given to candidates whose simulated and observed ranges were more similar.

The main difference for the ranging-based 3DMA GNSS methods lies in the estimation of the pseudorange error of NLOS-labeled satellites. In this paper, two state-of-the-art ranging-based 3DMA GNSSs are compared, like Sky-mask 3DMA (Ng et al., 2020a), and likelihood-based ranging (Groves et al., 2020).

The modeled pseudorange for satellite 𝑖 at candidate 𝑗 is notated as where ∗ represents the type of the ranging-based 3DMA GNSS. For example, the subscript SKY stands for Skymask 3DMA, while LBR stands for likelihood-based ranging.

The modeled and measured pseudoranges are used to find the pseudorange difference of a candidate:

3

The pseudorange difference is then used to obtain the weighted geometric mean square error of 𝛿𝜌_{𝑗,∗} and the score for candidate 𝑆_{𝑗,∗}, expressed as follows:

4

5

where 𝐐 refers to a diagonal matrix that contains the uncertainty for each satellite:

6

where 𝜏^{𝑖} represents the measurement uncertainty of the satellite 𝑖.

In this paper, the uncertainty is calculated based on the carrier-to-noise ratio 𝐶∕𝑁_{0}, and elevation angle 𝐸𝑙 (Realini & Reguzzoni, 2013), expressed as follows:

7

where 𝑇, 𝐹, 𝐴, and 𝑎 refer to the empirical constants to control the model performance, respectively. Their values are 𝑇 = 50, 𝐹 = 20, 𝐴 = 50, and 𝑎 = 30, determined historically.

In the case that the corresponding measurement is in L5-band, the weighting tunes have a tuning factor of 10 (Ng et al., 2020b):

8

The multipath effect on the L5-band measurement is 10 times smaller than the L1-band measurement to the maximum extent since the chip length of the L5-band signal is 10 times shorter than that of the L1-band signal. Therefore, the weighting for the L5-band measurement is 10 times larger than that of the L1-band measurement.

In theory, if the tracking algorithms are identical for both L1 and L5 signals, a commonly used non-coherent, early-late power delay lock loop (DLL) discriminator and thermal noise power for the code-based measurements should be positively proportional to the code length (Hegarty & Kaplan, 2005). The code length difference between L1 and L5 presents a relationship in the magnitude of 10 times; therefore, the maximum noise is assumed as a factor of 10.

#### 2.3.1 Pseudorange modeling

On each candidate, a modeled pseudorange is compared with the measured pseudorange . In theory, the candidate near ground truth should obtain a smaller difference between the modeled and measured pseudorange, that is, .

For satellite 𝑖 at candidate 𝑗, can be expressed as follows:

9

where refers to the geometry range calculated using a given satellite’s Earth-centered, Earth-fixed (ECEF) position 𝐩^{𝐢} and candidate ECEF position 𝐩_{𝐣}, given by ; c presents the speed of light; 𝛿𝑡_{𝑗} denotes the time delay of the receiver clock to corresponding satellite constellation; 𝛿𝑡^{𝑖} is the satellite clock delay given by ephemeris data; and refers to the modeled ionospheric delay for the L1-band signal (Klobuchar, 1987).

Due to the different wavelength for L1-band signals and L5-band signals, the value of this error varies. Therefore, a constant for tuning the ionospheric error is required, which is calculated based on the wavelengths of the L1-band 𝜆^{𝑖,𝐿1} and L5-band 𝜆^{𝑖,𝐿5} of a corresponding satellite, that is, (𝜆^{𝑖,𝐿5}∕𝜆^{𝑖,𝐿1})_{2}. of the modeled tropospheric delay (Saastamoinen, 1973).

It is worth noting that the correction provided by tropospheric and ionospheric model scans is replaced by a satellite-based augmentation system (SBAS) message, which can improve the modeling of pseudoranges (Walter et al., 2012). We use the conventional ones for simplicity.

𝜀𝑁𝐿𝑂𝑆 refers to the reflection delay for the NLOS reception, which can be estimated by Skymask 3DMA in Section 3b or likelihood-based ranging GNSS in Section 3c. denotes the thermal noise of the receiver.

In fact, the receiver clock bias 𝛿𝑡_{𝑗}gives a considerable uncertainty to the pseudorange modeling process. Two popular approaches are introduced in the following two sub-sections to deal with it.

##### Single difference (SD) to remove receiver clock bias

Single differencing (SD) of the pseudorange measurements between two satellites from the same constellation can theoretically remove the receiver clock bias from the SD measurement. In the case that the inter-constellation clock offset can be obtained, only one reference satellite is required.

The selection of a reference satellite is important for SD. If there are any errors, especially any errors with the multipath and NLOS in the urban area within the reference satellite, those errors will then be accumulated into the SD measurement.

Generally speaking, the best approach is to select the satellite with the highest elevation angle to avoid the multipath effect. However, in a deep urban environment, the LOS satellite with the highest elevation angle may also suffer from multipath effects.

Therefore, there are two popular approaches for the selection of the reference satellite: First, the reference satellite can be selected using an algorithm based on 𝐶∕𝑁_{0}, elevation angle, and any surrounding buildings (Groves et al., 2020). Only the satellite with the highest score would then be selected and used as the reference satellite for all sampling candidates.

Alternatively, the reference satellite can be selected based on each position’s hypothesized candidate, which therefore can be different across various candidates (Ng et al., 2020a). The performance of each approach is discussed in Section 4.

The SD pseudorange for the 𝑖-th satellite can be expressed as follows:

10

where ∗^{𝑟(𝑖)} refers to the reference satellite for the 𝑖-th satellite.

##### Receiver clock bias estimated by weighted least squares (WLS)

Compared to the SD approach, the general weighted least squares (WLS) approach does not rely on just one reference satellite. However, the error of all the measurements used may have an impact on the resulting estimation.

In Hsu et al. (2016), a common receiver clock bias was estimated and used for all positioning hypothesized candidates. However, in the case that the estimation was erroneous, the simulated pseudorange would not be accurate across all candidates. In other words, range-based 3DMA GNSS would not be robust if this approach is used.

In this paper, we propose estimating the receiver clock delay with WLS based on the geolocation of each candidate, and this WLS only estimates the receiver clock bias. The state vector that includes the receiver clock bias on different constellation 𝐱 can be calculated by the following formula:

11

where the receiver clock delay vector 𝐱_{𝑗} = [𝛿𝑡_{𝑗,𝐺𝑃𝑆},𝛿𝑡_{𝑗,Δ𝐺𝐴𝐿},𝛿𝑡_{𝑗,Δ𝐺𝐿𝑂},𝛿𝑡_{𝑗,Δ𝐵𝐷𝑆}]^{T}, and 𝛿𝑡_{𝑗,𝐺𝑃𝑆} represents the GPS receiver clock delay at candidate 𝑗; 𝛿𝑡_{𝑗,Δ𝐺𝐴𝐿}, 𝛿𝑡_{𝑗,Δ𝐺𝐿𝑂}, and 𝛿𝑡_{𝑗,Δ𝐵𝐷𝑆} denote the inter-constellation delay between GPS, Galileo, GLONASS, and Beidou, respectively.

It should be noted that we only use the predicted LOS satellites on each candidate and its pseudorange to estimate the receiver clock delay. 𝐀 refers to the design matrix consisting of the flag of the available satellites’ constellation:

12

where c refers to the speed of light; represents the flag to decide whether the satellite of 𝑖 belongs to ∗ constellation or not. In the case that satellite 𝑖 belongs to ∗ constellation, the value of one will be assigned; otherwise, zero will be assigned.

It is worth noting that if the inter-constellation error can be estimated, the satellites will be considered as one single constellation. 𝐀 is a single column vector with all speed of light of c, where the size will be the same as or .

This paper also compares the positioning result on estimating the receiver clock delay with inter-constellation correction.

#### 2.3.2 Skymask 3DMA GNSS

Skymask 3DMA (Ng et al., 2020a) is a simplified version of a ray-tracing GNSS (Hsu et al., 2016; Miura et al., 2015) which reduces its computational load while maintaining an almost identical accuracy.

Skymask 3DMA utilizes an enhanced skymask, which also consists of the building height information and azimuth angle of the reflecting planes (AARP) associated with each azimuth angle. Consequently, the reflecting point can be detected using the AARP value and building height information.

In the beginning, Skymask 3DMA required the NLOS satellite’s azimuth and elevation angles of and , which enhanced the skymask of the candidate of 𝑆𝑀_{𝑗}. Then, the reflecting point azimuth and elevation angles could be obtained by using the AARP value.

With the building height information on the corresponding azimuth and elevation angle of the reflecting point, the position of the reflecting point can theoretically be obtained as well. Therefore, the reflection delay of can be calculated by subtracting the distance between satellite and candidate from the total distance of the reflection path.

Therefore, the modeled pseudorange of Skymask 3DMA can be obtained, and it is worth noting that only the LOS and NLOS with reflection found satellites are going to score the candidate (i.e., {LOS, NLOS with reflection}).

The NLOS reflection delay can be obtained by Equation (2), notated as in Ng et al. (2020a) after locating the reflection point which is determined by the enhanced skymask at the candidate and AARP of the surrounding potential reflector.

The pseudorange difference of Δ𝝆_{𝑗,𝑆𝐾𝑌}, weighted root mean square error of 𝛿𝜌_{𝑗,SKY}, and the score for the candidate of 𝑆_{𝑗,𝑆𝐾𝑌} can be obtained by Equations (4) and (5), respectively.

#### 2.3.3 Likelihood-based ranging

Unlike the Skymask 3DMA or ray-tracing GNSS that correct the NLOS delay with building geometry, likelihood-based ranging (LBR) corrects the NLOS delay using a statistical approach.

First, a satellite is selected as the reference satellite to obtain the SD measurement for LBR based on the scoring algorithm (Groves et al., 2020). The implementation of the LBR used is based on the work of Groves et al. (2020).

At each positioning candidate, the satellite visibility will firstly be estimated by the candidate’s skymask. For the NLOS-classified satellite at a candidate, the pseudorange difference will be remapped.

For NLOS-labeled pseudoranges, the only unknown term is the in Equation (9). Therefore, the modeled pseudorange with NLOS error is given by the following formula:

13

And the pseudorange difference should be dominated by the NLOS error at this stage:

14

Consequently, the LBR remaps the NLOS pseudorange difference from to , and mitigates the NLOS error of . LBR first estimates the skew-normal distribution for NLOS measurements before substituting the NLOS pseudorange difference to the skew-normal distribution and obtaining the cumulative probability.

Then, this cumulative probability is fed back to the direct LOS cumulative probability function to get the corresponding pseudorange difference which is assumed to be NLOS error-free. The full algorithm can be found in Appendix A.

All pseudorange differences at a candidate Δ𝝆_{𝑗,𝐿𝐵𝑅} are obtained and employed to calculate the weighted root mean square error of 𝛿𝜌_{𝑗,𝐿𝐵𝑅} and score for the candidate 𝑆_{𝑗,𝐿𝐵𝑅} with Equations (4) and (5), respectively.

#### 2.3.4 Comparison of the ranging-based 3DMA GNSS methods

A comparison between two ranging-based 3DMA GNSS is given in Table 1.

The main difference between the two ranging-based methods used in this paper (i.e., Skymask 3DMA and likelihood-based ranging GNSS) presents the theory behind NLOS correction. The Skymask 3DMA based on geometry using the building model generated the enhanced skymask to propagate the NLOS reflection path, thereby estimating the NLOS reflection delay, which is similar to ray-tracing GNSS.

The uncertainty of Skymask 3DMA is mainly sourced from the provided building model for an enhanced sky-mask generation, and the unmodeled building structure has an impact on the transmission path propagation, resulting in the fact that some of the NLOS-labeled satellites are not available for candidate scoring if no reflecting point can be found.

The likelihood-based ranging GNSS employs the statistical model for the correction of the NLOS error. Therefore, the selected reference satellite becomes one of the main uncertainties of this algorithm. In the case that the multipath effect is observed on the master satellite, the NLOS error may be wrongly estimated for the NLOS measurement, which will deteriorate positioning performance.

However, unlike Skymask 3DMA, in theory, likelihood-based ranging GNSS is able to correct all NLOS-labeled satellites at a candidate. Therefore, based on the nature of LBR and SKY, the integration of the two ranging-based 3DMA GNSS methods is as follows.

Firstly, SKY is used if a reflection path can be found. Otherwise, the LBR will be used. In theory, this integration strategy can mitigate the error of the master satellite, contributing to all NLOS satellites while using all NLOS measurements for positioning.

### 2.4 Position solution and integrated solution

In this paper, we integrated the two ranging-based 3DMA GNSS methods and shadow matching. For the integration of the two ranging-based 3DMA GNSS algorithms, the simulated pseudorange is obtained by the following formula:

15

It should be noted that for the LOS-labeled satellite, the pseudorange difference is identical for both SKY and LBR at the candidate as no NLOS delay correction is required. Therefore, the pseudorange difference of integrated ranging-based 3DMA at a candidate Δ𝝆_{𝑗,𝑅𝑁𝐺} is obtained by Equation (3).

The weighted root means square error of 𝛿𝜌_{𝑗,RNG} and the score for the candidate of 𝑆_{𝑗,𝑅𝑁𝐺} can be calculated by using Equations (4) and (5), respectively.

Followed by the integration of shadow matching and ranging-based 3DMA algorithms, the candidate score is given by the following formula:

16

The position solution is calculated by the weighted average of the scored positioning candidates. And the score is given by corresponding algorithm or integration method, ∗:

17

where 𝐽 refers to the total number of the position-hypothesized candidates. The position solution calculations are identical for all 3DMA GNSS positioning methods, including the shadow matching in Section 2.2, the ranging-based 3DMA GNSS in Sections 2.3.2 and 2.3.3, and the integrated solution in Section 2.4.

## 3 DESIGNED EXPERIMENT RESULTS

### 3.1 Experiment Setup

Several designed experiments were conducted in urban canyons in Hong Kong using a Xiaomi Mi 8 smartphone that supports the L1- and L5-band measurements. The constellations used were GPS, Galileo, GLONASS, and Beidou, with additional experiment information summarized in Figure 2 and Table 2.

The ground truth is labeled manually based on Google Earth (with a working accuracy of 1 meter in our experience). The ratio of building height to street width is calculated by (𝑏𝑢𝑖𝑙𝑑𝑖𝑛𝑔_ℎ𝑒𝑖𝑔ℎ𝑡)∕(𝑠𝑡𝑟𝑒𝑒𝑡_𝑤𝑖𝑑𝑡ℎ), and if the value is higher, the street is narrower with taller buildings surrounded, and the environment is more challenging for GNSS positioning.

There are seven 3DMA GNSS methods evaluated in total, and the summarized information is given in Table 3. The score mentioned above can be utilized to calculate the positioning by replacing 𝑆_{𝑗,∗} in Equation (17).

The analysis of the pseudorange quality of the L1 and L5 measurements used in 3DMA GNSS is given in Section 3.2. The performance of the approaches on the estimation of receiver clock bias is discussed in Section 3.3, and the evaluation of the 3DMA GNSS positioning results is shown in Section 3.4.

The criteria for the evaluation of the positioning performance include root-mean-square (RMS) error, mean value, standard deviation (STD), maximum value (MAX), and minimum value (MIN).

### 3.2 Pseudorange quality

Experiment 3 is used to evaluate the pseudorange difference between the pseudorange measurement and simulated pseudorange with NLOS correction based on the ground truth (GT) location of Δ𝝆_{𝐺𝑇,SKY}, that is, (3) , and the result is shown in Figure 3. Only the received satellites with both L1- and L5-band measurements are shown in the figure.

In Figure 3, the upper and lower rows represent the L1- and L5-band pseudorange difference for each satellite, respectively. The first entry for “weighted AVG” represents the weighted average pseudorange difference of all available satellites. The upper row refers to the weighted average L1-band pseudorange difference, while the lower one defines the L5-band measurements replacing the L1-band ones.

The receiver clock bias is estimated using the weighted least-squares with inter-constellation correction, as shown in Equation (11). Therefore, 𝐀 represents a single column vector with c. The reference station gives the interconstellation correction, and the measurement from the reference station SatRef by the Lands Department of Hong Kong is employed.

It can be seen from Figure 3 that the pseudorange difference using an L5-band measurement is smaller than that of an L1-band, especially with a satellite with a low elevation angle, like satellites G03 and G08. The ray-tracing simulation is also performed on the ground truth location, as shown in Figure 4.

The ray-tracing simulation indicates that the satellite of G03 is affected by multipath. In fact, both G30 and G01 are similar (but the single reflection cannot be identified through ray-tracing) to G03. Thus, we illustrate the result mainly based on G03. The pseudorange error calculated based on the double-differencing method is shown in Figure 5(a) which was introduced in Section 2.3.3. of our previous work (Xu et al., 2019).

As mentioned, the L5 signal is at a higher chipping rate of 10.23 Mchip/s, and a shorter wavelength is available for the reduction of the waveform distortion. Therefore, the L5 signal can achieve a higher code measurement accuracy when both L1 and L5 correlator share the same early-late spacing.

In particular, for positioning in an environment featuring many reflective obstacles such as Experiment 3, the L5 can effectively reduce the multipath effect to provide a better positioning accuracy.

The average pseudorange error of L1 is 6.32 m with a standard deviation of 15.53 m, while the average error of L5 is 1.46 m with a standard deviation of 1.86 m. The maximum difference between L1 and L5 pseudorange error is at epoch 970 with 69.49 m, where the pseudorange error of L1 and L5 are 71.57 m and 2.08 m, respectively.

In the entire experiment, the curves of the L5 pseudorange error fluctuate more slightly than the one related to L1, showing that the L5 pseudorange quality performs better in the experiment, which is beneficial to the positioning.

Seen from the probability density distribution of L1 and L5 of G03 in Figure 5(b), the L5 signal shows a higher precision on the pseudorange than the L1. As shown in Figure 6, the cumulative distribution function (CDF) of L1- and L5-band pseudorange error based on the data we collected in all the experiments also supports this point of view.

The 1-sigma pseudorange error of L5- and L1-band are approximately 11 m and 20 m, respectively. If taking a comprehensive consideration of the data distribution, 90% of the L1 signal is within 58.3 m, while the L5 signal is within 26.8 m. In this manner, we experimentally show that the L5 range measurement is less affected by multipath compared with that of L1.

For the ranging-based 3DMA GNSS, this is beneficial to the improvement of the concentration of the likelihood distribution of candidates due to the high consistency between the modeled and the observed pseudorange.

On the other hand, the G08 satellite is an NLOS diffracted signal, shown on the right side of Figure 4, which explains why we observe a large pseudorange error for G08. It is worth noting that we only assume that all NLOS receptions are reflected signals, and only apply reflection correction to these satellites in this paper.

Diffraction and reflection delays are different, and different strategies should be applied separately to these effects (Zhang & Hsu, 2021). This diffraction path of G08 passes through the narrow gap between two buildings, and multiple reflections may occur on these two surfaces.

Lastly, G08 is a satellite with a low elevation angle of about 15^{◦}, and the average 𝐶∕𝑁_{0} for L1- and L5-band measurements are 12.8 dB-Hz and 10.6 dB-Hz, respectively. Due to the low-cost design approaches, the commercial receiver is usually designed to have the risk of receiving these complicated reflected/diffracted signals.

The range measurement with very weak 𝐶∕𝑁_{0} can partially reflect the actual reflection interference experienced by the signal. Therefore, we cannot conclude that the L5 measurement of G08 has a higher quality.

### 3.3 Evaluation of receiver clock bias estimation methods

The estimation methods compared are:

**SD-IC**: The single difference (SD) method introduced in Section 2.3.1.1. in which the inter-constellation (IC) offset is estimated, all available satellites are considered as one system. The reference satellite is selected by using a scoring algorithm [(14) in Groves et al. (2020)] and used for all candidates, same as Groves et al. (2020) proposed.**P-SD-IC**: The SD method introduced in Section 2.3.1.1. in which the inter-constellation (IC) offset is estimated; all available satellites are considered as one system. The reference satellite is selected for each candidate with the highest elevation angle (Ng et al., 2020a).**WLS-IC**: The WLS method introduced in Section 2.3.1.2. in which the inter-constellation (IC) offset is estimated; all available satellites are considered as one system. All available satellites are used to calculate the receiver clock error (Hsu et al., 2016).**P-WLS-IC**: The proposed particle-based WLS (P-WLS) introduced in Section 2.3.1.2. in which the interconstellation (IC) offset is estimated; all available satellites are considered as one system. Only LOS satellites are used to estimate the receiver clock error. The receiver clock offset is calculated by Equation (11). IC clock offsets are used, hence, 𝐀 refers to a single column vector with c.**P-WLS**: The proposed particle-based WLS (P-WLS) introduced in Section 2.3.1.2; only LOS satellites are used to estimate the receiver clock error. The receiver clock bias is estimated using Equations (11) and (12).

The clock offset of inter-constellations is very useful to improve the estimation accuracy of the receiver clock bias. In particular, when only receiving a few satellites for a specific constellation with a low elevation angle, an erroneous receiver clock bias of that constellation may be estimated.

In this paper, the inter-constellation offset (IC) is estimated using the measurements from the reference station of SatRef built by the Lands Department of Hong Kong. Table 4 shows the RMS positioning error of Experiment 3 on different receiver clock error estimation methods with L1- and L5-band measurements.

Method 3, WLS-IC, is available for the best positioning result for SKY, which is followed by Method 4 with an extra 0.46 m RMSE in which one can obtain the best positioning results for LBR. For SKY+LBR, Method 1 achieves the best positioning result, which is followed by Method 4 with an extra 0.1 m RMSE.

In terms of the average RMSE of SKY, LBR, and SKY+LBR, Methods 1, 3, and 4 are within 12 m. In general, the WLS is available for better positioning results compared with SD. These results also show that the IC clock offset is available for a better estimation of the receiver error and positioning.

The RMSE of Methods 1 and 4 are similar, but the STD and MAX behave differently, as shown in Table 5. The STD of SD is about 5 m, while the WLS value is about 4 m on SKY, indicating that the WLS can provide a more stable positioning performance. This is because the SD only employs one satellite to eliminate the receiver clock error by means of differencing across all other satellites.

In the case that an error exists in the selected reference satellite, the error will also accumulate on other satellites. However, WLS estimates the clock error with all LOS satellites, which indicates that the error within one satellite will be mitigated with the presence of other healthy measurements.

Furthermore, compared with SD on LBR, the STD error is found to be smaller by 2 m on WLS. A similar difference is found in MAX where the SKY with SD is about 57 m, while the WLS is 42 m.

Seen from the positioning heatmap, a larger positioning error can be found in SD at several epochs. We selected an epoch from Experiment 3 for investigation, as shown in Figure 7(a)–7(f). Besides, the skymask with a 𝐶∕𝑁_{0} value of each satellite is shown in Figure 7(g).

It can be observed from Figure 7 that higher score particles in SD heatmaps from (a)-(c) are more concentrated to their solution compared with the particle-based WLS ones in (d)-(f). Generally, the size of the candidate of SD in light blue to red is smaller than the WLS one, and the major axis length is about 10 m and 20 m for SD and WLS, respectively.

In general, this is a sign indicating a stable and accurate 3DMA GNSS solution (Ng et al., 2021). However, high-score particles drift off from the truth location compared to the WLS one. This error may come from the reference satellite when performing the SD.

In this paper, G01 is selected as the master satellite. Compared with other satellites, G01 is not near the building boundary and has a relatively high 𝐶∕𝑁_{0} value. However, G01 shows an elevation angle of more than 60^{◦}, but its 𝐶∕𝑁_{0} is only about 35 dB-Hz.

Although the ray-tracing for G01 is not available to find any reflected path, there could be some multipath with destructive interference occurred to the signal, resulting in an estimation error for Method 1. Moreover, this satellite doesn’t have a very high elevation angle. Another evaluation using the data obtained from Experiment 4 with an identical conclusion is included in Appendix B.1 as supplementary information.

Figure 8 shows the CDF on how the receiver clock error estimation methods contribute to the positioning error. It can be observed from the result that Methods 1, 3, and 4 outperform the others on most of the ranging-based 3DMA algorithms. With further investigation on the zoom-in plot [subfigure (d) in Figure 8], Methods 3 and 4 involving the WLS method with IC clock offset correction perform better in the position domain.

So far, short conclusions for Sections 3.2 and 3.3 can be made as follows:

The integration of L5-band measurements is available for the improvement of 3DMA GNSS positioning performance.

The utilization of IC clock corrections in the estimation of the receiver clock bias is beneficial to ranging-based 3DMA GNSS.

Method 4, as the proposed particle-based WLS approach, is more robust compared to other estimation methods on receiver clock bias.

### 3.4 Positioning results

Based on the short conclusions drawn above, the L5-band measurements plus the proposed P-WLS-IC (Method 4) should be able to obtain better performance. Therefore, the positioning results of 3DMA GNSS methods are evaluated based on this configuration. The positioning statistics of all experiments are summarized in Table 6, and a graphical view is illustrated in Figure 9. Besides, all different positioning result statistics are presented in Appendices B.2 and B.3.

It is found from the positioning results that the positioning RMSE of the proposed integrated solution (SKY+LBR+SDM) is within 10 m in a different scenario. Even in very deep urban canyons (i.e., Experiments 4 and 5), an integrated solution can still achieve the RMSE of 5 m. It is also observed that the algorithm integrated solution SKY+LBR+SDM outperforms the others in Experiments 1, 3, and 5 in terms of RMSE. While in Experiment 2 and 4, LBR+SDM obtains a slightly better positioning result compared to the proposed SKY+LBR+SDM.

In Experiment 5, the RMSE of positioning is about 8 m, 9 m, and 23 m for SKY, LBR, and SDM standalone positioning solution, respectively. After integration, both SKY+LBR and SKY+LBR+SDM RMSE are within 5m.

To better demonstrate the effectiveness of integration of SKY and LBR, we show a positioning heatmap of one of the epochs in Figure 10. The high-scoring particles of the SKY, LBR, and SDM have drifted away from the truth location toward the intersection center. After integration, only the common high score candidate remains.

The combination of SKY and LBR is shown in Figure 10(d). Although the reddish area of LBR is large, only the candidates near the truth location keep the high score after combining with the SKY. The positioning result is 4.7 m, which is a huge improvement compared with SKY and LBR alone.

Furthermore, the solution of SDM locates on the opposite side of the street with an error of 28.3 m. After the integration with the ranging-based 3DMA algorithm, the positioning result is decreased 4.6 m, as shown in Figure 10(e). The integration process based on different algorithms can keep the commonly high score candidate, and improve the positioning results.

However, an integrated solution can get worse if the standalone solution drifts away sometimes. Taking an epoch in Experiment 1 for example, as shown in Figure 11, we can observe that SKY heatmap in 11(a) has drifted away from the ground truth, causing the positioning error to exceed 20 m.

This is due to the fact that the SKY only uses the corrected NLOS satellites for positioning. However, there is no correction found for the NLOS satellites in this case. SKY can only use the LOS satellite in this paper for positioning.

While the LBR heatmap in 11(b) concentrates on the ground truth, the excellent performance of the LBR helps it successfully obtain a good positioning result when integrating with SKY [i.e., SKY+LBR in Figure 11(d)].

However, when integrating with SDM, the heatmap becomes more discrete, which is due to the fact that SDM heatmap in Figure 11(b) is scattered along the building edge. Consequently, although the positioning result of SKY+LBR+SDM reaches 2.7 m, we can see that the heatmap high score area is larger than that of the SKY+LBR one.

It can be seen from the heatmap positioning result that integration is not always available for a good result, or sometimes can add uncertainty to the integrated solution, which is also observed in Experiment 2. As the experiment is conducted at a single side building, no correction can be found for the NLOS satellites. As the result, the SKY can only use the LOS satellite for positioning.

The positioning RMSE achieved is approximately near 40 m. This worst positioning performance also affects the integration between other algorithms. It can be observed that with the integration between SKY and LBR, the positioning RMSE becomes larger than that of LBR itself.

Even if it combines with SDM, that is, SKY+LBR+SDM, the positioning RMSE is about 5.83 m. In comparison, LBR+SDM RMSE is 5.47 m which is about 0.4 m better, showing that the uncertainty is added by SKY when integrating its solution with others.

## 4 CONCLUSIONS AND FUTURE WORK

This paper integrated the L5-band measurements to the ranging-based 3DMA GNSS algorithms, and estimated the receiver clock bias using a particle-based WLS method. We integrated shadow matching and two ranging-based 3DMA GNSSs (i.e., the likelihood-based ranging and skymask-based approaches).

According to the experiments conducted in Hong Kong urban canyons, the L1-L5 3DMA GNSS is able to provide more stable positioning performance with an accuracy within 10 m on average. The L5-band measurements can reduce the pseudorange difference through a noticeable improvement in the positioning accuracy. With new smartphones starting to support the L5-band signal, it will be the future development in urban positioning.

One interesting aspect to improve the 3DMA GNSS is the generalization on the several heuristic models used. One suggested theme of future work is to determine the parameters of these heuristic models based on quantitative features that were extracted from the surrounding building environment.

In terms of the computational load, currently, the 3DMA GNSS is still processed in positioning candidate distribution. The positioning performance relies on the accuracy of initial positioning to distribute the candidates, and the distributed candidates have to cover the truth location to obtain the best solution.

In theory, the distributed candidates far from the truth location are unnecessary computations for positioning. Therefore, an optimal method is required to estimate the solution. Shortly, we will try to introduce factor graph optimization (FGO) into 3DMA GNSS to estimate the solution iteratively rather than distributing candidates.

Nevertheless, for positioning in grid-distributed candidate approaches, some candidates far from the truth location can still achieve a high score. These local minima issues can deteriorate the positioning performance. Therefore, the reliability examination and unreliable candidate exclusion may be required later.

In terms of accuracy, the integrated solution of 3DMA GNSS can achieve a positioning accuracy of 5 m in an urban environment. However, this is still meter-level accuracy positioning. If we want to enable the real-time-kinematic (RTK) GNSS, the integrated solution can be utilized as the initial position and visibility estimation.

## HOW TO CITE THIS ARTICLE

Ng H-F, Zhang G, Luo Y, Hsu L-T. Urban positioning: 3D mapping-aided GNSS using dual-frequency pseudorange measurements from smartphones. *NAVIGATION*. 2021;*68:*727–749. https://doi.org/10.1002/navi.448

## ACKNOWLEDGMENT

The authors would like to thank Dr. Taro Suzuki from the Chiba Institute of Technology in Japan for providing an algorithm to obtain the sky-pointing fisheye image from Google Earth Studio.

## APPENDIX A: 3DMA ALGORITHM IMPLEMENTATION

### A.1 Shadow matching (SDM)

The implementation of the shadow matching used is based on the work of Groves et al. (2020). After the construction of the position hypothesis candidates, the satellite visibility is estimated for each satellite on each candidate’s skymask. And the probability of satellite 𝑖 at candidate 𝑗 is predicted to be direct LOS using skymask, , where 𝐵𝐵 stands for building boundary, is given by the following formula:

A)

The applied values in this paper are determined empirically, which account for the prediction uncertainties, such as 3D building model accuracy or signal diffraction by obstacles not included in the model. The probability of each signal to be directed to LOS is predicted by using 𝐶∕𝑁_{0} measurement, so 𝑝(𝐿𝑂𝑆|𝐶∕𝑁_{0})^{𝑖} is given by the following formula:

A)

where p_{0‒min}, s_{min}, p_{0‒max}, s_{max}, a_{0}, a_{1}, and a_{2} refer to the empirically determined constant, and we set p_{0‒min} = 0.26 dB-Hz, s_{min} = 22 dB-Hz, p_{0‒max} = 0.9 dB-Hz, s_{max} = 32 dB-Hz, a_{0} = ‒2.252 dB-Hz, a_{1} = 0.1492 dB-Hz, and a_{2} = ‒0.001588 dB-Hz. These constants are determined based on our private communication with Dr. Paul Groves and experience. Therefore, the score of the matching between the predicted and measured satellite visibility can be calculated as:

A3

Finally, the score of shadow matching on a candidate can be obtained by the geometric mean of:

A4

where 𝐼 refers to the total number of the satellite. The geometric mean of all satellite matching is to normalize the score at each candidate. According to our experience, this is beneficial to the integration process. Heatmaps will not concentrate after normalizing the score. Consequently, when integrating with ranging-based 3DMA, shadow matching results will not dominate other methods.

### A.2 Likelihood-based ranging (LBR)

The implementation of likelihood-based ranging follows (Groves et al., 2020). LBR first selects a reference satellite with the scoring algorithm. Then, it calculates the standard deviation of all errors except for the NLOS path delay with 𝐶∕𝑁_{0}:

A5

where a and b refer to the empirically determined constant, and we set 205,700.0 m^{2} and 18.72 m^{2} in this paper, respectively.

At each positioning candidate, the satellite visibility will firstly be estimated using skymask. For the NLOS-classified satellite at a candidate, the pseudorange difference will be remapped. The NLOS pseudorange difference of is mapped using a skew-normal distribution to determine the cumulative probability of 𝐹_{𝑁𝐿𝑂𝑆}.

The estimation process first gets the parameters for the skew-normal distribution which is defined by three main parameters (i.e., location 𝜉, scale 𝜔, and shape 𝛼). These parameters are obtained by using the following formula:

A6

A7

A8

where 𝜇_{N} and 𝜎_{N} refer to the empirical constant on the mean and standard deviation of the NLOS path delay, respectively. We set 𝜇_{N} = 42 m and 𝜎_{N} = 36 m in this paper. 𝜇_{𝐿}represents the mean of LOS pseudorange difference, and 𝜎_{r} denotes the error standard deviation of the reference satellite. The cumulative probability of NLOS pseudorange difference 𝐹_{𝑁𝐿𝑂𝑆} can be obtained by substituting the values from Equation (A6) into Equation (A8):

A9

where integral of the normal distribution (also known as Gauss error function) erf(𝑥) (Groves et al., 2020) and Owen’s function T can be calculated by the following formula:

A10

A11

Finally, NLOS measurement can be remapped to the corresponding direct LOS error distribution of zero-mean Gaussian distribution by CDF. The remapped pseudorange difference can be obtained from 𝐹_{𝑁𝐿𝑂𝑆} by solving:

A12

where 𝐹_{𝑁𝐿𝑂𝑆} at the left-hand side refers to the NLOS CDF, while the right-hand side includes the corresponding LOS CDF, which makes the LOS and NLOS measurements comparable in the same hypothesis for positioning.

## APPENDIX B: DETAILED EXPERIMENT RESULTS

This appendix provides more details on the experiment results or the evaluation not included in the main content.

### B.1 Receiver clock bias estimation methods evaluation results on Experiment 4

To show the repeatability of the investigation described, we also compare different receiver clock estimation methods in Experiment 4 with the results shown in Table B1.

As shown in the table, the urban canyon of Experiment 4 is narrower than that of Experiment 3, in other words, near half of the sky is blocked by buildings. Reference satellite selection consequently becomes more challenging. The selected master satellite shows a higher possibility of being a multipath one that easily contributes error to the SD approached.

In a challenging environment, like that of Experiment 4, the performance for estimating receiver clock error by WLS outperforms that by SD. Due to the strong blockage offered by the surrounding buildings, it is challenging for SD approaches to select a healthy reference satellite for differencing.

For instance, Method 1 received 10 m+ of RMS error for all ranging-based 3DMA GNSS methods. In comparison, Method 4 constrains the RMS error to about 7 m.

Figure B1 shows an epoch of the particle scoring (the heatmap) in Experiment 4. It can be observed that Method 1 has two clusters with a high score, resulting in the solution driving to the wrong location.

Method 2 selects the master satellite on each candidate, and the high score cluster is near the ground truth. However, there is another small cluster away similar to the case in Method 1, which will increase the uncertainty.

In Method 4, only one cluster is found, which is near the ground truth, but the size is relatively larger than that of Method 2. This shows that Method 4 can still provide a more reliable performance here, as all high-scoring candidates are concentrated near the ground truth.

### B.2 Positioning statistic on only using L1-band measurements

The statistic of positioning results using only L1-band measurements is shown in Table B2, which is separated into different receiver clock modeling approaches.

In general, the integrated solution to SKY+LBR+SDM with inter-constellation correction shows the best performance among all methods.

### B.3 Positioning statistic on using L1- and L5-band measurements

The statistic of positioning results using L1- and L5-band measurements is shown in Table B3, which is separated into different receiver clock modeling approaches.

In general, the integrated solution to SKY+LBR+SDM with inter-constellation correction shows the best performance among all methods.

- Received December 9, 2020.
- Revision received August 2, 2021.

- © 2021 Institute of Navigation

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.