Urban positioning: 3D mapping-aided GNSS using dual-frequency pseudorange measurements from smartphones

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.

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.
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 opensourced building models are available on the market for free.
The nature of different 3DMA GNSS algorithms varies.Specifically, shadow matching can provide a better performance in the across-street direction, while the rangingbased 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 inte-gration 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 outof-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 L1band 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  (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 rangingbased 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 than 5 m on average.
To conclude, this paper makes three contributions: 1. Introduces the multi-constellation L5-band measurements into 3DMA GNSS using a smartphone 2. Compares the different approaches of the receiver clock estimation for 3DMA GNSS 3. 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.

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).

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: ,5 if satellite  receives both L1-and L5-bands (1) Similar to the carrier-to-noise ratio for satellite , (∕ 0 )  can be expressed as follows: 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.

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.

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-theart ranging-based 3DMA GNSSs are compared, like Skymask 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 rangingbased 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: The pseudorange difference is then used to obtain the weighted geometric mean square error of  , * and the score for candidate  , * , expressed as follows: where  refers to a diagonal matrix that contains the uncertainty for each satellite: 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: 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 L5band, the weighting tunes have a tuning factor of 10 (Ng et al., 2020b): 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.

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: )2 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  ,1  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 L1band  ,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: 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: 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: 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.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: And the pseudorange difference should be dominated by the NLOS error at this stage: 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 skewnormal 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 skymask 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, likelihoodbased 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.

Position solution and integrated solution
In this paper, we integrated the two ranging-based 3DMA GNSS methods and shadow matching.For the integration 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: 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, * : where  refers to the total number of the positionhypothesized 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.

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).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.

Pseudorange quality
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  (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.57m 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 recep-tions 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.

Evaluation of receiver clock bias estimation methods
The estimation methods compared are: 1. 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.
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).4. 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. 5. 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, highscore 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: 1.The integration of L5-band measurements is available for the improvement of 3DMA GNSS positioning performance.2. The utilization of IC clock corrections in the estimation of the receiver clock bias is beneficial to ranging-based 3DMA GNSS.3. Method 4, as the proposed particle-based WLS approach, is more robust compared to other estimation methods on receiver clock bias.

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 5 m.
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, 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.

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-timekinematic (RTK) GNSS, the integrated solution can be utilized as the initial position and visibility estimation.

A C K N O W L E D G M E N T
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.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.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.

F
Flowchart on the proposed 3DMA GNSS [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org]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

F
I G U R E 2 (a)-(d) Location of experiment 1 to 5, respectively; and (e)-(i) Sky-pointing fisheye camera image of experiment 1 to 5, respectively [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org] of the two ranging-based 3DMA GNSS algorithms, the simulated pseudorange is obtained by the following formula: ρ , = { ρ , if a reflection is found by SKY ρ , if no reflection is found by SKY (15) 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) Δ ,SKY = ρ − ρ,SKY , and the result is shown in Figure 3.Only the received satellites with both F I G U R E 3 L1-and L5-band pseudorange difference between simulated pseudoranges [e.g., (3), on Satellite G01, G03, G08, and G30] and the weighted average (Weighted AVG) of all available satellites at ground truth for Experiment 3.For each data set, L1-and L5-band pseudorange differences are separated into the upper row in the red and the lower row in the blue frame.The Y-axis label in the satellite's ID, the elevation angle inside the bracket, and the labeled satellite visibility, LOS and NLOS, are also presented.It should be noted that the color axis is limited within the range of 20 m [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org]F I G U R E 4 Ray-tracing results for satellite (a) G03, a multipath satellite in which the LOS path is not blocked, and the reflection path is found; and (b) G08, an NLOS diffracted satellite in which the LOS path is blocked, and the diffraction path is found (the red point illustrates the ground truth, the yellow line shows the LOS path, and the blue line represents the reflection/diffraction path) [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org]L1-and L5-band measurements are shown in the figure.

F
Multipath labeled satellite G03's (a) pseudorange error in the entire experiment, and (b) probability density function plot [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org]based on the double-differencing method is shown in Figure 5(a) which was introduced in Section 2.3.3. of our previous work

F
Cumulative distribution function (CDF) of L1-and L5-band pseudorange error by using the data collected in all experiments [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org]

F
I G U R E 7 (a)-(c) are heatmaps of SKY, LBR, and SKY+LBR using Method 1, respectively; (d)-(f) are heatmaps of SKY, LBR, and SKY+LBR using Method 4, respectively; (g) is the skymask with satellite, in which color points represent the ∕ 0 value for the corresponding satellites [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org]F I G U R E 8 CDF plots on 2D positioning error of different receiver clock error estimation methods contribute to (a) Skymask 3DMA (SKY), (b) likelihood-based ranging GNSS (LBR), (c) integration of SKY and LBR (SKY+LBR), and (d) a combination of all ranging-based 3DMA (SKY, LBR, and SKY+LBR) [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org]

F
Positioning results for Experiments 1 through 5 utilizing both L1-and L5-band measurements [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org]F I G U R E 1 0 Heatmap of Experiment 5 in (a) SKY, (b) LBR, (c) SDM, (d) SKY+LBR (RNG), and (e) SKY+LBR+SDM.The green star refers to the truth location, and the purple point represents the positioning solution of the corresponding algorithm.The value inside the bracket denotes the 2D positioning error [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org]F I G U R E 1 1 Heatmap of Experiment 1 in (a) SKY, (b) LBR, (c) SDM, (d) SKY+LBR (RNG), and (e) SKY+LBR+SDM.The green star represents the truth location, and the purple point refers to the positioning solution of the corresponding algorithm.The value inside the bracket denotes the 2D positioning error [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org]we can observe that SKY heatmap in 11(a) has drifted away from the ground truth, causing the positioning error to exceed 20 m.
Guohao Zhang https://orcid.org/0000-0002-4661-9058Li-Ta Hsu https://orcid.org/0000-0002-0352-741XR E F E R E N C E S F I G U R E B 1 LBR heatmaps based on different receiver clock error estimation methods: (a) Method 1, (b) Method 2, and (c) Method 4 [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org] Comparison of the two popular ranging-based 3DMA GNSS methods TA B L E 1 3DMA GNSS positioning RMS error of Experiment 3 based on different receiver clock error estimation methods TA B L E 4 Summary on the 3DMA GNSS results using L1-and L5-band measurements, in meters TA B L E 6 Positioning results using L1-band measurements only, in meters Positioning results using L1-and L5-band measurements, in meters TA B L E B 2