Monitoring of Water and Tillage Soil Erosion in Agricultural Basins, a Comparison of Measurements Acquired by Differential Interferometric Analysis of Sentinel TopSAR Images and a Terrestrial LIDAR System

Abstract

Agricultural production, the main pillar of food security, is highly dependent on soil quality, and threatened by erosion processes that degrade soil quality. This article is part of a research to verify the usefulness of differential interferometric analysis on TopSAR (Terrain Observation with Progressive Scans SAR, Synthetic Aperture Radar) images to measure water and tillage erosion in small agricultural basins. For this, images from the Sentinel 1 mission are used, analyzing the deformations on the earth’s surface. The purpose of this research is to verify the accuracy of the proposed method by comparing its measures with the ones taken with the gold standard laser terrestrial LIDAR (Light Detection and Ranging) system, as well as to establish a basic step period framework that guarantees an admissible loss of coherence. The results on a pilot plot in El Molar (north of Madrid, Spain) showed that the differences lay within the range of the error associated with the very LIDAR system and showed that coherence losses correspond with the deformations measured. Given the economic and labor advantages of the differential interferometric analysis, this method could be regarded as an excellent alternative to the use of LIDAR in large-scale studies for measuring ground deformation caused by water and tillage erosion.

1. Introduction

Food security depends upon the sustainability of agricultural production, which depends on several factors involving soil quality, water availability, climate, pests, etc. Soil erosion is one of the main problems concerning soil quality and, therefore, sustainable agriculture and food security throughout the world [1,2]. Water and tillage erosion, and the consequent loss of soil, affects much of Europe’s active agricultural land [3] as well as agricultural land now in disuse [4,5,6,7]. Water erosion is generally observed at lower concave slope positions following the drainage network, while tillage erosion tends to be found at higher convex slope positions following the direction of tillage. The rate of tillage erosion may be much greater than water erosion [8,9].

Several models have been used for erosion prediction, such as the Universal Soil Loss Equation USLE [10], the Modified Universal Soil Loss Equation MUSLE [11], the Revised Universal Soil Loss Equation RUSLE [12], and, more recently, models based upon the solution of the continuity equation [13,14]. Yet all of them are based on predictions and do not involve real measurements except for calibration.

Historically, there have been four fundamental ways of erosion measurement [15]: Change in weight (CW), sediment collection from erosion plots and watershed (SC), change in channel cross section (XS), and change in surface elevation (SE).

The first (CW) and second (SC) are direct measurements that do not allow us to know the distribution of the processes inside the area monitored with the control point where measurements are taken. The third (XS) is limited to the river channels, while the last (SE) is the one that really allows for the study of the distribution of erosion/sedimentation processes.

The tools traditionally used to measure erosion belong to the SC class and only allow for non-localized quantification. For example, the erosion of basins has long been measured in terms of the total material lost [15], but the actual places from where this material came were never mapped, and processes of erosion-sedimentation inside the basin were never considered.

More recently, mathematical models that make quantitative assessments of local erosion [16,17] based on the availability of geographical information, mainly digital terrain modelling, vegetation cover maps, and rainfall data, have become more important [18]. However, to provide reliable results, these methods need calibration, requiring field experiments and measurements to be taken on-site [19,20,21]. A system able to remotely monitor local erosion without such drawbacks is clearly needed. Monitoring and predicting erosion from remote sensing data acquired from airborne full-waveform lidar systems [22] comes to offer a solution, but it is expensive and has limited coverage and low frequency.

Finally, erosion study by using differential interferometry and remote sensing to measure changes in earth surface elevation [23] has recently come as a new solution that takes advantage of the freely available, frequently produced, high resolution, and wide coverage synthetic aperture radar (SAR) images provided by the Sentinel-1 satellite system (part of the European Space Agency’s [ESA] Copernicus Programme, in collaboration with the European Commission). For descriptions of the Sentinel-1 system see Peter et al. (2017) [24] and Berger and Aschbacher (2012) [25].

Differential interferometry has been successfully used to locally analyze the effect on land surfaces of subsidence [26], volcanic eruptions [27], coal mining [28], earthquakes [29,30], slope movements [31], surface mass movements, and landslips [32,33]. The effects of engineering constructions [34,35] and fracking [36] have also been analysed, and even those of a jökulhlaup, a type of glacial outbreak flood, have been so-examined [37].

Interferometry is based on making precise measurements on SAR images of the intensity and phase of electromagnetic radiation bouncing back from the surface of the Earth to a radar device equipped with a synthetic aperture. The comparison of SAR images of the same area allows the generation of an interferogram in which the phase difference information is strongly related to the topography of the ground [38]. Deformations can then be represented in the form of maps [39]. This method permits erosion to be monitored in a spatially localized manner, allowing an understanding of the processes involved. This could be used to measure erosion-sedimentation processes at the scale of just a few square meters [23], but it is still necessary to verify such aspects as:

• The accuracy of the measurements
• Coherence loss boundary
• The effect of soil moisture
• The effect of soil decompaction
• Type of vegetation

The present work tries to clarify the two first issues, the accuracy of measurements and coherence loss boundary, by using an experimental plot at El Molar (Madrid, Spain), comparing the results obtained with differential interferometry with those obtained using a terrestrial laser-based LIDAR system, a current gold standard that provides results widely recognised to be a true reflection of ground erosion [40,41,42].

The deformation of the earth surface of the study plot between different dates has been measured with both the reference system (LIDAR) and the assessed (InSAR) system, and the maximum possible error of the reference system has been adopted as the maximum permissible difference to validate the assessed system.

2. Materials and Methods

2.1. Study Area

The study area was an experimental plot in El Molar, just north of the city of Madrid (Figure 1) with centroid coordinates EPSG25830: 450548.37, 4510027.21.

The study area belongs to an urban plot of approximately 101 m × 62 m located on the northern outskirts of the village of El Molar. The soils in the plots are the typical haploxerepts (inceptisols) of this area, moderately deep and well-drained. It has never been built, by the building of high coherence [43] in its east, west, and south sides, while on the north side, there is a natural dehesa. It is also surrounded on its south and west sides by traffic lanes with parking slots. The north part of the plot had some vegetation (including a couple of trees); meanwhile, the south part of the plot was cleared of vegetation just before beginning the fieldwork of the present study and several erosion-favoring vehicle tracks were present.

Only the south part of the plot was used for the study, and it was chosen because it was:

• Accessible and monitorable: These enabled checks to be made during satellite overflights for the presence of walkers, livestock, or machines that might figure as variations in the ground surface.
• Clear: The scant and small vegetation in the plot would not affect the differential interferometry analysis of the SAR images [44]. However, it would affect the measurements made by the LIDAR system.
• Size: Large enough to contain several interferometry result raster cells (13.93 m side), but not so much to make the daily terrestrial LIDAR work unaffordable.
• An area with surrounding plots suitable for comparison: Plots used in this kind of study need to have a high coherence contrast with bordering plots [44,45]. This allows their localization in coherence maps produced by differential interferometry. Since the present study required a plot with almost naked, erodible ground (low coherence), it was desirable that it be surrounded by plots with high coherence (e.g., asphalted or built-up areas).
• A nearby control plot was available: A nearby high-coherence control plot is desirable since it can provide the reference for zero deformation [37].

The plot was monitored at the exact moments when a Sentinel-1 satellite was overhead (approximately 07:15 h and 19:05 h every six days). A nearby road education court (which was always empty during satellite passes) within the school’s facilities at the east was the preferred zero deformation control point.

2.2. The Sentinel-1 Mission

The Sentinel-1 mission involves two satellites equipped with SAR devices that work in the microwave C band, which allows images to be captured both day and night and irrespective of cloud conditions [24,46]. Their polar orbits with a phase shift of 180° in the same plane allow each of them to cover every point on Earth at least once every 12 days; each point is crossed by one of the two satellites every six days.

The sensors have four acquisition modes [46]: stripmap (SM), interferometric wide swath (IW), extra wide swath (EWS), and wave (W). The most suitable mode available for this kind of work is IW [23].

2.3. Data Acquisition

Data were collected using the methods below.

2.3.1. TopSAR Images for Differential Interferometric Analysis

TopSAR images are able to be downloaded from the Open Access Hub of the ESA (https://scihub.copernicus.eu/ (accessed on the sames dates that the acquisitions were taken)) just a few hours after they are obtained. It is important to select images with the same ascending/descending orbit and with the maximum overlap area possible, to minimise the quality loss during coregistration.

The images used were captured by either the Sentinel 1A or 1B satellite, but always in IW mode and single-look captures. The preferred polarisation was VV to minimises the effect of small vegetation [47,48].

Table 1. The TopSAR images used in the present work.

Date	Satellite Sentinel	Orbit	Local Time
28 April 2017	S-1A	Ascending	19:10–19:12
10 May 2017	S-1A	Ascending	19:10–19:12
22 May 2017	S-1A	Ascending	19:10–19:12
3 June 2017	S-1A	Ascending	19:10–19:12
15 June 2017	S-1A	Ascending	19:10–19:12
27 June 2017	S-1A	Ascending	19:10–19:12
9 July 2017	S-1A	Ascending	19:10–19:12
15 July 2017	S-1B	Ascending	19:10–19:11
30 December 2017	S-1B	Ascending	19:10–19:11
11 January 2018	S-1B	Ascending	19:10–19:11

records the characteristics of the TopSAR images used in the present work.

2.3.2. In Situ Surface Measurement Using the Terrestrial LIDAR System

Every day when one of the Sentel-1 satellites passed overhead, the study plot was modelled three-dimensionally by terrestrial LIDAR (model FARO 3D), taking data from two measuring positions (see Figure 2).

Four to six spheres located strategically in the plot were used to act as reference points during data acquisition and registering [49,50], although additional reference field points such as buildings and streetlights were used whenever it was necessary [51]. The same spheres later served to help compare the models for different dates and thus recreate the changes that occurred over time.

2.4. Analytical Methods. Differential Interferometry

Interferometry with SARs allows the accurate measurement of the intensity and the phase of electromagnetic radiation during its journey, as recorded in a SAR image (Figure 3). A SAR image is composed of a two-dimensional mosaic of elements termed pixels, each one referring to a small area of the ground, or cell. Each pixel contains the intensity and phase information of the reflected radiation in its cell.

$$\mathsf{\varphi }=\frac{2\mathsf{\pi }\mathrm{R}}{\mathsf{\lambda }}2\mathrm{R}=\frac{4\mathsf{\pi }}{\mathsf{\lambda }}\mathrm{R}$$

(1)

SAR interferometry computes the phase difference between two SAR observations of the same ground unit; these observations are taken from slightly different positions to allow the distance to the ground to be determined [52] (Figure 4).

$${\mathsf{\varphi }}_1=\frac{4\mathsf{\pi }}{\mathsf{\lambda }}\mathrm{R}$$

(2)

$${\mathsf{\varphi }}_2=\frac{4\mathsf{\pi }}{\mathsf{\lambda }}\left(\mathrm{R}+\mathrm{\Delta }\mathrm{R}\right)$$

(3)

$$\mathrm{\Delta }\mathsf{\varphi }={\mathsf{\varphi }}_2-{\mathsf{\varphi }}_1=\frac{4\mathsf{\pi }\mathrm{\Delta }\mathrm{R}}{\mathsf{\lambda }}$$

(4)

With the coregistering or combination of the differences in phase for the two observations, an interferogram can be constructed in which the phase difference information is strongly related to the topography of the terrain [38]. The deformations of the ground can then be mapped. The phase difference can be caused by five factors [52]:

$$\mathrm{\Delta }\mathsf{\varphi }=\mathrm{\Delta }{\mathsf{\varphi }}_{flat}+\mathrm{\Delta }\mathsf{\varphi }{\mathsf{\phi }}_{elevation}+\mathrm{\Delta }{\mathsf{\varphi }}_{displacement}+\mathrm{\Delta }{\mathsf{\varphi }}_{atmosphere}+\mathrm{\Delta }{\mathsf{\varphi }}_{noise}$$

(5)

• flat = the curvature of the Earth.
• elevation = topography.
• displacement = deformations of the terrain between the two acquisitions of information.
• atmosphere = differences in relative humidity, pressure, and temperature between information acquisitions.
• noise = changes over time in volume scattering and acquisition angles, etc.

The flat, atmosphere, and noise sources of error were removed during inteferometric analysis to be left with only those of interest: elevation (to generate digital models of the terrain) and displacement (to analyse changes occurring to the terrain). If the difference in phase caused by the topography is removed from the interferogram, the resulting difference must correspond to the pattern of deformation between the information acquisition dates.

A coherence band is also calculated. This shows how similar the two studied images are at the pixel level on a scale of 0 (low) to 1 (high). The loss of coherence between a pair of images is explained by the following equation [44,53]:

$$\mathsf{\gamma }={\mathsf{\gamma }}_{\mathrm{T}}\ast {\mathsf{\gamma }}_{\mathrm{G}}\ast {\mathsf{\gamma }}_{\mathrm{V}}\ast {\mathsf{\gamma }}_{\mathrm{P}}$$

(6)

where:

• ${\mathsf{\gamma }}_{\mathrm{T}}$ = Temporal factor. This cannot be avoided and is due to the differences in the ground occurring between the two information acquisitions; it is important to the present work.
• ${\mathsf{\gamma }}_{\mathrm{G}}$ = Geometric factor. This is caused by errors in the orbit of the satellites; it can be partially accounted for.
• ${\mathsf{\gamma }}_{\mathrm{V}}$ = Volumetric factor. Unavoidable due to the presence of vegetation.
• ${\mathsf{\gamma }}_{\mathrm{P}}$ = Processing factor. Caused by errors of calculation; this should be avoided.

In classic interferometry, with the objective of producing digital terrain models, coherence is very important, requiring a very high grade of similarity between acquisitions to produce reliable DTMs, but differential interferometry has the particularity that the very own object of study, the deformation of the earth surface, produces a loss of coherence, so the important thing is to reduce the loss of coherence from other sources, as can be observed in Equation (6).

2.4.1. Software

The ESA’s Sentinel Application Platform Zuhlke (SNAP) was used for the differential inteferometric analysis of pairs of TopSAR images [54,55], employing the SNAPHU algorithm produced at Stanford University for the unwrapping phase [56,57,58].

2.4.2. Coregistering and Interferogram Construction

To construct interferograms, interferometric pairs of SAR images were coregistered, aligning the images at the pixel level, using the elder of the two as the reference. This coregistering was performed using 21-point Bisinc interpolation, enhanced spectral display, and an orbital correction using the Sentinel precise polynomial degree 3 [59].

Once the images were aligned, their interferograms (intensity, phase, and coherence bands) were calculated, correcting the error caused by the curvature of the Earth using a fifth-degree polynomial and performing an interpolation of the capturing satellite’s orbit with a third-degree polynomial.

A peculiarity of TopSAR images is their acquisition by radar bursts [60]. Thus, to advance in the interferometric analysis, the spaces between the bursts were eliminated [61].

The image was then cropped to the area of study, drastically reducing the time required for the following operations.

Once the interferogram had been obtained, the dephasing owed to topography was eliminated using an SRTM 1 sec HGT digital terrain model. As the difference in phase caused by the topography has been removed from the interferogram, the resulting difference must correspond to the pattern of deformation between the information acquisition dates [38]. Finally, the effect of noise was reduced using Goldstein phase filtering [62].

Along with the difference phase band, the coherence band, used for the deformation raster horizontal and vertical adjustment, is obtained.

The unwrapping of the interferogram was performed using Chen and Zebker algorithm implemented in SNAPHU software [56,57,58] and then vertical displacement in each pixel was calculated using the following equation [63]:

$$\mathrm{d}=-\frac{\mathsf{\lambda }}{4\mathsf{\pi }\cos {\mathsf{ϴ}}_{\mathrm{inc}}}\mathrm{\Delta }{\mathsf{\varphi }}_{\mathrm{d}}$$

(7)

where:

• d = the deformation according to the vertical axis.
• λ = the wavelength used.
• $\mathrm{\Delta }{\mathsf{\varphi }}_{\mathrm{d}}$ = the phase displacement for each cell between data acquisitions.
• ϴ_inc = incident angle.

Differential interferometry can obtain deformation only in the direction of the line of sight (LOS), so vertical deformation is estimated as a vertical projection of the LOS in the vertical direction [64].

Finally, geometric correction of the terrain was performed to compensate for the distortions of the SAR image due to variations in the terrain with respect to the inclination of the sensor. The aim was to produce a two-dimensional representation that faithfully represented the real-world situation. This was performed using a Range-Doppler orthorectification method [65] to geode-codify the SAR images, considering information for the capturing satellite’s orbit, time records, conversion factors for the oblique distance to the Earth’s surface, and a reference digital terrain model to determine the exact location represented by each pixel.

The ETRS89 projection system was used to perform the above geometric correction. This provided two raster products—the terrain deformation between data acquisitions, and the coherence between the acquisitions.

The coherence band was used to make the horizontal adjustment to compensate for inaccuracies in the horizontal plane of the raster representation of the ground deformation. Such adjustment requires that the study plot has clear, high-contrast reference points surrounding it [44,53]. In the present work, it was necessary to perform for every pair of result rasters (coherence and deformation) a horizontal translation to match the areas of greatest coherence with the buildings around the plot.

The coherence band was also used to make the vertical adjustment to compensate for any vertical error in the estimation of local deformations. To make this adjustment, a cell with very high coherence, thus with little likelihood of being deformed, was chosen to provide the zero-deformation reference, then the value of its twin cell in the deformation raster was subtracted from the entire deformation raster. The values of control cells selected for the present work are presented in

Table 2. Values of the zero-deformation control points selected.

Vertical Adjustment (m)
Dates		Control Point Value	Control Point Coherence
28 April 2017	10 May 2017	−0.02574	0.98794
28 April 2017	22 May 2017	−0.05999	0.97597
28 April 2017	15 June 2017	−0.06752	0.96981
28 April 2017	9 July 2017	−0.16464	0.93012
10 May 2017	10 May 2017	−0.03508	0.97794
22 May 2017	3 June 2017	0.01051	0.96236
3 June 2017	15 June 2017	−0.03450	0.96150
15 June 2017	27 June 2017	0.03143	0.98654
27 June 2017	9 July 2017	−0.05201	0.98177
30 December 2017	11 January 2018	−0.00314	0.99392

.

This kind of “direct” correction is suitable only for studies of small areas. Studies of larger areas require the interpolation of multiple high-coherence cells [37]; this was not necessary for the present work.

2.4.3. Terrestrial LIDAR Measurements

Materials

The equipment used for terrestrial LIDAR measurements was a LIDAR FARO Focus 3D device. CloudCompare Software was used for image production.

Study Period

Observations were made every six days over three periods of time. Study period 1 lasted from 23 March 2017 until 10 April 2017. Data were recorded from one measuring station (providing 1 cloud of data points per model) and using four reference spheres. During this period a steady deformation was noticed moving away from the measuring station, but visual inspection suggested this did not appear to be caused by changes in the terrain. Indeed, this deformation was found to be the result of a non-perfect correspondence between the radial measuring directions used by the LIDAR system to take measurements on different dates. The data collected for this period were therefore rejected. In the next study period, two measuring stations (producing two clouds of data points per model) were used. These were separated by 12,3 m, such that their radial measuring directions crossed, minimising the effect of radial de-correlation (Figure 5).

Study period 2 lasted from 28 April 2017 until 15 July 2017. Study period 3 lasted from 30 December 2017 until 23 January 2018.

Table 3. Characteristics of LIDAR data collection used.

Date	Measured Using:
23 March 2017	1 measuring station and 4 spheres. Radial de-correlation. DATA REJECTED
29 March 2017	1 station and 4 spheres. Radial de-correlation. DATA REJECTED
4 April 2017	1 station and 4 spheres. Radial de-correlation. DATA REJECTED
10 April 2017	1 station and 4 spheres. Radial de-correlation. DATA REJECTED
16 April 2017	1 station and 4 spheres. Radial de-correlation. DATA REJECTED
28 April 2017	2 stations and 4 spheres (2 clouds per model)
4 May 2017	2 stations and 4 spheres (2 clouds per model)
10 May 2017	2 stations and 4 spheres (2 clouds per model)
16 May 2017	2 stations and 4 spheres (2 clouds per model)
22 May 2017	2 stations and 4 spheres (2 clouds per model)
28 May 2017	2 stations and 4 spheres (2 clouds per model)
3 June 2017	2 stations and 4 spheres (2 clouds per model)
9 June 2017	2 stations and 4 spheres (2 clouds per model)
15 June 2017	2 stations and 4 spheres (2 clouds per model)
21 June 2017	2 stations and 4 spheres (2 clouds per model)
27 June 2017	2 stations and 4 spheres (2 clouds per model)
3 July 2017	2 stations and 4 spheres (2 clouds per model)
9 July 2017	2 stations and 4 spheres (2 clouds per model)
15 July 2017	Repeated 16 July 2017 due to error in the memory card. 2 stations and 4 spheres (2 clouds per model)
30 December 2017	2 stations and 6 spheres (2 clouds per model)
5 January 2018	No data taken due to rain
11 January 2018	2 stations and 6 spheres (2 clouds per model)

summarises the characteristics of each data collection period.

Data Processing

Paired clouds of data points were finely registered by point pair alignment using at least four spheres complemented, when necessary, with fixed elements of the surrounding building until obtaining a root mean square error of <1 cm, a reasonable value for a natural surface [66]. Iterative Closest Point (ICP) algorithm was tested and discharged at the beginning of the study because in this specific case, the RMS obtained were higher than with the alignment with a pair of points. The resulting data clouds had approximately 35 million points for each measurement time (Figure 6).

Each model was then cropped to contain only the study area and processed using the CANUPO surface identification algorithm to eliminate vegetation [67,68]. The surface translator used by the algorithm was constructed using collected data points known to represent vegetation. The final data cloud consisted of 32 million points (Figure 7).

As a complementary measure, a manual analysis was performed to remove any remaining vegetation not detected by the CANUPO algorithm and to reinstate any ground that should not have been eliminated.

The Multiscale Model to Model Cloud Compare (M3C2) algorithm was then used to calculate the distances between the clouds of data points on different days, as shown in Figure 8, using calculation settings selected in function of the rugosity of the study area [66,69]. After testing several sets of parameters, the following set was selected:

• Normal diameter: 0.20 m
• Projected diameter: 0.05 m
• Maximum distance: 0.25–1.00 m
• Core points: 0.025 m

During the process, two possible sources of error were detected: (1) the turning mechanism of the LIDAR device, and (2) the assembly alignment process, both for clouds of data points for generating models, and between models when performing comparisons.

Errors of the same origin were combined using the root mean square technique [66]:

$${\mathrm{E}}_{\mathrm{T}}=\sqrt[2]{\frac{{\mathrm{E}}_1^2+{\mathrm{E}}_2^2+\dots +{\mathrm{E}}_{\mathrm{n}}^2}{\mathrm{n}}}$$

(8)

Errors of different origins were combined using the root sum squared technique [70]:

$${\mathrm{E}}_{\mathrm{T}}=\sqrt[2]{{\mathrm{E}}_1^2+{\mathrm{E}}_2^2+\dots +{\mathrm{E}}_{\mathrm{n}}^2}$$

(9)

This combination of errors is used as the allowable difference limit of the InSAR system in relation to the LIDAR system, but it does not mean that these errors are actual values, they are the maximum expected or possible errors as a result of the different maximum errors in each of the LIDAR measurements.

The distance calculation results were exported to a GIS using a raster image with 10 × 10 cm cell dimensions, providing a mean value of all the points for each cell. This is vital for normalising and homogenising values; not doing this would cause the different densities of data points in radial distributions to produce overestimations of any deformation values obtained close to the laser emitter, and underestimations of those obtained further away.

A raster image was finally obtained for comparison with the output data from the differential interferometry technique.

3. Results

3.1. Comparison of Clouds of Data Points Obtained Using the Terrestrial LIDAR System

Comparison of clouds of points obtained using the terrestrial LIDAR system allowed the production of raster images of ground deformation for each pair of dates contemplated. These raster images, which had 10 × 10 cm cell dimensions, were normalised against the 13.93 × 13.93 m cells to make possible the comparison with the results of the differential interferometry method: for each pair of monitoring dates, the mean of all cells for each dataset was determined, leaving four mean cells (A, B, C, and D) with which to make comparisons (Figure 9 and

Table 4. Deformation and error measurements by comparing LIDAR Clouds.

Lidar Estimated Errors and Deformation Measures (m)
Reference	Aligned	Error Position Reference Base			Error Position Aligned Base			EP Total	RMS				Max. Possible Error	Deformation Measures
		EP R1	EP R2	EP Ref.	EP A1	EP A4	EP Aligned		RMS CC Ref.	RMS CC Aligned	RMS CC Compare	RMS CC Total		AL	BL	CL	DL	Global (A B C D)L
Absolute		R1	R2	*	A1	A2	**	***	RMSR	RMSA	RMSC	****	*****	AL	BL	CL	DL
28 April 2017	10 May 2017	0.0021	0.0039	0.0031	0.0024	0.0022	0.0023	0.0027	0.0068	0.0073	0.0110	0.0085	0.0090	−0.0003	−0.0029	0.0017	0.0016	0.0000
28 April 2017	22 May 2017	0.0021	0.0039	0.0031	0.0022	0.0486	0.0344	0.0244	0.0068	0.0125	0.0093	0.0098	0.0263	−0.0122	−0.0136	−0.0043	−0.0018	−0.0080
28 April 2017	15 June 2017	0.0021	0.0039	0.0031	0.0019	0.0018	0.0018	0.0025	0.0068	0.0088	0.0083	0.0080	0.0084	0.0025	−0.0033	0.0132	−0.0012	0.0028
28 April 2017	9 July 2017	0.0021	0.0039	0.0031	0.0044	0.0022	0.0035	0.0033	0.0068	0.0083	0.0117	0.0091	0.0097	−0.0024	−0.0075	0.0101	0.0026	0.0007
Incremental
28 April 2017	10 May 2017	0.0021	0.0039	0.0031	0.0022	0.0024	0.0023	0.0027	0.0068	0.0073	0.0110	0.0085	0.0090	−0.0003	−0.0029	0.0017	0.0016	0.0000
10 May 2017	22 May 2017	0.0022	0.0024	0.0023	0.0022	0.0486	0.0344	0.0244	0.0073	0.0125	0.0098	0.0101	0.0264	−0.0119	−0.0112	−0.0057	−0.0032	−0.0080
22 May 2017	20170603	0.0022	0.0486	0.0344	0.0041	0.0019	0.0032	0.0244	0.0125	0.0090	0.0102	0.0107	0.0267	0.0191	0.0153	0.0170	−0.0008	0.0127
3 June 2017	15 June 2017	0.0041	0.0019	0.0032	0.0019	0.0018	0.0018	0.0026	0.0090	0.0088	0.0083	0.0087	0.0091	−0.0040	−0.0039	0.0012	0.0014	−0.0013
15 June 2017	27 June 2017	0.0019	0.0018	0.0018	0.0021	0.0021	0.0021	0.0020	0.0088	0.0109	0.0083	0.0094	0.0096	−0.0013	−0.0008	−0.0014	−0.0010	−0.0011
27 June 2017	9 July 2017	0.0021	0.0021	0.0021	0.0044	0.0022	0.0035	0.0029	0.0109	0.0124	0.0111	0.0115	0.0118	−0.0044	−0.0035	−0.0023	0.0048	−0.0014
30 December 2017	11 January 2018	0.0025	0.0021	0.0023	0.0020	0.0020	0.0020	0.0021	0.0101	0.0053	0.0061	0.0075	0.0078	−0.0016	0.0017	0.0027	−0.0018	0.0002

* EPR = [(R1² + R2²+ … +Rn²)/n]^1/n; ** EPA = [(A1² + A2² + … +An²)/n]^1/n; *** EPT = [(EPR² + EPA²)/2]^1/2; **** RMST = [(RMSR² + RMSA² + RMSC²)/3]^1/2; ***** MPE = (EPT² + RMST²)^1/2.

).

For each of the pairs of dates, a maximum potential error was calculated from the errors detected:

• The turning mechanism of the LIDAR system at measuring station 1 on the first (reference) date (Error station R1)
• The turning mechanism of the LIDAR system at measuring station 2 on the first (reference) date (Error station R2)
• The alignment of the clouds of points (measuring stations 1 and 2) on the first (reference) date (RMS CC Reference)
• The turning mechanism of the LIDAR system at measuring station 2 on the first (reference) date (Error station A1)
• The turning mechanism of the LIDAR system at measuring station 2 on the second (aligned) date (Error station A2)
• The alignment of the clouds of data points (measuring stations 1 and 2) on the second (aligned) date (RMS CC Aligned)
• The alignment of the resulting models for the first and second (references vs aligned) dates (RMS CC Compare)

All the errors recorded were of similar magnitude, except on 22 May 2017 when there was an unexpected error in the turning mechanism of the LIDAR device at the second measuring station; the differences with respect to the differential interferometry results were exceptionally large on this date. The deformation values obtained were of the order of magnitude of millimetres and tenths of a millimetre, which seemed coherent with the erosion observed in the field in the visual inspections (except for the results that involves the measure on 22 May 2017).

The resulting combination of errors for each date is used as the allowable difference limit of the InSAR system in relation to the LIDAR system, but it does not mean that these errors are actual values, they are the maximum expected or possible error as a result of the different maximum errors in each of the LIDAR measurements.

3.2. Differential Inteferometric Analysis of TopSAR Images

For each pair of dates, two products were obtained: a deformation band and a coherence band. The graphical information within the coherence band was essential for confirming the concordance of high coherence areas (white) with buildings, and low coherence (dark) with wasteland and pasture. For each interferometric pair of SAR images, the displacement in the horizontal plane was determined with the coherence band and used to correct the deformation band.

The numerical results revealed mild erosive processes. Some of the results, however, were due to other processes. For example, between 4 May and 16 May 2017, some gravel was dumped in the plot, affecting cell C (see the distribution of cells in Figure 10), and between 27 June and 9 July 2017, some soil was also dumped, affecting cells A, C, and D.

The deformation values obtained were of the order of magnitude of millimetres and tenths of a millimetre, which seemed coherent with the erosion observed in the field in the field inspections.

Figure 10 shows the graphical representation of the differential interferometry results:

• First column shows the dates between deformations.
• Second-column deformation obtained by the differential interferometry method between the cited dates. Green represents a positive deformation, and red a negative deformation. In this figure, the deformation raster has been already corrected in the horizontal plane and in the vertical axis and its value is indicated in the own cell.
• Third-column shows the coherence values obtained using the differential interferometry method between the dates mentioned, vital for making horizontal adjustments and for locating the zero-deformation control points for making vertical adjustments. In this figure the coherence band has been already corrected in the horizontal plane. The values vary between 0 and 1, where 1 (white) represents 100% coherence, and 0 (black) 0% coherence.
• Fourth-column contains the legends for deformation and coherence values.

Table 5. Deformation measurements by Differential Interferometry over IW TopSAR images, with and without vertical adjustment using zero-deformation control points (very high coherence cells).

INSAR TOPSAR Measures (m)
Reference	Aligned	Vertical Error Adjustment (m)		Original Measure (m)				Corrected Measure (m)
		Control Point Value	Control Point Coherence	A	B	C	D	A’	B’	C’	D’	Global A’B’C’D’
Absolute		x		A	B	C	D	A’ = (A − x)	B’ = (B − x)	C’ = (C − x)	D’ = (D − x)
28 April 2017	10 May 2017	−0.02574	0.98794	−0.02383	−0.02348	−0.02335	−0.02358	0.00192	0.00226	0.00239	0.00216	0.00218
28 April 2017	22 May 2017	−0.05999	0.97597	−0.05749	−0.05854	−0.05835	−0.05815	0.00250	0.00145	0.00164	0.00184	0.00186
28 April 2017	15 June 2017	−0.06752	0.96981	−0.06637	−0.06718	−0.06762	−0.06716	0.00114	0.00033	−0.00010	0.00035	0.00043
28 April 2017	9 July 2017	−0.16464	0.93012	−0.16202	−0.16218	−0.16165	−0.16452	0.00262	0.00246	0.00298	0.00012	0.00204
Incremental
28 April 2017	10 May 2017	−0.02574	0.98794	−0.02383	−0.02348	−0.02335	−0.02358	0.00192	0.00226	0.00239	0.00216	0.00218
10 May 2017	22 May 2017	−0.03508	0.97794	−0.03582	−0.03584	−0.03493	−0.03517	−0.00074	−0.00076	0.00014	−0.00010	−0.00036
22 May 2017	3 June 2017	0.01051	0.96236	0.01247	0.01243	0.01185	0.01163	0.00196	0.00191	0.00133	0.00112	0.00158
3 June 2017	15 June 2017	−0.03450	0.96150	−0.03636	−0.03106	−0.03679	−0.03574	−0.00186	0.00344	−0.00229	−0.00124	−0.00049
15 June 2017	27 June 2017	0.03143	0.98654	0.03267	0.03113	0.03258	0.03207	0.00125	−0.00030	0.00115	0.00064	0.00068
27 June 2017	9 July 2017	−0.05201	0.98177	−0.04979	−0.05145	−0.05039	−0.05052	0.00222	−0.00056	0.00162	0.00149	0.00147
30 December 2017	11 January 2018	−0.00314	0.99392	−0.00244	−0.00211	−0.00341	−0.00438	0.00070	0.00103	−0.00027	−0.00124	0.00006

shows the numerical results of the differential interferometry results, including the original measurements, the vertical correction, and the final values:

3.3. Comparison of Differential Interferometry and LIDAR Deformation Measures

As can be observed in

Table 6. Comparison of the results obtained by the two methods.

InSAR TopSAR Measurements vs. LIDAR Measurements (m)
Reference	Aligned	Permissible Range (LIDAR Max. Error)	Difference between LIDAR and InSAR
			ΔA	ΔB	ΔC	ΔD	Global Difference (A B C D)
Absolute			ΔA = (A’ − AL)	ΔB = (B’ − BL)	ΔC = (C’ − CL)	ΔD = (D’ − DL)
28 April 2017	10 May 2017	0.0090	0.0022	0.0052	0.0007	0.0006	0.0022
28 April 2017	22 May 2017	0.0263	0.0147	0.0151	0.0060	0.0036	0.0098
28 April 2017	15 June 2017	0.0084	0.0014	0.0036	0.0133	0.0015	0.0024
28 April 2017	09 July 2017	0.0097	0.0050	0.0100	0.0071	0.0025	0.0014
Incremental
28 April 2017	10 May 2017	0.0090	0.0022	0.0052	0.0007	0.0006	0.0022
10 May 2017	22 May 2017	0.0264	0.0111	0.0104	0.0059	0.0031	0.0076
22 May 2017	03 June 2017	0.0267	0.0171	0.0134	0.0157	0.0019	0.0111
03 June 2017	15 June 2017	0.0091	0.0021	0.0073	0.0035	0.0026	0.0008
15 June 2017	27 June 2017	0.0096	0.0025	0.0005	0.0026	0.0017	0.0018
27 June 2017	09 July 2017	0.0118	0.0067	0.0041	0.0039	0.0033	0.0028
30 December 2017	11 January 2018	0.0078	0.0023	0.0007	0.0029	0.0005	0.0002

, only two of a total of forty observations surpassed the range of maximum possible error associated with the LIDAR system, both were in periods longer than 12 days and it was for just 0.0003 and 0.0049 m. From 12 to 12 days, there were none exceeding one.

All the values were in the same range, and as the same errors are more important as lower is the erosion event, the potential error associated with the technique must be small.

4. Discussion

There are two objectives pursued by the current research, accuracy of measurements and coherence loss boundary:

• Concerning the accuracy of measurements, deformation values obtained with both systems were of the order of magnitude of millimetres and tenths of a millimetre, consistent with the erosion observed in the field in the visual inspections, except for the LIDAR results that involve the measurement on 22 May 2017, because of an error in the turning mechanism of LIDAR in position 2. Only two of a total of forty observations surpassed the range of maximum possible error associated with the LIDAR system, both were in periods longer than 12 days (48 and 86 days) and it was for just 0.0003 and 0.0049 m. From 12 to 12 days, there were none exceeding one.
• Concerning the loss of coherence, the results show how coherence is lost as the period between acquisitions increases, so deformation measurements lose accuracy. As exposed above, limiting the time between acquisition to 12 days guarantees that deformation measures are accurate, so the coherence loss corresponds directly to the deformation registered between acquisition, and not to other factors.

The numerical results revealed mild erosive processes. Some of the results, however, were due to other processes. For example, between 4 May and 16 May 2017, some gravel was dumped in the plot, affecting cell C, and then between 27 June and 9 July 2017, some soil was also dumped, also affecting cells A, C, and D.

Small plots such as those used in this study facilitate the work with the terrestrial LIDAR but make it difficult to ignore the impact of horizontal inaccuracies of TopSAR images in a small study area. The horizontal correction was made using coherence band information, thanks to a special characteristic of the chosen plot: surrounded by many reference points as buildings and pavemented areas. The graphical information within the coherence band of TopSAR differential interferometry was used to compare the concordance of high-coherence areas with buildings, and low-coherence areas with wasteland and pasture. For each interferometric pair of SAR images, the displacement in the horizontal plane was determined with the coherence band and was used to horizontally correct the deformation band.

The quality of the results in the raster deformation (vertical deformations) depends on adjusting its values by means of zero-vertical deformation control points, a nearby cell with a high coherence value that makes sure that no deformation has occurred between acquisitions.

Finally, the presence of vegetation was a problem for the gold standard LIDAR method and required using the CANUPO surface identification algorithm to eliminate vegetation, but not for the proposed interferometric method, in which C band microwave radiation signals are not affected by small vegetation (while dealing with herbaceous vegetation, not a shrub or tree vegetation).

5. Conclusions

Differential interferometry on TopSAR images of the Sentinel-1 mission has proven to be a real option to monitor tillage and water erosion with the following considerations.

• The study period of time between acquisition should not exceed 12 days, to guarantee that all the loss of coherence is because of the deformation of the earth’s surface and not from other sources.
• Longer periods can be studied by combining consecutive 12-day steps.
• Small vegetation such as cereal crops does not interfere with InSAR measurements, but larger vegetation such as sunflowers can do so. Verification of the interference of larger vegetation will be the subject of future research.
• The accuracy of the results is highly related to the existence of reliable zero-deformation control points (very high coherence points), used to adjust the deformation in the study area.
• The different soil moisture between acquisitions can alter the deformation measurements by affecting the expansion capability of the clay or the soil dielectric characteristic and its microwave reflecting capability. This effect will be the subject of future research, but meanwhile, it can be avoided by comparing acquisitions with similar soil moisture degrees.
• The effects of soil compaction/decompaction can also affect the measurements because this methodology is based on changes in the earth’s surface but does not consider changes in density. This effect will also be the subject of future research.

The spatial nature of this methodology, capable to analyse large basins using a grid with cells of only 13.93 m, becomes it a powerful tool not only to monitor the basin erosion but also to analyse the efficiency of different managements or anti-erosion measures proved in the field against short- to long- term rainy periods, making it possible to design and verify more efficient managing strategies and anti- erosion measures.