GNSS RUMS: GNSS Realistic Urban Multi-agent Simulator for Collaborative Positioning Research

: Accurate localization of road agents is the basis of intelligent transportation systems, which is still difficult to achieve for GNSS positioning in urban areas due to the signal interferences from buildings. Various collaborative positioning techniques are recently developed to improve the positioning performance by the aid from neighboring agents. However, it is still challenging to study their performances comprehensively. The GNSS measurement error behavior is complicated in urban areas and unable to be represented by naive models. On the other hand, real experiment requiring numbers of devices is hard to be conducted, especially for a large-scale test. Therefore, a GNSS realistic urban measurement simulator is developed to provide measurements for collaborative positioning studies. The proposed simulator employs a ray-tracing technique searching for all possible interferences in the urban area. Then, it categorizes them into direct, reflected, diffracted, and multipath signal to simulate the pseudorange, carrier-phase, 𝐶/𝑁 0 , and Doppler shift measurements correspondingly. The performance of the proposed simulator is validated through real experimental comparisons with different scenarios. The proposed simulator is also applied with different positioning algorithms, which verifies it is sophisticated enough for the collaborative positioning studies in the urban area.


Introduction
The development of intelligent transportation is one of the essential objectives of building a smart city, namely smart mobility [1].The intelligent transportation system (ITS) aims to alleviate traffic jams and enhance traffic efficiency with sufficient assistance and planning.However, as the foundation of ITS, the precise navigation of transportation participants can hardly be achieved, especially in the dense urban area [2].Most of the navigation algorithms rely on the global navigation satellite system (GNSS) for functionality or initialization, since it is the primary sensor directly providing the absolute positioning solution.Besides the systematic interferences that can be adequately modeled (e.g., atmospheric delay and ephemeris error), the buildings may obstruct or reflect the signals from the satellite.The reception of the reflection interfered signal or only the reflected signal is called the multipath or non-line-of-sight reception.Such interference leads to an enormous positioning error, possibly exceeding 30 meters, which is unqualified for ITS [2,3].Although various approaches, such as the integration of GNSS with INS, camera [4],or LiDAR-based simultaneous localization and mapping (SLAM) [5,6], are capable of achieving better positioning accuracy, they still require an accurate absolute positioning solution from GNSS for initialization.It is inevitable to improve the GNSS positioning accuracy in the urban area prior to ITS development.
Due to the rapid development of vehicle-to-everything (V2X) communication technology [7], the collaboration of transportation participants to achieve accurate positioning solution becomes possible, namely the collaborative positioning (CP).In most CP techniques, the neighboring road agents measure and share the inter-agent ranges as well as their positioning solutions, in order to optimize the overall positioning solution of each involving agents.By making use of the additional inter-agent information, the CP can achieve higher positioning accuracy compared with conventional standalone solution [8][9][10].Another straight-forward benefit of employing CP is to reduce the common noise by involving more road agents during positioning [11].Moreover, CP can make use of the network geometry formed by road agents to improve the overall positioning optimization [12].Unlike the transponder-based CP approach using DSRC or UWB, the GNSS-based CP can maintain the collaboration even when there exist obstacles between agents [13].Among various GNSS-based CP methods, a popular approach is to apply double difference (DD) on the pseudorange measurements from neighboring road agents to eliminate the systematic error during relative positioning, such as atmospheric delay and ephemeris error [14,15].Besides, based on the angular geometry between satellites and agents, the inter-agent range can be estimated for CP even when some of the satellite signals are blocked by buildings [16].However, these methods can be significantly degraded by the multipath and NLOS reception error.Recently, the 3D mapping aided (3DMA) GNSS technique has been developed, using the 3D building model to predict and mitigate those error in the urban area.Two representative approaches are the 3DMA GNSS shadow matching based on the satellite visibility prediction [17] and the 3DMA GNSS ray-tracing based on the LOS/NLOS pseudorange prediction [18].Then, the CP algorithms are extended by the 3DMA GNSS to achieve better performance in the urban area.The predicted satellite visibilities of neighboring road agents from the 3D building model are collaborated to achieve better position estimation [19,20].Moreover, a complementary integration of the 3DMA GNSS ray-tracing technique and the DD based GNSS CP is developed to eliminate the systematic error and mitigate the uncorrelated NLOS error simultaneously [21].The integration of 3DMA GNSS and CP has excellent potential to achieve satisfactory positioning accuracy in the urban area.
However, most of the CP algorithms are validated through naive simulated measurements, which behave very differently from the real experiments.A popular approach is to test the algorithm with additional simulated errors based on statistical models [15,22,23].Although the normal distribution can model many of the measurement noise, the GNSS measurement error in the urban environment is dominated by the enormous signal delay due to building reflection.This delay is uniquely related to the geometry between satellites, GNSS receiver, and the reflecting surface, which may differ between each road agent.Hence, it is inappropriate to validate the CP algorithm with the GNSS measurement error simulated by the normal distribution or other statistical distribution models.Some studies directly use real measurements collected from multiple receivers for validation [14,21], whereas the agent amount is limited to a few.It is tough to test large-scale performance with large numbers of agents or study the influence of network size.Therefore, it is necessary to develop a realistic GNSS measurement simulator to study the CP performance in the urban area comprehensively.
There exist various types of GNSS measurement simulator, from tracking-level to measurementlevel.The MUltipath Simulator Taking into Account Reflection and Diffraction (MUSTARD) developed by the Jet Propulsion Laboratory (JPL) is a famous GNSS measurement simulator [24], considering the signal reflection and diffraction effect based on the ray-tracing algorithm.The Satellite Navigation Radio Channel Signal Simulator (SNACKS) from the German Aerospace Center (DLR) further considers the multipath channel effect during simulation [25].Moreover, the commercial GNSS simulator SimGEN®+SE-NAV developed by the Spirent Communication not only use ray-tracing to track all valid reflection or diffraction signals, but also consider the interaction between those traced signals [26].To reduce the simulation complexity, another approach is developed, which simulates the tracking-level measurements in urban based on a reference data collected from software defined receiver (SDR) [27].The above simulators employ complicated algorithms to simulate the GNSS tracking-level measurements, whereas most CP studies only employs GNSS measurement-level data.For simplification, other GNSS simulators focusing on measurement-level data have been developed, such as the GNSS carrier-to-noise ratio ( / 0 ) simulation for positioning aiding [28], the GNSS LOS/NLOS pseudorange simulation by ray-tracing [18], or the Doppler shift modeling under multipath effect [29,30].However, many of them are not sophisticated enough with comprehensive error and noise modeling, for example, neglecting the multipath interaction error on pseudorange, or the tracking loop noise related to / 0 .Moreover, those GNSS simulators are designed for one specific kind of measurement of a single receiver.It is also necessary to develop a simulator supporting multi-agent full GNSS raw measurement simulation for CP studies.
In this study, a realistic GNSS measurement simulator is developed for multiple agents in the urban environment.Based on the satellite ephemeris, 3D building model, and the location of different road agents, the GNSS measurement-level data corresponding to each agent is simulated through the ray-tracing algorithm, including pseudorange, carrier-phase, / 0 , and Doppler frequency.To make the simulator suitable for CP studies, the locations of road agents are generated based on the dynamics and transportation mobility of road agents by SUMO (Simulation of Urban MObility) [31].By considering the direct, reflected, diffracted, and multipath GNSS signal through ray-tracing, the measurement-level GNSS data with sophisticated modeling noises are simulated for multiple road agents in the urban area, in order to supply realistic large scale data for CP research.The contributions of this study summarized as follows: 1) Providing multi-agent GNSS simulation for ITS and CP applications considering the transportation mobility from SUMO; 2) Holistically organizing and combing the important findings related to GNSS simulation from the related works and providing realistic measurement-level GNSS simulator based on comprehensive error models without referencing RF data from SDR; and 3) A Sophisticated GNSS simulator is developed for the urban scenario, in which the effects of signal reflection, diffraction and the interferences in-between on the pseudorange, / 0 , and Doppler frequency measurements are considered.
The rest of the paper is structured as follows: The overall structure of the proposed realistic multi-agent urban GNSS simulator is demonstrated in Section 2. The detailed procedures of raytracing simulation on direct, reflected, and diffracted GNSS signals are explained in Section 3. In Section 4, the GNSS measurement simulation will be elaborated based on different cases, including direct propagation, reflection, diffraction, and multipath.The modeling of GNSS measurement systematic error and noise are also explained in this section.Then, the performance of the proposed simulation is validated through different real experimental data and the positioning result based on popular positioning algorithms in Section 5.After discussing results and the future works in Section 6, the conclusion is drawn in Section 7.  The overall structure of the proposed GNSS realistic urban multi-agent simulator (RUMS) is shown in Figure 1.Table 1 shows the key models and their corresponding references being employed in the simulator developed in this paper.First, the locations or trajectories of all the involving road agents require to be generated before GNSS measurement simulation.For large scale CP studies, the transportation environment could influence the traffic flow and the available collaborator in a specific range.Therefore, it would be more realistic to evaluate the CP performance based on the simulated measurement considering the traffic condition.SUMO is an urban transportation simulator generating the trajectories of multiple road agents based on the transportation environment and facilities in a certain area.The traffic flow, road agent type (pedestrian or automobile), and other advanced settings can also be customized.In this study, we employ the simulated position   and velocity   of each road agent on GPS time from SUMO, namely the Road Agent PVT Data, to consider the urban mobility effect during the GNSS measurement simulation.Besides the information of road agents, the satellite positions   from ephemeris and the 3D building model with building corner positions   are also required before simulation.The last required information is the Receiver Parameter/Model representing the hardware characteristic of the GNSS receiver in the simulation, which will be interpreted with details in Section 4.5.

Simulator Structure
For each road agent, the ray-tracing simulation is firstly employed based on the agent and satellite position as well as the 3D building model, in order to search all valid GNSS signal ray via direct propagation, diffraction or reflection from each satellite.Based on the ray-tracing result, the measurement of each satellite will be categorized into four types following different simulation strategies: 1) LOS only: only the line-of-sight (LOS) signal propagation path is available and unobstructed; 2) Single diffraction: only one diffraction path is available; 3) Single reflection: only one reflection path is available; and 4) Multipath: more than one propagation paths are available, including any combination of the three cases above.
For the case of LOS only, the / 0 is simulated from an elevation-based open-sky / 0 regression model.The ranging measurement is the direct distance between satellite and agent without delay.The Doppler shift measurement is simulated from the Doppler effect model based on the location and velocity of the satellite and receiver.
For the single diffraction case, the / 0 measurement is simulated by the open-sky model and the uniform geometrical theory of diffraction (UTD) model [33].The corresponding ranging measurement is the summation of the direct range and the extra delay due to diffraction.The diffracted Doppler shift is modeled by the Doppler effect with an intermediate point [29].
For the single reflection case, the / 0 is simulated based on the open-sky model with an attenuation factor from the GNSS reflectometer (GNSS-R) model [37,47].The ranging measurement is modeled by adding the reflection delay on the direct distance.Similar to diffraction, the Doppler shift measurement of a reflected signal is simulated by the Doppler effect with an intermediate point.
For the case of multipath, the / 0 is simulated based on the superposition of fields from each valid signal, considering the interaction due to the differences in phase and amplitude.The corresponding ranging measurement is simulated by the multipath noise envelope based on the strength and delay of each signal, which takes the interferences in-between into account.The multipath Doppler shift is simulated based on the dominating signal path with the highest / 0 among all available paths.
Finally, the measurement systematic error and noise are simulated on top of the preceding C/N 0 , ranging and Doppler shift simulation result based on the GNSS receiver characteristics.The C/N 0 measurement noise is simulated based on the standard deviation of the narrow-to-wide power ratio evaluation method [40].The pseudorange measurement is simulated by incorporating the receiver clock bias, tropospheric delay, ionospheric delay, ephemeris error, satellite clock error, and a normally distributed noise based on the pseudorange error budget [39] for the above terms.The carrier phase measurement is simulated by adding the integer ambiguity on the pseudorange measurement.Moreover, the C/N 0 -related delay lock loop (DLL) and phase lock loop (PLL) noise are considered during the pseudorange and carrier phase simulation, respectively.The frequency lock loop (FLL) noise is also added on the preceding simulated Doppler shift.The simulated GNSS measurements of different agents are collected along GPS time and used for the CP algorithm test.

GNSS Ray-tracing Simulation
In an urban scenario, the GNSS measurement is mainly degraded due to three effects, the signal blockage by buildings, the diffraction effect when the direct signal path is closed to the building edge, and the reception of reflected signals from buildings.Therefore, three types of satellite signals can be received: the LOS signal, the diffracted signal, and the reflected signal, as shown in Figure 2, respectively.Based on the ray-tracing technique [48], all three types of signal can be traced with the geometrical relationship between buildings, satellite, and receiver.

LOS GNSS Signal
Assuming the GNSS signal propagation follows a straight line, the GNSS LOS signal of a specific satellite can be directly traced as a line from the satellite position to the receiver position, as Figure 2 (a) shows.Since there may exist buildings in urban that obstructing the GNSS signal, only the LOS signal without any intersection with the building surfaces from the 3D building model will be considered as the valid LOS signal.

Diffracted GNSS Signal
The GNSS signal is an electromagnetic wave broadcasting from the satellite to ground.The propagation of waves towards a specific receiver may be bent when its straight transmitting path is slightly obstructed by the building, namely the diffraction effect.More specifically, the GNSS diffraction can also be regarded as the reception of the secondary wave emitted from the building edge.A detailed explanation of the GNSS diffraction phenomenon can be found in [35].
Based on the UTD [34], the GNSS diffraction can be modeled as an effect dominated by a specific diffracted signal, as Figure 2 (b), where the signal propagates along a bent path via the point of diffraction on the building edge.A valid diffraction path follows the Fermat's principle that the angle  ′ between incidence and the building edge is equal to the angle  between the diffracted ray and the edge [32].By searching the point   on the edges of the 3D building model that fulfill such geometrical relationship between satellite and receiver, the ray-tracing technique can be extended to trace all valid path of the GNSS diffracted signal.In this study, we only consider the single diffraction effect and the diffracted signal being obstructed by other buildings is invalid.

Reflected GNSS Signal
Besides direct propagation and diffraction, the GNSS signal can also be reflected by the building surface and then received by the road agent, as Figure 2 (c) shows.The reflected signal follows the Fermat's principle that the incident angle  ′ is equal to the reflected angle  .Based on this geometrical relationship, the ray-tracing technique finds the reflected signal path by searching for valid point of reflection   on all the surfaces in the 3D building model.
In detail, the receiver (road agent) location is firstly mirrored with respect to a specific building surface and denoted as   .The point of reflection fulfills the geometrical relationship is located at the intersection between that surface and the line from the satellite to the mirrored location.If the point of reflection is within that mirroring surface and no obstruction from other building surfaces along the path, the corresponding path from the satellite passing the point of reflection to the receiver is a valid reflected signal.By applying this searching process to all the surfaces determined by the corner locations in the 3D building model, all valid reflected GNSS signal path can be traced.

GNSS Measurement Simulation
After searching all the valid direct, diffracted, and reflected GNSS signals from the ray-tracing, the measurement status of each satellite is categorized into four types: LOS only, single NLOS, single diffraction, and multipath.Different strategies are employed during the GNSS / 0 , pseudorange, carrier-phase, and Doppler frequency measurement simulation based on the measurement status.Finally, additional measurement error and noise based on comprehensive models are added on top of the preceding simulation result to make the final simulated measurement realistic.

Overall GNSS Measurement Simulation
To simulate the GNSS measurement realistically, all the error components need to be considered with appropriate models.For the / 0 measurement, besides the simulated value based on the signal strength, its hardware-related estimation noise also needs to be considered.The final simulated GNSS / 0 measurement is obtained by / 0, = / 0, + 10 log 10  / 0 (1) where / 0, is the simulated / 0 corresponding to the signal type and  ∈ {, , , ℎ}. / 0 ~(0,  / 0 2 ) is a random variable follows Gaussian distribution (denoted by ) with zero bias and the / 0 estimation variance  / 0 2 .The final simulated GNSS pseudorange measurement is given by where   is the simulated GNSS ranging measurement corresponding to different signal type,  is the speed of light,   is the satellite clock bias,   is the receiver clock bias.  and   are the ionospheric delay and the tropospheric delay, respectively.  ~(0,   2 ) and ) are random variable with the variance from the DLL tracking loop noise and the systematic error model noise, respectively.
The final carrier-phase measurement simulation result (in the unit of meters) is obtained from where the random variable   ~(0,   2 ) with the PLL tracking loop noise related variance,  is the wavelength corresponding to the simulated signal. is the integer ambiguity simulated by a random variable of integer that fixed for a specific satellite during operation, assuming no cycle-slip. ℎ is the multipath error on the carrier phase.
The final simulated Doppler shift measurement is obtained by where ∆  is the simulated GNSS Doppler shift corresponding to different signal types, and   is a random variable that   ~(0,   2 ) with the variance from the FLL tracking loop noise.Then, all four kinds of popular-used GNSS measurements are simulated and ready to be applied with different algorithms for evaluation.
The following Sections 4.2-4.5 demonstrate the procedures of obtaining the signal-type-related measurement error terms, corresponding to the cases of LOS-only, single diffraction, single reflection, and multipath.Section 4.6 introduces other systematic error terms related to the atmosphere, the hardware characteristics, etc.

LOS Only Case
The GNSS / 0 measurement relates to the strength of the signal, where a higher / 0 usually representing the signal being healthy with less interference.In an open-sky scenario without obstructions between a satellite and a receiver, the / 0 corresponding to a satellite signal is closely related to its elevation angle.Therefore, the / 0 of an unobstructed signal can be simulated based on a GNSS elevation-/ 0 regression model from long period open-sky data [28].In this study, the / 0 from GPS satellite, Beidou satellite on geosynchronous equatorial orbit (GEO), on inclined geosynchronous orbit (IGSO), and on medium Earth orbit (MEO) are simulated with different opensky models as Figure 3 shows, in order to consider the behavior varies between different systems and operating orbits.First-order polynomial fitting is employed between the elevation angle and the / 0 with the unit of Hz.The GNSS pseudorange measurement of the LOS only case is the direct range between satellite and road agent as follows.
The Doppler frequency measurement of the LOS only case can be simulated based on the Doppler effect between two dynamic objects [29], using where   and   are the velocity of the road agent and the satellite, respectively. is the lineof-sight unit vector from the road agent to the satellite. is the speed of light and  is the wavelength of the corresponding satellite signal. ̇ is the receiver clock drift as a tuning parameter of the simulation.

Single Diffraction Case
For the satellite with only one diffracted signal path valid from the ray-tracing, the corresponding / 0 , pseudorange, and Doppler measurements are simulated based on its geometrical behavior.Firstly, the / 0 is simulated based on the UTD, which describes the signal attenuation by the geometrical parameters of the diffracted signal and the building surfaces causing diffraction.The diffraction coefficient with respect to the strength of the LOS signal can be obtained from where   is the attenuation factor between the right-hand circular polarization (RHCP) electric field before and after diffraction on the point of diffraction. ∥ and  ⊥ are the soft and hard diffraction coefficient of a linear polarized field derived from the geometrical parameter of the diffracted signal and the building. ∥  and  ⊥  denote the orthogonal components of the GNSS signal electric field parallel and perpendicular to the incident plane, respectively. ∥  and  ⊥  denote the orthogonal components of the GNSS signal electric field parallel and perpendicular to the diffraction plane, respectively.Equation ( 4) extends the UTD from the linear polarized field to the RHCP field appropriate for the GNSS signal.  is the distance between the point of diffraction and the receiver,  is the GNSS signal wavenumber, and  is the imaginary unit. 1 √   ⁄ and  −  denotes the spreading factor and the phase shift of the signal emitted from the point of diffraction, respectively.The detailed procedures of computing the diffraction coefficient from the geometrical parameters can be found in [35].Then, the / 0 of the diffracted signal can be obtained using which is derived based on the relationships between the signal-to-noise ratio (SNR) with the diffraction coefficient and / 0 [49] as follows / 0 = 10  +  BW (11) where   and   denote the SNR peak power of the LOS and the diffracted signals, respectively.Equation (7) describes the relationship between / 0 (in the unit of dB-Hz) and SNR. BW is the Receiver front-end bandwidth which will be eliminated when only considers the power ratio between the LOS and diffracted signal.However, the / 0 of the original unobstructed LOS signal is hard to be obtained based on the receiver parameter that varies from different antennas.A convenient approach is to estimate the / 0, from an open-sky regression model, as the preceding LOS only case.After obtaining the / 0 of the diffracted signal, the corresponding GNSS ranging measurement can be simulated by its total traveling distance using where   ′ is the distance between the satellite and the point of diffraction.The pseudorange error due to diffraction can be denoted as   =   −   .Finally, the Doppler frequency of the diffracted signal can be derived based on the Doppler shift model extended with an intermediate point [29], using where   ′ is the line-of-sight unit vector from the point of diffraction   to the satellite and   is the unit vector from the road agent to the point of diffraction.Based on the position and velocity of the road agent and the satellite as well as the location of the point of diffraction, the Doppler frequency of the GNSS signal with a specific wavelength can be simulated correspondingly.

Single Reflection Case
The GNSS-R has been widely employed for geodetic survey application, which estimates the geometrical parameter during reflection based on the reflected / 0 measurement.The mechanism of GNSS-R can be, in turn, used to simulate the GNSS / 0 measurement based on the geometrical parameters during reflection.The / 0 relationship between the LOS signal and the reflected signal can be derived by based on Equation ( 7) and the relationship that where   ′ is the distance between the satellite and the point of reflection and   is the distance between the point of reflection and the receiver (road agent).Equation ( 10) is derived by Equation ( 11) and ( 12), which describe the reflected/LOS SNR peak power ratio [50] and the conversion from / 0 (in the unit of dB-Hz) to SNR, respectively. BW is the Receiver front-end bandwidth, which will be eliminated when only considers the power ratio between the LOS and diffracted signal.Therefore, the / 0 of a reflected signal can be represented by its LOS / 0 attenuated by a spreading factor and a reflection coefficient ℛ  .Noted that the GNSS signal is RHCP, whereas the reflected GNSS signal may be inverted into left-hand circular polarization (LHCP).By considering the polarization interference, the reflection coefficient ℛ  from RHCP to LHCP can be written as the combination of the reflection coefficient on horizontal and vertical polarization [36,51], as follows where ℛ ℎℎ and ℛ  are the Fresnel reflection coefficient of the horizontal and vertical linear polarization component with respect to the reflecting surface, derived from ′ is the incident angle of the reflected signal,  1 is the dielectric constant of the air, and  2 is the dielectric constant of the reflecting surface.In this study, the dielectric constant of glass ( 2 = 4.7) is used for the reflector, and the conduction due to the lossy medium is neglected.Noted that some articles denote ℛ ℎℎ by ℛ ⊥ and ℛ  by ℛ ∥ as the linear polarized component perpendicular and parallel to the incident plane (the plane determined by the incident and reflected signals), which are also representing the soft and hard reflections, respectively.Based on the geometrical parameter, including the incident angle  and the distances between satellite, road agent, and the point of reflection, the / 0 of the corresponding reflected GNSS signal can be simulated from its unobstructed / 0, by the open-sky model.The ranging measurement of single reflected GNSS signal can be simulated by the total travelling distance on the reflection path, as follows.
The corresponding pseudorange delay due to reflection can be written as The Doppler frequency of the reflected signal can be derived based on the Doppler shift model extended with an intermediate point [29], using where   ′ is the line-of-sight unit vector from the point of reflection   to the satellite, and   is the unit vector from the road agent to the point of reflection.Based on the position and velocity of the road agent and the satellite as well as the location of the point of reflection, the Doppler frequency of the GNSS signal with a specific wavelength can be simulated correspondingly.

Multipath Case
Besides the reception of a single GNSS signal from LOS, reflection or diffraction, the GNSS receiver can also receive multiple signals simultaneously, if available, namely the multipath case.The corresponding measurement will be degraded due to the interferences between each signal.In this study, we only consider the interaction between the shortest two signals with less than 20dB attenuation.Other available signals with larger delay or lower strength have fewer effects and they are neglected.Therefore, the multipath case contains the signal combination of: LOS and reflection, LOS and diffraction, reflection and diffraction, double reflections, and double diffractions.
Based on the expression of individual fields, the expression of the joint signal field with interferences between signals can be derived through the superposition of fields.The detailed derivation with different signal combination cases is demonstrated in Appendix A. Then, the strength ratio between the multipath joint field and the unobstructed field on the road agent location can be expressed by where   and   represent the stand-alone signal strength ratio compared to the unobstructed signal for different cases,   and   denote the extra traveling distance compared to the unobstructed signal for the diffracted signal or the reflected signal, respectively.Analogous to the single reflection case or single diffraction case with the approximation the extra delay due to building interference is negligible compared to the total range from satellite to receiver, the / 0 simulation of the multipath case can be obtained from / 0,ℎ = / 0, + 20| ℎ | On the other hand, the GNSS pseudorange measurement under the multipath case is also affected by the interaction between signals.The code pseudorange error due to this interaction needs to be evaluated by the multipath noise envelope model with the amplitude ratio and the delay between two signals.Based on the multipath noise envelope for the early-minus-late power discriminator [38], the multipath error on the pseudorange can be simulated by  ℎ =  ℎ   (, , ∆, ) where  ℎ   denote the multipath noise envelope model demonstrated in Figure 4,  denotes the signal amplitude ratio ( ≤ 1 for this case),  denotes the carrier phase offset and between two signals. is the time spacing between early and late correlator.∆ and ∆ are the multipath time and distance delay in-between, respectively. ℎ is the code chip width.In this study, the carrier phase shift induced by the reflector is neglected.Then, the GNSS code pseudorange measurement can be obtained by However, the multipath error from the interaction between signals on carrier phase measurement is much smaller compared to pseudorange measurement, since the baseband carrier has a shorter wavelength.The multipath error on the carrier phase can be simulated based on the PLL discriminator model, as follows [38].
For the Doppler effect, it is assumed to be dominated by the strongest signal [30], which has the highest / 0 .Hence, the Doppler shift for the multipath case can be simulated by the dominated signal through Equation ( 2), (9), or (19), as follows.

General Measurement Error and Noise
To ensure the GNSS measurement simulation realistic, other measurement errors and noises need to be considered besides the preceding building-related errors.In this study, sophisticated models are employed to simulate different kinds of error and noise generally appeared in GNSS measurements, including the tracking loop noise, / 0 estimation noise, satellite clock bias, receiver clock error, ionospheric delay, tropospheric delay, and pseudorange error modeling noise.

Tracking Loop Noise
The pseudorange measurement noise from DLL consists of the thermal noise and the dynamic stress [39].The thermal noise is related to the signal / 0 and the receiver parameters, whereas the dynamic stress is related to the loop order and bandwidth.This paper follows [39] to generate the tracking loop noise.Since the dynamic stress can be almost removed via the carrier aided code technique, it is neglected during the DLL noise simulation.As a result, the 1-sigma DLL noise of C/A code under BPSK-R modulation is obtained by where   is the 1-sigma thermal noise code tracing jitter,  is the early-to-late correlator spacing (unit of chip),   is the front-end bandwidth,   is the chip rate,   is the code loop noise bandwidth,   is the predetection integration time and   is the chip period.
The carrier phase measurement noise from a 3 rd order PLL consists of the thermal noise   , the vibration-induced oscillator jitter   , the Allan-variance-induced oscillator jitter  3 , and the dynamic stress  3 [39], as follows (in the unit of degrees).
where   is the carrier loop noise bandwidth,   is the L-band carrier frequency.  (  ) and (  ) are the oscillator vibration frequency sensitivity and the power curve, respectively, related to the random vibration modulation frequency   .  is the Allan deviation of the reference oscillator,  3   3 ⁄ is the maximum LOS jerk dynamics.The 1-sigma Doppler measurement noise   from a 2 nd order FLL is consist of the thermal noise frequency jitter   and the dynamic stress term   [39], which can be estimated as follows (in the unit of Hz).
where  equals to 1 for high / 0 and equals to 2 for the / 0 closed to the threshold 1/(12  ). 0 =   /0.53 is the loop filter natural radian frequency.Figure 5 shows an example of the 1-sigma noise distribution with respect to the / 0 value on measurements corresponding to DLL, PLL, and FLL tracking loop.Normally, a high / 0 indicates the signal is healthy, whereas a low / 0 indicates the signal has larger noise during the tracking loop.Therefore, the tracking loop noises on the GNSS pseudorange, carrier phase, and Doppler shift measurement are modeled by the random variables following zero-bias Gaussian distribution with the corresponding standard deviation.For real cases, the GNSS / 0 could also be noisy due to interferences or weak signal power.During the narrow-to-wide power method, a popular approach for / 0 estimation, the / 0 measurement noise can be evaluated statistically in the form of standard deviation [40].
The variable  and  denote that the prompt I and Q are divided into  intervals with  samples in each interval for the / 0 estimation. ℎ is the coherent integration time,  ̅ / and (  /  ) are the averaged value and standard deviation of the narrow and wide power measurements ratio on  intervals, respectively.However,  ̅ / and (  /  ) are obtained from the tracking-level measurements not available in the proposed simulator.An alternative approach is to model these two parameters by curve fitting models with large amounts of data in urban scenarios.In this study, based on the / 0 value (which is the / 0 with the unit of Hz), a rational fitting model is employed for  ̅ / , whereas an exponential fitting model is employed for (  /  ), as Figure 6 shows.Then,  ̅ / and (  /  ) can be modeled based on the simulated / 0 value, and further combined with the receiver parameter  ℎ ,  and  to simulate the noise on the / 0 measurement.

Satellite Clock Bias
The actual GNSS pseudorange measurement contains the satellite clock error due to imperfect synchronization with the corresponding system time.The correction for this error is estimated by the control segment and broadcast to users through the navigation message from satellite.Therefore, this error term can be obtained for the pseudorange and carrier phase simulation based on the correction parameters from the satellite ephemeris [39,52], as follows.

Receiver Clock Bias
Besides the satellite clock bias, the other term making the GNSS ranging measurement being 'pseudo-range', is the receiver clock bias, which is usually an unknown need to be estimated during positioning.In a real case, the receiver clock bias is not stable and drifts variously, related to the referencing oscillator of the receiver.A practical approach to simulate the receiver clock bias   with drift is through a 1 st order linear approximation model based on the receiver parameters [41,42], as follows.
=  ,0 + ∆  ( −  0 ) (45) ,0 is the initial clock bias of the receiver on time  0 ,  is the receiver current time.∆  is the frequency offset determined by the initial frequency offset ∆ 0 , random walk frequency modulation noise   , flicker frequency modulation noise   , white frequency modulation noise   , and the local oscillator frequency  0 .Figure 7 shows an example of the clock errors from   ,   , and   based on the TCXO parameters from [43].

Ionospheric Delay
The real GNSS pseudorange and carrier phase measurement contains the error due to ionospheric delay, which is usually corrected by the Klobuchar model with detailed procedures in [44], and summarized as below.
= 2( − 50400) ∑   Φ   3 =0 (48) where  is an elevation-related slant factor,   and   are the ionospheric parameters from ephemeris data, Φ   is the geomagnetic latitude,  is the local time.Based on the receiver latitude and longitude, satellite elevation and azimuth angles, receiver time, and the ionospheric parameters from the ephemeris, the corresponding ionospheric delay can be simulated via the Klobuchar model correspondingly.

Tropospheric Delay
Another error term always contained in the pseudorange and carrier phase measurements is the delay due to tropospheric refraction, described by the Saastamoinen model [45], as follows.
where  is the satellite elevation angle.  ,   ,   , and  are the total barometer pressure, the partial pressure of water vapor, absolute temperature (in Kelvin), and the coefficient term related to the ellipsoid height of the receiver.Based on the satellite elevation angle and the receiver ellipsoid height, the corresponding tropospheric delay   can be simulated correspondingly.

Pseudorange Error Modelling Noise
Besides the exact value of errors from different models, each error term may also have a certain variation during real operation.These variations of errors are caused by various unknown factors, which are hard to be modeled.A practical way to simulate these variations is according to the GNSS user-equivalent-range-error (UERE) budget from [39].The nominal 1-sigma error of the satellite clock   is given as 1.1 meters.The 1-sigma ephemeris prediction error  ℎ on pseudorange and carrier phase is given as 0.8 meters.The noise on the relativistic effect error is negligible.The 1-sigma error on ionospheric delay   is given by a higher bound as 7.0 meters, which is reduced in half in this study by assuming a good ionospheric correction.The 1-sigma error on Tropospheric delay is obtained from an elevation-angle-related model [46], as follows.
= 0.12 1.001 √0.002001 + sin 2  (50) As a result, all the error noise above can be combined into a total error noise factor, as the error modeling noise   shown below.
Then, the total modeling noise on the pseudorange measurements is simulated by a random variable from zero-bias Gaussian distribution with the standard deviation of   .

Simulation Performance Verification
In this section, the proposed realistic GNSS measurement simulator for multiple agents in the urban area is verified by two approaches.On the one hand, the simulated GNSS measurements are compared with the real experimental measurements collected in different scenarios to verify whether different interferences can be appropriately simulated.On the other hand, the simulated GNSS measurements of different agents will be applied with conventional positioning algorithms and collaborative positioning algorithms, aiming to verify that simulated measurements can reproduce realistic performances of current positioning algorithms.

Experiment Setup
The experimental verification of the proposed simulator consists of three scenarios: 1) One-hour static experiment on the intersection of the urban area; 2) dynamic vehicular experiment in the urban area; and 3) static experiment with multi-agent in different environments.For the verification using the one-hour static data, the experiment locations and the corresponding sky-plots with buildings are shown in Figure 8.The GNSS measurement on this location is expected to under severe degradation from reflection, diffraction, or multipath.For the dynamic experiment, the trajectory is shown in Figure 9, where the vehicle moving from the open-sky area to the urban area, and then back to the starting point.The receiver locations and sky-plots during the multi-agent static experiment are shown in Figure 10 with a summary of the surrounding environment of each receiver in Table 2.During all the experiments, the commercial-grade GNSS receiver ublox EVK-M8T with patch antenna is employed to collect GNSS (GPS and Beidou) raw measurements.The receiver true location of each static experiment is obtained based on the landmark on Google Earth.The receiver true location during the dynamic experiment is obtained based on RTK GNSS/INS integrated solution from the Novatel SPAN-CPT.The value of receiver parameters used in the proposed GNSS simulator are listed in Table 3.During the simulation, unlike the / 0 and Doppler frequency, the pseudorange measurement always contains the receiver clock bias, which is a receiver-related unknown.It is hard and unnecessary to simulate this bias exactly the same as the real measurement, since it is also estimated in the solution during positioning.Therefore, instead of comparing with the bias-contained pseudorange, we directly compare the pseudorange error between real measurement and simulation for verification.The pseudorange error in real measurement is labeled by applying the double different technique with the measurement from a nearby reference station [53].Similarly, the real Doppler frequency measurement also contains a receiver clock drift hard to be exactly simulated.Therefore, we apply the single difference approach to estimate the Doppler shift error for the consistency evaluation between real measurement and simulation.The Doppler shift error of the  ℎ satellite is estimated by where the superscript  denotes the Doppler shift of the selected master satellite (with the highest elevation angle).The subscript  denotes the Doppler shift from measurement or simulation.The sub-script  denotes the true Doppler shift labeled using the user velocity, satellite velocity, and the satellite line-of-sight vector.By applying the single difference on the Doppler shift between the target satellite and the master satellite, the shared receiver clock drift can be eliminated, which provides a clearer consistency evaluation between the real received and the simulated Doppler shift.During the static experiment on the urban intersection, the GNSS measurements from each satellite are simulated by the proposed simulator and compared with the real received measurement.Figure 11 shows the simulation result of / 0 , pseudorange error, and Doppler shift error from the B4 satellite, which is classified as a LOS satellite during the whole simulation.The simulated / 0 measurement is consistent with the real measurement, verifying the proposed open-sky regression model can be used to simulate the / 0 from the LOS satellites.For the pseudorange measurement, the simulation result is also consistent with the real measurement, which is a typical LOS measurement without large errors.Although the real measurement has lees fluctuation compared to the simulation due to some filtering technique in receivers, both simulated and real measurements have a similar magnitude of noise, representing the behavior of the LOS measurement.For the Doppler shift, the error from simulation has a slightly higher noise compared to that from the real measurement, but both of them are in a similar level.The simulation results from the G20 satellite experiencing all four measurement status are compared with the real measurements in Figure 12.In the first 350 seconds, only one valid reflected signal path is found by the ray-tracing, which corresponding to the single reflection case in the simulator.The / 0 is simulated with a large attenuation to around 28 dB-Hz, whereas the pseudorange measurement is simulated with around 30 meters delay.In this period, the real received measurement is consistent with the simulation, containing a similar / 0 attenuation and pseudorange delay.In the later period, the proposed simulator only finds one valid diffracted signal while the previous reflection becomes invalid, which turns the measurement status to single diffraction case.The / 0 is gradually increased, since the satellite elevation angle rises towards the building edge with a decreasing signal strength attenuation coefficient.On the other hand, the pseudorange errors in real measurement and simulation both are reduced since the diffraction delay is normally negligible comparing to other noise.After the 500 epoch, the satellite is raised from the building boundary and becomes visible to the receiver, making both diffracted and direct signal path valid to the receiver in the simulation.Therefore, the measurement is simulated with the multipath case.Consistent with the real measurement, the / 0 is continually increasing while the simulated pseudorange error fluctuates around zero.After the 1300 epoch, the diffracted signal path becomes invalid for the receiver, and only the direct signal remains.Similar to the real measurement, the / 0 is maintained in a high value and the simulated pseudorange error is dominated by zero-biased random noises.In the ending period, the real / 0 is attenuated compared to the simulation, but no large constant delay is found in the pseudorange error.This could possibly be due to the signal interferences from unfounded reflections.For the satellite G29 with large pseudorange errors probably due to reflection, the corresponding simulated and real received measurements are shown in Figure 13.During the experiment, its real received / 0 measurement is fluctuated with a large amplitude.This is a typical consequence of the multipath effect, where the phase shift between two signal paths is rapidly switching between inphase and out-of-phase due to the satellite movement.Moreover, the pseudorange error in the real received measurement also has a large fluctuation but contains an enormous bias from zero.This is also another evidence of occurring multipath effect, or more specifically, the multipath effect with a dominating reflected signal and no valid direct signal.The proposed simulator traces multiple valid propagation paths and correctly simulates the occurrence of multipath effect in the beginning (except the first 50 epochs) and ending period.Both the overall biases and the fluctuations in the pseudorange error are simulated consistent with the real measurement in these periods.However, during the period around 850 to 3200 epoch, the simulator cannot find another valid signal other than the dominating reflecting signal from ray-tracing.Therefore, only the enormous delay in the pseudorange is simulated, but the large fluctuation is missed.Despite this, the proposed simulator can still mimic the dominating pseudorange error due to reflection, which has the most impact on the positioning performance in the urban area.As the measurements having the greatest impact on the positioning performance, the GNSS / 0 and pseudorange simulation results for other satellites are shown in Figures 14 and 15, respectively.For most of the satellites, the simulated / 0 measurements are consistent with the real measurements, especially many of them appropriately simulate the attenuation on signal strength due to reflection or diffraction.On the other hand, the simulated pseudorange error is consistent with the real measurement in most epochs, even for those satellites with enormous delay due to reflection.However, for some of the satellites, the simulated / 0 and pseudorange are stable and healthy, while the real measurement is biased or noisy with fluctuations.This is probably due to the limitation of the ray-tracing technique, only considering the level-of-detail one (LoD-1) building model.The ray-tracing technique cannot trace those valid reflection or diffraction paths relate to the detailed structure of the building.As a result, the corresponding interferences are missed during the simulation.Despite this, those missed interferences are relatively small compared to those interferences introduced by severe reflections or diffractions, which are properly simulated.In general, the proposed simulator is able to simulate the measurement behavior of each satellite in the urban area in a realistic manner.The least squares positioning result using the simulated measurements are shown in Figure 16, compared with the positioning result using real measurements.Although the positioning result based on the real measurement has a larger variance than the simulation-based solutions, both the positioning solution has a similar error distribution mainly along the road direction.The corresponding positioning errors during the experiment are shown in Figure 17 with respect to East-West (E-W) and North-South (N-S) direction.The positioning errors from simulation have the same trend with the positioning errors from real measurements, except few epochs with large errors that unable to be traced, probably due to the limitation of ray-tracing and model accuracy.In summary, the proposed simulator can realistically simulate the GNSS measurement and corresponding positioning error behavior in the urban area during the long-period static experiment.

Dynamic Experimental Verification
During the dynamic experiment, the GNSS measurements from each satellite available in the ephemeris are simulated by the proposed simulator.The simulation result of the B1 satellite is shown in Figure 18, compared to the real received measurements.Most of the time, the satellite elevation angle is higher than the building boundary.The proposed simulator appropriately categorizes it as LOS or multipath dominated by the direct signal, and only a small pseudorange error is simulated, which matches the real measurement in the experiment.Besides, during the epoch around 430 and 700, the proposed simulation can appropriately simulate the attenuation of / 0 due to multipath (interference between two diffracted signals) and reflection.The slight increment of pseudorange error during 660-780 epoch due to reflection from the simulation is also consistent with the real measurements.Notice that, the real received / 0 has additional attenuation compared to the simulation in some epochs.This could possibly be due to the interferences from the signal transporting along an unexpected path that miss-detect by the ray-tracing algorithm.For example, during the first 170 epochs, the vehicle is nearby a tall building still under construction without a proper model.The interferences related to that building cannot be simulated correctly.In general, the simulated GNSS measurement show good consistency with real measurements.On the other hand, the simulation results of the satellite B6 is demonstrated in Figure 19.The GNSS measurement status is frequently changing between four cases during the simulation since the surrounding environment is rapidly changing during the experiment.Similar to the preceding result, the simulated measurement of / 0 and pseudorange error from this satellite are consistent with the real measurement, including the signal attenuation/delay due to the multipath effect.For the Doppler shift, the proposed simulator can simulate the error with a level and trend consistent with the real measurement under a rapidly changing environment.Noticed that there may still have few epochs that the real measurement has greater degradation than simulation, due to the receptions of unexpected interfered signals.From a simulation point of view, the proposed simulator can instantly determine the measurement status based on the surrounding environment.The measurement behavior for each status is comprehensively simulated with reasonable magnitude and fluctuation.The rapid change of measurement behavior also reflects the unstable measurement status and the surrounding environment during the dynamic operation.It is much more realistic than the statistical noise model that is normally used to validate the GNSS positioning algorithms in the urban area.The overall positioning solutions and corresponding errors from simulation and real measurements are shown in Figure 20.The proposed simulator appropriately simulates the increment of positioning error when entering a dense urban area, and the recovery of positioning accuracy when back to an open-sky environment.The positioning performance based on simulation behaves similar to that based on real measurements, except the beginning and ending periods.As mentioned before, the vehicle is nearby a constructing building without a model in the beginning period.Hence, the measurement errors from that building are underestimated, leading to a much smaller positioning error from the simulation.The second inconsistent period is probably because the vehicle is nearby the building with a complicated (ship-like) structure.The ray-tracing algorithm based on the LoD-1 building model cannot comprehensively predict all the valid interfered signals related to this building.As a result, the corresponding GNSS measurement error and the positioning error are underestimated.In summary, the proposed simulator can realistically simulate the GNSS measurement and positioning behavior of a vehicular agent in the urban area.During the multi-agent experimental verification, six receivers are set up at the designed locations (Figure 10) and collect the GNSS simultaneously.On the other hand, the proposed simulator simulates the GNSS measurements of those receivers during the same period based on the receiver locations, satellite ephemeris, and 3D building model.After that, the GNSS positioning behavior from the simulator is compared with that from real collected measurements to verify the performance of the proposed simulator.
The GNSS least squares positioning results of each agent based on real measurements and simulator are shown in Figure 21.The corresponding positioning 2D errors are shown in Figure 22.For Agent 1 and Agent 2 located in the open-sky scenario, the proposed simulator classifies most of the measurements are healthy.Only small errors are simulated in the measurements, which provides accurate positioning solutions consistent with the real measurement performance.Note that the positioning result from simulation contains higher noise than the real measurements.This could possibly be because of the measurement filtering technique that is usually embedded in the commercial GNSS receiver.For other receivers, the positioning error distribution from the proposed simulation is very similar to that from the real measurement.Especially for Agent 4, the increment of positioning error is also simulated at the epoch consistent with the real measurement.In summary, the proposed simulator can simultaneously simulate the GNSS measurements of multiple agents with realistic error behaviors in the urban area.

Simulation Setup
Besides validating the consistency between the proposed simulator and the real measurements, it is also important to evaluate the positioning performance after applying different algorithms.Many algorithms are developed based on a certain assumption; for example, the measurement error follows the normal distribution.If only the naive errors are simulated in the measurements, these errors can be easily mitigated by applying advanced positioning algorithms, resulting in an unrealistic outstanding performance.A realistic simulator needs to generate the measurements that maintain the difficulties in urban GNSS positioning, which can be employed to appropriately evaluate or improve the performance of various positioning algorithms.
To evaluate the positioning performance based on the measurements from the proposed simulator, a realistic multi-agent positioning scenario in the urban area is generated from SUMO [31] beforehand.The trajectory of each agent from SUMO considering urban transportation behaviors is shown in Figure 23.Five road agents are simulated, including one pedestrian in the open-sky area (Agent 1), one pedestrian in the dense urban area (Agent 2), and three vehicular agents operating with changing environments (Agent 3-5).After simulating the GNSS measurements corresponding to each agent trajectory by the proposed simulator, the simulated measurements are applied with six different urban GNSS positioning algorithms as follows: (1) LS: Conventional least squares positioning for a single agent; (2) CC: Least squares positioning with consistency check for a single agent [54]; (3) RT: 3DMA GNSS ray-tracing positioning for single agent [18]; (4) RT-CP: 3DMA GNSS ray-tracing based collaborative positioning with factor graph optimization [21].
Figure 23.The trajectory of each agent simulated from SUMO [31] for the proposed multi-agent GNSS measurement simulation and positioning performance evaluation.Agents 1 and 2 are simulated as pedestrians while Agents 3-5 are simulated as vehicles.The arrow on the trajectory denotes the moving direction.

Positioning Performance
The solutions of each agent by applying the above four different GNSS positioning algorithms to the simulated measurements are demonstrated in Figure 24.The corresponding 2D positioning error distributions and root-mean-square errors (RMSE) of each agent are shown in Figure 25 and Table 4, respectively.For Agent 1 in the open-sky area, all four positioning algorithms have similar accurate performances, since most of the measurements are simulated without interferences from buildings.
For Agent 2 in the urban area, the least-squares positioning RMSE is close to 30 meters, which is reasonable to a low-cost GNSS receiver sensitive with the reflection interferences from buildings.In this scenario, the agent moves along one side of the buildings, and most of the reflected signals are probably from the opposite buildings 70 meters away.Moreover, the sky-view of Agent 2 is very limited due to building blockage, which makes the multipath effect containing the direct signal less likely to occur.Hence, besides the healthy LOS measurements, the rest are more likely to be the reflected measurements with enormous delay.In this case, those degraded measurements are very inconsistent with other measurements.Therefore, the consistency check method can easily detect and isolate those outliers, achieving excellent positioning performance.On the other hand, the occurrence of multipath effect with two reflected signals is possible for this scenario, which reaches the limitation of the current ray-tracing positioning technique only considers one dominating signal.As a result, not much improvement is achieved by applying the ray-tracing positioning algorithm.Similarly, the ray-tracing based collaborative positioning only slightly improves the positioning accuracy by the aid from neighboring agents.
For vehicular Agent 3 and Agent 4, the surrounding environment is constantly changing, which has a higher chance of experiencing the scenario with complicated measurement error behavior, for example, the intersection with complicated multipath effects or the narrow street with very limited measurement number.Therefore, the consistency check only achieves limited improvements.On the other hand, the ray-tracing algorithm corrects the reflected measurement, which guarantees a sufficient amount of reliable measurement for better positioning accuracy.The ray-tracing based collaborative positioning method further employs the measurements from neighboring agents to eliminate the systematic errors, resulting in a slightly better performance.
For Agent 5 that is always in a narrow street, many of the measurements contain enormous errors, which makes both the LS method and CC method significantly degraded.The ray-tracing method is only capable of correcting the single reflection delay, which limits the positioning improvement.Its extension with collaborative positioning (RT-CP) can only reduce the RMSE to 18.6 meters.
In summary, the performances of different positioning algorithms applied to our simulated measurements are reasonable and consistent with the real performance of the low-cost GNSS receiver in the urban area.Therefore, the proposed simulator is capable of providing realistic GNSS measurements of multiple agents in the urban scenario for the study and analysis of collaborative positioning algorithms.

Discussion
The results from the proposed simulator show great consistency with the real GNSS measurement in the urban area, which validates its capability to provide realistic GNSS measurements for various urban positioning algorithm evaluations or developments.The proposed simulator is developed based on sophisticated models, covering most of the interferences in the urban area.As the results in Table 4, by applying different advanced positioning algorithms, the remaining positioning errors are consistent with the real experimental performances reported in [18,21,54].Therefore, the proposed simulator can appropriately reflect the challenges of urban GNSS positioning for future studies, especially the potential algorithms hard to conduct experimental verification, such as the large-scale collaborative positioning algorithms.
Compared to other existing algorithms, the proposed simulator is much less complicated, but still maintains a sufficient verisimilitude level.Besides the measurement availability prediction achieved by the simulator from UCL [55], our simulation further simulates all the basic GNSS raw measurements employed for positioning.Both the conventional GNSS raw measurement simulator, MUSTARD [24] and SNACS [25], simulates the interferences most likely from the ground, whereas our simulator focus on a more complex urban scenario considering severe interferences come from buildings.Another advanced GNSS simulator, SimGEN @ +SE-NAV [26], also employs the ray-tracing algorithm and sophisticated models for the interferences simulation in the urban area.However, it involves complicated radio channel modeling, which may not be useful for verifying the algorithms only using measurement-level data.Moreover, the detailed procedures of interference modeling in [26] are not explained.A recent study employs a reference RF data to simplify the channel modeling during GNSS measurement simulation, but still requires SDR for simulation.Different from the above approaches, our proposed simulator avoids the complicated channel modeling by using an open-sky / 0 model for reference, and directly simulates measurement-level GNSS data, which is convenient and adequate to verify various positioning algorithms.The modeling of different interferences, including reflection, diffraction, and multipath, is explained comprehensively, which can be extracted for individual study or evaluation by other scholars.
Besides providing realistic GNSS data, the proposed simulator can, in turn, improve the existing 3DMA GNSS positioning algorithms.Most of the 3DMA GNSS positioning algorithms determine the agent position by searching for a candidate location with the predicted measurement best matches the real measurement.However, the current methods only conduct matching on the satellite visibility or the direct/reflected pseudorange [17,18].The diffraction and multipath models in this study can be used to extend the 3DMA GNSS positioning algorithm by considering the matching on the / 0 involving diffraction or multipath.The 3DMA GNSS considering the diffraction will be our future work.
However, the proposed simulator still has two limitations.On the one hand, the measurement simulation is highly relying on the ray-tracing technique, which still unable to consider all the interferences from buildings comprehensively.An advanced ray-tracing technique needs to be developed with the consideration of detailed building models and material effects in the future.On the other hand, the proposed simulator only provides the measurements based on individual epoch, which neglects the GNSS receiver dynamics during the operation.Therefore, different filtering techniques need to be considered to model the measurement dynamics along time in the future.

Conclusions
In this study, a realistic GNSS measurement simulator for multiple agents in the urban area is developed, considering the measurement degradation due to surrounding buildings.Four types of GNSS measurements, including LOS, reflection, diffraction, and multipath, are considered during simulation, covering most of the interferences in the urban area.From each measurement point of view, both the one-hour static and vehicular dynamic experiments validate that the measurement from the proposed simulator has a consistent and reasonable error behavior compared to the real collected measurement.From the overall performance point of view, the multi-agent experiment validates that the positioning results based on the measurement from the proposed simulator can realistically reflect the GNSS positioning error behavior in the urban area.Moreover, the proposed simulator is integrated with SUMO to simulate the GNSS measurement of multiple agents considering the transportation behaviors.The simulated measurements are further applied with advanced positioning algorithms, which verifies that the proposed simulator can appropriately express the current difficulties of precise positioning as well as the bottleneck of different positioning algorithms.Therefore, the proposed simulator can provide realistic GNSS measurements of multiple agents to study and improve the state-of-art GNSS collaborative positioning algorithms in the urban area.
Supplemental files: This paper provides 3 examples of the GNSS observation data in the format of the RINEX 3.02.The ground truth of the data is also included.Please be noted the GNSS ephemeris data (RINEX nav file) is not included since the data can be easily downloaded online.Please follow the Readme.txtfor the detail information on downloading the ephemeris data via HK SatRef (https://www.geodetic.gov.hk/en/satref/satref.htm).

Figure 1 .
Figure 1.The flowchart of the proposed GNSS realistic urban multi-agent simulator, GNSS RUMS.

Figure 2 .
Figure 2. The valid GNSS signal tracked by ray-tracing algorithm, including: (a) the LOS signal; (b) the diffracted signal; and (c) the reflected signal.

Figure 4 .
Figure 4.The demonstration of the multipath noise envelope when the delayed signal has half amplitude attenuation compared to the early signal.

Figure 6 .
Figure 6.The curve fitting model for C/N 0 noise simulation based on large amount data in the urban scenario: (a)  ̅ / fitted by a rational model; (b) (  /  ) fitted by an exponential model.

Figure 7 .
Figure 7.An example of the simulated values of random walk frequency modulation noise   , flicker frequency modulation noise   and white frequency modulation noise   based on the receiver parameters with TCXO from [43].

Figure 8 .
Figure 8.The location and the corresponding sky-plot with building boundaries for the one-hour static experiment on an intersection in an urban area.The right sky-plot demonstrate the satellite distribution and the building blockage (shaded area) on the experiment location (red marker).

Figure 9 .
Figure 9.The trajectory during the dynamic vehicular experiment.

Figure 10 .
Figure 10.The receiver locations and the sky-plots during the multi-agent static experiment.

Figure 11 .
Figure 11.The comparison of the simulated and actual received GNSS measurements from B4 satellite during the one-hour static experiment, including: (a) / 0 ; (b) pseudorange error; and (c) Singledifferenced Doppler shift error.The black marker and red marker denote the real and simulated measurement, respectively.The green background denotes the measurement type as LOS during the simulation.

Figure 12 .
Figure 12.The comparison of the simulated and actual received GNSS measurements from G20 satellite during the one-hour static experiment, including: (a) / 0 ; (b) pseudorange error; and (c) Doppler shift error.The coloured background denotes the measurement status being categorized during the simulation, including: LOS, Single reflection, Single diffraction and Multipath.

Figure 13 .
Figure 13.The comparison of the simulated and actual received GNSS measurements from G29 satellite during the long-period static experiment, including: (a) / 0 ; (b) pseudorange error; and (c) Doppler shift error.The coloured background denotes the measurement status being categorized during the simulation, including: LOS, Single reflection, Single diffraction and Multipath.

Figure 14 .
Figure 14.The simulated / 0 and corresponding real measurements from each satellite during the long-period static experiment.The notation below each sub-figure denotes the index of each satellite from different constellation, where 'G' for GPS and 'B' for Beidou.

Figure 15 .
Figure 15.The simulated pseudorange error and corresponding labelled value in real measurements for each satellite during the long-period static experiment.

Figure 16 .
Figure 16.The least squares positioning solution based on (a) real received GNSS measurement; (b) simulated GNSS measurement from the proposed simulator.

Figure 17 .
Figure 17.The least squares positioning error from the real measurements and simulated measurements during the one-hour static experiment.

Figure 18 .
Figure 18.The comparison of the simulated and actual received GNSS measurements from B1 satellite during the dynamic vehicular experiment, including: (a) / 0 ; (b) pseudorange error; (c) Doppler shift error; and (d) the elevation angles of the satellite and the building boundary on satellite's azimuth direction.The coloured background denotes the measurement status being categorized during the simulation, including: LOS, Single reflection, Single diffraction and Multipath.

Figure 19 .
Figure 19.The comparison of the simulated and actual received GNSS measurements from B6 satellite during the dynamic vehicular experiment, including: (a) / 0 ; (b) pseudorange error; (c) Doppler shift error; and (d) the elevation angles of the satellite and the building boundary on satellite's azimuth direction.

Figure 20 .
Figure 20.(a) The positioning solutions from real measurements (blue) and simulated measurements (red) during the dynamic experiment; (b) the corresponding least squares positioning 2D errors from the simulation and real measurements.

Figure 21 .
Figure 21.The positioning solutions of different agents based on real measurements (blue) and simulated measurements (red) during the multi-agent experiment.The number denotes the agent index corresponding to Figure 10.

Figure 22 .
Figure 22.The least squares positioning 2D errors of different agents based on the simulated and real measurements during the multi-agent experiment.The number is corresponding to the agent index in Figure 21.

Figure 24 .
Figure24.The solutions of different positioning algorithms applied on the GNSS measurements from the proposed simulator on different scenarios, including the least squares positioning method (LS), the least squares method with consistency check (CC), the 3DMA GNSS ray-tracing positioning method (RT), and the 3DMA GNSS ray-tracing based collaborative positioning with factor graph optimization (RT-CP).

Figure 25 .
Figure 25.The 2D positioning error of different algorithms corresponding to each agent during the multi-agent simulation performance verification.

Table 1 .
The employed models and the corresponding references in the proposed simulator.

Table 2 .
Surrounding environment of each receiver during the multi-agent experiment.

Table 3 .
The value of receiver-related parameters during simulation.

Table 4 .
The root-mean-square error (in meter) of different GNSS positioning algorithms applied on the simulated measurements.