A New Concept of Soil Line Retrieval from Landsat 8 Images for Estimating Plant Biophysical Parameters

Extraction of vegetation information from remotely sensed images has remained a long-term challenge due to the influence of soil background. To reduce this effect, the slope and intercept of the soil line (SL) should be known to calculate SL-related vegetation indices (VIs). These VIs can be used to estimate the biophysical parameters of agricultural crops. However, it is a difficult task to retrieve the SL parameters under the vegetation canopy. A feasible method for retrieving these parameters involves extracting the bottom boundary line in two-dimensional spectral spaces (i.e., red and near-infrared bands). In this study, the slope and intercept of the SL was extracted from Landsat 8 OLI images of a test site in northeastern Germany. Different statistical methods, including the Red-NIRmin method, quantile regression method (using a floating tau with the smallest p-value), and a new approach proposed in this paper using a fixed quantile tau known as the diffuse non-interceptance (DIFN) value, were applied to retrieve the SL parameters. The DIFN value describes the amount of light visible below the canopy that reaches the soil surface. Therefore, this value can be used as a threshold for retrieving the bottom soil line. The simulated SLs were compared with actual ones extracted from ground truth data, as recorded by a handheld spectrometer, and were also compared with the SL retrieved from bare soil pixels of the Landsat 8 image collected after harvest. Subsequently, the SL parameters were used to separately estimate the dry biomasses of winter wheat (Triticum aestivum L.), barley (Hordeum vulgare L.), and canola (Brassica napus L.) at the local and field scales using different SL-related vegetation indices. The SL can be retrieved more accurately at the local scale compared with the field scale, and its simulation can be critical in the field due to significant differences from the actual SL. Moreover, the slope and intercept of the simulated SLs found using the floating and fixed quantile tau (slope ≈ 1.1 and intercept ≈ 0.05) show better agreement with the actual SL parameters (slope ≈ 1.2 and intercept ≈ 0.03) in the late growing stages (i.e., end of ripening and senescence stages) of crops. The slope and intercept of the soil line extracted from bare soil pixels of the Landsat 8 OLI data after harvest (slope = 1.3, intercept = 0.03, and R2 = 0.94) are similar to those of the simulated SL. The correlation coefficient (R2) of the simulated SLs are greater than 0.97 during different growing stage and all of the SL parameters are statistically significant (p < 0.05) at the local scale. The results also imply the need for different vegetation indices to best retrieve the crop biomass depending on the growing stage, but relatively small differences in performances were observed in this study. Remote Sens. 2016, 8, 738; doi:10.3390/rs8090738 www.mdpi.com/journal/remotesensing Remote Sens. 2016, 8, 738 2 of 23


Introduction
In optical remote sensing (RS) for agriculture, extraction of the vegetation canopy information has remained a long-term challenge due to background effects.Baret et al. (1993) [1] acquired important information related to this subject by creating the so-called soil line (SL).In fact, soil reflectance varies with time and place depending on many factors, such as the soil type [2], soil physical and chemical properties, soil moisture [3], soil mineral composition (i.e., inorganic matter or ash content), soil texture [4], surface roughness [5] and vegetation residuals [1].
Mineral soil components usually increase the reflectance from the visible to the shortwave infrared wavelengths.However, organic matter might also indirectly disturb the spectral response of soil based on its structure [2,[6][7][8][9].Additionally, increasing soil moisture decreases the reflectance intensity [9].Another factor that can alter the soil spectral properties is the surface roughness.Experimental results show that the soil reflectance increases with decreasing particle size for a given type of soil [10,11].
The SL can be extracted via regression analysis of the soil reflectance (in a two-dimensional space) of a group of samples in the lab.In 2004, Fox [8] developed an automated method for extracting the SLs from remotely sensed images of cropland using a digital aerial camera that fits a linear regression by deriving a set of minimum near-infrared digital numbers across the red and NIR bands from bare soil pixels as the initial SL.Subsequently, they used an iterative process to delete the points (i.e., contradicting SL pixels) that differ greatly from the initial SL to obtain the final one.Later, Xu and Guo [12] extended Fox's methods to grassland using two different statistical methods (Red-NIR min ) and quantile regression.
To calculate the SL-related VIs in RS applications, the presence of bare soil pixels or other prior knowledge of the soil is necessary to calculate the slope and intercept of the SL.However, this information is not always available in real applications.Signals in the optical domain of the electromagnetic spectrum have rather shallow penetration ability, and therefore, extracting the soil reflectance spectra under the vegetation canopy has always been challenging.In these situations, extraction of the bottom boundary line in two-dimensional spectral spaces can be a feasible solution for extraction of the SL and its parameters.Nonetheless, few studies thus far have retrieved this bottom boundary line and provide the parameters (i.e., slope and intercept) for the SL concept, e.g., [8,12].The SL and its parameters are highly important for many applications in RS science.
In 1976, Kauth and Thomas [13] described the famous triangular tasseled cap in Red-NIR space for simulation of the SL using MSS (multi-spectral scanner) data.The SL is a hypothetical line that can explain the variation of the soil spectrum in multi-or hyper-spectral RS images.Usually, this line can be designated using the spectral reflectance of the soil background in the red and near-infrared bands [13,14].In 1977, Richardson and Wiegand [15] introduced the perpendicular vegetation index (PVI) to reduce the disturbance of the soil background when extracting vegetation signals from Landsat MSS images based on the concept of SL.Since then, additional vegetation indices, that is, the so-called SL related vegetation indices, were developed based on the slope and intercept of the SL to reduce the background influences when retrieving the canopy characteristics [16][17][18].These vegetation indices can evaluate different characteristics of the target using the Euclidean or angular distance from the SL [16,19,20].
To this end, this study addresses the hypothesis that the bottom boundary line of the R-NIR scatter plot can be extracted more accurately using the new proposed concept of tau (fixed and floating tau), especially at the end of the growing season and these methods can be used as a complementary tool for SL parameters retrieval.Therefore, this study develops the new concept of the floating and fixed tau method using quantile regression.In addition, this study compares the aforementioned method with the (Red-NIR min ) method and the soil spectra collected in the field (i.e., ground truth data) to extract SL from Landsat 8 images in three croplands (barley, winter wheat, and canola).The simulated SLs are also compared with the SL extracted from the bare soil pixels of Landsat 8 data after harvest.More specifically, the objectives of this research are: (1) to simulate the SL using the aforementioned methods from a two-dimensional spectral space using the NIR and Red bands of Landsat 8 OLI data and (2) to evaluate the efficiency of the SL parameters (i.e., slope and intercept) extracted by all of these approaches to retrieve the dry biomass of different crops using different SL-related vegetation indices.

Test Site and Study Area
The study area is located in northeastern Germany, Durable Environmental Multidisciplinary Monitoring Information Network (DEMMIN).Different crops are planted in this area (e.g., barley, wheat, and canola).The dominant crops are winter types covering almost 60% of the fields in the area [33].Soil substrates are dominated by loamy sands and sandy loams alternation with pure sand patches or clayey areas [33].Moreover, the region belongs to the temperate zone.Maximum temperature was reached with approximately 32 • C at the end of July and the minimum value with approximately 0 • C were reached in the beginning of May 2013 during the field campaign.In 2013, the mean annual air and soil surface temperature in the study area were 8.88 • C and 8.95 • C, respectively, and the mean precipitation was 41.55 mm, as recorded by the closest weather station (approximately 2 km away).Within this test site, a study area of 32 km 2 that includes three winter wheat (225 ha), barley (117 ha) and canola fields (25 ha) (Figure 1) was chosen for ground truth data collection.For additional info, please refer to [34].To this end, this study addresses the hypothesis that the bottom boundary line of the R-NIR scatter plot can be extracted more accurately using the new proposed concept of tau (fixed and floating tau), especially at the end of the growing season and these methods can be used as a complementary tool for SL parameters retrieval.Therefore, this study develops the new concept of the floating and fixed tau method using quantile regression.In addition, this study compares the aforementioned method with the (Red-NIRmin) method and the soil spectra collected in the field (i.e., ground truth data) to extract SL from Landsat 8 images in three croplands (barley, winter wheat, and canola).The simulated SLs are also compared with the SL extracted from the bare soil pixels of Landsat 8 data after harvest.More specifically, the objectives of this research are: (1) to simulate the SL using the aforementioned methods from a two-dimensional spectral space using the NIR and Red bands of Landsat 8 OLI data and (2) to evaluate the efficiency of the SL parameters (i.e., slope and intercept) extracted by all of these approaches to retrieve the dry biomass of different crops using different SL-related vegetation indices.

Test Site and Study Area
The study area is located in northeastern Germany, Durable Environmental Multidisciplinary Monitoring Information Network (DEMMIN).Different crops are planted in this area (e.g., barley, wheat, and canola).The dominant crops are winter types covering almost 60% of the fields in the area [33].Soil substrates are dominated by loamy sands and sandy loams alternation with pure sand patches or clayey areas [33].Moreover, the region belongs to the temperate zone.Maximum temperature was reached with approximately 32 °C at the end of July and the minimum value with approximately 0 °C were reached in the beginning of May 2013 during the field campaign.In 2013, the mean annual air and soil surface temperature in the study area were 8.88 °C and 8.95 °C, respectively, and the mean precipitation was 41.55 mm, as recorded by the closest weather station (approximately 2 km away).Within this test site, a study area of 32 km 2 that includes three winter wheat (225 ha), barley (117 ha) and canola fields (25 ha) (Figure 1) was chosen for ground truth data collection.For additional info, please refer to [34].

Field Data and Laboratory Analysis
The whole monitoring campaign for the crops was conducted for 31 weeks from 17 April until 13 November 2013 at weekly intervals.During this period, winter wheat, barley, and canola were monitored throughout their entire growing stages.At each field trip, two random center on the field and five sampling locations using five quadrats (50 cm × 50 cm) around each center were chosen.The distance to the center was listed and the GPS coordinates were recorded with a handheld Trimble device for mapping of data in a Geographic Information System (GIS).
The vegetation parameters such as fresh plant biomass (kg/m 2 ) were collected from inside of the quadrats.The crops were cut off as close to the ground as possible (approximately 10-12 samples per field during each field trips).These vegetation samples were taken to the laboratory in sealable plastic bags and then the fresh biomass was measured for all components (stems, leaves, and shoots if applicable).Dry biomass was determined after drying at 75 • C for at least 48 h.
During each field expedition, after removing the biomass, the reflectance of the bare soil surface was examined using a handheld JAZ Spectrometer (Ocean Optics) (see Figure 2) (i.e., approximately 10-12 measurements collected in each crop field during each field trip).The wavelength range of this device is approximately 340-1100 nm with an optical resolution of approximately 0.3-10.0nm.We used a WS-1 white spectralon (i.e., reflectance standard) for 99% reflectivity.Soil samples were also collected using a hand sledge, and 5.6-cm-diameter cylinders/rings (i.e., approximately 10-12 samples were collected in each crop field during each field trip) from the top 0-10 cm of the soil surface inside the quadrats (Figure 2).To estimate the ash content (mineral residues, value between 0 and 1) of the soil, the samples were taken to the laboratory and burned at 550 • C for 4 to 5 h in a furnace, and the residual minerals (inorganic material) were determined.This analysis was performed based on A5TM D 2974-87 Standard Test Methods for Moisture, Ash, and Organic Matter of Peat and Other Organic Soils.Soil samples were also used in soil surface texture analyses, which refer to the size of the particles that make up the soil.Soil texture analysis for each crop field was conducted in the laboratory with a laser particle sizer ("ANALYSETTE 22", MicroTec plus, FRITSCH GmbH, Idar-Oberstein, Germany).Because the soil surface texture does not change throughout the growing season unless it is altered by tillage, texture analysis was performed once at the beginning of the growing season.
The leaf area index (LAI) was measured using a LAI-2200 plant canopy analyzer (LI-CORE, Lincoln, NE, USA) close to the quadrats.The 270 • view cap was used to hide the operator from the sensor.Normally, the LAI can be computed using two readings, the first above (A) and the second below the canopy (B).To determine the number of below (B) readings necessary for 95% confidence that the true LAI mean is within ±10% of the measured LAI, it is necessary to obtain a LAI reading based on 6 below (B) readings that include both the thinnest and densest portions of the canopy.By applying this procedure, the standard error (SEL) of the LAI was calculated and was divided by the LAI (i.e., SEL/LAI), and a table from LI-CORE was used to determine the number of B readings for all consequent readings.After collecting the LAI values in the field, post-processing of data was performed using the software from LI-CORE (i.e., FV2200).Moreover, the DIFN was calculated by integrating the gap fraction (GAPS) to yield a value that indicates the fraction of the sky that is not blocked by foliage.This value ranges from 0 (no sky visible to the sensor) to 1 (no foliage visible to the sensor) and can be obtained from the leaf area index (LAI) file.To obtain the maximum foliage of the crop in the field of view of the LAI sensor, the device was placed on the soil surface for each and every measurement.

Satellite Data
Landsat 8 Operational Land Imager (OLI) consists of nine spectral bands with a spatial resolution of 30 m for bands 1 to 7 and 9 and 15 m for panchromatic band (band 8) and temporal resolution of 16 days.A series of Landsat 8 OLI images of the study area and the different fields was acquired in 2013 on 6 May, 7 June, 14 June, 9 July, 16 July, 10 August, and 26 August, with no or small amounts (less than 10%) of cloud content.Because the images had already been geometrically corrected, radiometric and atmospheric correction was performed in ENVI 5.1 (Exelis) using the FLAASH module.FLAASH uses a standard equation for spectral radiance at a sensor pixel L that applies to the solar wavelength range (thermal emission is neglected) and flat Lambertian materials or their equivalents.The equation is given as follows: where ρ is the pixel surface reflectance, ρ e is the average surface reflectance for the pixel and a surrounding region, S is the spherical albedo of the atmosphere, L a is the radiance backscattered by the atmosphere, and A and B are coefficients that depend on atmospheric and geometric conditions but not on the surface [35].The values of A, B, S, and L a are determined from MODTRAN4 calculations that use the viewing and solar angles as well as the mean surface elevation of the measurement and assume a certain model atmosphere, aerosol type, and visible range [36].
After atmospheric correction, a subset of the study area and the winter wheat, barley and canola fields were extracted from airborne data with notably high spatial resolution (30 cm) using a shapefile created in ARCGIS software.The clouds, shadows (if applicable), and kettle holes were visually identified and manually masked if necessary.The projection of all the images was WGS_1984_UTM_Zone_33N.For a better visualization and overview of the ground truth and satellite data, Table 1 summarizes the variable names, number of samples for ground truth data collection, and satellite images used in this study.

(Red-NIR min ) Regression Method
The common method for simulating the bottom line of the scatter plot is to fit the simple linear regression of a set of points characterized by the minimum near-infrared (NIR) reflectance value within the certain red (R) reflectance values, which is known as the (Red-NIR min ) method [12].To fit the line in the two-dimensional space of the red and near infrared bands, the entire range of reflectance values in red were divided into several portions with equal intervals (e.g., 0 < R ≤ 0.001, 0.001 < R ≤ 0.002, etc.).Subsequently, one point with a minimum NIR value was selected in each portion.Finally, a general linear model was fitted to those points, and the best simulated SL was retrieved with the maximum coefficient of determination (R 2 ) and minimum p-value.The intervals might be changed (i.e., decreased or increased) gradually to retrieve the best-fit line.

Quantile Regression Method
With developments in linear programming [37], the quantile regression method is able to compete with the traditional least squares methods in terms of computational effort [38].Quantile regression is a method for estimating the functional relations between variables for all portions of a probability distribution [38].Instead of focusing on the mean, quantile regression gives the relationship between at least one covariate and the conditional median or other quantiles of the distribution of the response variable, which is usually between zero and one [39][40][41][42].Because this method fits the regression line to a portion of the distribution [43], it can be used to retrieve boundary lines in two-dimensional spectral space [44,45].
To retrieve the bottom boundary (i.e., SL) by using the quantile regression technique in NIR versus Red scatter plots, the quantile tau is initially set to 0.00001 to extract the bottom SL, and the tau is gradually changed to depict the best bottom line with the minimum p-value and maximum t-stat value.The quantile is set by "the number of points below or on the bottom line" divided by "the total number of points".It should be mentioned that by setting the quantile to zero, all of the points are above the simulated line.In the second approach, the DIFN value was used as a fixed tau to reduce the difficulties in choosing the best tau.As mentioned previously, the DIFN is calculated by integrating the gap fractions and represents the amount of diffuse short-wave radiation that is not intercepted or blocked by the plant canopy.The DIFN is an accurate and direct measure of porosity [46] and thus represents the plant canopy structure because it is influenced by the canopy LAI and mean tilt angle (MTA) [47].In other words, this value measures the amount of light that is visible below the canopy [48] that reaches the soil surface.Therefore, this value was used as the threshold (i.e., quantile tau) for specifying the bottom boundary line of the scatter plot using the quantile regression method.In this work, it can be assumed that the signals that progress upward from the LAI sensor (soil surface) to the satellite sensor (i.e., Landsat 8) are the same as the signals that travel downward to the LAI device.This value was used as a proxy for vegetation cover in this study.In other words, one minus DIFN is equal to the vegetation cover because DIFN shows the fraction of the soil that is not blocked by foliage and can be observed from the remote sensor (Figure 3).For example, a DIFN value of 0.2 indicates that 20% of the signal reaches the soil, and therefore, we can assume that the lower 20% of the data distribution in the R-NIR scatter plot represents the soil.

Quantile Regression Method
With developments in linear programming [37], the quantile regression method is able to compete with the traditional least squares methods in terms of computational effort [38].Quantile regression is a method for estimating the functional relations between variables for all portions of a probability distribution [38].Instead of focusing on the mean, quantile regression gives the relationship between at least one covariate and the conditional median or other quantiles of the distribution of the response variable, which is usually between zero and one [39][40][41][42].Because this method fits the regression line to a portion of the distribution [43], it can be used to retrieve boundary lines in two-dimensional spectral space [44,45].
To retrieve the bottom boundary (i.e., SL) by using the quantile regression technique in NIR versus Red scatter plots, the quantile tau is initially set to 0.00001 to extract the bottom SL, and the tau is gradually changed to depict the best bottom line with the minimum p-value and maximum t-stat value.The quantile is set by "the number of points below or on the bottom line" divided by "the total number of points".It should be mentioned that by setting the quantile to zero, all of the points are above the simulated line.In the second approach, the DIFN value was used as a fixed tau to reduce the difficulties in choosing the best tau.As mentioned previously, the DIFN is calculated by integrating the gap fractions and represents the amount of diffuse short-wave radiation that is not intercepted or blocked by the plant canopy.The DIFN is an accurate and direct measure of porosity [46] and thus represents the plant canopy structure because it is influenced by the canopy LAI and mean tilt angle (MTA) [47].In other words, this value measures the amount of light that is visible below the canopy [48] that reaches the soil surface.Therefore, this value was used as the threshold (i.e., quantile tau) for specifying the bottom boundary line of the scatter plot using the quantile regression method.In this work, it can be assumed that the signals that progress upward from the LAI sensor (soil surface) to the satellite sensor (i.e., Landsat 8) are the same as the signals that travel downward to the LAI device.This value was used as a proxy for vegetation cover in this study.In other words, one minus DIFN is equal to the vegetation cover because DIFN shows the fraction of the soil that is not blocked by foliage and can be observed from the remote sensor (Figure 3).For example, a DIFN value of 0.2 indicates that 20% of the signal reaches the soil, and therefore, we can assume that the lower 20% of the data distribution in the R-NIR scatter plot represents the soil.
A simple linear relationship was constructed between the ground truth data and the corresponding SL-related VI pixel values from the satellite images.These values were used to evaluate the relationships between the SL-related VIs and the dry biomass of crops and to calibrate the models.To examine the autocorrelation between the data, the Durbin-Watson test for autocorrelated errors was performed on the dataset.Moreover, four diagnostic plots, including residuals vs. fitted, normal Q-Q, scale location, and residuals vs. leverage, were studied to evaluate the assumptions of the regression models.Based on the aforementioned methods, the assumptions of the regression models appear to be upheld, including the non-existence of any autocorrelation between data.
SL-related vegetation indices (Table 2) [15,16,32,[51][52][53] were used to evaluate the efficiency of the slopes and intercepts of the SLs extracted in this study using different methods by comparing their ability to estimate the dry biomass from the field campaign.The correlation coefficient (R 2 ) and root mean squared error (RMSE) were also reported for assessment of goodness of fit.Although R 2 (i.e., correlation coefficient) is often used to assess model accuracy, RMSE is often used to quantify model precision [54].

Results
Ash content analysis shows high levels of inorganic matter (i.e., ash content) with relatively constant amounts and small fluctuations throughout the winter wheat, barley, and canola life cycle.The amounts of ash content are approximately 0.96 and 0.97, 0.94 (corresponding to 4%, 3%, and 6% organic matter in the soil) for the winter wheat, barley, and canola fields, respectively.As observed from Figure 4, soil from the canola field has the highest organic matter (i.e., lowest ash content), whereas soil from the barley field shows the highest ash content (i.e., approximately 0.97 ash or 0.03 organic matter).
The DIFN values are presented in Figure 5 for winter wheat, barley, canola and their combinations using pooled data (using all measurements).In the beginning of the growing season, the DIFN values are high, corresponding to a low LAI.As the plants grow, the LAI values increase, leading to decreasing DIFN values because the portion of sky that is visible to the sensor becomes smaller, and vice versa.At the end of the growing season, the LAI again decreases, which leads to a slight increase in DIFN values.The dry biomass also shows a relatively increasing monotonic trend during the entire growing season.The error bars of the dry biomass depict the standard deviation of the data.As mentioned, the bare soil reflectances were collected using a handheld spectrometer in the fields after biomass removal (see Figure 2).These reflectances were spectrally re-sampled according to the Landsat 8 images using the Automated Radiative Transfer Models Operator (ARTMO) package [55] (Figure 6).The reflectance values of all fields were combined and used as pooled data to define a global SL in the study area for all fields.As observed from Table 3, the global SL can be defined with high accuracy (R 2 = 0.95) in this study.As mentioned, the bare soil reflectances were collected using a handheld spectrometer in the fields after biomass removal (see Figure 2).These reflectances were spectrally re-sampled according to the Landsat 8 images using the Automated Radiative Transfer Models Operator (ARTMO) package [55] (Figure 6).The reflectance values of all fields were combined and used as pooled data to define a global SL in the study area for all fields.As observed from Table 3, the global SL can be defined with high accuracy (R 2 = 0.95) in this study.As mentioned, the bare soil reflectances were collected using a handheld spectrometer in the fields after biomass removal (see Figure 2).These reflectances were spectrally re-sampled according to the Landsat 8 images using the Automated Radiative Transfer Models Operator (ARTMO) package [55] (Figure 6).The reflectance values of all fields were combined and used as pooled data to define a global SL in the study area for all fields.As observed from Table 3, the global SL can be defined with high accuracy (R 2 = 0.95) in this study.The averaged bare soil reflectances recorded by the handheld spectrometer from different crop fields are also presented in Figure 7 for better visualization.As observed from the figure, a slight concave shape appears between 525 and 627 nm for the barley field soil.Moreover, a slight convex to flat shape was observed between 525 and 693 nm for the canola field soil.Soils with high OM have a convex or flat shape from the visible to NIR wavelengths with low reflectance intensity.However, low OM contents have a higher reflectance intensity.In addition, barley has a higher sand content and low clay content, as revealed by the slight concave shape and value of the table in Figure 7.The action of these components is also discriminated by the angle at which the spectral curve initiates.We observe that the canola spectra have a low angle and that barley spectra display a higher angle.Due to OM, both samples do not indicate the absorption features of hematite, as ratified by its soil color.
Different characteristics and statistical information for the SL parameters extracted from Landsat 8 OLI images using the (Red-NIRmin) and quantile regression using the floating tau and fixed tau (i.e., DIFN) methods are listed in Tables 4 and 5 for the study area and field scale, respectively.The p-values of SLs extracted by the (Red-NIRmin) method are statistically significant except for the intercept on the 10th of August, whereas the quantile regression using floating tau and fixed tau (i.e., DIFN) methods have statistically significant slopes and intercept values for all dates.
The slopes of the SLs extracted by the quantile regression using the floating tau method in the heading stage of winter wheat and the slopes of the SLs extracted by the (Red-NIRmin) method for canola in the heading stage are not statistically significant (Table 5) at the field scale.The averaged bare soil reflectances recorded by the handheld spectrometer from different crop fields are also presented in Figure 7 for better visualization.As observed from the figure, a slight concave shape appears between 525 and 627 nm for the barley field soil.Moreover, a slight convex to flat shape was observed between 525 and 693 nm for the canola field soil.Soils with high OM have a convex or flat shape from the visible to NIR wavelengths with low reflectance intensity.However, low OM contents have a higher reflectance intensity.In addition, barley has a higher sand content and low clay content, as revealed by the slight concave shape and value of the table in Figure 7.The action of these components is also discriminated by the angle at which the spectral curve initiates.We observe that the canola spectra have a low angle and that barley spectra display a higher angle.Due to OM, both samples do not indicate the absorption features of hematite, as ratified by its soil color.

Soils of Different
Different characteristics and statistical information for the SL parameters extracted from Landsat 8 OLI images using the (Red-NIR min ) and quantile regression using the floating tau and fixed tau (i.e., DIFN) methods are listed in Tables 4 and 5 for the study area and field scale, respectively.The p-values of SLs extracted by the (Red-NIR min ) method are statistically significant except for the intercept on the 10th of August, whereas the quantile regression using floating tau and fixed tau (i.e., DIFN) methods have statistically significant slopes and intercept values for all dates.
The slopes of the SLs extracted by the quantile regression using the floating tau method in the heading stage of winter wheat and the slopes of the SLs extracted by the (Red-NIR min ) method for canola in the heading stage are not statistically significant (Table 5) at the field scale.
According to the vegetation growing stages, winter wheat started to green in the middle of April, became mature in July and reached senescence in middle to late August, whereas barley and canola started to green in the middle of April, became mature in June, and reached senescence in middle to late July.The xand y-axes in the scatter plots show the surface reflectance in the red and near-infrared bands, respectively, of Landsat 8 OLI (Figures 8 and 9).The slope of the SLs extracted from the (Red-NIR min ) method is higher than that extracted from the quantile regression methods (Figure 8), except for the 10th of August in the study area, when the canola and barley were already harvested.The slopes of the SLs extracted by the (Red-NIR min ) method are closer to the actual SL for all dates except for the 10th of August, when the quantile regression methods show better agreement with the actual SL.The SLs retrieved from both quantile regression methods have lower slopes and higher intercepts than those extracted from the (Red-NIR min ) method and the actual soil.They do not precisely fit the bottom line except for the last date, which means that sparse points at the bottom of scatter plots cause difficulty in quantile regression for extraction of SL, however, with a relatively straight bottom line on the 10th of August, these methods could estimate the actual SL better than the (Red-NIR min ) method.In other words, the density and/or distribution of the surface reflectance values have an influence on the quantile regression results.Therefore, the sparse points at the bottom of the scatter plot affect the quantile regression results on different dates.
Field Soils OM (%) Ash Content (%) Sand (%) Silt (%) Clay (%)        Based on Figure 9, at the field scale, the scatter plots contain certain sparse points during both the heading and end of the growing stages.These sparse points lead to a certain amount of uncertainty in capturing the bottom boundary line for all of the methods except for the situation in which the scatter plot shows a relatively straight bottom boundary line, such as ripening of wheat Based on Figure 9, at the field scale, the scatter plots contain certain sparse points during both the heading and end of the growing stages.These sparse points lead to a certain amount of uncertainty in capturing the bottom boundary line for all of the methods except for the situation in which the scatter plot shows a relatively straight bottom boundary line, such as ripening of wheat and senescence of barley.As mentioned previously, the overall distribution of data influences SLs extraction via the quantile methods, particularly the scatter plots in the heading stage of the crops (Figure 9a,c,e).The slope of the bottom boundary line extracted from the (Red-NIRmin) method is higher than that extracted from the quantile regression methods (Figure 9a,b and Table 5) for winter wheat.Furthermore, the same phenomenon can be observed for barley, except for the senescence stage.The slopes of the bottom boundary line of canola extracted using the quantile regression methods are higher than that of the (Red-NIRmin method) in the heading stage.The scatter plots in the heading stage of all crops (Figure 9a,c,e) show an accumulated region with relatively higher NIR reflectance in the low red values.
The winter wheat had already started to ripen on the 10th of August, and the kernels were hard on this date.At this stage, the bottom boundary lines extracted using the quantile regression methods form two relatively similar straight bottom lines for winter wheat (Figure 9b), and it should be noted that according to Table 5, although they have the same value for the slope, their intercepts show slightly different values.Additionally, the scatterplot for the senescence period of barley The slope of the bottom boundary line extracted from the (Red-NIR min ) method is higher than that extracted from the quantile regression methods (Figure 9a,b and Table 5) for winter wheat.Furthermore, the same phenomenon can be observed for barley, except for the senescence stage.The slopes of the bottom boundary line of canola extracted using the quantile regression methods are higher than that of the (Red-NIR min method) in the heading stage.The scatter plots in the heading stage of all crops (Figure 9a,c,e) show an accumulated region with relatively higher NIR reflectance in the low red values.
The winter wheat had already started to ripen on the 10th of August, and the kernels were hard on this date.At this stage, the bottom boundary lines extracted using the quantile regression methods form two relatively similar straight bottom lines for winter wheat (Figure 9b), and it should be noted that according to Table 5, although they have the same value for the slope, their intercepts show slightly different values.Additionally, the scatterplot for the senescence period of barley growth reveals the ability of the quantile regression to capture the lower bottom line.Moreover, the scatter plot has a relatively clear straight bottom line the same as the ripening period of winter wheat.For the scatter plots in the ripening and senescence stages of winter wheat and barley with a clear straight bottom line, all three lines simulate the bottom lines quite well.In the red near-infrared spectral space, the surface reflectances in the heading stages and at the end of the growing stages of the crops have scatter plots with higher NIR reflectance values compared with the actual SL (Figure 9).Therefore, the SL cannot be observed in any of the aforementioned subfigures except for the ripening of winter wheat (Figure 9b).
To validate the soil lines from the ground truth data collection and the SL simulated with different methods, a cloud-free Landsat 8 image of the barley field was acquired after harvest.This image shows the reflectance of the bare soil of the barley field.The field photo of the corresponding image is shown in Figure 2d.
As observed in Figure 10a, the slope and intercept of the scatter plot are rather close to those measured on bare soil using a handheld spectrometer for both the barley field and global SL.Furthermore, this slope and intercept are also close to those from the simulated SL using quantile regression methods (fixed and floating tau) for the study area, especially at the end of growing season (i.e., Table 4, 2013/08/10).Indeed, these data are in agreement with [56,57], for whom soils with different textures appear at different positions in the SL extracted from Landsat TM (Figure 10b,c).The pixels in Figure 10 are taken from bare soil, and it is interesting to note that they are spread all over the soil line, which indicates that the pixels have different textures.The pixels at the top are sandier than the ones at the bottom.According to Demattê et al., [56] the pattern (i.e., Figure 10b) values on the lower left of the scatter plot (i.e., black circles) represent clay soil.The results of granulometric analysis for barley soil in Figure 7 (bottom table) confirm that barley soil has a notably low percentage of clay (approximately 2%).The results of Landsat 8 OLI for the barley field shown in Figure 10a and the results of field spectrometry shown in Figure 6 confirms this fact, with fewer points on the lower left of the scatter plot.However, Figure 10a also shows that most of the points are located in the middle, which represents loamy soil, according to the Demattê et al. [56] pattern and Figure 10c, but our granulometric analysis shows that the amount of silt is approximately 17% and the barley soil is mostly sandy.This result might be because our granulometric analysis was based on a few samples that were taken from a small area within the fields or could be attributed to a direct response to the variability in soil moisture.As observed from Figure 10d [58], differences in soil moisture cause a shift of the values along the SL with moister values in the lower portion and drier values in the higher red and NIR regions.The Demattê et al., [56] pattern in Figure 10b also shows the scattered sandy soil points along the SL due to spread albedos.The differences in albedo can be attributed again to differences in soil moisture.Therefore, the quantile regression methods were successful in simulating the SL with high accuracy, especially at the end of the growing season.Unfortunately, due to cloud cover and rapid re-cultivation, no bare soil images were available for the winter wheat and canola field.
Figure 11 shows the correlation coefficients and RMSE values for the dry biomass of winter wheat, barley, and canola with the indices based on the extracted SLs from this study and the slope and intercept from the global SL using the pooled reflectance data for the study area.All of the coefficients are at a significance level of p-value less than 0.05.In the local area, for the winter wheat in this study, based on the SLs of the 6th of May from all of the methods except quantile regression using DIFN as the fixed tau, TSAVI slightly improved the dry biomass extraction over NDVI and other VIs (Figure 11a).The same phenomenon was observed for PVI2 on the 9th of July using all of the methods by which this index shows its superiority against NDVI and the other VIs.The NDVI performs better compared with the other indices on the 7th of June and 10th of August for winter wheat.Estimation of dry biomass in the barley fields shows that GESAVI is a better estimator of dry biomass than NDVI and the other VIs for the 06th of May.PVI2 can estimate the dry biomass of barley better than NDVI and the other indices on the 7th of June and 9th of July.For canola, ATSAVI and PVI2 show stronger relationships than NDVI and the other SL-related indices based on all methods on the 6th of May and 7th of June, respectively.However, NDVI shows its superiority against the other SL VIs on the 9th of July for the assessment of dry biomass of canola.
Figure 12 presents the correlation coefficients and RMSE values for dry biomass of the three crops with the indices based on the extracted SLs from this study and the slope and intercept of the actual SL at the field scale.It is worth mentioning that different indices produce the highest coefficient correlation with different methods for different dates and growing stages, but most of the indices show the same relative performance throughout the entire growing stages of the crops except for TSAVI, which shows relatively different R 2 values, especially during the heading stage of all three crops (Figure 12).A decreasing trend in the correlation coefficients is observed for all of the studied indices from the heading stage to the end of the life cycle for winter wheat and barley, whereas R 2 increases during the ripening period of canola.Figure 12 presents the correlation coefficients and RMSE values for dry biomass of the three crops with the indices based on the extracted SLs from this study and the slope and intercept of the actual SL at the field scale.It is worth mentioning that different indices produce the highest coefficient correlation with different methods for different dates and growing stages, but most of the indices show the same relative performance throughout the entire growing stages of the crops except for TSAVI, which shows relatively different R2 values, especially during the heading stage of all three crops (Figure 12).A decreasing trend in the correlation coefficients is observed for all of the studied indices from the heading stage to the end of the life cycle for winter wheat and barley, whereas R 2 increases during the ripening period of canola.

Discussion
As mentioned previously, the small amount of organic matters (OM) in the field (i.e., high ash content) throughout the growing season leads to the small influence of this factor on the soil surface reflectance [6][7][8][9]57].Therefore, the slightly higher slope and intercept of the SL extracted from the barley field can be attributed to a lower OM because the soil reflectance increases as this attribute decreases.This result occurs because according to the soil texture measurements, both fields have a relatively similar texture, as corroborated by the barley soil spectra in Figure 7, which has a slightly concave shape, indicating low OM, in accordance with [59].The higher slope of the barley field might be related to the higher amount of silt in the soil texture, causing increasing reflectance with decreasing particle size, which also agrees with the findings of Demattê et al. (2014) [59].Previous studies, such as [60], have suggested that for a given spectral measuring instrument (a spectrometer or imager), an unique SL (global SL) appears to exist for soil reflectance measurements.This observation is confirmed by the results of Figure 6 and Table 3, which show that using the pooled reflectance data, a SL can be defined for the entire study area with high accuracy (high R 2 ), and an SL with approximately similar slopes and intercepts can be found [61].
Pixels that contain the pure bare-soil spectral signature were not available in this study at the field scale.Thus, it was much more difficult to determine the SL parameters as mentioned by [8].According to Equation (2) [52]: The slope (α) of the vegetation isoline in the Red-NIR space is dependent on the slope of SL " ," LAI, and difference in the red and NIR canopy extinction coefficient.The red and NIR transmittance (extinction) through a photosynthetically active canopy differ significantly with a much higher optical thickness in the red band due to the highly absorptive properties of leaf pigments and relatively low NIR optical extinction due to the highly scattered signal [62].If LAI = 0 (i.e., bare soil), (α) is the bare soil slope.In the same manner, if kred, is notably close to knir or kred = knir, the slope of vegetation isolines should be similar to the SL slope "a" for a relatively constant LAI value.This is the particular case for senescent vegetation and corroborates the results for which the "senescent line" [16,51] has the same slope as the corresponding bare SL.This fact can be confirmed by Figure 8d and Figure 9b (ripening stage of winter wheat), when all three methods were able to retrieve the hypothetical bottom line without any challenges, and the slope and intercepts of the extracted SL are relatively close to the slopes and intercepts of actual SLs.If kred > knir the slope of the vegetation

Discussion
As mentioned previously, the small amount of organic matters (OM) in the field (i.e., high ash content) throughout the growing season leads to the small influence of this factor on the soil surface reflectance [6][7][8][9]57].Therefore, the slightly higher slope and intercept of the SL extracted from the barley field can be attributed to a lower OM because the soil reflectance increases as this attribute decreases.This result occurs because according to the soil texture measurements, both fields have a relatively similar texture, as corroborated by the barley soil spectra in Figure 7, which has a slightly concave shape, indicating low OM, in accordance with [59].The higher slope of the barley field might be related to the higher amount of silt in the soil texture, causing increasing reflectance with decreasing particle size, which also agrees with the findings of Demattê et al. (2014) [59].Previous studies, such as [60], have suggested that for a given spectral measuring instrument (a spectrometer or imager), an unique SL (global SL) appears to exist for soil reflectance measurements.This observation is confirmed by the results of Figure 6 and Table 3, which show that using the pooled reflectance data, a SL can be defined for the entire study area with high accuracy (high R 2 ), and an SL with approximately similar slopes and intercepts can be found [61].
Pixels that contain the pure bare-soil spectral signature were not available in this study at the field scale.Thus, it was much more difficult to determine the SL parameters as mentioned by [8].
According to Equation (2) [52]: The slope (α) of the vegetation isoline in the Red-NIR space is dependent on the slope of SL "a," LAI, and difference in the red and NIR canopy extinction coefficient.The red and NIR transmittance (extinction) through a photosynthetically active canopy differ significantly with a much higher optical thickness in the red band due to the highly absorptive properties of leaf pigments and relatively low NIR optical extinction due to the highly scattered signal [62].If LAI = 0 (i.e., bare soil), (α) is the bare soil slope.In the same manner, if k red , is notably close to k nir or k red = k nir , the slope of vegetation isolines should be similar to the SL slope "a" for a relatively constant LAI value.This is the particular case for senescent vegetation and corroborates the results for which the "senescent line" [16,51] has the same slope as the corresponding bare SL.This fact can be confirmed by Figures 8d and 9b (ripening stage of winter wheat), when all three methods were able to retrieve the hypothetical bottom line without any challenges, and the slope and intercepts of the extracted SL are relatively close to the slopes and intercepts of actual SLs.If k red > k nir the slope of the vegetation isoline becomes greater than the SL slope [51,52].Thus, the difference between red and NIR extinction through a vegetated canopy and the LAI value [63] determines the slope behavior of a vegetation isoline and the subsequent soil influence.This phenomenon can also explain the difficulties in extracting the bottom SL at the field scale during the heading stages of crops using all three methods.Furthermore, the clusters with little variation of the surface reflectance ranges along both the red and NIR bands during the heading stage of the crops create a challenge in simulating the bottom boundary line parameters.This phenomenon was also observed by Fox et al. [8].Additionally, the relatively smaller influence of the soil on the vegetation isolines at larger LAI (lower DIFN) values due to the lower transmittance [64] can explain the difficulties in extracting the bottom lines and the significant difference in the bottom boundary lines with the actual SL, especially at the field scale (Figure 9).In the other words, it can be assumed that the spectral reflectance of the crop canopy is a mixture of the reflectance spectra of the crop and the soil beneath it.As the vegetation grows, the contribution of the soil considerably decreases (decreasing DIFN values in Figure 5), but might still remain significant, depending on many factors, such as the crop density, row distance, canopy geometry, wind effects [6].Even at the end of the growing season, when the crop is likely to be completely senescent, the bottom boundary line has a higher slope and intercept and still lies above the actual SL.Therefore, these observations suggest that the calculated bottom boundary line is actually a "senescent vegetation line" (Figure 9b,d).
Thompson and Wehmanen [65] proposed a technique for extracting the SL by calculating the Kauth-Thomas vectors, rejecting any improper vector for agricultural data, assuming the remaining pixels as good pixels, and discarding 1% of the pixels with the lowest greenness values.These researchers also hypothesized that the remaining lowest greenness value becomes the SL.Rejection of a further 1% of the total data creates further protection against low outliers [65].Nonetheless, rejection of 1% of the total data lacks a rigid scientific background compared with quantile regression using DIFN as a fixed tau.Hence, the latter can be used as a substitute approach for similar applications.
As discussed by Galvão and Vitorello [66], SLs obtained with broad or narrow NIR bands were positioned at shorter wavelengths; consequently, at smaller distance, differences between the pair of bands (i.e., center difference), presented better fitted lines than those obtained with NIR bands located at longer wavelengths and at larger distances.These researchers concluded that by the same reasoning, the LANDSAT 5-TM might present better results than MSS because the NIR band from the latter is positioned at longer wavelengths.However, compared with the broad bands, spectrally better-positioned narrower NIR/R bands (e.g., Landsat 8 OLI) can also produce similar results.
The study conducted by Galvão and Vitorello [66] showed that alterations in the slope and intercept of the SLs indicate that the reflectance segments that compose them tend to have less divergence and scattering if the NIR band position is shifted to shorter wavelengths, which could be the reason why the SL simulated by Xu and Guo [12] with different methods and Landsat TM showed slope differences on the order of approximately 0.1 or less, but we generally observed slope differences on the order of approximate 0.3 in this study.
Finally, the narrower-band versions of the vegetation indices (e.g., Landsat 8 OLI) produce similar or only slightly better accuracy than their broad-band counterparts (e.g., Landsat ETM+ and TM) [67,68].Because the SL parameters depend on soil characteristics and band configuration, it is important to understand the possible performance of VIs against changes in the SL parameters under an identical vegetation canopy [63].As mentioned previously, as the amount of vegetation increases, the difference between the red and NIR penetration to the soil surface increases, causing the SL-related VIs lose their ability to estimate crop parameters properly [52].This fact can explain the lower performance of VIs during the heading of winter wheat and barley crops.As shown by [16], for a LAI greater than 4 (heading stages in this study), the worst index is TSAVI.This phenomenon can be observed in Figure 12 for winter wheat and barley.This result might be due to the fact that this VI is near its saturation level and perhaps reaches this level before the other VIs.However, this index shows a notably good performance during the heading stage of canola, and the reason for this outcome still remains unknown.The lower correlation coefficient and higher RMSE values at the end of the growing stage of winter wheat and barley crops can be attributed to the significantly higher standard deviation of the data (Figure 5).
The TSAVI index was defined for a low LAI or low vegetation coverage according to [51]. Figure 11 shows this phenomenon, and this index performs quite well for the lower LAI values (10 August) at the end of the growing season.However, this index shows worse performance for higher LAI values, especially during the heading stage of winter wheat and barley crops.
Although the value of the Z-coefficient for GESAVI was selected to minimize the variation of the canopy reflectance with the soil background for intermediate LAI levels [32], this value was nevertheless well suited for the wide range of canopy LAI [32] values in this study.A fixed value of L = 0.5 was used in this study because the determination of this factor requires information on the vegetation isoline at LAI = 1 prior to the optimization process, which is most likely not available in real applications, as mentioned by [64].It appears that the good results of SAVI are in fact due to the almost constant sensitivity to soil at all LAI values [6].Furthermore, Baret and Guyot [16] mentioned that ATSAVI and SAVI are closely correlated, in accordance with Figures 11 and 12.
The optimal choice of a vegetation index is indeed related to the purpose of the study and the type of vegetation considered [6].According to the performances of different indices in this study over a single soil type and single field, small differences are noted between the indices studied based on different crop types throughout the crop life cycle.
In evaluating the SL from RS data, such as those from Landsat 8, use of a larger scale, such as local and/or regional data (rather than data solely from within an agricultural field), appear to produce better results because the presence of pure bare soil pixels is more likely in regional image data.Furthermore, although these three methods can retrieve the bottom boundary line of the red-NIR spectral space on the field scale, due to the penetration ability of optical wavelengths, these bottom boundary lines are significantly different from the actual SL.Therefore, if retrieving the SL when the ground is highly vegetated and if no information is available on the soil, use of these approaches could be critical at the field scale.

Conclusions
The ability to extract the SL parameters without manual investigation of bare soil pixels within the images allows the SL-related vegetation indices or other applications to be more easily investigated.Simulation of the SL parameters in the presence of vegetation using remote sensing data is a critical issue because a major shift in the reflectance values occurs for red and near-infrared radiation impinging on a crop canopy during the season.However, the actual SL does not change dramatically in croplands during one growing season because the soil attributes usually remain the same during one growing season (if the soil has not been altered by tillage).
The capability of the three statistical methods in SL simulation varies in different crop growth stages because the bottom lines change through crop growth, and it is difficult to say whether (Red-NIR min ) or quantile regression is the better approach.We found that the slopes of the SLs extracted by the (Red-NIR min ) method are closer to the actual SL during the crop life cycle.However, the quantile regression methods show better agreement with the actual SL at the end of growing season and if more patches of soil are available in the image (certain fields are already harvested) on the local scale.The better agreement of the quantile methods at the end of growing stage was also confirmed by the Landsat 8 image acquired on the bare soil of the barley field.We also found that the spectral reflectance of a crop canopy can be assumed to be a mixture of the reflectance spectra of the crop and soil beneath it and the contribution of soil considerably decreases as the vegetation grows, but might still remain significant (depending on the DIFN value).Furthermore, extraction of the hypothetical SL is highly critical at the field scale because even at the end of the growing season, the bottom boundary line has a higher slope and intercept and still lies above the actual SL.
Over a single soil type and field, the optimal choice of a proper VI is difficult.This conclusion can be stated because a small difference is observed among VIs based on their abilities to retrieve the biophysical parameters of crops during a special phenological stage using the approaches studied in this research.
The proposed method (quantile regression using DIFN as the fixed tau) is not intended as a replacement for the conventional methods, but rather as a complementary tool for extraction, analysis and interpretation of SL concepts in future studies.

Figure 1 .
Figure 1.Study area including barley (A), winter wheat (B), and canola (C) at the (DEMMIN) test site.As an example of ground truth data collection, the light-orange squares show the location of the soil samples, and the blue squares show the quadrats.The red circle represents the center.The background data were taken from GeoBasis (M-V) (DOP40), GeoBasis (M-V) (ATKIS), and GeoBasis (DE).

Figure 1 .
Figure 1.Study area including barley (A), winter wheat (B), and canola (C) at the (DEMMIN) test site.As an example of ground truth data collection, the light-orange squares show the location of the soil samples, and the blue squares show the quadrats.The red circle represents the center.The background data were taken from GeoBasis (M-V) (DOP40), GeoBasis (M-V) (ATKIS), and GeoBasis (DE).

Figure 2 .
Figure 2. Field data collection photos and flowchart of methods used in this study.Photos of the fields during different growing season: (a) Beginning of Growing Season; (b) Middle of Growing Season; (c) End of Growing Season; (d) After Harvest; (e) Soil sample and soil reflectance measurement; (f) Flowchart of methods used to retrieve the SL parameters and estimating biomass in this study.

Figure 2 .
Figure 2. Field data collection photos and flowchart of methods used in this study.Photos of the fields during different growing season: (a) Beginning of Growing Season; (b) Middle of Growing Season; (c) End of Growing Season; (d) After Harvest; (e) Soil sample and soil reflectance measurement; (f) Flowchart of methods used to retrieve the SL parameters and estimating biomass in this study.

Figure 3 .
Figure 3. Schematic figure of quantile regression using the fixed tau (DIFN) approach in this study.The background photo was taken from: [49,50].

Figure 3 .
Figure 3. Schematic figure of quantile regression using the fixed tau (DIFN) approach in this study.The background photo was taken from: [49,50].

Table 2 .
SL-related vegetation indices based on soil line (SL) parameters.Vegetation Index Algorithm Description PVI [15] N IR − aRED − b √ a 2 + 1 a = slope of the SL, b = intercept of the SL TSAVI 2 or ATSAVI [16] a(N IR − aRED − b) RED + aN IR − ab + X(1 + a 2 ) a = slope of the SL, b = intercept of the SL, X = 0.08 GESAVI [32] N IR − bRED − a RED + Z Z = 0.35, b = slope of the SL, a = intercept of the SL TSAVI [51] a(N IR − aRED − b) RED + aN IR − ab a = slope of the SL, b = intercept of the SL SAVI [52] Remote Sens. 2016, 8, 738 9 of 23 during the entire growing season.The error bars of the dry biomass depict the standard deviation of the data.

Figure 4 .
Figure 4. Analysis of the soil ash content (average value) for different crop fields derived from soil samples during the entire growing season.

Figure 5 .
Figure 5. In-situ measurements of dry biomass and the DIFN value for three crop fields during the entire growing season.

Figure 4 .
Figure 4. Analysis of the soil ash content (average value) for different crop fields derived from soil samples during the entire growing season.

Figure 4 .
Figure 4. Analysis of the soil ash content (average value) for different crop fields derived from soil samples during the entire growing season.

Figure 5 .
Figure 5. In-situ measurements of dry biomass and the DIFN value for three crop fields during the entire growing season.

Figure 5 .
Figure 5. In-situ measurements of dry biomass and the DIFN value for three crop fields during the entire growing season.

Figure 6 .
Figure 6.Resampling of measured bare soil reflectance collected from the fields after biomass removal using a handheld spectrometer based on the red and NIR bands of Landsat 8 OLI.The plus sign (+), circle (•), and upward-pointing triangle represent the bare soil reflectances of the canola, barley and wheat fields, respectively.Table 3. Slopes, intercepts and statistical information for SLs extracted from Figure 6 for different soils of crop fields and pooled data for the study area from field data collection.a and b represent the slope and intercept of the SL, respectively.Soils of Different Fields a b R 2 RMSE Wheat 1.187 0.02862 0.9821 0.0089 Barley 1.188 0.03929 0.9632 0.0125 Canola 1.362 0.02875 0.9636 0.0134 Global SL (Pool Data) 1.22 0.03425 0.9555 0.0147

Figure 6 .
Figure 6.Resampling of measured bare soil reflectance collected from the fields after biomass removal using a handheld spectrometer based on the red and NIR bands of Landsat 8 OLI.The plus sign (+), circle ( ), and upward-pointing triangle represent the bare soil reflectances of the canola, barley and wheat fields, respectively.

Figure 7 .
Figure 7.Comparison of the average bare soil reflectances recorded on different crop fields after biomass removal using a hand-held spectrometer, texture analysis, amount of OM, and ash content in their soils.

Figure 7 .
Figure 7.Comparison of the average bare soil reflectances recorded on different crop fields after biomass removal using a hand-held spectrometer, texture analysis, amount of OM, and ash content in their soils.

Figure 8 .
Figure 8. SLs extracted by the (Red-NIRmin) method, quantile regression method using floating tau, and quantile regression method using fixed tau equal to DIFN from Landsat 8 OLI images, and the actual SL for the study area retrieved from resampling of the spectrometer on different dates: (a) 6 May; (b) 7 June; (c) 9 July; and (d) 10 August 2013.

Figure 8 .
Figure 8. SLs extracted by the (Red-NIRmin) method, quantile regression method using floating tau, and quantile regression method using fixed tau equal to DIFN from Landsat 8 OLI images, and the actual SL for the study area retrieved from resampling of the spectrometer on different dates: (a) 6 May; (b) 7 June; (c) 9 July; and (d) 10 August 2013.
Remote Sens. 2016, 8, 738 13 of 23 and senescence of barley.As mentioned previously, the overall distribution of data influences SLs extraction via the quantile methods, particularly the scatter plots in the heading stage of the crops (Figure 9a,c,e).

Figure 9 .
Figure 9. SLs extracted by the (Red-NIRmin) method, quantile regression method using floating tau, and quantile regression method using fixed tau equal to DIFN from Landsat 8 OLI images, and the actual SL retrieved from resampling of spectrometer for (a) heading of wheat; (b) ripening of wheat (c) heading of barley; (d) senescence of barley; and (e) heading of canola (f) ripening of canola.

Figure 9 .
Figure 9. SLs extracted by the (Red-NIR min ) method, quantile regression method using floating tau, and quantile regression method using fixed tau equal to DIFN from Landsat 8 OLI images, and the actual SL retrieved from resampling of spectrometer for (a) heading of wheat; (b) ripening of wheat (c) heading of barley; (d) senescence of barley; and (e) heading of canola (f) ripening of canola.

Figure 10 .
Figure 10.(a) Scatter plot of the bare soil reflectance for the barley field soil retrieved from the Landsat 8 OLI image after harvest; (b) scatter plot developed by [56] indicating the SL of soils with different textures; (c) SL pixels from Landsat TM [57], where the upper circle represents sandy soils.The bottom represents clay and the middle indicates loam; (d) Soil line evaluated from multi-temporal image data for georeferenced bare soil surfaces with varying degrees of moisture [58].

Figure 10 .Figure 10 .Figure 11 .
Figure 10.(a) Scatter plot of the bare soil reflectance for the barley field soil retrieved from the Landsat 8 OLI image after harvest; (b) scatter plot developed by [56] indicating the SL of soils with different textures; (c) SL pixels from Landsat TM [57], where the upper circle represents sandy soils.The bottom represents clay and the middle indicates loam; (d) Soil line evaluated from multi-temporal image data for georeferenced bare soil surfaces with varying degrees of moisture [58].

Figure 11 .
Figure 11.Statistical information of the indices based on the slope and intercepts of SLs extracted by different methods and global SL for (a) dry biomass of winter wheat; (b) dry biomass of barley; (c) dry biomass of canola for the study area.

Figure 12 .
Figure 12.Statistical information of the indices for estimation of the dry biomass of crops based on the slope and intercepts of SLs extracted by different methods and the actual SL of each crop at heading and end of growing stage for different crops.

Table 1 .
Variable names, number of samples for ground truth data collection and satellite images collected during the entire crop-growing season.

Table 4 .
Slope and intercept of extracted SLs from different methods using Landsat 8 images for the study area.Date format: yy/mm/dd.

Table 5 .
Slope and intercept of SLs extracted by different methods for different crops at the field scale using Landsat 8 images.Date format: yy/mm/dd.

Table 4 .
Slope and intercept of extracted SLs from different methods using Landsat 8 images for the study area.Date format: yyyy/mm/dd.

Table 5 .
Slope and intercept of SLs extracted by different methods for different crops at the field scale using Landsat 8 images.Date format: yyyy/mm/dd.

Table 5 .
Slope and intercept of SLs extracted by different methods for different crops at the field scale using Landsat 8 images.Date format: yy/mm/dd.
Statistical information of the indices for estimation of the dry biomass of crops based on the slope and intercepts of SLs extracted by different methods and the actual SL of each crop at heading and end of growing stage for different crops.