GNSS Urban Positioning with Vision-Aided NLOS Identification

Abstract

The global navigation satellite system (GNSS) has played an important role in a broad range of consumer and industrial applications. In particular, cities have become GNSS major application scenarios; however, GNSS signals suffer from blocking, reflection and attenuation in harsh urban environments, resulting in diverse received signals, e.g., non-line-of-sight (NLOS) and multipath signals. NLOS signals often cause severe deterioration in positioning, navigation, and timing (PNT) solutions, which should be identified and excluded. In this paper, we propose a vision-aided NLOS identification method to augment GNSS urban positioning. A skyward omnidirectional camera is installed on a GNSS antenna to collect omnidirectional images of the sky region. After being rectified, these images are processed for sky region segmentation, which is improved by leveraging gradient information and energy function optimization. Image morphology processing is further employed to smooth slender boundaries. After sky region segmentation, the satellites are projected onto the omnidirectional image, from which NLOS satellites are identified. Finally, the identified NLOS satellites are excluded from GNSS PNT estimation, promoting accuracy and stability. Practical test results show that the proposed sky region segmentation module achieves over 96% accuracy, and that completely accurate NLOS identification is achieved for the experimental images. We validate the performance of our method on public datasets. Compared with the raw measurements without screening, the vision-aided NLOS identification method enables improvements of 60.3%, 12.4% and 63.3% in the E, N, and U directions, respectively, as well as an improvement of 58.5% in 3D accuracy.

1. Introduction

The global navigation satellite system (GNSS) can provide convenient and efficient positioning services, and has become an essential part of a variety of products [1]. With the development of the GNSS, the requirements for positioning services in urban environments have increased, presenting steep challenge [2]. GNSS signals are essentially electromagnetic wave signals, which are reflected by glass or metal surfaces and blocked by traffic road bridges and tall buildings. Thus, pseudorange errors are introduced due to both the multipath effect and non-line-of-sight (NLOS) reception. Multipath signals usually include line-of-sight (LOS) and reflection/diffraction signals, while NLOS reception includes only reflection/diffraction signals. The pseudorange observation error caused by the multipath effect is generally at the meter level, and the carrier phase observation error is at the centimeter level. NLOS reception can lead to positioning errors of more than 100 m in deep urban canyons [3]. Clearly, the requirements for high-precision positioning information often cannot be met in complex urban environments, such as unmanned driving and pedestrian navigation. To improve the accuracy of GNSS positioning, NLOS signals should be effectively excluded. Given that satellite navigation still occupies a dominant status for vehicle and pedestrian navigation, several promising solutions for NLOS elimination have been proposed by scholars and engineers.

Detecting NLOS signals based on raw GNSS measurements is straightforward. Ref. [4] proposed a method to mitigate the multipath effect and select the proper signals based on signal-to-noise ratio (SNR) measurements. Ref. [5] presented an improved algorithm to identify and discard measurements that are corrupted by multipath and NLOS errors. In addition, determining whether a GNSS measurement is an NLOS or LOS signal has been regarded as a typical classification problem. Hence, machine learning methods are very popular for NLOS signal detection. Ref. [6] developed a LOS/NLOS classifier algorithm based on support vector machine (SVM), which makes full use of GNSS satellites’ raw measurements, such as SNR, normalized pseudorange residual, and elevation angle, to discuss satellite visibility. Shadow matching was integrated with the proposed SVM LOS/NLOS classifier to eliminate NLOS signals and improve the robustness and accuracy of single point positioning (SPP). However, in complex urban areas, methods based on raw GNSS measurements do not significantly improve the positioning accuracy. Moreover, some methods using machine learning require a large number of training samples, and the generalization ability of a model is also a problem worth considering. In addition to NLOS signal detection methods from raw GNSS measurements, some scholars have proposed light detection and ranging (LiDAR)-based methods. LiDAR-based three-dimensional reconstruction and mapping is an environmental perception strategy. Accordingly, if we can collect prior information about the space in which a vehicle drives, the NLOS signals can be effectively rejected by combining the satellite visibility analysis [7], or raw measurements can be acquired to correct the NLOS delay in a pseudorange model [8]. However, accurate 3D building models are expensive to establish and maintain, and the use of such models places high demands on the performance of computing devices. Furthermore, an accurate 3D model must be created in advance. Some moving objects with a tall height are not among the buildings modeled in the 3D city map, which can also lead to NLOS reception [9]. Therefore, the accuracy and timeliness limit the popularity of this technology to a certain extent [10].

Compared to using LiDAR-based methods to detect NLOS signals, a vision-aided NLOS elimination method has a lower sensor cost. A vision-aided method uses a skyward camera to detect the visible sky and other objects present in the field of view. Then, the satellites are projected onto the processed image to identify whether signals can be received directly (in a sky region) or not (in a non-sky region). An important problem in vision-aided methods is the segmentation of sky and non-sky regions of the images captured by the camera. Sets of unsupervised and supervised clustering algorithms have been compared to determine the best sky region classifier in terms of good classification rate and processing time [11]. The method of [12] first divides the image into small regions. Then, the feature points in the image are tracked, and the distance traveled in the image is computed to identify those small clusters that correspond to obstacles and extract the sky regions. In addition, methods based on the hue–saturation–value (HSV) color model [13] and deep learning [14] are commonly used to identify sky regions. However, these methods cannot address the reflective phenomena on the upper floors and surfaces of buildings in urban environments, as shown in Figure 1, leading to the misidentification of sky regions. In terms of positioning results, both [14,15] conducted only experiments in static scenes, which cannot well simulate the actual effect when driving a vehicle in an urban environment.

In this paper, we propose a vision-aided method to exclude the NLOS signals caused by tall buildings in urban scenarios. Specifically, images taken by a skyward omnidirectional camera are used for sky region segmentation. Sky regions are detected by leveraging the gradient information and energy function optimization. For slender sky boundaries caused by strong light, we utilize an image morphology processing method to restore the misdetected sky regions. Then, we project satellites on the omnidirectional image after sky region segmentation, with the intent of determining which satellites are located in visible sky regions and which are hidden behind buildings and sending NLOS signals. Finally, the GNSS positioning result is calculated based on the remaining visible satellites. The method proposed in this paper significantly improves positioning accuracy in urban environments with the aid of a low-cost sensor, and the calculation is simple. The proposed method does not require a large number of samples for training and does not need to consider the generalization ability of the model.

The remainder of this paper is structured as follows. Section 2.1 provides a detailed description of the sky region segmentation algorithm. Section 2.2 discusses how NLOS signals are rejected. In Section 2.3, a code-based kinematic positioning algorithm is implemented with an extended Kalman filter for the estimation of positioning. In Section 3, we evaluate the effectiveness of the proposed method by means of experiments. Finally, Section 4 summarizes the work by presenting several conclusions and provides an outlook on future research.

2. Methodology

The vision-aided GNSS urban positioning system proposed in this paper is composed of three sensors: a GNSS receiver, an inertial measurement unit (IMU), and an omnidirectional camera. The images captured by the omnidirectional camera are first rectified, and then sky region segmentation is performed using an improved method based on gradient information and energy function optimization, which can efficiently extract the sky regions in urban canyons, under flyovers and in other challenging environments. The proposed approach using gradient information and energy function optimization is different from the above-mentioned sky region segmentation methods. The proposed method is able to correctly identify the reflective phenomena generated by the surfaces and windows of tall buildings. Therefore, the method in [16] is more suitable for the image characteristics of omnidirectional camera images after rectification and has better performance in solving the problem of reflection from building surfaces in urban environments. In our work, we developed improvements based on [16] for the problem of sky region segmentation when buildings, viaducts and flyovers block the field of view of the camera. After sky region segmentation has been completed, the omnidirectional image is recovered in accordance with the remapping parameters. The IMU is used to obtain the angle between the system and the true north direction to align the satellite skyplot with the omnidirectional image. The imaging principle of the omnidirectional camera is used to project the satellites onto the omnidirectional image after sky region segmentation for visibility analysis. Finally, the global position of the system is obtained using code-based kinematic positioning after the elimination of the identified NLOS satellites. A flow chart of the system is shown in Figure 2.

2.1. Sky Region Segmentation

The sky region segmentation module proposed in this paper first rectifies the input omnidirectional image. We perform sky region segmentation on the image after rectification using a method based on improved gradient information and energy function optimization [16]. The processed images are converted back to omnidirectional images using remapping parameters after sky region segmentation is completed. Then, image morphology processing is applied to the resulting slender boundaries. The whole process is shown in Figure 3.

2.1.1. Omnidirectional Image Rectification

For vision-aided NLOS identification, the selection of a suitable camera must be considered. Excellent details in images are not what we are seeking; instead, our needs are mainly based on the shooting angle. The simplest choice is a normal-scale camera but with a limited shooting angle, while a wide-range camera can slightly improve the shooting angle but at the cost of inevitable distortion. Therefore, an omnidirectional camera can be considered the best selection because of its affordability, zero blind spot features and storage optimization.

The first step after the camera captures an image is to rectify the omnidirectional image. Omnidirectional images have very large distortions, so they are not compatible with the sky region segmentation algorithm proposed in this paper. To rectify omnidirectional images, the imaging principle of the omnidirectional camera should be specified. An imaging system with a single effective viewpoint is called a central projection system. A conventional perspective camera is an example of a central projection system. Most commercial cameras can be described as pinhole cameras, modeled by the well-known perspective projection. However, to widen the camera field of view, an omnidirectional camera uses a combination of mirrors (such as parabolic, hyperbolic, or elliptical mirrors) and a standard perspective camera that provides a 360-degree field of view in the horizontal plane. A unified projection model for a central catadioptric camera was proposed in [17,18]. All omnidirectional cameras are isomorphic to project mappings from a sphere to a plane with a projection center perpendicular to the plane. In other words, a point in three-dimensional space is first projected onto a unit sphere and then onto the imaging plane of a pinhole camera, with the center of the unit sphere at distance $\xi$ from the imaging plane. Thus, the projection model of an omnidirectional camera [19] can be represented as shown in Figure 4.

A point in three-dimensional space ${\left(\chi \right)}_{O_s}={\left(X,Y,Z\right)}^T$(where the subscript $O_s$ denotes the unit spherical coordinate system with $O_s$ as the coordinate origin) is projected onto the unit sphere to obtain the point denoted by ${\chi }_s$.

$${\chi }_s=\frac{\chi }{\parallel \chi \parallel }={\left(X_s,Y_s,Z_s\right)}^T$$

(1)

Then, ${\chi }_s$ on the unit sphere is transformed into the imaging plane coordinate system, with $O_p$ as the coordinate origin.

$${\chi }_s={\left(X_s,Y_s,Z_s+\xi \right)}^T$$

(2)

Subsequently, ${\chi }_s$ is projected onto the normalized imaging plane ${\pi }_u$ to obtain the point $m_u$.

$$m_u={\left(\frac{X_s}{Z_s+\xi },\frac{Y_s}{Z_s+\xi },1\right)}^T={\left(x_u,y_u,1\right)}^T=\Gamma \left({\chi }_s\right)$$

(3)

Taking into account the radial and tangential distortion in the projection process, the point $m_d={\left(x_d,y_d,1\right)}^T$ is obtained after the distortion of the point $m_u$.

$$\begin{array}{c}{x}_d=x_u\left(1+k_1{r}^2+k_2{r}^4\right)+2p_1{x}_u{y}_u+p_2\left(r^2+2{x_u}^2\right)\hfill \\ y_d=y_u\left(1+k_1{r}^2+k_2{r}^4\right)+p_1\left(r^2+2{y_u}^2\right)+2p_2{x}_u{y}_u\hfill \end{array}$$

(4)

where $r^2={x_u}^2+{y_u}^2$, $k_1$ and $k_2$ are the radial distortion coefficients, and $p_1$ and $p_2$ are the tangential coefficients. The distortion matrix D can be expressed as $D=\left[k_1,k_2,p_1,p_2\right]$.

$m_d$ is transformed into the pixel plane ${\pi }_p$ using the camera intrinsic matrix K to obtain the point p.

$$p={\left(u_p,v_p,1\right)}^T=K{m}_d$$

(5)

$$K=\left[\begin{array}{ccc}{f}_x& 0& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right]$$

(6)

The above is the omnidirectional camera image processing procedure under the uniform projection model. To rectify an omnidirectional image to obtain a normal image, it is necessary to invert $m_u$ to ${\chi }_s$.

$${\Gamma }^{-1}\left(m_u\right)={\chi }_s=\left[\begin{array}{c}{X}_s\\ Y_s\\ Z_s\end{array}\right]=\left[\begin{array}{c}\frac{\xi +\sqrt{1+\left(1-{\xi }^2\right)\left({x_u}^2+{y_u}^2\right)}}{{x_u}^2+{y_u}^2+1}{x}_u\\ \frac{\xi +\sqrt{1+\left(1-{\xi }^2\right)\left({x_u}^2+{y_u}^2\right)}}{{x_u}^2+{y_u}^2+1}{y}_u\\ \frac{\xi +\sqrt{1+\left(1-{\xi }^2\right)\left({x_u}^2+{y_u}^2\right)}}{{x_u}^2+{y_u}^2+1}-\xi \end{array}\right]$$

(7)

The rectified coordinates can be obtained by projecting ${\chi }_s$ onto the normalized imaging plane.

$$\left[\begin{array}{c}\frac{X_s}{Z_s+\xi }\\ \frac{Y_s}{Z_s+\xi }\\ \frac{Z_s}{Z_s+\xi }\end{array}\right]=\left[\begin{array}{c}{x}_u\\ y_u\\ 1-\xi \frac{{x_u}^2+{y_u}^2+1}{\xi +\sqrt{1+\left(1-{\xi }^2\right)\left({x_u}^2+{y_u}^2\right)}}\end{array}\right]$$

(8)

The rectification process outputs the remapping parameters of the prerectification and postrectification images. Using these remapping parameters, the rectified image of the omnidirectional image can be obtained. The result of the omnidirectional image correction is shown in Figure 5. Similarly, we can transform the image back into an omnidirectional image after sky region detection.

2.1.2. Sky Region Detection

There is a clear visual difference between the sky and building parts in the images captured by a camera in outdoor scenes. The gradient information of the image can accurately reflect the difference between the sky and the non-sky regions in most cases. Therefore, the operation of the image with the sky region starts from converting the color image into a greyscale image. After the input image has been converted to grayscale image, the corresponding gradient image is calculated using the Sobel operator.

The sky border position function $\mathrm{border}\left(x\right)$ is defined as

$$1\le \mathrm{border}\left(x\right)\le H(1\le x\le W)$$

(9)

where W and H represent the width and height of the image, respectively, and the sky and non-sky regions are defined as follows:

$$\mathrm{sky}=\left\{\right(x,y)\mid 1\le x\le W,1\le y\le \mathrm{border}(x\left)\right\}$$

(10)

$$\mathrm{non-sky}=\left\{\right(x,y)\mid 1\le x\le W,\mathrm{border}(x)\le y\le H\}$$

(11)

Threshold t is used to calculate the sky boundaries in accordance with Formulas (10) and (11). In the gradient image, the first pixel in each column with a gradient value greater than threshold t is defined as the dividing point between the sky and non-sky regions in that column. The sky boundary function $\mathrm{border}$ (x) is preliminarily obtained by traversing each column of the gradient image. The Algorithm 1 shows the algorithm for calculating the sky boundary function border (x) from threshold t. However, choosing threshold t to obtain the optimal sky boundary function is a complex problem.

Algorithm 1 Calculating the sky boundary function border (x)

Input:t, $grad$ Output:$\mathrm{border}\left(x\right)$ 1: forx = 1 to W do 2:    $\mathrm{border}\left(x\right)$ = H 3:    for y = 1 to H do 4:      if $grad(y,x)>t$ then 5:         $\mathrm{border}\left(x\right)$ = y 6:         break 7:      end if 8:    end for 9: end for

In [20], an energy function was first defined to optimize the boundary. In this paper, an improved energy function is used, which is as follows:

$$J_n=\frac{1}{\gamma ·\left|{\Sigma }_s\right|+\left|{\Sigma }_g\right|+\gamma ·\left|{\lambda }_1^s\right|+\left|{\lambda }_1^g\right|}$$

(12)

where ${\Sigma }_s$ and ${\Sigma }_g$ are the covariance matrices of the pixels that are described by the red-green-blue (RGB) values in sky and non-sky regions, respectively. $\gamma$ is a parameter for the homogeneity in the sky region, which is set to 2. ${\lambda }_1^s$ and ${\lambda }_1^g$ correspond to the largest eigenvalues of the two matrices. ${\Sigma }_s$ and ${\Sigma }_g$ are defined as follows:

$${\Sigma }_s=\frac{1}{N_s}\sum _{(y,x)\in sky}\left(I^s(y,x)-{\mu }_s\right){\left(I^s(y,x)-{\mu }_s\right)}^T$$

(13)

$${\Sigma }_g=\frac{1}{N_s}\sum _{(y,x)\in non-sky}\left(I^g(y,x)-{\mu }_g\right){\left(I^g(y,x)-{\mu }_g\right)}^T$$

(14)

where $N_s$ and $N_g$ are the numbers of pixels in the sky region and non-sky region, respectively; I is the input RGB image; and ${\mu }_s$ and ${\mu }_g$ are vectors that represent the average RGB values in the sky region and non-sky region, respectively.

For a given threshold t, a corresponding energy function $J_n\left(t\right)$ can be obtained. The sky boundary that maximizes the energy function is the optimal boundary. However, the relationship between the threshold t and the energy function $J_n\left(t\right)$ is nonlinear, and selecting threshold t so as to maximize the energy function $J_n\left(t\right)$ is a complex problem. Therefore, we use a step search method to obtain the optimal threshold t and calculate the sky boundary function $\mathrm{border}\left(x\right)$. Using the sky boundary function, we can then extract the sky region, as shown in Figure 6b, in which the red part is the extracted sky region. For subsequent NLOS signal rejection, we binarize the image, as shown in Figure 6c, in which the sky region is represented in white, and the non-sky region is represented in black.

However, this sky region detection algorithm cannot solve the problem that flyovers block the field of view of the camera in urban environments, because Formula (10) assumes that the sky region is always present at the top of the image. A rectified omnidirectional image taken when the vehicle is passing under a flyover is shown in Figure 7a. It can be seen that in this case, the assumption that the sky region lies at the top of the image is no longer satisfied. The results of the sky region detection using the algorithm described above are shown in Figure 7b,c. It can be seen that the algorithm mistakes the boundary between the flyover and the sky for the boundary between the sky and non-sky regions, resulting in completely incorrect sky region detection, which will lead to errors in NLOS signal rejection.

In this paper, an improved algorithm is introduced to solve this problem. The HSV color model is used to detect whether the extracted sky region is instead a flyover region. The HSV model is another popular color model in addition to BGR, where H is the hue, S is the saturation, and V is the value, which is more in line with people’s description of color. The range of black in the HSV color model is $0\le H\le 180$, $0\le S\le 255$, and $0\le V\le 46$. After detecting all columns in row $H/6$ of the image, if at least fifty percent of the pixels are considered to be black, it can be assumed that the image contains a flyover region. When a flyover region is detected in the image, the sky region is still determined in accordance with the procedure of the above algorithm; however, the detected sky region is actually a flyover. The BGR values in the flyover region are all set to [255, 255, 255], and the flyover region is temporarily treated as part of the sky region, as shown in Figure 8a. Then, the sky region is detected again, and the results are shown in Figure 8b. Using the boundary function of the flyover region detected during the first execution of the algorithm, the BGR values in the flyover region are set to [0, 0, 0] to obtain the real sky region, as shown in Figure 8c.

2.1.3. Sky Boundary Smoothing

In an urban canyon environment, the close buildings are relatively flat in the omnidirectional camera shots, and the corresponding results of the sky region segmentation algorithm are satisfactory. In contrast, the edges of distant buildings and high-rise buildings are greatly affected by sky reflections, resulting in less clear boundaries. The sky region segmentation algorithm tends to misidentify these reflective regions as sky regions, producing narrow interruptions and slender gullies, resulting in discontinuous sky boundaries, as shown in Figure 9a. The result of binarization after converting back to the omnidirectional image is shown in Figure 9b.

To solve this problem, an image morphology algorithm is employed to process the omnidirectional image after sky region segmentation to make the sky boundary smoother. First, we conduct an erosion operation to restore the parts of the image misidentified as belonging to the sky region to the non-sky region and eliminate the thin interruptions between buildings; the result is shown in Figure 10a. After the erosion operation, the narrow interruptions and slender gullies caused by the influence of strong light are restored to the non-sky regions. However, the erosion operation also results in some loss of the sky regions. Therefore, we subsequently perform a dilation operation to restore the part of the sky region lost due to the erosion operation, with the result shown in Figure 10b. Compared with the original image, as shown in Figure 10c, we can see that after smoothing, our algorithm can well identify the parts of the non-sky region that are affected by strong light.

2.2. NLOS Signal Rejection

The proposed system involves four reference frames: the camera frame, image frame, IMU frame, and GNSS receiver frame. Since the GNSS satellites are far from the ground user, the distances between the GNSS receiver antenna, omnidirectional camera, and IMU are negligible from the perspective of the satellites; therefore, we consider these three frames to be the same frame, and use the omnidirectional camera frame to uniformly represent them.

The pseudorange $pr$, elevation angle $el$ and azimuth angle $az$ of a satellite can be calculated from the received GNSS observations. A satellite S can be represented by $\left(pr,el,az\right)$. Formula (15) is used to represent the satellite in the ENU frame, $S_e=\left(x,y,z\right)$. The x, y, and z axes of the ENU frame point in the east, north, and up directions, respectively.

$$\begin{array}{cc}\hfill & x=prcos\left(el\right)cos\left(az\right)\hfill \\ \hfill & y=prcos\left(el\right)sin\left(az\right)\hfill \\ \hfill & z=prsin\left(el\right)\hfill \end{array}$$

(15)

Then, the satellite is transformed from the ENU frame into the omnidirectional camera frame using the Euler angles $\left(\psi ,\phi ,\theta \right)$ collected from the IMU, where $\psi$ is the roll angle, which represents the rotation angle around the x axis; $\phi$ is the pitch angle, which represents the rotation angle around the y axis; and $\theta$ is the yaw angle, which represents the rotation angle around the z axis. These angles are used to define the rotation matrices around each axis of the system. Formulas (16)–(18) describe the rotations around the x axis, y axis, and z axis, respectively:

$$R\left(\psi \right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& cos\psi & -sin\psi \\ 0& sin\psi & cos\psi \end{array}\right]$$

(16)

$$R\left(\phi \right)=\left[\begin{array}{ccc}cos\phi & 0& sin\phi \\ 0& 1& 0\\ -sin\phi & 0& cos\phi \end{array}\right]$$

(17)

$$R\left(\theta \right)=\left[\begin{array}{ccc}cos\theta & -sin\theta & 0\\ sin\theta & cos\theta & 0\\ 0& 0& 1\end{array}\right]$$

(18)

Thus, the satellite can be represented in the omnidirectional camera frame using Formula (19):

$$S_c=R\left(\theta \right)R\left(\phi \right)R\left(\psi \right)S_e$$

(19)

Since small changes in the roll and pitch angles can be considered constant values to the satellite when the vehicle is moving, we use only the yaw angle for the transformation from the ENU frame to the omnidirectional camera frame. Therefore, Formula (19) can be simplified to Formula (20):

$$S_c=R\left(\theta \right)S_e$$

(20)

After representing the satellites in the omnidirectional camera frame, each satellite in the skyplot at the time of omnidirectional image acquisition can be projected onto the omnidirectional image using the imaging principle of the omnidirectional camera. Then, we check whether each satellite is located in a sky region, and we consider a satellite to be a LOS satellite if it is in a sky region or an NLOS satellite if it is in a non-sky region.

2.3. GNSS Kinematic Positioning Based on Extended Kalman Filter

2.3.1. GNSS Measurement Model

GNSS kinematic positioning in urban environments mainly relies on the pseudorange, which is the basic measurement of a GNSS receiver. The pseudorange is modeled as the sum of the Euclidean distance between the receiver and the satellite, the receiver clock error, the ionospheric and tropospheric errors, and the multipath error [21], as follows:

$$p_{r,j}^s={\rho }_s^r+c·(\delta t_r-\delta t^s)+\delta T_r^s+\delta I_{r,j}^s+DC{B}_{r,j}^s+M_{r,j}^s+{ϵ}_{p,j}$$

(21)

where we have the following:

• ${\rho }_r^s$ is the real distance between the satellite’s antenna and the receiver’s antenna (in meters);
• $\delta t_r$ is the receiver clock error (in seconds);
• $\delta t^s$ is the satellite clock error (in seconds);
• $\delta T_r^s$ is the tropospheric delay (in meters);
• $\delta I_{r,j}^s$ is the ionospheric delay (in meters);
• $DC{B}_{r,j}^s$ is the differential code bias (in seconds);
• $M_{r,j}^s$ is the multipath error for code observation (in meters);
• ${ϵ}_{p,j}$ is the measurement noise for code observation (in meters);
• c is the speed of light in vacuum (in m/s).

The superscript ‘s’ refers to the s-th satellite, and the subscript ‘j’ refers to the j-th frequency band.

We can deduce from Formula (21) that the raw measurement contains various errors that need to be eliminated for an accurate geometric distance measurement. The effect of the tropospheric delay will be eliminated based on ephemeris data from the International GNSS Service (IGS) according to the Hopfield/Saastamoinen model, while the Klobuchar model will be applied to correct the effect of the ionospheric delay, also based on ephemeris data from the IGS. Clock products from the IGS are applied in the correction of the satellite clock bias, while the receiver clock bias is estimated using an extended Kalman filter. The multipath errors and measurement noise are still contained in the pseudorange, affecting the ranging accuracy, which is one of the error sources in the pseudorange measurement.

2.3.2. Extended Kalman Filter

An extended Kalman filter (EKF) [22] is employed to process the nonlinear GNSS measurement model. Two equations can be used to represent the state space model of the EKF: the state equation,

$$x_k=f\left(x_{k-1}\right)+w_k$$

(22)

and the measurement equation:

$$z_k=h\left(x_k\right)+v_k$$

(23)

where ${\mathbf{x}}_k\in {\mathbb{R}}^{n\times 1}$ is the state vector at time $t_k$, ${\mathbf{z}}_k\in {\mathbb{R}}^{m\times 1}$ is the measurement vector, ${\mathbf{F}}_k\in {\mathbb{R}}^{n\times n}$ is the system transition matrix from epoch $t_k$ to $t_{k+1}$, ${\mathbf{H}}_k\in {\mathbb{R}}^{m\times n}$ is the measurement matrix, ${\mathbf{w}}_k\in {\mathbb{R}}^{n\times 1}$ is the system noise vector, and ${\mathbf{v}}_k\in {\mathbb{R}}^{n\times 1}$ is the measurement noise vector. The extended Kalman filter assumes that the system noise and measurement noise are Gaussian white noise and uncorrelated. That is,

$$\begin{array}{cc}\hfill E\left(w_k{w}_l^T\right)& =Q_k{\delta }_{kl}\hfill \\ \hfill E\left(v_k{v}_l^T\right)& =R_k{\delta }_{kl}\hfill \\ \hfill E\left(w_k{v}_l^T\right)& =0\hfill \end{array}$$

(24)

In the GNSS kinematic positioning problem, only the measurement equation is nonlinear, while the state equation remains linear. Therefore, the state equation can be expressed as follows:

$$x_k=F_k{x}_{k-1}+w_k$$

(25)

Then, observation ${\mathbf{z}}_k$ is incorporated to refine the state prediction, updating the state estimate and covariance, as follows.

Prediction is

$$\begin{array}{cc}\hfill {\widehat{x}}_k^{-}& =f\left({\widehat{x}}_{k-1}\right)\hfill \\ \hfill P_k^{-}& =F_{k-1}{P}_{k-1}{F}_{k-1}^T+W_k\hfill \end{array}$$

(26)

Updating is

$$\begin{array}{cc}\hfill K_k& =P_k^{-}{H}_k^T{[V_k+H_k{P}_k^{-}{H}_k^T]}^{-1}\hfill \\ \hfill {\widehat{x}}_k& ={\widehat{x}}_k^{-}+K_k[z_k-h\left({\widehat{x}}_k^{-}\right)]\hfill \\ \hfill P_k& =(I-K_k{H}_K)P_k^{-}\hfill \end{array}$$

(27)

2.3.3. Receiver Motion Model

The accuracy of the motion trajectory estimation using the EKF depends on the accuracy of the description of the vehicle’s motion state. Based on this consideration, we adopt a position–velocity (PV) model. In the PV model, the state vector can be expressed as follows:

$$x_k\triangleq {\left[x,\dot{x},y,\dot{y},z,\dot{z},r_{t_r}\right]}_k^T$$

(28)

A proper discrete time model can be found in [23], and the state transition equation is as follows:

$$F_k=\left[\begin{array}{ccccccc}1& T& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 1& T& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 1& T& 0\\ 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 1\end{array}\right]$$

(29)

The measurement matrix is derived from the Jacobian matrix of $h\left(x\right)$. For sequential Kalman filtering, the measurement matrix reflects the relationship between the receiver position and the individual satellite observations, and the Jacobian matrix can be expressed as follows:

$$H_k={\left[\begin{array}{ccccccc}\frac{\partial p_{r,j}^S}{\partial x}& 0& \frac{\partial p_{r,j}^S}{\partial y}& 0& \frac{\partial p_{r,j}^S}{\partial z}& 0& \frac{\partial p_{r,j}^S}{\partial \delta t_r}\end{array}\right]}_k$$

(30)

The partial derivatives that appear above are shown in detail below:

$$\begin{array}{cc}\hfill \frac{\partial p_{r,j}^S}{\partial x}& =\frac{-\left(X_i-x\right)}{{\left[{\left(X_i-x\right)}^2+{\left(Y_i-y\right)}^2+{\left(Z_i-z\right)}^2\right]}^{1/2}}\hfill \\ \hfill \frac{\partial p_{r,j}^S}{\partial y}& =\frac{-\left(Y_i-y\right)}{{\left[{\left(X_i-x\right)}^2+{\left(Y_i-y\right)}^2+{\left(Z_i-z\right)}^2\right]}^{1/2}}\hfill \\ \hfill \frac{\partial p_{r,j}^S}{\partial z}& =\frac{-\left(Z_i-z\right)}{{\left[{\left(X_i-x\right)}^2+{\left(Y_i-y\right)}^2+{\left(Z_i-z\right)}^2\right]}^{1/2}}\hfill \\ \hfill \frac{\partial p_{r,j}^S}{\partial \delta t_r}& =1\hfill \end{array}$$

(31)

where $(x,y,z)$ is the position of the receiver and $(X_i,Y_i,Z_i)$ is the position of the satellite in the Earth-centered, Earth-fixed coordinate.

More details of the process noise covariance matrix can be found in [23].

3. Experimental Results and Discussions

We evaluated our system on the public UrbanLoco dataset [24]. This dataset was chosen because it is a relatively new public dataset, and other public datasets do not specifically target at the mapping and localization tasks in dense urban applications. The UrbanLoco dataset is a mapping/localization dataset collected in highly urbanized environments with a full sensor suite. The dataset includes a wide variety of urban terrains: urban canyons, bridges, sharp turns, etc. More importantly, the dataset collected in Hong Kong provides information from an omnidirectional camera, an IMU, and a GNSS receiver. In addition, this dataset provides an accurate ground truth that can be easily utilized to evaluate our method. The data collection platform used in Hong Kong, which is shown in Figure 11, was equipped with a Grasshopper3 5.0 MP (GS3-U3-51S5C-C) camera with a Fujinon FE185C057HA-1 fisheye lens (185 HFOV, 185 V-FOV, 10 Hz), an Xsens MTi 10 IMU (100 Hz), and a u-blox M8T GNSS receiver (GPS/BeiDou, 1 Hz). Based on the Hong Kong UrbanLoco dataset, we evaluated the proposed method at three levels, i.e., sky region segmentation, NLOS signal rejection, and urban positioning.

3.1. Evaluation of Sky Region Segmentation

Figure 12 shows the results of the sky region segmentation algorithm when applied to omnidirectional images in the UrbanLoco dataset, which depict several typical urban environments. The first column shows the original omnidirectional images, the second column shows the rectified images, the third column shows the images processed by the sky region segmentation algorithm, the fourth column shows the images binarized from the third column, and the fifth column shows the results of image morphology processing after conversion to omnidirectional images. As shown in Figure 12, the proposed algorithm can effectively segment omnidirectional images into sky and non-sky regions. In the first row, the boundary between the sky and non-sky regions is obvious, but the presence of streetlights and billboards produces some minor misidentifications. In the second and third rows, there are a large number of bright areas on the surfaces of buildings with colors similar to the sky, and reflective windows are present on the surfaces of buildings in the fourth row, but our algorithm can correctly detect these difficulties and achieve satisfactory segmentation results. In the fifth row, the algorithm works successfully in the case of some thin buildings. There is a flyover in the field of view of the camera in the sixth row, and our algorithm can also accurately detect the sky region in this case.

Since the true values of the sky regions are not provided in the UrbanLoco dataset, to quantitatively evaluate the accuracy of the sky region segmentation algorithm, we manually labeled the sky regions in the experimental images. The sky region of an image is colored in green, as shown in the first column of Figure 13. Because the boundaries could not be definitely distinguished, we referred to the corresponding locations on the Google Street View map to ensure the accuracy of manual sky region segmentation. After manual sky region segmentation, the images were binarized and used as benchmarks to evaluate the proposed sky region segmentation algorithm, as shown in the second column. The third column shows the results of our sky region segmentation algorithm, where the red pixels represent the parts of the sky region that are incorrectly detected by our algorithm. The results were obtained by comparing each pixel in the results of our algorithm with the corresponding benchmark image. To calculate the accuracy of the sky region segmentation algorithm more accurately, we did not include the black areas in the original omnidirectional images in the calculation. The sky region segmentation accuracy $P_c$ was calculated using Formula (32), where $cnt\_correct$ is the number of correct pixels; r is the radius of the omnidirectional image, which is 220; and the image size is 612 × 512. We refer to the four images in Figure 13, from top to bottom, as image 1, image 2, image 3, and image 4. The accuracy on these four images is summarized in

Table 1. Accuracy of sky region segmentation on four images.

Image	${\mathit{P}}_{\mathit{c}}$
Image 1	97.05%
Image 2	98.68%
Image 3	97.57%
Image 4	96.89%

.

$$P_c=\frac{cnt\_correct}{\pi \ast r^2}\times 100\%$$

(32)

Figure 13 and Table 1 show that the proposed sky region segmentation algorithm achieves over 96% accuracy in sky region segmentation and can effectively extract sky regions from omnidirectional images acquired in various urban environments without requiring image samples for training.

3.2. Evaluation of NLOS Signal Rejection

The satellite skyplot at a certain time is shown in Figure 14a. Here, ’G’ stands for GPS satellites, and ’C’ stands for BDS satellites. We used the method proposed in Section 2.2 to project each satellite onto the omnidirectional image after sky region segmentation to obtain the corresponding NLOS signal rejection result, as shown in Figure 14b, where green circles represent LOS satellites and red circles represent NLOS satellites.

To verify the accuracy of NLOS signal rejection, we used a 3D photo-realistic model of Hong Kong to derive the satellite visibility. We imported the downloaded 3D photo-realistic model of typical urban environments into Rhino for 3D modeling. From this 3D photo-realistic model, we could easily obtain the true north direction, and each satellite was plotted in the 3D realistic model based on its elevation and azimuth angles. Instead of using the actual distance of the satellite from the ground, we used a fixed distance for observation purposes. Once the satellites had been modeled, the LOS satellites could be clearly distinguished from the NLOS satellites based on the 3D photo-realistic model. The results of the modeling are shown in Figure 15. In the rendering mode, as shown in Figure 15a, we could see whether the 3D photo-realistic model we imported matched the scene captured by the omnidirectional camera. In the coloring mode, as shown in Figure 15b, we marked the LOS satellites in green and the NLOS satellites in red. A comparison with the results in Figure 14b shows that our NLOS signals rejection results are exactly correct. The red lines are discontinuous because the NLOS satellites are obscured by buildings.

Similarly, we verified the NLOS signal rejection performance for the other two images. The results are shown in Figure 16. Among them, Figure 16a,b show the results of the image taken when the vehicle was at a standstill waiting for a turn. Figure 16c,d show the results of the image taken when it was driving through a narrow intersection. From these results, it can be concluded that the NLOS signal rejection algorithm is quite effective in various urban environments.

3.3. Evaluation of Positioning Results

To further validate the vision-aided NLOS rejection method, we selected the HK-Data20190426-2 UrbanLoco dataset, as it contains typical urban environments, such as narrow lanes, half-open environments, and intersections. The complete driving environment and trajectory are shown in Figure 17, where the vehicle started in a narrow lane in area A, and the vehicle in areas C and D stopped waiting for a turn. When collecting the HK-Data20190426-2 dataset, the vehicle traveled a distance of 0.7 km and had a travel time of 260 s, including time spent waiting for traffic lights.

We conducted code-based kinematic positioning experiments with three data screening methods. Specifically, GNSS/Raw estimates the coordinates straightforwardly without any data screening. GNSS/SNR and ELE employs the classic data screening strategy based on the SNR and elevation. In the classic method, thresholds are set on the SNR and elevation to eliminate suspicious NLOS measurements because NLOS signals usually have a small SNR and elevation [25]. In our experiments, the SNR and elevation thresholds of the classic screening method were set to 20 dB and ${15}^{\circ }$ respectively. GNSS/Vision adopts the proposed vision-aided NLOS rejection algorithm. The positioning results are shown in Figure 18.

As shown in Figure 18, the GNSS positioning accuracy exhibits significant deterioration in the harsh urban environment, which is mostly caused by NLOS signals. The positioning results estimated from the raw measurements without screening contain large errors and deviate from the true trajectory in areas B, C, and D, and the positioning results are quite unreliable in highly urbanized environments. Compared with the direct use of the raw measurements without screening, the classic method enables a certain improvement in positioning accuracy. However, in the harsh urban environment, fixed thresholds on the SNR and elevation cannot ensure the effective separation of NLOS signals from LOS signals. When the thresholds are set too high, some valid satellite signals are eliminated, resulting in an insufficient number of satellites. When the thresholds are set too low, a large number of NLOS signals are still retained in the measurements, leading to considerable errors. For example, in area A, incorrectly eliminating NLOS satellites leads to larger positioning errors. In contrast, the vision-aided screening method effectively eliminates the NLOS signals and avoids large-scale mutations. The positioning results are closer to the real trajectory and offer higher reliability.

Based on the reference trajectory provided in the dataset, the positioning errors in the east (E), north (N), and up (U) directions as well as the 3D positioning error are plotted in Figure 19 for the three positioning solutions. In parts of the urban canyon, the obscuration was too serious to observe adequate normal satellites. Conversely, the number of NLOS signals was greatly increased, resulting in incorrect results in some epochs. The vision-aided screening method avoided some large fluctuations occurring in the positioning solutions without screening, contributing to smoother and more reliable solutions. We also generated boxplots to represent the error distributions for the comparison of the three methods, as shown in Figure 20. It is obvious from this figure that the vision-aided screening method has great advantages in improving the results in all directions. The vision-aided screening method effectively reduces the occurrence of large errors while keeping the median error at a low value. It can also be seen from this figure that the classic screening method leads to larger errors at certain times due to the incorrect removal of LOS satellites and retention of NLOS satellites. The vision-aided method can more accurately remove NLOS satellites, thus demonstrating the superiority of our method.

Table 2. RMS positioning results with different methods.

Method	East (m)	North (m)	Up (m)	3D (m)
GNSS/Raw	15.16	12.55	49.48	53.25
GNSS/SNR&ELE	9.39	13.26	31.58	35.51
GNSS/Vision	6.02	11.00	18.17	22.08

lists the root mean square (RMS) values in the E, N, and U directions and the 3D RMS for each solution. Compared with the raw measurements without screening, the vision-aided method enables significant improvements in all directions. In contrast, due to the incorrect removal of LOS satellites at certain times, the classic screening method actually increases the RMS error in the northward direction. The positioning accuracy improvements relative to the raw measurements without screening, expressed as percentages, are shown in

Table 3. Improvements in positioning accuracy expressed as percentages.

Method	East (%)	North (%)	Up (%)	3D (%)
GNSS/SNR&ELE	38.1	−5.6	36.2	33.3
GNSS/Vision	60.3	12.4	63.3	58.5

. Although the classic screening method achieves a certain improvement in the positioning accuracy, the proposed vision-aided screening method improves the positioning accuracy more significantly. This is because the omnidirectional camera captures real information about the obscuration of the surroundings, allowing NLOS signals to be more effectively identified. Specifically, the classic screening method exhibits improvements of 38.1%, −5.6% and 36.2% in the E, N, and U directions, respectively, as well as a 3D accuracy improvement of 33.3%. In comparison, the vision-aided screening method exhibits improvements of 60.3%, 12.4% and 63.3% in the E, N, and U directions, respectively, as well as a 3D accuracy improvement of 58.5%.

Figure 21 shows the time-series changes in the number of satellites available for each of the three methods. The numbers of available satellites are also listed in

Table 4. Numbers of available satellites for the HK-Data20190426-2 dataset.

Satellites	Mean Number	Max Number	Min Number
GNSS/Raw	13.6	18	7
GNSS/SNR&ELE	11.4	16	5
GNSS/Vision	9.2	13	2

. Since the GNSS receiver can receive both GPS and BDS satellite signals, the mean number of satellites received is 13.6. After the exclusion of NLOS satellites using the classic screening method or the vision-aided method, the mean number of satellites drops to 11.4 or 9.2, respectively. In some extreme cases, such as the vehicle passing under a flyover, the vision-aided method culls more satellites and reduces the number of available satellites to fewer than 4. If more satellite navigation systems were to be used, there would be more satellites available that could be used for positioning.

4. Conclusions

As urbanization is becoming increasingly rapid, it is becoming increasingly important to perform accurate positioning in urban environments. In this paper, we propose a GNSS urban positioning method with omnidirectional camera-aided NLOS identification. Our method adopts a novel solution to perform accurate sky region segmentation for omnidirectional images acquired in various urban environments. When a satellite is located in a non-sky region, it is identified as an NLOS satellite. Our proposed method was tested using code-based kinematic positioning on the HK-Data20190426-2 UrbanLoco dataset. The proposed method shows a clear advantage over the classic screening method in terms of the positioning results. The proposed method exhibits improvements of 60.3%, 12.4% and 63.3% in the E, N, and U directions, respectively, as well as a 3D accuracy improvement of 58.5%. In addition, the vision-aided method is simple to calculate and does not require a large number of samples for training. Finally, we conclude that the exclusion of NLOS signals is necessary for achieving more accurate and reliable GNSS urban positioning.

In our work, when NLOS satellites are detected and eliminated, they are not used for positioning, which does not make good use of all the signals. In some extreme cases, this may result in fewer than four usable satellites, and positioning may not be possible. In future work, we will investigate methods of making full use of NLOS signals for more extreme cases.