Research ArticleOriginal Article
Open Access

Multi-Epoch 3D-Mapping-Aided Positioning using Bayesian Filtering Techniques

Qiming Zhong and Paul D. Groves
NAVIGATION: Journal of the Institute of Navigation June 2022, 69 (2) navi.515; DOI: https://doi.org/10.33012/navi.515

Summary

The performance of different filtering algorithms combined with 3D mapping-aided (3DMA) techniques is investigated in this paper. Several single- and multi-epoch filtering algorithms were implemented and then tested on static pedestrian navigation data collected in the City of London using a u-blox EVK M8T GNSS receiver and vehicle navigation data collected in Canary Wharf, London, by a trial van with a Racelogic Labsat 3 GNSS front-end. The results show that filtering has a greater impact on mobile positioning than static positioning, while 3DMA GNSS brings more significant improvements to positioning accuracy in denser environments than in more open areas. Thus, multi-epoch 3DMA GNSS filtering should bring the maximum benefit to mobile positioning in dense environments. In vehicle tests at Canary Wharf, 3DMA GNSS filtering reduced the RMS horizontal position error by approximately 68% and 57% compared to the single-epoch 3DMA GNSS and filtered conventional GNSS, respectively.

Keywords

1 INTRODUCTION

Navigation and positioning are indispensable in modern life, and GNSS is one of the most widely used technologies. The demand for and potential of navigation and positioning services in urban canyons are enormous. Location-based services and applications have penetrated many aspects of people’s lives, such as in travel, entertainment, and health (Usman et al., 2018). In recent years, positioning modules have been built into a wide range of consumer products such as smartwatches and smartphones, making it easy for people to access their location. The involvement of location information has made many existing applications smarter and easier to use, while also giving rise to many new applications. Some sports apps on mobile phones, for example, record the trajectory of the runner and allow photos taken along the way to be attached. These services and applications, in turn, promote demand for positioning methods for consumer-level equipment that achieves meter-level horizontal accuracy without additional computing resources and power consumption (European GSA, 2018). However, the performance of GNSS in cities is not satisfactory. According to a survey of Android customers in 2018, the low performance of positioning in cities has become a primary concern for users (van Diggelen, 2021b).

There are several reasons for the poor performance of GNSS in cities. In an environment known as an urban canyon, traffic circulates in the streets and many buildings stand on each side of the street. The complex environment makes the full capabilities of GNSS almost impossible to realize. Tall buildings, large vehicles, and leafy trees all affect satellite signals to varying degrees. Blocking and reflection by buildings are considered to be the main sources of deterioration in positioning accuracy (Groves et al., 2013; Ji et al., 2010; Nur et al., 2013). Blocking reduces the number of line-of-sight (LOS) satellites available. In a typical urban canyon, obstructions on each side of the road can block signals perpendicular to the street direction, leaving only a small portion of the sky overhead unobstructed, resulting in the majority of LOS satellites being distributed along the street. This unhealthy distribution often leads to a sharp increase in errors perpendicular to the street (Wang, 2015).

In addition to blocking, reflection effects interfere with GNSS positioning. In cities, especially emerging ones, large sheets of metallized glass are widely used on the outer surface of buildings. These are powerful radio reflectors that can easily reflect GNSS satellite signals. Users may receive a mixture of LOS satellite signals and corresponding reflected replicas from smooth buildings and the ground at the same time, which is known as multipath interference (Groves, 2013; McGraw et al., 2020; Misra & Enge, 2010; Ward et al., 2017). When the direct and reflected signals are mixed, the resultant signal received by the receiver is distorted, resulting in a ranging error. The carrier phase, Doppler shift, and signal strength measurements may also be impacted. The phenomenon of non-line-of-sight (NLOS) reception is the simultaneous occurrence of reflection and blocking. In other words, the LOS signal is blocked, and the user can only receive the reflected replicas. Since reflection may significantly increase the signal propagation path, NLOS reception would introduce considerable positive errors of tens of meters, or even hundreds of meters in some extreme cases. NLOS signal strength varies drastically. Some signals are weak, while others are as strong as ordinary LOS signals (Groves & Jiang, 2013; Groves et al., 2013; Wang, 2015). Multipath interference and NLOS reception, which occur frequently, are the two key causes of restricted conventional GNSS positioning performance in urban canyons.

The emergence of 3D-mapping data provides additional information for navigation and positioning, and also enriches methods for dealing with issues such as multipath and NLOS reception that often occur in urban environments. Those 3D models are often used to predict, at any given location, which satellite signals are directly visible and which are blocked by obstacles, and even, in some implementations, to estimate path delays. It certainly provides additional useful information for positioning. The methods for implementing 3DMA GNSS are flexible and varied. In the past few years, many different 3D mapping-aided (3DMA) GNSS techniques (Groves et al., 2020; Hsu et al., 2016; Ng et al., 2020; Suzuki, 2016; Suzuki & Kubo, 2013; Wang et al., 2015; Ziedan, 2017, 2019) have been demonstrated to significantly improve the performance of GNSS in cities.

Location-based services typically use single-epoch positioning, while pedestrian and vehicle navigation applications use filtered solutions. Compared to single-epoch solutions, filtered solutions are less affected by noise-like errors that occur during the measurement process. Currently, filtering algorithms commonly used in positioning and navigation include the extended Kalman filter (EKF) and the particle filter (PF), both of which can handle the non-linear distribution of state estimates. The involvement of 3D-mapping data may be able to further improve positioning performance. Based on the characteristics of University College of London’s (UCL) 3DMA GNSS core algorithm (Groves et al., 2020), this paper demonstrates two GNSS filters embedded within the 3DMA algorithms for a multi-epoch case, namely the 3DMA particle filter (PF) and the 3DMA grid filter (GF). The main difference between them is the method used to represent the state estimates. PFs use a fixed number of unevenly distributed particles, while GFs use a variable number of uniformly distributed particles.

Two GNSS data sets are used to test the algorithms mentioned in this paper. The signals in the first data set were collected with a u-blox EVK M8T GNSS receiver in the City of London, which represents a traditional European city with predominantly masonry buildings. The other data set was recorded with a Racelogic Labsat 3 GNSS front-end in the Canary Wharf area of London, which has many glass-covered tall buildings similar to those found in North American and Asian cities.

This paper starts with the problem of GNSS in urban positioning and briefly reviews the various 3DMA GNSS positioning techniques. Then, detailed descriptions of multi-epoch 3DMA GNSS algorithms based on particle filtering and grid filtering are given in Section 3. Section 4 shows the positioning results from, respectively, single-epoch conventional GNSS, single-epoch 3DMA GNSS, multi-epoch conventional filtering, and multi-epoch 3DMA filtering in two test environments, followed by a comparison and analysis of the results of these algorithms. Finally, conclusions are summarized in Section 5.

2 BACKGROUND

2.1 3D-Mapping-Aided GNSS

The additional information that 3D-mapping data can provide for positioning and navigation (including but not limited to terrain height, building location, building orientation, and roof height) enriches and enables new methods for dealing with issues such as multipath and NLOS reception that often occur in urban environments.

Some studies (Amt & Raquet, 2006; Groves & Jiang, 2013) show that terrain height is a helpful supplement to the geometric distribution of satellites to improve the accuracy of positioning. The prerequisite for terrain height aid to be effective is that the accuracy of the terrain height obtained from the digital terrain model (DTM) must be higher than that of the pseudorange measurement.

In addition to terrain height data, 3D building models are often used to assist in positioning. These models can be used to not only predict the visibility of satellites at a given position (Bradbury et al., 2007; Suzuki & Kubo, 2015; Wang, 2015; Wang et al., 2012), but also estimate the distance traveled by the signal (Betaille et al., 2013; Ercek et al., 2006; Gu & Kamijo, 2017; Hsu et al., 2015; Ng et al., 2020; Zhang, Ng, et al., 2020). According to the different principles of solution determination, positioning and navigation algorithms using 3D-mapping data can be divided into shadow matching and 3DMA ranging. The former uses signal strength measurements, while the latter uses pseudoranges similar to the conventional GNSS method.

Shadow matching adopts the idea of pattern matching and determines a solution by comparing received signal strength with satellite visibility prediction at a series of candidate positions. In early research (Wang et al., 2013), the measured satellite visibility was seen as a binary value obtained by a hard threshold on the carrier-to-noise density power ratio (C/N0) value. The degree of matching between the visibility prediction and measurement uses the exclusive-not or (XNOR) logical operation that returns true if its inputs are the same, otherwise false. The position solution is determined by a weighted average of the coordinates of several candidate points with high scores. Some subsequent studies (Isaacs et al., 2014; Wang, 2015; Wang et al., 2015; Zhang et al., 2020) have shown better performance using probability-based satellite visibility and Bayesian theory-based matching determination. Some different research groups (Isaacs et al., 2014; Yozevitch & ben Moshe, 2015) have also demonstrated continuous positioning based on particle filtering.

There are many different approaches to 3DMA ranging. One of the most intuitive approaches is to exclude NLOS signals detected by using 3D models from the calculation (Ng & Hsu, 2021; Obst et al., 2012; Peyraud et al., 2013), which usually requires a fairly accurate initial position to enable subsequent NLOS detection algorithms to confidently predict satellite visibility without much time. These methods are mostly used in continuous positioning, as most positioning applications cannot provide a sufficiently accurate solution within a few seconds of launching.

Many research groups tend to use NLOS measurements instead of simply deleting them. Hypothesis testing is one of the most commonly used methods. At a series of candidate positions generated around a rough position solution, the path delay of the NLOS signal can be estimated by its 3D building model. These candidate positions are then scored based on path prediction and the corresponding actual measurement (Gu & Kamijo, 2017; Hsu etal., 2015, 2016; Zhang et al., 2020). However, the primary limitation of using path delays is that the propagation path calculation requires a large amount of computing resources. A less computationally intensive method (Groves et al., 2020) exploits the difference in symmetry between the pseudorange error distributions of LOS and NLOS signals. It uses a different combination of error distributions at each candidate position based on the visibility predictions from 3D-mapping data, which enables those NLOS pseudoranges to participate in the position calculation without explicitly computing the additional distance traveled by them. Recently, Google has incorporated 3DMA GNSS into their Android system to provide better location services for smartphone users in many cities (van Diggelen, 2021a).

2.2 UCL’s 3DMA GNSS Core Algorithms

UCL’s 3DMA GNSS algorithms consist mainly of shadow matching, likelihood-based ranging, and an integration algorithm, as shown in Figure 1. Both shadow matching and likelihood-based ranging are performed in the way of hypothesis testing on candidate positions. The candidate positions are a set of three-dimensional coordinates. For land positioning and navigation applications, in order to reduce the complexity of the problem, the height dimension is set to the sum of the terrain height at that horizontal position and the height of the user device above the ground.

FIGURE 1
FIGURE 1

Components of the 3DMA GNSS core algorithms

3D-mapping data is used to predict the visibility of each satellite signal (i.e., LOS or NLOS) at each candidate position. This step is relatively computationally intensive and time consuming. Therefore, an intermediate step called building boundaries has been introduced to achieve the goal of being able to operate in real time over a large number of candidate positions. The building boundary refers to the maximum elevation of all buildings within a certain distance at a given azimuth. In other words, it is the minimum elevation of a satellite that allows users to receive its signals directly in that direction. The building boundary is pre-computed and stored for each candidate position. When required, the signal can be classified as LOS or NLOS by simply comparing the satellite elevation with the building boundary at the corresponding azimuth.

The shadow-matching algorithm (Wang et al., 2015) compares satellite visibility predictions with a counterpart determined by the received signal strength to calculate the degree of matching at different candidate positions, thereby giving the optimal solution. The algorithm is comprised of the following steps (Groves et al., 2020; Wang et al., 2015):

  • The predicted visibility of each satellite signal at each candidate position is obtained.

  • For each received signal, the probability that it is LOS is determined from the measurement of the carrier-to-noise density power ratio, C/N0, using an appropriate statistical model.

  • A matching score is obtained by evaluating each satellite at each candidate position based on the match between its predicted visibility and measured C/N0.

  • The final score for each candidate position is a combination of the matching scores for each satellite at that position.

The likelihood-based ranging algorithm (Groves et al., 2020) applies different statistical distributions to pseudorange errors according to satellite visibility predictions, and then evaluates the correspondence between the measured and predicted pseudoranges to give the positioning solution. The algorithm is comprised of the following steps (Groves et al., 2020):

  • The predicted visibility of each satellite signal at each candidate position is obtained.

  • At each candidate position, one of the satellites predicted to be LOS is selected as the reference.

  • At each candidate position, the measurement innovation for each satellite is obtained by subtracting the LOS range and known errors, such as satellite clock errors, atmospheric delays, and inter-constellation offsets, from the measured pseudorange, and then differencing with respect to the reference satellite to remove receiver clock offset.

  • At each candidate position, the cumulative probability of the measurement innovation on a skew-normal distribution is determined for each satellite predicted to be NLOS. These NLOS innovations are then replaced by corresponding LOS innovations with the same cumulative probability.

  • The final score for each candidate position is calculated using the modified measurement innovations and their error covariance matrix.

The intention of shadow matching is to improve the accuracy in the direction perpendicular to the street, whereas likelihood-based ranging is considered to be more accurate in the direction along the street. Therefore, a hypothesis-domain integration algorithm is executed to give a comprehensive single score for each candidate position based on the scoring surfaces from shadow matching and likelihood-based ranging. Finally, the position solution is obtained using the combined scores to weight the candidate positions.

Single-epoch positioning using the 3DMA GNSS core algorithms described above has been demonstrated in the high-density central area of Canary Wharf (Groves et al., 2020).

2.3 Conventional Multi-Epoch GNSS

All road and pedestrian navigation applications use filtered solutions. Filtering algorithms use new measurements to correct navigation solutions predicted from previous information. Specifically, the previous clock drift and drift rate are used to predict the current counterparts, and the previous position and velocity solutions are used to give predictions of current position and velocity. Finally, the current measurements are used to correct the predictions to obtain the final solution.

Filtered solutions mainly have the following three advantages. First, the code tracking noise can be smoothed, which also reduces the negative effects of multipath errors when moving. Second, since more information is available to compare each measurement with, the sensitivity of outlier detection can be improved. Finally, when the number of satellite signals is insufficient to determine the solution, or even when all are shielded, the solution from the previous epoch can also maintain navigation with lower accuracy.

The extended Kalman filter and particle filter are the two most popular multi-epoch GNSS positioning techniques. The extended Kalman filter is a non-linear version of the Kalman filter, which linearizes the state transition and observation models using Taylor’s theorem (Groves, 2013). A particle filter is a sequential Monte Carlo estimation algorithm (Groves, 2013). It uses a set of particles to represent the estimated probability distribution of a set of states, regardless of the form of the distribution. Therefore, we expect filtering to also benefit 3D-mapping-aided GNSS over multiple epochs.

3 3D-MAPPING-AIDED MULTI-EPOCH GNSS

As with conventional multi-epoch GNSS filters, 3DMA GNSS filters need to be able to handle non-linear state estimation. In addition, the 3DMA GNSS algorithm from UCL is built on a hypothesis testing approach that requires the filter to provide a range of position candidates. We therefore propose two different schemes, namely the 3DMA GNSS particle filter and the 3DMA GNSS grid filter. Compared to the extended Kalman filter, which assumes a Gaussian position distribution, the particle and grid filters offer better flexibility. In addition, they are more convenient for generating position candidates.

To represent the estimates, the particle filter and the grid filter use different strategies, which is one of their most significant differences. A set of particles are used by the particle filter in which the total weight of a small cluster of particles determines the likelihood of its location. Conversely, the particles in the grid filter are distributed uniformly in a grid form, which directly describes the probability of their corresponding position being the solution. Figure 2 visualizes the difference in the representation of position estimates between the particle filter and the grid filter, in which the colors of the particles represent their likelihood values and the blank areas represent buildings and locations outside the search area.

FIGURE 2
FIGURE 2

Different representations of position estimates by particle filtering and grid filtering

3.1 3DMA GNSS Particle Filter

As the name suggests, the 3DMA GNSS particle filter is an application of the conventional particle filter described in Groves (2013). The conventional particle-filtering scheme is naturally able to incorporate 3DMA GNSS algorithms as its particles can play the role of position candidates. Figure 3 illustrates the six stages of the multi-epoch 3DMA particle filter implemented in this paper. Three of the components (i.e., system propagation, position measurement update, and resampling) remain consistent with conventional GNSS particle filtering, while the initialization, particle probability update, and velocity filtering are modified to use 3D mapping data.

FIGURE 3
FIGURE 3

Components of a 3DMA GNSS particle filter

3.1.1 Initialization

Initialization is divided into two steps to initialize velocity and position, respectively. For position initialization, a standard single-epoch 3DMA GNSS algorithm set (as detailed in Groves et al. [2020]) is used. The output of the 3DMA GNSS algorithm is a likelihood surface and does not need to be further converted to a position solution with covariance, as the particle filter is capable of representing any shape of the probability distribution. The likelihood surface is then directly sampled to generate a set of particles with equal weight. The number of particles, N, should be determined as appropriate, and the initial weight of each particle is assigned as 1/N.

Velocity initialization takes place after the initialization of the position. The least squares algorithm is used to estimate the velocity and clock drift from the pseudorange rate measurements. It weights the pseudorange rate measurements based on the measured signal strength and the average probability of the corresponding satellite being predicted to be LOS at all candidate positions, weighted according to the single-epoch 3DMA GNSS position solution. The velocity initialization proceeds as follows:

  • At each candidate position, satellite visibility is predicted.

  • The LOS probability for each satellite is estimated based on the satellite visibility predictions for each particle and the corresponding 3DMA likelihoods.

  • The satellite signals are weighted according to their signal strength and LOS probabilities.

  • The velocity and clock offset solution are estimated by a least squares method.

Full details are given in Section 2.1 of Appendix A.

3.1.2 System Propagation

In the system propagation stage, the state estimate of each particle changes while its probability remains unchanged. The process is similar to conventional GNSS particle filtering. The particle state used here only contains the position (in easting and northing form), while the velocity is considered separately in Section 3.1.5. First, the sampling of the system noise is performed independently for each particle. Each particle, xp,k, is then propagated separately through the system propagation model given by:

Embedded Image 1

where p denotes the p-th particle, k denotes the k-th epoch, ϕk−1 is the transition function, and wp,k−1 is the randomly generated system noise vector based on the known probability density function (PDF) of wk−1. A simple system transition model that uses the velocity solution from the previous epoch to calculate the displacement is adopted in this paper. The displacement prediction is given by Equation (A30) and its corresponding error covariance is given by Equation (A31).

3.1.3 3DMA GNSS Scoring

The next step is to apply the 3DMA GNSS core algorithms. To reduce the computational load, particle coordinates rounded to the nearest integer multiple of the 1-meter grid spacing are used for satellite visibility prediction. Next, shadow matching and likelihood-based ranging algorithms are executed on each particle (whose coordinates are not rounded here). The scores from these two algorithms are then combined by hypothesis-domain integration to give the 3DMA likelihood for each particle, Embedded Image. Note that the particle here is equivalent to the candidate position in the 3DMA GNSS algorithm description. Full details are given in Section 1 of Appendix A.

3.1.4 Position Measurement Update

In contrast to system propagation, the measurement update phase changes the probabilities but not the state estimates. This step is processed in the same way as its counterpart in conventional GNSS particle filtering. However, the observation likelihoods, Embedded Image, are different from those in a conventional measurement update. Full details are given in Section 2.4 of Appendix A.

3.1.5 3DMA GNSS Velocity Filter

Similar to the initialization step, the position solution and the velocity solution are obtained separately. The velocity and clock drift are maintained by an extended Kalman filter. It is a two-step process. In the prediction step, the velocity state estimates between epochs remain unchanged as there is no information to propagate them, while their error covariance is increased to model the unknown acceleration between epochs and the unknown change in receiver clock drift. In the update step, the measurement errors are weighted based on the measured signal strength and the average probability of the corresponding satellite being predicted to be LOS, with smaller standard deviation being given to measurement errors with higher LOS probabilities and C/N0 values. In static positioning, the velocity filter can be omitted. Full details are given in Section 2.5 of Appendix A.

3.1.6 Resampling

The final phase of 3DMA particle filtering is resampling. As with conventional GNSS particle filters, many of the particle probabilities may degenerate to zero after several consecutive recursions, leaving only a few particles with relatively large weights. Those low-weighted particles crowd computational resources while reducing the number of particles used to represent the core features of the state estimate distribution, resulting in a decrease in estimate performance and efficiency. Therefore, a resampling step is introduced to remove particles with small weights and duplicate particles with significant weights to mitigate the particle degeneracy problem. Resampling can be performed on every fixed number of epochs, or based on the degree of degeneracy, Neff, given by:

Embedded Image 2

falling below a threshold. This paper adopts the latter option with a threshold of 4N/5, below which resampling is performed. Finally, resampled particles are reallocated with equal weight, 1/N.

3.2 3DMA GNSS Grid Filter

Grid filtering uses a set of particles (also called samples) with different likelihoods, uniformly distributed over a certain search space to represent the posterior distribution of some stochastic process given noisy and/or partial measurements. Focusing on the 3DMA grid filter used in this paper, the search space is a horizontal plane represented by the Easting-Northing position coordinates, and the likelihoods of the particles are determined by the 3DMA GNSS core algorithm. A multi-epoch 3DMA grid filter has six phases, shown in Figure 4. The initialization, 3DMA GNSS scoring, and 3DMA velocity filter phases are equivalent to their 3DMA particle filter counterparts.

FIGURE 4
FIGURE 4

Components of a 3DMA GNSS grid filter

3.2.1 Initialization

The initialization of the 3DMA grid filter is almost identical to the corresponding step in the 3DMA particle filter. A single-epoch 3DMA GNSS algorithm is executed to obtain a series of candidate positions and their likelihoods. The pseudorange rate measurements are then weighted by the satellite LOS probabilities derived from the likelihoods and satellite visibilities at these candidate positions to initialize the overall velocity using the least squares method. Full details are given in Section 2.1 of Appendix A.

3.2.2 System Propagation

Position propagation is the process of predicting the state estimate of the next epoch based on currently known information, which moves position candidates and reassigns their likelihoods. This is conducted in three steps as follows:

  • Search area extension

  • Motion-based likelihood redistribution

  • Confidence-based likelihood redistribution

The candidates with high scores can be located anywhere in the search area. Updated position solutions for the current epoch are more likely to appear in and around these high scoring regions than others. Therefore, the search area needs to be expanded large enough to allow high-scoring candidates that fall far from the center and their surrounding areas to be taken into account in subsequent processing.

In motion-based likelihood redistribution, the movement of the candidate positions is carried out according to the system transition model, reflecting the motion of the receiver. It is less likely that the position translation of the receivers exactly conforms to the grid. The translation between epochs is therefore split into an integer grid space and a remainder in the east and north directions, respectively. For the integer part, the likelihood of each candidate position is inherited directly from its likelihood before the translation. For the remainder part, at each candidate position, the likelihood is collectively determined by the likelihoods within a certain range around it.

Confidence-based likelihood redistribution aims to add a small base likelihood for all candidates in the search area, as all candidates have the potential to be the final solution. The base likelihood, in fact, models the error of the 3DMA GNSS core algorithms, with larger values indicating less confidence in the solution given by the 3DMA GNSS algorithm. Full details are given in Section 2.2 of Appendix A.

3.2.3 Search Area Determination

In the previous system propagation step (Section 3.2.2), the search area was expanded and the number of candidate positions grew to about four times its original size. In this step, the search area is redefined according to the likelihood distribution in the expanded search area, which not only reduces the number of candidates, but also centers the search area around the high-scoring candidates. Full details are given in Section 2.3 of Appendix A.

3.2.4 3DMA GNSS Scoring

The 3DMA GNSS scoring process is the same as the corresponding process of the 3DMA GNSS particle filter. The standard single-epoch likelihood-based 3DMA GNSS ranging, shadow matching, and hypothesis-domain integration algorithms are used over the search area determined in the last step. The output is a set of likelihoods for each candidate position in the search area, Embedded Image, where p denotes the candidate position p, and k denotes the k-th epoch. Full details are given in Section 1 of Appendix A.

3.2.5 Position Measurement Update

The position measurement is updated in a similar way to the equivalent part of the 3DMA GNSS particle filter. The position likelihood distribution is updated to incorporate new measurements from the current epoch, k. At each candidate position, the propagated likelihood, Embedded Image is multiplied by the likelihood, Embedded Image derived by the 3DMA GNSS algorithm based on the measurements of the current epoch k to obtain a composite likelihood. Either of the two likelihood surfaces with a higher confidence level will have a higher and narrower peak, making it more likely to dominate the composite distribution. In addition, a constant weighting parameter is introduced to adjust the dominance of the likelihood from the measurement in the combined distribution, regulating the receptivity of the filter to the new measurements. Finally, the position solution is obtained using the combined scores to weight the candidate positions. Full details are given in Section 2.4 of Appendix A.

3.2.6 3DMA GNSS Velocity Filter

The receiver velocities in the 3DMA grid filter are also maintained by a Kalman filter, which is implemented in the same way as in the 3DMA particle filter. Full details are available in Section 2.5 of Appendix A.

4 EXPERIMENTAL TESTS

The first data set consisted of a number of 2-minute GNSS records from the three constellations of GPS, GLONASS, and Galileo. They were collected at various locations in the City of London in July 2017 using a u-blox EVK M8T GNSS receiver at a recording frequency of 1 Hz. The antenna was maintained at a height of 1.1 meters above the ground. The experimental locations are marked in Figure 5. Note that at each test site, there were two sets of data collected in the morning and afternoon, which gave the satellites enough time in orbit to reach significantly different positions, allowing the morning and afternoon data to be independent of each other. Therefore, the afternoon data set was used for tuning the configurable parameters of the positioning and filtering algorithms, while the morning one was used for testing. The City of London is a typical European city. The roads in such cities are generally narrow, and the walls of the buildings are mainly made of masonry.

FIGURE 5
FIGURE 5

True positions in City of London. Background map ©Google Maps

The other data set was collected in Canary Wharf in July 2019 by a van equipped with a Racelogic Labsat 3 GNSS front-end. The measurement data was made up of intermediate frequency samples and was subsequently processed by Focal Point Positioning Ltd. The data set consisted of 1,602 epochs of GPS and Galileo measurements from conventional code tracking. Figure 6 illustrates the true trajectory of the trial vehicle provided by a Novatel iMAR INS/GNSS system. The vehicle departed from the lower-right corner of the map, moved clockwise on the path in the direction indicated by the orange arrows, and traversed the northern area twice before returning to the vicinity of the starting point from the road on the right of the map. The central area of Canary Wharf is marked by a red rectangle in Figure 6. There are many high-rise buildings with glass and steel surfaces in this area, which are common in some large cities in North America and Asia. The environmental characteristics in non-central areas of Canary Wharf are more open than those in the City of London.

FIGURE 6
FIGURE 6

True trajectory of van trial in Canary Wharf, London. Background map ©Google Maps

In order to simulate navigation from different locations, the Canary Wharf data set was divided equally into eight segments of 200 one-second epochs each, marked with different colored dots in Figure 6. Three of the eight segments occured throughout the non-central areas, while the remaining five segments spanned both central and non-central areas, with two segments starting in the central area.

The statistics of the satellites in the two data sets are given in Table B4.

To evaluate the performance of multi-epoch 3DMA GNSS, the following six positioning algorithms were implemented, tested, and compared in the two test data sets mentioned above:

  • Single-epoch conventional GNSS with outlier detection and terrain-height aid

  • Single-epoch 3DMA GNSS

  • Conventional extended Kalman filter (EKF) with terrain-height aid

  • Conventional particle filter (PF) with terrain-height aid

  • 3DMA GNSS particle filter (PF)

  • 3DMA GNSS grid filter (GF)

The single-epoch conventional GNSS algorithm is described in Groves et al. (2020). The single-epoch 3DMA GNSS algorithm is described in both Groves et al. (2020) and Appendix A, and the parameters used can also be found in Appendix A. The implementation of the conventional extended Kalman and particle filters can be found in Groves (2013). The algorithms for 3DMA GNSS filtering are summarized in Section 3 and detailed in Section 2 of Appendix A. Note that the 3DMA GNSS core algorithm already includes an implicit terrain-height aiding algorithm.

In general, it is better to separate the tuning and testing data. For the data set collected in the City of London, the data collected in the afternoon was used to tune the parameters of the algorithms, while the morning data was used to test and generate the results presented in this section. However, for the Canary Wharf data set, the same data was used to tune and test the algorithms due to the limited duration of the data set.

Figures 7 and 9 show the root-mean-square (RMS) position errors obtained from tests on the City of London and Canary Wharf data sets, respectively. As all methods include a terrain-height aiding technique, only errors in the horizontal radial direction were assessed. Figures 8 and 10 illustrate the maximum positioning error at different confidence levels for the solutions of the City of London and Canary Wharf data sets, respectively. The blue bars represent the maximum position error within the 90% confidence interval (i.e., the maximum error after excluding the 10% of solutions with the largest absolute errors). Similarly, the orange bars show the maximum error within the 50% confidence interval. It is worth noting that a very small number of faulty epochs were excluded from the single-epoch positioning in the vehicle test because the number of observed satellites did not meet the minimum requirements for running the least squares ranging algorithm. Detailed results are presented in Appendix B.

FIGURE 7
FIGURE 7

Horizontal radial position root-mean-square error in City of London (stationary)

FIGURE 8
FIGURE 8

Maximum horizontal radial position error at various confidence levels in City of London (stationary)

FIGURE 9
FIGURE 9

Horizontal radial position root-mean-square error in Canary Wharf (vehicle)

FIGURE 10
FIGURE 10

Maximum horizontal radial position error at various confidence levels in Canary Wharf (vehicle)

The results show that the conventional GNSS filtering algorithms performed better than the conventional single-epoch GNSS algorithm in most cases. In vehicle tests, in particular, the overall root-mean-square (RMS) error in the position solution from even the worst performing multi-epoch GNSS algorithm (i.e., the conventional extended Kalman filter) was approximately 40% lower than that of the single-epoch conventional least squares. However, in the static positioning of the City of London data set, the RMS error of the filtered solution was just slightly reduced compared to the single-epoch GNSS, which is not as significant as in the vehicle tests. It is clear that filtering has a greater impact on the vehicle results. One possible reason for this is that, for static positioning, the NLOS and multipath errors are largely correlated over successive epochs, whereas for mobile positioning, they vary more significantly and can thus be more easily mitigated by the filtering process.

With the use of 3DMA GNSS techniques, the single-epoch positioning algorithm improved significantly in vehicle tests. The overall RMS error of the single-epoch 3DMA GNSS position solution was reduced by about a quarter compared to the conventional one, while its maximum position error at 90% and 50% confidence levels was even lower than the conventional filtered results. It can be found from Table B6 in Appendix B that the single-epoch 3DMA GNSS performed better in segments with more epochs in the central region of Canary Wharf. This indicates that 3DMA GNSS is more advantageous in denser environments than in the more open areas, which is also corroborated by the City of London results. In the City of London, the 3DMA GNSS techniques provided little improvement, with only an obvious reduction in error at the 50% confidence level compared to the conventional.

3DMA GNSS techniques and filtering algorithms benefit from each other. With the filtered solution providing a better initial position for the 3DMA GNSS algorithms in the Canary Wharf data set, the position error of the filtered 3DMA GNSS solution was further reduced by more than 60% compared to that of the single-epoch 3DMA GNSS. Compared to the results of the conventional filtering solutions, the filtered solutions with the 3DMA GNSS techniques showed a reduction of approximately 57% and 67% in, respectively, the overall RMS error and the maximum error at 90% confidence. Particularly in denser environments, such as epochs 1,001–1,200, 3DMA GNSS reduced the error in the filtered solution to approximately one third of the conventional filtered solution. In the City of London data set, where 3DMA GNSS performed unsatisfactorily, the introduction of 3DMA GNSS provided little benefit to the filtering.

To gain insight into the reasons behind the poor performance of the 3DMA filter in static positioning, a case study is presented using location C14_N in static positioning and epochs 401–600 in dynamic positioning. Figures 11 and 12 show the position errors in the east and north directions. In static positioning, building geometry remains almost constant and the NLOS and multipath errors are largely correlated over successive epochs, resulting in position solutions not being improved by filtering. However, in dynamic positioning in which the building geometry changes with vehicle movement, the NLOS and multipath errors vary more significantly, allowing the filtering to mitigate the errors more easily.

FIGURE 11
FIGURE 11

Horizontal position error at Site C14_N, City of London (stationary)

FIGURE 12
FIGURE 12

Horizontal position error of Epochs 401-600 in Canary Wharf (vehicle)

Comparing the results from the filtered 3DMA GNSS algorithms, it can be seen that the overall accuracy of the filtered solutions from grid filtering and particle filtering are similar to each other for both static positioning in the City of London and vehicle navigation in Canary Wharf. This is in line with expectations that grid filtering essentially uses the likelihood of a grid of candidate points to represent the distribution density of particles at these candidate positions, which would not change the performance in terms of positional accuracy. Figures 8 and 10 show that the maximum position error of the grid filter solution at the 90% confidence level was slightly lower than that of the particle filter in both data sets, but not significantly different at the 50% confidence level. Specifically, according to the solution error tables listed in the Appendix B, the difference in RMS error between the filtered solutions from these two filters remained within 1 meter in approximately 50% of the segments in the data sets. The largest differences were found in the C10_E and C2_E data segments for the City of London, and in Epochs 200 to 600 for Canary Wharf.

The difference between 3DMA GNSS grid filtering and particle filtering is more significant in terms of computational load. The computational load of the 3DMA GNSS filters depends mainly on the number of candidates to be computed in each epoch. The count of candidate positions can therefore be used to roughly evaluate the efficiency of these filters. For the 3DMA particle filter, the number of candidate positions is the number of particles maintained. Whereas for the 3DMA grid filter, the number of candidate points may vary, depending on the size and spacing of the grid and the proportion of buildings in the region. The shape and size of the grid can be specibed in advance or based on the uncertainty of the position solution from the previous epoch. For example, a circular grid with a radius of 20 meters and a spacing of 1 meter has approximately 1,256 candidate points. If the area occupied by buildings is removed, then the actual number of candidates to be computed is even less. Therefore, the average of the counts is suitable for comparison. The number of particles in the particle filter is determined by the minimum number (in steps of 500) that completely covers all the candidate positions in a search area spaced at 1 m and of the same radius.

The particle filter can benefit to some extent from an increase in the number of particles, while the grid filter can also benefit from a decrease in the spacing. In order to facilitate comparison of their computational loads, the number of particles in the particle filter was adjusted to achieve an accuracy similar to that of the 1-meter interval grid filter. For the Canary Wharf data set, the particle filter processed 5,000 candidate positions per epoch, while the grid filter processed an average of 3,604 candidate positions per epoch. For the City of London data set, the particle filter computed 1,000 candidates per epoch, while the grid filter computed an average of only 523 candidates per epoch. Therefore, extrapolating from the number of candidate positions in each epoch, 3DMA GNSS grid filtering should be approximately 40% to 50% faster than 3DMA GNSS particle filtering. However, this is not always true in a practical implementation. In the static positioning test (i.e., the City of London data set), the grid filtering was indeed almost twice as efficient as the particle filtering. However, in the vehicle navigation test at Canary Wharf, grid filtering consumed as much time as particle filtering, as it traversed all candidate locations (i.e. Equation [A36]) when performing the motion-based likelihood redistribution in the system propagation step, which took a significant amount of time and resources.

The problem of getting lost can happen with any filter, and the 3DMA GNSS filters are no exception. The error in the position solution given by a lost filter increases with time, while the uncertainty in the solution remains consistently smaller than the position error. For the 3DMA GNSS filters, the true position lying outside the search area can be defined as getting lost. The lost problem usually starts with an incorrect state estimate. In the specific case of the 3DMA GNSS filters shown in this paper, getting lost can originate from a position and/or velocity solution with large errors. The true position may be too far away from the position prediction given in the system propagation step to lie within the 3DMA GNSS search area in the current epoch. Due to the limitations of the 3DMA GNSS core algorithm itself, the true position must lie within the search area. Otherwise, it is impossible for the 3DMA GNSS algorithm to provide the correct position solution. Once the true position wanders outside the search area, it is difficult for the lost filter to get the correct position and velocity solution again, and thus difficult to recover in subsequent epochs. Therefore, improving the positioning accuracy of the 3DMA GNSS core algorithm can reduce the chances of the filter getting lost. However, any algorithm can fail in some situations. Artificially expanding the search area to include the true position retains the opportunity for a lost filter to recover in subsequent epochs. Other methods of detecting and recovering lost 3DMA GNSS filters are worth investigating in future research.

5 CONCLUSION

The results show that filtering has a greater impact on the results of mobile positioning with significant movement compared to static positioning. In vehicle tests, the conventional GNSS filtering algorithms improved positioning accuracy by more than 40% compared to conventional single-epoch GNSS, while in static positioning they offered only a slight improvement. The advantages of 3DMA GNSS are more apparent in denser environments than in more open areas. In the single-epoch positioning of the Canary Wharf data set, 3DMA GNSS improved the overall RMS position error by about a quarter. But in more open areas, such as the City of London and the non-central areas of Canary Wharf, it doesn’t bring much benefit.

3DMA GNSS techniques and filtering algorithms benefit from each other. The former provides the latter with a better position solution in the measurement update step, while the latter, in turn, rewards the former with a better initial position and a smaller search area. In the Canary Wharf data set, the filtered 3DMA GNSS solution showed a further reduction in position error of approximately 68% compared to the single-epoch 3DMA GNSS. The filtered results with the 3DMA GNSS techniques reduced the overall RMS error and the maximum error at 90% confidence level by approximately 57% and 67%, respectively, compared to the conventional filtering solution. It can be inferred that multi-epoch 3DMA GNSS filtering should maximize the benefits for mobile positioning in dense environments.

3DMA GNSS grid filtering and particle filtering show similar performance in position accuracy. In terms of efficiency, 3DMA GNSS grid filtering is able to achieve solutions with fewer particles with errors comparable to those of 3DMA GNSS particle filtering. Theoretically, grid filtering should consume less resources and run faster than particle filtering. However, the grid filtering in the practical implementation may be slowed down by traversing all candidate positions in the system propagation.

There is still much room for improvement in the 3DMA GNSS filtering algorithm. For example, any improvements to the 3DMA GNSS core algorithm could improve the performance of the 3DMA GNSS filtering algorithm to some extent. Potential topics include optimizing satellite visibility predictions, improving the scoring model for likelihood-based ranging and shadow matching, adding outlier detection, etc. The system propagation step in grid filtering could also be optimized to improve operational efficiency.

HOW TO CITE THIS ARTICLE

Zhong, Q., & Groves, P. D. (2022) Multi-epoch 3D-mapping-aided positioning using Bayesian filtering techniques. NAVIGATION, 69(2). https://doi.org/10.33012/navi.515

ACKNOWLEDGMENTS

Qiming Zhong is funded by the China Scholarship Council and a UCL Engineering Faculty Scholarship. The authors also thank Focal Point Positioning for the conventional pseudoranges generated by their software receiver and Imperial College London for use of their trials van and reference system.

APPENDIX A: A DETAILED DESCRIPTION OF ALGORITHMS

This appendix provides more details of the 3DMA positioning algorithms used to generate the results presented in the paper. As the equipment used differs in the two data sets mentioned in this paper, a set of tuning parameters is applied to each data set individually.

1 3DMA CORE ALGORITHMS

The input to the core 3DMA algorithm is a 3D city model of a region and a series of candidate positions within the region, and the output is a scoring surface (also known as likelihood distribution) over these candidates. The score for each candidate position is obtained by the following four algorithms.

1.1 Visibility Prediction

Once a grid of candidate positions has been set up, the visibility of the satellites at these candidates can be predicted. The same set of predictions is used in all subsequent algorithms including shadow matching and likelihood-based ranging. For each candidate position, the probability that the satellite is predicted to be direct line-of-sight (LOS), p(LOS|BB), is determined using the building boundaries pre-computed from the 3D-mapping data. For the City of London data set, p(LOS|BB) was set to 0.85 where the satellite elevation was above the building boundary and 0.15 otherwise. For the Canary Wharf data set, p(LOS|BB) was set to 0.8 and 0.2 for satellites predicted to be LOS and NLOS, respectively. These values account for the possibility that predictions may be wrong due to errors and resolution limitations in the 3D city model and unpredictable factors such as passing vehicles.

1.2 Shadow Matching

Once a set of candidate positions and the satellite visibility predictions at them have been determined, shadow matching is comprised of the following steps:

  1. The probability, p(LOS|C/N0), that each received signal is LOS is determined from the GNSS receiver’s C/N0 measurements using the following statistical model:

    Embedded Image A1

    where:

    Embedded Image A2

    and the coefficients are given in Tables A1 and A2.

    View this table:
    TABLE A1

    Tuning Parameters for Determining LOS Probability from Measured C/N0 in City of London

    View this table:
    TABLE A2

    Tuning Parameters for Determining LOS Probability from Measured C/N0 in Canary Wharf

  2. Each candidate position and satellite is then scored according to the match between the predicted and measured satellite visibility:

    Embedded Image A3

  3. An overall likelihood score, ΛSp, for each position, p, is obtained by multiplying the scores for each satellite.

1.3 Likelihood-Based Ranging

Once a set of candidate positions and the satellite visibility predictions at them have been determined, the likelihood-based ranging algorithm is comprised of the following steps:

  1. For each candidate position, the satellites are ranked based on the elevation angle, measured C/N0, and surrounding buildings, and the one with the highest score is selected as the reference. The score is calculated by:

    Embedded Image A4

    where j denotes the j-th satellite, C/N0 values are rounded to the nearest multiple of five; Embedded Image is the average value of the difference between the satellite elevation; θj, and the corresponding building boundary, BBp, at the candidate point, p = 1, and its immediate neighbors, p = 2,3, …, np. The neighbors whose distance to the candidate point is less than 1.5 times the grid spacing are considered.

  2. At each candidate position, measurement innovations are obtained by subtracting the computed LOS ranges (refer to Groves [2013] for calculations) from the measured pseudoranges and then differencing with respect to the reference satellite to eliminate the receiver clock offset. Note that some known errors, such as atmospheric delays and satellite clock offsets, have been modeled in the calculation of the range to each satellite at the candidate positions.

  3. The error standard deviation of all errors except for the NLOS path delay is computed as a function of C/N0 using:

    Embedded Image A5

    where a and b are empirically determined constants and are given in Table A3.

    View this table:
    TABLE A3

    Parameters Used in the Likelihood-Based Ranging (LBR) Algorithm

  4. At each candidate position, the measurement innovations for satellites predicted to be LOS are modeled by a normal distribution with a mean of μL, while those for satellites predicted to be NLOS are modeled by a skew-normal distribution with the location, ξN, scale, ωN, and shape, αN, as follows:

    Embedded Image A6

    Embedded Image A7

    Embedded Image A8

    where μN and σN are the mean and standard deviation of the NLOS path delay, respectively, and σr is the error standard deviation of the reference satellite, given in Table A3. The cumulative probability, F, of the NLOS measurement innovation, δzpj, is then computed using:

    Embedded Image A9

    where erf(x) is the integral of the normal distribution and T(x, α) is Owen’s T-function.

    The modified measurement innovation, Embedded Image, is then obtained by solving:

    Embedded Image A10

    For measurements predicted to be LOS, the measurement innovation remains unchanged (i.e., Embedded Image).

  5. To prevent excessively large innovations from producing very low likelihood scores, limiting is applied to each innovation:

    Embedded Image A11

    where ⌊•⌋ and ⌈•⌉ are floor and ceiling functions, respectively, and δzmax is given in the Table A3.

  6. A likelihood score for each candidate position, p, is finally computed using:

    Embedded Image A12

    where Embedded Image is the vector of re-mapped measurement innovations and Cδz, p is the measurement error covariance matrix, given by:

    Embedded Image A13

1.4 Hypothesis Domain Integration

Hypothesis domain integration combines the shadow matching and likelihood-based ranging scores to give a single score for each candidate position:

Embedded Image A14

where Wp is the weighting factor, which is a value related to the number of LOS signals in the current epoch.

Embedded Image A15

where nLOS,p and nNLOS,p are the number of satellites predicted to be LOS and NLOS at candidate position p, respectively, and α is an empirically determined constant. For the City of London data set, α takes a value of 4.6, while for the Canary Wharf data set, α is 3.6.

2 3DMA FILTERING ALGORITHMS

This appendix shows several important components of the 3DMA filtering algorithm excluding the 3DMA core algorithm (described separately in Section 1). Some of these components (i.e., initialization, position measurement update, and velocity filtering) are applicable in both particle filtering and grid filtering, while others are only applicable in grid filtering as the equivalent components are available in conventional particle filtering. For example, 3DMA particle filtering can use the conventional system propagation scheme directly. Note that the velocity initialization and velocity filtering can be omitted in static positioning.

2.1 Particle and Grid Filter Initialization

The initialization is divided into position and velocity components, where the position is initialized by a standard singleepoch 3DMA GNSS algorithm set and the velocity is initialized by a weighted least squares method. For static positioning applications, the velocity is considered known and does not need to be initialized.

The single-epoch 3DMA GNSS for position initialization consists of three steps: the weighted least squares ranging is used to determine a rough position solution from which a grid of candidate positions are generated, and finally the core 3DMA GNSS algorithm scores these candidates to derive a position solution. The detailed algorithm is described in Groves et al. (2020), and the 3DMA GNSS core algorithm and parameters can also be found in Section 1 of the appendix.

The velocity initialization is comprised of the following steps:

  1. For each candidate position, predict the satellite visibility using building boundaries. Note that the predictions are already made in the 3DMA GNSS positioning algorithm and can, therefore, be used directly here without recalculation.

  2. The overall LOS probability of a satellite over the search area is calculated by:

    Embedded Image A16

    where j denotes the satellite, Embedded Image is the likelihood of candidate position p, and Embedded Image is a Boolean value representing the predicted visibility of satellite j at candidate position p.

  3. The variance of the measurement error for each satellite is determined using:

    Embedded Image A17

    where the LOS probability-based part is given by:

    Embedded Image A18

    and the C/N0-based part is given by:

    Embedded Image A19

    where the coefficients σr0 = 0.2, β = 10, ar = 500m2/s2, and br = 0.015m2/s2 are determined empirically.

  4. The measurement innovation, Embedded Image, is:

    Embedded Image A20

    where Embedded Image is the measured pseudorange rate from satellite j, and Embedded Image is the predicted range rate for satellite j after eliminating known errors such as the effect of Earth rotation, calculated using the method in Groves (2013).

    The measurement matrix, Embedded Image, is computed using:

    Embedded Image A21

    where Embedded Image is the LOS unit vector of satellite j, calculated using the method in Groves (2013).

    The measurement error covariance matrix, Cr, is given by:

    Embedded Image A22

    The velocity state vector consists of the velocity, Embedded Image, and the clock drift, Embedded Image, and is obtained by:

    Embedded Image A23

    and the state estimation error covariance is initialized using:

    Embedded Image A24

  5. The east and north velocities are then given by:

    Embedded Image A25

    and their error covariance is:

    Embedded Image A26

    where Embedded Image is the first three columns and rows of Embedded Image.

    The north-east-down (NED) to Earth-centered, Earth-fixed (ECEF) coordinate transformation matrix, Embedded Image, in Equations (A25) and (A26) is given by:

    Embedded Image A27

    where the latitude, La, and longitude, λa, are obtained from the position solution, while Embedded Image is the ECEF to NED coordinate transformation matrix.

2.2 Grid Filter System Propagation

System propagation predicts the state estimate of the next epoch based on currently known information, which moves position candidates and redistributes their likelihoods. The process is divided into the following three steps:

  1. The search area is expanded to ensure that candidates with high scores from the previous epoch (or initialization) fall in the center of the search area in subsequent calculations. In the case of a circular search area, for example, the radius rp is increased to rq = 2rp with the same center.

    Thus, likelihoods for the candidate positions within the expanded area, indexed by q, are initialized using:

    Embedded Image A28

    where Embedded Image is the likelihood at the preceding epoch, k–1, following the position measurement update, of the point pk− that has the same coordinates as point qk, and Pk−1 is the set of candidate positions used for the position measurement update at epoch k – 1. The likelihood of the newly added candidates is temporarily assigned zero, pending further processing in subsequent steps. The coordinate notation of the candidate position is also updated after the search area extension to avoid confusion.

    Embedded Image A29

  2. The position estimate between epochs can be propagated by a velocity solution and its associated error covariance. The predicted displacement between epochs k – 1 and k is:

    Embedded Image A30

    where Embedded Image and Embedded Image are, respectively, the estimated east and north velocities at epoch k – 1, and τs is the time interval between epochs.

    The error covariance of the predicted displacement is:

    Embedded Image A31

    where Embedded Image is the error covariance of the east and north velocities at the preceding epoch k – 1, SaE,k and SaN,k are the power spectral densities (PSDs) of the east and north accelerations, respectively, and SaEN,k is the cross-spectral density of the east and north accelerations. The same values as in the velocity filter, given by (A54), are used.

    To facilitate alignment with the grid in the search area, the displacement is partitioned into integer grid space Embedded Image and remainder Embedded Image components:

    Embedded Image A32

    where Δpg is the grid spacing, and the integer part is given by:

    Embedded Image A33

    The position likelihood is redistributed to account for the integer displacement simply by changing the position associated with each likelihood. Thus:

    Embedded Image A34

    Embedded Image A35

    For the remainder part, the likelihoods are then redistributed using:

    Embedded Image A36

    where the weighting for each epoch is computed by:

    Embedded Image A37

    where:

    Embedded Image A38

    Embedded Image A39

    Embedded Image A40

    In static positioning, since the velocity of the receiver is known to be zero, the motion-based likelihood redistribution does not need to be performed. Thus, the coordinates and likelihoods of the candidate positions are inherited directly from the expanded search area:

    Embedded Image A41

    Embedded Image A42

  3. The confidence-based likelihood redistribution adds a minimum likelihood to all candidates in the search area. Let CP be the level of confidence in the position solution from the preceding epoch (or initialization), where 0 <CP < 1. The likelihoods are redistributed using:

    Embedded Image A43

    where Ap is the area of the search area calculated by:

    Embedded Image A44

    for a circular area with radius, rp, or by:

    Embedded Image A45

    for an elliptical area with semi-major axis, αp, and semiminor axis, bp. CP takes a value of 0.99 in the City of London data set and 0.75 in the Canary Wharf data set.

2.3 Grid Filter Search Area Determination

The center of the search area for the 3DMA GNSS position measurement corresponds to the weighted average position solution obtained from the propagated likelihood grid, giving by:

Embedded Image A46

The likelihood of the candidate positions remains unchanged within the newly defined search area. Thus, the notation is updated by:

Embedded Image A47

The size of the search area can be fixed or vary with the error covariance. In the City of London and Canary Wharf tests, the search area was fixed to a circular area of a 20-m and 40-m radius, respectively.

2.4 Particle and Grid Filter Position Measurement Update

At each candidate position, the likelihood is updated with the propagated likelihood, Embedded Image, and the 3DMA scoring, Embedded Image, by:

Embedded Image A48

where Wm is the empirically determined measurement weighting factor. A value greater than one gives more weighting to new measurement data while a value less than one gives less. Wm takes a value of one for both the City of London and Canary Wharf data sets. The overall likelihood is then normalized by:

Embedded Image A49

Finally, the updated position solution is obtained using:

Embedded Image A50

where Ep,k and Np,k are the easting and northing coordinates of candidate position p.

2.5 Velocity Filter (Used with Particle and Grid Filters)

The velocity filter includes prediction and update steps. In the propagation step, with no additional information, the velocity state estimates remain unchanged, while the error covariance increases. Thus:

Embedded Image A51

Embedded Image A52

where Embedded Image is the receiver clock frequency drift PSD, τs is the time interval between epochs, and Embedded Image is the acceleration PSD matrix that is given by:

Embedded Image A53

where Embedded Image and Embedded Image are transformation matrices between NED and ECEF coordinates as shown in Equation (A27), SαE,k is the east acceleration PSD, SαN,k is the north acceleration PSD, SαEN,k is the cross-spectral density of the east and north accelerations, and SαD,k is the vertical acceleration PSD.

For a typical vehicle, the acceleration is similar in the easting and northing directions, while the vertical acceleration is relatively small. The horizontal acceleration tends to be greater at lower speeds. Therefore, the acceleration PSD model is proposed as follows:

Embedded Image A54

where the constant, Ka = 1.5m2/s3, is determined empirically.

In the measurement update step, the weighting model applied to the measurement errors is similar to the one used in the velocity initialization shown in Equations (A16), (A17), (A18), and (A19).

The measurement noise covariance, Rk, is then obtained using:

Embedded Image A55

where σr,j is given by Equation (A17). The measurement matrix, Embedded Image, is given by:

Embedded Image A56

where Embedded Image is the line-of-sight unit vector of satellite j.

The Kalman gain is then computed using:

Embedded Image A57

The measurement innovation, Embedded Image, is computed using:

Embedded Image A58

where Embedded Image is the measured pseudorange rate from satellite j, Embedded Image is the modeled range rate from satellite j, and Embedded Image is the predicted receiver clock drift.

The state estimate consisting of the velocity, Embedded Image, and clock drift, Embedded Image, is then updated using:

Embedded Image A59

The state estimation error covariance is updated using:

Embedded Image A60

The east and north velocity and their error covariance can then be obtained using Equations (A25–A27).

APPENDIX B: DETAILED EXPERIMENTAL RESULTS

View this table:
TABLE B1

Statistics of the Satellites

View this table:
TABLE B2

Horizontal Radial Position RMS Errors (in Meters) for Different Algorithms in City of London

View this table:
TABLE B3

Horizontal Radial Position RMS Errors (in Meters) for Different Algorithms in Canary Wharf

View this table:
TABLE B4

Maximum Horizontal Radial Position Error (in Meters) for Different Algorithms in City of London, 90% of Confidence

View this table:
TABLE B5

Maximum Horizontal Radial Position Error (in Meters) for Different Algorithms in Canary Wharf, 90% of Confidence

View this table:
TABLE B6

Maximum Horizontal Radial Position Error (in Meters) for Different Algorithms in City of London, 50% of Confidence

View this table:
TABLE B7

Maximum Horizontal Radial Position Error (in Meters) for Different Algorithms in Canary Wharf, 50% of Confidence

Abbreviations

3DMA
3D mapping-aided
EKF
extended Kalman filter
PF
particle filter
GF
grid filter
LOS
line-of-sight
NLOS
non-line-of-sight
SM
shadow matching
LBR
likelihood-based ranging

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.

REFERENCES

  1. ↵
    1. Amt, J. H., &
    2. Raquet, J. F.
    (2006). Positioning for range-based land navigation systems using surface topography. Proc. of the 19th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS 2006), Fort Worth, TX, 1494–1505. https://www.ion.org/publications/abstract.cfm?articleID=6773
  2. ↵
    1. Betaille, D.,
    2. Peyret, F.,
    3. Ortiz, M.,
    4. Miquel, S., &
    5. Fontenay, L.
    (2013). A new modeling based on urban trenches to improve GNSS positioning quality of service in cities. IEEE Intelligent Transportation Systems Magazine, 5(3), 59–70. https://doi.org/10.1109/MITS.2013.2263460
  3. ↵
    1. Bradbury, J.,
    2. Ziebart, M.,
    3. Cross, P. A.,
    4. Boulton, P., &
    5. Read, A.
    (2007). Code multipath modelling in the urban environment using large virtual reality city models: Determining the local environment. The Journal of Navigation, 60(1), 95–105. https://doi.org/10.1017/S0373463307004079
  4. ↵
    1. Ercek, R.,
    2. De Doncker, P., &
    3. Grenez, F.
    (2006). NLOS-multipath effects on pseudo-range estimation in urban canyons for GNSS applications. 2006 First European Conference on Antennas and Propagation, Nice, France. https://doi.org/10.1109/eucap.2006.4584889
  5. ↵
    European GSA. (2018). GNSS User Technology Report Issue 2. European Global Navigation Satellite Systems Agency. https://www.gsa.europa.eu/system/files/reports/gnss_user_tech_report_2018.pdf
  6. ↵
    1. Groves, P. D.
    (2013). Principles of GNSS, inertial, and multisensor integrated navigation systems. Artech House.
  7. ↵
    1. Groves, P. D., &
    2. Jiang, Z.
    (2013). Height aiding, C/N0 weighting and consistency checking for GNSS NLOS and multipath mitigation in urban areas. The Journal of Navigation, 66(5), 653–669. https://doi.org/10.1017/S0373463313000350
  8. ↵
    1. Groves, P. D.,
    2. Jiang, Z.,
    3. Rudi, M., &
    4. Strode, P.
    (2013). A portfolio approach to NLOS and multipath mitigation in dense urban areas. Proc. of the 26th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2013), Nashville, TN, 3231–3247. https://www.ion.org/publications/abstract.cfm?articleID=11264
  9. ↵
    1. Groves, P. D.,
    2. Zhong, Q.,
    3. Faragher, R., &
    4. Esteves, P.
    (2020). Combining inertially-aided extended coherent integration (supercorrelation) with 3D-mapping-aided GNSS. Proc. of the 33rd International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2020), 2327–2346. https://doi.org/10.33012/2020.17767
  10. ↵
    1. Gu, Y., &
    2. Kamijo, S.
    (2017). GNSS positioning in deep urban city with 3D map and double reflection. 2017 European Navigation Conference (ENC), Lausanne, Switzerland. https://doi.org/10.1109/EURONAV.2017.7954196
  11. ↵
    1. Hsu, L. -T.,
    2. Gu, Y., &
    3. Kamijo, S.
    (2015). NLOS correction/exclusion for GNSS measurement using RAIM and city building models. Sensors, 15(7), 17329–17349. https://doi.org/10.3390/s150717329
  12. ↵
    1. Hsu, L. -T.,
    2. Gu, Y., &
    3. Kamijo, S.
    (2016). 3D building model-based pedestrian positioning method using GPS/GLONASS/QZSS and its reliability calculation. GPS Solutions, 20(3), 413–428. https://doi.org/10.1007/s10291-015-0451-7
  13. ↵
    1. Isaacs, J. T.,
    2. Irish, A. T.,
    3. Quitin, F.,
    4. Madhow, U., &
    5. Hespanha, J. P.
    (2014). Bayesian localization and mapping using GNSS SNR measurements. 2014 IEEE/ION Position, Location, and Navigation Symposium (PLANS 2014), Monterey, CA. https://doi.org/10.1109/PLANS.2014.6851402
  14. ↵
    1. Ji, S.,
    2. Chen, W.,
    3. Ding, X.,
    4. Chen, Y.,
    5. Zhao, C., &
    6. Hu, C.
    (2010). Potential benefits of GPS/GLONASS/GALILEO integration in an urban canyon – Hong Kong. The Journal of Navigation, 63(4), 681–693. https://doi.org/10.1017/S0373463310000081
  15. ↵
    1. McGraw, G. A.,
    2. Groves, P. D., &
    3. Ashman, B. W.
    (2020). Robust positioning in the presence of multipath and NLOS GNSS signals. In Y. T. Jade Morton, F. van Diggelen, J. J. Spilker Jr.., B. W. Parkinson, S. Lo, & G. Gao (Eds.), Position, navigation, and timing technologies in the 21st century: Integrated satellite navigation, sensor systems, and civil applications, volume 1 (pp. 551–589). John Wiley & Sons, Ltd. https://doi.org/10.1002/9781119458449.ch22
  16. ↵
    1. Misra, P., &
    2. Enge, P.
    (2010). Global positioning system: Signals, measurements, and performance. Ganga-Jamuna Press.
  17. ↵
    1. Ng, H. -F., &
    2. Hsu, L. -T.
    (2021). 3D mapping database-aided GNSS RTK and its assessments in urban canyons. IEEE Transactions on Aerospace and Electronic Systems, 57(5), 3150–3166. https://doi.org/10.1109/TAES.2021.3069271
  18. ↵
    1. Ng, H. -F.,
    2. Zhang, G., &
    3. Hsu, L. -T.
    (2020). A computation effective range-based 3D mapping aided GNSS with NLOS correction method. The Journal of Navigation, 73(6), 1202–1222. https://doi.org/10.1017/S037346332000003X
  19. ↵
    1. Nur, K.,
    2. Feng, S.,
    3. Ling, C., &
    4. Ochieng, W.
    (2013). Integration of GPS with a WiFi high accuracy ranging functionality. Geo-Spatial Information Science, 16(3), 155–168. https://doi.org/10.1080/10095020.2013.817106
  20. ↵
    1. Obst, M.,
    2. Bauer, S., &
    3. Wanielik, G.
    (2012). Urban multipath detection and mitigation with dynamic 3D maps for reliable land vehicle localization. Proc. of the 2012 IEEE/ION Position, Location, and Navigation Symposium, Myrtle Beach, SC, 685–691. https://doi.org/10.1109/PLANS.2012.6236944
  21. ↵
    1. Peyraud, S.,
    2. Bétaille, D.,
    3. Renault, S.,
    4. Ortiz, M.,
    5. Mougel, F.,
    6. Meizel, D., &
    7. Peyret, F.
    (2013). About non-line-of-sight satellite detection and exclusion in a 3D map-aided localization algorithm. Sensors, 13(1), 829–847. https://doi.org/10.3390/s130100829
  22. ↵
    1. Suzuki, T.
    (2016). Integration of GNSS positioning and 3D map using particle filter. Proc. of the 29th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2016), Portland, OR, 1296–1304. https://doi.org/10.33012/2016.14857
  23. ↵
    1. Suzuki, T., &
    2. Kubo, N.
    (2013). Correcting GNSS multipath errors using a 3D surface model and particle filter. Proc. of the 26th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2013), Nashville, TN, 1583–1595. https://www.ion.org/publications/abstract.cfm?articleID=11171
  24. ↵
    1. Suzuki, T., &
    2. Kubo, N.
    (2015). Simulation of GNSS satellite availability in urban environments using Google Earth. Proc. of the ION 2015 Pacific PNT Meeting, Honolulu, HI, 1069–1079. https://www.ion.org/publications/abstract.cfm?articleID=12782
  25. ↵
    1. Usman, M.,
    2. Asghar, M. R.,
    3. Ansari, I. S.,
    4. Granelli, F., &
    5. Qaraqe, K. A.
    (2018). Technologies and solutions for location-based services in smart cities: Past, present, and future. IEEE Access, 6, 22240–22248. https://doi.org/10.1109/ACCESS.2018.2826041
  26. ↵
    1. van Diggelen, F.
    (2021a). End game for urban GNSS: Google’s use of 3D building models. Inside GNSS. https://insidegnss.com/end-game-for-urban-gnss-googles-use-of-3d-building-models
  27. ↵
    1. van Diggelen, F.
    (2021b). Google’s use of 3D building models to solve urban GNSS. ION ITM/PTTI 2021.
  28. ↵
    1. Wang, L.
    (2015). Investigation of shadow matching for GNSS positioning in urban canyons [Doctoral dissertation, University College London]. UCL Discovery. https://discovery.ucl.ac.uk/id/eprint/1464060
  29. ↵
    1. Wang, L.,
    2. Groves, P. D., &
    3. Ziebart, M. K.
    (2012). Multi-constellation GNSS performance evaluation for urban canyons using large virtual reality city models. The Journal of Navigation, 65(3), 459–476. https://doi.org/10.1017/S0373463312000082
  30. ↵
    1. Wang, L.,
    2. Groves, P. D., &
    3. Ziebart, M. K.
    (2013). Urban positioning on a smartphone: Real-time shadow matching using GNSS and 3D city models. Proc. of the 26th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2013), Nashville, TN, 1606–1619. https://www.ion.org/publications/abstract.cfm?articleID=11339
  31. ↵
    1. Wang, L.,
    2. Groves, P. D., &
    3. Ziebart, M. K.
    (2015). Smartphone shadow matching for better cross-street GNSS positioning in urban environments. The Journal of Navigation, 68(3), 411–433. https://doi.org/10.1017/S0373463314000836
  32. ↵
    1. Ward, P. W.,
    2. Betz, J. W., &
    3. Hegarty, C. J.
    (2017). GNSS Disruptions. In E. D. Kaplan & C. Hegarty (Eds.), Understanding gps/gnss - principles and applications (3rd ed., pp. 549–617). Artech House.
  33. ↵
    1. Yozevitch, R., &
    2. ben Moshe, B.
    (2015). A robust shadow matching algorithm for GNSS positioning. NAVIGATION, 62(2), 95–109. https://doi.org/10.1002/navi.85
    1. Zhang, G.,
    2. Ng, H. -F.,
    3. Wen, W., &
    4. Hsu, L. -T.
    (2020). 3D mapping database aided GNSS based collaborative positioning using factor graph optimization. IEEE Transactions on Intelligent Transportation Systems, 22(10), 6157–6187. https://doi.org/10.1109/tits.2020.2988531
    1. Zhang, G.,
    2. Wen, W.,
    3. Xu, B., &
    4. Hsu, L. -T.
    (2020). Extending shadow matching to tightly-coupled GNSS/INS integration system. IEEE Transactions on Vehicular Technology, 69(5), 4979–4991. https://doi.org/10.1109/TVT.2020.2981093
  36. ↵
    1. Ziedan, N. I.
    (2017). Urban positioning accuracy enhancement utilizing 3D buildings model and accelerated ray tracing algorithm. Proc. of the 30th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2017), Portland, OR, 3253–3268. https://doi.org/10.33012/2017.15366
  37. ↵
    1. Ziedan, N. I.
    (2019). Enhancing GNSS mobile positioning in urban environments through utilization of multipath prediction and consistency analysis. Proc. of the 32nd International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2019), Miami, FL, 3469–3483. https://doi.org/10.33012/2019.16929
Back to top

In this issue

Print
Download PDF
Article Alerts
Email Article
Citation Tools
Multi-Epoch 3D-Mapping-Aided Positioning using Bayesian Filtering Techniques
Qiming Zhong, Paul D. Groves
NAVIGATION: Journal of the Institute of Navigation Jun 2022, 69 (2) navi.515; DOI: 10.33012/navi.515
Twitter logo Facebook logo Mendeley logo
Bookmark this article

Jump to section

Keywords