Investigating GNSS multipath effects induced by co-located Radar Corner Reflectors

1 Introduction

Interferometric synthetic aperture radar (InSAR) technology has been widely applied to measure and map the Earth’s land surface deformation arising from natural or anthropogenic phenomena [1], [2], [3]. In contrast to other geodetic techniques, InSAR can remotely image large areas at high spatial resolution and at regular repeat intervals (several days). InSAR is complementary to other ground-based geodetic networks such as Global Navigation Satellite Systems (GNSS) sites, levelling or gravity benchmarks, which provide point-wise access to stable, absolute geodetic terrestrial reference frames such as the International Terrestrial Reference Frame (ITRF, [4]).

In order to make precise comparisons between InSAR and other geodetic measurement techniques and to fully integrate InSAR into geodetic reference frames, a network of ground-based artificial radar targets can be deployed [5]. Such radar targets provide a high-quality, temporally stable phase response in SAR imagery, which is required for reliable time series analysis of InSAR phase measurements to derive surface deformation [6].

Openly available Sentinel-1 data creates the opportunity to use InSAR for deformation mapping on national and continental scales [7], [8], [9]. To improve the generation of large-scale and seamless deformation maps and to integrate Sentinel-1 products into geodetic reference frames, there is an increasing trend for radar targets suitable for InSAR analysis to be installed in geodetic networks [10], [11], [12], [13], [14], [15]. While passive radar Corner Reflectors (CR) have been deployed in the field for this purpose since the early 2000s [16] and successfully used with various SAR missions [17], active radar transponder devices are also being increasingly used [18], [19]. Co-locating the radar targets with other geodetic sensors such as GNSS has the advantage that InSAR displacement measurements can be directly compared and validated with the GNSS measurements [20] which helps integrating the relative InSAR displacements into an absolute reference frame defined by GNSS [14], [21]. However, there are concerns about degrading GNSS signals when co-locating CR at continuously operating GNSS sites [13]. Metal structures such as CR situated in the direct vicinity of a GNSS antenna may induce multipath effects by reflecting additional GNSS signals, particularly if located in direct vicinity of the GNSS receiver antenna. Such multipath effects result from interference of the direct path of GNSS observations from GNSS satellite to receiver antenna with indirect paths caused by near and far field reflections in the surroundings of the antenna [22]. Multipath effects can be classified depending on the distance between GNSS antenna and the source of reflection. While far-field effects show short-periodic properties (up to half an hour [23]) and can be averaged out over long observation periods, near-field effects possess nonzero mean and long-period characteristics of up to several hours [24]. Strong multipath effects may result in oscillating patterns in the GNSS observations [25]. In the context of GNSS reflectometry experiments (e. g. used to measure soil moisture, snow depth or water height [26]), Fresnel zones are used to identify and describe the regions of active scattering [27]. However, using Fresnel zones to describe the full multipath environment of a GNSS site becomes increasingly complex when considering real-world scenarios with non-horizontal surfaces and multipath effects induced from several reflectors around the site [28], [29].

A number of strategies to reduce GNSS multipath effects have been proposed in the past decades. The most straightforward option for multipath mitigation is to prevent reflections in the vicinity of the GNSS antenna by appropriate site selection or usage of microwave-absorbing materials [30]. Other mitigation options focus on the hardware by improving the antenna design (e. g. polarisation-related antenna gain pattern, choke rings, ground planes, antenna arrays [31]) or the receiver architecture (e. g. correlator techniques [32]). In addition to site selection and hardware improvements, most other published mitigation options make use of GNSS data processing strategies to model and mitigate multipath effects. These approaches can be classified into in-situ calibration [33], [34], [35] including the use of Fresnel zones [29], wavelet-based multipath reduction [36], multipath mitigation using the signal-to-noise ratio (SNR) [25], [37], least mean square adaptive filtering [38], improved stochastic modelling [39], [40], computational electromagnetic modelling [41], sidereal stacking and filtering [42], [43], [44] and empirical site modelling using observation residuals [45], [46], [47]. The latter approach has several advantages including that it is applicable to all GNSS (in contrast for example to the sidereal filtering approaches), capable of detecting both near-field and far-field multipath signals, independent of GNSS equipment and data sampling rates and suitable for easy assessment of site quality. We will therefore use this approach to map multipath patterns in the surroundings of continuously operating GNSS sites in Section 3.3. We also perform sidereal filtering in Section 3.2 to investigate changes to SNR at a GNSS site equipped with CR.

There is limited published literature on operating radar targets such as CR at co-located GNSS sites. In a scenario where a CR was deployed 30 m from a continuously operating GNSS site, [13] observed no detectable increase in the average root mean square (RMS) of GNSS phase residuals but suggested carrying out more comprehensive multipath testing for different deployment scenarios. Despite this recommendation and to the authors’ knowledge, there is no published study investigating the multipath effects of radar CR on GNSS measurements for a coupled GNSS/InSAR co-location scenario. In this contribution we intend to fill this gap by assessing the effect on GNSS observations and resulting site coordinates of two metallic CR mounted directly to a GNSS monument. Furthermore, we will visualise the spatial pattern of GNSS multipath effects induced by the CR using GNSS observation residuals. Our investigations make use of a regional network of geodetic monitoring sites located in the Greater Sydney area, Australia, since 2016 and consisting of five continuously operating GNSS sites co-located with CR (see Section 2.1).

The paper is structured as follows. Section 2 presents the study area and data sets and summarises the data analysis techniques used for the different investigations. Section 3 consists of three subsections presenting the results of coordinate time series analysis, SNR value comparison and multipath stacking map generation, while the focus of investigations is on the latter subsection. In Section 4, we discuss our findings again with a focus on multipath stacking maps and their interpretation. The discussion includes a theoretical modelling of GNSS signal reflections for the specific CR scenario using Fresnel zones. Section 5 gives the major conclusions of this study.

2 Material & methods

2.1 Study area and GNSS data sets

In this study, we analyse GNSS data observed at 18 continuously operating sites in the Greater Sydney area, Australia. Four of the sites (namely CA07, CA10, CA17, CA19, referred to as CAxx in the following) were constructed in June 2016 and are equipped with two CR (see Figure 1 and [14]). One additional site, MENA (operating since 2010), was retrofitted with two CR mounted to the antenna pole on 25 May 2016 (DOY 146, see Figure 1b). This site is therefore particularly well suited for detecting differences in coordinates, signal-to-noise-ratio (SNR) or post-fit residuals caused by the CR deployment.

The CR used in this experiment are made of three identical 60 cm × 60 cm square plates composed of a welded stainless steel frame and aluminium mesh. The two CR at each site are aligned for ascending and descending orbital passes of SAR satellites and the azimuth and baseplate elevation is given in Table 1. The CR are mounted on the GNSS antenna pole with a vertical distance of approximately 50 cm from the top of the CR to the bottom of the GNSS antenna (at both the CAxx sites and site MENA). In addition to the five continuously operating sites equipped with CR, we perform GNSS analyses for several other sites in the area with the same antenna as the four CAxx sites and site MENA. Figure 2 gives an overview of the location of all sites used within this study, most of them forming part of the CORSnet-NSW network [48] operated by New South Wales (NSW) Spatial Services.

Figure 1

GNSS sites in the Sydney area with CR attached to the antenna pole. a) site MENA at Menangle, NSW, without CR (before 25 May 2016); b) site MENA with CR (since 25 May 2016); c) site CA07 at Wilton Park, NSW (operating since July 2016).

Table 1

Location of continuously operating GNSS sites with co-located Corner Reflectors (CR) and orientation parameters of CR aligned for ascending and descending Radarsat-2 passes. Orientation parameters were rounded to 0.5 deg for alignment in the field with the True Azimuth adjusted by magnetic declination to enable usage of a sighting compass as described by [49].

Site name	Longitude [deg]	Latitude [deg]	Ellipsoidal height [m]	True Azimuth [deg]		Baseplate Elevation [deg]
				Ascending	Descending	Ascending	Descending
CA07	150.639152	−34.221610	233.9	253.46	106.96	17.00	15.50
CA10	150.760197	−34.071943	157.4	253.46	106.96	17.00	15.50
CA17	150.691522	−34.046270	87.5	253.46	106.96	17.00	15.50
CA19	150.691241	−34.178878	182.7	253.46	106.96	17.00	15.50
MENA	150.743758	−34.126087	111.3	253.46	106.96	17.00	15.50

Figure 2

Location of 18 continuously operating GNSS sites in the Sydney area used in this study. The inset shows the location of the map extent with respect to the Australian continent (red box). Further details on the antenna types are given in Table 2.

Table 2

Overview of different GNSS analyses at continuously operating sites in the Sydney area including data periods and processing strategy.

Start DOY	Site(s) of interest	Sites for validation	Antenna type	Processing strategy	Comments
End DOY
2013:019	MENA	–	ASH701945E_M	Differential	Split into 1222 days before and 1222 days after CR deployment
2019:273			SCIS
2015:146	MENA	GLBN, MGRV, NWRA, SPWD, VLWD, WFAL	ASH701945E_M	PPP	Split into two one-year periods before and after CR deployment at site MENA, other sites have the same antenna as MENA
2017:145			SCIS
2017:121	CA07, CA10, CA17, CA19	BANK, FTDN, GERR, HNSB, KTMB, OBRN, PCTN	TRM57971.00	PPP	One-year period minimising the outage periods at CA19 and full coverage for CA07, CA10, CA17. Other sites have the same antenna as CAxx
2018:120			NONE

At site MENA, we perform three different investigations: (i) comparison of daily site coordinates resulting from differential GNSS processing of two periods consisting of 1222 days before and after CR deployment, (ii) comparison of stacked SNR values using sidereal lags for a 7-days period before and after CR deployment and (iii) comparison of stacked PPP ionosphere-free linear combination residuals (L3 residuals) for a one-year period before and after CR deployment. The latter analysis is also performed for six other continuously operating sites located in the Sydney area for validation. These sites are equipped with the same Ashtech antenna as site MENA, but without CR being attached to the monument (see Table 2).

At the CAxx sites, we perform residual stacking for a one-year period to visualise patterns related to multipath. Again, a set of other continuously operating sites located in the Sydney area with the same Trimble Zephyr antenna as the CAxx sites, but without CR, are used for validation. A different data period is used compared to the periods used for the multipath stacking at MENA in order to minimise outage periods present at site CA19. Table 2 gives an overview of the analysed periods used for the different investigations at the different sites. Note that the receiver types are not consistent for the different sites and the chosen data periods. The manifold receiver changes which have occurred at the various sites analysed in this study are summarised in Table 3.

Table 3

Overview of GNSS receivers used at the study sites during the analysed data periods.

Start DOY	Sites	Receiver type	Secondary receiver type (start date)
End DOY
2013:019	MENA	LEICA GRX1200GGPRO	LEICA GR30 (since 2017:326)
2019:273
2015:146	MENA, GLBN, NWRA	LEICA GRX1200GGPRO
2017:145	MGRV	LEICA GRX1200GGPRO	TRIMBLE NETR9 (since 2017:118)
	SPWD	LEICA GRX1200+GNSS
	VLWD,	TRIMBLE NETR9
	WFAL	TRIMBLE NETR5	TRIMBLE NETR9 (since 2017:044)
2017:121	CA07, CA10, CA17, CA19, FTDN, GERR, HNSB, KTMB, OBRN, PCTN	TRIMBLE NETR9
2018:120
	BANK	TRIMBLE NETR5	TRIMBLE NETR9 (since 2017:201)

2.2 Data analysis

In the following, we describe our processing strategies and data analysis tools. Note that we analyse GPS (Global Positioning System) data only even if data from other GNSS have been observed at some sites. For the comparison of daily site coordinates, we use the double difference differential data processing methodology implemented in the Bernese GNSS software (version 5.2, [50]). High-quality sites belonging to the International GNSS Service (IGS) and Asia Pacific Reference Frame (APREF) networks were used as reference sites [51], [52]. The most relevant processing parameters and models used for the differential processing are listed in Table 4.

Table 4

Overview of major processing parameters and models used for the differential and PPP-mode GPS processing at site MENA.

GPS processing parameters	Differential processing	PPP processing
Geodetic Datum	IGS14, epoch 2016.0	IGS14, epoch 2016.0
GPS observation data	30-second phase observations (L1, L2)	30-second phase observations (L1, L2)
Observation weighting model	Elevation-dependent ( sin 2 E)	Elevation-dependent ( sin 2 E)
Elevation cut-off angle	7 degrees	7 degrees
EOP / satellite orbit / satellite clock	Final IGS products (24 h / 15 min / 5 min)	Final IGS products (24 h / 15 min / 30 s)
Differencing level	Double	Zero
Ambiguity resolution	Yes (Sigma and QIF strategies)	No (float ambiguities estimated per epoch)
Ionosphere handling	Ionosphere-free linear combination L3	Ionosphere-free linear combination L3
Troposphere a priori model	Dry GMF	Dry GMF
Tropospheric mapping function	Wet GMF	Wet GMF
Site-specific troposphere parameters	Zenith delay (2 h) + gradient estimation (24 h, model: Chen and Herring [53])	Zenith delay (1 h) + gradient estimation (24 h, model: Chen and Herring [53])
Ocean / atmospheric loading models	FES2004 (Le Provost [54]) / Ray and Ponte 2003 [55]	FES2004 (Le Provost [54]) / Ray and Ponte 2003 [55]
Nutation / pole model	IERS convention 2010 [56]	IERS convention 2010 [56]
Satellite antenna correction	Absolute antenna model (IGS)	Absolute antenna model (IGS)
Receiver antenna correction	Absolute calibration (IGS type-specific)	Absolute calibration (IGS type-specific)

The SNR analysis at site MENA uses SNR values for L1 and L2 as given in the corresponding RINEX files. We then compare SNR for each visible satellite in the seven days before and after CR deployment. To reduce noise we take advantage of the ground track repeat period of GPS satellites and stack SNR values by applying the time lag of the sidereal day [42], [57].

For the generation of multipath stacking maps we use the post-fit residuals obtained from a PPP processing completed using the Bernese PPP processing engine. The most relevant processing parameters and models used for the PPP processing are listed in Table 4. The variable phase centre characteristics of the antennas installed at the CAxx sites were individually calibrated on the Geoscience Australia robotic antenna calibration facility [58]. However, we find that the differences in the PPP L3 post-fit residuals for a processing using IGS type-specific antenna correction values and using individually calibrated antenna corrections are small. The mean absolute difference of L3 residuals resulting from PPP processing of one day at site CA07 when using IGS or individual antenna correction values is 0.4 mm with a maximum difference of 9.2 mm, while the L3 residuals themselves reach absolute values of up to 126.5 mm on the analysed day. Since the antennas at the other sites used for validation have not been individually calibrated, we use IGS type-mean antenna corrections for all sites for consistency.

The spatial stacking is performed in congruent cells with a resolution of 1 deg × 1 deg (elevation × azimuth) at the local horizon, but lower resolution in azimuth with increasing elevation angles. This is done to enable approximately equal cell size at all elevation levels resulting in a more equal number of observations in each cell [47]. We also apply the statistical assessment of residuals within each stacking cell as described by [47]. To obtain a good spatial coverage of cells with valid stacked values, the stacking is performed for periods of a full year. Using shorter periods results in gaps in the spatial coverage of stacking cells due to the geometry of GPS satellite tracks (see Figure 3). The enhanced spatial coverage of grid cells containing a stacked value in Figure 3c enables a detailed interpretation of effects and provides information on the location of obstacles surrounding the site. For example, the view on the hemispheres in Figure 3 is from a northeasterly direction revealing oscillating patterns which indicate strong multipath effects in that direction [25]. Using a full year of data also reduces the impact of potential seasonal effects in the residuals on the stacked values. Residual stacking is performed for site MENA (and other sites with the same antenna) separately for the two one-year periods before and after CR deployment at MENA, as well as for a one-year period at the CAxx sites (and other sites with the same antenna, cf. Table 2).

Figure 3

3D view of stacked residuals at site SPWD (Springwood, NSW) showing the stacking cells on a hemisphere above the local horizon. Different time periods are used for residual stacking in this example: a) 10 days, b) 30 days, c) 365 days. The stacking is performed at congruent cells with a resolution of 1 deg in elevation and 1 deg in azimuth (at the local horizon, lower resolution for increasing elevation angles).

3 Results

In this section, we present the results of the described data analysis: (i) coordinate time series at site MENA (Section 3.1), (ii) SNR values at site MENA (Section 3.2), (iii) multipath stacking maps at all analysed sites (Section 3.3).

3.1 Daily site coordinates

The full time series of daily site coordinates at MENA is given in Figure 4 for the East component with respect to the first analysed day (DOY 2013:019). Note that the motion of the Australian continental plate has already been subtracted using a plate motion model [59]. All coordinate differences are within +/− 8 mm for the East component. From this time series, no effect of the CR deployment on 25 May 2016 (dashed line in Figure 4a) is evident. Similar to the East component, no visible change in coordinates is detected for the North and Up components. To further investigate potential changes to the coordinate quality before and after CR deployment at the site, we visualise the 2-sigma coordinate standard deviations resulting from differential GPS processing for each daily site coordinate (Figure 4b). Note that the absolute values of formal standard deviations resulting from processing of daily GPS data are overly optimistic which can be attributed to un-modelled errors sources and serial correlations present in GPS observations [60], [61]. However, the standard deviations can still be used to investigate relative changes in the uncertainty level. A clear seasonal effect is visible in the standard deviations, particularly for the Up component. But as for the coordinate components, no significant change is visible when comparing the periods before and after the CR deployment.

Figure 4

Daily coordinate solutions resulting from GPS processing with Bernese GNSS Software at site MENA: a) coordinate differences for all three coordinate components (East, North, Up) with respect to the first epoch on 19 January 2013 and b) standard deviation of coordinate estimates for all three components. The date when CR were attached to the monument is marked by a grey, dashed line.

In order to statistically assess potential effects of the CR deployment, we calculate the average coordinate standard deviation for the period before and after deployment, as well as a measure of coordinate variability defined by the absolute difference to a smoothed model of the daily coordinates. The smoothing model uses a generalised moving average spanning 100 data points and filter coefficients determined by an unweighted linear least-squares regression with a polynomial model of degree two [62]. The corresponding numbers for average coordinate variability and standard deviation are listed in Table 5 along with the difference between the two periods (after CR deployment minus before CR deployment). A slight increase is observed for both statistical measures (i. e. coordinate variability and standard deviation) and in all coordinate components. However, the difference is small with respect to the absolute numbers of both coordinate variability and standard deviation, and could also be explained by other site-specific effects in the surroundings of site MENA (e. g. vegetation growth). To rigorously test the statistical significance of the difference between the two periods, we perform hypothesis testing for both standard deviation and the coordinate variability samples. Since both samples are not normally distributed (i. e. non-gaussian) we use a nonparametric test for equality of population medians of two independent samples known both as the Wilcoxon Rank Sum test or the Mann-Whitney U-test [63], [64]. We test the null hypothesis that the data sets of the two periods are samples from continuous distributions with equal medians against the alternative that they are not. The Rank Sum test statistics reveal that the difference in standard deviations are significant for the East and Up components at a 99 % confidence level (with a p-value of 0.001 and 0.004, respectively, which is less than the chosen significance level of 1 %), but not significant for the North component (p-value of 0.1). However, for all three components the relative increase in the mean or median standard deviation values is 2 % of absolute value only (cf. Table 5). The difference in coordinate variability is not significant at a 99 % confidence level for all three components (East, North and Up) according to our hypothesis testing (with the p-value of 0.58, 0.18 and 0.66, respectively, being above the 1 % significance level). As a third measure to quantify the difference in daily site coordinates before and after CR deployment, we calculate the RMS level of coordinate values for the two periods by referencing the two data sets to the first coordinate in each period. The RMS values of the period after CR deployment are smaller for all three coordinate components (see Table 5). While the difference between the two periods is negligible for the East component, a larger decrease in RMS of 1 mm and 3 mm is observed for the North and Up component, respectively. We can therefore conclude that the impact of the CR deployment on daily site coordinates is below the level of variability related to other effects at a continuously operating GNSS site (e. g. hydrological or thermal variations).

Table 5

Statistics of daily site coordinates at site MENA before and after CR deployment.

Period	# of days	Coordinate variability [mm]			Mean standard deviation [mm]			Root mean square level [mm]
		East	North	Up	East	North	Up	East	North	Up
Before	1222	0.70	0.63	3.20	1.04	1.09	3.15	2.64	5.46	8.54
After	1222	0.74	0.67	3.23	1.06	1.11	3.21	2.55	4.51	5.33
Difference		0.04	0.04	0.03	0.02	0.02	0.06	−0.09	−0.95	−3.21

3.2 SNR values

The signal-to-noise ratio (SNR) is a measurement of the strength of GNSS signals received. SNR values are measured in the GNSS receiver, however different algorithms and definitions are used by the various manufacturers making it difficult to compare SNR values measured with different equipment (e. g. different quantisation levels [65], [66], [67]). It is therefore recommended to calculate effective SNR from the raw SNR values saved in GNSS RINEX data for signal interference detection [65], [68]. Many multipath studies make use of the assumption that SNR values follow a nominal, elevation-dependent curve mainly determined by the antenna gain pattern. Some manufacturers and research groups have published receiver antenna gain (or radiation) patterns for right-hand and left-hand circular polarisation (e. g. [69], [70]), which can be used for comparison when the specific antenna is operated in a real-world scenario. The signal strength pattern of the Trimble antennas used at the CAxx sites was measured during individual antenna calibrations revealing a circular-shaped, almost purely elevation-dependent pattern with maximum signal strength at high elevation angles (approximately 51 dB-Hz and 44 dB-Hz at 90 degrees elevation for L1 and L2, respectively) decreasing to smaller values for low elevation angles (approximately 39 dB-Hz and 21 dB-Hz at 0 degrees elevation for L1 and L2, respectively). For the other antenna/receiver combinations used in this study (see Table 3), the gain patterns are unknown. At site MENA, however, we can investigate the relative changes in SNR values before and after CR deployment by applying the sidereal lag. Sidereal filtering approaches work best using high-rate SNR data (e. g. 1 Hz). For this study, we only have access to SNR values recorded at 30 s intervals, which makes a more sophisticated SNR multipath analysis difficult and may result in aliasing of high-rate multipath effects.

To compare SNR values before and after CR deployment, we perform stacking of epoch-wise SNR values for seven consecutive days by applying a sidereal transformation (shift by the nominal lag of 3 min 56 s for each day). Note that the multiples of the sideral lag ( 236 s ∗ k, with k = 1 … 7) are rounded to the closest 30 s interval in order to match the sampling of the input data. For the same reason (i. e. the 30 s sampling interval of the input data), we neglect variations from the nominal lag for individual GPS satellites which can be up to 10 s [42], [57]. The comparison of stacked, epoch-wise SNR values for each GPS satellite for a seven-days period before and after CR deployment is shown in Figure 5 for the L1 and L2 frequencies. Note that GPS satellite G02 was not observed on all seven days of the second period. Also, a data gap at DOY 143 (first period, 1 h 39 min) results in no data for the first 200 epochs. The fit between the SNR curves for the two periods (before and after) is relatively good with an overall RMS difference value of 0.25 dB-Hz and 0.40 dB-Hz for L1 and L2, respectively (note that the difference of two values in the logarithmic dB scale equals the logarithm of the ratio of the two values in a linear scale). The mean/minimum/maximum SNR differences (after minus before) are 0.07/−1.64/2.25 dB-Hz and 0.06/−3.18/4.36 dB-Hz for L1 and L2, respectively, indicating a slight increase in SNR for the period after CR deployment. Hypothesis testing of the SNR before and after CR deployment (using the Wilcoxon Rank Sum test as outlined in Section 3.1) reveals that the difference of the two samples is not significant, neither for L1 nor for L2 SNR values. However, the low temporal sampling used for this investigation (30 s) inhibits detecting potential high-rate multipath effects present in the SNR data and might result in aliasing effects. A dedicated test setup using high-rate SNR observations (1 Hz) would be required to detect potential multipath effects acting at a shorter temporal sampling and to test the statistical significance of SNR differences more rigorously.

Figure 5

Comparison of stacked SNR values before and after CR deployment at site MENA. Stacking of seven days for both periods (before and after) has been performed by applying a sidereal transformation. a) SNR values for L1 frequency, b) SNR values for L2 frequency.

To further investigate potential changes to SNR before and after CR deployment, we have a closer look at differences in stacked SNR values for observations from elevation angles above 60 degrees. The rationale behind this investigation relates to multipath patterns observed at high elevation angles (above 60 degrees), particularly for the sites equipped with Trimble Zephyr antennas (see Section 3.3.2). Again, the seven-day period after CR deployment has a slightly higher SNR level for both L1 and L2. The mean differences considering SNR values for observations from elevation angles above 60 degrees are +0.03 dB-Hz and +0.07 dB-Hz for L1 and L2, respectively. Almost all differences for individual satellite passes are positive (see Figure 6). Particularly for L2 a relatively large increase in SNR of up to 0.15 dB-Hz is observed for the period after CR deployment. The increase in SNR for observations from high elevation angles might be a result of constructive interference (e. g. [71]) of signals reflected at the CR with direct signals.

Figure 6

Difference in stacked SNR values before and after CR deployment for observations from elevation angles above 60 degrees (mean value for individual GPS satellites). The differences are plotted with respect to the mean azimuth angle of the corresponding GPS satellite pass. The grey lines give the centre orientation of ascending (253 degrees) and descending (107 degrees) pass CR.

3.3 Multipath stacking maps

The results of multipath stacking of a full year of PPP L3 (Lc) residuals are presented in this Section separately for site MENA (Section 3.3.1) and the four CAxx sites (Section 3.3.2). In both subsections, other sites equipped with the same antenna are used for comparison and validation.

3.3.1 Site MENA

At site MENA we perform multipath stacking (as described in Section 2.2) separately for two one-year periods before and after CR deployment. The resulting stacking maps are shown in Figure 7 along with their difference (after minus before) as skyplots. Larger positive and negative residuals are observed at low elevation angles which can be related to higher noise in the original GPS observations (e. g. due to atmospheric path delays and interaction with vegetation). An oscillating pattern is visible in the West for elevation angles between 20 and 40 degrees for both stacking periods (Figure 7a and Figure 7b). We relate this to multipath effects caused by a large metal shed located in that direction, approximately five meters away from the GNSS antenna and exceeding the GNSS site in height. The difference in stacked residuals (Figure 7c) reveals changes of up to 53 mm. A pattern of large positive change is visible in the ESE direction for elevation angles around 10 degrees, and may be related to dynamic changes in the vegetation (trees and bushes growing in that area over the course of a one year time period).

Figure 7

Multipath stacking maps at site MENA for different one-year periods. The local azimuth is given in polar coordinates with North, East, South and West directions marked by solid lines. Concentric rings are shown for 0, 30 and 60 degree elevation angles, with a dotted line indicating the elevation cut-off of 7 degrees. a) before CR deployment (DOY 2015:146 to DOY 2016:145), b) after CR deployment (DOY 2016:147 to DOY 2017:145), c) difference between a) and b), i. e. “after” minus “before”.

A slight change in the multipath pattern is visible at higher elevation angles shown in a zoomed view of the stacked residuals before and after CR deployment (Figure 8). A pattern of oscillating residual values is visible after CR deployment at elevation angles above 60 degrees, aligned approximately in North-South direction (Figure 8b). This oscillating pattern is not visible before the CR deployment but is small in magnitude with maximum values of around ±5 mm. The location of the oscillating pattern in azimuth and elevation suggests that it is related to small multipath effects induced by the ascending and descending pass CR. Since similar patterns are observed at the CAxx sites at the same location, but at a greater magnitude (see Section 3.3.2), the pattern in Figure 8b can be addressed to the CR deployment at site MENA. However, other effects visible in the multipath stacking maps (e. g. effects related to particular satellite tracks, particularly in the period before CR deployment) make the interpretation of the small, CR-related multipath at site MENA difficult.

Figure 8

Zoomed view of the multipath stacking maps given in Figure 7 at site MENA for high elevation angles. a) before CR deployment, b) after CR deployment, c) difference between a) and b), i. e. “after” minus “before”.

To validate the effects visible in the stacking maps of site MENA, we perform residual stacking for the same periods as for MENA for a set of other sites in the greater Sydney area equipped with the same antenna (i. e. ASH701945E_M, SCIS, see Table 2). Similar as at site MENA, large positive and negative residuals are visible at low elevation angles for all of the six analysed validation sites. Some sites show oscillating patterns related to multipath induced by structures in the vicinity of a site (for example site SPWD, cf. Figure 3c). However, none of the six sites shows a similar pattern at elevation angles above 60 degrees as visible in Figure 8b at site MENA. Comparing the variation of stacked residuals and their change in the two periods, we find that site MENA has a small increase in both standard deviation and mean absolute deviation of stacked residuals for the second period compared to the first period, whereas the six validation sites rather show a slight decrease or no change in these statistical parameters (Table 6). However, parts of this increase might be related to changes in the surroundings of site MENA not related to the CR deployment (e. g. vegetation growth). Comparing the absolute numbers at the different sites, site MENA is still performing better than most of the other sites, independent of the analysed period (i. e. second smallest standard deviation and mean absolute deviation in both periods, see Table 6).

Table 6

Comparison of standard deviation (STD) and mean absolute deviation (MAD) of stacked residuals at site MENA with other sites in the Sydney area equipped with the same antenna. STD and MAD are given in millimetres for both one-year periods (before and after CR deployment at site MENA) as well as the change between the periods (in percent).

Site	1st period		2nd period		Change [%]
	STD	MAD	STD	MAD	STD	MAD
MENA	4.4	2.8	4.8	2.9	+9 %	+6 %
GLBN	5.1	3.6	4.9	3.4	−4 %	−6 %
MGRV	4.1	2.7	3.9	2.6	−5 %	−5 %
NWRA	6.7	4.6	6.7	4.5	−0 %	−1 %
SPWD	8.6	5.7	8.5	5.6	−1 %	−1 %
VLWD	5.6	3.5	5.4	3.4	−2 %	−4 %
WFAL	6.5	4.2	6.4	4.1	−1 %	−2 %
all others	6.5	4.3	6.1	4.0	−2 %	−3 %

3.3.2 CAxx sites

At the four CAxx sites, multipath stacking can only be performed for periods with the CR attached to the antenna pole since these sites were newly constructed in June 2016 and used for geodetic monitoring with GNSS and InSAR since installation. We select a one-year period that minimises data outages, particularly at site CA19, which had several battery-related outage periods during its lifetime. The resulting stacking maps for this period for site CA07, CA10, CA17 and CA19 are shown in Figure 9. Large residuals are visible at low elevation angles. While site CA07 reveals smaller magnitudes for the stacked residual values at low elevations, the other three sites show some larger effects (particularly site CA19), which can be related to the individual setting of the four sites with regards to obstacles in the surroundings (particularly trees and vegetation). At high elevation angles, a similar oscillating pattern is visible at all four sites. This pattern was already faintly visible at site MENA (Figure 8b) and can be related to multipath interaction with the ascending and descending pass CR. The magnitude at the CAxx sites is larger compared to site MENA, reaching values of around ±15 mm (see Figure 10).

Figure 9

Multipath stacking maps at CAxx sites for a one-year period (DOY 2017:121 to DOY 2018:120) post-installation. Same axis and colour scaling as in Figure 7 (colour range ±40 mm). a) site CA07, b) site CA10, c) site CA17, d) site CA19.

Figure 10

Zoomed view of the multipath stacking maps given in Figure 9 at the CAxx sites. Same colour scale as in Figure 8 (colour range ±20 mm). a) site CA07, b) site CA10, c) site CA17, d) site CA19.

As for site MENA, we validate the effects visible in the stacking maps of the CAxx sites with a set of other sites in the greater Sydney area equipped with the same Trimble Zephyr antenna (see Table 2). Again, large positive and negative residuals are visible at low elevation angles for all seven analysed validation sites as well as oscillating patterns indicating multipath effects at some sites (particularly at site HNSB). However, none of the seven validation sites shows a similar pattern at elevation angles above 60 degrees as visible in Figure 10 for the CAxx sites which confirms that this pattern is exclusive to the CR being attached to these sites. Comparing the variation of stacked residuals for the same data period (see Table 7), we find that the CAxx sites are performing better than most of the other sites which are mostly mounted on or located close to buildings. The average standard deviation and mean absolute deviation of the CAxx sites is 6.4 mm and 4.1 mm, respectively, which is significantly lower than the corresponding average values at the validation sites of 7.9 mm and 5.4 mm, respectively. The values in Table 6 and Table 7 are used to assess the long-term quality of the different sites. High-frequency multipath effects at a site may be averaged out and not be reflected in these average deviation estimates. The reasonability of these quality indicators is evidenced by checking and comparing the environment of the different sites (see [72] for images of most of the sites). For example, the site revealing the largest deviations (i. e. site Hornsby, HNSB) is located on a metal roof with metal poles and rails surrounding the antenna, while site Oberon, OBRN, is attached to a concrete wall at the edge of a building and is therefore showing a much smaller standard deviation compared to HNSB. When comparing the values in Table 6 and Table 7, it becomes obvious that the sites equipped with the higher quality choke ring antennas (i. e. the sites listed in Table 6) have smaller deviations on average compared to the sites equipped with the lower quality Trimble Zephyr antennas. One exception is site Springwood, SPWD, which is located next to a railway line and surrounded by metal fences, pillars and horizontal bars, some of them well above the local horizon of the antenna, inducing multipath at various azimuth and elevation locations (see Figure 3c).

Table 7

Comparison of standard deviation (STD) and mean absolute deviation (MAD) of stacked residuals at CAxx sites to stacked residuals of other sites in the Sydney area equipped with the same antenna. STD and MAD are given in millimetres.

Site	STD	MAD
CA07	5.6	3.8
CA10	5.9	3.8
CA17	6.4	4.0
CA19	7.9	4.7
all CAxx	6.4	4.1
BANK	9.6	6.9
FTDN	8.4	5.2
GERR	7.2	5.4
HNSB	10.1	7.2
KTMB	6.5	4.1
OBRN	6.4	4.3
PCTN	7.0	4.8
all others	7.9	5.4

4 Discussion

In Section 3, we presented the impact of CR deployment at GNSS sites based on detailed analyses of (i) daily coordinates, (ii) SNR and (iii) observation residuals. For (i) and (ii) we have investigated the impact at a site which has a continuous time series with and without CR attached to the GNSS monument. For (iii) we could make use of four additional sites which had CR attached to the monument since their installation. In the following, we will summarise and discuss our findings.

A slight increase in coordinate variation and uncertainty (< 0.1 mm) is visible when comparing daily site coordinate estimates for a period of 1222 days before and after CR deployment at site MENA (see Section 3.1). Hypothesis testing of these results indicates the change in coordinate variation is not significant for all three components (East, North and Up), while the change in standard deviation is slightly significant for the East and Up components (at a 1 % significance level). A comparison of RMS values for the same periods (before and after CR deployment) shows lower RMS for the period with CR deployed at the site, particularly for the East and Up components. This observation indicates that other site-specific effects at continuously operating GNSS sites (e. g. related to hydrology or vegetation growth) have a larger impact on daily coordinate estimates compared to the deployment of CR.

A comparison of SNR values before and after CR deployment by applying the time lag of the sidereal day to the 30s SNR observations was presented in Section 3.2. The change in SNR values at site MENA for a seven-day period before and after CR deployment is small with RMS values of 0.25 dB-Hz and 0.40 dB-Hz for L1 and L2, respectively, which is less than 1 % of the absolute SNR values. We find a slight increase of SNR values for the period after CR deployment for elevation angles above 60 degrees which could be related to constructive interference of GPS signals reflected at the CR with direct signals. The SNR results indicate that the CR produce more constructive reflected signals at high elevation angles compared to low elevation angles. This can be explained by the location of the CR with respect to the GNSS antenna and the geometry of incoming GPS signals. According to Snell’s law of reflection, most of the reflected signals from the CR resulting from low elevation angle incoming signals would miss the bottom of the antenna, while signals from high elevation angles would be able to get reflected onto the GNSS antenna. The SNR increase at high elevation angles is small however (maximum of 0.15 dB-Hz averaged for one satellite track), and may equally well be related to effects other than interaction with the CR.

Our major assessment of potential multipath induced by CR focussed on empirical site models resulting from stacking of observation residuals. This spatial visualisation of site-specific residuals enabled multipath from objects in the surroundings to be identified. We find that multipath patterns are visible at sites with and without CR. At site MENA, we compare the location and magnitude of site-specific effects before and after CR deployment and find that differences which can be attributed to the CR deployment are only visible at high elevation angles (> 60°), revealing a low-magnitude oscillating pattern along an East-West directed axis. This observation is in agreement with theoretical considerations on the alignment and orientation of the CR and the angles of incidence and reflection of GPS signals (Snell’s law, see below). Multipath patterns with a similar shape and spatial location are also visible at the CAxx sites, but larger in magnitude (approximately ±15 mm compared to ±5 mm at MENA). Compared to site MENA, the CAxx sites are equipped with lower-quality Trimble Zephyr antennas (no choke ring) and therefore more susceptible to receiving reflected signals.

To check the plausibility of the patterns observed in Figure 10, we perform a theoretical valuation of signal reflections using the concept of Fresnel ellipsoids for the specific GNSS antenna and CR geometry. Significant signal reflection only occurs for objects which are large enough in size to span a significant fraction of a cross-section of the first Fresnel zone. For example, objects less than a wavelength λ in size ( λ = 19 cm for the GPS L1 signal) are not significant sources of reflected energy. The cross-section of the first Fresnel zone is an ellipse located on the reflecting surface. For satellite observations from zenith directions (i. e. elevation angle E = 90 °), the cross-sectional area results in a circle [28]. For a non-horizontal reflector, the tilt angle of the reflecting surface needs to be considered. We apply the equations for reflections at arbitrarily oriented surfaces given by [29] to the specific case of our CR (tilted by an angle α = 16 ° against the horizontal, antenna at an orthogonal distance d of approximately 1.1 m to the reflecting surface, see also sketch in Figure 11) to calculate the size of the Fresnel zone:

The resulting elevation-dependent Fresnel zones for this case (for the GPS L1 frequency) are still almost circular for an elevation angle of 90° (semi major and minor axis of the ellipse of a = 50 cm and b = 48 cm, respectively), but increasing for smaller elevation angles (e. g. semi major and minor axis of the ellipse of a = 81 cm and b = 57 cm, respectively for an elevation angle of 60°). The CR plates are square shaped with a size of 60 cm × 60 cm. Considering the calculated Fresnel zone sizes, signal reflections from the CR are more likely to occur for high elevation angles and for areas close to the centre of the CR baseplate.

Apart from the size of the reflector, a second criterion relating to the roughness of the reflecting surface needs to be fulfilled for specular multipath to occur (in contrast to diffuse reflection which is less critical due to its random nature). If the reflecting surface is smooth compared to the signal wavelength, specular reflection will be the dominating type of reflection. Roughness can be described as a deviation from the mean height of the (tilted) surface. In the case of CR, the curvature of the plates and other surface irregularities need to be small to be useful radar targets (e. g. within 1 mm for X-band radar data [73]). For the CR used in this study, flatness of the baseplates resulting from the manufacturing process of the material is deemed to be less than a few millimetres across the face. For the case of electromagnetic signal reflections, the elevation-dependent Raleigh criterion can be used to decide whether specular or diffuse reflections will dominate (e. g. [29]). The smallest threshold Δ H m i n results from elevation angles of 90° yielding a value of 2.4 cm for the GPS L1 frequency:

The roughness of our CR is much smaller than Δ H m i n , meaning that mainly specular reflection will occur at the baseplate of the CR. Figure 11 displays a model of the vertical angles of multipath signals at the CR baseplates for the case of specular reflections at the centre of the baseplate. The sketch shows a vertical profile of the setting viewed from a location perpendicular to the CR orientation axis (i. e. at a cross section showing the baseplate diagonal with a length of 84.5 cm). Using the given CR dimensions, the CR tilt angle ( α = 16 °), the width of the CR mount (ca. 10 cm) and its vertical distance from the GNSS antenna (ca. 50 cm), we can calculate the length dimensions given in Figure 11 from trigonometry. The angle β between the CR baseplate and the reflected GPS signal can then be calculated for a signal reaching the GNSS antenna centre (red dashed line in Figure 11) from:

Figure 11

Sketch of GPS signal reflections for the specific CR scenario presented in this study. The view is from an angle perpendicular to the CR orientation axis.

For locations reaching the GNSS antenna at its borders (yellow and blue dashed lines in Figure 11), we add or subtract 17 cm which corresponds to the radius of the Trimble Zephyr antennas used at the CAxx sites (34.3 cm diameter according to specifications). To derive the elevation angle E of the incoming GPS signals, we make use of the identical incoming and outgoing angles of reflection θ relative to the normal of the CR baseplate:

For this scenario the resulting elevations are 78.7°, 84.9° and 92.3° (i. e. 87.7° on the other side of the zenith direction), for incoming signals at the far, centre and near locations of the GNSS antenna, respectively. The results of this basic model of signal reflections agree with the observations in Figure 10 with large absolute values of stacked residuals occurring mostly at elevation areas greater than 75°.

In Figure 12, we model the horizontal direction of reflected signals by looking at the geometry of reflected signals at the two CR at each site from a top-down view. We consider two cases, (i) a Fresnel zone with the same diameter as the CR (shown for the ascending pass CR on the left in Figure 12) and (ii) Fresnel zones with a smaller diameter (i. e. 50 cm) resulting in different centre locations of reflection on the CR baseplate (shown for the descending pass CR on the right in Figure 12). The azimuth angles of incoming GPS signals can be calculated from trigonometry using the known distances to the reflection points on the CR and the given CR orientations. For case (i), the azimuth of incoming GPS signals varies between 93° and 121° for a descending pass CR and between 239° and 268° for an ascending pass CR. In the case of smaller Fresnel zones with multiple reflection locations on the baseplate, the azimuths vary between 85° and 129° for the descending pass CR and between 231° and 276° for the ascending pass CR. The observed azimuth regions with large absolute values of stacked residuals in Figure 10 indicate that multipath is occurring for azimuths below 90° and above 125° as well as below 235° and above 270°. This means that multipath is likely occurring from reflection points offset from the centre of the baseplate and Fresnel zones smaller than 60 cm in diameter. The relatively small sizes of potential reflection zones compared to the theoretical numbers resulting from Eq. (1) might be explained with the flatness of the CR resulting in almost purely specular reflection of incoming GPS signal. Note that we have not considered the gain pattern of the GNSS antenna in these simple reflection models which may explain some of the irregularities visible in the observed patterns in Figure 10. However, the model does in general confirm the azimuth and elevation zones visible in the multipath stacking maps.

Figure 12

Top-down view of the CR scenario presented in this study modelling horizontal angles of incoming GPS signals. The ascending pass CR left of the GNSS antenna shows an example of one Fresnel zone with a 60 cm diameter (i. e. covering the entire CR baseplate). The descending pass CR right of the GNSS antenna shows an example of four Fresnel zones with diameters of 50 cm with centre-points off from the centre of the CR baseplate.

When comparing the impact of detected multipath patterns at the CAxx sites to effects at other sites in the region (see Table 7), we find that the CR-induced multipath is small compared to the effects at the other sites. In fact, a comparison of the overall magnitude and variation of stacked residuals reveals a smaller level of site-specific effects at the sites equipped with CR than at most of the other sites under investigation. In other words, the quality of the five sites equipped with CR is still better than the average quality of surrounding sites with the same antenna. This observation can be explained by the fact that most of the validation sites are located on top of a building or close to other metal structures which induce multipath effects to the GPS signals received at the antenna.

5 Conclusions

The potentially detrimental effect that the presence of CR co-located at GNSS sites has on the GNSS signal is an understandable concern for GNSS network operators. In this study we aimed to establish whether there is any validity to these concerns. Note that the results of this study are specific to our design of CR and to the location of the CR relative to the GNSS antenna. Different sizes, plate shapes or volumes of metal used for the CR will probably change the results. In our set-up, we quantified the effects for a situation where the CR are located very close to the GNSS antenna (i. e. two CR attached to the antenna pole with a vertical distance of approximately 50 cm). In our set-up we find that CR directly attached to a GNSS antenna pole:

1. have no significant effect on daily site coordinates: for the site tested in this study the effect on site coordinates is less than 0.1 mm;

2. may induce small multipath effects to GNSS phase observations: mainly for observations from elevation angles above 60 degrees and depending on the GNSS antenna type used;

3. induce smaller multipath effects compared to multipath observed at GNSS sites on or at a building.

Furthermore, our investigations have shown that multipath effects induced by CR are significantly smaller when the co-located GNSS site is equipped with a choke-ring antenna (approximately by a factor of three compared to a ground plane antenna without choke rings). However, a quality assessment for the specific equipment and scenario used at a co-location site is recommended. As an example, we have applied a geometric reflection model using Fresnel zones to describe the regions of impact for the specific CR setting confirming the observed multipath patterns at high elevation angles and for specific azimuth zones. The multipath maps resulting from empirical site models using residual stacking can be used in the future to investigate multipath effects for other CR/GNSS co-location scenarios. Quantifying site-specific effects using the presented residual stacking approach can also be used to generate azimuth- and elevation-dependent correction data sets for GNSS sites (similar to GNSS antenna correction files), which allow for the correction of GNSS signals at a site at the observation level, prior to the parameter estimation process.