Mapping GPS Multipath : a Case Study for the Lunar Laser Ranger Timing Antenna at HartRAO

Accounting for multipath in Global Navigation Satellite Systems (GNSS) is a difficult task and an important one, especially during the pre-investigation phase for the installation of a permanent GNSS station for positioning or timing applications. Sites with a high level of multipath can cause positioning errors or timing errors resulting in the quality of GNSS products (position or timing) becoming degraded by several metres or nanoseconds. We investigate and attempt to map multipath as part of the site investigation for the installation of the timing antenna for lunar laser ranging applications at the Hartebeesthoek Radio Astronomy Observatory (HartRAO). A high-resolution wavelet power spectrum and a standard deviation parameter are used to map multipath in both the time and frequency domain as well as spatial variations on the sky plot. The high standard deviation values on the sky map are attributed to reflections due to shrubs or trees on the site, while smaller standard deviation areas are attributed to bare soil or less vegetated as this would give constant reflection over time provided the ground has constant moisture. We conclude that the site is suitable for installation of the timing antenna and that a mask of 15°-20° elevation angle will be applied to the timing antenna to minimise multipath at lower elevations.


Introduction
Multipath signals in most Global Navigation Satellite Systems (GNSS) applications are considered noise.Multipath occurs when an electromagnetic signal arrives at an antenna along an indirect path due to the deflection and reflection of the signal by nearby objects or the surface near the antenna (Larson et al., 2008).Multipath signals are caused by surrounding features such as trees or soil reflecting transmitted signals from the satellites to the antenna.Trees affect the frequencies in the L-band, the leaves and the trunk cause shadowing and attenuation, scattering originates from the tree branches, reflection and diffraction are caused mainly by the tree trunk (Schubert et al., 2010).Soil moisture changes affect the phase of the Signal-to-noise ratio (SNR) modulation pattern and its magnitude when a GPS satellite decends or ascends at low elevations (Zavorotny et al., 2010).The Lband signals penetrate further in dry soil as compared to wet soil (Njoku and Entekhabi, 1996).
The errors caused by multipath on GNSS data adversely affect GNSS products, these products include positioning (Closas et al., 2009), timing (Ray and Senior, 2003) and Integrated Precipitable Water Vapour (IPWV) (Tregoning et al., 1998).However, recent studies have shown that multipath can be used in the GNSS-Reflectometry (GNSS-R) technique to study soil moisture content around GNSS antennas (Larson et al., 2010;Mironov and Muzalecsjiy, 2012), snow depth estimation (Botteron et al., 2013), and in the estimation of height and moisture content of the vegetation around GNSS antennas (Wan et al., 2015).Application of GNSS-R is also important in characterising sites for geodetic applications, where accuracy is of high priority such as time transfer techniques using GNSS (Ray and Senior, 2003).
Radio signals transmitted by the GNSS satellites are also affected by the ionosphere.The ionosphere is the upper layer of the atmosphere about 60km from the Earth's surface and it contains free ions and electrons that are produced by solar radiation (Ratcliffe, 1972).The ionosphere causes a delay of radio signals.The ionospheric effect can result in a delay of up to 30ns (Giffard, 1999).However, the error introduced by the ionospheric effect can be minimised by using dual-frequency receivers and applying the double differencing method (Pullen et al., 2009).The troposphere also introduces a delay of up to metres of the radio wave.The delay introduced by water vapour in the troposphere increases with a decrease in elevation angle.This delay is usually corrected by applying tropospheric delay models such as the GPT2w developed by Böhm et al. (2015).
HartRAO is currently building a Lunar Laser Ranging (LLR) station.The LLR station requires an accurate timing system, which must be steered by the Global Positioning System (GPS).The GPS satellites broadcast position and time information to the end-user.Satellites are constantly updated with accurate timing information from the ground-based network of Caesium and Hydrogen master clocks (Allan, 1980), which is referred to as time transfer and GPS plays an integral part in this regard.A geodetic observatory that requires precise accurate station time.Station time is required in order to be accurate enough (within nanoseconds of UTC) to be able to track earth-orbiting satellites or retro-reflectors located on the surface of the Moon with high accuracy.
The technique of LLR involves transmitting short laser pulses to the retro-reflectors that are located on the surface of the Moon.The round-trip time-of-flight of the laser pulses is measured using an event timer.The event timer requires a stable frequency (which can range from 5 to 100MHz) from the local timing system (e.g.rubidium clock).The local timing system drifts slowly over time and this requires accurate time updates from GPS to be updated and monitored.The measurement data are used to support and maintain the International Terrestrial Reference Frame (ITRF), determination of the Moon's orbit, and to study geodynamics and Earth-Moon dynamics.HartRAO aims to achieve millimetre-ranging accuracy.This will require all sub-systems (these include: pointing model, timing, laser, transmit and receive path efficiencies, etc.) to be accurate, for the minimisation of station error budget to better than centimetre level.A detailed review of the LLR technique is given in Munghemezulu et al., (2016) and the references therein.In this paper, we present an experimental study, which aims to map multipath around the new site for the timing antenna for the new LLR station located at the Hartebeesthoek Radio Astronomy Observatory (HartRAO).

Study Area and Data
HartRAO is located in a valley of Magaliesberg mountain ranges, 50km northwest of Johannesburg.HartRAO is collocated with four instruments that are fundamental for space geodetic techniques, namely: Very Long Baseline Interferometry (VLBI), Satellite Laser Ranging (SLR), Global Navigation Satellite Systems (GNSS) antenna and Doppler Orbitography and Radio-positioning Integrated by Satellite (DORIS).Currently, the observatory is building a Lunar Laser Ranging system and a new 15 m VLBI2010 Global Observing System (VGOS) antenna (Figure 1).
The site for the timing antenna was chosen according to the following considerations: (i) the site must be close to the control room to minimise signal delay due to cable length, as this reduces the systematic errors introduced by the cable and variations of the ambient air temperature, (ii) the site must have minimal objects that can introduce multipath such as trees and reflective objects such as the LLR enclosure and the 26m VLBI antenna.However, since some of the objects cannot be removed, an observational strategy of observing above certain elevation angle has to be established to minimise the effects of multipath.Figure 2 illustrates the visibility of the satellites at the HartRAO site, the centre point represents the location of the proposed site for the timing antenna.
An experiment was carried out using a Topcon antenna and receiver (GB1000) for a period of 10 days.Observations were scheduled for days when there was a clear sky as clouds often introduce random multipath errors (Rao et al., 2013).The raw observations files were converted from Topcon Positioning System (TPS) format to the Receiver Independent Exchange (RINEX) format as recommended by Gurtner (2007).

Methodology
A method for computing pseudo-range multipath errors (MP1 and MP2) was devised to assess the level of multipath in the vicinity of the GNSS antenna by making use of Translation Editing and Quality Check (TEQC) software (Estey and Meertens, 1999).Bilich and Larson (2007) indicated certain limitations associated with MP1 and MP2 derived from TEQC software to map multipath environments.One limitation is the use of Root Mean Square (RMS) error values to represent multipath.In the current work the RMS values that resulted are shown in Table 1.The RMS error values lack spatial components and as a result, multipath around the GPS station cannot be visualised.The second limitation involves MP1 and MP2 being a measure of pseudo-range precision, which is highly dependent on smoothing algorithms, the antenna used, and firmware changes of the receiver (Ray and Senior, 2003).A good example to visualise multipath is given by Ogaja and Hedfors (2007), where they developed a MATLAB program for creating colourised polar maps of high-frequency multipath using TEQC report files.(2007) for the derivation of Equation 1.However, this method works well when high data rates (1 s sampling interval) are available.
We make use of the MP1 and MP2 values derived from the TEQC software to further enhance multipath mapping by creating a gridded sky plot from the derived multipath values using a linear interpolation method and computing the wavelet power spectrum of the multipath time-series data to understand the magnitude in relation to the spatial location of the objects inducing multipath.
Multipath signals from GNSS are non-stationary signals since the frequency of the signals varies in space and time.This makes Wavelet Analysis (WA) suitable for studying multipath behaviour over time as compared to the traditional Fourier Transform (FT) method.We use the Morlet wavelet, which is given by Where  is the time and o  is the wavenumber (Torrence and Compo, 1998).To construct a picture depicting the amplitude of any feature with its variations over time and scale, a continuous wavelet transform ( n W ) defined as the convolution of n x  (MP1 and MP2) time series with a scaled and translated version of  is used to derive the wavelet power spectrum, which is given in Equation 1, Where N is the number of points in the time series at sampling time interval δt (30 seconds) and ( * ) represents the complex conjugate (Torrence and Compo, 1998).
The standard deviation is calculated for each pixel of the linear interpolated multipath matrix as given by Equation 4, 2 1
Where i x is the linear interpolated MP1 or MP2 values and x is the mean value of MP1 or MP2.These data are projected to azimuth ( ) and elevation ( ) values as a function of time i.e. sky or dry bare soil or less vegetated areas).However, it must be noted that this method is limited to the time span of the data since  will change depending on the number of observations available.Therefore, multipath analysis using this method requires long timeseries (e.g. more than a month) to capture all the seasonal variations due to changes in seasons, which will affect e.g.growth patterns of the trees and this will affect reflected patterns of the L1 and L2 signals.

Results and Discussion
The GPS satellites provide many data points, but these points are only concentrated at specific satellite tracks and some parts of the sky are not fully covered (e.g.areas between 90 and 270 azimuths have a poor coverage in Figure 3).The gridded plots in Figure 3 (C, B and D) provide a complete picture of the distribution of multipath as well as signal-to-noise ratio (SNR) around the GPS station.As illustrated in Figure 3, most multipath errors occur below 30 elevation angle and this is associated with low SNR ratio values.signals from satellites at low elevation angles travel longer in the atmosphere (i.e.troposphere) than when the satellites are at high elevations (Degnan, 1993).Hence, the observations below 30 elevations are noisier as compared to those at a higher elevation.The "spaghetti" features in MP1 and MP2 plots of Figure 3 at both low and high elevations indicating a variable atmosphere and these features could be attributed to clouds or cirrus clouds that vary in space and time.Infrared satellite cloud cover data indicated variable cloud cover during the time of the experiment.The data for the cloud cover can be retrieved from (http://www2.sat24.com/history.aspx?culture=en) for the given dates of the experiment.

Conclusion
Mapping multipath from GNSS signals is a difficult task given the many variables that change over a short and long period of time.These variables include moisture content of the soil or vegetation, canopy structure and clouds.The technique of multipath mapping can be very useful in preliminary site investigation for geodetic applications such as the installation of a new GNSS base station for timing or positional applications.We used RINEX observations for 10 days to map spatial location of multipath and its variation over time to aid in understanding the distribution of multipath before an installation of the timing antenna at HartRAO for LLR timing reference purposes.
High-resolution time-frequency wavelet analysis demonstrates that the MP1 and MP2 signals peak at different frequencies with different intensities when reflected by the same feature.This analysis allows for extraction of information such as the time of day that the GNSS station receive maximum multipath (altered by changing geometry of the satellites constellation), how long this condition lasts, and whether an observational strategy can be implemented.Through standard deviation maps, we conclude that the multipath on the site is mainly due to reflections from the mountains surrounding the geodetic site and that the site offers clear sky view of the satellites with minimal multipath from 30 and above in elevation.
The high standard deviation is attributed to reflections due to shrubs, vegetation and bare moist soil around the site.An observational strategy limiting the use of satellite below an elevation of 15° to 20° will be implemented as this will offer minimal multipath and still enable adequate satellite visibility.

Figure 1 .
Figure 1.Space geodetic techniques collocated at HartRAO.The timing antenna will be mounted on top of the control room.The continuous operating GNSS station (HRAO) is important as it acts as a base/control station in solving positions of other points

Figure 3 .
Figure 3. Enhanced visualisation of multipath and SNR data, the normal sky plot indicate satellite tracks across the sky while gridded plots indicate a complete picture distribution of multipath around the GPS station.

Figure 4 .
Figure 4.The multipath signals derived from the L1 band and its associated wavelet power spectrum

Table 1 .
Averaged Root Mean Square (RMS) calculated for all the satellites in view during specific days of the experiment.Each observation lasted for 24 hours.SNR) data that is reported in the standard RINEX file format, where SN1 and SN2 SNR data represents raw signal strength values as given by the receiver for the L1 and L2 phase observations.In the multipath environment, the SNR values measured by the tracking algorithm are a composite of direct and indirect signals.Equation [1]gives SNR (in an environment where only one indirect signal is introduced) as a function of multipath amplitude Am, direct amplitude Ad and multipath relative phase .Refer toBilich and Larson