Is GNSS real-time positioning a reliable option to validate erosion studies at olive grove environments?

Aim of study: Soil degradation in agricultural areas is a widespread problem. In this framework, a data validation methodology is presented, including a study of the spatial resolution of Global Navigation Satellite System (GNSS) measurements, the calculation of erosion/deposition models, and the contribution of dual frequency and low-cost single frequency GNSS receivers.Area of study: A test olive grove in SE Spain.Material and methods: The study is based on three observation campaigns, between 2016 and 2018, using different GNSS receivers and working modes. The comparison between different surveys provide the volumetric variation over the analyzed period.Main results: Considering the dual-frequency receiver, there was no statistically significant difference between the means and the variances from 1.5 m and from 4.5 m data resolution at the 0.05 significance level. In order to estimate vertical differences from successive GNSS campaigns a differential digital elevation approach was applied. Although the differences depended on the zone of the test area and they changed along the monitoring period, the erosion rate could be catalogued as very low. The dual-frequency receiver satisfied the vertical centimetric precision limits for high accurate Digital Elevation Model (DEM), making it a reliable and accurate option to validate erosion studies in small areas.Research highlights: The results have allowed the characterization of multi-annual spatial redistribution of the topsoil at local scale, being of great help to design future prevention actions for the “tillage erosion” in olive grove environments. However, more tests are needed to guarantee the feasibility of low-cost receivers.


Introduction
Soil erosion is a problem that produces an important impact on the landscape, especially in olive (Olea europaea L.) grove environment (Vanwalleghem et al., 2011;Gómez et al., 2014). Sloping olive orchards are widespread, occupying large parts of SE Spain landscape where soil erosion is one of the most important environmental risks. The national government estimates erosion rates of above 50 t ha-1 yr-1 on olive orchards located in mountainous areas in Andalusia, which illustrates the severity of the environmental degradation risk (Taguas & Gómez, 2015). Here, the surface of olive groves with low erosion (0 to 12 t ha -1 yr -1 ) represent 47.2% of the total area of the Andalusian olive grove. Surfaces with moderate erosion (12-50 t ha -1 yr -1 ), high (50-100 t ha -1 yr -1 ) and very high (more than 100 t ha -1 yr -1 ), represent respectively 29.7%, 11.8% and 11.2%. In absolute values, the province of Jaén is the largest area of olive groves with high or very high erosion (310,258 ha) followed by Córdoba (159,818 ha) (Junta de Andalucía, 2015). Soil erosion is a highly variable process in space and time and requires precise and accurate spatial representation of surface associated with each time epoch (Ramos et al., 2008). The accuracy of altitudes affects the quality of digital elevation models (DEMs), so this requirement must be added to high resolution DEMs. The validation of soil erosion modelling using field data is usually limited (Žížala et al., 2019) and cost-effective methods are necessary. The current terrestrial and airborne laser scanning collect elevation data to generate a 3D surface model, however, they involve high cost so they are only suitable for large area and not for small project with low budget (Tahar et al., 2013). Others sources of spatial information, as DEM from UAS (Unmanned Aerial System) images (Stöcker et al., 2015;Pineux et al., 2017) and geodetic positioning methods, are also being widely used (Abd Aziz et al., 2012;Garrido et al., 2013). All of them offer different results in terms of accuracy, resolution and cost. Real-Time Global Navigation Satellite System (GNSS) positioning using high precision receivers is currently a commonly approach used in precision agriculture (Álamo et al., 2012;Guo et al., 2018), however, the use of low-cost receivers has gained special interest in recent years (Kabir et al., 2016;Keskin et al., 2017) due to their multiple advantages: small size, ease of use, low-cost, and high precision even in real-time. Recent studies show that these receivers can achieve competitive positioning performance to survey-grade receivers (Dabove & Manzino, 2014;Odolinski & Teunissen, 2016Garrido et al., 2019a). The use of Real Time Kinematic (RTK) positioning has remarkable benefits when capturing elevation data at local scale in agricultural fields. In addition to low-cost post-processing efforts, local elevations and depressions are well represented, the horizontally and vertically accuracy is adequate for the representation of local terrain features and GNSS surveys may be applied with a temporal resolution adapted to the needs of a specific application.
In order to examine the reliability of RTK positioning to validate and/or to complement erosion studies at olive grove environments, three GNSS surveys, between 2016 and 2018, using dual-frequency and low-cost GNSS receivers and applying different RTK working modes, have been performed at a test olive grove in SE Spain. This study seeks to achieve the following objectives: first, to analyze the resolution of spatial data surveyed with GNSS receivers in order to optimize, in terms of reliability and cost, the data collection. Secondly, to evaluate the erosion or deposition rates in the test area throughout the analyzed period and finally, to compare the dual frequency and low cost single frequency GNSS receivers in terms of suitability for local erosion/deposition studies.

Real time GNSS positioning: RTK and NRTK solutions
In real time differential single-base GNSS positioning, the rover calculates its position relative to a reference station. The main limitation using single-base RTK solution is that some errors (orbital errors and ionospheric and tropospheric refraction) are decorrelated with distance (Hofmann-Wellenhof et al., 2008). The RTK positioning requires solving integer carrier phase ambiguities. This task may be relatively easy over short baselines (few kilometers); however, it becomes increasingly difficult as the baseline increases. The distance-dependent errors can be modeled using GNSS observations from several reference stations around the rover position, extending the RTK solution to a Network-based RTK positioning (NRTK). Different approaches coexist to generate the NRTK corrections (Takac & Zelzer, 2008), including the MAC (Master Auxiliary Concept), VRS (Virtual Reference Station), PRS (Pseudo-reference Station) and FKP (Flächen Korrektur Parameter) approaches. Currently, the first two approaches are the most used. In the MAC approach, the network RTK solution is calculated considering a main reference station (master) and several auxiliary reference stations (sub-network). This approach provides ambiguity-level observation data as correction differences of dispersive (ionospheric delay) and non-dispersive data (tropospheric delay and orbit errors) for each satellite-receiver pair. The main challenge of the network processing software is to reduce the ambiguities for the phase ranges from reference stations in the sub-network to a common ambiguity level (Euler et al., 2001). In the VRS approach, the rover sends a navigation solution to the network control center, which accepts this current position as the location of a "virtual" reference station, calculates the corrections relating to the VRS, and transmits them to the rover (Landau et al., 2002). From the user perspective, both approaches deliver positioning results at the same accuracy-level and are supported by the main GNSS receivers, however, Janssen (2009) indicated that the MAC approach appears to be superior, since it supplies raw correction data related to real reference stations, rather than modeled data. Today, RTK networks are an indispensable complement to GNSS positioning. They provide 3D positions in real-time with high precision and accuracy and define a common reference frame for different users and applications. In multi-temporal studies, as erosion processes, it is an essential requirement to have a stable reference frame. The RTK solutions used in this study considers as geodetic reference frame the Andalusian Positioning Network (RAP, http://www. ideandalucia.es/portal/web/portal-posicionamiento/rap).
This active network is a local geodetic infrastructure in the Andalusian Community (South of Spain) managed by the Institute of Statistics and Cartography of Andalusia (Páez et al., 2017). More than twenty permanent reference stations are distributed over the region of Andalusia with a maximum inter-stations distance of 70 km. RAP provides RINEX (Receiver Independent EXchange) observation files for post-processing solutions and real time differential correction services. This active network guarantees that the RTK solutions for each time epoch are in a compatible realization of the official reference system ETRS89 (European Terrestrial Reference System 1989).

Field site
The study region is located within the Northeast part of the Autonomous Community of Andalusia (S Spain), in the Province of Jaén (37º57´3´´N, 3º 58´54´´W) (Figs. 1A and 1B). The olive grove selected is representative of the majority of the olive plots in the South of Spain. It has a typically Mediterranean climate with a continental character and a mean annual rainfall of 500 mm, most of which falls between October and April. The test olive grove selected is 0.5 ha, unirrigated and the slope gradient varies between 2% and 25%. The relief is undulating with slopes with a decreasing elevation from north to south ( Fig. 1C and 1D). According the World Reference Base for Soil Resources (WRB) Soil Classification System (IUSS Working Group WRB, 2015), the soil is qualified as vertisol. The olive orchard has trees of similar characteristics -three trunks and about 40 years old. The mean diameter of the canopy of the olive trees is 5.5 m and tree spacing is 10 m × 10 m. In this olive grove, traditional tillage is the soil management system applied. Thus, sediment transport is significant due to "tillage erosion", moving soil from the top of the field downward where the test olive grove is selected. The test area includes a geodetic pillar denominated OLIV, which is composed by a concrete pillar anchored to the ground that incorporates an embedded reinforced centering system (Figs. 1D and 1E).

GNSS receivers: types and configurations
Two different GNSS receivers, configurations and working modes were applied in this study: -Geodetic dual-frequency GNSS receiver and antenna, to receive network-based RTK solution provided by RAP network, via Internet Protocol (NTRIP). The GNSS geodetic equipment used in 2016, 2017 and 2018 observation campaigns include a dual-frequency Leica GS10 receiver (equipment number 3660839), a triple frequency (GPS/GLONASS/Galileo) Leica AS10 antenna, and a radio field controller CS10 (Leica Geosystems AG. Heerbrugg, Switzerland). This receiver is equipped with a Siemens MC75 GSM/GPRS module and it is configured to receive NRTK MAC corrections from the Andalusian active positioning network (RAP) in RTCM format through a NTRIP connection. The RTK setting were 1-sec observation rate, maximum position dilution of precision (PDOP) = 6 and measurements number = 1, horizontal coordinate quality (CQ2D) ≤ 30 mm and vertical coordinate quality (CQ1D) ≤ 50 mm. The standard deviations indicated by the manufacturer specifications considering the network RTK solution are σxy = ± (8 mm + 0.5 ppm) for horizontal components and σh = ± (15 mm + 0.5 ppm) for vertical component. The instrumentation includes a telescopic carbon survey pole (locks at 2.00 m) with a circular bubble.
-Low-cost single-frequency GNSS receivers and patch antennas, considering a single-base RTK solution. A lowcost single-frequency GNSS u-blox NEO-M8P GNSS RTK module, type number C94-M8P-3-11 (u-blox AG. Thalwil, Switzerland) was introduced in the 2018 survey campaign. U-blox M8P module can receive GPS L1 C/A, GLONASS L1OF, and BeiDou B1 observables. To get an optimal positioning accuracy, the accuracy of the non-permanent reference station position, located at geodetic pillar OLIV, must be guaranteed. For that, previously we carried out static GNSS measurements at OLIV reference station using a dual-frequency GNSS receiver. The GNSS measurements were post-processed relative to the nearest RAP reference stations to obtain its coordinates in the same reference system and frame than the RAP network (ETRS89). The base station module is placed on the upper end of a 2 m pole on the pillar in order to reduce obstructions due to the canopy of the olive trees ( Fig. 2A). The rover module receives RTK corrections from temporal base station using RTCM protocol via a radio communication link enabling the rover to output its position relative to the base station down to centimeter-level precision. The manufacturer specifies the circular probable error for RTK positioning as ± 25 mm + 1 ppm (ppm limited to baselines up to 10 km). In addition, the patch antennas must be coupled to a metallic ground plane in order to mitigate the multipath effect (Fig. 2B). The u-center V8.29 GNSS evaluation software is used for the configuration of u-blox GNSS receivers, collecting the data and the performance analysis.

Data acquisition
The detection of slight changes at terrain morphology requires precise spatial information that must be acquired by repeated surveys. The comparison between different GNSS surveys performed at different years provide the volumetric variation over the analyzed period. The applied measurement technique considers a set of parallel altimetric profile lines established on the test olive grove. The profiles surveyed are oriented mostly Northeast-Southwest with an approximate azimuth of 10 degrees according to the maximum slope. According to previous studies in this test olive grove (Ramos et al., 2008;Garrido et al., 2019b), the distance between adjacent profiles is stablished in 1.5 m. This interval was maintained between consecutive points surveyed in the same profile. The reference GNSS survey was performed on June 2016. In 2017, the profiles were surveyed so that the start and end points of each profile were coincident with the first GNSS campaign. The RAP NRTK solution (MAC approach) was considered. The framework was established by the RAP network, assuring the same reference frame for all GNSS surveys. The collected data were referred to the ETRS89 -UTM projection (Zone 30) considering the ellipsoidal height (h).
On September 2018, the test olive grove was surveyed using two GNSS receivers, adding a low-cost GNSS unit to survey-grade receiver (Fig. 2B) this GNSS survey was to determine quantitatively how the low-cost single-base GNSS altimetric solution differs from the NRTK altimetric solution. Altimetric profiles were surveyed simultaneously with a dual frequency GNSS receiver, applying NRTK corrections from RAP active network, and with a low-cost single frequency one, receiving RTK corrections from a temporal reference station installed at OLIV geodetic pillar. In order to reduce the time and the cost of the GNSS survey, the distance between adjacent profiles was fixed at 5 m. It was also the approximate distance between consecutive points in a profile. Details of the different GNSS surveys can be found in Table 1. Figure 3 shows a detailed methodological scheme followed in the analysis process.

Analysis of data resolution
Fine DEMs derived from high-resolution spatial data will generate accurate estimations of soil erosion at local scale. However, researchers cannot always acquire such information due to the expensive cost. The optimal resolution must achieve a balance between the resource  optimization, the survey accuracy, the data processing and the need for interpretable outputs. The resolution of spatial data has been analyzed comparing two samples of elevation from dual frequency GNSS receivers. Spatial data with a mean resolution of 1.5 m were compared with a sample considering higher sparse resolution (Fig.  3). Summary statistics of elevation data sets for 2016 and 2017 surveys are given in Table 2. For each GNSS observation campaign, two independent sets of elevation data were compared through "F-test" (comparison of variances) and "t-test" (comparison of means) (Teunissen, 2000) to determine if the differences were significant (0.05 significance level) (Table 3). These tests have been used by Kumhálová et al. (2011) to compare DEMs and models of slope derived from topographic data. The results of the "F-test" showed that there was no statistically significant difference between the variances of both data sets at the 0.05 significance level, even taking into account the large difference in number of measurements of the two data sets. In addition, the results of "t-test" showed that there was no statistically significant difference between the means of the data derived from 1.5 m data resolution and from 4.5 m data resolution. These results suggest that both data sets described elevation similarly. This is true for both survey campaigns.

Analysis of erosion and deposition
In order to analyze the evolution of the erosion process at local scale in the test olive grove, for each GNSS survey an interpolation process was applied. The initial high density of points per hectare allowed obtaining high-  resolution DEMs (Fig. 3). As shown in Liu et al. (2011), DEM accuracy is very sensitive to horizontal resolution, vertical precision and density of sample points. DEMs with a grid resolution of 1 m × 1 m were generated from 2016 and 2017 survey campaigns considering original elevation data. Although different interpolation methods were available (Guo et al., 2010), the algorithm used in this study was the triangulation. This method is most suited to spatial data distributed over the area and produces a regular gridded raster from a set of input data points by creating a Delaunay triangulation. Additionally, for each survey, a model is generated considering the standard deviations of the ellipsoidal height (σh) as well as the corresponding histogram (Fig. 4). The spatial data were treated using MapInfo Pro Advanced V.16 software, which provides an ideal environment for geo-referencing, visualization, DEMs generation, subtraction and interpretation. After DEMs were obtained, a first differential model between 2016 and 2017 DEMs was calculated in order to quantify the annual processes of  erosion/deposition. High-resolution differential digital elevation model (DDEM) was created subtracting the elevation grid values from 2017 and 2016 with MapInfo Pro Advanced V.16 using a SQL mathematical sentence. In order not to overestimate the results, it is very important to remove those cells with null values in either of the two models considered. A quantitative approach was applied in order to estimate vertical differences between DEMs derived from successive GNSS survey campaigns. This methodology allows to calculate vertical differences and also identify loss and accumulation zones inside the test olive grove. The grid of the elevation differences (Fig.  5A) presents negative and positive values, but low in absolute values. That is, some zones experiment slight loss of soil (erosion) while others have sediment deposition (accretion). Besides, according to Ramos et al. (2008) and Wheaton et al. (2010), individual height uncertainties in DEMs were propagated in a DEM difference as, where σ hiDif is the propagated standard deviation in the differenced DEM, and σhi2016 and σhi2017 are the individual standard deviation at each point in the two DEMs considered. Although the observed elevation differences reached maximum values in the order of one decimeter in the vertical component, larger than the accuracy limit of GNSS positioning, only 60% of the differences calculated were greater than σ hiDif (Fig. 5B) and therefore could be considered significant differences. To this, we should add that the mean value for the elevation differences (-0.005 m) was very close to zero. In addition, we computed the erosion and deposition volumes between a later DEM (2017) and an earlier DEM dataset (2016). We highlight that the lost soil volume (96.96 m 3 ) was very similar to the deposition volume (90.65 m 3 ) estimated on the test area. Considering that the timespan between two GNSS surveys was one year, the erosion rate can be catalogued as very low. Nevertheless, the observed differences depended on the zone of the test area and they may change along the monitoring period.
Once the use of spatial data with lower resolution was justified (see previous section), and based on models from successive survey campaigns, a new quantitative approach was applied. To generate DDEM between 2016, 2017 and 2018 surveys, the spatial data originally at 1.5 m were resampled at 4.5 m in order to smooth the bias that eventually DEMs could have using different cell size. The DDEMs between campaigns could be calculated based on the same resolution level. If a resolution discrepancy exists between the two periods, the higher resolution scale should be converted to the lower one. After that, DEMs with a grid resolution of 3 m × 3 m were used in order to quantify the inter-annual processes of erosion/deposition on the test area (Fig. 3). The results derived from differences between DEMs corresponding to 2016, 2017 and 2018 surveys are shown in Table 4. We highlight that vertical differences in the range of a few decimeters were derived. Standard deviations for elevation differences extracted from 2017-2016 and 2018-2017 DEMs were very similar. The erosion and deposition volumes between 2016 and 2018 are also presented. The erosion and deposition volumes between 2016 and 2017 were according to the estimated values with high-resolution DEMs on the test area. The estimated values, especially the erosion volume, increased in the second period analyzed (Jun'2017-Sept'2018. Such differences may be due in part to the fact that the tillage is done during the summer season. It should be noted that the last period was slightly higher than the first one, including a greater number of months of tillage. However, new survey campaigns will allow a better analysis of this fact. The results obtained have allowed the characterization of multi-annual spatial redistribution of the topsoil at local scale in the test olive grove. This information is great help to design future prevention actions for the so-called "tillage erosion" in olive grove environments. However, further investigation is necessary to analyze the accordance of these results with erosion/deposition values derived from other techniques of massive acquisition of spatial data (for example, LIDAR system or UAS) and the compliance of the short time erosion/deposition rates with long-term rates on this environment.

Analysis of GNSS receivers
Complementing the above analysis and in order to reduce the cost of acquiring accurate spatial data, altitudes derived from two different GNSS receivers applied in 2018 GNSS survey were analyzed (Fig. 3). The objective was to assess the performance of a low-cost single-frequency receiver (u-blox NEO-M8P) using a single-base RTK solution and contrasts it with a geodetic dual-frequency one (Leica GS10 receiver with AS10 antenna) considering network-based RTK solution.
Two sets of elevation data, from geodetic dual frequency GNSS receiver and from low-cost single-frequency receiver, were available with a similar mean standard deviation of the ellipsoidal height (0.015m). In this case, the size of data sets was very similar: 244 from dual frequency and 236 from single-frequency receiver. The difference in the size of data sets was given by specific problems in the resolution of phase ambiguities with the low-cost receiver. Again, two data sets were compared applying the "F-test" and "t-test" in order to determine whether they were significantly different at the 0.05 significance level (Table 5). The "F-test" (0.9822) and corresponding p-value (0.5554) showed that there was no statistically significant difference between the variances of both data sets. The results of "t-test" (0.0644) and p-value (0.9486) also showed that there was no statistically significant difference between the means of the data sets at the 0.05 significance level. Therefore, after this preliminary analysis, it could be inferred that the low cost GNSS receiver could be used in order to reduce costs without affecting significantly the altimetric results. However, if we compare DEMs derived from both sets of elevation data, even though the mean value for vertical differences was equal to 2 cm, the standard deviation doubled this value (0.039 m), being these differences on the interval from -0.090 m to 0.153 m. Thus, we cannot conclude that the differences between both DEMs were not significant, especially when the spatial data corresponded exactly to the same period.
The results indicate that the geodetic dual-frequency receiver satisfies the vertical centimetric precision limits for high accurate DEM, making it a feasible and accurate option to validate erosion studies in small areas.
However, this study also shows that, even though the single-frequency u-blox NEO-M8P receiver combined with a low-cost patch antenna can achieve centimeter-level precision in real-time, the comparison of DEMs derived from geodetic and low-cost GNSS receivers shows significant differences. A disadvantage during the surveying campaign is that two low-cost GNSS receivers (reference and rover) must have a constant radio link. This may cause problems in the process of ambiguity resolution in an olive grove environment. In addition, the stability of the GNSS receiver used as reference station may affect the accuracy of altimetric results. More tests are needed to guarantee the quality of the altitude values obtained with a low-cost GNSS receiver in this environment. The use of low cost receivers in this type of applications should be questioned based on the preliminary results obtained, at least with single frequency low-cost receivers. It is necessary to advance in this line of research in order to guarantee the quality of the altimetry obtained with low cost GNSS receivers in an olive grove environment, for example, improving the stability of the GNSS receiver used as reference station and analyzing the results obtained with the new low-cost dual frequency receivers. However, this research demonstrates that GNSS geodetic tools have a great potential for advancing knowledge on on-site soil processes and could efficiently support current research related to soil redistribution and soil conservation strategies at olive grove environment.