Coordinate Frames and Transformations in GNSS Ray-Tracing for Autonomous Driving in Urban Areas

: 3D Mapping-Aided (3DMA) Global Navigation Satellite System (GNSS) is a widely used method to mitigate multipath errors. Various research has been presented which utilizes 3D building model data in conjunction with ray-tracing algorithms to compute and predict satellites’ visibility conditions and compute delays caused by signal reﬂection. To simulate, model and potentially correct multipath errors in highly dynamic applications, such as, e.g., autonomous driving, the satellite–receiver–reﬂector geometry has to be known precisely in a common reference frame. Three-dimensional building models are often provided by regional public or private services and the coordinate information is usually given in a coordinate system of a map projection. Inconsistencies in the coordinate frames used to express the satellite and user coordinates, as well as the reﬂector surfaces, lead to falsely determined multipath errors and, thus, reduce the performance of 3DMA GNSS. This paper aims to provide the needed transformation steps to consider when integrating 3D building model data, user position, and GNSS orbit information. The impact of frame inconsistencies on the computed extra path delay is quantiﬁed based on a simulation study in a local 3D building model; they can easily amount to several meters. Differences between the extra path-delay computations in a metric system and a map projection are evaluated and corrections are proposed to both variants depending on the accuracy needs and the intended use.


Introduction
Quality requirements for positioning in terms of accuracy, continuity and integrity are very stringent for urban navigation applications, such as autonomous driving [1][2][3].The GNSS sensor is the only navigation system providing absolute positioning, and is, therefore, indispensable in multi-sensor navigation systems.To enable accurate and precise GNSS positioning, various error sources need to be modeled and corrected.The GNSS error sources, such as orbit errors, satellite clock errors or errors from ionospheric and tropospheric refraction, as well as relativistic effects are well-known and can be modeled or estimated accurately in relative positioning approaches [4] and Precise Point Positioning (PPP) [5], respectively.However, in urban areas, multipath and non-line-of-sight (NLOS) signal propagation are remaining error sources.NLOS delays can reach up to twice the orthogonal distance to the reflector plane (e.g., tens of meters for small urban trenches in European inner cities [6,7] and up to hundreds of meters for urban canyons in metropolitan cities [8,9]).
The delay caused by a multipath depends on the extra path delay and the receiver settings and is typically assessed by the multipath envelope function for code observations yielding up to 30-50 m.For the carrier phase, the impact ranges from millimeters to a maximum of a quarter of the wavelength (typically around 5 cm) [10].
In order to mitigate multipath errors, 3DMA GNSS is a widely adopted method [11].Furthermore, 3DMA GNSS can be divided into three specific categories: terrain height aiding, shadow matching and 3DMA ranging techniques.Firstly, terrain height aiding increases the horizontal accuracy in dense urban areas by constraining the positioning solution to a known surface height using a digital terrain model [12,13].Secondly, shadow matching is a GNSS positioning technique which uses 3D building model data to determine the user position by exploiting the satellite visibility prediction [14] and comparing the predicted and measured signal-to-noise ratios [15,16].
In terms of 3DMA ranging techniques, various research has been presented which utilizes 3D building model data in conjunction with a ray-tracing algorithm to predict and detect satellites' visibility conditions (line-of-sight (LOS), multipath, diffraction, NLOS) [17,18] and subsequent delays caused by signal reflection for GNSS multipath error correction [19][20][21].This method is also used to realistically simulate GNSS signals in challenging environments [22].The preceding approaches can be combined to further improve GNSS positioning in urban environments through machine learning [23].
In order to determine accurate extra path delays and subsequent multipath information, the satellite-user-reflector geometry must be given in a common metric coordinate system.However, in practice, this requirement can be violated since different coordinate systems are in use and the 3D building models are often given in a map projection such as Universal Transverse Mercator (UTM), which leads to distortions of distances [24,25].
Several cities across the world provide building models in CityGML using projected coordinate systems which are comparable to the 3D building model data provided by the city of Hannover [26].Examples are the cities of Berlin [27], Amsterdam [28], or New York City [29].In addition, building model data is becoming more accurate, including level of detail (LoD) 3 building information and photo-realistic textures [30] or point clouds from LiDAR sensors [31,32].
Using the precise height information of the GNSS antenna locations, the accurate multipath modeling of ground reflections can be performed.In [33,34], GNSS multipath errors are simulated and validated with simple, geometry-controlled experiments.They present that multipath signal amplitudes are modeled accurately when the antenna gain is given, and the geometry of the antenna environment is known precisely.To simulate, model and potentially correct multipath errors in dynamic applications, e.g., in autonomous driving, the surrounding geometry has to be known precisely even for GNSS code observations, since the error is sensitive w.r.t. the phase shift to the cm-level [35].Hence, inconsistencies between the various coordinate frames lead to falsely determined multipath errors.A sound description of how ray-tracing techniques deal with different coordinate frames, i.e., the coordinate frames of given building model data and user and satellite coordinates, is still lacking.
The remainder of this paper is structured as follows.Firstly, we briefly introduce our ray classification algorithm and the computation of the extra path delays.Secondly, we introduce the coordinate systems and general transformation steps that are needed when using data from different coordinate representations.We illustrate their application in a case study from Hannover, Germany.We present the different data types as well as their coordinate systems and frames used for 3DMA ray-tracing, in our case, and describe the application of the aforementioned steps.Thirdly, we quantify the impact of the inconsistencies between the involved coordinate representations on the computed extra path delay based on a simulation study.In this context, minimizing the differences between the computations of the extra path delay in a metric and a projected representation is attempted.

Ray Classification Algorithm
In 3DMA applications, it is tested whether the ray between a satellite and antenna is free of obstructions by buildings or blocked, reflected or diffracted, as introduced in [36].

LOS/NLOS Classification
The LOS/NLOS classification needs a common coordinate basis, as the ray is spanned between the antenna position A and the satellite position S.The resulting ray is then checked for intersections with the 3D building model.More specifically, each building's polygon needs to be tested for an intersection with the ray.
To increase the algorithm's speed, we firstly determine which buildings are in the vicinity of the ray and, secondly, divide all polygons into triangles to perform a ray-triangle intersection check [7].If any of the polygons are intersected by the ray, the satellite is classified as NLOS.Otherwise, it is classified as LOS.

Determination of Signal Reflections
In order to investigate if a signal is under the influence of multipath or experiences reflections on building surfaces, reflection points need to be determined.This, again, requires the vector of the ray between an antenna and satellite, as well as nearby building surfaces, to be in a common coordinate frame.
In Figure 1, the geometry of the reflection calculation is displayed.To find the reflection point R on the reflecting surface P, the antenna position A is mirrored w.r.t this plane.The projection of the antenna onto the reflection surface K is calculated using the normalized outer normal vector n of the reflection plane and any point P 0 on the plane.The point P 0 can be any vertex, and n can be calculated from three vertices of the surface.Next, the mirrored antenna point A is calculated.
The intersection between the ray from the mirrored antenna point A to the satellite S will give the reflection point R.As a result, the ray can be classified into one of the four categories: 1. LOS: rays without obstruction and without any possible reflection points; 2.
Multipath: satellite is in LOS and additionally a signal reflection is determined; 3.
NLOS: the direct LOS is obstructed by a building but the signal can be received via reflection on a building surface; 4.
Blocked: all other cases in which the direct LOS is obstructed and a single reflection does not suffice for the signal to reach the antenna.
An example of the satellite ray classification of a kinematic trajectory is shown in Figure 2, where the ray-tracing results are displayed in a GNSS Feature Map [21].The y-axis displays the distance travelled.At each epoch, i.e., travelled distance, the respective skyplot is projected onto one line and represented as such, i.e., the eastern part of the hemisphere is represented to the right of the center line while the western part is shown to the left.Figure 2a shows the classification of the visible satellites while travelling six times through the trajectory depicted in Figure 2b.Each almost vertical line is such a recorded and classified track.It can be seen that ray-tracing in combination with the GNSS Feature Map is a useful tool to determine critical parts of a trajectory (street segments are separated by horizontal black lines) in terms of signal propagation characteristics.As an example, in segment B, a multipath is indicated for all satellites in the western part of the hemisphere which are reflected in the buildings to the right of the trajectory.

Extra Path-Delay Computation and Multipath Error
The calculation of the extra path delay δ in Equation ( 4) is based on information about the geometric conditions of a satellite-reflector-antenna scenario published in [37] and demands knowledge about the outer normal vector n of the reflection surface, the position of the antenna A and an arbitrary point P 0 on the reflection surface P, δ = 2 • ∆x t cos α p cos p + ∆y t sin α p cos p + ∆z t sin p • cos (α p − α) cos p cos + sin p sin .(4) The vector ∆x t represents the vector between the antenna and the arbitrary point P 0 in the topocentric system; α p and p denote the azimuth and elevation of the outer normal vector n, respectively; and α and are the azimuth and elevation of the satellite, respectively.
For NLOS conditions, the extra path delay represents the ranging error.In multipath conditions, the extra path delay is input to compute the multipath error.In case of the superposition of only one reflected signal path to the direct path, the multipath error can be expressed as for code ρ MP and carrier phase Φ MP , respectively.In the Equations ( 5), r denotes the amplitude ratio between the direct and reflected ray, which depends mainly on the material properties and antenna gain [8,38,39], and λ the wavelength of the respective carrier.In Ref. [35], it was shown that small variations, by 1 cm of the horizontal distance to the reflection plane, can cause variations in the code multipath error of up to 20 cm for distances to the reflector of up to 3 m.This is due to the oscillating character of the multipath error with periods of the wavelength λ.Thus, small errors in the relative geometry will distort the computed multipath error significantly.

Coordinate Frames
For classifying the GNSS signal propagation condition (such as LOS, multipath, NLOS and blocked) and for computing the extra path delays, the satellite and user position as well as the reflector surface must be given in a common Cartesian coordinate frame; we summarize the respective relations in the following.

Coordinate Frames and Datum Definitions for GNSS Satellite Orbits
Various GNSS use specific Earth-Centered, Earth-Fixed (ECEF) 3D Cartesian reference frames to provide orbit information as broadcast ephemerides to the user.Examples are the realizations of the World Geodetic System 1984 (WGS84) for GPS, the Parametry Zemli-90 (PZ-90) for Global'naja Nawigatsionnaja Sputnikowaya Sistema (GLONASS), the Galileo Terrestrial Reference Frame (GTRF) for Galileo, or the China Geodetic Coordinate System 2000 (CGCS20) for BeiDou.The system providers are making continuous efforts to align the respective frame realization with the current version of the International Terrestrial Reference Frame (ITRF) [40,41].Corresponding transformations can be found in, e.g., Ref. [42].The ITRF reference frame is represented by a set of station coordinates and respective velocities given at a reference epoch, which allows to consider the station movements, e.g., by plate motion, which occurred between the epoch of coordinate definition and the epoch of measurement.The Geodetic reference system 1980 (GRS80) ellipsoid is associated with this frame.
In addition, for post-processing purposes, the precise ephemerides from the IGS Multi-GNSS Experiment (MGEX) project provide multi-GNSS satellite orbits in a common frame, the IGS frame is currently IGS14 [43], whose parameters are aligned to the ITRF.
Although the International GNSS Service (IGS) 14 frame and the ITRF2014 are not technically identical, the differences between them are in the mm range [44].

Coordinate Frames and Datum Definitions for 3D Building Models
Building models are often provided by regional public or private services.The coordinate frames are based on specific regional standards, such as the plate-fixed 3D Cartesian European Terrestrial Reference System 1989 (ETRS89) in Europe.It is derived from the ITRF global reference system, and was consistent with the ITRF at epoch 1989.0.Since it is fixed to the Eurasian tectonic plate and is co-moving with that plate (about 3 cm per year in the north-eastern direction), a drift compared to realizations of the International Terrestrial Reference System (ITRS) is introduced [45].Associated ellipsoid heights, longitudes and latitudes refer to the GRS80.A similar situation can occur with the different realizations of WGS84 [42].
The horizontal coordinates of a 3D building model are provided in the coordinates of a map projection of the ellipsoidal coordinates.This is due to the frequent or even mandatory use of map projections in planning, surveying or cadastre [25], and the subsequent use of projected 2D coordinates, such as Northing and Easting, in UTM.For this purpose, a reference ellipsoid has to be selected.The height component is generally a physical height and, thus, refers to a different reference surface [46].Consequently, for each application, it must be considered whether the height components are physical or ellipsoidal heights and whether respective conversions using suitable geoid undulations are required; see, also, the discussion in [25].
The 3D building model applied in this study is provided by the city of Hannover in LoD 2 using the CityGML format [26,47].The data basis is the building contours from cadastral real-estate records combined with laser-scan point clouds to recognize the roof shape.Therefore, the 2D accuracy results from the accuracy of the building contours (cm-level), whereas the height accuracy of the roofs is given at mostly 1 m, depending on the complexity of the roof [48].The 2D x and y coordinates of the 3D building model are given as the projected UTM values Easting and Northing.The UTM projection relies on the ETRS89 and the European Terrestrial Reference Frame (ETRF) 2000 with its realization of December 1st 2016, which corresponds to the official definition of the German cadastre.The physical height component is given according to the German Main Height Network 2016 (DHHN2016).

Coordinate Frames for User Positions
The coordinates of the user antenna position can be given in different coordinate systems.Typical situations are the determination of the user position by GNSS.In case of absolute position by Single Point Positioning (SPP) or PPP, the user position is given in the frame of the orbits introduced during the analysis.In case of relative positioning, the coordinates can refer to local frames, e.g., when using network real-time kinematic positioning (RTK).For user convenience, some receivers apply coordinate transformations internally.

Transformation between Different Frames
To resolve the geometric relations between the building model, the satellite position and the user position, it is required to convert all data to a common coordinate system at the same epoch.The most intuitive forms of representation would be the representation of the situation in local topocentric 3D Cartesian East, North, Up (ENU) coordinates or the map projecting environment of the building model at the required accuracy level.However, only the metric 3D Cartesian system, which represents the real positional relations between the objects, allows reality-true computations.
The required transformation process is exemplified in Figure 3.

Conversion from Map Projection into a 3D Cartesian System
The steps that have to be considered during the conversion of the available data into one metric system are shown in Figure 3 following the blue colored arrows.Taking the ellipsoid dimensions (i.e., semi-major axis a and flattening f ) into account, the back projection of the UTM Easting and Northing coordinates into ellipsoidal coordinates (ϕ, λ) is applied as described in [49]: In parallel, the particular height properties of the building model have to be considered and eventually converted into an ellipsoidal height.In our case, the heights of the building model are physical heights given in DHHN2016.The latest German official geoid, the German Combined Quasigeoid 2016, allows the conversion between ellipsoidal heights h in the ETRS89 realization of 2016 and physical heights H in the DHHN2016 [50].The model's accuracy is 1 cm in the lowlands, 2 cm in the high mountains and 2 to 6 cm on the shore.Due to the position of Hannover, we can consider an accuracy of 1 cm.
Using the geoid undulation ξ, the conversion can be computed Next, using the dimensional parameters of the ellipsoid, the ellipsoidal coordinate triplet (latitude, longitude, ellipsoidal height) is converted into 3D ECEF Cartesian coordinates in the ETRF frame at the epoch of definition of the building model, where is the radius of curvature in the prime vertical [51].The eccentricity e reads e 2 = 2 f − f 2 .

Transformation between Plate-Fixed and Earth-Fixed Frames, Datum Transformation
Now, it is possible to apply a transformation between the reference frames, considering the realizations and the respective epochs.The Reference Frame Sub-Commission for Europe (EUREF) provides transformation parameters and transformation rules between the ITRF and the ETRF (cf.Equation ( 12)) as well as between single realizations of the ITRF (cf.Equation ( 13)) [52].Thus, the transformation from the ETRF2000 to the ITRF2014 reads x with the translation vectors t and t yy , the scale factor D, the respective rotation matrix R and the rotation rate matrix Ṙ.The index t refers to the year of the measurement in Gregorian years.Using decimal notation and decimal digits, plate motions within one year can be considered.The respective values are dependent on the realization and the reference frame of the present and the desired datum.The values are determined and provided by the EUREF [45].To realize an accurate transformation between coordinates, it is required to apply first the transformation between ETRF2000 and ITRF2000, and then the transformation between ITRF2000 and ITRF2014 [53].The calculated coordinates are rounded to the fourth decimal place to prevent problems with the precision after the transformation.In addition, the truncation is to prevent inaccuracies when calculating the reflection points by a precision, which is not given in the original 3D building model.For the inverse transformation, first, the frame in the respective realizations of the ITRS has to be matched to the one of the 3D building model.The transformation can be performed inverse to Equations ( 12) and (13).Similar considerations have to be made for the different realizations of the WGS84, see [42].

Representation of Coordinate Differences in a Topocentric System
For the representation of 3D Cartesian-coordinate differences in a local topocentric frame, it is necessary to select the global position of the topocenter, which could be the user position.Depending on the coordinate representation of the user position, it could be required to perform the same conversion steps as for the 3D building model.
Using ellipsoid parameters, the 3D Cartesian coordinates of the origin of the topocenter are converted into ellipsoidal longitude, latitude and height.For this, the iterative formula λ = arctan y e x e , ( 14) − N (j) , ( 15) 1 − e 2 N (j) can be used with p = x 2 + y 2 and the starting value for the iteration Then, 3D ECEF Cartesian-coordinate differences ∆x e , ∆y e , ∆z e are computed w.r.t. the position of the origin of the topocenter.This can be the differences between the user antenna position and the satellite position, as well as the user antenna position and the vertices of the facade of the building modeled transformed to the same ECEF frame.
Finally, the transformation from ECEF to topocentric coordinates reads:

Conversion from 3D ECEF Cartesian Coordinates to a Map Projection
Although only the metric Cartesian 3D system allows realistic calculations of the extra path delay, the calculation of the extra path delay in another system, such as the map projection of the building model, can be useful, e.g., to save computational effort, or when integrating local measurement sensors, such as laser scanners or cameras.Several characteristics have to be taken into consideration to achieve an accuracy that closely approximates the results of the 3D Cartesian calculation, such as the meridian convergence γ, the mapping distortion and a height correction.Distortions due to the curvature of the Earth may be neglected in the calculation of the extra path delay because of the comparatively small extra path delay compared to the Earth curvature.
First, the respective realizations, epochs and reference systems of the respective representations must be considered.Analogous to Section 4.2, the 3D ECEF Cartesian coordinates must be converted from the given realizations into the realization of the map projection.To start the map projection, the 3D coordinates are usually converted into ellipsoidal coordinates first, cf.Equations ( 9) to (11).The height components must be considered separately, and conversion into the reference height of the objects in the map projection must be applied if necessary.The relationship is given in Equation (8).
The satellite position is usually given in 3D ECEF Cartesian coordinates.Converted into topocentric coordinates and taking the meridian convergence γ [54] into account, the lineof-sight unit vector to the satellite can be represented by the elevation and the corrected azimuth, then used in Equation ( 4): with (λ A − λ 0 ) being the longitudinal difference of the investigated position A to the central meridian and ϕ A being the latitude of the investigated position [54].Inaccuracies in the representation of the elevation are neglected.
A mapping correction k is required, due to a horizontal distance distortion resulting from the urge to minimize scale variations in a specific zone of a map projection.For UTM, each zone is leveled with the most favorable transverse Mercator projection.The projection cylinder intersects the surface, and a part of the surface, thus, emerges from the projection.As a result, the meridian in the center of the circles of contact between the cylinder and the Earth's surface is shortened by a factor of 0.9996.This scale factor k is dependent on the geodetic latitude ϕ A and longitude λ A and increases by the distance from the reference meridian.It can be calculated by where k 0 is the scale factor at the reference meridian [55].
Distance distortion due to the height over the reference ellipsoid can be adjusted using the scale factor k sl [56].
where h ell denotes the mean ellipsoidal height of two respective measurement points and R denotes the radius of the Earth.The scale factors can be considered by dividing a respective length in the projection environment by the corresponding scale factors.The distortion correction is applied within the computation of the extra path delay in Equation ( 4), which depends mainly on the relative geometry.The first parenthesis term denotes the horizontal distance between reflection surface and antenna and can be corrected by the scale factors.

Simulation Set Up
The simulation study aims to show the impact of the above-described conversion and transformation steps onto the computed extra path delays: namely, the impacts of (i) frame transformation between plate fixed and ECEF frames, (ii) meridian convergence, and (iii) distortion due to map projection.
The extra path delay is calculated according to Equation ( 4) and in the respective system.The conversion of the user position to the UTM environment and the computation of the extra path delay in the map projection is performed step wise to show the impact of the single steps.The involved coordinate frames in the Hannover use case are showcased in Table 1.In order to cover all incidence angles, the simulated satellites are equally distributed over a 360 • × 90 • grid in 1 • increments from 1 • to 360 • in azimuth and from 1 • to 90 • in elevation.The simulation was conducted on the 274th day of the year 2022.The reflection points are calculated for four antenna positions in front of a 240 m long and 64 m high wall.Its representation in the UTM environment can be seen in Figure 4a.The artificial wall is generated at the location of the university campus and the ellipsoidal coordinates of the closest antenna correspond to 52.2888°North and 9.7127°East.All four antennas are set to be at two meters height and are located orthogonal to the center of the wall.The antenna reflector spacing of the individual antennas is chosen so that they could represent typical situations in densely populated areas.For example, positioning on a pavement, the position of a car on roads of different widths, or relatively distant reflections along an urban canyon.In order to perform the calculations in the metric system, the wall is brought from the UTM coordinates into the metric environment.The steps in Figure 3 are followed.The second computation is performed in UTM in order to highlight the impact of the different corrections and transformation steps, cf. Figure 3. Figure 4b shows the origin of all signals that result in a reflection.

Impact of Neglecting Different Reference Frames and Meridian Convergence
Figure 5 shows the differences in the extra path-delay computation for the scenario where neither the transformation between the different reference systems nor the meridian convergence is applied to the computation in the UTM environment.This indicates that the user position, given in the ITRF coordinate frame, and the 3D building model, given in the ETRF coordinate frame, are not aligned in one common frame.Additionally, the deviation in their north direction of the 3D building model from the geographic north direction is not considered.The differences are shown against the azimuth values of the satellites computed in the metric system and can be compared to the skyplot in Figure 4b.The color scale indicates the incidence angle, which refers to the angle between the ray incident on the surface and the line perpendicular to the reflection plane in the metric system.Higher incidence angles imply more distant reflection points, as the signals are reflected alongside the wall.The metric distances indicate the corresponding antenna reflector spacing.The single Root Mean Square Error (RMSE) for the different antenna reflector spacings can be seen in Table 2 and reach values up to 128 cm.The magnitude of the differences reaches maximum values of 2.13 m.
Table 2. RMSE of ∆δ for different distances between the antenna position and the reflection point after the application of the meridian convergence and the system transformation as well as the additional application of the distortion correction.A dependence of the difference on the incidence angle and, thus, the distance to the reflection point can be seen, especially for a small antenna reflector spacing.However, the symmetrical dependence is skewed, especially at large antenna reflector spacings.The behavior of the differences concerning the azimuth angle of the satellite, as well as the positive offset of the differences, can be explained in comparison with Figure 6a,b.

Transformation to a Common Reference System
Figure 6a shows the differences between both computed extra path delays after the data given in the ITRF reference frame are aligned to the reference frame of the 3D building model.However, the difference between the grid and geographic north direction is not considered.
The impact of the transformation into a common reference system can be recognized by the negative offset visible in Figure 5 compared to Figure 6a.The offset visualizes that δ utm is too small.The difference evolves from the temporal evolution of the ITRS2014 relative to the ETRS89.The distance between the reflection plane and antenna position decreases since the reflection plane is located north of the antenna positions, and the relative shift between both systems results in a north-east shift in the antenna position towards the reflection plane.The conversion to conformal systems ensured the uniformity of the differences.The differences are small for a signal origin perpendicular to the reflection wall, here an azimuth of 180 degrees.For azimuth angles smaller or larger than 180 degrees, i.e., opening incidence angles, the magnitudes of the differences grow approximately uniformly for the respective antenna reflector spacings.The distribution of negative and positive magnitudes over the azimuth angle indicates that the calculated δ utm of signals originating from the west are calculated as too small.At the same time, the δ utm of signals originating from the east are computed as too long.Responsible for that behavior is the neglected meridian convergence.The orientation of the reflection surface with respect to the direction vector of the satellite ray is misaligned.The misalignment results in extended or shortened reflection paths, depending on the relative geometry between the directional signal vector and the orientation of the reflection surface, as can be seen in Equation (4).Subsequently, reflection points that are close to the antenna projection on the reflection surface are not as strongly affected by the meridian convergence and its resulting shift as reflection points at the more eastern or western ends of the surface.

Different Frames but Consideration of the Meridian Convergence
Considering the respective meridian convergence at the satellite's azimuth in the UTM representation, the results are shown in Figure 6b.Based on the origin of the satellite on the azimuth plane, the differences ∆δ are now symmetrically distributed around an orthogonal origin of the signal source.The differences increase in accordance with the decreasing incidence angles.This is due to the additional travel distance of the respective signals being up to twice the antenna reflector spacing, and, subsequently, the bigger impact of errors evolving from the neglected temporal evolution of the ITRS2014 relative to the ETRS89.The magnitude fits to the magnitude of the temporal drift between the reference frames, and is, therefore, independent of the horizontal antenna reflector spacing.This leads to the overlap of the data sets of the different antenna reflector spacings.The negative offset, resulting from the different reference frames, is clearly visible.
Both corrections are able to significantly reduce the overall RMSE, as shown in Table 2.However, the impact of both corrections varies with the antenna reflection spacing.The effect of the aligned reference systems is greatest at shorter spacing, while the effect of the considered meridian convergence increases with the spacing.This is very well visualized by the row-by-row evolution of the single RMSEs for the different antenna reflector spacing in Table 2.

Transformation to Common Reference Systems and Consideration of the Meridian Convergence
Figure 6c shows the results of the combined application of the meridian convergence and frame transformations.The magnitude of the differences can be drastically reduced and the RMSE values improve drastically compared to the use of different frames and no consideration of the meridian convergence.However, a general dependence on the antenna reflector spacing can be seen as ∆δ increases, the further the antenna positions are away from the reflection surface.Similar to Figure 6b, ∆δ increases with decreasing incidence angles for a fixed antenna position.This is due to the varying order of magnitude of the extra path, which depends on the incidence angle.For signals with a low incidence angle, the additional paths are longer, causing the error induced by the horizontal distortion to have a greater effect.Additionally, the distance distortion in the map projection leads to shortened horizontal distances, which is visible by the negative offset of the differences in Figure 6c.

Transformation to Common Reference Systems, Consideration of the Meridian Convergence and Distortion Correction
Figure 6d shows the differences ∆δ after the distortion correction.The magnitude of the difference decreases significantly for the different antenna positions compared to Figure 6c.The maximum magnitude is smaller than 2.3 mm for an antenna position 90 m distant from the reflection surface and smaller than 0.2 mm for closer antenna positions at a maximum distance of 20 m from the reflection surface.However, the positive differences for close antenna reflector spacings show that the distortion correction is too big in case of small antenna reflector spacings.In fact, the individual RMSE value for an antenna reflector spacing of 2.2 m is bigger than for an antenna reflector spacing of 11.2 m.Nevertheless, the individual RMSE values are in the lower single-digit cm range or in the mm range, as displayed in Table 2. Another error source becomes visible when introducing the impact of the component-wise distortion correction mentioned in Section 4. The color scaling in Figure 6d indicates the relation between the vertical and horizontal difference between user position and reflection point in the metric system.A high value corresponds to a large share of the vertical component in the total distance.It coincides with larger differences in ∆δ for the same antenna reflector spacing.The reason is that the distortion correction is applied to the horizontal distance only.Thus, a LOS vector between the user position and the reflection point in the projected environment contains a vertical component which is affected by vertical distortion if they are not parallel to the horizontal plane.

Conclusions
In this paper, 3DMA methods are used to improve GNSS-based positioning and navigation in severe multipath environments such as urban areas.The accuracy of the underlying ray-tracing in classifying the satellite signal (e.g., LOS or NLOS) and estimating extra path delays depends on the consistency of the relative geometry between satellite, receiver antenna and reflector.Due to the different coordinate systems of the 3D building models involved, the GNSS receiver antenna positions and the satellite geometry, the peculiarities may not be obvious in all applications or to all users.
We investigated and clarified the typical frames and definitions associated with the three quantities (3D building models, GNSS receiver antenna positions and satellite geometry) needed for ray tracing.The respective frame and datum transformations as well as corrections are summarized for the transformation from a building model in map coordinates such as UTM to 3D Cartesian topocentric coordinates.As a result, all three quantities are provided in the same Cartesian coordinate frame and the same datum, which is a prerequisite for error-free ray-tracing computations.Since it is computationally expensive to transform the building model in a 3D Cartesian frame, we developed a correction scheme for frame inconsistencies, meridian convergence and map distortions.
Using a simulation study, the impact of the different transformation and correction steps on the computation of the extra path delay is revealed.We showed that the application of the meridian convergence to the directions of the satellite and the frame transformations in case of a plate-fixed frame such as ETRF with a long accumulation time of plate motion w.r.t. an ECEF frame such as ITRF, have a large impact.Depending on the antenna reflector spacing, the error in the extra path delay can reach more than 200 cm for 120 m and more distant reflection points.The impact of distortion due to the UTM projection is at the cm-level.
Taking all frame definitions into account and applying all corrections accordingly, we showed that the extra path delay can also be accurately determined in 2D UTM coordinates completed with an independent height, as typically performed in 3D building models.

Figure 1 .
Figure 1.Reflection of a satellite signal on a building surface P at the reflection point R. The user position is denoted by A, its image w.r.t.P as A .The extra path delay is depicted as the difference in the distances between the satellite S and the mirrored antenna point A , as well as the satellite and antenna A.

Figure 2 .
Figure 2. Visualization of satellite ray classification results.(a) Ray classification results of Global Positioning System (GPS) satellites in a GNSS Feature Map.The specific trajectory is shown together with used building model data in (b).

Figure 3 .
Figure 3. General conversion scheme for building model data referring to a map projection, satellite-position data referring to ECEF coordinates and user-position data referring to a map projection or ECEF coordinates.Blue arrows indicates the transformation steps into a common topocentric system, orange ones into the coordinate system of the map projection.

Figure 4 .
Figure 4. Simulation set up: Antenna positions and reflection surface, as well as distribution of signals affected by multipath.(a) Reflection surface and the four antenna positions indicated by a red cross, and (b) computed extra path delays in metric coordinates for the closest antenna position (2.2 m distance).

Figure 5 .
Figure 5. ∆δ = δ utm − δ enu vs. satellite azimuth, without any transformation to a common reference frame and no consideration of the meridian convergence to the δ utm computation.The metric values refer to the horizontal antenna reflector spacing and thus identify the antenna to which the corresponding data set belongs.

Figure 6 .
Figure 6.Differences in the computation of the extra path delay when applying different correction steps.(a) Section 6.2 transformation to a common reference system, (b) Section 6.3 different frames but consideration of the meridian convergence, (c) Section 6.4 transformation to common reference systems and consideration of the meridian convergence, (d) Section 6.5 transformation to common reference systems, consideration of the meridian convergence and distortion correction.

Table 1 .
Input data for 3DMA ray-tracing and their corresponding coordinate system and frame.