# Spatially Explicit Estimation of Clay and Organic Carbon Content in Agricultural Soils Using Multi-Annual Imaging Spectroscopy Data

Authors: Heike Gerighausen; Gunter Menz; Hermann Kaufmann
Journal: Applied and Environmental Soil Science
 (2012)

---

## Abstract

Information on soil clay and organic carbon content on a regional to local scale is vital for a multitude of reasons such as soil conservation, precision agriculture, and possibly also in the context of global environmental change. The objective of this study was to evaluate the potential of multi-annual hyperspectral images acquired with the HyMap sensor (450–2480 nm) during three flight campaigns in 2004, 2005, and 2008 for the prediction of clay and organic carbon content on croplands by means of partial least squares regression (PLSR). Supplementary, laboratory reflectance measurements were acquired under standardized conditions. Laboratory spectroscopy yielded prediction errors between 19.48 and 35.55 g kg−1 for clay and 1.92 and 2.46 g kg−1 for organic carbon. Estimation errors with HyMap image spectra ranged from 15.99 to 23.39 g kg−1 for clay and 1.61 to 2.13 g kg−1 for organic carbon. A comparison of parameter predictions from different years confirmed the predictive ability of the models. BRDF effects increased model errors in the overlap of neighboring flight strips up to 3 times, but an appropriated preprocessing method can mitigate these negative influences. Using multi-annual image data, soil parameter maps could be successively complemented. They are exemplarily shown providing field specific information on prediction accuracy and image data source.


---

## Body

## 1. Introduction

Soil is the uppermost, weathered layer of the earth’s crust which forms the interface of the lithosphere, biosphere, hydrosphere, and atmosphere. As such, it acts as one of the major resources available to man whose conservation must be given high priority [1]. Clay and organic carbon are two key soil attributes as they contribute many benefits to its physical, chemical, and biological properties such as soil structure, soil water holding capacity, and soil fertility. Furthermore, the question whether CO2 can be sequestered in agricultural soils, to what degree and under which circumstances [2], has dramatically increased the interest in estimating soil organic carbon stocks (e.g., [3–5]). Yet, monitoring soil nutrient supply, degradation status, or organic carbon stocks requires information on soil properties over large areas with a high spatial resolution. The rate of repetition will depend on the issue considered. In any case, conventional soil sampling strategies suffer from the lack of providing data with a high spatial and temporal resolution as they are costly, labor intensive, and time consuming.As an alternative to standard analytical methods, visible and near infrared reflectance spectroscopy (VNIRS) has been exploited successfully for the quantitative analysis of various soil properties in the laboratory since many years (e.g., [6–11]). Because of the spectral characteristics of clay minerals [12, 13] and organic carbon [14], there is a good chance for their prediction from reflectance spectra with a reasonable to good accuracy. Yet, interactions with other soil properties such as iron or sand can diminish the prediction accuracy [15].Imaging spectroscopy measures surface reflectance with a given spatial resolution on a regional scale instead of a single point, thus offering an opportunity for spatially explicit estimations of soil characteristics. Early studies on soil clay and organic carbon content were published by Krüger et al. [16], Chabrillat et al. [17], and Ben-Dor et al. [18]. Meantime, other studies have proven the potential of imaging spectroscopy to retrieve quantitative information on soil clay or organic carbon content (e.g., [19–21]). But the number of studies is still limited and many problems remain unsolved including radiometric calibration of the sensor and atmospheric effects [22], variable water content [23], scattering effects of wide field of view sensors [20], or vegetation [24, 25]. The latter is one of the main factors hampering the estimation of soil properties on agricultural fields as arable land is covered by crops most time of the year due to cultivation practices. Siegal and Goetz [24] reported that already 10% of green vegetation can mask mineral absorption characteristics beyond recognition. According to their results, dry vegetation does not change spectral characteristics of minerals but only affects albedo. Murphy and Wadge [26] concluded the reverse for the shortwave infrared region in a study of different soil types in arid terrain. These findings were approved by Murphy [27] who investigated the impact of vegetation on soil mineral absorption features near 2200 nm. Bartholomeus et al. [28] and Ouerghemmi et al. [29] stated that prediction errors for clay and organic carbon content decreased progressively as vegetation cover exceeded 5 to 10%. Different strategies have been developed to account for the negative impact of vegetation cover on soil parameter estimations. Bartholomeus et al. [30] showed that a combination of several spectral indices could reduce the estimation errors for soil iron content in partially vegetated areas in their specific case. Rodger and Cudahy [31] successfully introduced a vegetation corrected continuum depth (VCCD) method to remove, or negate, the obscuring effect of the vegetation. But research by Lagacherie et al. [22] and Gomez et al. [32] disclosed that continuum removal analysis is clearly inferior to the PLSR approach for VNIRS quantitative analysis. Apart from that, VCCD only applies to the absorption band at 2200 nm typical for Al-OH bearing minerals such as clay but cannot be adopted for the prediction of organic carbon content. A spectral unmixing based approach was utilized by Bartholomeus et al. [28] to eliminate the influence of vegetation from the spectral reflectance of mixed pixels over partially covered maize fields to estimate the soil organic carbon concentrations. The success of their approach essentially depends on the quality of the unmixing results which may deteriorate if the spectral variability of the studied area increases. Ouerghemmi et al. [29] explored a “blind source separation” (BSS) algorithm to isolate the soil signal from mixed surfaces of vegetation and soil. They yielded accurate predictions with this procedure in 40% of their test cases. Problems were encountered if vegetation cover exceeded 50%. Furthermore, they subsumed that a better understanding of the behaviour of the BSS algorithm is required to improve the selection process for the extracted soil signals. A completely different approach was accomplished by Kooistra et al. [33] who used grass reflectance as a proxy for the estimation of soil metal concentrations. However, the most commonly applied strategy, and probably the most reliable in terms of prediction accuracy, is the restriction of parameter estimates to bare soils (e.g., [22, 32, 34]). The major disadvantage of this approach for croplands is the loss of a noticeable percentage of the cultivated area for parameter estimation. This may be compensated if several images of one study area acquired in different years are analyzed. Parameter estimation would hence no longer be restricted to a few fields free of vegetation but permits a gradual completion of soil maps of cropland on a regional scale. Applying such an approach, particular attention has to be paid to the temporal variability of the soil parameters of interest. While the clay content or soil particle size distribution does not change a lot within a single decade [35] except for transport and deposition of large amounts of soil particles during major erosion events, soil organic carbon is subject to both seasonal and interannual variation mainly controlled by temperature and soil moisture (e.g., [36–38]). Leinweber et al. [36] reported pronounced seasonal variations in topsoil organic carbon concentrations on arable land between 2.40 and 4.30 g kg−1, measuring weekly from April to September. Most distinctive was a decline between June and August. Rogasik et al. [39] revealed variations between 0.40 g kg−1 and 1.15 g kg−1 organic carbon though their measurements were carried out monthly from November to April. Yearly measurements of topsoil organic carbon contents over more than 20 years suppose that the year-to-year variability is of the same order than seasonal variations and range between 1.94 g kg−1 and 4.40 g kg−1 [40, 41]. The absolute amount of change varies with treatment, tillage intensity, and type of crop litter.The aim of this study is to investigate the potential of imaging spectroscopy for spatially explicit estimations of two important soil constituents, clay and organic carbon, on agricultural fields by evaluating multi-annual hyperspectral image data using multivariate calibration. This has, to our knowledge, not been done so far. Three image data sets of the Australian sensor HyMap, obtained during HyEurope campaigns in 2004, 2005, and 2008, were explored by means of partial least squares regression (PLSR) and reference data collected during a field campaign in 2006. The paper addresses in particular (i) the gain in percentage area with soil parameter estimates, (ii) the verification of model accuracies by comparing predictions from different years and in the overlap of adjacent flight strips additional to the common practice of validation with mostly a limited number of reference samples, and (iii) the potential impact of BRDF effects on model predictions and possibilities to mitigate these effects. Prior to image data analysis, laboratory spectra measured under standardized conditions were investigated to assure reliable relationships between chemical and spectral properties. This not only allows the direct comparison between the laboratory data and the image data but also between the laboratory models and the image models. Furthermore, knowledge could be gained on the likely predictive ability and model behavior outside the calibration data range.


## 2. Materials and Methods

### 2.1. Study Site

The study area is the test site of DEMMIN (Durable Environmental Multidisciplinary Monitoring Information Network) located in the lowlands of northeastern Germany (upper left corner: 54°2′54.29′′N, 12°52′17.98′′E, lower right corner: 53°45′40.42′′N, 13°27′49.45′′E). The test site is an intensively used agricultural ecosystem mainly grown with winter grains (about 60%) and root crops (about 13%). About 25% are cultivated as grassland and pasture. The topography is rather flat in the north but hilly to undulate in the south with an altitudinal range of 120 m. Considerable differences in the parent material and relief caused a high spatial variability of soil types in some parts. The flat and slightly undulating plains are characterized by sand rich regions and extensive areas with glacial till. The sand rich regions are dominated by Cambisols, Luvisols, and Albeluvisols. Luvisols, Albeluvisols, and Stagnosols evolved on glacial tills poor in sand but abound in loam and clay. Soils on glacial till in hilly terrain are often truncated and colluvial sediments have accumulated on the bottom of slopes. The floodplains which are mainly used as grassland are characterized by peaty soils [42]. The long-term mean temperature of the area is about 8°C (281 K) with an annual precipitation from 500 to 600 mm increasing from south-east to north-west [43]. The research herein focuses on the eastern part of the test site due to the coverage of the hyperspectral image data (Figure 1).Figure 1
Location of the DEMMIN test site and its agricultural fields in the northeast of Germany near the city of Demmin. The black polygon outlines the maximum extent of the HyMap data acquired in the years 2004, 2005, and 2008. Sampling points are pictured as black circles. The red flags indicate the location of the two fields shown in Figure13 (eastern flag) and Figure 14 (western flag).


### 2.2. Soil Sampling and Analysis

135 soil samples were collected from the agricultural fields in DEMMIN during a field campaign in September 2006 (Figure1). To obtain a representative set of soil samples of the area, sampling sites were preselected with the help of existing soil maps and a digital terrain model. These were then reconciled with the hyperspectral image data available at this time (2004, 2005) to assure the occurrence of vegetation-free fields with reference samples in these images. Finally, they were adapted to field conditions if appropriate. In total, 14 fields were sampled which stretched across approximately 22 km in north-south direction and 8 km in east-west direction. 8 fields were intensively sampled (6 to 43 points per field). From 6 further fields, samples from 2 to 5 locations were taken. The mean distance between adjacent sampling points within the fields was about 138 m reaching 743 m at maximum. The described spatial distribution of the samples was mainly determined by the pedological setting of the test site which is characterized by extensive areas of sandy and loamy soils, and distinct but smaller regions of loamy and heavy clay soils. In conjunction with the flight windows, soils abound in clay and organic carbon could only be encompassed in the reference data by sampling single fields.Soil samples were taken from the top 5 cm of soil within a radius of 2 m. The geographical coordinates of the sampling points were measured using differential GPS with an accuracy of less than 1 m. At 18 sampling points, five samples were taken within the radius of 2 m and analyzed in the laboratory to assess the error of the applied sampling strategy and to get an idea of the possible impacts of subsampling [44]. We estimated this standard error of sampling (SES) using the following equation from Lozán and Kausch [45]:
(1)SES=∑i=1m(∑j=1kxij2-(∑j=1kxij)2/k)n-m,
where x is the measured value, k is the number of samples taken at each point, m is the number of sampling points, and n is the total number of samples analyzed.Particle size distribution [g kg−1], carbon content [g kg−1], total iron [g kg−1], and soil acidity (pH) were analyzed by the Central Laboratory of the Centre for Agricultural Landscape Research (ZALF) using air-dried samples. Particle size distribution was determined using a combined sieve (fraction 0.063–2 mm) and pipette (fraction <0.063 mm) analysis. The total amount of carbon (TC) was determined by dry combustion using an element analyzer (CNS2000 from LECO). Inorganic carbon content (IC) was analyzed measuring the change in volume during decomposition with phosphoric acid using a Carmhomat 12D. Organic carbon content (OC) was then calculated as the difference between TC and IC. Iron contents were determined by dissolution with aqua regia. Soil acidity was measured by radiometer analytics with a TitraMaster85. The standard error of laboratory (SEL) for soil clay content was 5.0 g kg−1. Analytical techniques for the determination of TC and IC had a SEL of 0.17 g kg−1 and 2.88 g kg−1, respectively.


### 2.3. Spectral Measurements

#### 2.3.1. Laboratory Spectroscopy

Spectral measurements were acquired from all soil samples using an ASD spectroradiometer FieldSpec Pro covering the full spectral range from 350 to 2500 nm. The ASD has a fixed spectral sampling interval of 1.4 nm in the visible and near infrared portion of the spectrum (350–1000 nm) and 2 nm in the shortwave infrared (1000–2500 nm) with a full-width half maximum (FWHM) of 3 nm and 10 nm, respectively. All measurements were taken with bare fibre and without foreoptics (25° FOV) relative to a standardized spectralon panel. A 1000 W tungsten lamp was used as light source to illuminate the target at an incident angle of 30°. Before the measurements, all soil samples were grounded with a mortar and passed through a 2 mm sieve to reduce anisotropic light scattering. To standardize the moisture level, they were oven-dried for 48 h at 105°C (378 K).A spectrum of a soil probe was averaged out of ten measurements with a system integration time of 100 scans per single collected spectrum. Noise removal and smoothing was accomplished using the Savitzky Golay filter [46] implemented in IDL [47].


#### 2.3.2. Imaging Spectroscopy

Hyperspectral image data from the test site DEMMIN was acquired by the Australian sensor HyMap mounted on a Dornier Do 228 aircraft and a Cessna 208 B during three HyEurope flight campaigns. The HyMap airborne imaging spectrometer records reflected radiances in 128 bands in the wavelength range from 450 to 2480 nm with a spectral bandwidth (FWHM) of 15-16 nm between 450 and 1900 nm and 18–20 nm between 1950 and 2480 nm. The sensor’s field of view (FOV) is 60° with 512 pixel in across track direction. The instantaneous field of view (IFOV) is 2.5 mrad along track and 2.0 mrad across track. Overflights took place on August 9, 2004, May 27, 2005, and July 29, 2008 at an average flight altitude of about 2000 m resulting in a spatial resolution of4×4 m² (Table 1). At each of the campaigns, seven to ten strips of about 16 km length were acquired under clear-sky conditions covering an area of approximately 200 km² and 260 km², respectively.Table 1
Characteristics of available HyMap data sets and its acquisition parameters.
DoA
ToA (UTC)*
Number of strips
Solar angles (°)
Heading** (°) (±2°)
MFA MSL (m)
PS  (m)
August 09, 2004
13:30–14:20
7
Azimuth: 117.2–130Zenith: 46.3–52.2
0/180
~1990
4
May 27, 2005
9:50–10:50
7
Azimuth: 152.7–178.3Zenith: 31.7–34.0
1/180
~1940
4
July 29, 2008
9:50–11:20
10
Azimuth: 147.1–183.1Zenith: 35.3–38.7
0/180
~1934
4
*Local time is given by UTC+2. **Neighboring strips have the opposite flight direction. DoA: day of acquisition, ToA: time of acquisition, MFA ASL: mean flight altitude above mean sea level, and PS: pixel size.Radiometric calibration of the HyMap data was performed inflight using reference ground measurements [48]. Atmospheric correction and conversion of at-sensor radiances in reflectance values was accomplished using the ATCOR-4 software [48, 49]. In a final step, the data was geometrically corrected by means of a generic ortho image processor [50].For each of the sampling points in the field, HyMap spectra were extracted using a window of3×3 pixels centered on the pixel with the coordinates closest to the plot location. The average of these nine spectra was then used for further statistical analysis.


### 2.4. Spectral Preprocessing

Laboratory spectra were resampled to fit the spectral resolution of the HyMap sensor. Due to the annual recalibration of the HyMap sensor by HyVista, minor differences in the centre wavelengths and bandwidth (FWHM) occur in between years. To obtain spectral consistency, laboratory spectra and image spectra from the 2004 and 2008 flight campaigns were adjusted to the configuration of the HyMap sensor in 2005. Several spectral transformations were applied. Absorbance was computed by taking the logarithm to the base 10 of the inverse reflectance. First and second derivative of absorbance were computed after the Savitzky Golay method with a fifth-order polynomial including 5 data points either side of each point in the filter. Bands prone to noise were eliminated before statistical analysis. Wavelengths removed were located either side of the water absorption bands near 1400 nm and 1900 nm and below 500 nm and above 2400 nm.


### 2.5. Statistical Analysis

A multitude of statistical approaches comprising principal component analysis, multiple linear regression, artificial neural networks, or regression trees, has been utilized for the quantification of soil parameters from spectroscopic data with varying success (e.g., [10, 19, 51–53]). Partial least squares regression (PLSR) has been exploited by many authors and proofed a good efficiency in quantitative spectroscopy ([11], Table 1). PLSR sets up a linear regression model by projecting the predictor variables onto a few so-called latent variables or factors taking into account both, the response and predictor variation. This way, two problems of conventional multiple linear regression modeling are handled effectively: multicollinearity and variable selection. However, choosing the optimum number of factors is a crucial point in PLS modeling. They were determined applying a leave-one-out cross-validation approach to the calibration data set, choosing the model with the lowest Akaike information criterion (AIC). Rather than just using the RMSE as a selection criterion, AIC also integrates the number of samples and the number of PLS factors avoiding a model overfit more effectively. Statistical analysis of the spectra was carried out using the ParLes Software [54] which implements the orthogonal PLS regression for one y-variable after Martens and Næs [55]. The datasets were randomly split into a calibration set and a validation set. About two third was used for calibration purposes. The rest of the samples were retained from model calibration and used as an independent validation dataset. If the total number of available samples was less than 30, no separate validation was conducted but model accuracy was solely assessed via cross-validation. Calibration and validation data schemes are described in detail in Section 3.5.1.The performance of the models was evaluated by calculating the coefficient of determination (R2), the standard error of calibration (SEC), and the standard error of validation (SEP). The latter is given by the following equation:
(2)SEP=1n-1·∑i=1n(y^i-yi)2,
where y^i and yi are the predicted and observed values of the sample i and n is the total number of samples in the validation data set. SEC is calculated the same way but additionally takes into account the number of PLS factors used in the model. As an independent measure of model accuracy, the ratio of performance to deviation (RPD) was computed by dividing the standard deviation of the reference samples through the SEP [56]. If no independent validation is performed, information on the predictive ability of the model can be retrieved by calculating the ratio of the SEC and the standard deviation (SD) of the calibration data set (SD/SEC). For each model variable, importance in the projection (VIP) and PLS-regression coefficients (b-coefficients) were studied to identify important wavelengths used during model calibration and evaluate their selection with respect to known absorption features of the parameter of interest.To compare parameter predictions between years, root mean squared error was calculated. It is given by the equation(3)RMSE=1n·∑k=1n(y^k,A(ij)-y^k,B(ij))2.
where y^k,A(ij) and y^k,B(ij) are the predicted values of the pixel ij in data set A and in data set B, and n is the total number of pixels.


## 2.1. Study Site

The study area is the test site of DEMMIN (Durable Environmental Multidisciplinary Monitoring Information Network) located in the lowlands of northeastern Germany (upper left corner: 54°2′54.29′′N, 12°52′17.98′′E, lower right corner: 53°45′40.42′′N, 13°27′49.45′′E). The test site is an intensively used agricultural ecosystem mainly grown with winter grains (about 60%) and root crops (about 13%). About 25% are cultivated as grassland and pasture. The topography is rather flat in the north but hilly to undulate in the south with an altitudinal range of 120 m. Considerable differences in the parent material and relief caused a high spatial variability of soil types in some parts. The flat and slightly undulating plains are characterized by sand rich regions and extensive areas with glacial till. The sand rich regions are dominated by Cambisols, Luvisols, and Albeluvisols. Luvisols, Albeluvisols, and Stagnosols evolved on glacial tills poor in sand but abound in loam and clay. Soils on glacial till in hilly terrain are often truncated and colluvial sediments have accumulated on the bottom of slopes. The floodplains which are mainly used as grassland are characterized by peaty soils [42]. The long-term mean temperature of the area is about 8°C (281 K) with an annual precipitation from 500 to 600 mm increasing from south-east to north-west [43]. The research herein focuses on the eastern part of the test site due to the coverage of the hyperspectral image data (Figure 1).Figure 1
Location of the DEMMIN test site and its agricultural fields in the northeast of Germany near the city of Demmin. The black polygon outlines the maximum extent of the HyMap data acquired in the years 2004, 2005, and 2008. Sampling points are pictured as black circles. The red flags indicate the location of the two fields shown in Figure13 (eastern flag) and Figure 14 (western flag).


## 2.2. Soil Sampling and Analysis

135 soil samples were collected from the agricultural fields in DEMMIN during a field campaign in September 2006 (Figure1). To obtain a representative set of soil samples of the area, sampling sites were preselected with the help of existing soil maps and a digital terrain model. These were then reconciled with the hyperspectral image data available at this time (2004, 2005) to assure the occurrence of vegetation-free fields with reference samples in these images. Finally, they were adapted to field conditions if appropriate. In total, 14 fields were sampled which stretched across approximately 22 km in north-south direction and 8 km in east-west direction. 8 fields were intensively sampled (6 to 43 points per field). From 6 further fields, samples from 2 to 5 locations were taken. The mean distance between adjacent sampling points within the fields was about 138 m reaching 743 m at maximum. The described spatial distribution of the samples was mainly determined by the pedological setting of the test site which is characterized by extensive areas of sandy and loamy soils, and distinct but smaller regions of loamy and heavy clay soils. In conjunction with the flight windows, soils abound in clay and organic carbon could only be encompassed in the reference data by sampling single fields.Soil samples were taken from the top 5 cm of soil within a radius of 2 m. The geographical coordinates of the sampling points were measured using differential GPS with an accuracy of less than 1 m. At 18 sampling points, five samples were taken within the radius of 2 m and analyzed in the laboratory to assess the error of the applied sampling strategy and to get an idea of the possible impacts of subsampling [44]. We estimated this standard error of sampling (SES) using the following equation from Lozán and Kausch [45]:
(1)SES=∑i=1m(∑j=1kxij2-(∑j=1kxij)2/k)n-m,
where x is the measured value, k is the number of samples taken at each point, m is the number of sampling points, and n is the total number of samples analyzed.Particle size distribution [g kg−1], carbon content [g kg−1], total iron [g kg−1], and soil acidity (pH) were analyzed by the Central Laboratory of the Centre for Agricultural Landscape Research (ZALF) using air-dried samples. Particle size distribution was determined using a combined sieve (fraction 0.063–2 mm) and pipette (fraction <0.063 mm) analysis. The total amount of carbon (TC) was determined by dry combustion using an element analyzer (CNS2000 from LECO). Inorganic carbon content (IC) was analyzed measuring the change in volume during decomposition with phosphoric acid using a Carmhomat 12D. Organic carbon content (OC) was then calculated as the difference between TC and IC. Iron contents were determined by dissolution with aqua regia. Soil acidity was measured by radiometer analytics with a TitraMaster85. The standard error of laboratory (SEL) for soil clay content was 5.0 g kg−1. Analytical techniques for the determination of TC and IC had a SEL of 0.17 g kg−1 and 2.88 g kg−1, respectively.


## 2.3. Spectral Measurements

### 2.3.1. Laboratory Spectroscopy

Spectral measurements were acquired from all soil samples using an ASD spectroradiometer FieldSpec Pro covering the full spectral range from 350 to 2500 nm. The ASD has a fixed spectral sampling interval of 1.4 nm in the visible and near infrared portion of the spectrum (350–1000 nm) and 2 nm in the shortwave infrared (1000–2500 nm) with a full-width half maximum (FWHM) of 3 nm and 10 nm, respectively. All measurements were taken with bare fibre and without foreoptics (25° FOV) relative to a standardized spectralon panel. A 1000 W tungsten lamp was used as light source to illuminate the target at an incident angle of 30°. Before the measurements, all soil samples were grounded with a mortar and passed through a 2 mm sieve to reduce anisotropic light scattering. To standardize the moisture level, they were oven-dried for 48 h at 105°C (378 K).A spectrum of a soil probe was averaged out of ten measurements with a system integration time of 100 scans per single collected spectrum. Noise removal and smoothing was accomplished using the Savitzky Golay filter [46] implemented in IDL [47].


### 2.3.2. Imaging Spectroscopy

Hyperspectral image data from the test site DEMMIN was acquired by the Australian sensor HyMap mounted on a Dornier Do 228 aircraft and a Cessna 208 B during three HyEurope flight campaigns. The HyMap airborne imaging spectrometer records reflected radiances in 128 bands in the wavelength range from 450 to 2480 nm with a spectral bandwidth (FWHM) of 15-16 nm between 450 and 1900 nm and 18–20 nm between 1950 and 2480 nm. The sensor’s field of view (FOV) is 60° with 512 pixel in across track direction. The instantaneous field of view (IFOV) is 2.5 mrad along track and 2.0 mrad across track. Overflights took place on August 9, 2004, May 27, 2005, and July 29, 2008 at an average flight altitude of about 2000 m resulting in a spatial resolution of4×4 m² (Table 1). At each of the campaigns, seven to ten strips of about 16 km length were acquired under clear-sky conditions covering an area of approximately 200 km² and 260 km², respectively.Table 1
Characteristics of available HyMap data sets and its acquisition parameters.
DoA
ToA (UTC)*
Number of strips
Solar angles (°)
Heading** (°) (±2°)
MFA MSL (m)
PS  (m)
August 09, 2004
13:30–14:20
7
Azimuth: 117.2–130Zenith: 46.3–52.2
0/180
~1990
4
May 27, 2005
9:50–10:50
7
Azimuth: 152.7–178.3Zenith: 31.7–34.0
1/180
~1940
4
July 29, 2008
9:50–11:20
10
Azimuth: 147.1–183.1Zenith: 35.3–38.7
0/180
~1934
4
*Local time is given by UTC+2. **Neighboring strips have the opposite flight direction. DoA: day of acquisition, ToA: time of acquisition, MFA ASL: mean flight altitude above mean sea level, and PS: pixel size.Radiometric calibration of the HyMap data was performed inflight using reference ground measurements [48]. Atmospheric correction and conversion of at-sensor radiances in reflectance values was accomplished using the ATCOR-4 software [48, 49]. In a final step, the data was geometrically corrected by means of a generic ortho image processor [50].For each of the sampling points in the field, HyMap spectra were extracted using a window of3×3 pixels centered on the pixel with the coordinates closest to the plot location. The average of these nine spectra was then used for further statistical analysis.


## 2.3.1. Laboratory Spectroscopy

Spectral measurements were acquired from all soil samples using an ASD spectroradiometer FieldSpec Pro covering the full spectral range from 350 to 2500 nm. The ASD has a fixed spectral sampling interval of 1.4 nm in the visible and near infrared portion of the spectrum (350–1000 nm) and 2 nm in the shortwave infrared (1000–2500 nm) with a full-width half maximum (FWHM) of 3 nm and 10 nm, respectively. All measurements were taken with bare fibre and without foreoptics (25° FOV) relative to a standardized spectralon panel. A 1000 W tungsten lamp was used as light source to illuminate the target at an incident angle of 30°. Before the measurements, all soil samples were grounded with a mortar and passed through a 2 mm sieve to reduce anisotropic light scattering. To standardize the moisture level, they were oven-dried for 48 h at 105°C (378 K).A spectrum of a soil probe was averaged out of ten measurements with a system integration time of 100 scans per single collected spectrum. Noise removal and smoothing was accomplished using the Savitzky Golay filter [46] implemented in IDL [47].


## 2.3.2. Imaging Spectroscopy

Hyperspectral image data from the test site DEMMIN was acquired by the Australian sensor HyMap mounted on a Dornier Do 228 aircraft and a Cessna 208 B during three HyEurope flight campaigns. The HyMap airborne imaging spectrometer records reflected radiances in 128 bands in the wavelength range from 450 to 2480 nm with a spectral bandwidth (FWHM) of 15-16 nm between 450 and 1900 nm and 18–20 nm between 1950 and 2480 nm. The sensor’s field of view (FOV) is 60° with 512 pixel in across track direction. The instantaneous field of view (IFOV) is 2.5 mrad along track and 2.0 mrad across track. Overflights took place on August 9, 2004, May 27, 2005, and July 29, 2008 at an average flight altitude of about 2000 m resulting in a spatial resolution of4×4 m² (Table 1). At each of the campaigns, seven to ten strips of about 16 km length were acquired under clear-sky conditions covering an area of approximately 200 km² and 260 km², respectively.Table 1
Characteristics of available HyMap data sets and its acquisition parameters.
DoA
ToA (UTC)*
Number of strips
Solar angles (°)
Heading** (°) (±2°)
MFA MSL (m)
PS  (m)
August 09, 2004
13:30–14:20
7
Azimuth: 117.2–130Zenith: 46.3–52.2
0/180
~1990
4
May 27, 2005
9:50–10:50
7
Azimuth: 152.7–178.3Zenith: 31.7–34.0
1/180
~1940
4
July 29, 2008
9:50–11:20
10
Azimuth: 147.1–183.1Zenith: 35.3–38.7
0/180
~1934
4
*Local time is given by UTC+2. **Neighboring strips have the opposite flight direction. DoA: day of acquisition, ToA: time of acquisition, MFA ASL: mean flight altitude above mean sea level, and PS: pixel size.Radiometric calibration of the HyMap data was performed inflight using reference ground measurements [48]. Atmospheric correction and conversion of at-sensor radiances in reflectance values was accomplished using the ATCOR-4 software [48, 49]. In a final step, the data was geometrically corrected by means of a generic ortho image processor [50].For each of the sampling points in the field, HyMap spectra were extracted using a window of3×3 pixels centered on the pixel with the coordinates closest to the plot location. The average of these nine spectra was then used for further statistical analysis.


## 2.4. Spectral Preprocessing

Laboratory spectra were resampled to fit the spectral resolution of the HyMap sensor. Due to the annual recalibration of the HyMap sensor by HyVista, minor differences in the centre wavelengths and bandwidth (FWHM) occur in between years. To obtain spectral consistency, laboratory spectra and image spectra from the 2004 and 2008 flight campaigns were adjusted to the configuration of the HyMap sensor in 2005. Several spectral transformations were applied. Absorbance was computed by taking the logarithm to the base 10 of the inverse reflectance. First and second derivative of absorbance were computed after the Savitzky Golay method with a fifth-order polynomial including 5 data points either side of each point in the filter. Bands prone to noise were eliminated before statistical analysis. Wavelengths removed were located either side of the water absorption bands near 1400 nm and 1900 nm and below 500 nm and above 2400 nm.


## 2.5. Statistical Analysis

A multitude of statistical approaches comprising principal component analysis, multiple linear regression, artificial neural networks, or regression trees, has been utilized for the quantification of soil parameters from spectroscopic data with varying success (e.g., [10, 19, 51–53]). Partial least squares regression (PLSR) has been exploited by many authors and proofed a good efficiency in quantitative spectroscopy ([11], Table 1). PLSR sets up a linear regression model by projecting the predictor variables onto a few so-called latent variables or factors taking into account both, the response and predictor variation. This way, two problems of conventional multiple linear regression modeling are handled effectively: multicollinearity and variable selection. However, choosing the optimum number of factors is a crucial point in PLS modeling. They were determined applying a leave-one-out cross-validation approach to the calibration data set, choosing the model with the lowest Akaike information criterion (AIC). Rather than just using the RMSE as a selection criterion, AIC also integrates the number of samples and the number of PLS factors avoiding a model overfit more effectively. Statistical analysis of the spectra was carried out using the ParLes Software [54] which implements the orthogonal PLS regression for one y-variable after Martens and Næs [55]. The datasets were randomly split into a calibration set and a validation set. About two third was used for calibration purposes. The rest of the samples were retained from model calibration and used as an independent validation dataset. If the total number of available samples was less than 30, no separate validation was conducted but model accuracy was solely assessed via cross-validation. Calibration and validation data schemes are described in detail in Section 3.5.1.The performance of the models was evaluated by calculating the coefficient of determination (R2), the standard error of calibration (SEC), and the standard error of validation (SEP). The latter is given by the following equation:
(2)SEP=1n-1·∑i=1n(y^i-yi)2,
where y^i and yi are the predicted and observed values of the sample i and n is the total number of samples in the validation data set. SEC is calculated the same way but additionally takes into account the number of PLS factors used in the model. As an independent measure of model accuracy, the ratio of performance to deviation (RPD) was computed by dividing the standard deviation of the reference samples through the SEP [56]. If no independent validation is performed, information on the predictive ability of the model can be retrieved by calculating the ratio of the SEC and the standard deviation (SD) of the calibration data set (SD/SEC). For each model variable, importance in the projection (VIP) and PLS-regression coefficients (b-coefficients) were studied to identify important wavelengths used during model calibration and evaluate their selection with respect to known absorption features of the parameter of interest.To compare parameter predictions between years, root mean squared error was calculated. It is given by the equation(3)RMSE=1n·∑k=1n(y^k,A(ij)-y^k,B(ij))2.
where y^k,A(ij) and y^k,B(ij) are the predicted values of the pixel ij in data set A and in data set B, and n is the total number of pixels.


## 3. Results and Discussion

### 3.1. Soil Samples Descriptive Statistics

Soil texture analysis revealed a prevalence of silty sand, loamy sand, and sandy loam for the soils from the test site DEMMIN with a number of samples representing clayey loam and pure sands. The minimum amount of clay is 36.53 g kg−1 and the maximum is 329.47 g kg−1. Sand and silt contents ranged from 380.17 to 899.62 g kg−1 and from 61.93 to 395.67 g kg−1, respectively. The organic carbon range of all samples was 4.31 to 24.79 g kg−1 but 75% of the samples show OC concentrations less than 11.62 g kg−1. Inorganic carbon content is generally very low but a few samples feature IC contents of up to 10.79 g kg−1 (Figure 2). This may be due to local geology, that is, the occurrence of carbonate rich till, or to liming on agricultural fields. The correlation matrix reveals partly high interrelations between the soil chemical and physical parameters. A high positive correlation exists between clay and iron content (R=0.93). A reasonable relationship was detected for clay and OC. High negative dependencies exist between the sand fraction of the soil and the amount of silt or clay. Negative correlations were also encountered between sand and OC and sand and iron (Table 2).Table 2
Correlation matrix of the physical and chemical soil properties.
Sand
Silt
Clay
IC
OC
Fe2O3
pH
Sand
1
Silt
-
0.81
1
Clay
-
0.85
0.38
1
IC
-
0.03
-
0.15
0.18
1
OC
-
0.71
0.42
0.73
-
0.11
1
Fe2O3
-
0.77
0.32
0.93
0.26
0.57
1
pH
-
0.21
-
0.09
0.42
0.40
0.28
0.42
1Physical and chemical characteristics of the soil samples (n=135) from the test site DEMMIN: particle size fractions sand (63–2000 μm), silt (2–63 μm), and clay (0–2 μm) (a) and total carbon (TC), organic carbon (OC), inorganic carbon (IC), and iron content (Fe2O3) (b).
(a)
(b)Physical and chemical characteristics of the soil samples from the test site DEMMIN exhibit typical properties of agricultural soils on arable land in north eastern Germany [57]. With respect to OC, they furthermore resemble humus conditions on a majority of cropland in Germany [58] but also parts of the central-European agricultural region (e.g., [34]).Findings from the five-fold sampling at 18 points showed that the applied sampling strategy and errors due to subsampling cause a slightly higher uncertainty of the analytical results for clay (7.97 g kg−1) and TC (0.69 g kg−1) than is to be expected by the standard error of laboratory. Measurement errors for IC account for 0.43 gkg−1. This result is in conflict with the SEL for IC which is rather high (cf. Section 2.2). It allows the presumption that the effective standard analytical error for OC for the samples of DEMMIN is not as high as the sum of the standard errors of TC and IC, that is, 3.05 g kg−1 but approximates about 1.12 g kg−1.


### 3.2. Spectral Data Quality

Quality of the spectral data as acquired by laboratory measurements and the HyMap imaging spectrometer was evaluated visually by plotting the reflectance and first derivative of the reflectance of the soil probes and its corresponding sampling points. Figure3 displays the spectral signal of both measurement techniques for a sandy loam with 10.4 g OC kg−1. Maximum reflectance of the sampling point varied between 24% and 34% in case of the HyMap imaging spectrometer and the FieldSpec Pro, respectively. Higher reflectance of the laboratory measurements is probably due to oven-drying of the soil samples as soil moisture reduces the reflectance in the entire wavelength range (e.g., [59]). Soil roughness may also reduce reflectance in the field and is mitigated by sieving the samples before laboratory measurements. Despite the difference in absolute reflectance, the curves are similar in shape. Both show a pronounced absorption feature near 2200 nm. A less distinct and broader absorption feature can be observed near 500 nm and 650 nm which is clearly displayed by the first derivative. The first derivative of the reflectance discloses some differences between the spectra near 880 nm, around 1400 nm and 1900 nm, and at 2350 nm. Deviations around 1400 nm and 1900 nm are probably caused by atmospheric water absorption. It is presumed that the disparity near 880 nm is also related to atmospheric disturbances since a minor water absorption band occurs near 900 nm [14]. Minor differences near 2350 nm are likely a result of different surface conditions in the laboratory and the field where such weak absorption features are less pronounced due to interference with other factors.Reflectance (a) and first derivative of the reflectance (b) of one sampling point as measured under laboratory conditions (black solid line) and by the HyMap imaging spectrometer on May 27, 2005 (red solid line with diamonds indicating centre wavelengths). The soil was classified by soil analysis as a sandy loam with 10.4 g OC kg−1.
(a)
(b)


### 3.3. Image Data Quality

The quality of the atmospheric corrected HyMap images from 2004, 2005, and 2008 was assessed in two ways. First, the general consistency between the HyMap acquisitions of the different years was evaluated. This was done by comparing the mean image spectra of a relatively homogenous concrete area among each other and with field spectral measurements from August 2004. The overall shape of the reflectance curves of the three different years is in good agreement. A slight decline in reflectivity in the shortwave infrared region is observed for the spectra from 2005 and 2008 compared to 2004, whereas reflectivity in the visible and near-infrared region is highest for the spectra from 2008. But apart from a few minor deviations, all spectra range within the standard deviation of the ASD field spectral measurements (Figure4). Second, the relative accuracy of atmospheric correction between image strips of each year was assessed. Spectra of identical image pixels (within the given geometric accuracy) in the overlap of neighboring images were extracted. The ratio of the corresponding spectra of a soil surface and a former runway as a function of wavelength is shown in Figure 5. There is a good match for the runway in all years and for the soil surface in 2005 and 2008. The mean factor of 0.95 is within the expected accuracy of the atmospheric correction including minor BRDF effects during optimal illumination conditions. Increasing noise below 0.5 μm and above 2.4 μm is attributed to the lower signal-to-noise ratio of the HyMap sensor in these spectral regions. Artifacts near 1.4 μm and 1.7 to 1.9 μm are owed to atmospheric water absorption bands. The low ratio of the soil spectra in 2004 indicates significant differences in the overall reflectance of this surface type and occurs systematically in the overlap of all flight strips. Image acquisition parameters in Table 1 suggest that a high solar zenith and a small solar azimuth across the flight heading promote the appearance of BRDF effects on rough soil surfaces in the 2004 data set. This is confirmed by an across-track brightness gradient that occurs after atmospheric correction in all scenes of this year [48]. A comparison to a topographic map with a scale of 1 : 25.000 revealed an absolute geometric accuracy of less than 25 m. The relative geometric accuracy of the image data between the years amounts to two pixels in x or y direction, that is, 8 m.Figure 4
Comparison of mean reflectance spectra of a concrete area on a farm premise near the city of Kruckow acquired with an ASD Fieldspec Pro in 2004 and the HyMap sensor in 2004, 2005, and 2008.Ratio of spectra of identical image pixels (within the given geometric accuracy) in the overlap of neighboring scenes for a former runway and a soil surface.
(a)
(b)


### 3.4. Preparing Image Data for Spatially Explicit Parameter Estimation

In the central European agricultural region, bare soil is only encountered in early spring while fields are prepared for sowing of sugar beet, maize, potatoes, or summer grain, and in late summer after harvesting and seedbed preparation for winter grains (e.g., winter barley, winter wheat). Mitigating the impact of vegetation or excluding vegetated fields from the analysis is therefore a prerequisite for a successful application of spectroscopic approaches. In a preliminary study, Daughtry et al. [60] showed that photosynthetic active vegetation, nonphotosynthetic vegetation, and bare soils could be distinguished by a combination of vegetation indices, the cellulose absorption index (CAI) [61, 62], and the normalized difference vegetation index (NDVI). They obtained better results on AVIRIS data with the combination of indices than with spectral unmixing. In the present study, this approach was adopted to identify fields with bare soils and mask pixels with vegetation signal before image data analysis. Indices were computed for all HyMap images using wavelengths near 680 and 800 nm for the NDVI and wavelengths near 2000, 2100, and 2200 nm for CAI. Thresholds for green and dry vegetation were determined interactively and knowledge based and were set to 0.2 for NDVI and −1 to −2.5 for CAI. The quality of the soil mask was assessed using 200 ground truth regions of interest randomly distributed in equal parts over the agricultural area classified as vegetation and bare soil. The classification of the pixels was evaluated knowledge based with the help of the corresponding image spectra. There was almost no misinterpretation for pixels classified as vegetation. Disagreements in the soil class were mostly related to emerging green vegetation as some of the ground truth regions showed a small but discrete absorption feature near 680 nm but NDVI values below the assumed threshold. Accuracy statistics revealed a producer’s accuracy of 99% and 84% and a user’s accuracy of 86.1% and 98.8% for vegetation and soil, respectively. The kappa coefficient was 0.83 indicating a good overall agreement for the two observed classes.A maximum of 8.21% of the total area covered by the HyMap scenes from May 2005 was detected free of vegetation applying the combined threshold of CAI and NDVI. In the data sets from August 2004 and July 2008, only 2.2% and 1.29% of the total image area were identified as bare soils and made available for topsoil parameter estimation (Table3). Since the total image area covered by the HyMap data in 2008 exceeds the total area covered in 2004 by roughly one-third, the absolute number of image pixels available for image data analysis is within the same dimension.Table 3
Area free of vegetation in the HyMap images from 2004, 2005, and 2008 given in absolute values and percentage of the total image coverage.
Time of acquisition
Total number of image pixels
Image pixels with NDVI<0.2 and-2.5<CAI<-1.0
Proportion of image pixels with NDVI<0.2 and -2.5<CAI<-1.0(%)
09.08.2004
12.076.995
266.019
2.20
27.05.2005
12.495.804
1.026.187
8.21
29.07.2008
16.125.125
208.374
1.29In Figure6, CAI is plotted against NDVI for all 135 sampling points derived from HyMap images of the three years. In 2005, sixty-seven sampling points reside on fields without vegetation. The remaining sixty-eight points are situated on fields with green vegetation indicated by NDVI values above 0.6. In 2004 and 2008, respectively, the spectral signal of only 20 and 11 sampling points suggested the appearance of bare soil. Apart from a few data points with higher NDVI values, the majority of sampling points exhibit NDVI less than 0.2 but CAI values above zero which is clear evidence for nonphotosynthetic, dry vegetation. Due to the time of image data acquisition, these fields are either still grown with senescent crops or covered with stubbles which remain one the fields after harvesting.Figure 6
CAI versus NDVI derived from HyMap spectra of the years 2004, 2005, and 2008 of the 135 sampling points in the field.


### 3.5. Model Calibration and Validation

#### 3.5.1. Calibration and Validation Schemes

Calibration and validation schemes were set up with respect to the availability of vegetation-free reference points in the HyMap data from 2004, 2005, and 2008. In total, three calibration data sets were generated. The first one comprises the spectra of the 20 sampling points identified in the image data from 2004 (Lab_2004, HyMap_2004). The second calibration set contains 38 out of 67 vegetation-free spectra in the image data from 2005 (Lab_2005, HyMap_2005). The third-one is made up of 11 spectra corresponding with the vegetation-free sampling points in the image data acquired during the flight campaign in 2008 (Lab_2008, HyMap_2008). With the laboratory spectroscopy, each of the calibration models were validated with the rest of the laboratory spectra withheld during calibration that are 115, 97 and 124 in number. For the calibration data set with 38 soil spectra, an additional validation was performed utilizing only those 29 samples which are located on fields with bare soil in the image data from 2005 (Lab_2005b). Validation of the HyMap models was confined to cross-validation in case of the HyMap_2004 and HyMap_2008 models due to the lack of vegetation-free reference points. The HyMap_2005 model was additionally validated with the 29 samples not used during calibration. Descriptive statistics of all calibration and validation data sets are summarized in Table4. As sampling has been carried out field-by-field with a mean within-field sampling distance of adjacent sampling points of 138 m, validation cannot be considered to be fully independent from the calibration data because of spatial autocorrelation. Regardless of how the validation subsample is selected, this can lead to an overestimation of the predictive accuracy of the models in particular for unsampled fields [63].Table 4
Descriptive statistics of the calibration and validation data sets for the prediction of clay and organic carbon content with laboratory and HyMap image spectra.
Calibration data sets
Laboratory
HyMap
Number of samples
Clay (g kg−1)
Organic carbon  (g kg−1)
Mean
SD
Min
Max
Mean
SD
Min
Max
Lab_2004
HyMap_2004
20
108.07
48.38
48.03
186.64
9.60
3.17
5.25
19.59
Lab_2005
HyMap_2005
38
99.12
50.05
38.45
216.18
9.16
2.95
5.01
17.48
Lab_2008
HyMap_2008
11
97.31
45.68
52.22
197.20
6.66
1.47
4.31
10.34
Validation data sets
Laboratory
HyMap
Number of samples
Clay (g kg−1)
Organic carbon (g kg−1)
Mean
SD
Min
Max
Mean
SD
Min
Max
Lab_2004
—
115
130.22
67.45
36.53
329.47
10.42
4.02
4.31
24.79
Lab_2005a
—
97
137.83
67.53
36.53
329.47
10.74
4.16
4.31
24.79
Lab_2005b
HyMap_2005
29
92.06
52.36
36.53
237.52
8.89
3.73
5.55
24.79
Lab_2004
—
124
129.57
66.26
36.53
329.47
10.62
3.90
5.01
24.79Except for the cal/val configuration with 38 and 29 samples, data range and standard deviation of the calibration data sets were considerably smaller than in the validation data sets. For the most part of spectroscopic studies, validation is performed either on data resembling the data range of the calibration data set or solely by cross-validation as model validity is restricted to the data range covered during calibration. Here, the described procedure was chosen to evaluate model performance outside the calibration data range since soil characteristics of the test site DEMMIN are not adequately represented by the small number of vegetation-free sampling points in the image data sets.


#### 3.5.2. Laboratory Spectra

Table5 summarizes calibration and validation statistics of the PLSR models for the prediction of clay and organic carbon content. Tests with varying preprocessing methods showed that either absorbance or first derivative of the absorbance gave the best prediction results. Calibration data sets with a greater number of samples generally yield higher R2 and lower standard errors of calibration (SEC) than models with a fewer number of samples. For clay, R2 ranged from 0.70 to 0.90. SEC was lowest in the Lab_2005_Clay model (17.42 g kg−1) and highest in the Lab_2008_Clay model (35.68 g kg−1). Standard error of validation (SEP) varied from 19.48 to 35.55 g kg−1. Plotting measured against predicted clay contents reveals a close relationship and a homogeneous scatter along the 1 : 1 line for the Lab_2005_Clay model validated with 29 and 97 soil samples (Figure 7). For the Lab_2004_Clay model, measured and predicted values are generally in good agreement but soils with more than 200 g kg−1 clay, of which none was enclosed in the calibration data set, are slightly overestimated. Samples show a systematic, linear offset from the 1 : 1 line. However, deviations are within the error variance of the whole validation data set. The calibration model with the smallest number of samples (Lab_2008_Clay) caused large deviations for some validation samples with less than 200 g kg−1 clay but predictions of clay rich soils are less erroneous (Figure 7).Table 5
Calibration and validation performance statistics of laboratory spectra.
Parameter
Clay
Model name
Number of samples(cal/val)
Preprocessing
Number of LV
SDcal†
SDval†
R
²
SEC
SD/SEC
SEP
RPD
Lab_2004_Clay
20/115
Absorption
6
48.21
75.21
0.76
29.30
1.65
26.59
2.54
Lab_2005_Clay
38/97
First derivative of absorption
5
49.32
70.06
0.90
17.42
2.87
20.77
3.25
Lab_2005_Clay
38/29
First derivative of absorption
5
49.32
48.08
0.90
17.42
2.87
19.48
2.69
Lab_2008_Clay
11/124
Absorption
5
39.92
69.18
0.70
35.68
1.28
35.55
1.86
Parameter
Organic carbon
Model name
Number of samples(cal/val)
Preprocessing
Number of LV
SDcal†
SDval†
R
2
SEC
SD/SEC
SEP
RPD
Lab_2004_OC
20/115
Absorption
2
2.53
3.77
0.51
2.36
1.34
2.36
1.71
Lab_2005_OC
38/97
Absorption
6
2.87
3.90
0.88
1.11
2.67
1.92
2.16
Lab_2005_OC
38/29
Absorption
6
2.87
2.63
0.88
1.11
2.67
1.94
1.93
Lab_2008_OC
11/124
First derivative of absorption
4
1.49
4.12
0.78
0.97
1.52
2.46
1.58
†SD: standard deviation of estimated parameter values in (g kg−1), R²: coefficient of determination, SEC: standard error of calibration in (g kg−1), SEP: standard error of validation in (g kg−1), RPD: ratio of performance to deviation.Plots of measured versus predicted values for clay (gkg−1) in the calibration (black symbol, ×) and validation (red symbol, +) data sets based on laboratory spectra. The respective calibration and validation data schemes (ncal./nval.) are marked on the plots. For further explanation see text.
(a)
(b)
(c)
(d)Calibrations based on at least 20 samples had RPD values clearly above 2, indicating excellent prediction models. They can be assigned to group A in the classification of Chang et al. [56]. Reducing the number of samples in the calibration data set to just 11 spectra, results in less predictive power with an RPD of 1.86. Such models belong to group B and may be improved by different calibration strategies [56]. Similar results on the prediction of clay from laboratory spectra are reported, for example, from Waiser et al. [64], Chang et al. [56], and Cozzolino and Morón [65]. The latter performed a cross-validation with a modified PLSR approach on spectra from Uruguayan soils and achieved a value of 2.7 (ratio of standard deviation by standard error of cross validation). Waiser et al. [64] tested four different soil preparation techniques including air-dried ground samples from Texas and achieved RPD values between 1.95 and 3.51 with PLSR. Principal component regression and cross validation were applied by Chang et al. [56]. They obtained an RPD of 1.7. In contrast to our work, the range of clay contents in the calibration and validation data sets in those studies were similar. The small range of clay concentrations in the soil samples of the calibration data sets varying from about 40 to 220 g kg−1 clay did hardly affect the validity of our models even for samples with higher clay contents. But a very small number of samples hamper the generation of a solid calibration model for the prediction of clay content using VNIR spectroscopy.R
2 values of the PLSR models on organic carbon ranged from 0.51 to 0.88 with SEC between 1.11 and 2.36 g kg−1. SEP varied between 1.92 and 2.46 g kg−1. The best results were again achieved providing 38 samples in the calibration data set (Lab_2005_OC) having the highest R2 and lowest standard errors. Calibration models for organic carbon exhibit a comparable predictive ability as those for clay. The maximum ratio of standard deviation to SEC is 2.67 for the Lab_2005_OC model and minimum 1.34 for the Lab_2004_OC model. In contrast to the prediction of clay, the distribution of data points in the scatter plots for the validation data sets indicated less precise models (Figure 8). A good agreement is observed between measured and predicted values in the validation data sets for the Lab_2005_OC model for samples with less than 15.0 g kg−1 OC. Samples with higher concentrations show a trend of underestimation. Validating the models with 20 and 11 calibration samples reveals a great variety in the measured versus predicted plots which seems to increase with higher organic carbon contents. This decrease in prediction accuracy is reflected by lower RPD values. With the exception of the Lab_2005_OC model which has a RPD close to 2 or above 2, all models must be categorized in group B following the classification of Chang et al. [56]. The precision obtained for OC in our study is within the precision obtained by previous authors (e.g., [34, 51, 66]). Stevens et al. [34], for example, generated a PLSR model with 117 soil samples from agricultural fields in Belgian Lorraine with a RPD of 2.11. Volkan Bilgili et al. [66] predicted several soil properties of a 32 ha field in the north of Turkey among them soil organic matter (SOM). In their study, a higher data range in the validation data set (3.9–68.7 g kg−1 SOM) than in the calibration data set (5.0–33.8 g kg−1 SOM) did not result in a decrease of the models predictive ability, neither for PLSR (RPD: 1.93) nor with multivariate adaptive regression splines (MARS) (RPD: 1.90). This may be due to the great number of samples in their cal/val sets (153/359).Plots of measured versus predicted values for organic carbon (g kg−1) in the calibration (black symbol, ×) and validation (red symbol, +) data sets based on laboratory spectra. The respective calibration and validation data schemes (ncal./nval.) are marked on the plots. For further explanation see text.
(a)
(b)
(c)
(d)Resampling laboratory spectra to HyMap, spectral resolution had only little affect on the prediction accuracy of clay and none on the prediction accuracy of organic carbon (results not shown). Similarly, Gomez et al. [21] reported little influence on prediction accuracy while resampling spectra to Hyperion resolution.Variable importance in the projection (VIP) and PLS-regression coefficients (b-coefficients) provide information about the weight of the x variables (wavelengths) in the PLSR model. Large b-coefficients indicate wavelengths which are relevant in the modeling of Y (soil parameter). High VIP values point out wavelengths which are essential for the modeling of both, Y and X [67]. Thus, x variables with a high VIP and a high b-coefficient indicate important wavelengths and allow insights on the physical basis of the predictions. VIP values and b-coefficients of all laboratory models are plotted in Figure 9. Curve shapes of the Lab_2005_Clay model and the Lab_2008_OC model have a different appearance compared to the rest due to the applied spectral preprocessing (first derivative of absorbance). For clay, VIP and b-coefficients indicate an emphasis of wavelengths between 2000 and 2500 nm and in particular around 2200 nm across all models. In fact, clay minerals are spectrally active in this part of the spectrum [12]. Absorption features near 2200 to 2300 nm produced by combination vibrations involving an OH stretch and metal-OH bend are diagnostic for all clay minerals. The exact position of the absorption features varies depending on the metal which is bond to the hydroxyl group [13, 14]. Furthermore, VIP and b-coefficients suggest that wavelengths in the visible part of the spectrum near 500 nm and 650 nm contain important information on clay minerals or rather clay related features as the spectral activity of clays is confined to the infrared portion of the spectrum. Iron, either as main constituent or associated with the structure of clay minerals, shows spectral features in the VNIR due to electronic transitions of iron cations [14]. Hunt and Salisbury [12], for example, attributed weak absorption bands near 500 nm to ferrous iron in montmorillonites. In illites, absorption near 650 nm was observed by Clark et al. [13] who ascribed this feature to ferric iron. The very close correlation of iron and clay (Table 2) in the soil samples of the test site DEMMIN encourages the assumption that iron has an influence on the PLS regression for clay.Variable importance in the projection (VIP) andb-coefficients of PLSR calibration models based on laboratory spectra using 20, 38, and 11 soil samples for clay (left side) and organic carbon (right side).
(a)
(b)
(c)
(d)
(e)
(f)VIP andb-coefficient curves of the organic carbon models exhibit less consistency than for clay. Whereas high loadings can be observed for the Lab_2005_OC and Lab_2008_OC model between 2000 and 2500 nm similar to clay, this region seems of little importance in the Lab_2004_OC model although a small accentuation occurs near 2200 nm. In the latter, generally high VIP values and b-coefficients prevail in the 500 to 900 nm region. In the Lab_2005_OC and Lab_2008_OC models wavelengths near 650 nm and 500 to 600 nm, respectively, are emphasized. As organic carbon in soils consists of very complex chemical compounds having numerous, frequently overlapping absorption features [68, 69], a spectral feature assignment to specific wavelengths is difficult [7, 70]. In addition, good correlations between clay and organic carbon or iron and organic carbon hamper the elucidation of the causal relationships between model parameters and its physical background.


#### 3.5.3. HyMap Image Spectra

In a first attempt, it was tried to set up a stable calibration model using the image data from 2005 and validating this model with the remaining sampling points from 2004, 2005, and 2008. Model calibration and validation with image spectra from 2005 yielded good results with RPD of 2.70 and 1.80 for clay and organic carbon, respectively (Table7). Validation with image spectra from different flight campaigns, that is, 2004 and 2008, reached RPD below 1.0 indicating poor results with no or little predictive ability especially for clay (Table 6). This is the case although soil clay and OC contents of the respective sampling points are within the data range of the calibration model and spectral data quality assessment showed a good overall agreement between years (cf. Section 3.2). Experiments with laboratory spectra revealed that the composition of the validation data set did only slightly impair the model’s predictive ability. The problem with empirical calibration is that coefficients are precisely adjusted to the spectral properties of the calibration input data. Even small changes in spectral data properties hamper a reasonable application. Laboratory spectra are measured under standard conditions where variations in spectral properties are solely a consequence of variations in physical and chemical characteristics. Under field conditions, soil status (e.g., roughness, crusting, moisture) varies with time. Further variations may be introduced through radiometric and atmospheric correction (see Figure 4). Thus, it was not possible to set up one stable calibration model for the estimation of soil parameters from image data but data set specific calibrations seem to be crucial to retrieve reasonable prediction accuracies.Table 6
Validation performance statistics of HyMap_2005 model transfers to HyMap spectra from different acquisition years.
Transfer of HyMap_2005 model to
Clay
Organic carbon
SEP (g kg−1)
RPD
SEP (g kg−1)
RPD
HyMap spectra 2004
175.04
0.28
4.82
0.66
HyMap spectra 2008
239.66
0.19
1.58
0.93Table 7
Calibration and validation performance statistics of HyMap image spectra.
Parameter
Clay
Model name
Number of samples(cal/val)
Preprocessing
Number of LV
SDcal†
SDval†
R
²
SEC
SD/SEC
SEP
RPD
HyMap_2004_Clay
20
Absorption
3
44.72
—
0.80
23.39
2.07
—
—
HyMap_2005_Clay
38/29
First derivative of absorption
3
47.18
50.49
0.85
20.34
2.46
19.41
2.70
HyMap_2008_Clay
11
Absorption
3
42.36
—
0.92
15.99
2.86
—
—
Parameter
Organic carbon
Model name
Number of samples(cal/val)
Preprocessing
Number of LV
SDcal†
SDval†
R
²
SEC
SD/SEC
SEP
RPD
HyMap_2004_OC
20
First derivative of absorption
3
2.70
—
0.62
2.13
1.48
—
—
HyMap_2005_OC
38/29
First derivative of absorption
2
2.58
2.20
0.71
1.64
1.80
2.07
1.80
HyMap_2008_OC
11
First derivative of absorption
5
1.20
—
0.45
1.61
0.91
—
—
†SD: standard deviation of estimated parameter values in (g kg−1), R2: coefficient of determination, SEC: standard error of calibration in (g kg−1), SEP: standard error of validation in (g kg−1), RPD: ratio of performance to deviation.Results of data set specific PLS regressions with the HyMap spectra from 2004, 2005, and 2008 are listed in Table7. Calibration and validation were performed as described in Section 3.5.1 (Table 4). As with laboratory spectra prediction accuracy was higher for clay than for organic carbon but a high number of samples in the calibration data set did not necessarily yield the best results. For clay, R2 ranged from 0.80 to 0.92 with SEC as low as 15.99 to 23.39 g kg−1. SD divided by SEC is highest for the HyMap spectra from 2008 reaching a value of 2.86 (HyMap_2008_clay). The calibration model for the image spectra from 2005 is characterized by a SD/SEC ratio of 2.46 (HyMap_2005_clay). Independent validation with the remaining 29 data points resulted in a RPD of 2.70 confirming a very high prediction accuracy of this model similar to the image spectra from 2008. The lowest prediction ability is achieved with the calibration model based on image spectra from 2004 with a SD/SEC ratio of 2.07 (HyMap_2004_clay). Measured versus predicted values suggest that samples with low clay contents (<80 g kg−1) may not be adequately predicted applying this model as variability of estimations is very high for samples with more or less identical clay contents (Figure 10). According to Chang et al. [56], all three models can be classified as group A with an excellent ability to predict soil clay content from hyperspectral image data. A direct comparison of our results to the prediction accuracy obtained by other authors is difficult since none of them reports RPD or SD/SEC values. Standard errors may only be compared between studies if soils have an analogous range of clay content. This is the case with the work of Selige et al. [19] who reported a cross-validated R2 of 0.65 with a RMSE of 38 g kg−1 for MLR and a R2 of 0.71 with a RMSE of 42 g kg−1 for PLSR for estimations of clay from HyMap data in an agricultural area in the German federal state Saxony-Anhalt. Although their samples reached a maximum of 260 g kg−1 clay, these error measures clearly indicate that prediction accuracy for clay from HyMap image data in DEMMIN is superior.Plots of measured versus predicted values for clay (gkg−1) in the calibration (black symbol, ×) and validation (red symbol, +) data sets based on HyMap image spectra. The respective calibration and validation data schemes (ncal./nval.) are marked on the plots. For further explanation see text.
(a)
(b)
(c)R
2 obtained from model calibration for the prediction of OC content from HyMap spectra are 0.71 at most and 0.45 at least with SEC between 1.61 and 2.13 g kg−1. Prediction ability of the PLSR model based on the 2005 image data is acceptable with a SD/SEC ratio and an RPD of 1.80. Being restricted to 29 validation samples, no clear conclusion can be drawn on the model’s behaviour above 20 g kg−1 OC. Similar to the corresponding laboratory cal/val configuration (Lab_2005_Clay, nval=29), a strong underestimation of one data point with about 25 g kg−1 OC is observed (Figures 8 and 11). This gives reason to interpret that the HyMap_2005_OC model is subject to a nonlinear trend just as examined in the validation of the laboratory model with 97 samples. Calibration models based on image data from 2004 and 2008 exhibit SD/SEC values of 1.49 and 0.91, respectively. Thus, models for the prediction of OC based on image spectra belong to group B (HyMap_2004_OC, HyMap_2005_OC) and group C (HyMap_2008_OC) after Chang et al. [56]. Prediction accuracy obtained by PLS regression on HyMap images from 2004 and 2005 match results from Stevens et al. [34] and Stevens et al. [23] applying image data of the sensors AHS-160 (RPD: 1.40) and CASI (RPD: 1.86). Stevens et al. [20] showed that regrouping a global cal/val data set by soil type into several smaller cal/val data sets could improve prediction accuracy about two times with an RPD as high as 3.15. However, this approach is only feasible if the number of reference points is large and if auxiliary information on soil type is provided for delineation of soil types in the image data in order to correctly apply the different calibration models.Plots of measured versus predicted values for organic carbon (gkg−1) in the calibration (black symbol, ×) and validation (red symbol, +) data sets based on HyMap image spectra. The respective calibration and validation data schemes (ncal./nval.) are marked on the plots. For further explanation see text.
(a)
(b)
(c)Since laboratory measurements are acquired under standardized conditions, it was expected that calibrations with HyMap image spectra results in a deterioration of prediction accuracy. Lagacherie et al. [22], for example, attributed a decrease in model performance from laboratory to imaging spectroscopy to radiometric and wavelength calibration uncertainties of the HyMap sensor as well as atmospheric effects. Gomez et al. [21] assumed that a higher signal-to-noise ratio and the large spatial support area (30×30 m²) of Hyperion images caused a decline in prediction ability. Other factors degrading the model performance may be variable soil moisture content [23] or soil surface conditions [71] as vegetation affected spectra were removed before the image data analysis. In fact, higher standard errors were observed in this study for the HyMap_2005_OC and the HyMap_2008_OC models. Apart from the factors mentioned above, the interannual variability of soil organic carbon [40, 41] may have influenced the spectral calibrations for these two years. Since image data was analyzed with reference data from 2006, possible variations of soil organic carbon caused by the specific temperature and soil moisture regime of each year could not be taken into account. Prediction ability of the HyMap_2005_Clay model remained more or less constant (Tables 5 and 7). Judging from the SD/SEC values model performance increased when applying the image data from 2004 (HyMap_2004_OC, HyMap_2004_Clay). Accuracy also improved when estimating clay using the data from 2008 although it decreased when estimating organic carbon with the same data set. A reasonable explanation for this behavior was not found. In any case, it must be stated that only a reverification of the predictive accuracy with totally independent reference samples can reveal further insight and the true predictive ability of the presented PLSR models. This is of particular importance for fields not sampled during field work. To obtain a first idea of possible implications caused by spatial autocorrelation, a calibration and validation scheme for the HyMAP spectra from 2005 was tested in which all samples from three fields were retained from calibration at a time and used for model validation only. This was done for five different cal/val constellations. Standard errors of validation (SEP) retrieved this way for clay were less than the SEP of the HyMAP_2005_Clay model in three cases. In two of these cal/val schemes, SEP increased by about 11% and 15%, that is, a maximum SEP of 22.3 g kg−1 was obtained. For OC, SEP increased by about 9% and 29% for two cal/val schemes. The top most SEP was 2.68 g kg−1. In the remaining test cases considered, SEP was less than the SEP of the HyMAP_2005_OC model. These results corroborate the urgent need for independent validation. Although an increase in standard errors by about 30% may still seem acceptable, parameter estimations for soils with low clay or OC contents would become very inaccurate.VIP andb-coefficient curves of the PLSR models for clay based on HyMap spectra presented in Table 7 signalize an importance of wavelength between 2000 and 2500 nm just as observed for the laboratory models. Again special emphasis is given to the 2200 nm band. This is less pronounced in the HyMap_2008_Clay model where wavelengths in the visible and near-infrared part of the spectrum seem to play a major role. Other than with the laboratory spectra, wavelengths between 900 and 1300 nm contain high values in this model with a slight accentuation near 1100 nm. This inflection can also be seen in the Lab_2008_Clay model and may be associated with absorption caused by ferrous iron in illite minerals [13]. Increased VIP and b-coefficients in all three models can also be found in the 500 to 650 nm region and are probably linked to iron as discussed with the laboratory spectra. For organic carbon, several wavelengths in the visible part of the spectrum are stressed. Peaks occur at 550, 680, and 770 nm. Likewise, bands near 2000, 2200, and 2350 nm carry high loadings in all three models similar to the laboratory models. As has already been stated in Section 3.5.2, there are noticeable analogies to the clay model parameters.


### 3.6. Spatially Explicit Parameter Estimation

In Figure12 results of the multi-annual soil parameter estimation using hyperspectral image data and PLS regression is exemplarily displayed for the parameter clay. Flight strips of each flight campaign were processed individually. Mosaics were created subsequently. Since statistical analysis was done for soils from arable land only, quantitative mapping was restricted to agricultural crop land by means of a vector layer with field boundaries. White regions represent nonarable land (e.g., settlements, forest, and grassland). Grey regions indicate arable land completely covered with vegetation which was masked using CAI and NDVI (cf. Section 3.4). Supplementary information is given on the field specific model predictive ability classified after Chang et al. [56] including details on the respective calibration data range (Figure 12(b)). Careful documentation of error measures and calibration data range is considered essential as it differs with multi-annual image data. Figure 12(c) provides direct access to the field-specific image data source. As expected, clay content is less than 125 g kg−1 within most of the fields. Despite this predominance of sandy loam and loamy sand, the map reveals a typical characteristic of the DEMMIN area which is the occurrence of partly large local heterogeneities of the soil type. A single field of 40 ha near the village Heydenhof is an example (Figure 13) with clay contents close to 300 g kg−1. Similarly, organic carbon contents are above average in the eastern part of the field reaching up to 17.5 g kg−1. The spatial pattern of clay and organic carbon content in this field shows little congruence. This is in agreement with a moderate correlation coefficient (Table 2) between the two soil properties and considerable deviations that occur for a few sandy soils with very high organic carbon contents but very little clay. Spatial pattern of clay and OC in this field is most likely not related to topography as the area is very flat.Figure 12
Map of soil clay content on agricultural fields (a) derived by multi-annual hyperspectral image data analysis. Map of RPD values (b) indicate field-specific model predictive ability after Chang et al. [56]. A field assignment to the year of acquisition is illustrated in Map (c).Figure 13
Map of soil clay and organic carbon content derived from HyMap image data (May 27, 2005) and PLSR for a field near the village of Heydenhof. The location of this field on the test site DEMMIN is marked by a red flag in Figure1.Performing multi-annual estimations of soil parameters offers the opportunity to compare parameter estimations from different years which provides an additional measure of the model’s predictive ability. Because of the time of image data acquisition, only two fields could be identified where parameters were estimated in either 2004 and 2005 or 2005 and 2008. No overlap was detected for all three years. Figure14 displays the difference image of estimated clay and organic carbon contents based on HyMap images from 2004 and 2005. During the overflight in 2004, this field was partly covered with crop residues which are more or less arranged in line according to the tillage direction. Furthermore, the state of tillage differed within the field. In the eastern part, vegetation residues almost completely cover the soil surface. These pixels are masked (white). On average, predicted clay contents in 2005 are 3.42 g kg−1 below the estimates in 2004. 50% of all data points show differences between −14.19 and 9.11 g kg−1 clay. Mean differences of predicted OC contents are close to 0 whereby half the number of data points differ between −0.11 and 0.57 g kg−1 OC. The spatial pattern of both difference images is similar. Higher negative differences appear with both parameters in the upper left of the field (dark blue). A small area in the centre of the field exhibits higher positive deviations (dark orange-red). Differences reach extreme values for a few isolated pixels. They are attributed to artifacts or problems caused by mixed pixels. RMSE values calculated for both parameters are equal to or less than standard errors of model calibration (Table 7). This result is encouraging as it confirms statements made on model accuracies and emphasizes the potential of imaging spectroscopy to map soil parameters. However, clay and OC contents on this particular field are quite low. As laboratory models showed an increase in variance with increasing OC concentrations, the actual RMSE calculated for OC may be higher. This was not the case for clay but only further investigations may provide additional evidence.Figure 14
Difference image of estimated clay and organic carbon content based on HyMap images from 2004 and 2005 for an agricultural field near the city of Marienfelde. Q1:  lower quartile, Q3: upper quartile. Mean deviation (RMSE) calculated from estimations in 2005 and 2004. The location of this field on the test site DEMMIN is marked by a red flag in Figure1.Reflectance of most natural surfaces varies depending on the viewing and solar illumination conditions. This fact is described by the bidirectional reflectance distribution function (BDRF) [72]. Evaluating prediction accuracy in the overlap area of the flight strips per year may give an indication for inaccuracies related to these reflectance anisotropies. Stevens et al. [20] presumed that difficulties to produce robust calibrations were partly caused by strong backward- or forward-scattering at the image edges due to the large field of view of the AHS sensor (90°). Computing RMSE revealed substantially higher errors in the overlap area for the 2004 data set than was expected from the validation, respectively, cross-validation of the models (Table 8). For clay, RMSE was about 3 times larger than SEC using reflectance transformed to absorbance. Then applying reflectance transformed to the first derivative of absorbance RMSE was about 1.4 times larger than SEC. Effects were less pronounced for organic carbon but RMSE still exceeded SEC by a factor of 1.6 using absorbance. In contrast, RMSE were either only slightly increased or even less than standard errors of validation for HyMap spectra from 2005 and 2008 independent from the preprocessing method applied. These results provide strong evidence that BRDF effects have negative influence on the performance of the models based on image data from 2004. As no negative impact was observed in the two other data sets, this effect seems solely due to the unfavorable viewing and illumination conditions during the overflight in 2004, that is, a large zenith and a low azimuth angle (Table 1), but not to the sensors field of view in principle. Negative effects in the off-nadir area can be diminished applying the first derivate even though a somewhat higher standard error must be accepted.Table 8
Mean deviation (RMSE) of predicted clay and organic carbon contents in the overlap area of the HyMap flight strips and standard errors of PLSR models generated using the spectra transformed to absorbance (Abs) and spectra transformed to first derivative of absorbance (Abs_Deriv1). Bold numbers indicate model parameters presented in Table7.
Day of acquisition
Number of pixels estimated
Clay
Organic carbon
RMSEoverlap
SEP* or SEC**
RMSEoverlap
SEP* or SEC**
Abs
Abs_Deriv1
Abs
Abs_Deriv1
Abs
Abs_Deriv1
Abs
Abs_Deriv1
August 9, 2004
36015
69.63
38.99
23.3
9
*
*
28.35**
3.86
2.20
2.35**
2.1
3
*
*
May 27, 2005
260206
22.46
17.03
20.96*
19.4
1
*
2.34
0.60
3.13*
2.0
7
*
July 29, 2008
60308
12.95
15.39
15.99
*
*
22.89**
2.21
0.79
1.88**
1.6
1
*
*
SEC: standard error of calibration in (g kg−1) and SEP: standard error of validation in (g kg−1).


## 3.1. Soil Samples Descriptive Statistics

Soil texture analysis revealed a prevalence of silty sand, loamy sand, and sandy loam for the soils from the test site DEMMIN with a number of samples representing clayey loam and pure sands. The minimum amount of clay is 36.53 g kg−1 and the maximum is 329.47 g kg−1. Sand and silt contents ranged from 380.17 to 899.62 g kg−1 and from 61.93 to 395.67 g kg−1, respectively. The organic carbon range of all samples was 4.31 to 24.79 g kg−1 but 75% of the samples show OC concentrations less than 11.62 g kg−1. Inorganic carbon content is generally very low but a few samples feature IC contents of up to 10.79 g kg−1 (Figure 2). This may be due to local geology, that is, the occurrence of carbonate rich till, or to liming on agricultural fields. The correlation matrix reveals partly high interrelations between the soil chemical and physical parameters. A high positive correlation exists between clay and iron content (R=0.93). A reasonable relationship was detected for clay and OC. High negative dependencies exist between the sand fraction of the soil and the amount of silt or clay. Negative correlations were also encountered between sand and OC and sand and iron (Table 2).Table 2
Correlation matrix of the physical and chemical soil properties.
Sand
Silt
Clay
IC
OC
Fe2O3
pH
Sand
1
Silt
-
0.81
1
Clay
-
0.85
0.38
1
IC
-
0.03
-
0.15
0.18
1
OC
-
0.71
0.42
0.73
-
0.11
1
Fe2O3
-
0.77
0.32
0.93
0.26
0.57
1
pH
-
0.21
-
0.09
0.42
0.40
0.28
0.42
1Physical and chemical characteristics of the soil samples (n=135) from the test site DEMMIN: particle size fractions sand (63–2000 μm), silt (2–63 μm), and clay (0–2 μm) (a) and total carbon (TC), organic carbon (OC), inorganic carbon (IC), and iron content (Fe2O3) (b).
(a)
(b)Physical and chemical characteristics of the soil samples from the test site DEMMIN exhibit typical properties of agricultural soils on arable land in north eastern Germany [57]. With respect to OC, they furthermore resemble humus conditions on a majority of cropland in Germany [58] but also parts of the central-European agricultural region (e.g., [34]).Findings from the five-fold sampling at 18 points showed that the applied sampling strategy and errors due to subsampling cause a slightly higher uncertainty of the analytical results for clay (7.97 g kg−1) and TC (0.69 g kg−1) than is to be expected by the standard error of laboratory. Measurement errors for IC account for 0.43 gkg−1. This result is in conflict with the SEL for IC which is rather high (cf. Section 2.2). It allows the presumption that the effective standard analytical error for OC for the samples of DEMMIN is not as high as the sum of the standard errors of TC and IC, that is, 3.05 g kg−1 but approximates about 1.12 g kg−1.


## 3.2. Spectral Data Quality

Quality of the spectral data as acquired by laboratory measurements and the HyMap imaging spectrometer was evaluated visually by plotting the reflectance and first derivative of the reflectance of the soil probes and its corresponding sampling points. Figure3 displays the spectral signal of both measurement techniques for a sandy loam with 10.4 g OC kg−1. Maximum reflectance of the sampling point varied between 24% and 34% in case of the HyMap imaging spectrometer and the FieldSpec Pro, respectively. Higher reflectance of the laboratory measurements is probably due to oven-drying of the soil samples as soil moisture reduces the reflectance in the entire wavelength range (e.g., [59]). Soil roughness may also reduce reflectance in the field and is mitigated by sieving the samples before laboratory measurements. Despite the difference in absolute reflectance, the curves are similar in shape. Both show a pronounced absorption feature near 2200 nm. A less distinct and broader absorption feature can be observed near 500 nm and 650 nm which is clearly displayed by the first derivative. The first derivative of the reflectance discloses some differences between the spectra near 880 nm, around 1400 nm and 1900 nm, and at 2350 nm. Deviations around 1400 nm and 1900 nm are probably caused by atmospheric water absorption. It is presumed that the disparity near 880 nm is also related to atmospheric disturbances since a minor water absorption band occurs near 900 nm [14]. Minor differences near 2350 nm are likely a result of different surface conditions in the laboratory and the field where such weak absorption features are less pronounced due to interference with other factors.Reflectance (a) and first derivative of the reflectance (b) of one sampling point as measured under laboratory conditions (black solid line) and by the HyMap imaging spectrometer on May 27, 2005 (red solid line with diamonds indicating centre wavelengths). The soil was classified by soil analysis as a sandy loam with 10.4 g OC kg−1.
(a)
(b)


## 3.3. Image Data Quality

The quality of the atmospheric corrected HyMap images from 2004, 2005, and 2008 was assessed in two ways. First, the general consistency between the HyMap acquisitions of the different years was evaluated. This was done by comparing the mean image spectra of a relatively homogenous concrete area among each other and with field spectral measurements from August 2004. The overall shape of the reflectance curves of the three different years is in good agreement. A slight decline in reflectivity in the shortwave infrared region is observed for the spectra from 2005 and 2008 compared to 2004, whereas reflectivity in the visible and near-infrared region is highest for the spectra from 2008. But apart from a few minor deviations, all spectra range within the standard deviation of the ASD field spectral measurements (Figure4). Second, the relative accuracy of atmospheric correction between image strips of each year was assessed. Spectra of identical image pixels (within the given geometric accuracy) in the overlap of neighboring images were extracted. The ratio of the corresponding spectra of a soil surface and a former runway as a function of wavelength is shown in Figure 5. There is a good match for the runway in all years and for the soil surface in 2005 and 2008. The mean factor of 0.95 is within the expected accuracy of the atmospheric correction including minor BRDF effects during optimal illumination conditions. Increasing noise below 0.5 μm and above 2.4 μm is attributed to the lower signal-to-noise ratio of the HyMap sensor in these spectral regions. Artifacts near 1.4 μm and 1.7 to 1.9 μm are owed to atmospheric water absorption bands. The low ratio of the soil spectra in 2004 indicates significant differences in the overall reflectance of this surface type and occurs systematically in the overlap of all flight strips. Image acquisition parameters in Table 1 suggest that a high solar zenith and a small solar azimuth across the flight heading promote the appearance of BRDF effects on rough soil surfaces in the 2004 data set. This is confirmed by an across-track brightness gradient that occurs after atmospheric correction in all scenes of this year [48]. A comparison to a topographic map with a scale of 1 : 25.000 revealed an absolute geometric accuracy of less than 25 m. The relative geometric accuracy of the image data between the years amounts to two pixels in x or y direction, that is, 8 m.Figure 4
Comparison of mean reflectance spectra of a concrete area on a farm premise near the city of Kruckow acquired with an ASD Fieldspec Pro in 2004 and the HyMap sensor in 2004, 2005, and 2008.Ratio of spectra of identical image pixels (within the given geometric accuracy) in the overlap of neighboring scenes for a former runway and a soil surface.
(a)
(b)


## 3.4. Preparing Image Data for Spatially Explicit Parameter Estimation

In the central European agricultural region, bare soil is only encountered in early spring while fields are prepared for sowing of sugar beet, maize, potatoes, or summer grain, and in late summer after harvesting and seedbed preparation for winter grains (e.g., winter barley, winter wheat). Mitigating the impact of vegetation or excluding vegetated fields from the analysis is therefore a prerequisite for a successful application of spectroscopic approaches. In a preliminary study, Daughtry et al. [60] showed that photosynthetic active vegetation, nonphotosynthetic vegetation, and bare soils could be distinguished by a combination of vegetation indices, the cellulose absorption index (CAI) [61, 62], and the normalized difference vegetation index (NDVI). They obtained better results on AVIRIS data with the combination of indices than with spectral unmixing. In the present study, this approach was adopted to identify fields with bare soils and mask pixels with vegetation signal before image data analysis. Indices were computed for all HyMap images using wavelengths near 680 and 800 nm for the NDVI and wavelengths near 2000, 2100, and 2200 nm for CAI. Thresholds for green and dry vegetation were determined interactively and knowledge based and were set to 0.2 for NDVI and −1 to −2.5 for CAI. The quality of the soil mask was assessed using 200 ground truth regions of interest randomly distributed in equal parts over the agricultural area classified as vegetation and bare soil. The classification of the pixels was evaluated knowledge based with the help of the corresponding image spectra. There was almost no misinterpretation for pixels classified as vegetation. Disagreements in the soil class were mostly related to emerging green vegetation as some of the ground truth regions showed a small but discrete absorption feature near 680 nm but NDVI values below the assumed threshold. Accuracy statistics revealed a producer’s accuracy of 99% and 84% and a user’s accuracy of 86.1% and 98.8% for vegetation and soil, respectively. The kappa coefficient was 0.83 indicating a good overall agreement for the two observed classes.A maximum of 8.21% of the total area covered by the HyMap scenes from May 2005 was detected free of vegetation applying the combined threshold of CAI and NDVI. In the data sets from August 2004 and July 2008, only 2.2% and 1.29% of the total image area were identified as bare soils and made available for topsoil parameter estimation (Table3). Since the total image area covered by the HyMap data in 2008 exceeds the total area covered in 2004 by roughly one-third, the absolute number of image pixels available for image data analysis is within the same dimension.Table 3
Area free of vegetation in the HyMap images from 2004, 2005, and 2008 given in absolute values and percentage of the total image coverage.
Time of acquisition
Total number of image pixels
Image pixels with NDVI<0.2 and-2.5<CAI<-1.0
Proportion of image pixels with NDVI<0.2 and -2.5<CAI<-1.0(%)
09.08.2004
12.076.995
266.019
2.20
27.05.2005
12.495.804
1.026.187
8.21
29.07.2008
16.125.125
208.374
1.29In Figure6, CAI is plotted against NDVI for all 135 sampling points derived from HyMap images of the three years. In 2005, sixty-seven sampling points reside on fields without vegetation. The remaining sixty-eight points are situated on fields with green vegetation indicated by NDVI values above 0.6. In 2004 and 2008, respectively, the spectral signal of only 20 and 11 sampling points suggested the appearance of bare soil. Apart from a few data points with higher NDVI values, the majority of sampling points exhibit NDVI less than 0.2 but CAI values above zero which is clear evidence for nonphotosynthetic, dry vegetation. Due to the time of image data acquisition, these fields are either still grown with senescent crops or covered with stubbles which remain one the fields after harvesting.Figure 6
CAI versus NDVI derived from HyMap spectra of the years 2004, 2005, and 2008 of the 135 sampling points in the field.


## 3.5. Model Calibration and Validation

### 3.5.1. Calibration and Validation Schemes

Calibration and validation schemes were set up with respect to the availability of vegetation-free reference points in the HyMap data from 2004, 2005, and 2008. In total, three calibration data sets were generated. The first one comprises the spectra of the 20 sampling points identified in the image data from 2004 (Lab_2004, HyMap_2004). The second calibration set contains 38 out of 67 vegetation-free spectra in the image data from 2005 (Lab_2005, HyMap_2005). The third-one is made up of 11 spectra corresponding with the vegetation-free sampling points in the image data acquired during the flight campaign in 2008 (Lab_2008, HyMap_2008). With the laboratory spectroscopy, each of the calibration models were validated with the rest of the laboratory spectra withheld during calibration that are 115, 97 and 124 in number. For the calibration data set with 38 soil spectra, an additional validation was performed utilizing only those 29 samples which are located on fields with bare soil in the image data from 2005 (Lab_2005b). Validation of the HyMap models was confined to cross-validation in case of the HyMap_2004 and HyMap_2008 models due to the lack of vegetation-free reference points. The HyMap_2005 model was additionally validated with the 29 samples not used during calibration. Descriptive statistics of all calibration and validation data sets are summarized in Table4. As sampling has been carried out field-by-field with a mean within-field sampling distance of adjacent sampling points of 138 m, validation cannot be considered to be fully independent from the calibration data because of spatial autocorrelation. Regardless of how the validation subsample is selected, this can lead to an overestimation of the predictive accuracy of the models in particular for unsampled fields [63].Table 4
Descriptive statistics of the calibration and validation data sets for the prediction of clay and organic carbon content with laboratory and HyMap image spectra.
Calibration data sets
Laboratory
HyMap
Number of samples
Clay (g kg−1)
Organic carbon  (g kg−1)
Mean
SD
Min
Max
Mean
SD
Min
Max
Lab_2004
HyMap_2004
20
108.07
48.38
48.03
186.64
9.60
3.17
5.25
19.59
Lab_2005
HyMap_2005
38
99.12
50.05
38.45
216.18
9.16
2.95
5.01
17.48
Lab_2008
HyMap_2008
11
97.31
45.68
52.22
197.20
6.66
1.47
4.31
10.34
Validation data sets
Laboratory
HyMap
Number of samples
Clay (g kg−1)
Organic carbon (g kg−1)
Mean
SD
Min
Max
Mean
SD
Min
Max
Lab_2004
—
115
130.22
67.45
36.53
329.47
10.42
4.02
4.31
24.79
Lab_2005a
—
97
137.83
67.53
36.53
329.47
10.74
4.16
4.31
24.79
Lab_2005b
HyMap_2005
29
92.06
52.36
36.53
237.52
8.89
3.73
5.55
24.79
Lab_2004
—
124
129.57
66.26
36.53
329.47
10.62
3.90
5.01
24.79Except for the cal/val configuration with 38 and 29 samples, data range and standard deviation of the calibration data sets were considerably smaller than in the validation data sets. For the most part of spectroscopic studies, validation is performed either on data resembling the data range of the calibration data set or solely by cross-validation as model validity is restricted to the data range covered during calibration. Here, the described procedure was chosen to evaluate model performance outside the calibration data range since soil characteristics of the test site DEMMIN are not adequately represented by the small number of vegetation-free sampling points in the image data sets.


### 3.5.2. Laboratory Spectra

Table5 summarizes calibration and validation statistics of the PLSR models for the prediction of clay and organic carbon content. Tests with varying preprocessing methods showed that either absorbance or first derivative of the absorbance gave the best prediction results. Calibration data sets with a greater number of samples generally yield higher R2 and lower standard errors of calibration (SEC) than models with a fewer number of samples. For clay, R2 ranged from 0.70 to 0.90. SEC was lowest in the Lab_2005_Clay model (17.42 g kg−1) and highest in the Lab_2008_Clay model (35.68 g kg−1). Standard error of validation (SEP) varied from 19.48 to 35.55 g kg−1. Plotting measured against predicted clay contents reveals a close relationship and a homogeneous scatter along the 1 : 1 line for the Lab_2005_Clay model validated with 29 and 97 soil samples (Figure 7). For the Lab_2004_Clay model, measured and predicted values are generally in good agreement but soils with more than 200 g kg−1 clay, of which none was enclosed in the calibration data set, are slightly overestimated. Samples show a systematic, linear offset from the 1 : 1 line. However, deviations are within the error variance of the whole validation data set. The calibration model with the smallest number of samples (Lab_2008_Clay) caused large deviations for some validation samples with less than 200 g kg−1 clay but predictions of clay rich soils are less erroneous (Figure 7).Table 5
Calibration and validation performance statistics of laboratory spectra.
Parameter
Clay
Model name
Number of samples(cal/val)
Preprocessing
Number of LV
SDcal†
SDval†
R
²
SEC
SD/SEC
SEP
RPD
Lab_2004_Clay
20/115
Absorption
6
48.21
75.21
0.76
29.30
1.65
26.59
2.54
Lab_2005_Clay
38/97
First derivative of absorption
5
49.32
70.06
0.90
17.42
2.87
20.77
3.25
Lab_2005_Clay
38/29
First derivative of absorption
5
49.32
48.08
0.90
17.42
2.87
19.48
2.69
Lab_2008_Clay
11/124
Absorption
5
39.92
69.18
0.70
35.68
1.28
35.55
1.86
Parameter
Organic carbon
Model name
Number of samples(cal/val)
Preprocessing
Number of LV
SDcal†
SDval†
R
2
SEC
SD/SEC
SEP
RPD
Lab_2004_OC
20/115
Absorption
2
2.53
3.77
0.51
2.36
1.34
2.36
1.71
Lab_2005_OC
38/97
Absorption
6
2.87
3.90
0.88
1.11
2.67
1.92
2.16
Lab_2005_OC
38/29
Absorption
6
2.87
2.63
0.88
1.11
2.67
1.94
1.93
Lab_2008_OC
11/124
First derivative of absorption
4
1.49
4.12
0.78
0.97
1.52
2.46
1.58
†SD: standard deviation of estimated parameter values in (g kg−1), R²: coefficient of determination, SEC: standard error of calibration in (g kg−1), SEP: standard error of validation in (g kg−1), RPD: ratio of performance to deviation.Plots of measured versus predicted values for clay (gkg−1) in the calibration (black symbol, ×) and validation (red symbol, +) data sets based on laboratory spectra. The respective calibration and validation data schemes (ncal./nval.) are marked on the plots. For further explanation see text.
(a)
(b)
(c)
(d)Calibrations based on at least 20 samples had RPD values clearly above 2, indicating excellent prediction models. They can be assigned to group A in the classification of Chang et al. [56]. Reducing the number of samples in the calibration data set to just 11 spectra, results in less predictive power with an RPD of 1.86. Such models belong to group B and may be improved by different calibration strategies [56]. Similar results on the prediction of clay from laboratory spectra are reported, for example, from Waiser et al. [64], Chang et al. [56], and Cozzolino and Morón [65]. The latter performed a cross-validation with a modified PLSR approach on spectra from Uruguayan soils and achieved a value of 2.7 (ratio of standard deviation by standard error of cross validation). Waiser et al. [64] tested four different soil preparation techniques including air-dried ground samples from Texas and achieved RPD values between 1.95 and 3.51 with PLSR. Principal component regression and cross validation were applied by Chang et al. [56]. They obtained an RPD of 1.7. In contrast to our work, the range of clay contents in the calibration and validation data sets in those studies were similar. The small range of clay concentrations in the soil samples of the calibration data sets varying from about 40 to 220 g kg−1 clay did hardly affect the validity of our models even for samples with higher clay contents. But a very small number of samples hamper the generation of a solid calibration model for the prediction of clay content using VNIR spectroscopy.R
2 values of the PLSR models on organic carbon ranged from 0.51 to 0.88 with SEC between 1.11 and 2.36 g kg−1. SEP varied between 1.92 and 2.46 g kg−1. The best results were again achieved providing 38 samples in the calibration data set (Lab_2005_OC) having the highest R2 and lowest standard errors. Calibration models for organic carbon exhibit a comparable predictive ability as those for clay. The maximum ratio of standard deviation to SEC is 2.67 for the Lab_2005_OC model and minimum 1.34 for the Lab_2004_OC model. In contrast to the prediction of clay, the distribution of data points in the scatter plots for the validation data sets indicated less precise models (Figure 8). A good agreement is observed between measured and predicted values in the validation data sets for the Lab_2005_OC model for samples with less than 15.0 g kg−1 OC. Samples with higher concentrations show a trend of underestimation. Validating the models with 20 and 11 calibration samples reveals a great variety in the measured versus predicted plots which seems to increase with higher organic carbon contents. This decrease in prediction accuracy is reflected by lower RPD values. With the exception of the Lab_2005_OC model which has a RPD close to 2 or above 2, all models must be categorized in group B following the classification of Chang et al. [56]. The precision obtained for OC in our study is within the precision obtained by previous authors (e.g., [34, 51, 66]). Stevens et al. [34], for example, generated a PLSR model with 117 soil samples from agricultural fields in Belgian Lorraine with a RPD of 2.11. Volkan Bilgili et al. [66] predicted several soil properties of a 32 ha field in the north of Turkey among them soil organic matter (SOM). In their study, a higher data range in the validation data set (3.9–68.7 g kg−1 SOM) than in the calibration data set (5.0–33.8 g kg−1 SOM) did not result in a decrease of the models predictive ability, neither for PLSR (RPD: 1.93) nor with multivariate adaptive regression splines (MARS) (RPD: 1.90). This may be due to the great number of samples in their cal/val sets (153/359).Plots of measured versus predicted values for organic carbon (g kg−1) in the calibration (black symbol, ×) and validation (red symbol, +) data sets based on laboratory spectra. The respective calibration and validation data schemes (ncal./nval.) are marked on the plots. For further explanation see text.
(a)
(b)
(c)
(d)Resampling laboratory spectra to HyMap, spectral resolution had only little affect on the prediction accuracy of clay and none on the prediction accuracy of organic carbon (results not shown). Similarly, Gomez et al. [21] reported little influence on prediction accuracy while resampling spectra to Hyperion resolution.Variable importance in the projection (VIP) and PLS-regression coefficients (b-coefficients) provide information about the weight of the x variables (wavelengths) in the PLSR model. Large b-coefficients indicate wavelengths which are relevant in the modeling of Y (soil parameter). High VIP values point out wavelengths which are essential for the modeling of both, Y and X [67]. Thus, x variables with a high VIP and a high b-coefficient indicate important wavelengths and allow insights on the physical basis of the predictions. VIP values and b-coefficients of all laboratory models are plotted in Figure 9. Curve shapes of the Lab_2005_Clay model and the Lab_2008_OC model have a different appearance compared to the rest due to the applied spectral preprocessing (first derivative of absorbance). For clay, VIP and b-coefficients indicate an emphasis of wavelengths between 2000 and 2500 nm and in particular around 2200 nm across all models. In fact, clay minerals are spectrally active in this part of the spectrum [12]. Absorption features near 2200 to 2300 nm produced by combination vibrations involving an OH stretch and metal-OH bend are diagnostic for all clay minerals. The exact position of the absorption features varies depending on the metal which is bond to the hydroxyl group [13, 14]. Furthermore, VIP and b-coefficients suggest that wavelengths in the visible part of the spectrum near 500 nm and 650 nm contain important information on clay minerals or rather clay related features as the spectral activity of clays is confined to the infrared portion of the spectrum. Iron, either as main constituent or associated with the structure of clay minerals, shows spectral features in the VNIR due to electronic transitions of iron cations [14]. Hunt and Salisbury [12], for example, attributed weak absorption bands near 500 nm to ferrous iron in montmorillonites. In illites, absorption near 650 nm was observed by Clark et al. [13] who ascribed this feature to ferric iron. The very close correlation of iron and clay (Table 2) in the soil samples of the test site DEMMIN encourages the assumption that iron has an influence on the PLS regression for clay.Variable importance in the projection (VIP) andb-coefficients of PLSR calibration models based on laboratory spectra using 20, 38, and 11 soil samples for clay (left side) and organic carbon (right side).
(a)
(b)
(c)
(d)
(e)
(f)VIP andb-coefficient curves of the organic carbon models exhibit less consistency than for clay. Whereas high loadings can be observed for the Lab_2005_OC and Lab_2008_OC model between 2000 and 2500 nm similar to clay, this region seems of little importance in the Lab_2004_OC model although a small accentuation occurs near 2200 nm. In the latter, generally high VIP values and b-coefficients prevail in the 500 to 900 nm region. In the Lab_2005_OC and Lab_2008_OC models wavelengths near 650 nm and 500 to 600 nm, respectively, are emphasized. As organic carbon in soils consists of very complex chemical compounds having numerous, frequently overlapping absorption features [68, 69], a spectral feature assignment to specific wavelengths is difficult [7, 70]. In addition, good correlations between clay and organic carbon or iron and organic carbon hamper the elucidation of the causal relationships between model parameters and its physical background.


### 3.5.3. HyMap Image Spectra

In a first attempt, it was tried to set up a stable calibration model using the image data from 2005 and validating this model with the remaining sampling points from 2004, 2005, and 2008. Model calibration and validation with image spectra from 2005 yielded good results with RPD of 2.70 and 1.80 for clay and organic carbon, respectively (Table7). Validation with image spectra from different flight campaigns, that is, 2004 and 2008, reached RPD below 1.0 indicating poor results with no or little predictive ability especially for clay (Table 6). This is the case although soil clay and OC contents of the respective sampling points are within the data range of the calibration model and spectral data quality assessment showed a good overall agreement between years (cf. Section 3.2). Experiments with laboratory spectra revealed that the composition of the validation data set did only slightly impair the model’s predictive ability. The problem with empirical calibration is that coefficients are precisely adjusted to the spectral properties of the calibration input data. Even small changes in spectral data properties hamper a reasonable application. Laboratory spectra are measured under standard conditions where variations in spectral properties are solely a consequence of variations in physical and chemical characteristics. Under field conditions, soil status (e.g., roughness, crusting, moisture) varies with time. Further variations may be introduced through radiometric and atmospheric correction (see Figure 4). Thus, it was not possible to set up one stable calibration model for the estimation of soil parameters from image data but data set specific calibrations seem to be crucial to retrieve reasonable prediction accuracies.Table 6
Validation performance statistics of HyMap_2005 model transfers to HyMap spectra from different acquisition years.
Transfer of HyMap_2005 model to
Clay
Organic carbon
SEP (g kg−1)
RPD
SEP (g kg−1)
RPD
HyMap spectra 2004
175.04
0.28
4.82
0.66
HyMap spectra 2008
239.66
0.19
1.58
0.93Table 7
Calibration and validation performance statistics of HyMap image spectra.
Parameter
Clay
Model name
Number of samples(cal/val)
Preprocessing
Number of LV
SDcal†
SDval†
R
²
SEC
SD/SEC
SEP
RPD
HyMap_2004_Clay
20
Absorption
3
44.72
—
0.80
23.39
2.07
—
—
HyMap_2005_Clay
38/29
First derivative of absorption
3
47.18
50.49
0.85
20.34
2.46
19.41
2.70
HyMap_2008_Clay
11
Absorption
3
42.36
—
0.92
15.99
2.86
—
—
Parameter
Organic carbon
Model name
Number of samples(cal/val)
Preprocessing
Number of LV
SDcal†
SDval†
R
²
SEC
SD/SEC
SEP
RPD
HyMap_2004_OC
20
First derivative of absorption
3
2.70
—
0.62
2.13
1.48
—
—
HyMap_2005_OC
38/29
First derivative of absorption
2
2.58
2.20
0.71
1.64
1.80
2.07
1.80
HyMap_2008_OC
11
First derivative of absorption
5
1.20
—
0.45
1.61
0.91
—
—
†SD: standard deviation of estimated parameter values in (g kg−1), R2: coefficient of determination, SEC: standard error of calibration in (g kg−1), SEP: standard error of validation in (g kg−1), RPD: ratio of performance to deviation.Results of data set specific PLS regressions with the HyMap spectra from 2004, 2005, and 2008 are listed in Table7. Calibration and validation were performed as described in Section 3.5.1 (Table 4). As with laboratory spectra prediction accuracy was higher for clay than for organic carbon but a high number of samples in the calibration data set did not necessarily yield the best results. For clay, R2 ranged from 0.80 to 0.92 with SEC as low as 15.99 to 23.39 g kg−1. SD divided by SEC is highest for the HyMap spectra from 2008 reaching a value of 2.86 (HyMap_2008_clay). The calibration model for the image spectra from 2005 is characterized by a SD/SEC ratio of 2.46 (HyMap_2005_clay). Independent validation with the remaining 29 data points resulted in a RPD of 2.70 confirming a very high prediction accuracy of this model similar to the image spectra from 2008. The lowest prediction ability is achieved with the calibration model based on image spectra from 2004 with a SD/SEC ratio of 2.07 (HyMap_2004_clay). Measured versus predicted values suggest that samples with low clay contents (<80 g kg−1) may not be adequately predicted applying this model as variability of estimations is very high for samples with more or less identical clay contents (Figure 10). According to Chang et al. [56], all three models can be classified as group A with an excellent ability to predict soil clay content from hyperspectral image data. A direct comparison of our results to the prediction accuracy obtained by other authors is difficult since none of them reports RPD or SD/SEC values. Standard errors may only be compared between studies if soils have an analogous range of clay content. This is the case with the work of Selige et al. [19] who reported a cross-validated R2 of 0.65 with a RMSE of 38 g kg−1 for MLR and a R2 of 0.71 with a RMSE of 42 g kg−1 for PLSR for estimations of clay from HyMap data in an agricultural area in the German federal state Saxony-Anhalt. Although their samples reached a maximum of 260 g kg−1 clay, these error measures clearly indicate that prediction accuracy for clay from HyMap image data in DEMMIN is superior.Plots of measured versus predicted values for clay (gkg−1) in the calibration (black symbol, ×) and validation (red symbol, +) data sets based on HyMap image spectra. The respective calibration and validation data schemes (ncal./nval.) are marked on the plots. For further explanation see text.
(a)
(b)
(c)R
2 obtained from model calibration for the prediction of OC content from HyMap spectra are 0.71 at most and 0.45 at least with SEC between 1.61 and 2.13 g kg−1. Prediction ability of the PLSR model based on the 2005 image data is acceptable with a SD/SEC ratio and an RPD of 1.80. Being restricted to 29 validation samples, no clear conclusion can be drawn on the model’s behaviour above 20 g kg−1 OC. Similar to the corresponding laboratory cal/val configuration (Lab_2005_Clay, nval=29), a strong underestimation of one data point with about 25 g kg−1 OC is observed (Figures 8 and 11). This gives reason to interpret that the HyMap_2005_OC model is subject to a nonlinear trend just as examined in the validation of the laboratory model with 97 samples. Calibration models based on image data from 2004 and 2008 exhibit SD/SEC values of 1.49 and 0.91, respectively. Thus, models for the prediction of OC based on image spectra belong to group B (HyMap_2004_OC, HyMap_2005_OC) and group C (HyMap_2008_OC) after Chang et al. [56]. Prediction accuracy obtained by PLS regression on HyMap images from 2004 and 2005 match results from Stevens et al. [34] and Stevens et al. [23] applying image data of the sensors AHS-160 (RPD: 1.40) and CASI (RPD: 1.86). Stevens et al. [20] showed that regrouping a global cal/val data set by soil type into several smaller cal/val data sets could improve prediction accuracy about two times with an RPD as high as 3.15. However, this approach is only feasible if the number of reference points is large and if auxiliary information on soil type is provided for delineation of soil types in the image data in order to correctly apply the different calibration models.Plots of measured versus predicted values for organic carbon (gkg−1) in the calibration (black symbol, ×) and validation (red symbol, +) data sets based on HyMap image spectra. The respective calibration and validation data schemes (ncal./nval.) are marked on the plots. For further explanation see text.
(a)
(b)
(c)Since laboratory measurements are acquired under standardized conditions, it was expected that calibrations with HyMap image spectra results in a deterioration of prediction accuracy. Lagacherie et al. [22], for example, attributed a decrease in model performance from laboratory to imaging spectroscopy to radiometric and wavelength calibration uncertainties of the HyMap sensor as well as atmospheric effects. Gomez et al. [21] assumed that a higher signal-to-noise ratio and the large spatial support area (30×30 m²) of Hyperion images caused a decline in prediction ability. Other factors degrading the model performance may be variable soil moisture content [23] or soil surface conditions [71] as vegetation affected spectra were removed before the image data analysis. In fact, higher standard errors were observed in this study for the HyMap_2005_OC and the HyMap_2008_OC models. Apart from the factors mentioned above, the interannual variability of soil organic carbon [40, 41] may have influenced the spectral calibrations for these two years. Since image data was analyzed with reference data from 2006, possible variations of soil organic carbon caused by the specific temperature and soil moisture regime of each year could not be taken into account. Prediction ability of the HyMap_2005_Clay model remained more or less constant (Tables 5 and 7). Judging from the SD/SEC values model performance increased when applying the image data from 2004 (HyMap_2004_OC, HyMap_2004_Clay). Accuracy also improved when estimating clay using the data from 2008 although it decreased when estimating organic carbon with the same data set. A reasonable explanation for this behavior was not found. In any case, it must be stated that only a reverification of the predictive accuracy with totally independent reference samples can reveal further insight and the true predictive ability of the presented PLSR models. This is of particular importance for fields not sampled during field work. To obtain a first idea of possible implications caused by spatial autocorrelation, a calibration and validation scheme for the HyMAP spectra from 2005 was tested in which all samples from three fields were retained from calibration at a time and used for model validation only. This was done for five different cal/val constellations. Standard errors of validation (SEP) retrieved this way for clay were less than the SEP of the HyMAP_2005_Clay model in three cases. In two of these cal/val schemes, SEP increased by about 11% and 15%, that is, a maximum SEP of 22.3 g kg−1 was obtained. For OC, SEP increased by about 9% and 29% for two cal/val schemes. The top most SEP was 2.68 g kg−1. In the remaining test cases considered, SEP was less than the SEP of the HyMAP_2005_OC model. These results corroborate the urgent need for independent validation. Although an increase in standard errors by about 30% may still seem acceptable, parameter estimations for soils with low clay or OC contents would become very inaccurate.VIP andb-coefficient curves of the PLSR models for clay based on HyMap spectra presented in Table 7 signalize an importance of wavelength between 2000 and 2500 nm just as observed for the laboratory models. Again special emphasis is given to the 2200 nm band. This is less pronounced in the HyMap_2008_Clay model where wavelengths in the visible and near-infrared part of the spectrum seem to play a major role. Other than with the laboratory spectra, wavelengths between 900 and 1300 nm contain high values in this model with a slight accentuation near 1100 nm. This inflection can also be seen in the Lab_2008_Clay model and may be associated with absorption caused by ferrous iron in illite minerals [13]. Increased VIP and b-coefficients in all three models can also be found in the 500 to 650 nm region and are probably linked to iron as discussed with the laboratory spectra. For organic carbon, several wavelengths in the visible part of the spectrum are stressed. Peaks occur at 550, 680, and 770 nm. Likewise, bands near 2000, 2200, and 2350 nm carry high loadings in all three models similar to the laboratory models. As has already been stated in Section 3.5.2, there are noticeable analogies to the clay model parameters.


## 3.5.1. Calibration and Validation Schemes

Calibration and validation schemes were set up with respect to the availability of vegetation-free reference points in the HyMap data from 2004, 2005, and 2008. In total, three calibration data sets were generated. The first one comprises the spectra of the 20 sampling points identified in the image data from 2004 (Lab_2004, HyMap_2004). The second calibration set contains 38 out of 67 vegetation-free spectra in the image data from 2005 (Lab_2005, HyMap_2005). The third-one is made up of 11 spectra corresponding with the vegetation-free sampling points in the image data acquired during the flight campaign in 2008 (Lab_2008, HyMap_2008). With the laboratory spectroscopy, each of the calibration models were validated with the rest of the laboratory spectra withheld during calibration that are 115, 97 and 124 in number. For the calibration data set with 38 soil spectra, an additional validation was performed utilizing only those 29 samples which are located on fields with bare soil in the image data from 2005 (Lab_2005b). Validation of the HyMap models was confined to cross-validation in case of the HyMap_2004 and HyMap_2008 models due to the lack of vegetation-free reference points. The HyMap_2005 model was additionally validated with the 29 samples not used during calibration. Descriptive statistics of all calibration and validation data sets are summarized in Table4. As sampling has been carried out field-by-field with a mean within-field sampling distance of adjacent sampling points of 138 m, validation cannot be considered to be fully independent from the calibration data because of spatial autocorrelation. Regardless of how the validation subsample is selected, this can lead to an overestimation of the predictive accuracy of the models in particular for unsampled fields [63].Table 4
Descriptive statistics of the calibration and validation data sets for the prediction of clay and organic carbon content with laboratory and HyMap image spectra.
Calibration data sets
Laboratory
HyMap
Number of samples
Clay (g kg−1)
Organic carbon  (g kg−1)
Mean
SD
Min
Max
Mean
SD
Min
Max
Lab_2004
HyMap_2004
20
108.07
48.38
48.03
186.64
9.60
3.17
5.25
19.59
Lab_2005
HyMap_2005
38
99.12
50.05
38.45
216.18
9.16
2.95
5.01
17.48
Lab_2008
HyMap_2008
11
97.31
45.68
52.22
197.20
6.66
1.47
4.31
10.34
Validation data sets
Laboratory
HyMap
Number of samples
Clay (g kg−1)
Organic carbon (g kg−1)
Mean
SD
Min
Max
Mean
SD
Min
Max
Lab_2004
—
115
130.22
67.45
36.53
329.47
10.42
4.02
4.31
24.79
Lab_2005a
—
97
137.83
67.53
36.53
329.47
10.74
4.16
4.31
24.79
Lab_2005b
HyMap_2005
29
92.06
52.36
36.53
237.52
8.89
3.73
5.55
24.79
Lab_2004
—
124
129.57
66.26
36.53
329.47
10.62
3.90
5.01
24.79Except for the cal/val configuration with 38 and 29 samples, data range and standard deviation of the calibration data sets were considerably smaller than in the validation data sets. For the most part of spectroscopic studies, validation is performed either on data resembling the data range of the calibration data set or solely by cross-validation as model validity is restricted to the data range covered during calibration. Here, the described procedure was chosen to evaluate model performance outside the calibration data range since soil characteristics of the test site DEMMIN are not adequately represented by the small number of vegetation-free sampling points in the image data sets.


## 3.5.2. Laboratory Spectra

Table5 summarizes calibration and validation statistics of the PLSR models for the prediction of clay and organic carbon content. Tests with varying preprocessing methods showed that either absorbance or first derivative of the absorbance gave the best prediction results. Calibration data sets with a greater number of samples generally yield higher R2 and lower standard errors of calibration (SEC) than models with a fewer number of samples. For clay, R2 ranged from 0.70 to 0.90. SEC was lowest in the Lab_2005_Clay model (17.42 g kg−1) and highest in the Lab_2008_Clay model (35.68 g kg−1). Standard error of validation (SEP) varied from 19.48 to 35.55 g kg−1. Plotting measured against predicted clay contents reveals a close relationship and a homogeneous scatter along the 1 : 1 line for the Lab_2005_Clay model validated with 29 and 97 soil samples (Figure 7). For the Lab_2004_Clay model, measured and predicted values are generally in good agreement but soils with more than 200 g kg−1 clay, of which none was enclosed in the calibration data set, are slightly overestimated. Samples show a systematic, linear offset from the 1 : 1 line. However, deviations are within the error variance of the whole validation data set. The calibration model with the smallest number of samples (Lab_2008_Clay) caused large deviations for some validation samples with less than 200 g kg−1 clay but predictions of clay rich soils are less erroneous (Figure 7).Table 5
Calibration and validation performance statistics of laboratory spectra.
Parameter
Clay
Model name
Number of samples(cal/val)
Preprocessing
Number of LV
SDcal†
SDval†
R
²
SEC
SD/SEC
SEP
RPD
Lab_2004_Clay
20/115
Absorption
6
48.21
75.21
0.76
29.30
1.65
26.59
2.54
Lab_2005_Clay
38/97
First derivative of absorption
5
49.32
70.06
0.90
17.42
2.87
20.77
3.25
Lab_2005_Clay
38/29
First derivative of absorption
5
49.32
48.08
0.90
17.42
2.87
19.48
2.69
Lab_2008_Clay
11/124
Absorption
5
39.92
69.18
0.70
35.68
1.28
35.55
1.86
Parameter
Organic carbon
Model name
Number of samples(cal/val)
Preprocessing
Number of LV
SDcal†
SDval†
R
2
SEC
SD/SEC
SEP
RPD
Lab_2004_OC
20/115
Absorption
2
2.53
3.77
0.51
2.36
1.34
2.36
1.71
Lab_2005_OC
38/97
Absorption
6
2.87
3.90
0.88
1.11
2.67
1.92
2.16
Lab_2005_OC
38/29
Absorption
6
2.87
2.63
0.88
1.11
2.67
1.94
1.93
Lab_2008_OC
11/124
First derivative of absorption
4
1.49
4.12
0.78
0.97
1.52
2.46
1.58
†SD: standard deviation of estimated parameter values in (g kg−1), R²: coefficient of determination, SEC: standard error of calibration in (g kg−1), SEP: standard error of validation in (g kg−1), RPD: ratio of performance to deviation.Plots of measured versus predicted values for clay (gkg−1) in the calibration (black symbol, ×) and validation (red symbol, +) data sets based on laboratory spectra. The respective calibration and validation data schemes (ncal./nval.) are marked on the plots. For further explanation see text.
(a)
(b)
(c)
(d)Calibrations based on at least 20 samples had RPD values clearly above 2, indicating excellent prediction models. They can be assigned to group A in the classification of Chang et al. [56]. Reducing the number of samples in the calibration data set to just 11 spectra, results in less predictive power with an RPD of 1.86. Such models belong to group B and may be improved by different calibration strategies [56]. Similar results on the prediction of clay from laboratory spectra are reported, for example, from Waiser et al. [64], Chang et al. [56], and Cozzolino and Morón [65]. The latter performed a cross-validation with a modified PLSR approach on spectra from Uruguayan soils and achieved a value of 2.7 (ratio of standard deviation by standard error of cross validation). Waiser et al. [64] tested four different soil preparation techniques including air-dried ground samples from Texas and achieved RPD values between 1.95 and 3.51 with PLSR. Principal component regression and cross validation were applied by Chang et al. [56]. They obtained an RPD of 1.7. In contrast to our work, the range of clay contents in the calibration and validation data sets in those studies were similar. The small range of clay concentrations in the soil samples of the calibration data sets varying from about 40 to 220 g kg−1 clay did hardly affect the validity of our models even for samples with higher clay contents. But a very small number of samples hamper the generation of a solid calibration model for the prediction of clay content using VNIR spectroscopy.R
2 values of the PLSR models on organic carbon ranged from 0.51 to 0.88 with SEC between 1.11 and 2.36 g kg−1. SEP varied between 1.92 and 2.46 g kg−1. The best results were again achieved providing 38 samples in the calibration data set (Lab_2005_OC) having the highest R2 and lowest standard errors. Calibration models for organic carbon exhibit a comparable predictive ability as those for clay. The maximum ratio of standard deviation to SEC is 2.67 for the Lab_2005_OC model and minimum 1.34 for the Lab_2004_OC model. In contrast to the prediction of clay, the distribution of data points in the scatter plots for the validation data sets indicated less precise models (Figure 8). A good agreement is observed between measured and predicted values in the validation data sets for the Lab_2005_OC model for samples with less than 15.0 g kg−1 OC. Samples with higher concentrations show a trend of underestimation. Validating the models with 20 and 11 calibration samples reveals a great variety in the measured versus predicted plots which seems to increase with higher organic carbon contents. This decrease in prediction accuracy is reflected by lower RPD values. With the exception of the Lab_2005_OC model which has a RPD close to 2 or above 2, all models must be categorized in group B following the classification of Chang et al. [56]. The precision obtained for OC in our study is within the precision obtained by previous authors (e.g., [34, 51, 66]). Stevens et al. [34], for example, generated a PLSR model with 117 soil samples from agricultural fields in Belgian Lorraine with a RPD of 2.11. Volkan Bilgili et al. [66] predicted several soil properties of a 32 ha field in the north of Turkey among them soil organic matter (SOM). In their study, a higher data range in the validation data set (3.9–68.7 g kg−1 SOM) than in the calibration data set (5.0–33.8 g kg−1 SOM) did not result in a decrease of the models predictive ability, neither for PLSR (RPD: 1.93) nor with multivariate adaptive regression splines (MARS) (RPD: 1.90). This may be due to the great number of samples in their cal/val sets (153/359).Plots of measured versus predicted values for organic carbon (g kg−1) in the calibration (black symbol, ×) and validation (red symbol, +) data sets based on laboratory spectra. The respective calibration and validation data schemes (ncal./nval.) are marked on the plots. For further explanation see text.
(a)
(b)
(c)
(d)Resampling laboratory spectra to HyMap, spectral resolution had only little affect on the prediction accuracy of clay and none on the prediction accuracy of organic carbon (results not shown). Similarly, Gomez et al. [21] reported little influence on prediction accuracy while resampling spectra to Hyperion resolution.Variable importance in the projection (VIP) and PLS-regression coefficients (b-coefficients) provide information about the weight of the x variables (wavelengths) in the PLSR model. Large b-coefficients indicate wavelengths which are relevant in the modeling of Y (soil parameter). High VIP values point out wavelengths which are essential for the modeling of both, Y and X [67]. Thus, x variables with a high VIP and a high b-coefficient indicate important wavelengths and allow insights on the physical basis of the predictions. VIP values and b-coefficients of all laboratory models are plotted in Figure 9. Curve shapes of the Lab_2005_Clay model and the Lab_2008_OC model have a different appearance compared to the rest due to the applied spectral preprocessing (first derivative of absorbance). For clay, VIP and b-coefficients indicate an emphasis of wavelengths between 2000 and 2500 nm and in particular around 2200 nm across all models. In fact, clay minerals are spectrally active in this part of the spectrum [12]. Absorption features near 2200 to 2300 nm produced by combination vibrations involving an OH stretch and metal-OH bend are diagnostic for all clay minerals. The exact position of the absorption features varies depending on the metal which is bond to the hydroxyl group [13, 14]. Furthermore, VIP and b-coefficients suggest that wavelengths in the visible part of the spectrum near 500 nm and 650 nm contain important information on clay minerals or rather clay related features as the spectral activity of clays is confined to the infrared portion of the spectrum. Iron, either as main constituent or associated with the structure of clay minerals, shows spectral features in the VNIR due to electronic transitions of iron cations [14]. Hunt and Salisbury [12], for example, attributed weak absorption bands near 500 nm to ferrous iron in montmorillonites. In illites, absorption near 650 nm was observed by Clark et al. [13] who ascribed this feature to ferric iron. The very close correlation of iron and clay (Table 2) in the soil samples of the test site DEMMIN encourages the assumption that iron has an influence on the PLS regression for clay.Variable importance in the projection (VIP) andb-coefficients of PLSR calibration models based on laboratory spectra using 20, 38, and 11 soil samples for clay (left side) and organic carbon (right side).
(a)
(b)
(c)
(d)
(e)
(f)VIP andb-coefficient curves of the organic carbon models exhibit less consistency than for clay. Whereas high loadings can be observed for the Lab_2005_OC and Lab_2008_OC model between 2000 and 2500 nm similar to clay, this region seems of little importance in the Lab_2004_OC model although a small accentuation occurs near 2200 nm. In the latter, generally high VIP values and b-coefficients prevail in the 500 to 900 nm region. In the Lab_2005_OC and Lab_2008_OC models wavelengths near 650 nm and 500 to 600 nm, respectively, are emphasized. As organic carbon in soils consists of very complex chemical compounds having numerous, frequently overlapping absorption features [68, 69], a spectral feature assignment to specific wavelengths is difficult [7, 70]. In addition, good correlations between clay and organic carbon or iron and organic carbon hamper the elucidation of the causal relationships between model parameters and its physical background.


## 3.5.3. HyMap Image Spectra

In a first attempt, it was tried to set up a stable calibration model using the image data from 2005 and validating this model with the remaining sampling points from 2004, 2005, and 2008. Model calibration and validation with image spectra from 2005 yielded good results with RPD of 2.70 and 1.80 for clay and organic carbon, respectively (Table7). Validation with image spectra from different flight campaigns, that is, 2004 and 2008, reached RPD below 1.0 indicating poor results with no or little predictive ability especially for clay (Table 6). This is the case although soil clay and OC contents of the respective sampling points are within the data range of the calibration model and spectral data quality assessment showed a good overall agreement between years (cf. Section 3.2). Experiments with laboratory spectra revealed that the composition of the validation data set did only slightly impair the model’s predictive ability. The problem with empirical calibration is that coefficients are precisely adjusted to the spectral properties of the calibration input data. Even small changes in spectral data properties hamper a reasonable application. Laboratory spectra are measured under standard conditions where variations in spectral properties are solely a consequence of variations in physical and chemical characteristics. Under field conditions, soil status (e.g., roughness, crusting, moisture) varies with time. Further variations may be introduced through radiometric and atmospheric correction (see Figure 4). Thus, it was not possible to set up one stable calibration model for the estimation of soil parameters from image data but data set specific calibrations seem to be crucial to retrieve reasonable prediction accuracies.Table 6
Validation performance statistics of HyMap_2005 model transfers to HyMap spectra from different acquisition years.
Transfer of HyMap_2005 model to
Clay
Organic carbon
SEP (g kg−1)
RPD
SEP (g kg−1)
RPD
HyMap spectra 2004
175.04
0.28
4.82
0.66
HyMap spectra 2008
239.66
0.19
1.58
0.93Table 7
Calibration and validation performance statistics of HyMap image spectra.
Parameter
Clay
Model name
Number of samples(cal/val)
Preprocessing
Number of LV
SDcal†
SDval†
R
²
SEC
SD/SEC
SEP
RPD
HyMap_2004_Clay
20
Absorption
3
44.72
—
0.80
23.39
2.07
—
—
HyMap_2005_Clay
38/29
First derivative of absorption
3
47.18
50.49
0.85
20.34
2.46
19.41
2.70
HyMap_2008_Clay
11
Absorption
3
42.36
—
0.92
15.99
2.86
—
—
Parameter
Organic carbon
Model name
Number of samples(cal/val)
Preprocessing
Number of LV
SDcal†
SDval†
R
²
SEC
SD/SEC
SEP
RPD
HyMap_2004_OC
20
First derivative of absorption
3
2.70
—
0.62
2.13
1.48
—
—
HyMap_2005_OC
38/29
First derivative of absorption
2
2.58
2.20
0.71
1.64
1.80
2.07
1.80
HyMap_2008_OC
11
First derivative of absorption
5
1.20
—
0.45
1.61
0.91
—
—
†SD: standard deviation of estimated parameter values in (g kg−1), R2: coefficient of determination, SEC: standard error of calibration in (g kg−1), SEP: standard error of validation in (g kg−1), RPD: ratio of performance to deviation.Results of data set specific PLS regressions with the HyMap spectra from 2004, 2005, and 2008 are listed in Table7. Calibration and validation were performed as described in Section 3.5.1 (Table 4). As with laboratory spectra prediction accuracy was higher for clay than for organic carbon but a high number of samples in the calibration data set did not necessarily yield the best results. For clay, R2 ranged from 0.80 to 0.92 with SEC as low as 15.99 to 23.39 g kg−1. SD divided by SEC is highest for the HyMap spectra from 2008 reaching a value of 2.86 (HyMap_2008_clay). The calibration model for the image spectra from 2005 is characterized by a SD/SEC ratio of 2.46 (HyMap_2005_clay). Independent validation with the remaining 29 data points resulted in a RPD of 2.70 confirming a very high prediction accuracy of this model similar to the image spectra from 2008. The lowest prediction ability is achieved with the calibration model based on image spectra from 2004 with a SD/SEC ratio of 2.07 (HyMap_2004_clay). Measured versus predicted values suggest that samples with low clay contents (<80 g kg−1) may not be adequately predicted applying this model as variability of estimations is very high for samples with more or less identical clay contents (Figure 10). According to Chang et al. [56], all three models can be classified as group A with an excellent ability to predict soil clay content from hyperspectral image data. A direct comparison of our results to the prediction accuracy obtained by other authors is difficult since none of them reports RPD or SD/SEC values. Standard errors may only be compared between studies if soils have an analogous range of clay content. This is the case with the work of Selige et al. [19] who reported a cross-validated R2 of 0.65 with a RMSE of 38 g kg−1 for MLR and a R2 of 0.71 with a RMSE of 42 g kg−1 for PLSR for estimations of clay from HyMap data in an agricultural area in the German federal state Saxony-Anhalt. Although their samples reached a maximum of 260 g kg−1 clay, these error measures clearly indicate that prediction accuracy for clay from HyMap image data in DEMMIN is superior.Plots of measured versus predicted values for clay (gkg−1) in the calibration (black symbol, ×) and validation (red symbol, +) data sets based on HyMap image spectra. The respective calibration and validation data schemes (ncal./nval.) are marked on the plots. For further explanation see text.
(a)
(b)
(c)R
2 obtained from model calibration for the prediction of OC content from HyMap spectra are 0.71 at most and 0.45 at least with SEC between 1.61 and 2.13 g kg−1. Prediction ability of the PLSR model based on the 2005 image data is acceptable with a SD/SEC ratio and an RPD of 1.80. Being restricted to 29 validation samples, no clear conclusion can be drawn on the model’s behaviour above 20 g kg−1 OC. Similar to the corresponding laboratory cal/val configuration (Lab_2005_Clay, nval=29), a strong underestimation of one data point with about 25 g kg−1 OC is observed (Figures 8 and 11). This gives reason to interpret that the HyMap_2005_OC model is subject to a nonlinear trend just as examined in the validation of the laboratory model with 97 samples. Calibration models based on image data from 2004 and 2008 exhibit SD/SEC values of 1.49 and 0.91, respectively. Thus, models for the prediction of OC based on image spectra belong to group B (HyMap_2004_OC, HyMap_2005_OC) and group C (HyMap_2008_OC) after Chang et al. [56]. Prediction accuracy obtained by PLS regression on HyMap images from 2004 and 2005 match results from Stevens et al. [34] and Stevens et al. [23] applying image data of the sensors AHS-160 (RPD: 1.40) and CASI (RPD: 1.86). Stevens et al. [20] showed that regrouping a global cal/val data set by soil type into several smaller cal/val data sets could improve prediction accuracy about two times with an RPD as high as 3.15. However, this approach is only feasible if the number of reference points is large and if auxiliary information on soil type is provided for delineation of soil types in the image data in order to correctly apply the different calibration models.Plots of measured versus predicted values for organic carbon (gkg−1) in the calibration (black symbol, ×) and validation (red symbol, +) data sets based on HyMap image spectra. The respective calibration and validation data schemes (ncal./nval.) are marked on the plots. For further explanation see text.
(a)
(b)
(c)Since laboratory measurements are acquired under standardized conditions, it was expected that calibrations with HyMap image spectra results in a deterioration of prediction accuracy. Lagacherie et al. [22], for example, attributed a decrease in model performance from laboratory to imaging spectroscopy to radiometric and wavelength calibration uncertainties of the HyMap sensor as well as atmospheric effects. Gomez et al. [21] assumed that a higher signal-to-noise ratio and the large spatial support area (30×30 m²) of Hyperion images caused a decline in prediction ability. Other factors degrading the model performance may be variable soil moisture content [23] or soil surface conditions [71] as vegetation affected spectra were removed before the image data analysis. In fact, higher standard errors were observed in this study for the HyMap_2005_OC and the HyMap_2008_OC models. Apart from the factors mentioned above, the interannual variability of soil organic carbon [40, 41] may have influenced the spectral calibrations for these two years. Since image data was analyzed with reference data from 2006, possible variations of soil organic carbon caused by the specific temperature and soil moisture regime of each year could not be taken into account. Prediction ability of the HyMap_2005_Clay model remained more or less constant (Tables 5 and 7). Judging from the SD/SEC values model performance increased when applying the image data from 2004 (HyMap_2004_OC, HyMap_2004_Clay). Accuracy also improved when estimating clay using the data from 2008 although it decreased when estimating organic carbon with the same data set. A reasonable explanation for this behavior was not found. In any case, it must be stated that only a reverification of the predictive accuracy with totally independent reference samples can reveal further insight and the true predictive ability of the presented PLSR models. This is of particular importance for fields not sampled during field work. To obtain a first idea of possible implications caused by spatial autocorrelation, a calibration and validation scheme for the HyMAP spectra from 2005 was tested in which all samples from three fields were retained from calibration at a time and used for model validation only. This was done for five different cal/val constellations. Standard errors of validation (SEP) retrieved this way for clay were less than the SEP of the HyMAP_2005_Clay model in three cases. In two of these cal/val schemes, SEP increased by about 11% and 15%, that is, a maximum SEP of 22.3 g kg−1 was obtained. For OC, SEP increased by about 9% and 29% for two cal/val schemes. The top most SEP was 2.68 g kg−1. In the remaining test cases considered, SEP was less than the SEP of the HyMAP_2005_OC model. These results corroborate the urgent need for independent validation. Although an increase in standard errors by about 30% may still seem acceptable, parameter estimations for soils with low clay or OC contents would become very inaccurate.VIP andb-coefficient curves of the PLSR models for clay based on HyMap spectra presented in Table 7 signalize an importance of wavelength between 2000 and 2500 nm just as observed for the laboratory models. Again special emphasis is given to the 2200 nm band. This is less pronounced in the HyMap_2008_Clay model where wavelengths in the visible and near-infrared part of the spectrum seem to play a major role. Other than with the laboratory spectra, wavelengths between 900 and 1300 nm contain high values in this model with a slight accentuation near 1100 nm. This inflection can also be seen in the Lab_2008_Clay model and may be associated with absorption caused by ferrous iron in illite minerals [13]. Increased VIP and b-coefficients in all three models can also be found in the 500 to 650 nm region and are probably linked to iron as discussed with the laboratory spectra. For organic carbon, several wavelengths in the visible part of the spectrum are stressed. Peaks occur at 550, 680, and 770 nm. Likewise, bands near 2000, 2200, and 2350 nm carry high loadings in all three models similar to the laboratory models. As has already been stated in Section 3.5.2, there are noticeable analogies to the clay model parameters.


## 3.6. Spatially Explicit Parameter Estimation

In Figure12 results of the multi-annual soil parameter estimation using hyperspectral image data and PLS regression is exemplarily displayed for the parameter clay. Flight strips of each flight campaign were processed individually. Mosaics were created subsequently. Since statistical analysis was done for soils from arable land only, quantitative mapping was restricted to agricultural crop land by means of a vector layer with field boundaries. White regions represent nonarable land (e.g., settlements, forest, and grassland). Grey regions indicate arable land completely covered with vegetation which was masked using CAI and NDVI (cf. Section 3.4). Supplementary information is given on the field specific model predictive ability classified after Chang et al. [56] including details on the respective calibration data range (Figure 12(b)). Careful documentation of error measures and calibration data range is considered essential as it differs with multi-annual image data. Figure 12(c) provides direct access to the field-specific image data source. As expected, clay content is less than 125 g kg−1 within most of the fields. Despite this predominance of sandy loam and loamy sand, the map reveals a typical characteristic of the DEMMIN area which is the occurrence of partly large local heterogeneities of the soil type. A single field of 40 ha near the village Heydenhof is an example (Figure 13) with clay contents close to 300 g kg−1. Similarly, organic carbon contents are above average in the eastern part of the field reaching up to 17.5 g kg−1. The spatial pattern of clay and organic carbon content in this field shows little congruence. This is in agreement with a moderate correlation coefficient (Table 2) between the two soil properties and considerable deviations that occur for a few sandy soils with very high organic carbon contents but very little clay. Spatial pattern of clay and OC in this field is most likely not related to topography as the area is very flat.Figure 12
Map of soil clay content on agricultural fields (a) derived by multi-annual hyperspectral image data analysis. Map of RPD values (b) indicate field-specific model predictive ability after Chang et al. [56]. A field assignment to the year of acquisition is illustrated in Map (c).Figure 13
Map of soil clay and organic carbon content derived from HyMap image data (May 27, 2005) and PLSR for a field near the village of Heydenhof. The location of this field on the test site DEMMIN is marked by a red flag in Figure1.Performing multi-annual estimations of soil parameters offers the opportunity to compare parameter estimations from different years which provides an additional measure of the model’s predictive ability. Because of the time of image data acquisition, only two fields could be identified where parameters were estimated in either 2004 and 2005 or 2005 and 2008. No overlap was detected for all three years. Figure14 displays the difference image of estimated clay and organic carbon contents based on HyMap images from 2004 and 2005. During the overflight in 2004, this field was partly covered with crop residues which are more or less arranged in line according to the tillage direction. Furthermore, the state of tillage differed within the field. In the eastern part, vegetation residues almost completely cover the soil surface. These pixels are masked (white). On average, predicted clay contents in 2005 are 3.42 g kg−1 below the estimates in 2004. 50% of all data points show differences between −14.19 and 9.11 g kg−1 clay. Mean differences of predicted OC contents are close to 0 whereby half the number of data points differ between −0.11 and 0.57 g kg−1 OC. The spatial pattern of both difference images is similar. Higher negative differences appear with both parameters in the upper left of the field (dark blue). A small area in the centre of the field exhibits higher positive deviations (dark orange-red). Differences reach extreme values for a few isolated pixels. They are attributed to artifacts or problems caused by mixed pixels. RMSE values calculated for both parameters are equal to or less than standard errors of model calibration (Table 7). This result is encouraging as it confirms statements made on model accuracies and emphasizes the potential of imaging spectroscopy to map soil parameters. However, clay and OC contents on this particular field are quite low. As laboratory models showed an increase in variance with increasing OC concentrations, the actual RMSE calculated for OC may be higher. This was not the case for clay but only further investigations may provide additional evidence.Figure 14
Difference image of estimated clay and organic carbon content based on HyMap images from 2004 and 2005 for an agricultural field near the city of Marienfelde. Q1:  lower quartile, Q3: upper quartile. Mean deviation (RMSE) calculated from estimations in 2005 and 2004. The location of this field on the test site DEMMIN is marked by a red flag in Figure1.Reflectance of most natural surfaces varies depending on the viewing and solar illumination conditions. This fact is described by the bidirectional reflectance distribution function (BDRF) [72]. Evaluating prediction accuracy in the overlap area of the flight strips per year may give an indication for inaccuracies related to these reflectance anisotropies. Stevens et al. [20] presumed that difficulties to produce robust calibrations were partly caused by strong backward- or forward-scattering at the image edges due to the large field of view of the AHS sensor (90°). Computing RMSE revealed substantially higher errors in the overlap area for the 2004 data set than was expected from the validation, respectively, cross-validation of the models (Table 8). For clay, RMSE was about 3 times larger than SEC using reflectance transformed to absorbance. Then applying reflectance transformed to the first derivative of absorbance RMSE was about 1.4 times larger than SEC. Effects were less pronounced for organic carbon but RMSE still exceeded SEC by a factor of 1.6 using absorbance. In contrast, RMSE were either only slightly increased or even less than standard errors of validation for HyMap spectra from 2005 and 2008 independent from the preprocessing method applied. These results provide strong evidence that BRDF effects have negative influence on the performance of the models based on image data from 2004. As no negative impact was observed in the two other data sets, this effect seems solely due to the unfavorable viewing and illumination conditions during the overflight in 2004, that is, a large zenith and a low azimuth angle (Table 1), but not to the sensors field of view in principle. Negative effects in the off-nadir area can be diminished applying the first derivate even though a somewhat higher standard error must be accepted.Table 8
Mean deviation (RMSE) of predicted clay and organic carbon contents in the overlap area of the HyMap flight strips and standard errors of PLSR models generated using the spectra transformed to absorbance (Abs) and spectra transformed to first derivative of absorbance (Abs_Deriv1). Bold numbers indicate model parameters presented in Table7.
Day of acquisition
Number of pixels estimated
Clay
Organic carbon
RMSEoverlap
SEP* or SEC**
RMSEoverlap
SEP* or SEC**
Abs
Abs_Deriv1
Abs
Abs_Deriv1
Abs
Abs_Deriv1
Abs
Abs_Deriv1
August 9, 2004
36015
69.63
38.99
23.3
9
*
*
28.35**
3.86
2.20
2.35**
2.1
3
*
*
May 27, 2005
260206
22.46
17.03
20.96*
19.4
1
*
2.34
0.60
3.13*
2.0
7
*
July 29, 2008
60308
12.95
15.39
15.99
*
*
22.89**
2.21
0.79
1.88**
1.6
1
*
*
SEC: standard error of calibration in (g kg−1) and SEP: standard error of validation in (g kg−1).


## 4. Conclusions

Multi-annual spatially explicit soil parameter estimation was performed using three HyMap images of the test site DEMMIN and PLSR. It was shown that this approach provides an opportunity to subsequently map soil parameters on agricultural fields despite long periods of vegetation coverage. Image data applied within this study was acquired in May, July, and August. Whereas the spring image offered a good chance to estimate soil parameters on fields grown with root crops or summer grains, images acquired in summer were less suited. More than 97.8% of the images were covered with dry or green vegetation at this time of the year. Having three hyperspectral images at our disposal, estimations on clay and organic carbon content could be performed for about 10% of the total area covered. The presented approach will be facilitated by future hyperspectral satellite missions such as EnMAP [73] since there is no restriction to a very few if not only one overflight per year as is the case with airborne sensors. Thus, requirements by researchers interested in vegetation as well as soil parameters will be met. Due to the high seasonal and interannual variations of soil organic carbon [36, 41], model calibrations for this parameter may be negatively affected if the time of acquisition of the image data and the reference data deviates. Problems could be solved by sampling the soil more or less at the same time as the flight campaign. For soil parameters with little temporal variability such as clay, this will not be required.Calibration models for clay content based on image spectra showed excellent prediction ability with a maximum standard error of 23.39 g kg−1 clay. Clay estimations were generally more precise than predictions of organic carbon content whose standard error was 2.13 g kg−1 OC at most. For both soil parameters, validation with soil samples withheld from calibration approved prediction errors of the models although for two image data sets validation was confined to leave-one-out cross-validation due to the little number of reference samples. A comparison of predictions from different years suggests that there is a very good agreement between estimations of different years and thus confirms validation results. Investigations under laboratory conditions provided an additional indication that standard errors for clay may also apply to soils with concentrations outside the calibration data range. However, despite these promising results, further validation with truly independent samples will be crucial to reveal the model’s predictive ability for the entire data set including fields without reference information. If this is known, soil parameter maps derived from imaging spectroscopy may serve as a basis for precision agriculture, environmental modeling, or carbon sequestration studies, for example.Furthermore, it was shown that off-nadir parameter estimations were not affected by BRDF effects if image data was acquired under optimal viewing and illumination conditions and in flat terrain. But if acquisition took place while zenith angle was high and azimuth angle was low inaccuracies enlarged up to about 3 times. This negative effect could be mitigated using an appropriate spectral preprocessing such as the first derivative even though this resulted in a certain loss of the overall prediction accuracy.Results presented here demonstrate that high resolution soil parameter maps of agricultural fields can be derived on a regional scale with multi-annual imaging spectroscopy data. Prediction accuracy both for laboratory and imaging spectroscopy does not come up to conventional analytical techniques in the laboratory. Model errors are at least about two times higher than SES and SEL taking into account that model errors include measurement errors of reference values. The great benefit of imaging spectroscopy is the provision of information on soil parameters as a two-dimensional array rather than single data points whereby inaccuracies associated with spatial interpolation are avoided. Because of frequent plowing on croplands, surface reflectance measured by hyperspectral instruments does not only reflect the characteristics of the upper five centimeters but does very likely represent the soil characteristics of the upper 20 to 30 cm. As this layer is most prone to degradation, imaging spectroscopy can provide valuable information for soil conservation measures on agricultural land.



---

*Source: 868090-2012-10-10.xml*