Automatic Grassland Cutting Status Detection in the Context of Spatiotemporal Sentinel-1 Imagery Analysis and Artificial Neural Networks

Grassland contributes to carbon storage and animal feed production. Its yield is largely determined by the cutting times of grassland. Previous studies have used remote sensing data for grassland biomass estimation, but only a few studies have focused on SAR remote sensing approaches for automatic grassland cutting status detection. Due to the occurrence of multiple cuttings in a year, it is crucial to effectively monitor grassland cutting events in order to achieve accurate biomass estimations of a whole season. In this study, we examined the capabilities of multilayer perceptron neural networks for automatic grassland cutting status detection using SAR imagery. The proposed model inputs are a time series dataset of VV and VH Sentinel-1 C-band SAR and second-order texture metrics (homogeneity, entropy, contrast and dissimilarity). The proposed approach has been successfully tested on a dataset collected from several fields in Germany in 2016, with an overall accuracy of 85.71% for the validation set.


Introduction
Grassland covers 40% of the global landmass and is a primary contributor to the production of forage for dairy and meat industries [1,2].Further, it stores about 34% of the global carbon stock [3].Therefore, grassland can be considered a key component in regulating the global carbon cycle [4] and animal biodiversity [5,6].In order to preserve grassland where regular cutting is not economically feasible, the European Union (EU) has decided to pay subsidies to farmers [7].These payments are scaled according to cuttings in a precisely mapped area in a given year, which have highlighted the need for accurate mapping of cutting events.Moreover, the monitoring of grassland cutting events is crucial considering their ecological and agricultural importance as well as their key role in modeling yearly grassland yield [8].As in-situ observations are time-consuming and costly, it would be advantageous to find a more cost-and time-efficient method to determine whether, when and to what extent a grassland has been mown.To address the aforementioned issues, remote sensing may be used as a more cost-and time-effective technique that provides information on cutting status (cut or uncut fields) at a high frequency over large spatial scales [9].
In the remote sensing-related literature, grassland has not been studied as thoroughly as other agricultural crops [10].Most of the remote sensing approaches have focused on the agronomy of grassland [1,4], such as growth rate modeling [11] or biophysical parameter retrieval, including chlorophyll content [12] and Leaf Area Index (LAI) [13].These approaches rely on the use of optical remote sensing.Synthetic Aperture Radar (SAR) images have not been explored extensively in grassland studies [12,14].Since SAR data remain relatively independent of perturbances in the region of Middle Franconia in northern Bavaria at an elevation of approximately 420 m above sea level on relatively flat lowland terrain.The typical soil type of the region is brown earth, which originated from weathered sand-and claystone [34].
The second site (Kempten) is located in southern Bavaria, east of Lake Constance and close to the Alps.This site is found in the Alpine Foreland at an elevation of around 750 m above sea level on relatively flat terrain of a former river basin in an otherwise hilly region.The typical soil type of the region is brown and para brown earth, which originated from limestone and loamy moraine sediments [34].
The mean temperatures of 2016 were above the long-term average at both test sites (8.9 and 8.6 °C).The total precipitation was significantly below average in Triesdorf (596 mm) and slightly above average in Kempten (1330 mm).The 30-year average (1981-2010) temperatures were 8.3 and 7.6 °C while the average yearly precipitations were 751 and 1261 mm for Triesdorf and Kempten, respectively [35].According to the Koeppen Geiger climate classification, the climate for the study areas can be classified as a temperate, oceanic type (warm summer, wet winter; Climate zone Cfb).

Remotely Sensed Data
The dataset that we used consists of 50 Sentinel-1 Ground Range Detected (GRD), Interferometric Wide swath mode (IW) images captured in 2016.Focused Sentinel-1 C-band SAR data was detected, multi-looked and projected to ground range using an Earth ellipsoid model in order to create Level-1 GRD products.The resulting products have approximately square pixel spacing (~10 m in range and azimuth) with reduced speckles.We used the images from acquisition dates every nine days (with some exceptions) from three geometries defined by the Relative Orbit Number (RON) (Table 1).Figure 2 summarizes the acquisition dates for each geometry.The second site (Kempten) is located in southern Bavaria, east of Lake Constance and close to the Alps.This site is found in the Alpine Foreland at an elevation of around 750 m above sea level on relatively flat terrain of a former river basin in an otherwise hilly region.The typical soil type of the region is brown and para brown earth, which originated from limestone and loamy moraine sediments [34].
The mean temperatures of 2016 were above the long-term average at both test sites (8.9 and 8.6 • C).The total precipitation was significantly below average in Triesdorf (596 mm) and slightly above average in Kempten (1330 mm).The 30-year average (1981-2010) temperatures were 8.3 and 7.6 • C while the average yearly precipitations were 751 and 1261 mm for Triesdorf and Kempten, respectively [35].According to the Koeppen Geiger climate classification, the climate for the study areas can be classified as a temperate, oceanic type (warm summer, wet winter; Climate zone Cfb).

Remotely Sensed Data
The dataset that we used consists of 50 Sentinel-1 Ground Range Detected (GRD), Interferometric Wide swath mode (IW) images captured in 2016.Focused Sentinel-1 C-band SAR data was detected, multi-looked and projected to ground range using an Earth ellipsoid model in order to create Level-1 GRD products.The resulting products have approximately square pixel spacing (~10 m in range and azimuth) with reduced speckles.We used the images from acquisition dates every nine days (with some exceptions) from three geometries defined by the Relative Orbit Number (RON) (Table 1).Figure 2 summarizes the acquisition dates for each geometry.

Ground Truth Campaign
We have used a total of five fields in each of our study sites (Triesdorf and Kempten).The investigated grasslands are intensively used and cut four to five times throughout a season.The typical cutting times were early May, mid of June, late July and October while there was sometimes a fifth cutting in early September.Having been run in parallel to the Sentinel-1 acquisitions, the field surveys of the Bavarian state research center for agriculture (LfL) provided the exact dates of cutting events for each test field (e.g., on 25 August 2017, field one in Kempten is cut), enabling us to determine the cutting status in each image.

Data Preprocessing
For preprocessing of SAR data, we used the open source Sentinel-1 toolbox (S1TBX) that is provided by the European Space Agency (ESA).Figure 3 shows a schematic view of the preprocessing flow.Since the proposed model relies on highly accurate radiometric values, we removed thermal noise.To improve the co-registration accuracy, we applied the precise orbit model to each image.This provides information about the geometric location of the SAR images that improves the position accuracy during co-registration.In the next phase, we calibrated all images to the backscatter coefficient.After this, the radiometric terrain flattening using Shuttle Radar Topography Mission (SRTM) 3 sec data and bilinear interpolation was applied to the dataset.For speckle filtering, a 3 × 3 Lee filter was applied and finally, the images were georeferenced using a digital terrain model (Range Doppler Terrain Correction, SRTM) and re-projected to Universal Transverse Mercator (UTM), WGS84.For further analysis, we transformed the corrected intensity data to the logarithmic scale (unit dB).

Ground Truth Campaign
We have used a total of five fields in each of our study sites (Triesdorf and Kempten).The investigated grasslands are intensively used and cut four to five times throughout a season.The typical cutting times were early May, mid of June, late July and October while there was sometimes a fifth cutting in early September.Having been run in parallel to the Sentinel-1 acquisitions, the field surveys of the Bavarian state research center for agriculture (LfL) provided the exact dates of cutting events for each test field (e.g., on 25 August 2017, field one in Kempten is cut), enabling us to determine the cutting status in each image.

Data Preprocessing
For preprocessing of SAR data, we used the open source Sentinel-1 toolbox (S1TBX) that is provided by the European Space Agency (ESA).Figure 3 shows a schematic view of the preprocessing flow.Since the proposed model relies on highly accurate radiometric values, we removed thermal noise.To improve the co-registration accuracy, we applied the precise orbit model to each image.This provides information about the geometric location of the SAR images that improves the position accuracy during co-registration.In the next phase, we calibrated all images to the backscatter coefficient.After this, the radiometric terrain flattening using Shuttle Radar Topography Mission (SRTM) 3 sec data and bilinear interpolation was applied to the dataset.For speckle filtering, a 3 × 3 Lee filter was applied and finally, the images were georeferenced using a digital terrain model (Range Doppler Terrain Correction, SRTM) and re-projected to Universal Transverse Mercator (UTM), WGS84.For further analysis, we transformed the corrected intensity data to the logarithmic scale (unit dB).

ANN Model Description
Artificial Neural Networks (ANNs) are useful in modeling a variety of nonlinear behaviors [36].ANNs are composed of neurons that are nonlinear computational elements.The neurons establish inverse mapping, which discriminates the relations between inputs and outputs during the training phase.Several studies point out that ANNs are robust in terms of computational speed, stability and accuracy with respect to other investigated algorithms [36,37].
One of the most popular ANNs topologies is the feedforward topology in which each neuron of a layer is connected to all neurons of the successive layer but has no feedback to neurons in the previous layers [38].In feedforward ANNs, the weighted sum of input datasets is computed with an activation function by each neuron in the network.Equation 1 presents the mathematical representation of ANNs [38]: where d1, . . ., dn are inputs; w1, . . ., wn are associated connection weights; and b is the bias value.The most commonly used activation function is the sigmoid function.An example of the sigmoid function is the logistic function that is in the range of [0, 1]: where a is the curve steepness for v.The hyperbolic tangent function provides better results than those obtained with the sigmoid function of Equation 2. Therefore, in this study, we used an antisymmetric form of sigmoid function with the range from −1 to 1 (hyperbolic tangent function) (Equation (3)).
Figure 4 provides the backpropagation algorithm for the training of the network, which is very common in remote sensing studies [39,40].The backpropagation algorithm minimizes the Mean Square Error (MSE) between the multilayered feedforward perception output and the desired output by using a gradient search technique [38].

ANN Model Description
Artificial Neural Networks (ANNs) are useful in modeling a variety of nonlinear behaviors [36].ANNs are composed of neurons that are nonlinear computational elements.The neurons establish inverse mapping, which discriminates the relations between inputs and outputs during the training phase.Several studies point out that ANNs are robust in terms of computational speed, stability and accuracy with respect to other investigated algorithms [36,37].
One of the most popular ANNs topologies is the feedforward topology in which each neuron of a layer is connected to all neurons of the successive layer but has no feedback to neurons in the previous layers [38].In feedforward ANNs, the weighted sum of input datasets is computed with an activation function by each neuron in the network.Equation 1 presents the mathematical representation of ANNs [38]: where d 1 , . . ., d n are inputs; w 1 , . . ., w n are associated connection weights; and b is the bias value.The most commonly used activation function is the sigmoid function.An example of the sigmoid function is the logistic function that is in the range of [0, 1]: where a is the curve steepness for v.The hyperbolic tangent function provides better results than those obtained with the sigmoid function of Equation 2. Therefore, in this study, we used an anti-symmetric form of sigmoid function with the range from −1 to 1 (hyperbolic tangent function) (Equation ( 3)).
Figure 4 provides the backpropagation algorithm for the training of the network, which is very common in remote sensing studies [39,40].The backpropagation algorithm minimizes the Mean Square Error (MSE) between the multilayered feedforward perception output and the desired output by using a gradient search technique [38].

Model Implementation
In this study, we used multilayer perceptron neural networks, which are the most suitable topology for information retrieval [38,41].We have used VV and VH polarization backscattering intensities that are expected to be more suitable for grassland cutting status monitoring compared to HH and HV polarization [42,43].Therefore, six of the input parameters include VV and VH polarization backscattering intensity for the checking date; VV and VH polarization backscattering intensity for nine days before the checking date; and VV and VH polarization backscattering intensity for eighteen days before the checking date.
Further input variables are the second-order texture metrics (Table 2), including Homogeneity, Entropy, Contrast and Dissimilarity (for the checking date), that are calculated for each polarization (eight parameters in total for both VV and VH channels).Second-order texture metrics show the number of occurrences of the relationship between a pixel and its specified neighboring pixels [44].We used a 3 × 3 kernel as it is considered to be suitable for computing the texture matrix and extracting enough spatial information to distinguish among different land features [45].All the input parameters were normalized to a range [−1, 1] using the following equation [40]: where Dnew is the normalized value, Dmin is the minimum value and Dmax is the maximum value in the parameter dataset.

Model Implementation
In this study, we used multilayer perceptron neural networks, which are the most suitable topology for information retrieval [38,41].We have used VV and VH polarization backscattering intensities that are expected to be more suitable for grassland cutting status monitoring compared to HH and HV polarization [42,43].Therefore, six of the input parameters include VV and VH polarization backscattering intensity for the checking date; VV and VH polarization backscattering intensity for nine days before the checking date; and VV and VH polarization backscattering intensity for eighteen days before the checking date.
Further input variables are the second-order texture metrics (Table 2), including Homogeneity, Entropy, Contrast and Dissimilarity (for the checking date), that are calculated for each polarization (eight parameters in total for both VV and VH channels).Second-order texture metrics show the number of occurrences of the relationship between a pixel and its specified neighboring pixels [44].We used a 3 × 3 kernel as it is considered to be suitable for computing the texture matrix and extracting enough spatial information to distinguish among different land features [45].All the input parameters were normalized to a range [−1, 1] using the following equation [40]: where D new is the normalized value, D min is the minimum value and D max is the maximum value in the parameter dataset.Overall, the dataset consists of 70 samples (checking dates), which were utilized for training/testing and validating the model.We randomly selected 60% and 10% of the data set for training and testing, respectively.The remaining samples (10 and 11 samples for Triesdorf and Kempten, respectively) were used for final model validation.The network training was repeated ten times for each topology.Subsequently, we selected the network with the highest accuracy, which implies that the retained model resides within the best 10% of the distribution of all possible models [46].The early stopping procedure was used to determine the stop point for the training since it enables the optimization of the generalization properties of the network to avoid overfitting [38].The early stopping procedure constantly evaluates the performance of the network during the training process.The number of hidden nodes was optimized by trial and error and ranged between 5 and 35 depending on the number of input parameters.The output function of the network returns 0 for a detected cutting and 1 if no cutting occurred.The simulation of the ANNs was conducted with the Stuttgart Neural Network Simulator (SNNS) [47].

The Results
In the final step of the model development, we analyzed different topologies by changing the input parameters of the model in order to assess the influence of different sets of inputs on model accuracy.We started by considering only VV and VH parameters for the checking date as the input layer (called "a" in Table 3) of the models.After this, we increased the input parameters to six and ten by adding the VV and VH parameters of 9 days and 18 days before the checking dates (called "b" in Table 3) and the homogeneity and entropy parameters (for VV and VH polarizations) for the checking date (called "c" in Table 3), respectively.Finally, we increased the input parameters to fourteen by adding contrast and dissimilarity parameters (for VV and VH polarizations) for the checking date (called "d" in Table 3).For each input set, the topology and calibration with the highest performance in terms of overall accuracy and MSE has been selected and presented in Table 3. Table 3 also highlights the model architecture (no. of input parameters -no. of nodes in the hidden layer -output) of the best models.Cohen's Kappa for the input set "d" (Table 4) shows that the proposed algorithm achieves acceptable detection results for the validation dataset (with an overall accuracy of 85.71%).For cutting events, the commission and omission errors are 25.00% and 14.29%, respectively.Growing (no cutting) events show a commission error of 7.7% and an omission error of 14.29%.

Discussion
The results presented in Table 3 demonstrate that increasing the number of components in the input layer progressively decreased the MSE.We observed a significant improvement when the second-order texture metrics of the checking date were added to the input layer (MSE reduced by 3.14), which indicates that the ANN is able to constructively combine the information of each input node.The best overall accuracy of 85.71% is obtained for the model generated by the input set "d", which results in an increase of 71.43% in the overall accuracy compared to the best model using VV and VH parameters for the checking date (input set a).The overall accuracies increased by 52.38% and 57.14% when the homogeneity/entropy (input set b) and contrast/dissimilarity parameters (for VV and VH polarizations), respectively, were used along with the other parameters (input set c).
For model evaluation, it is essential to identify the conditions in which the proposed approach generates poor accuracies.The highest commission error (25%) occurred at a cutting event at Triesdorf on October 27 th , 2016.For this date, no significant correlation could be established between the backscattering coefficient and the vegetation biomass.Since various factors, i.e., systemdependent (e.g., radar wavelength, mean incidence angle, polarization) and target/weatherdependent (e.g., wind speed, temperature and precipitation) influence the microwave scattering of grassland, it is difficult to isolate a single factor.
Weather clearly affects both structure and water content of grass significantly.For example, Wooding et al. proved that rainfall of about 10 mm prior to SAR acquisition has led to increases of 1-4 dB in the backscattering of grass, bare soil and wheat-bare soil cover types [48] while Wind affects plant microclimates.Temperature changes are important due to the highly differing dielectric properties of ice and liquid water.Moreover, different species associations showed differences in the backscatter response that were partially related to herbage height.Furthermore, these differences  Cohen's Kappa for the input set "d" (Table 4) shows that the proposed algorithm achieves acceptable detection results for the validation dataset (with an overall accuracy of 85.71%).For cutting events, the commission and omission errors are 25.00% and 14.29%, respectively.Growing (no cutting) events show a commission error of 7.7% and an omission error of 14.29%.

Discussion
The results presented in Table 3 demonstrate that increasing the number of components in the input layer progressively decreased the MSE.We observed a significant improvement when the second-order texture metrics of the checking date were added to the input layer (MSE reduced by 3.14), which indicates that the ANN is able to constructively combine the information of each input node.The best overall accuracy of 85.71% is obtained for the model generated by the input set "d", which results in an increase of 71.43% in the overall accuracy compared to the best model using VV and VH parameters for the checking date (input set a).The overall accuracies increased by 52.38% and 57.14% when the homogeneity/entropy (input set b) and contrast/dissimilarity parameters (for VV and VH polarizations), respectively, were used along with the other parameters (input set c).
For model evaluation, it is essential to identify the conditions in which the proposed approach generates poor accuracies.The highest commission error (25%) occurred at a cutting event at Triesdorf on October 27th, 2016.For this date, no significant correlation could be established between the backscattering coefficient and the vegetation biomass.Since various factors, i.e., system-dependent (e.g., radar wavelength, mean incidence angle, polarization) and target/weather-dependent (e.g., wind speed, temperature and precipitation) influence the microwave scattering of grassland, it is difficult to isolate a single factor.
Weather clearly affects both structure and water content of grass significantly.For example, Wooding et al. proved that rainfall of about 10 mm prior to SAR acquisition has led to increases of 1-4 dB in the backscattering of grass, bare soil and wheat-bare soil cover types [48] while Wind affects plant microclimates.Temperature changes are important due to the highly differing dielectric properties of ice and liquid water.Moreover, different species associations showed differences in the backscatter response that were partially related to herbage height.Furthermore, these differences were significantly influenced by canopy structure and underlying soil condition and biomass water content [48].
Therefore, we checked the weather data before and during October 27.Since the weather stations in the study site did not record any rainfall or wind speed five days before Sentinel-1 overpass, a direct influence of rainfall or wind can be neglected.The species composition of the test fields is unknown.However, we assume that the species composition does not significantly change between the third/fourth and fifth cut.In spring and early summer, early cutting may positively influence the occurrence of rapidly growing grasses while smaller growing grasses and herbs may be reduced.Therefore, the species composition may be of special interest when analyzing different grassland types.
Furthermore, soil moisture influences the backscattering signal.Schieche et al. proved that the soil water content has little influence on C-Band backscatter in the presence of dense vegetation [25], which may explain why the proposed approach failed to detect a cutting event at Triesdorf on October 27th, 2016.In Germany, the grassland fields that were cut 4-5 times a year show a low amount of biomass in autumn (the field measurements indicate less than 100 grams/sqm) in comparison with the biomass amount in April-May (more than 3.0 kg/sqm of biomass on the study sites).We assume that there are different contributions of soil backscatter in October and late spring/summer [49,50].Therefore, a sparsely vegetated field with a high amount of soil contributing to the total backscatter probably caused the outlier event on October 27.
We assume that model performance can be improved by adding further input parameters, such as the average incidence angle.A shallow incidence angle data would allow us to retrieve more vegetation signals.Furthermore, the underlying soil is more evident at steeper incidence angles as the SAR pulse travel length through the vegetation layer increases with an increase in the incidence angle.Moreover, Ulaby et al. demonstrated the effect of varying incidence angles on backscattering coefficients along the crop cycle [51].They mentioned that the difference can reach 0.35 dB 0−1 in C-band over bare soils and decreases to a value close to zero when the surface is fully vegetated (about 0.05 dB 0−1 ).Tamm et al. demonstrated that median VV and VH polarization coherence values were significantly higher after a mowing event [24].Therefore, the use of interferometric coherence could also be beneficial for increasing the performance of the model.

Conclusions
The aim of this study was to demonstrate the capability of multilayer perceptron neural networks as an automated method for grassland cutting status detection.There were four different sets of input parameters, which include (a) VV and VH parameters for the checking date; (b) VV and VH parameters of 9 days and 18 days before the checking dates plus set 'a'; (c) homogeneity and entropy parameters (for VV and VH polarizations) for the checking date plus set 'b'; and (d) contrast and dissimilarity parameters (for VV and VH polarizations) for the checking date plus 'c'.These sets of parameters were used for model development in order to assess the influence of different inputs on model accuracies.
We used four different setups with an increasing number of input parameters (2 to 14, including VV and VH parameters for the checking date as well as 9 days and 18 days before the checking dates, homogeneity, entropy, contrast and dissimilarity).
Twenty-one different models with varying numbers of hidden nodes (for each input set) were separately trained.The network with 14 input parameters and 31 hidden nodes performed best in terms of generalization, time load and accuracy (overall accuracy of 85.71% for the validation dataset).The results demonstrate that the machine learning inversion model is capable of detecting grassland cutting status with a good level of precision.Nevertheless, multilayer perceptron models can be enhanced by a larger and more diverse training/testing dataset.
Concerning polarization as a system-dependent influencing factor, several studies have reported stronger backscatter from grass in VV polarization compared to HH polarization, which interacts more strongly with broad-leaved canopies [52,53].The geometric and dielectric properties of different grassland types change rapidly over the year, which alters the backscatter signal response.For each

Figure 1 .
Figure 1.Map of Germany with locations of in-situ study areas highlighted in white.The colors represent relief elevation.

Figure 1 .
Figure 1.Map of Germany with locations of in-situ study areas highlighted in white.The colors represent relief elevation.

Figure 2 .
Figure 2. Sentinel-1 acquisitions used in the study.The squares and circles represent Triesdorf and Kempten sites, respectively.

Figure 2 .
Figure 2. Sentinel-1 acquisitions used in the study.The squares and circles represent Triesdorf and Kempten sites, respectively.

Figure 4 .
Figure 4.The backpropagation ANN model used in this study.The input layer (d) is connected to a layer of hidden nodes (n), which are in turn connected to the output layer (x).

Figure 4 .
Figure 4.The backpropagation ANN model used in this study.The input layer (d) is connected to a layer of hidden nodes (n), which are in turn connected to the output layer (x).

Figure 5
Figure 5 exemplarily presents the MSE of the test set "d" as a function of the number of hidden nodes.Considering generalization, time load and accuracy, the network with 31 nodes provided the lowest MSE compared with the other possibilities.In total, approximately 9000 training cycles were sufficient for training the network.The results of the accuracy assessment based on our validation dataset are shown in Table4.

Figure 5 .
Figure 5. Training MSE as a function of the number of hidden nodes for input set "d".

Figure 5 .
Figure 5. Training MSE as a function of the number of hidden nodes for input set "d".

Table 1 .
RONs and their parameters used in this study.The values are given for the study areas.

Table 1 .
RONs and their parameters used in this study.The values are given for the study areas.

Table 2 .
List of input variables used.Gr stands for the distinct grey level numbers.

Table 2 .
List of input variables used.Gr stands for the distinct grey level numbers.

Table 3 .
Comparison of results using different input parameter combinations; the model architecture [x-y-z] summarizes the number of input parameters [x], the number of nodes in the hidden layer [y] and the number of output nodes [z].

Table 4 .
Confusion matrix and accuracies of the validation set on input set "d".

Table 4 .
Confusion matrix and accuracies of the validation set on input set "d".