Detection of GNSS signals propagation in urban canyos using 3D city models

INFORMATION AND COMMUNICATION TECHNOLOGIES AND SERVICES VOLUME: 13 |NUMBER: 1 |2015 |MARCH Detection of GNSS Signals Propagation in Urban Canyos Using 3D City Models Petra PISOVA, Jiri CHOD Department of Telecommunication Engineering, Faculty of Electrical Engineering, Czech Technical University in Prague, Technicka 2, 166 27 Prague, Czech Republic piso[email protected], cho[email protected] DOI: 10.15598/aeee.v13i1.1291 Abstract. This paper presents one of the solutions to the problem of multipath propagation and effects on Global Navigation Satellite Systems (GNSS) signals in urban canyons. GNSS signals may reach a receiver not only through Line-of-Sight (LOS) paths, but they are often blocked, reflected or diffracted from tall buildings, leading to unmodelled GNSS errors in position estimation. Therefore in order to detect and mitigate the impact of multipath, a new ray-tracing model for simulation of GNSS signals reception in urban canyons is proposed - based on digital 3D maps information, known positions of GNSS satellites and an assumed position of a receiver. The model is established and validated using experimental, as well as real data. It is specially designed for complex environments and situations where positioning with highest accuracy is required - a typical example is navigation for blind people. Keywords 3D buildings models, blind navigation, GNSS signal propagation, LOS, multipath, path delay estimation, urban canyons. 1. Introduction The aim of this paper is to describe the calculation of GNSS signal parameters in a complex environment with multiple reflections taking the knowledge of the receiver position in account. This approach can be used to estimate the subject’s unknown position while moving in the same environment. It is well known that multipath is a phenomenon that represents a dominant error source in precise GNSS positioning. In dynamic urban environment, it changes rapidly and therefore it is difficult to detect, predict or control. Generally, the direct (or Line-of-Sight, LOS) signal is normally the most wanted signal. Nevertheless signals can also arrive at a receiver via a number of different paths that may occur between a satellite and a receiver. These paths are results of reflections and diffractions from buildings, water, ground etc. Fig. 1: The real path of the user in urban canyon. Typically, an antenna receives the direct signal and one or more of its reflections. A reflected signal is usually a weaker version of the direct signal and takes more time to reach the receiver than the direct signal. This path delay, or better this difference between length of the path taken by the reflected signal and the direct signal between a satellite and a receiver, causes an important pseudorange measurement error [1], [2] and [3], as shown in Fig. 1 and Fig. 2. Reception of reflected signals also causes distortion in the code correlation peak within the receiver therefore the code phase of the direct LOS signal cannot be precisely determined. c 2015 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 22 INFORMATION AND COMMUNICATION TECHNOLOGIES AND SERVICES VOLUME: 13 |NUMBER: 1 |2015 |MARCH Fig. 2: Illustration of errors in GNSS positioning as computed by standard receivers: in red Mapfactor VTU009 - CSD mode, in yellow Mapfactor VTU009 - GPRS mode, in green MiO 710. Different multipath mitigation techniques were developed to improve positioning performance in built up areas. One promising solution is to use multiple satellite constellations with combination of 3D building models to calculate the path length of direct signals or better path delay of reflected signals and thus to eliminate multipath impact on GNSS signals. There are several satellite systems in operation today - the US NAVSTAR Global Positioning System (GPS), or the Russian GLONASS (Global’naya Navigatsionnaya Sputnikovaya Sistema). Chinese Compass navigation system and the European Union’s Galileo navigation system are currently under construction. 2. Background In order to understand multipath effects it is necessary to understand the electromagnetic wave properties and changes a GNSS signal goes through on the path from a satellite to a receiver. How the position is computed using pseudoranges and the typical multipath error in pseudorange measurements is also explained in this section. 2.1. Electromagnetic Wave Properties A GNSS signal is a Right-Hand Circularly Polarized (RHCP) transverse electromagnetic wave - it means a wave where the electric and magnetic fields oscillate perpendicularly to and in phase with one another, and is perpendicular to the direction of propagation [1], [2] and [4]. Mathematical representation of a plane wave propagating in the direction xcan be written as [5]: ~ E=~ E0cos(kx −ωt +ϕ0),(1) where ~ E0is the amplitude, kis the propagation constant, ωis the circular frequency and ϕ0is the phase constant: k=2π λ,(2) ω=kc =2πc λ,(3) ϕ= (kx −ω+ϕ0),(4) and ϕis the phase of the plane wave. The electric field vector ~ Emay be decomposed into the parallel Eland perpendicular Ercomponents as: El=El0cos (kx −ωt +ϕl0),(5) Er=Er0cos (kx −ωt +ϕr0),(6) with a simple manipulation we get: El El0 = cos (v) cos (ϕl0)−sin (v) sin (ϕl0),(7) Er Er0 = cos (v) cos (ϕr0)−sin (v) sin (ϕr0),(8) where: v=kx −ωt, (9) then we obtain: El El02 +Er Er02 − −El El0Er Er0cos (∆ϕ) = sin2(∆ϕ), (10) when ∆ϕ=nπ/2 (n=±1,±3, ....),El0=Er0=E0, sin (∆ϕ) = ±1and cos (∆ϕ) = 0 then Eq. (10) becomes: E2 l+E2 r=E2 0,(11) and Eq. (11) defines a circle or circularly polarized wave. Electromagnetic waves don’t propagate in straight lines (expect in a vacuum where they travel at 3·10−8m·s−1), they are bent as they pass through the different layers of the earth’s atmosphere to the earth - this increases the GNSS signal travel time from the satellite to the receiver [4]. c 2015 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 23 INFORMATION AND COMMUNICATION TECHNOLOGIES AND SERVICES VOLUME: 13 |NUMBER: 1 |2015 |MARCH Fig. 3: Illustration of GNSS signals propagation in urban environment. 2.2. Propagation of GNSS Signals in Urban Environment The propagation of GNSS signals is mainly impacted of three possible examples - reflection, diffraction and refraction [4], as illustrated in Fig. 3. For our research we expect the GNSS signal to be a straight line. Therefore the proposed algorithm considers reflection and diffraction examples only. Reflection is a process where a wave enters a different medium and its parameters can change at the boundary of the two media (i.e. speed, wavelength, and energy). The Law of Reflection states that the angle of incidence equals to the angle of reflection. Furthermore, the incident ray, the normal to the surface, and the reflected ray all lie in the same plane. As a special case is considered scattering - a wave reflects at all angles (i.e. reflection that doesn’t obey the Law of Reflection). Diffraction is bending of a wave process as it travels past the edge of or around sharps corners (for more look at Huygens’ principle of wave diffraction), [4]. 2.3. Position Computation The position is computed from pseudoranges and satellite positions (Xk,Yk,Zk) found from ephemerides data. The most commonly used algorithm for position computations from pseudoranges is based on the method of least-squares. This method is used when there are more observations than unknowns. A pseudorange is the time the signal travelled from the satellite to the receiver, multiplied by the speed of light [6], [7], [8]. As in [6], let the true geometrical range between satellite kand receiver ibe denoted pk i: pk i=q(Xk−Xi)2+ (Yk−Yi)2+ (Zk−Zi)2,(12) Let cdenote the speed of light, dtibe the receiver clock offset, dtkbe the satellite clock offset, be the tropospheric delay, be the ionospheric delay and be the observational error of the pseudorange. Then the basic equation for the pseudorange Pk iis: Pk i=pk i+cdti−dtk+Tk i+Ik i+ek i.(13) The position of the satellite (Xk,Yk,Zk) and satellite clock offset dtkare found in the ephemerides, the tropospheric delay and the ionospheric delay Tk iand the ionospheric delay Ik iare computed from different a priori models, the error ek iis minimized by using the least-squares method. There are four unknowns in Eq. (13), (X, Y, Z and dt), therefore to compute the position at least four pseudoranges are needed. Before using the leastsquares method, Eq. (13) has to be linearized. 2.4. Multipath Error in Pseudorange Typical multipath error in pseudorange measurements differs from 1 m (in a fair environment) to more than 5 m (in highly reflective environment) and the phase measurement error doesn’t exceed a quarter cycles (if the reflected signal has smaller amplitude than the direct signal), [1]: Consider an antenna receives two signals: a direct signal and a delayed reflected signal with phase shift ∆φand amplitude attenuation α. received_signal =Acos φ+αA cos (φ+ ∆φ).(14) The error in the carrier phase measurement due to the multipath is as: δφ = arctan sin ∆φ α−1+cos ∆φ,(15) in the worst case δφ = 90◦, for α > 1. 3. Urban Modelling Methodology In this section, the urban modelling methodology for creating a 3D model is introduced. At first, we describe how to extract building geometry of a given environment with coordinate data into a file to be able to read data in Matlab. Afterwards, the known coordinates of satellites are obtained from navigation messages. Finally, a conversion between the satellite coordinate system into the mapping coordinate system used for our calculations is presented. c 2015 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 24 INFORMATION AND COMMUNICATION TECHNOLOGIES AND SERVICES VOLUME: 13 |NUMBER: 1 |2015 |MARCH 3.1. 3D Environment Modelling 1) 3D Geometry Extraction To capture and extract 3D building geometry with textures and coordinate data from Google Earth into a .3dr file, 3D Ripper DX was used. The conversion to .obj file was done in Autodesk 3ds Max 2011, and then changed in Mathematica to .csv file. 2) 3D Environment for Plotting and Calculations The environment was simulated in Matlab using .cvs geometry coordinates data, considering the building layout to be made up of triangles (vertices A, B, C) put together. 3.2. Satellite Coordinates The known positions of satellites are calculated from ephemeris data and are found in the navigation message, i.e. in [9]. 1) Conversion of Satellite Coordinates into Mapping Coordinates The obtained satellites coordinates in earthcentred earth-fixed (ECEF) WGS-84 system (Xs,Ys,Zs) are transformed into the mapping coordinates (Xm,Ym,Zm) using the following matrix MECEF _MC : MECEF _MC =   −sin (β) cos β0 −cos (β)·sin (α)−sin (β)·sin (α) cos (α) cos (α)·cos (β) cos (α)·sin (β) sin (α)  , (16) where αand βrepresent latitude (N) and longitude (E) coordinates converted to radians. 4. Mathematical Model In the following, the proposed multipath detection model for dynamic propagation environments is presented. To investigate whether GNSS signals arrived at a receiver directly or through reflections, the raytracing algorithm is introduced based on information obtained from 3D buildings model, assumed position of the receiver and known positions of GNSS satellites extracted from navigation message. In section 4.1, calculation of direct path between a satellite and a receiver is described. In section 4.2, possible reflections of a GNSS signal are calculated, the coordinates of reflection points are determined and thus the reflection length (or the signal path delay) is computed. 4.1. Calculation of Direct Path Consider a ray (vector # » SP) from S(satellite) to P(receiver), and a triangle with vertices A,B,C(obtained from the 3D map), as shown in Fig. 4. To avoid all points (S,P,A,B,C) to be on one line, the algorithm provides an estimate of the number of linearly independent vectors # » AB, # » AC, and # » SP. To get the intersection of the vector # » SP with the ∆ABC, we first determine the intersection of the vector # » SP and a plane (∆ABC ∈ρ)as mentioned in section 3.1. If it does not intersect, then it also does not intersect the ∆ABC and the time and distance from a given satellite to the receiver are calculated from pseudoranges, as in Eq. (12). However, if they intersect in the point R(reflection), we need to determine if this point is inside the ∆ABC. A limitation is done to exclude the intersection point Ras a part of the ray # » SP (R /∈ # » SP). Fig. 4: Illustration of GNSS signal - direct path. 1) Intersection of a Ray with a Plane A plane is defined by three non-collinear points (three points not on a line) A,B,C. These three points define two distinct vectors # » AB and # » AC. The ∆ABC lies in the plane ρthrough Awith normal vector ~n [4], [10]: ~n = # » AB × # » AC. (17) Then the parametric plane equation through a point is: R=A+# » BA·s+# » CA·u, (18) where s,u∈R. c 2015 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 25 INFORMATION AND COMMUNICATION TECHNOLOGIES AND SERVICES VOLUME: 13 |NUMBER: 1 |2015 |MARCH The parametric equation of a ray that passes through a point is: R=P+# » PS·t, (19) where t∈R. By inserting Eq. (19) to Eq. (18), we get the solution to t,s,u: Mt,s,u =   Sx−PxBx−AxCx−Ax Sy−PyBy−AyCy−Ay Sz−PzBz−AzCz−Az  Px−Ax Py−Ay Pz−Az  . (20) By substituting the solution of Eq. (20) into Eq. (18) or Eq. (19), we get the parametric coordinates of the intersection point Rin the plane. 4.2. Calculation of Multipath Reflections 1) Single Reflection - Calculation of Reflection Length Areas with buildings on one side are considered for this model, as shown in Fig. 5. The basic criteria to determine whether a point lies in a triangle or on the edge of a triangle are as [4], [9]: α1+α2+α3≤1,(21) α1, α2, α3≥1,(22) where α1,α2,α3∈R. Fig. 5: Illustration of single reflection on a planar reflector. Given the point R(Rx,Ry,Rz) and ∆ABC, we find the solution to α1,α2,α3: Mα1,α2,α3=    1 1 1 AxBxCx AyByCy AzBzCz  1 Rx Ry Rz     .(23) To determine whether the intersection point Rlies inside or on an edge of the ∆ABC, the algorithm verifies the conditions described in Eq. (21) and Eq. (22). Given by known coordinates of S,R1, and P, the reflection length for a single reflection can be simply calculated as: ldelay =d1+d2,(24) where d1is |PR1|, and d2is |SR1|in meters. From the knowledge that GNSS signals travel at the speed of light and the known distance to the receiver, we can simply calculate the time the signal took to reach the receiver. The benefit of a single reflection is that after reflection the circular polarization of a GNSS signal is changed to left-handed (LHCP). An RHCP antenna can quite effective suppress the LHCP reflection and therefore minimize multipath reflection error. Naturally, a second reflection will cause the RHCP polarization again [11], [12], [13]. 2) Multiple Reflections - Calculation of Reflection Length GNSS availability in urban environment is seriously degraded by buildings blocking direct signals from more than one side, as in Fig. 6. From the knowledge of the single reflection point coordinates (Rx,Ry,Rz), the matrix describing reflections is found. As in Eq. (19), the reflection point Rcan be expressed by: S+~n1·t1=P+~n2·t2,(25) where t1,t2∈R,~n1is the directional vector of the line |SR|,~n2is the directional vector of the line |RP|. The relation between these two directional vectors is given by: ~n2=~n1×MREF ,(26) where MREF is the matrix describing reflection. Let a triangle be defined with vertices Ax,Bx,Cx, (x∈R) and vector ~w to be a normal vector to this plane: ~w = # » AB × # » AC, (27) then the matrix MREF can be expressed as: MREF = −inv   ~w # » AB # » AC  ·  −100 0 1 0 0 0 1  ·  ~w # » AB # » AC  . (28) c 2015 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 26 INFORMATION AND COMMUNICATION TECHNOLOGIES AND SERVICES VOLUME: 13 |NUMBER: 1 |2015 |MARCH Now when we know the reflection matrix, we can simply calculate if a GNSS signal arrived at the receiver via more than one reflection. For the illustration in Fig. 6: Illustration of multiple reflections on planar reflectors. Fig. 6, given by known coordinates of S,R1,R2,R3 and P, the reflection length for multiple reflections can be simply calculated as: ldelay =d1+d2+d3+d4,(29) where d1is |PR3|,d2is |R2R3|,d3is |R1R2|, d4is |PR3|in meters. From the knowledge that GNSS signals travel at the speed of light and the known distance to the receiver, we can calculate the time the signal took to reach the receiver. 5. Experimental Results In this section, the evaluation methodology and results are presented. GNSS availability, integrity and precision were evaluated over pedestrian routes using different combinations of GNSS constellations. Measurements were performed in Europe and Asia for GPS, GLONASS, Galileo and Beidou. In this case, we consider a building block in city model of Prague with variable parameters for street width, street length and building block parameters and known position of the receiver (50◦06010.000 N, 14◦23024.1800 E). Figure 7 shows a street captured in Google Earth used for our measurements as an example, Fig. 8 shows user view of the same street. Fig. 7: Street model of Prague in Google Earth. Fig. 8: Snapshot of user view in Google Earth. Figure 9 shows the same street after 3D geometry extraction into .obj file. Fig. 9: Street model of Prague in .obj file. After extracting the navigation message and converting the satellites coordinates into the mapping coordinates Eq. (16), the algorithm described above is applied. From the well-known position of the GNSS receiver and positions of currently available GNSS satellites in the view, the algorithm determines whether there is an obstacle in direct signal path or not. As a result, we get a list of all direct paths and possible reflection examples (with mapping coordinates of c 2015 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 27 INFORMATION AND COMMUNICATION TECHNOLOGIES AND SERVICES VOLUME: 13 |NUMBER: 1 |2015 |MARCH points of reflection) for each satellite and the receiver, as shown in Tab. 1. and Tab. 2. Tab. 1: Results shown for GPS, satellite Nr. 5. Satellite Nr. 5 Direct Path: No Number of signals with a single reflection: 1 Point of reflection - coordinates (Xm,Ym,Zm): 1.0·10+2 ·[−1.8066,−0.6410,−0.0499] Tab. 2: Results shown for GPS, satellite Nr. 7. Satellite Nr. 7 Direct Path: No Number of signals with a single reflection: 0 Number of signals with two reflections: 3 Signal 1 Point of reflection - coordinates (X1m,Y1m,Z1m): 1.0·10+11 ·[3.7131, −1.9701, 0.1399] Point of reflection - coordinates (X2m,Y2m,Z2m): 1.0·10+2 ·[−1.5763,−0.7414,−0.0099] Signal 2 Point of reflection - coordinates (X1m,Y1m,Z1m): 1.0·10+2 ·[−1.4315,−0.6138,−0.0076] Point of reflection - coordinates (X2m,Y2m,Z2m): 1.0·10+2 ·[−1.4427,−0.5883,−0.0084] Based on known coordinates of points of reflection, the time and the total propagation path length from a given satellite to the receiver are computed using Eq. (24) and Eq. (29). Computation time is extensively increasing with the amount of 3D map information, satellites, and receivers. It is obvious that these facts can be used to compute the subject’s unknown position in a complex environment in reverse order. This will be shown in the following article in the future. The principle of one of the possible methods can be described as a sequence of following steps: •Determine the number of satellites visible using acombination of intelligent maps and known position of the subject. •Perform calculations to determine the position (respectively the position in apresumed region) if the number of satellites that are in direct line-of-sight is greater than or equal to three. •Based on combination of intelligent maps data, their calculations (including terrain and buildings height) and measured GNSS data (with known positions of satellites obtained from navigation message) - i.e. the interval halving method can be used to iterate data in so many steps until an agreement between measured and calculated parameters is reached and thus the determination of the actual position (with sufficient accuracy) is achieved. This method can be combined with other options - i.e. in the first approximation, using an antenna capable of receiving only phase undistorted GNSS signals [14], or cameras able to determine which part of the sky is uninterrupted by objects etc. 6. Conclusion The presence of reflected and diffracted signals in urban environments seriously degrades GNSS positioning in terms of accuracy, integrity and precision. Different combinations of GNSS constellations promise considerable improvements with respect to GPS only. Although, the current situation allows using just GPS and GLONASS fully, in the near future additional satellites of Galileo and Compass will offer significantly improved GNSS performance and availability. Within this paper, a new ray-tracing model for simulation of GNSS signals reception in urban canyons was designed. The proposed multipath detection algorithm is based on digital 3D maps information, known positions of GNSS satellites and an assumed position of a receiver. For each satellite, the algorithm determines whether a signal arrived at a receiver through a direct path or through multipath reflections. Furthermore, the algorithm can estimate the number of multipath reflections and their coordinate data within the proposed simulation system. With this information, the distance a signal travelled to the receiver, as well the transit time, is calculated. Based on the research results presented in this paper, detection of GNSS signals propagation in urban environments can significantly improve accuracy of the pedestrian positioning. In the future work, an investigation in the multi-constellation GNSS area and right selection of direct/multipath signals with one or more reflections for position calculation using the method presented seems to be interesting. It is also expected the data will be sent to the cloud for faster computation. Acknowledgment This work was supported by the TACR grant TA03011396 called Advanced Navigation of Blind, by CTU in Prague, FEE. c 2015 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 28 INFORMATION AND COMMUNICATION TECHNOLOGIES AND SERVICES VOLUME: 13 |NUMBER: 1 |2015 |MARCH References [1] MISRA, P. and P. ENGE. Global Positioning System: Signals, Measurements and Performance. 2nd ed. Lincoln, MA: Ganga-Jamuna Press, 2006. ISBN 0-9709544-1-7. [2] BORRE, M. A., N. BERTELSEN, P. RINDER and S. H. JENSEN. A Software-Defined GPS and Galileo Receiver. A Single-Frequency Approach (Applied and Numerical Harmonic Analysis). Boston: Birkhauser, 2007. ISBN 978-0-81764540-3. [3] GROVES, P. D. Principles of GNSS, Inertial, and Multi-Sensor Integrated Navigation Systems (GNSS Technology and Applications). 2nd ed. Boston: Artech House, 2013. ISBN 978-1-60807005-3. [4] GOTTLIEB, M. A. and R. PFEIFFER. The Feynman Lectures on Physics, Volume I. In: The Feynman Lectures on Physics. [online]. 2013. Available at: http://www.feynmanlectures. caltech.edu/I_toc.html. [5] SOKOLIK, I. N. Properties of Electromagnetic Radiation. Polarization. Stokes’ Parameters. Main Radiation Laws. Brightness Temperature. Emission from Ocean and Land Surfaces. In: Georgia Institute of Technology [online]. 2008. Available at: http://irina.eas.gatech.edu/ EAS_Fall2008/Lecture3.pdf. [6] PARKINSON B. W. and J. J. SPILKER. Global Positioning System: Theory and Applications Volume 1 (Progress in Astronautics and Aeronautics). Washington, DC: American Institute of Aeronautics and Astronautics, 1996. ISBN 156347-106-X. [7] PARKINSON, B. W. and J. J. SPILKER. Global Positioning System: Theory and Applications Volume 2 (Progress in Astronautics and Aeronautics). Washington, DC: American Institute of Aeronautics and Astronautics, 1996. ISBN 156347-107-8. [8] HARPER, N. Server-Side GPS and AssistedGPS in Java. Norwood, MA: Artech House, 2010. ISBN 978-1-60783-985-9. [9] IGS Product Availability. In: International GNSS Service [online]. 2014. Available at: http://igscb.jpl.nasa.gov/ components/prods_cb.html. [10] SUNDAY, D. Intersection of a Ray/Segment with a Triangle. Available In: Geometry Algorithms [online]. 2014. Available at: http://geomalgorithms.com/a06-_ intersect-2.html. [11] PANY, T. Navigation Signal Processing for GNSS Software Receivers. Norwood, MA: Artech House, 2010. ISBN 978-1-60807-027-5. [12] PRASAD, R. and M. RUGGIERI. Applied Satellite Navigation Using GPS, GALILEO and Augmentation Systems. Boston: Artech House, 2005. ISBN 1-58053-814-2. [13] EL-RABBANY, A. Introduction to GPS. The Global Positioning System. 2nd ed. Norwood, MA: Artech House, 2006. ISBN 978-1-59693-016-2. [14] RAMA, R. B., W. KUNYSZ, R. FANTE and K. MCDONALD. GPS/GNSS Antennas (GNSS Technology and Applications). Boston: Artech House, 2013. ISBN 978-1596931503. About Authors Petra PISOVA is a Ph.D. candidate in the Department of Telecommunication Engineering at the Czech Technical University in Prague, Faculty of Electrical Engineering. She received her Master degree in Telecommunication Engineering in 2008 from the CTU in Prague. Her current research interests include the algorithm research and mathematical modelling with a GNSS background. Jiri CHOD joined the Department of Telecommunication Engineering, Faculty of Electrical Engineering, CTU in Prague in 1969. He obtained his Ph.D. in 1984 and became an associate professor in 1991. He is a long time supervisor of doctoral students and an investigator of a number of grants and projects. His scientific interests include the mobile communications, microprocessors and e-health technologies. c 2015 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 29