Utilization of Multi-Temporal Microwave Remote Sensing Data within a Geostatistical Regionalization Approach for the Derivation of Soil Texture

Abstract

Land Surface Models (LSM) have become indispensable tools to quantify water and nutrient fluxes in support of land management strategies or the prediction of climate change impacts. However, the utilization of LSM requires soil and vegetation parameters, which are seldom available in high spatial distribution or in an appropriate temporal frequency. As shown in recent studies, the quality of these model input parameters, especially the spatial heterogeneity and temporal variability of soil parameters, has a strong effect on LSM simulations. This paper assesses the potential of microwave remote sensing data for retrieving soil physical properties such as soil texture. Microwave remote sensing is able to penetrate in an imaged media (soil, vegetation), thus being capable of retrieving information beneath such a surface. In this study, airborne remote sensing data acquired at 1.3 GHz and in different polarization is utilized in conjunction with geostatistics to retrieve information about soil texture. The developed approach is validated with in-situ data from different field campaigns carried out over the TERENO test-site “North-Eastern German Lowland Observatorium”. With the proposed approach a high accuracy of the retrieved soil texture with a mean RMSE of 2.42 (Mass-%) could be achieved outperforming classical deterministic and geostatistical approaches.

1. Introduction

Climate Change is present and affects many ecological, social and economic developments. A new challenge to solve in climate change research is the prediction of large scale climate change to develop suitable adaption strategies [1]. The development of management and adaption strategies implies the utilization of robust Land Surface Models (LSM), as well as reliable model input parameters [2,3]. The robustness of the utilized LSM is defined by the physics behind it, but still, it also depends directly on the quality and accuracy of the model input parameters [4,5,6]. While standard vegetation parameters can nowadays be derived reliably from remote sensing data, the availability of high-quality soil data as model input parameters on different scales is often becoming the limiting factor in pedo-hydrological modeling attempts. Mostly, conventional soil maps cannot feasibly satisfy those requirements and are not available in the desired scale or resolution [6,7,8].

However, conventional measurements of soil characteristics (texture, bulk density, moisture) remain time-consuming and non-cost effective and are therefore continuously reduced respectively discontinued [7]. Due to this problem, different strategies have been recently developed to solve model input data uncertainties (Kalman Filter; [9], GLUE; [10], PEST; [11]), while the focus on an efficient retrieval of input parameters, especially soil parameters, is diminishing. Different approaches for retrieving soil parameters have been applied, which can be summarized under the item of pedometrics [7,12,13,14,15]).

Classical geostatistical approaches, based on the work of [16,17], utilize the theory of regionalized variables. In this respect, several soil parameters are connected due to their scale-dependent spatial correlation and allow them to reproduce and transfer spatial patterns of various other soil properties. It is to criticize that these approaches require a large amount of in-situ measurements and data and are only based on mathematical descriptions, while soil functional relationships are not considered. Thus, transferable regionalization models cannot be generated. The integration of available spatial distributed information of secondary parameters derived from digital elevation models (DEM) (such as relief and its derivatives), geological and vegetation data into regionalization models have shown a high potential for the regionalization of soil properties [18]. The so-called ‘CLORPT’ techniques are based on the ‘Factors of Soil Formation’, formulated by [19] according to the formulation of [20].

$$S=f(Cl,O,R,P,T)$$

(1)

where soil (S) is a function of climate (Cl), organisms (O), relief (R), parent rock (P) and time (T). These techniques include simple, multiple and polynomial regressions, generalized linear regionalization models, generalized additive models as well as neural networks. A promising and often adopted approach was introduced by [21], who utilized a hybrid technique of Regression Kriging by incorporating soil-relief and soil-vegetation relationships into the regionalization model. Underlying data are relief and vegetation information derived from high-resolution digital elevation models (DEM) and high resolution optical (Landsat) satellite data. It was found that these parameters are usually well correlated with physical and chemical soil properties. The major advantage of this hybrid model approach is the significant reduction of in-situ soil measurements [13,22,23]. While the incorporation of DEM information is only useful in areas with strong topographic gradients [24], the usefulness of remote sensing information as co-variables in smooth areas is as yet not sufficiently investigated. In various studies, several authors tackle the problem of the missing relief as a soil-forming factor using different remote sensing based co-variables. These can be either yield estimates [24], short term changes in spatial patterns due to rainfall events [25] or general vegetation characteristics [26].

It is understood that the spatial pattern of crop yield or the health status and phenology of vegetation is affected by the underlying soil properties, which govern mechanical stability as well as the water and nutrient retention capacities and exchange functions [27]. Consequently, it must be assumed that these parameters, which can be retrieved accurately from remote sensing data, can be utilized to invert soil property distribution [18,22]. Especially in the field of Precision Agriculture, these methods for the derivation of soil yield potential from remote sensing data are subject to many investigations and have often been applied successfully [26,28,29]. Numerous spectral indices are derived from the ratio of vegetation reflectance in the near-infrared and red spectral range, among which the simple Normalized Difference Vegetation Index (NDVI) is mostly applied. Due to high correlation, the spatial distribution of many essential vegetation parameters such as biomass, crop coverage, leaf area index (LAI) and vegetation height have been retrieved from NDVI products using linear or non-linear regression or other inversion techniques [30,31,32].

The work in [33] combined selected relief parameters with NDVI information retrieved from Landsat-TM data and achieved coefficients of determination from 42% to 78% for the prediction of phosphate, respectively 54% for the prediction of carbon content in different soil depths. The authors of [34] utilized NDVI information derived from Landsat-TM and ASTER data in combination with nine topographical parameters determined from a DEM to regionalize soil properties in a paddy soil area in China. Multiple linear regressions yielded in coefficients of determination for carbon content and silt content between 0.41 and 0.29. In a further attempt, Regression Kriging (RK) was performed and improved the prediction of carbon, nitrogen and silt content by 14% (carbon), 13% (nitrogen) and 10% (silt), as compared to Simple Kriging (SK). Besides, [25] used change patterns (dynamic feedback patterns) of the land surface, such as those captured daily by remote sensing images during a short period (6–7 d) after a significant rain event to differentiate soil types. However, the utilization of optical remote sensing data as co-variables in regionalization models is limited to cloud free acquisition dates which reduces the reliability and availability drastically.

Due to their independence from weather conditions, sensors in the microwave spectral range have been developed and tested in recent years to derive vegetation and soil parameters as well as their spatial distribution. From a remote sensing perspective, a microwave remote sensing system operating as an active system is the most appropriate and straight forward way to monitor earth system dynamics at high spatial resolutions [35]. The principle of an active microwave remote sensing system with synthetic-aperture radar (SAR) is based on the transmission and receiving of an electromagnetic wave at given frequency and polarization (e.g., H = horizontal, V = vertical) by an antenna (aperture) on board a carrier, for instance, satellite or airplane [35,36]. Due to its longer wavelength compared to (classical) systems operating in the optical domain, microwave remote sensing permits retrieval of information from the inner of an imaged media. In particular, over vegetated soil surfaces, depending on the utilized wavelength and object parameters, the backscatter can contain information from the soil surface as well as vegetation structure content. While imaging a bare soil surface, the backscatter is the result of the surface and subsurface conditions of the soil. Thus being well suited for soil properties characterization.

In general, a microwave signal is sensitive to the dielectric properties as well as the geometry and orientation of the imaged object. Concerning soil properties, the dielectric properties are described by the complex dielectric constant and are mainly driven by the soil moisture content as well as the bulk density of a given bare soil surface. In the presence of a vegetation layer above the soil, the dielectric constant is a composite of the vegetation and soil water content at high frequencies. While at lower frequencies, the impact of the vegetation water content on the backscattered signal is not that dominant. On the other hand, the geometry of an imaged bare soil surface is a composite of the distribution of the soil aggregates, particularly for an agricultural field, the composition of soil aggregates, seedbed rows and local topography summarized as soil surface roughness. Thus, a direct retrieval of soil properties is possible [37]. Even though changes in backscatter due to differences in soil texture are small compared to the impact of a change in soil moisture. In [38], an inversion scheme of an LSM model to retrieve soil texture based on spatially distributed soil moisture inputs retrieved from SAR data was used.

In the recent past, [39,40,41,42] used direct empirical relations between SAR backscatter values and soil texture values. All of the studies mentioned above observed a change of backscatter values due to a change of clay content. In particular, [39] demonstrated that for a soil containing more clay, backscatter values could be up to 3 dB weaker. Nevertheless, all of the studies retrieved soil texture over sparsely vegetated or bare soil crops, thus making short wave SAR data applicable, while for a more sophisticated (vegetated) agricultural scenery such an approach is more challenging as it might be biased by vegetation water content as well as differences in soil moisture. At long wavelengths (e.g., 23 cm), information from the sub-canopy layer or sub-soil is retrieved. For instance, [43] proposed an approach for the direct retrieval of soil bulk density over an agricultural test-site by using the spatial distribution of soil surface roughness retrieved from fully polarimetric SAR data as a proxy of bulk density. Nevertheless, an approach using a relative index in changes of backscatter from multi-temporal SAR data would solve the inherent dependency of the backscattered signal to soil moisture and soil surface roughness to allow retrieval of soil texture.

In this paper, we investigate the usefulness of microwave remote sensing data acquired at a longer wavelength for an improved regionalization of soil properties at field scale. Thus, the hypothesizes of the study is defined as the following:

• Microwave remote sensing can improve the prediction of soil texture by means of Regression Kriging compared to classical interpolation methods.

Therefore, soil characteristics in terms of soil texture were measured over agricultural fields in the framework of two different campaigns carried out over a test site in North-East Germany. Besides and simultaneous to the measurements of soil texture, microwave remote sensing data were acquired using an airborne synthetic aperture radar system.

2. Methods

The study was carried out with data sets from two different campaigns acquired in 2006 and 2008 over the Demmin test-site as part of the TERENO site “North-Eastern German Lowland Observatorium”. The campaign was part of the AgriSAR2006 campaign, which was funded by the European Space Agency.

2.1. Study Area

The test-site is located in the young moraine area, Mecklenburg-West Pomerania, North-East Germany, approximately 150 km north of Berlin. The altitudinal range within the test-site is about 60 m. The works in [44], as well as [45], classifies the area of the test-site as a gently undulating ground moraine generated within the Mecklenburger Stadium of the last ice age. Thus, the landscape is characterized by smooth topography, kettle holes as well as small regional meltwater valleys. Within the Mecklenburger Stadium calcareous marly till has been deposited, which was further overlayed with partially thick glacial sands [45,46]. During pedogenesis within the Holocene, those glacial sands have been weathered to sandy loam and loamy sand [47], which are according to [45] the main soil texture in that area. Due to these different parental substrates and the different topographic conditions, several soil types can be differentiated. Within the valleys, Planosols and Gleysols are dominant due to a thin glacial sand cover layer. In contrast, fen and peat can be found within kettle holes. At locations with a glacial sand layer larger than 80 cm Luvisols and Cambisols/Brunic Arenosols [45]. Information from the state soils evaluation (German: Bodenschaetzung) confirm the above-mentioned findings. The soils within the area are highly fertile and productive. As a consequence, they are used intensively for agricultural cultivation. The main crop rotation is winter wheat, winter rape and winter barley. Additionally, maize and sugar beet is sown in spring. A small percentage of forest and sporadic trees along roads and field paths result from intense agricultural cultivation. The mean field size in the area is 225 ha. Due to the very large fields and intensive cultivation, small and local wind, water and tillage erosion evidence can be observed.

2.2. In-Field Data

In-field data were sampled during two single campaigns in 2006 and 2008. In 2006 the focus of the sampling strategy was to sample dynamic soil properties. In 2008 only static soil parameters like texture were sampled. Besides, the sampling design of the campaign in 2008 was focused on a raster-based sampling scheme on two big fields, while in 2006, the sampling design for dynamic soil parameters was based on expert knowledge in terms of maximum differences. Figure 1 shows the sampling locations for both campaigns.

Table 1. Acquired soil data during different campaigns.

Year	n	Soil Texture	Ctot	N	Bulk Density	Soil Moisture	Soil Surface Roughness
2006	27	x			x	x	x
2008	400	x	x	x

summarizes the target variables which were sampled. The overall layout of the sampling design (dynamic/stable soil properties) was chosen due to the limited availability of microwave remote sensing data in 2008. As only for 2006 microwave remote sensing images were available weekly, only dynamic soil parameters were measured in 2006. Thus only static soil parameters (e.g., soil texture) were used for geostatistical analysis, as it can be assumed that a change of soil texture is not very likely within the period of two years. Soil texture was measured by sieving and sedimentation of the mineral soil following DIN ISO 11277:2002-08.

2.3. Remote Sensing Data

Microwave remote sensing data were acquired for a whole vegetation period (April–August) on a nearly weekly basis. Data was acquired by the experimental-SAR (E-SAR—now F-SAR) system operated by the German Aerospace Center [48]. Data were acquired at different frequencies and polarizations. As the penetration depth of a microwave imaging system into a medium is generally dependent on the deployed wavelength, only low-frequency acquisitions were used for this study. In particular, images acquired at 1.3 GHz (L-Band) were used.

2.3.1. Processing of Remote Sensing Data Sets

The used SAR data was processed and delivered by the German Aerospace Center in the framework of the AgriSAR2006 campaign in two different data formats: In terrain and radiometrically corrected backscatter values (GTC) and in raw single look complex image (SLC) formats. While the latter contains information about the phase and the amplitude of the received electromagnetic wave in radar geometry as a complex data type, the terrain and radiometrically corrected backscatter values (GTC) only contain information on the amplitude in integer format. As mentioned earlier, the backscattered signal’s sensitivity is mainly dependent on the dielectric constant (mainly soil and vegetation moisture) as well as geometric properties (surface soil and vegetation roughness) of the imaged surface. Therefore, using a statistical measure that allows to characterize the long term variance of the backscattered signal is essential to map potential differences in soil texture. It is assumed that the long term variations in backscatter contain information on the water retention curve, which is a function of soil texture and water content. To characterize the above mentioned long term variance, the E-SAR GTC data set containing the amplitude information was used for a simple calculation of the coefficient of variation over all consecutive acquisition dates on a weekly basis from April to August. Thus, indicating the normalized standard deviation over pixel and time by the mean over pixel and time as a proxy of the change in soil moisture over time.

$$v_{pp}=\frac{\sqrt{\frac{{\sum }_{i=1}^N{\sigma }_i^0-\overline{x}}{N}}}{\overline{x}}$$

(2)

where $\overline{x}$ is the mean pp-polarized backscatter values (${\sigma }^0$) per pixel for all acquisitions (N) normalizing the standard deviation of ${\sigma }^0$ per pixel i over all acquisitions. ${\sigma }^0$ is defined as:

$${\sigma }^0={\beta }^0·sin\left(\theta \right)$$

(3)

where ${\beta }^0$ is the calibrated pp-polarized backscatter value delivered by the German Aerospace Center in terms of the GTC data. $\theta$ is the local incidence angle of the electromagnetic wave with respect to the surface normal in radians calculated from a DEM.

To reduce the effect of speckle noise in the imagery, a 7 × 7 pixels moving window filter was applied during the calculation of the backscatter values [49].

2.4. Additional Co-Variables

Besides the remote sensing data, different co-variables have been calculated and used for the regression analysis. These co-variables were derived from a digital elevation model (DEM). The DEM was derived from an X-Band airborne interferometric SAR system and was resampled to a spatial resolution of 2 m by 2 m. Out of this DEM, primary and secondary terrain attributes where calculated, using the terrain analysis module of the SAGA GIS software package [50]. Slope ($\beta$), analytical hill-shading (ACN), plan curvature (Ch), profile curvature (Cv), divergence/convergence index (CDI), catchment network base level (CNBL), aspect (A), topographic wetness index (TWI) following [51], SAGA wetness index (SWI) [52], length-slope (LS) factor of the USLE and the stream power index (Ω) (following [53]). Detailed information regarding the terrain analysis procedure of SAGA GIS can be found in [14,54,55].

Each calculated variable represents terrain attributes that can be linked to the intensity of natural processes that influence the soil-forming processes. A good explanation of these variables and their physical meanings can be obtained from [6]. As the DEM showed significant inconsistencies in terms of abrupt changes, the DEM was low-pass filtered with a filter size of 25 by 25 pixels. Besides, all co-variables were cropped in its extension to the agricultural fields where the majority of samples have been taken from.

2.5. Geostatistical Approaches

Geostatistical techniques are commonly used tools in soil mapping and have been extensively used to study patterns of environmental variables [56,57,58]. A geostatistical approach was carried out in the applied study based on variogram and regression analysis between the target variables and different co-variables. For the co-variables different spatial grid-based information layers were used as explained in Section 2.3 and Section 2.4. In this study, a hybrid regionalization approach was used, widely known as Regression Kriging C, as defined by [21] and from hereon referred to as Regression Kriging (RK). In the literature, several other synonyms exist, such as kriging with external drift. However, from a mathematical point of view, they are all the same [14].

RK is one of the most popular, practical and robust hybrid spatial interpolation techniques in the digital soil mapper’s toolbox and is defined as the sum of the regressed values and kriged residuals from the regression [59].

The deterministic part of the RK model was calculated using a multiple-linear regression (MLR) as follows [15]:

$$\widehat{Z}\left(S_0\right)=\sum _{k=0}^p{\widehat{\beta }}_k·q_k\left(S_0\right);q_0\left(S_0\right)\equiv 1,$$

(4)

where $\widehat{Z}\left(S_0\right)$ is the predicted soil property, $q_k\left(S_0\right)$ are the values of the auxiliary variables at the target location, ${\widehat{\beta }}_k$ are the estimated regression coefficients and p is the number of regressors.

The Kriging of the residuals was done using an Ordinary Kriging approach. Predictions using Ordinary Kriging (OK) are commonly made by calculating some weigthed average of the observations [15,58]. With this the predicted value $\widehat{Z}\left(S_0\right)$ of the target variable at unvisited locations $S_0$ is calculated as follows:

$$\widehat{Z}\left(S_0\right)=\sum _{i=0}^n{\lambda }_i·z\left(s_i\right)$$

(5)

where $z\left(s_i\right),\dots ,z\left(s_n\right)$, is the sample data with coordinates. The weights ${\lambda }_i$ are optimized in a way that the prediction error is zero, the variance of the prediction error is minimized and the sum of weights is one [6,15]. The estimation of the weights ${\lambda }_i$ strongly depends on the variogram $\widehat{\gamma }$. Empirical variograms have been calculated using the standard method-of- moments estimator according to [60]:

$$\widehat{\gamma }=\frac{1}{2}\frac{1}{N\left(h\right)}\sum _{i=1}^{N\left(h\right)}{[Z(S_i+h)-Z\left(S_i\right)]}^2$$

(6)

Based on the Equations (4) and (5) the hybrid model is defined as follows:

$${\widehat{Z}}_{OK}\left(S_0\right)=\sum _{k=0}^p{\widehat{\beta }}_k·q_k\left(S_0\right)+\sum _{i=1}^n{\lambda }_i·ϵ\left(S_i\right)$$

(7)

where $\widehat{m}\left(S_O\right)$ is the fitted deterministic part, $\widehat{ϵ}\left(S_O\right)$ is the interpolated residual, ${\beta }_k$ are estimated deterministic model coefficients (${\widehat{\beta }}_k$ is the estimated intercept), ${\lambda }_i$ are kriging weights determined by the spatial dependence structure of the residual and where $ϵ\left(S_i\right)$ is the residual at locations $S_i$ [14].

The regression coefficients ${\beta }_k$ were estimated by applying the extracted pixel values of the co-variable grids Section 2.3 and Section 2.4 at sampling locations for each target variable into the MLR model and using the ordinary least squares (OLS) method to derive the best fit. A backward stepwise procedure (of the R-library “MASS” [61]), was used to reduce the number of co-variables and get the final regression equation. Only those auxiliary variables that were significantly correlated with the particular target variable regarding Pearson’s correlation coefficients and avoiding overfitting those without multi-collinearity were considered as initial regressors. The final MLR equations were applied on the co-variable grids to produce the grid $Zpr\ast$ (deterministic part of RK) for every target variable. The MLR residuals of the final regression models were applied into a variogram analysis (Equation (6)) and finally interpolated using the OK interpolation. With this procedure, a grid of the MLR residual error ($ϵ\ast$) was interpolated for every target variable. The final RK prediction map was calculated by the sum of the $Zpr\ast$ and $ϵ\ast$ grid for every soil texture.

$$Soi{l}_{map}=Zpr\ast +ϵ\ast$$

(8)

All geostatistical modelling was performed using the statistical software environment of R [62].

2.6. Validation of the Model Performance

The prediction performance of applied interpolation models was tested against the sampled soil physical properties. The model performance was assessed using a leave-one-out cross-validation (LOOCV) [63]. The Root Mean Square Error (RMSE) and the mean absolute error (MAE) are regularly employed in model evaluation studies. However, if the error distribution is expected to be Gaussian than RMSE is more appropriate to represent model performance than the MAE [64]. Due to that, the accuracy of each model run was compared with the Root Mean Square Error:

$$RMSE=\sqrt{\frac{1}{n}{\Sigma }_{i=1}^n({\widehat{z}}_i-z_i{)}^2}$$

(9)

with ${\widehat{z}}_i$ and $z_i$ being the predicted and measured values, respectively.

Results of the RK approach were also compared against Inverse-Distance-Weighted (IDW) interpolation and the basic Ordinary-Kriging (OK) approach. For both OK and IDW approaches, readers are requested to refer to standard literature such as [58].

3. Results

Figure 2, Figure 3, Figure 4 and Figure 5 show as an example several of the calculated co-variables for the site. While the DEM and the derived topographic indices reflect the obvious topographical field situation with smooth terrain, the coefficient of variation of the backscatter values reveals several additional features that do not follow the topographic features in that area. These are in particular areas with either different soil texture or soil aggregates, respectively, porosity [43]. The soil sampling locations are shown in Figure 1 and the lab-analyzed results summarized in

Table 2. Descriptive statistics of soil samples.

Texture Class	Mean [Mass-%]	Min [Mass-%]	Max [Mass-%]	std [Mass-%]	cv
sand	63	44	76	5.1	8.1
silt	27	16	36	3.6	13.4
clay	9	3	26	4.1	41.5

. The sand fraction usually dominates the texture composition of northern German topsoil. In the soil samples of the test site, sand is also the dominating fraction with a mean of 63% and a range between 44% and 76%. The clay content shows lowest values with a mean of 9% and a range between 3% to 26% and the silt fraction reaches values between 16% and 36% and a mean content of 27%. The coefficient of variation (cv), on the other hand, shows a reverse picture. The clay fraction of the samples show the highest cv with 41.5, the silt fraction shows a lower cv = 13.4 and sand has the lowest cv = 8.1. In Figure 6, soil samples are classified based on the clay and silt content following the German KA5 and the USDA classification system.

Table 3. Fitting results of the deterministic term of Equation (7) in terms of regression equations. Meanings of the abbreviations are explained in Section 2.4 and Section 2.3.1.

Texture Class	Best Fit MLR-Model	RMSE
Sand	−041 ∗ DEM + 4179 ∗ $\beta$ + 0.88 ∗ CNBL + 85.32	3.24
Silt	0.37 ∗ DEM + 8.92 ∗ LS − 2.70 ∗ $v_{HH}$ − 190.75 ∗ $\beta$ − 1.14 ∗ CNBL − 0.48 ∗ TWI + 14.92	2.65
Clay	6.23 ∗ ACN − 11.12 ∗ LS + 155.47 ∗ $\beta$ + 2.43 ∗ $v_{HH}$ + 0.4 ∗ TWI − 1.49	2.19

summarizes the outcome of the multiple-linear and stepwise regression approach as well as the resulting RMSE of the models. To select the best co-variables in terms of target variable prediction, the backward stepwise procedure of the R-package "stats" has been applied. The “stats” package uses the “AICstep” algorithm of the R-package “MASS” for the selection of the final model without any over- or under-fittings [61]. As can be seen, different models have been obtained for certain texture classes based on the stepwise regression approach. For all the texture classes, the slope ($\beta$) is an important estimator in all models, DEM, LS, CNBL, TWI and $v_{HH}$ are only relevant in two of the three texture classes. Having different deterministic models for single texture classes is also observed by other studies. Besides, the higher sensitivity of the HH polarization $v_{HH}$ compared to other polarization like VV or VH is also in good accordance with other studies. For instance, [65] reports a higher sensitivity of HH to soil properties than VV or VH, especially under vegetation. It is also to notice that the models for silt and clay both contain TWI. This leads to a certain pattern visible in both resulting maps but not in the one for sand.

Figure 7 shows the sample variogram and the fitted spherical variograms for the residual term of Equation (7). The best fit was reached by using a spherical function to describe the sample variogram and to estimate the weights of ${\lambda }_i$ in Equation (7). Using these weights in conjunction with the best MLR estimates, the spatial distribution of the three soil texture classes (sand, silt, clay) have been estimated. Figure 8, Figure 9 and Figure 10 show the outcome of the regionalized soil texture classes by using the RK approach. In addition to the interpolated texture values, the fraction of the corresponding texture class of the sampling campaign is also displayed with dots given in the same color range. By doing so, mismatches of the interpolated values compared to the sampled values become visible. As can be seen, the majority of the different sample points are well assessed by the used approach. However, several mismatched locations exist. These are, in general, the areas along the valley of the test-site. As can been seen, there is an evident shift from lower silt values in the samples to higher silt values in the interpolated values. The same, but vice-versa can be observed for the texture class sand. Consequently, the amount of sand in such areas is significantly underestimated by the used approach. Besides, a much smoother appearance of the interpolated sand map is visible (Figure 8). This missing small scale heterogeneity can be related to the absence of the TWI and LS-factor in the MLR model for sand. Both co-variables add significant structures to the interpolated results for clay and silt. However, the general performance of the used approach is very high.

Table 4. Root mean square error (RMSE) [Mass-%] for retrieved soil texture classes using the different approaches compared to in-situ data. RK: Regression Kriging, OK: Ordinary Kriging, IDW: Inverse Distance Weighted, * denotes a significant difference of the soil texture distributions from the different interpolation methods at 0.999 confidence level according to a Kolmogorov–Smirnov (K–S) test.

Texture Class	RK	OK	IDW
Sand	2.85 *	2.86 *	2.93 *
Silt	2.31 *	2.42 *	2.54 *
Clay	2.12 *	2.22 *	2.26 *
mean	2.43	2.51	2.57

shows the accuracy in terms of the RMSE and the performance in comparison to standard deterministic and geostatistical approaches such as IDW and ordinary-kriging (OK). With a mean RMSE over all texture classes of 2.4 mass.-% the approach is very accurate and outperforms the other approaches. Even though the performance of the other approaches is only slightly worse with slight advantages for OK. However, a two-sample Kolmogorov–Smirnov (K–S) test revealed a significant statistical difference between all the derived approaches with a level of confidence of 0.999.

Figure 11, Figure 12 and Figure 13 show the outcome of the regionalized soil texture classes by using the OK approach. In general, a much smoother appearance, especially for the silt and clay texture classes can be observed, as there are no co-variables involved in the interpolation process. In addition, due to the data and variogram driven interpolation range, a more natural pattern, compared the IDW approach (Figure 14, Figure 15 and Figure 16) is given. However, the variability along the slopes of the meltwater valley is not reproduced.

Figure 14, Figure 15 and Figure 16 show the outcome of the regionalized soil texture classes by using the IDW approach. In general, an IDW inherent bulls-eye appearance can be observed for all three texture classes. Even by changing the distance weighting factor, the local extreme values do not disappear or result in a strong smoothing of the interpolated values. Thus, a very unnatural pattern is obtained by the approach.

4. Discussion

In this paper, we presented an approach that combines geostatistics and microwave remote sensing based analysis to predict soil texture at field scale. Here, the remote sensing-based information is obtained by airborne microwave remote sensing SAR data in terms of the coefficient of variation of an image time series reflecting the long-term dynamics over the agricultural scene. It serves as a co-variable in a Regression Kriging approach as defined by [21]. In particular, type ‘C’ of the RK mappers toolbox is used [21]. Additional co-variables are used in terms of topographic indices such as slope, LS-factor, TWI, Analytical Hillshade and other to model soil genesis based on topographical features. In a statistical analysis, three independent deterministic models have been developed to model the single texture classes using the co-variables and Regression Kriging approaches. While in general, the co-variables based on the DEM reflect the dependency of the different texture classes from topography, the remote sensing based information is used to reflect the spatial distribution of texture independent from topography.

In the literature, authors used optical remote sensing data sets [24,25,26,34], which have clear disadvantages in cloud prone areas. Here microwave remote sensing has a clear advantage, as it is able to see through clouds. Besides, it also adds further information from inside the soil, as it can penetrate the soil due to its longer wavelength. In this study, a higher sensitivity of the HH polarization of the SAR data was observed compared to the VV polarization. Similar observations were made by other authors using, for instance, microwave remote sensing for the retrieval of soil moisture [65,66]. It is also to note that the HH polarized microwave data is only included in the MLR-models for silt and clay. Here, a good dependency of the backscatter signals due to the better water retention in those texture classes can be observed. In addition, the Slope and TWI were also observed as good descriptors of the spatial distribution of soil texture classes. Similar observations had been made by [6].

The general performance of the applied approach is excellent. With a mean RMSE of 2.43 mass-%, the RK approach slightly outperforms the other methods. This might have several reasons. First, the variability of the different texture classes is not very high, showing a homogeneity over the study site with only slight variations. This results in the comparable performance of the deterministic and geostatistical approaches compared to the used hybrid RK approach. Second, the study site’s extent is not very large, with a large number of sample points resulting in a dense sampling pattern, which is advantageous for the deterministic IDW approach. However, as revealed by a two-sample K-S test, differences in the retrieved soil texture maps for the different approaches are significant. Consequently, the RK approach can be considered very suitable for the retrieval of soil texture by incorporating microwave backscatter information as a co-variable.

Having a closer look at the retrieved soil texture patterns, the RK approach results show much more spatial details and heterogeneity than the OK and IDW approach. This is especially true for the clay and silt maps and results from the co-variables. Besides, the IDW approach results in a typical bulls-eye pattern, which is very common for IDW and looks very unnatural. The OK approach results in a much smoother appearance compared to other approaches. However, the local topographic effects in soil texture distributions are not well reflected.

It is notable that the smallest RMSE values are obtained for the clay fraction. This is also in good accordance with the studies of [40,41], which retrieved the highest sensitivity of backscatter values to clay content. The mean RMSE over all texture classes of the RK approach is also in good accordance with other studies. The authors of [40] retrieved soil texture with an RMSE of 1.2 mass-%, while [67] stated, that the overall accuracy of the global SoilGrids product is at 20–30% relative error. With respect to the German classification scheme (see Figure 6), the smallest texture class has a range of 5%. Thus, a target accuracy below 5 % is envisaged. With the retrieved RMSE value of 2.45, the proposed approach nicely fits this target accuracy.

As can been seen from Figure 9, there is an evident shift in from lower silt values in the samples to higher silt fraction in the interpolated values resulting in an underestimation of the sand fraction. This is especially evident in the areas of small scale slopes in the developed meltwater valleys. This shift might be due to the weaker performance of the approach in the area of gentle slopes or the inherent higher assumed soil moisture values in these areas. A higher soil moisture values would consequently lead to a higher backscatter value of the microwave acquisition, thus leading to a higher fraction of silt in that area.

5. Conclusions

In this study, an approach was shown to assess spatial soil texture estimates using microwave remote sensing and topographic indices jointly used in a hybrid geostatistical approach. In particular, a Regression Kriging (RK) approach was used with airborne L-Band microwave backscatter values and topographic indices derived from an airborne digital elevation model. With a performance of a mean RMSE of 2.45 mass-%, the approach shows good potential for future surface soil texture assessments even at high resolution or for global applications. This is especially true as the approach is applicable over any area, as the used microwave data is independent of clouds or illumination problems. In particular, with the upcoming remote sensing sensors and acquisition plans, the usability of such an approach increases significantly. With the availability of the Copernicus Sentinel mission [68,69], a continuous observation with free accessibility of the land surface with a SAR system is real. However, to overcome the shortcoming of a C-Band SAR system in terms of vegetation–soil surface interaction [35,70], longer wavelengths such as the deployed L-Band data are preferred. Here the already existing ALOS2 and SAOCOM, as well as the upcoming NiSAR or TanDEM-L missions, show high potential for future usage of the developed approach over complex terrain such as agricultural fields. With there acquisition plans, spatial and temporal high-resolution imaging is possible, allowing for the creation of a database for the prediction of soil texture in the proposed framework. Here the ALOS2 and upcoming TanDEM-L missions have to be highlighted, as they have the same imaging characteristics as the deployed E-SAR data.