Simulations of direct and reflected waves trajectories for in situ GNSS-R experiments

Simulations of direct and reflected waves trajectories for in situ GNSS-R experiments N. Roussel, F. Frappart, G. Ramillien, C. Desjardins, P. Gegout, F. Pérosanz, and R. Biancale Université de Toulouse, CNRS, IRD, GET-OMP, Toulouse, France Groupe de Recherche en Géodésie Spatiale, Toulouse, France Collecte Localisation Satellites, Ramonville Saint Agne, France Centre National d’Etudes Spatiales, Toulouse, France


Abstract
The detection of Global Navigation Satellite System (GNSS) signals that are reflected off the surface, together with the reception of direct GNSS signals offers a unique opportunity to monitor water level variations over land and ocean.The time delay between the reception of the direct and the reflected signal gives access to the altitude of the receiver over the reflecting surface.The field of view of the receiver is highly dependent on both the orbits of the GNSS satellites and the configuration of the study site geometries.A simulator has been developed to determine the accurate location of the reflection points on the surface by modelling the trajectories of GNSS electromagnetic waves that are reflected on the surface of the Earth.Only the geometric problem have been considered using a specular reflection assumption.The orbit of the GNSS constellations satellite (mainly GPS, GLONASS and Galileo), and the position of a fixed receiver are used as input.Three different simulation modes are proposed depending on the choice of the Earth surface (local sphere or ellipsoid) and the consideration of topography likely to cause masking effects.Atmospheric delay effects derived from adaptive mapping functions are also taken into account.This simulator was developed to determine where the GNSS-R receivers should be located to monitor efficiently a given study area.In this study, two test sites were considered.The first one at the top of the Cordouan lighthouse (45 • 35 11 N; 1 • 10 24 W; 65 m) and the second one in the shore of the Geneva lake (46 • 24 30 N; 6 • 43 6 E, with a 50 m receiver height).This site is hidden by mountains in the South (altitude up to 2000 m), and overlooking the lake in the North (altitude of 370 m).For this second test site configuration, reflections occur until 560 m from the receiver.The geometric differences between the positions of the specular reflection points obtained considering the Earth as a sphere or as an ellipsoid were found to be on average 44 cm for satellites elevation angle greater than 10 • and 1 m for satellite elevation angle between 5 • and 10 • .The simulations highlight the importance of the DEM integration: differences with and without integrating the DEM were found to be about 3.80 m with the minimum elevation angle equal to 5 • and

Introduction
The Global Navigation Satellite System (GNSS), which includes the American GPS, the Russian GLONASS, and the European Galileo (which is getting more and more denser) uses L-band microwave signals to provide accurate 3-D positioning on any point of the Earth surface or close to it.Along with the space segment development, the processing techniques have also improved considerably, with a better consideration of the various sources of error in the processing.Among them, multipaths still remain a major problem, and the mitigation of their influence has been widely investigated (Bilich, 2004).ESA (European Space Agency) first proposed the idea of taking advantage of the multipaths phenomenon in order to assess different parameters of the reflecting surface (Martin-Neira, 1993).This opportunistic remote sensing technique, known as GNSS-Reflectometry (GNSS-R), is based on the analysis of the electromagnetic signals emitted continuously by the GNSS satellites and detected by a receiver after reflection on the Earth's surface.Several parameters of the Earth surface can be retrieved either using the time-delay between the signals received by the upper (direct signal) and the lower (reflected signal) antennas, or by analyzing the waveforms (temporal evolution of the signal power) corresponding to the reflected signal.This Figures

Back Close
Full technique offers a wide-range of applications in Earth sciences.The time-delay can be interpreted in terms of altimetry as the difference of height between the receiver and the surface.Temporal variations of sea (Lowe et al., 2002;Ruffini et al., 2004;Löfgren et al., 2011;Semmling et al., 2011;Rius et al., 2012) and lakes level (Treuhaft et al., 2004;Helm, 2008) were recorded with an accuracy of a few centimeters using in situ and on-board aircraft antennas.Surface roughness can be estimated from the analysis of the Delay-Doppler Maps (DDM) derived from the waveforms of the reflected signals.They can be related to parameters such as soil moisture (Katzberg et al., 2006;Rodriguez-Alvarez et al., 2009, 2011) over land, or wave heights and wind speed (Komjathy et al., 2000;Zavorotny and Voronovich, 2000;Rius et al., 2002;Soulat et al., 2004) over the ocean, or ice properties (Gleason, 2006;Cardellach et al., 2012).GNSS-R technique presents two main advantages: (1) a dense spatial and temporal coverage, not only limited to a single measurement point or a non repetitive transect as using classical GNSS buoys, (2) a guarantee of service for the next decades (because of the strategic role played by these systems).GNSS-R altimetric accuracy is today at the level of few centimeters.But this technique will benefit, in the future, from improved processing technique and from the densification of the GNSS constellation.
The commonly-used GNSS-R system consist of two antennas (Fig. 1): the first one is right-hand circular polarized (RHCP) and zenith-facing to receive the direct waves.
The second one, left-hand circular polarized (LHCP) and nadir-facing to receive the reflected waves.These reflected waves will change their polarization from RCHP to LHCP by reflecting.The reflected signals have an additional path delay with respect to the direct ones.The analysis of the path difference between these direct and reflected signals is used to estimate the relative height difference between the two antennas.In order to anticipate the impact of the geometric configuration of the experiment, a simulator has been developed to estimate the positions of reflection points using a specular reflection point assumption.Three different methods were implemented: approximating the Earth as a sphere, as an ellipsoid or integrating a Digital Elevation Model (DEM).In addition, the signal bending due to the neutral part of atmosphere is taken into account Figures

Back Close
Full using the Adaptive Mapping Functions (AMF) from Gegout et al. (2011).Simulations were performed for different configurations: variations in the reflectometer height, mask effects due to terrain, satellite network geometry.This article is composed by three main parts following the logical structure of the Fig. 2. The first part presents the datasets used for initiating simulations, the second one concerns the methodologies for the determination of the reflection points while the last one deals with the simulator performances and simulation results.

Design of the simulator
The simulator has been developed in the GNU R language, generally used for data treatment and statistical analysis.A user manual and a description of the R language can be found on the website http://www.r-project.org/.The main interest of such a language remains in that it is distributed under GNU GPL license which does R routines an open source program, available on various platforms (i.e.GNU/Linux, FreeBSD, NetBSD, OpenBSD, Mac OS and Windows).
The simulator is composed by different blocks (Fig. 2): an input block which contains the different elements mandatory for the processing; a processing block where the user can choose which algorithm to be used, and an output block containing the different results of the simulation.As inputs, this simulator requires the receiver coordinates, the satellite ephemeris and a set of optional environmental parameters such as a DEM in order to take the possible masking of the terrestrial topography into account, as well as adaptive mapping functions to integrate atmospheric delays and bending effects.
As outputs, the simulator provides the time-varying reflection point coordinates, but also various KML files (Keyhole Markup Language -standard format used by Google Earth) such as satellites positions and ground paths, waves trajectories and Fresnel first surfaces which can be opened using the Google Earth visualization tool.

GNSS orbit parameters
The simulations are based on the determination of the positions of the specular reflection points, once the receiver and the satellites positions are known.Satellites coordinates can be obtained from the International GNSS Service (IGS) ephemeris final products which provide GNSS orbit and clock offset data with a temporal resolution of 15 min in the SP3 format for the past epochs, or derived from the Keplerian parameters (semi-major axis, inclination, and argument of perigee) to predict GNSS satellite positions.These products are available on the IGS website: http://igs.org/.

Radio-electric mask
Simulations are performed for a given receiver position in the WGS84 coordinates system and height above the ground.It is possible to apply an elevation or azimuthal mask to the simulations to avoid low elevation satellites for instance.The elevation mask commonly used is set to (10 • ; 90 • ) and no mask is set in azimuth.

SRTM Digital Elevation Model
The most realistic simulation needs the integration of a Digital Elevation Model (DEM) in order not to only take the possible masking of satellites into account, but to get more accurate and exact positions of the specular reflection points as well.The holefilled version 4 of the Shuttle Radar Topography Mission (SRTM) DEM, with a spatial resolution of 90 m at the equator is used (Jarvis et al., 2008).The altitudes are given with reference to the EGM2008 geoid model.Uncertainty on altitude is around 16 m over mountainous areas (Rodriguez et al., 2005).It is made available by files of 5

Earth Gravitational Model EGM2008
In order to be able to convert between ellipsoidal heights (with respect to the WGS84 ellipsoid) and altitudes (with respect to the EGM2008 geoid model) when producing KML files or when integrating a DEM, the knowledge of the geoid undulation is mandatory.The Earth Gravitational Model EGM2008 has a spatial resolution of 10 km (Pavlis et al., 2012).In this study, we used 2.5 min × 2.5 min Geoid Undulation Grid file derived from EGM 2008 model in a tide-free system released by the US National Geospatial-Intelligence Agency (NGA) EGM Development Team: http: //earth-info.nga.mil/GandG/wgs84/gravitymod/.

Adaptive mapping functions
The neutral atmosphere changes the propagation direction and bends the propagation path of the GNSS signal.The range between the satellite and the tracking site is neither the geometric distance nor the length of the propagation path, but the radio range of the propagation path (Marini, 1972).
For GNSS-R measurements, the tropospheric delays induced by the neutral part of the atmosphere are an important source of error.Indeed, GNSS-R measurements are often done at low elevation where the tropospheric effects are maximal.Accurate models of tropospheric delays have to be used to mitigate signal speed decrease and path bending.It is commonly accepted to model tropospheric delays by calculating the zenith tropospheric delay and obtaining the slant tropospheric delays with a mapping function.New mapping functions have been developed in the 2000's (Boehm et al., 2006;Niell, 2001) and significantly improve the geodetic positioning.Although modern mapping functions like VMF1 and GPT2/VMF1 are derived from numerical weather models (NWM), most of these mapping functions ignore the azimuth dependency which is usually introduced by two horizontal gradient parameters -in north-south and eastwest directions -estimated directly from observations (Chen et al., 1997).More recently, the use of ray-traced delays through NWM directly at observation level has Figures

Back Close
Full shown an improvement on geodetic results (Hobiger et al., 2008;Nafisi et al., 2012;Zus et al., 2012).The Adaptive Mapping Functions (AMF) are designed to fit the most information available in NWM -especially the azimuth dependency -preserving the classical mapping function strategy.AMF are thus used to approximate thousands of atmospheric ray-traced delays using a few tens of coefficients with millimetre accuracy at low elevation (Gegout et al., 2011).AMF have a classical form with terms which are function of the elevation.But, they also include coefficients which depend on the azimuth to represent the azimuthal dependency of ray-traced delays.In addition, AMF are suitable to adapt to complex weather by changing the truncation of the successive fractions.Therefore, the AMF are espacially suited to correct propagation of low elevation GNSS-R signals.

Data used for validation
In order to assess the simulator performance and the ocean tide influence on the positions of the reflection points estimated at an offshore experimental site located at the top of the Cordouan lighthouse (45 • 35 11 N ; 1 • 10 24 W), we use 24 h of REFMAR (Réseau de Référence des Observations Marégraphiques) tide gauge observations, with a sampling frequency of 5 min.The tide gauge records of the station of Royan (45 • 37 14.07 N; 1 • 01 40.12 , located 12 km from the lighthouse) are the property of MEDDE (Ministère de l'Ecologie, du Développement Durable et de l'Energie), and they are available on the REFMAR website (http://refmar.shom.fr).

Methodology: determination of the positions of reflection points
The difference of phase between the two antennas (A-RHCP and B-LHCP on Fig. 1) at an epoch t for the i th GNSS satellite can be seen as a classical simple difference between two receivers used for relative positioning as follows: Tables Figures

Back Close
Full where λ is the wavelength of the GNSS wavelength carrier, ∆φ i AB the measured carrier phase difference between the direct and received signals expressed in cycles, ∆δ i AB the difference in distance between the direct and received signals, ∆N i AB is the difference of phase ambiguity between the direct and received signals, c the speed of light in vacuum, ∆t AB the receivers clock bias difference.As the baseline between the two receivers is short (a few centimeters to a few tenth of centimeters), and in the case of low altitude of the receivers (typically from a few meters to a few hundred meters in the case of in situ experiments), both tropospheric and ionospheric effects can be neglected as they are cancelled out by single difference.Besides, when both antennas are connected to same receiver, the receiver clock bias difference is also cancelled out.
In this study, we only consider the difference in distance between direct and reflected signals as illustrated in Fig. 1.
The processing block contains three algorithms for determining the positions of the specular reflection points: the first considering the Earth as a sphere, the second as an ellipsoid, and the third one takes the Earth's topography into account.As it will be discussed in the Sect.4.2, the three algorithms have different characteristics, in terms of calculation time and accuracy of the positions determination.
All of them are based on iterative approaches to solve the Snell-Descartes law for reflection: the angle of incidence is equal to the angle of reflection on a plane interface separating two half-space media.

Local sphere approximation
Let us consider the vertical plane formed by the transmitter (GNSS) satellite (T ), the receiver (R) and O, the centre of the Earth (Fig. 3a).We assume that the specular reflection point (S) will be included in that plane.Let us consider the following orthonormal reference systems of coordinates: -R2(O, x, y): a local two-dimensional system, obtained by the rotation of the (O, X , Y , Z) system around the z axis, in such a way that x r = 0.
-R3(S, x , y ): a local two-dimensional system, obtained by a rotation around the z axis and a r E translation of the (O, x, y) system in such a way that x and the local vertical are colinear, and that the system origin coincides with the specular reflection point S.
If H is the height of the receiver above the ground, the position of the receiver is: (2) with 1 − e 2 sin 2 (ϕ r ) (3) the radius of the Earth at the latitude of the receiver ϕ r .e being the eccentricity, and a the semi-major axis of the WGS84 ellipsoid.
The position of the GNSS satellite transmitter considering ε the elevation angle of the satellite and τ the angle RT O is given by (Fig. 3b): Using the trigonometric sine formula in the R-T -O triangle (Fig. 3b):

Conclusions References
Tables Figures

Back Close
Full We finally obtain: The Snell-Descartes law for reflection can be expressed as the ratios of the coordinates of the receiver and the transmitter in (S, x , y ): The coordinates in R3 can be derived from the coordinates in (O, x, y) from: where γ is the rotation angle between the two systems (Fig. 3a).So Eq. ( 7) becomes: 2(x t x r − y t y r ) sin(γ) cos(γ) − (x t y r + y t x r )(cos 2 (γ) − sin 2 (γ)) − r e (x t + x r ) sin(γ) + r e (y t + y r ) cos(γ) = 0 (9) Following (Helm, 2008), we proceed to the substitution t = tan( γ 2 ), and Eq. ( 9) becomes: And finally becomes: with: Equation ( 11) is solved to determine the roots of this polynom using an iterative scheme based on the Newton method (Nocedal and Wright, 2006).

Local ellipsoid approximation
By knowing the locations of the transmitter and the receiver on the local ellipsoid included in the plane defined by the centre of the Earth, the receiver and the transmitter, let us consider the two normalized vectors between the specular reflection point and the transmitter, and the specular reflection point and the receiver.When the Snell-Descartes law is verified, the sum of the two vectors coincides with the local vertical (Fig. 4).The determination of the location of the reflection point is based on the following iterative process proposed earlier by Gleason et al. (2009): with K a coefficient and the correction in direction.The iterative process stops when the criteria: α is a parameter of tolerance.After several tests, we chose α = 1 mm, which is enough is most cases.To ensure a faster convergence, K can vary as a function of the difference of the distance between two consecutive steps.Following Gleason et al. (2009), we chose K equals 10 000 when r s (t + 1) − r s (t) ≥ 10 m and K equals 1000 afterwards.

Local ellipsoid approximation taking a DEM into account
The two first approaches presented above are well adapted in the case of an isolated receiver, located on the top of a light house, for instance.In most of the cases, the receiver is located on a cliff, a sand dune, or a building overhanging the sea surface or a lake.It can however be really judicious to incorporate a Digital Elevation Model (DEM) into the simulations, in order not to only take the mask effects (e.g., a mountain occulting a GNSS satellite) into account, but also to get more accurate and realistic positions of specular reflection points.The method we propose here consists of three steps later detailed in Sects.3.3.1-3.3.3: 1. a "visibility" determination approach to determine if the receiver is in sight of each GNSS satellite.
2. A transformation from 3-D coordinates to 2-D coordinates to increase the computation velocity.
3. A determination of the specular reflection point position.
We have to keep in mind that a DEM gives altitudes above a reference geoid.For consistency purpose, the positions of the receiver and the transmitter, and the DEM 1013 Introduction

Conclusions References
Tables Figures

Back Close
Full grid points have all to be in the same reference system.So it is absolutely mandatory to convert the altitudes of the DEM grid points into ellipsoidal heights by removing the geoid undulation.To do so, a global grid from the EGM2008 gravity field model with respect to the WGS84 ellipsoid was removed from SRTM DEM grid points.

Visibility of the GNSS satellite from the receiver
This algorithm aims to determine the presence of mask between the receiver and the satellite.The coordinates of the DEM points, the satellite and the receiver locations are already known in a 3-D Earth Centred Earth fixed (ECEF coordinate system).To simplify computation, the first step is to convert these 3-D coordinates in a local 3-D East North Up (ENU) system (Fig. 5).A similar approach in 3-D ECEF coordinates could have been possible (i.e., without projecting the positions of the transmitter and the receiver in the ECEF system), at the expense of longer and more complicated computations for finally a small gain in accuracy.As an example of that coordinates change, we consider a topography formed by the A, B, C, D and E points (Fig. 6).In the case, a satellite located in T is not visible from the receiver R due to the elevation of the topography (the direct signal encounter the topography in P at an elevation lower than the elevation of B).Points A, B, C, D and E , and also R (receiver) and T (transmitter) are known in the 3-D ECEF system.In order to transform these coordinate into a local 3-D ENU system, we proceed as follow: 1. an ENU system centred in the receiver projection on the ellipsoid is defined, such as x ENU and y ENU be included in the tangent plane to the ellipsoid, and z ENU be up-looking (the east, north and up components of the receiver in the ENU system).Introduction

Conclusions References
Tables Figures

Back Close
Full 4. The ellipsoidal height is added to these new points as attribute.We thus obtain 3-D ENU coordinates of the initial points.These points are noted T , A , B , etc. in Fig. 6.Obviously, the satellite in T still remains hidden from the receiver in R.
This coordinate transformation done, we can now begin to focus on the visibility problem.The algorithm consists of the following steps: 1. a buffer zone where to find the location of the reflection point is first defined to reduce both the computation time and the allocation of the memory occupied by the DEM.First, all the DEM grid points outside the zone delineated by the perpendiculars lines to the emitter-receiver segment are eliminated (regions 1 and 2 in Fig. 7).Then, a buffer zone is defined around segment [T R] (region 3 in Fig. 7).Here we chose 1.5 times the DEM resolution (i.e.1.5 m × 90 m) which assures to have enough points to interpolate altitudes.
2. The segment [T R] is run through according to a predefined sampling rate (here chosen equal to the DEM resolution, i.e. 90 m).For each step, the ellipsoidal heights interpolated from the DEM using the nearest neighbours method and derived from the satellite and receiver locations are compared.If the one interpolated from the DEM is higher than the one calculated from the satellite and receiver, it means that the satellite is not visible from the receiver.Conversely, if any point of the segment has an interpolated ellipsoidal height higher than the one calculated from the satellite and receiver coordinates.

Transformation from 3-D coordinates to 2-D coordinates
Once system origin (i.e., the receiver), and the corresponding ordinate will be determined by interpolating the ellipsoidal height from the nearest neighbours in the previous 3-D system.The topography is represented by a broken line between the receiver which abscissa equals to zero, and the satellite which abscissa equals to the planimetric distance between the receiver and the satellite.Consequently, the previous 3-D ENU system is now used as an intermediate between the 3-D ECEF system and this new 2-D system.

Position of the specular point
Once the satellite visibility from the receiver is confirmed, the last step consists in determining the location of the specular reflection point in the DEM along the 2-D profile.We suppose that the specular point is located into the plane formed by the satellite, the receiver and the centre of the Earth.The method is based on a predetermined part of the DEM at a sampling resolution of 10 cm on first estimation, and 1 mm afterwards for a better determination of the position according to the first law of Snell-Descartes.The main issue here is to pre-determine judiciously the part of the DEM susceptible to contain a specular reflection point.First, we compute the specular point disregarding the DEM, with sphere or ellipsoid approximation algorithm.We thus obtain the position of a point noted S1 (Fig. 8a).Given that the distance between the satellite and the receiver is really huge (about 20 000 km), we can approximate that waves emitted by that satellite placed at infinity are parallel between them (i.e. blue straight lines in Fig. 8a).In order to respect the Snell-Descartes law, locations of potential specular points are the intersection points between the line (R-S1) and the DEM.Following the example of the Fig. 8a, this would correspond to the point P 1.Let us consider now the other extreme, basing ourselves on the point of view of the direct wave (T -R).We can also suppose, on first approximation, the waves to be parallel (i.e.red straight lines in Fig. 8b).Disregarding the DEM, we thus obtain a second specular point S2.As previously done, according to the Snell-Descartes law, the sole places where reflected and incident angles would be equal will be the 1016 Figures

Back Close
Full intersection points between the DEM and the line (R-S2); in other words, points P 2, P 3 and P 4 in Fig. 8b.The potential specular reflection points will thus be located between these two extremes that correspond to green hatched area in Fig. 8c.Hence it will only be necessary to run through the part of DEM included in that area to determine the specular points.It will obviously be mandatory to check if these specular points are visible from the satellite and the receiver before, thanks to a visibility determination algorithm using 2-D coordinates.

Tropospheric corrections
In order to correct the anisotropy of propagation of radio waves used by the GNSS satellites, we use AMF calculated from the ECMWF model-level data.Given the geometric specificities of the specular reflection point, two paths have to be checked for propagation error: the first one from the satellite to the ground, and the second from the ground to the receiver.

Correction of the satellite-ground path
First and foremost, we solve the parallax problem for the wave emitted by a known GNSS satellite.At first sight, we consider the position of the specular reflection point calculated without any tropospheric correction, given by the algorithm approximating the Earth as a sphere given in Sect.3.1.We use here AMF calculated from a corresponding receiver ground position (i.e. the receiver position minus the receiver height), considering that the AMF planimetric variations are negligible for ground-based observations (i.e.we consider that we can use the same AMF for every specular reflection points, which is valid only if the specular reflection points are not too far from the receiver).We thus obtain the corrected elevation of the incoming wave.Considering the law of Snell-Descartes, the reflecting angle must be equal to the corrected elevation, for the specular reflection point position.

Correction of the ground-receiver path
With the corrected reflection angle, we can deduce the corrected geometric distance between the reflection point and the receiver, using this time AMF calculated from the receiver, assuming that the AMF altimetric variations are non-negligible (i.e. the part of the troposphere corresponding to the receiver height will have a non-negligible impact on the AMF).Knowing the corrected geometric distance between the reflection point and the receiver, the corrected position of the reflection point is obviously determined.It is indeed obtained by intersection between a circle whose radius is equal to the correct geometric distance, and the surface of the Earth assimilated as a sphere, an ellipsoid, or with a DEM, depending on which approximation of the Earth is taken into account.
We iterate the whole process until convergence to reach a better accuracy of the reflection point position.In fact, the first corrections were not accurate since calculated from an initially false reflection point position, and each iteration brings the point closer to the correct position, diminishing each time the correction to apply.Figure 9 shows an example of elevation corrections to apply as a function of the satellite elevations.This figure has been computed from simulations done on a receiver placed on the Geneva Lake shore (46 • 24 30 N; 6 • 43 6 E; 471 m): see Sect.4.2.

Footprint size of the reflected signal
The footprint of the reflected signal, referred as the glistening zone, corresponds to the area around the specular reflection point where the signal power is being scattered towards the receiver.It is defined by the intersection of iso-range (i.e.ellipses of equal delay) and iso-Doppler (i.e.paraboles of constant Doppler frequency) contours (Gleason, 2006;Helm, 2008).The signal power received is mostly due to coherent reflection and most of scattering is coming from the first Fresnel zone (Peckmann and Spizzichino, 1987).The first Fresnel surface can be described in the horizontal plane Introduction

Conclusions References
Tables Figures

Back Close
Full as an ellipse of semi-major axis (a) and semi-minor axis (b) equal to: With λ the wave length (m), h the receiver height (m) and the satellite elevation seen from the specular reflection point (rad) (i.e.corresponds to the reflection angle).
4 Simulator performance and results

Plot of the specular reflection points and recap text files
The simulator provides the position of the reflection points estimated during the selected time period of the simulation for each satellite, with a time-step of 15 min.These successive positions are mapped gradually on a pop-up window of the R software and their coordinates are contained in a text file which summarizes the different selected parameters of the simulation, as well.

KML files
The coordinates of the simulated specular reflection points are provided as KML files too: it is possible to use Google Earth to visualize them.This allows us to use the Google Earth time-selection cursor to visualize the simulation results either at every pre-step ∆t (i.e., every 15 min), or cumulated over longer timer period ∆T = n i =1 ∆t i .The different KML files created at the end of each simulation and viewable in Google Earth are the following: 1019 Introduction

Conclusions References
Tables Figures

Back Close
Full first Fresnel surface.

Simulation experiments
Simulations and tests of parameters have been performed on two main sites: the Cordouan lighthouse (45 • 35 11 N ; 1 • 10 24 W), in the Gironde Estuary, France.This lighthouse is about 60 m high, and it is surrounded by the sea.
-The shore of the Geneva lake ( 46• 24 30 N; 6 • 43 6 E).This site is hidden by mountains in the South (altitude up to 2000 m), and overlooks the lake in the North (altitude of 370 m).
For both sites, precise GPS and GLONASS ephemeris have been taken from the IGS website (15 min sampling interval).

Outputs
Examples of visualization of outputs for simulations in the case of the Cordouan lighthouse are respectively presented in Figs.10-14b.These simulations have been done considering the sphere approximation algorithm and a 15 min time-step.Introduction

Conclusions References
Tables Figures

Back Close
Full Figure 15 shows the variation of the distance between reflected points and receiver, as a function of the satellite elevation, and for several receiver heights and Fig. 16 shows the variation of the area of the first Fresnel surface.Such figures have been produced by doing simulations on the Cordouan lighthouse and varying the receiver height.The map of the reflected points obtained for an important receiver height will in fact be the same as the one obtained for a smaller receiver height, but more stretched.Henceforth, the higher the receiver height, the bigger the "measurable" area, but the less dense the ground coverage of the data (less reflection points per surface unit).

Assessment of the ocean tide influence
Simulations in the Cordouan lighthouse have been achieved integrating ocean tide from the tide gauge in Royan, by time-varying the receiver height in order to simulate the tide.The vertical visibility mask was set to 10-90 • , in order to avoid the weaker accuracy of determination of the specular reflection points positions for low elevation satellites, as highlighted in Sect.4.3.2.By comparing the results with simulations made with a fixedreceiver height of 60 m above the sea surface, it appears that the 3-D gaps reach values higher than 12 m for the maximum tide values (< 3 m) (Fig. 17).We can expect even higher discrepancies by taking into account satellites whose elevation would be lower than 10 • .

Geneva Lake
Three sets of simulation have been performed in the case of the Geneva Lake shore, for a 24 h experiment, on 14 October 2012: first configuration considering a receiver height of 5 m a.s.l.; second configuration considering a receiver height of 50 m a.s.l.; third configuration considering a receiver height of 300 m a.s.l. as for an airborne experiment (e.g.hovering helicopter).Introduction

Conclusions References
Tables Figures

Back Close
Full Each series has been computed using the three algorithms of determination of the reflection points (sphere approximation, ellipsoid approximation and the algorithm taking a DEM into account).Results are presented in Tables 1-6.They show the distances between the specular points and the receiver (arc lengths), and the differences between the positions given by each algorithm.

Influence of the receiver height
It appears that the differences between the methods used increase with the receiver height.Indeed, for a 5 m receiver height and satellite elevation greater than 5 • , the mean difference between the reflection points positions is 1.8 m whereas for a receiver height of 300 m, we reach 18 m.This is explainable by the fact that the higher the receiver is, the farther the reflection points will be from the receiver, and the bigger the impact of the Earth approximation will be.For a 5 m receiver height, reflection occurs until 57.1 m from the receiver, whereas for a 300 m receiver height, it occurs until 3302.5 m.It means that, in the second case, reflections occur in the mountains in the South of the receiver hence big differences between ellipsoid or sphere algorithms and the algorithm taking the DEM into account.

Influence of the satellite elevation angle
Secondly, by plotting the differences as functions of the satellite elevations, we can observe that the lapses between the different algorithms vary in an inversely proportional way than the satellite elevation (and so, proportionally to the point distance from the receiver).That is why we re-ran the simulations, putting a more restrictive mask of visibility, tolerating only satellites whose elevation is between 10 • and 90 • .Tables 4-6 show results we obtain by applying such a mask.By doing so, we get arc lengths smaller: about 1 m for a height of 5 m and 13.9 m for a height of 300 m (respectively 1.8 m and 18 m applying a 5 • elevation mask).The lower the satellite elevation is, the farther the specular reflection points from the receiver and the bigger the impact of Introduction

Conclusions References
Tables Figures

Back Close
Full the Earth approximation is.The choice of the algorithm used to perform the simulations becomes thus really important for the farthest reflection points (i.e for low satellite elevations, and high receiver height).

Influence of the DEM integration
Integrating a DEM has deleted 162 specular reflection points out of the 905 points determined during 24 h the 21 October 2012 with the sphere approximation algorithm (Fig. 18a).These 162 points came from a wave emitted by a satellite hidden by a mountain located in the south part of the area.In the north part, any reflection point is valid when taking a DEM into account, because in that direction, the relief is flat over the Geneva Lake, and so, satellites are all visible and reflections are possible (Fig. 18b).
Moreover, the points positions have been rectified while taking a DEM into account, since the others algorithms consider that reflections occur (in first approximation) in a plane around the projection of the receiver and without integrating the problem of the presence of relief.

Comparison between algorithms
For a 5 m receiver height, and for satellite elevations greater than 10 • , the mean difference between the ellipsoid and the sphere algorithm is equal to 8 cm whereas for a 300 m receiver height it is equal to 2.24 m.The approximation done by considering the Earth as a sphere or as an ellipsoid does not really affect the precision of the specular reflection point determination when reflection does not occur too far from the receiver (maximum equal to 60 cm for a distance inferior to 27 m) i.e. for low receiver height and high satellite elevation.When reflections occur far from the receiver, the choice of the approximation begins to be important: maximal differences of 6.8 m for distances inferior to 1672 m.Concerning the algorithm taking the DEM into account, the differences obtained with respect to the sphere or ellipsoid algorithms are quite big even if the specular reflection Introduction

Conclusions References
Tables Figures

Back Close
Full point is close enough from the receiver.For instance, the mean difference between the sphere or ellipsoid algorithm and the one integrating the DEM is bigger than 2.2 m for a 5 m receiver height, and bigger than 24 m for a 300 m receiver height, and with satellite elevation above 10 • .
It is also worth noticing that, globally, the ellipsoid approximation gives smaller differences with the algorithm integrating the DEM than the sphere approximation.

Tropospheric error
Given the geometric configuration of the satellite, the reflection point and the receiver, the same elevation correction will have a different effect according to the receiver height.It turns out that considering a same satellite at a given time, the corresponding reflection point will be farther for a big receiver height than for a smaller one.Consequently, for the same elevation correction, the resulting correction of the reflection point position will be higher in the first case than in the second one.Figure 19 shows the differences, in terms of geometric distances, between the reflection points positions obtained with and without taking the tropospheric correction into account (delay and bending) and for different receiver heights (with the ellipsoid approximation of the Earth).It appears that for low satellite elevation and high receiver height, the tropospheric error has a non-negligible influence on the specular point positions (103 m for a 300 m receiver height, satellites elevation inferior to 10 • ).

Calculation time
An assessment of the simulator performance has been achieved in terms of computation time from runs computed with the following two computers characteristics: -Computer 1: 2 GB ram, dual core processor Quad 2.66 GHz, using Windows 7, 64 bits operating system.Introduction

Conclusions References
Tables Figures
The different series of simulations have been processed with receiver heights of respectively 5, 10, 30, 50, 100, 300 and 500 m and during 24 h, the 21 October 2012.Each series has been processed 10 times and averaged for both computers, and with the three different algorithms.The results of such analysis are visible on Fig. 20.It is worth reminding that such a factor will highly be influenced by both the capacities of the processor used to do the calculations, and the chosen parameters to reach a precise estimate of position (notably in terms of convergence criteria and tolerances).
The major part of the calculation time is due to the conversion from ellipsoidal heights to altitudes (interpolation from a grid).Regarding the results, we can firstly notice that the sphere approximation algorithm is the fastest, followed by the ellipsoid approximation and then the algorithm taking a DEM into account, which is logical, given the operating mode of each algorithm (equation storage for the first one, iterative process for the two others).The receiver height does not affect calculation time for the sphere approximation, whereas it increases it while integrating a DEM.This is explained by the fact that when the receiver position becomes higher, the reflection area increases and so the DEM part to analyse becomes bigger.Regarding the ellipsoid approximation algorithm, calculation time is high for very small receiver heights, and reaches a minimum for a 30 m receiver height, and then slowly increases with the receiver height.The peak for small receiver heights can be explained by the fact that the coefficient K (shifting factor of the temporary position between two iterations) used during the iterative process (see Sect. 3.2) varies proportionally to the correction applied to the temporary position.For small receiver height, this correction will be small (because the reflection point will not be far away from the receiver), and consequently also the K coefficient, hence a slower (but more precise) convergence.Introduction

Conclusions References
Tables Figures

Back Close
Full

Conclusions
In this paper, we presented a simulator based on real GNSS satellite ephemeris, as a user-friendly tool, for modelling the trajectories of GNSS electromagnetic waves that are reflected on the surface of the Earth and therefore preparing GNSS-R campaigns more efficiently.The originality of this simulator remains mainly in the integration of a DEM and of the tropospheric error correction.The results of simulations led us to a better understanding of the influence of some parameters on the reflection geometry, namely by quantifying the impact of the receiver height but also the influence of the satellite elevations, the natural relief (DEM), and the tropospheric bending and delay.
The different simulations realized near to quite rugged topography lead us to the following conclusions: the DEM integration is really important for mountainous areas (differences up to 544 m for a 300 m receiver height, elevation superior to 10 • ).
-Differences between sphere and ellipsoid approximation are negligible for specular reflection point close from the receiver (inferior to 60 cm for a 5 m receiverheight, elevation superior to 10 • ) i.e. small receiver height and/or high satellites elevations.
-The tropospheric error correction is absolutely mandatory for reflection points farther than 55 m from the receiver (i.e.receiver height inferior to 5 m, satellite elevation superior to 5 • ) but can be negligible otherwise.The correction to apply is exponential for low elevation satellites.
Globally, it is worth reminding that the farther the specular reflection point from the receiver is, the more important the influence of the different error sources will be: Earth approximation, DEM integration, tropospheric error correction.The farthest specular reflection points will be obtained for high receiver height and low satellite elevation.Introduction

Conclusions References
Tables Figures

Back Close
Full implemented, such as receiver installed on a moving platform in order to map the area covered by airborne GNSS-R measurements campaigns and on-board a LEO satellite.Introduction Conclusions Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

-
R1(O, X , Y , Z): WGS84 Cartesian system (NIMA, 1997), with O the centre of the Earth.The receiver and transmitter coordinates are known in this system.Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | the visibility of the satellite is checked, a change of coordinates from 3-D to 2-D is achieved to decrease the computation time of the determination of the location of the specular points.The origin of this 2-D system is defined by the receiver, the abscissa axis is formed by the [RT ] planimetric segment and the ordinate is the ellipsoid height.The abscissa of the profile are thus given by the planimetric distance from the Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | positions of the specular reflection points; positions of the receiver and satellites; ground paths of the satellites; direct and reflected waves; Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |