Next Article in Journal
Use of Radarsat-2 and Landsat TM Images for Spatial Parameterization of Manning’s Roughness Coefficient in Hydraulic Modeling
Previous Article in Journal
Modeling Aboveground Biomass in Dense Tropical Submontane Rainforest Using Airborne Laser Scanner Data
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Developing in situ Non-Destructive Estimates of Crop Biomass to Address Issues of Scale in Remote Sensing

1
Climate Research Group, World Agroforestry Centre, United Nations Ave, Gigiri, P.O. Box 30677-00100, Nairobi, Kenya
2
Southwestern Geographic Center, United States Geological Survey, 2255 N. Gemini Dr, Flagstaff, AZ 86001, USA
*
Author to whom correspondence should be addressed.
Remote Sens. 2015, 7(1), 808-835; https://0-doi-org.brum.beds.ac.uk/10.3390/rs70100808
Submission received: 11 November 2014 / Accepted: 7 January 2015 / Published: 14 January 2015

Abstract

:
Ground-based estimates of aboveground wet (fresh) biomass (AWB) are an important input for crop growth models. In this study, we developed empirical equations of AWB for rice, maize, cotton, and alfalfa, by combining several in situ non-spectral and spectral predictors. The non-spectral predictors included: crop height (H), fraction of absorbed photosynthetically active radiation (FAPAR), leaf area index (LAI), and fraction of vegetation cover (FVC). The spectral predictors included 196 hyperspectral narrowbands (HNBs) from 350 to 2500 nm. The models for rice, maize, cotton, and alfalfa included H and HNBs in the near infrared (NIR); H, FAPAR, and HNBs in the NIR; H and HNBs in the visible and NIR; and FVC and HNBs in the visible; respectively. In each case, the non-spectral predictors were the most important, while the HNBs explained additional and statistically significant predictors, but with lower variance. The final models selected for validation yielded an R2 of 0.84, 0.59, 0.91, and 0.86 for rice, maize, cotton, and alfalfa, which when compared to models using HNBs alone from a previous study using the same spectral data, explained an additional 12%, 29%, 14%, and 6% in AWB variance. These integrated models will be used in an up-coming study to extrapolate AWB over 60 × 60 m transects to evaluate spaceborne multispectral broad bands and hyperspectral narrowbands.

Graphical Abstract

1. Introduction

The quantification of the terrestrial carbon balance is important for understanding anthropogenic climate change and costing ecosystem services, but remains a challenge, due to uncertainties in simulations of carbon sources and sinks over space and time [1]. Agriculture accounts for approximately 25% of the global greenhouse gas budget [2]. The two main sources of agriculture-related emissions—nitrous oxide and methane—come mostly from fertilizer/manure use and production, livestock ruminants and manure management, bush fires, and rice cultivation, while carbon dioxide (CO2) emissions are second to emissions from fossil fuel consumption and are primarily due to land use/cover change [3]. Given its importance in the global greenhouse gas budget and climate forcing, methods for estimating CO2 emissions from agriculture are important for designing appropriate mitigation strategies [4]. Plant biomass, defined here as the total aboveground dry-weight biologically-produced matter, is an important surrogate and relatively easy measure of carbon assimilation [5].
Physiological-based crop growth models are used to monitor and forecast crop condition and production to inform farm-level management and local, national, and global policy-making [6]. These models typically use the crop or variety-specific fraction of harvestable crop biomass (harvest index: [7]) to estimate yield [8]. Airborne/spaceborne remote sensing technology has been used since the early 1980s to estimate crop biomass and other biophysical components of crop growth models, because it offers a near real-time snapshot or long-term record of crop condition over multiple spatial scales [9]. The accuracy of remote sensing methods varies and depends largely on the ability of modelers to collect and scale up in situ data for model calibration/validation [10]. In order for remote sensing methods to accurately capture biomass, therefore, rapid and cost-effective ground-based techniques are needed that account for spatial heterogeneity.
Ground-based estimates of crop biomass can be made using either in situ destructive, non-spectral or spectral methods [11]. Due to the difficulty of collecting below-ground biomass, biomass estimation tends to focus on the above-ground portion only. In situ destructive methods are time-consuming and costly, because they involve harvesting, pruning, and weighing to estimate aboveground wet (fresh) biomass (AWB) and further drying to estimate crop biomass. In addition, this technique may be prohibited by land owners or managers, because it is destructive. In situ non-spectral and spectral techniques, on the other hand, facilitate intensive and extensive sampling campaigns with minimal financial impact on growers. In addition, with a high density of biomass estimates generated by these techniques, transects can be made that account for uncertainties over heterogeneous canopies for remote sensing applications [12]. Reviews of in situ non-spectral methods pertaining to crops and rangelands can be found in [13,14,15]. These techniques involve: visual assessment (e.g., [16]), crop height (H) (e.g., [17]), fraction of absorbed photosynthetically active radiation (FAPAR) derived from a ceptometer, fraction of vegetation cover (FVC) derived from red-green-blue (RGB) band photographs (e.g., [18]), electronic capacitance probes (e.g., [19]), weighted discs (e.g., [20]) or pendulum sensors (e.g., [21]).
Spectral techniques rely on the unique spectral absorption and scattering characteristics of crops and varieties. Ollinger [22] reviews these relationships and concludes based on previous studies that the visual and near infrared (NIR) regions of the spectrum are the most important for simulating biophysical and biochemical plant processes. In the visible range, preferential absorption occurs in the blue around 450 nm, scatter in the green around 550 nm, and absorption of red around 630–690 nm, due to chlorophyll and accessory pigments in plant leaves. In the NIR (700–1000 nm), plants scatter largely from the lack of absorption by chlorophyll and accessory pigments and alignment of plant cell walls [23]. The spectral region where a rapid change in absorption in the visible red and scatter in the NIR occurs (“red-edge”: 700–740 nm) is effective at detecting differences in crop stage [24] and leaf nitrogen content [25]. Ustin et al. and Goetz [26,27] review these relationships and extend the discussion to the Short-Wave Infrared 1 (SWIR1: 1000–1700 nm) and Short-Wave Infrared 2 (SWIR2: 1700–2500 nm) in the context of hyperspectral narrow bands (HNBs) retrieved from airborne/spaceborne sensors. Light around 1450, 1940, and 2500 nm is especially sensitive to leaf water content [28], while light around 2000 and 2200 nm is sensitive to structural carbon (lignin-cellulose) [29].
The primary goal of this manuscript is to convey a non-destructive method for calibrating/validating airborne/spaceborne remote sensing crop biomass estimates that integrates in situ non-spectral and spectral methods. The data were collected over several 60 × 60 m transects to identify the optimal blend of high, medium, and coarse broadband and hyper-spectral spaceborne sensors for estimating crop biomass over the Central Valley of California in an upcoming study. It is anticipated that the method can be used to calibrate/validate crop biomass models derived from upcoming missions, such as the Hyperspectral Infrared Imager (HyspIRI: [30]) as well. The dataset was calibrated empirically using in situ destructive measurements of AWB. We used AWB, because it is simple to implement in the field and highly correlated with crop biomass used in productivity modeling [31,32]. The dataset consisted of California’s leading water-intensive crops (alfalfa, cotton, maize, and rice). The in situ non-destructive techniques involved measuring or collecting: (a) H; (b) FAPAR; (c) FVC; and (d) canopy HNBs from 350 to 2500 nm at 1–10 nm intervals, depending on the spectral region. Given the large number of HNBs, data mining techniques discussed in [33,34] were used to develop indices of AWB, thus eliminating redundant bands and utilizing the non-redundant optimal wavebands to extract the best possible information.

2. Material and Methods

2.1. Study Area

In the summer of 2011 and 2012, spectral data, canopy structure, and AWB were collected during the major phenological phases (sprouting, flowering/silking, and bud/grain-filling) of alfalfa, cotton, maize, and rice varieties in the Central Valley of California. California is the most agriculturally productive state in the United States (US), representing approximately 12% of US agriculture revenues in 2012 [35]. The Central Valley of California spans a distance north to south of more than 700 km and covers an area of approximately 60,000 km2 between the Coastal and Sierra Nevada ranges of California (Figure 1). It includes seven of the most productive counties of the state (Fresno, Tulare, Kern, Merced, Stanislaus, San Joaquin, and King). Agricultural irrigation accounts for 75%–80% of the state’s annual water budget [36]. Climate variability and change, land use change, and other conspiring factors are expected to put further strain on water resources [37]. Alfalfa, cotton, maize, and rice are California’s largest water users by crop type [38] yet the tradeoffs between crop productivity and water conservation for these crops remains uncertain, necessitating improved modelling inputs and techniques. In this study, more than 500 AWB samples were collected evenly across crop type and phenological phase on 18 privately owned and University of California, Davis research farms to (1) develop accurate non-destructive estimates of AWB and (2) create 60 × 60 m transects for airborne/spaceborne remote sensing model calibration/validation at Landsat resolution. The transects where measurements were taken are shown in Figure 1 and are stratified across three of the four sub-regional basins of the Central Valley: the Sacramento Basin, the Delta and Eastside Streams, and the Tulare Basin [39].
Figure 1. The spatial extent of non-spectroradiometric, spectroradiometric, and aboveground wet biomass (AWB) samples taken in the Central Valley of California in 2011 and 2012. More than 500+ samples were taken for alfalfa, cotton, maize, and rice. The 60 × 60 m transects where AWB and non-destructive predictors were taken (▪) are overlaid with the National Agricultural Statistics Service Cropland Data Layer for 2012 [40].
Figure 1. The spatial extent of non-spectroradiometric, spectroradiometric, and aboveground wet biomass (AWB) samples taken in the Central Valley of California in 2011 and 2012. More than 500+ samples were taken for alfalfa, cotton, maize, and rice. The 60 × 60 m transects where AWB and non-destructive predictors were taken (▪) are overlaid with the National Agricultural Statistics Service Cropland Data Layer for 2012 [40].
Remotesensing 07 00808 g001

2.2. Non-Spectral Measurements and Processing

Non-spectral measurements collected for each crop included: AWB (g∙m−2), visible canopy and background RGB light intensity, H (cm), and the above and below canopy irradiance. The fraction of vegetation cover was calculated using RGB intensity, while above and below canopy irradiance was used to estimate the canopy gap fraction (β) expressed as a decimal percent and Leaf Area Index (LAI in m2∙m−2). Two to four replicate samples were taken for each measurement over a square-meter quadrat and averaged, so that the empirical AWB models driven by these predictors yielded square-meter units. The location of each quadrat was stored in a Garmin Geographic Information System® (precision = 1 m), so that the measurements could be repeated at each major phenological phase during the growing season. Ten quadrats were taken per 60 × 60 m transect following [41]. Approximately 10 transects were randomly distributed on each farm. The size of the transects were designed to account for geolocation error for the anticipated Landsat (30 m resolution) study. In order to minimize damage to the canopy and prevent the degradation of spectral retrievals, (1) AWB was only taken at one quadrat per transect and visit and (2) AWB samples were taken after all other measurements were taken per visit. The mean and standard deviation of each non-spectral plant measurement at each phenological phase for those quadrats where AWB was measured are shown in Table 1.
Table 1. Summary of statistics of alfalfa, cotton, maize, and rice AWB and non-spectral samples taken for each major phenological crop stage during the 2011 and 2012 growing seasons. Where N is the number of samples, H is height (cm), LAI is the leaf area index (m2∙m−2), β is the gap fraction (%), EXG is the excess greenness index, EXGR is the excess green minus excess red index, NDI is the normalized difference index, meanG is the mean chromatic greenness, and medianG is the median chromatic greenness.
Table 1. Summary of statistics of alfalfa, cotton, maize, and rice AWB and non-spectral samples taken for each major phenological crop stage during the 2011 and 2012 growing seasons. Where N is the number of samples, H is height (cm), LAI is the leaf area index (m2∙m−2), β is the gap fraction (%), EXG is the excess greenness index, EXGR is the excess green minus excess red index, NDI is the normalized difference index, meanG is the mean chromatic greenness, and medianG is the median chromatic greenness.
ParameterCrop StageStatisticAlfalfa (N = 136)Cotton (N = 147)Maize (N = 151)Rice (N = 106)
AWBSproutingμ19847767558774
(g∙m−2) σ25177573882589
Floweringμ161328134139592703
σ12945402528101089
Senescenceμ-9447118572727
σ-577034601133
HSproutingμ20.9231.30173.7434.46
σ13.6510.3887.0714.50
Floweringμ51.4198.21292.1877.36
σ14.0016.6729.4712.72
Senescenceμ-103.13316.3979.82
σ-20.0322.5615.71
LAISproutingμ0.770.983.770.93
σ1.090.881.861.43
Floweringμ3.894.185.493.86
σ1.631.781.631.68
Senescenceμ-3.623.994.86
σ-1.761.011.13
βSproutingμ0.740.720.190.78
σ0.260.180.180.25
Floweringμ0.150.150.050.15
σ0.140.160.060.13
Senescenceμ-0.210.080.05
σ-0.190.070.06
EXGSproutingμ0.610.670.580.54
σ0.120.140.250.29
Floweringμ0.810.710.540.82
σ0.260.100.320.09
Senescenceμ-0.590.530.60
σ-0.220.320.06
EXGRSproutingμ0.240.210.510.34
σ0.200.090.170.17
Floweringμ0.720.450.400.58
σ0.140.120.160.12
Senescenceμ-0.330.250.13
σ-0.150.170.10
NDISproutingμ0.480.560.510.45
σ0.140.140.190.21
Floweringμ0.790.580.520.68
σ0.090.100.130.12
Senescenceμ-0.530.400.34
σ-0.230.210.08
meanGSproutingμ0.340.340.370.35
σ0.030.030.030.03
Floweringμ0.410.340.350.39
σ0.030.030.030.05
Senescenceμ-0.330.330.30
σ-0.030.030.02
medianGSproutingμ0.380.380.410.37
σ0.030.040.040.05
Floweringμ0.450.390.380.44
σ0.020.040.050.04
Senescenceμ-0.360.360.34
σ-0.030.040.04
Individual plants were taken from the ground up and averaged in two different ways, because alfalfa is a perennial, while cotton, maize, and rice are annuals. Cotton, maize, and rice AWB was determined by weighing the replicate plants with a spring scale (precision = 0.1 g) in the field and multiplying the average by the number of plants in the quadrat from which the samples were taken. On the same day, a Nikon DX-90 digital camera® was used to collect top-of canopy RGB light intensity at nadir during ±2 h of solar noon. As with all light-sensitive measurements, this reduced errors in shadowing and other solar illumination effects. The digital photos were taken at a fixed height (1.5 m for alfalfa, cotton, and rice and 2.5 m for maize) and cropped to fit the boundaries of the quadrat. These photos were then used to estimate the number of cotton, maize, and rice plants during the sprouting phase, when it was easy to distinguish between plants. This number was decreased over the season, according to the number of samples removed from the quadrat. Since alfalfa grows in clusters, replicate clusters were taken, along with their length and width, assuming that the clusters grow in ellipses. The clusters were weighed and averaged, as above, however, the number of clusters in a quadrat was estimated by dividing a square-meter times FVC (estimated visually from the camera independent of the methods described in Section 2.2.1) by the average area of the ellipse.
Crop height was measured with a measuring tape for short canopies (alfalfa, cotton, and rice) and a telescoping measuring rod for maize (Figure 2A).
Figure 2. Non-spectroradiometric and spectroradiometric measurements used to predict aboveground wet biomass: (A) height (cm) using a measuring stick or telescoping rod; (B) fraction of vegetation cover derived from red-green-blue photos (in this case the excess green minus excess red index); (C) gap fraction and leaf area index derived from a ACCUPAR LP-80 ceptometer; and (D) hypersepctral narrowbands retrieved from an Analytical Spectral Devices Field Spec Pro 3 and white reflectance mounted to a tripod.
Figure 2. Non-spectroradiometric and spectroradiometric measurements used to predict aboveground wet biomass: (A) height (cm) using a measuring stick or telescoping rod; (B) fraction of vegetation cover derived from red-green-blue photos (in this case the excess green minus excess red index); (C) gap fraction and leaf area index derived from a ACCUPAR LP-80 ceptometer; and (D) hypersepctral narrowbands retrieved from an Analytical Spectral Devices Field Spec Pro 3 and white reflectance mounted to a tripod.
Remotesensing 07 00808 g002

2.2.1. Fraction of Vegetation Cover (FVC)

The fraction of absorbed photosynthetically active radiation was estimated using indices derived from RGB digital numbers. An automated threshold-based approach described in [42] was first used to derive FVC, but showed poor performance compared to the indices discussed below. The three indices evaluated in this study were:
NDI = (G − R)/(G + R)
EXG = 2G − R − B
EXGR = EXG − (1.4R − G)
where NDI is the normalized difference index [43], EXG is the excess greenness index [44], and EXGR is the excess green minus excess red index [45]. Indices were calculated for each pixel in the cropped images. The green and red variables shown in Equations (1)–(3) are expressed as chromatic coordinates. Binary images were generated for each index: vegetation pixels received a value of one and background pixels received a value of zero. The Otsu threshold method [46] was used to create the binary images for NDI and EXG. The threshold technique assumes the distribution of the image is bimodal and uses the index histogram to optimize discrimination between vegetation and background. In order to better distinguish image background from vegetation and to reduce camera bias, [47] developed EXGR. By subtracting the redness index from EXG shown in Equation (3), a threshold technique is not required: values less than zero are background, while values greater than zero are vegetation (Figure 2B). The mean (meanG) and median (medianG) chromatic greenness of each image was also calculated and included as relative predictors, because the former has shown to be superior to EXG in forested canopies [48], while the latter is less sensitive to outliers than the mean.

2.2.2. Gap Fraction (β) and Leaf Area Index (LAI)

The gap fraction and LAI were estimated using an ACCUPAR LP-80® ceptometer [49]. The LP-80 is comparable to other commonly used ceptometers, such as LAI-2000® and SunScan®, when evaluated against destructive measurements of LAI in maize fields under a wide array of sky conditions [50]. The LP-80 is connected to a datalogger and consists of a probe containing 80 independent photosensors, 1 cm apart (Figure 2C). The photosensors measure PAR from 400 to 700 nm in units of μmol∙m−2∙s−1 [51]. An external PAR sensor was mounted to a pole that extended above the canopy and connected to the datalogger, so that below and above canopy irradiance readings could be taken simultaneously. Ten above and below canopy measurements were taken for each replicate sample and averaged internally. The leaf area index or the area of leaves per unit area of soil surface was measured by inverting the below and above canopy irradiance. Solar zenith angle and fractional beam readings were internally computed and used along with a leaf area distribution parameter (χ) for the inversion. χ was computed at the beginning of daily ceptometer readings. χ describes how leaves of a plant are oriented. Planophile crops have χ > 1, while erectophile crops have χ < 1. χ was approximated as the ratio of the vertical gap fraction to the horizontal gap fraction determined from shadows of a representative canopy segment projected on a white gridded screen. The gap fraction was computed as follows:
β = Ib/Ia
where β (1 − FAPAR) was determined using ceptometer readings of above canopy irradiance (Ia) and below canopy irradiance (Ib).

2.3. Spectral Measurements and Processing

Spectral data were collected using an Analytical Spectral Devices (ASD) Field Spec Pro 3® [52]. A full description of spectroradiometric energy retrieval procedures and analysis can be found in [53]. Spectra consisted of energy between 350 and 2500 nm, were collected using an 18° field of view foreoptic connected to a pistol grip, and were transmitted to the spectroradiometer witha fiber optic cable (Figure 2D). The Field Spec Pro 3 yielded spectra at 1 nm intervals resampled from 3 nm beyond 700 nm and 10 nm beyond 1400 nm. The foreoptic was mounted to (at) the same pole (height) as the camera, so that spectra and RGB photos were taken simultaneously. To improve spatial and spectral uniformity across quadrats (see [54]), five replicate spectra were taken at random positions over each quadrat. Spectral data were collected for each quadrat in which in situ non-spectral measurements were made. To further increase the signal-to-noise ratio (S/N), each replicate spectrum consisted of internally averaged spectra: 30 spectra for clear sky conditions and 40 spectra for sub-optimal sky conditions. The Field Spec Pro 3 collected raw radiance (W∙sr−1∙m−2), which was standardized as surface reflectance (0–100%) using reference or “white” spectra and accompanying software (ViewSpec Pro®). Depending on sky conditions, white reflectance was collected every 2–10 min using a panel composed of BaSO4. Re-optimization was frequently performed as well, because current generated from the spectroradiometer (i.e., “dark current”) decreased the S/N over the day.
The ASD Field Spec Pro 3 fiber optic cable collected light and diverted it to three detectors, so inconsistencies in overlapping spectral regions did occur. To reduce these minor inconsistencies, the visible and NIR and SWIR2 portions of the replicate spectra were normalized using simple ratios with the endpoints of the SWIR1 portion of each replicate spectra. Each replicate spectrum group of five was then averaged after spectra exceeding one standard deviation of the group were omitted. Generally, one outlier was removed in every group. These processing steps yielded 576 spectra and AWB pairs for model-building. Regions of the spectra where atmospheric O3, H2O, and CO2 seriously degrade S/N (350–390 nm, 1350–1450 nm, 1790–2000 nm, and 2300–2500 nm) were masked before the 1 nm spectra were aggregated by simple averaging to 10 nm intervals. Previous studies [32,33] have performed this operation to reduce the number of redundant predictors. These two processes yielded a total of 196 bands, 10 nm wide, and over 350–2500 nm, which will be identified by their median wavelengths for the remainder of the paper.

2.3.1. Hyperspectral Data Transformations

The spectra were initially analyzed in four formats: untransformed, inverse-log, 1st order derivative transformed and 2nd order derivative transformed. Transformations and all further analytical procedures discussed were performed in R® [55]. The inverse-log transformation creates a pseudo-absorbance curve and has been used in a limited number of hyperspectral studies (e.g., [56]). First and 2nd order derivative transformations [57] are typically applied to spectra to amplify the vegetation signal, as changes in soil and other background spectra are more gradual than changes in the vegetation signal. The derivative signals are particularly strong at inflection points, such as the red-edge. Cubic spline functions in the “stats” package were used to smooth the spectra before 1st and 2nd order derivatives were computed for each HNB. The untransformed and 1st order derivative transformed spectra are presented here, because the other transformations were evaluated and performed worse.

2.3.2. Data Reduction for Model-building

The stepwise regression approach discussed in Section 2.4.1 requires that the number of predictors (p) is less than the number of samples (N), so a selection process was used a priori to remove predictors with relatively limited explanatory power. Other data reduction techniques were tried such as Principal Components Analysis (PCA), but yielded poorer results during validation (not shown). In addition, [53] found that PCA regression which uses the components of PCA for predictive purposes performed more poorly than other empirical methods using the same spectral data. In our procedure, a cutoff threshold was determined based on the highest ranked predictors: p values were ranked from highest to lowest based on their correlation with AWB and only the N – 1 highest ranked p was selected for further analysis.

2.4. Model Building

2.4.1. Multiple-Band Vegetation Index (MBVI)

The reduced list of untransformed or transformed HNBs and non-spectral (H, β, EXG, EXGR, NDI, meanG, and medianG) predictors were evaluated using stepwise regression to identify the optimal MBVI for each crop:
MBVI = i = 1 j m i X i + b
where i is the number of predictors used to build the model, j is the maximum number of predictors, Xi is the untransformed or transformed reflectance or non-spectral model input, mi is the respective predictor slope, and b is the model intercept.
Seventy percent of the data were selected at random and stratified by crop type to build the models, while the remaining 30% of the data were held back for validation. The biomass data were inherently non-linear, so they were log (base 2) transformed to facilitate linear model-building. The models were developed for each crop, as an “all crops” model showed significantly lower correlations with AWB. The linear models were developed using the “regsubsets” function in the “leaps” package. The regsubsets function performs a stepwise search for the linear model that minimizes the Bayesian Information Criterion (BIC). The BIC, as with the Akaike Information Criterion, penalizes over-fitting when new predictors are introduced [58]. The regsubset function identifies the optimal (minimum BIC) model using one or more predictors. The results are displayed in a convenient tabular format up to j, which is defined a priori. The maximum number of predictors was set using the ratio of the number of bands (M) to the number of samples (N) or M/N criterion [32]. When the M/N ratio exceeds a threshold of 0.15, over-fitting is likely to occur. The sample size for the calibration subset was 71, 103, 100, and 53 for alfalfa, cotton, maize, and rice, respectively. If a sample had any missing predictor values, it was omitted to maintain consistency across samples. Based on the M/N criterion, j = 9, 11, 11, and 6 for alfalfa, cotton, maize, and rice were used in the regsubsets function, respectively.

2.4.2. Two-Band Vegetation Index (TBVI)

The two-band vegetation index approach was handled differently from the single band approach. Two channels were selected and combined analogous to the Normalized Difference Vegetation Index [59]:
TBVI = R j R i R j + R i
where Ri is one untransformed or transformed HNB and Rj is another untransformed or transformed HNB. The TBVIs were ranked from highest to lowest correlated with AWB, again to eliminate relatively low explanatory predictors, so that the sample size requirement was met. The optimal combination of TBVIs and non-spectral predictors were determined using the regsubsets function as above.

2.4.3. Hyperspectral Narrowband versus Broadband Predictors

There are spectroradiometers that are considerably less expensive than the ASD Field Spec Pro, because they collect energy over broad ranges commensurate with multi-spectral broadband spaceborne sensors. Skye® [60], for example, designs sensors that measure ground reflectance for comparison with MODIS bands 1–7 (620–670, 841–876, 459–479, 545–565, 1230–1250, 1628–1652, and 2105–2155 nm). In order to evaluate the potential of multispectral broadband detectors such as this, the sample spectra were convolved over the MODIS ranges in ViewSpec Pro. These bands were combined with the non-spectral predictors to build MBVIs and TBVIs as above. In each case, however, the MODIS bands added little (∆R2 < 0.01) and statistically insignificant explanatory power when combined with non-spectral predictors and were therefore not included in the comparison with HNBs.

2.5. Model Validation

The final model selected from the calibration subset was used to estimate AWB for the validation subsets. The final model was selected that met the M/N criterion. The M/N criterion tended to yield models with fewer degrees of freedom that did not explain additional AWB variability, so a ΔR2 criterion was developed after visual inspection of the standardized residual plots (not shown) for models using the M/N criterion alone. The ΔR2 criterion, defined here as the first model to explain less than an additional one percent of AWB variability from the previous step of the stepwise regression, was selected for the validation phase. The ΔR2 criterion yielded models with high explanatory power and lower risk of over-fitting. The models developed from the calibration dataset in this manner were evaluated graphically and numerically using the coefficient of determination (R2) and the root mean square error (RMSE) with the validation subset.

3. Results

3.1. Non-Spectral Measurement Summary

Aboveground wet biomass across the full range of sample values varied from 1984–16,132, 776–9447, 7558–11,857, and 774–2727 g∙m−2 for alfalfa, cotton, maize, and rice, respectively (Table 1). Alfalfa is primarily used for cattle feed and is cut at the beginning of the flowering stage to maximize nutrient content, so it does not have a senescence phase. Even so, alfalfa yielded the largest AWB (16,132 g∙m−2) of any of the crops during senescence. Rice yielded the lowest end-of-season AWB (2727 g∙m−2). The standard deviation in AWB for the other crops increased substantially from the sprouting to flowering/tasseling and then bud/grain-filling stages. Conversely, maize sprouting means were closer to flowering (tasseling) than the other crops, because samples were taken for each crop around the same time and it has a longer growing season and is planted relatively earlier than the others. Maize overall was the tallest crop measured (mean = 316.39 cm during grain-filling), while rice was the shortest (mean = 79.82 cm during grain-filling). Height consistently increased with crop age over the period measured. Furthermore, with the exception of maize, variability increased with crop age, though not as dramatically as AWB. The leaf area index was less consistent with crop age, particularly for cotton and maize. For cotton and maize, LAI decreased from 4.18 and 5.49 to 3.62 and 3.99 from flowering/tasseling to bud/grain-filling, respectively. Similarly, β decreased for all crops between sprouting and flowering/tasseling phases, but increased from 0.15–0.21 and 0.05–0.08 for flowering/tasseling and grain/bud-filling cotton and maize, respectively. The results of the RGB derived FVC were less consistent than LAI and β. As with LAI and β estimates, cotton and maize FVC increased between sprouting and flowering/tasseling, but decreased between flowering/tasseling and bud/grain-filling. With maize, however, this was only the case for NDI—the other two indices showed a decrease between the first two stages as well. Rice also showed a decrease in FVC between flowering and grain-filling.

3.2. Spectral Measurement Summary

Figure 3 shows the strength and direction of relationships between AWB and non-spectral and untransformed spectral predictors using Pearson R correlations. The dotted lines show the cutoff for predictors that had high Pearson R correlations and met the sample size requirement for model building. Non-spectral predictors and HNBs that fall between the two lines were not included in the model-building process. The cutoff for untransformed rice (R = |0.73|) constrained spectral predictors to the NIR (783–1286 nm) and included H, LAI, and β as non-spectral predictors (Figure 3A). Height (R = 0.93) showed the strongest positive correlation with AWB, while β (R = −0.90) showed the strongest negative correlation. The highest and lowest correlated HNBs occurred at 1094 nm (R = 0.81) and 428 nm (R = −0.67), respectively. The cutoff for untransformed alfalfa (R = |0.66|) similarly included NIR HNBs, but due to its larger sample size, also included visible (428–682 nm) and SWIR1 (1437–1538 nm) (Figure 3B). Height, LAI, and β, as with rice, were included in the model-building phase; however, several FVCs were also included. The excess green minus excess red index (R = 0.87) showed the highest positive correlation, while the HNB at 672 nm (R = −0.81) showed the highest negative correlation with AWB. The HNB with the highest alfalfa biomass correlation occurred at 763 nm (R = 0.76). Cotton (cutoff R = |0.58|) included similar visible, NIR, and SWIR1 HNBs to alfalfa, and a very narrow range of significant correlations in the SWIR2 around 1992 nm (Figure 3C). As with the other crops, strong correlations were found between AWB and H, LAI, and β. Like alfalfa, EXGR also showed strong correlations with AWB. Height (R = 0.89) had the strongest positive correlation with biomass, while β had the strongest negative correlation (R = −0.90). The highest positively correlated HNB occurred at 1084 nm (R = 0.69), while the highest negatively correlated HNB occurred at 672 nm (R = −0.82). Maize (cutoff R = |0.13|) overall showed the lowest correlations and included predictors primarily in the NIR range (Figure 3D). Height, LAI, and β showed much stronger correlations with AWB than the HNBs, with H showing the highest positive correlation (R = 0.81) and β showing the highest negative correlation (R = −0.74). The best positive and negatively correlated HNB predictors were at 1084 nm (R = 0.42) and 428 nm (R = −0.16), respectively.
Figure 4 shows the strength and direction of relationships between biomass and non-spectral and 1st order derivative transformed spectral predictors using Pearson R correlations. Sample sizes for cutoff thresholds were the same, while Pearson thresholds (R = |0.70|, |0.72|, |0.53|, and |0.20|) for rice, alfalfa, cotton, and maize, respectively, were comparable. Correlations in general were over much narrower wavelength ranges than the untransformed HNBs. Correlations between rice AWB and HNBs were primarily at inflection points in the NIR, though important inflection points in the SWIR1 were above the cutoff threshold as well (Figure 4A). The maximum correlation occurred at 1245 nm (R = 0.85) and the minimum correlation occurred at 1336 nm (R = −0.88). As with the untransformed HNBs, alfalfa (Figure 4B) and cotton (Figure 4C) exhibited strong relationships in the visible, NIR, and SWIR1. The maximum and minimum correlation for alfalfa were at 712 nm (R = 0.84) and 560 nm (R = −0.86), while the maximum and minimum correlation for cotton were at 1235 nm (R = 0.86) and 1145 nm (R = −0.85). For maize, correlations improved significantly after the transformation, but were still considerably lower than the untransformed HNBs. Overall, SWIR2 remained the least significant spectral region, while NIR remained the most significant spectral region for AWB estimation (Figure 4D). The maximum and minimum correlations occurred at 773 nm (R = 0.66) and 1336 nm (R = −0.52).
Figure 3. Pearson correlation between aboveground wet biomass and candidate non-spectral and untransformed spectral predictors for (A) rice, (B) alfalfa, (C) cotton, and (D) maize. The dashed line shows the threshold used to include predictors in the model-building phase. LAI is the leaf area index, β is the gap fraction, EXG is the excess greenness index, EXGR is the excess green minus excess red index, NDI is the normalized difference index, meanG is the mean chromatic greenness index, and medianG is the median chromatic greenness index. The hyperspectral narrowbands are in nm.
Figure 3. Pearson correlation between aboveground wet biomass and candidate non-spectral and untransformed spectral predictors for (A) rice, (B) alfalfa, (C) cotton, and (D) maize. The dashed line shows the threshold used to include predictors in the model-building phase. LAI is the leaf area index, β is the gap fraction, EXG is the excess greenness index, EXGR is the excess green minus excess red index, NDI is the normalized difference index, meanG is the mean chromatic greenness index, and medianG is the median chromatic greenness index. The hyperspectral narrowbands are in nm.
Remotesensing 07 00808 g003aRemotesensing 07 00808 g003b
Figure 4. Pearson correlation between aboveground wet biomass and candidate non-spectral and 1st order derivative transformed spectral predictors for (A) rice, (B) alfalfa, (C) cotton, and (D) maize. The dashed line shows the threshold used to include predictors in the model-building phase. LAI is the leaf area index, β is the gap fraction, EXG is the excess greenness index, EXGR is the excess green minus excess red index, NDI is the normalized difference index, meanG is the mean chromatic greenness index, and medianG is the median chromatic greenness index. The hyperspectral narrowbands are in nm.
Figure 4. Pearson correlation between aboveground wet biomass and candidate non-spectral and 1st order derivative transformed spectral predictors for (A) rice, (B) alfalfa, (C) cotton, and (D) maize. The dashed line shows the threshold used to include predictors in the model-building phase. LAI is the leaf area index, β is the gap fraction, EXG is the excess greenness index, EXGR is the excess green minus excess red index, NDI is the normalized difference index, meanG is the mean chromatic greenness index, and medianG is the median chromatic greenness index. The hyperspectral narrowbands are in nm.
Remotesensing 07 00808 g004aRemotesensing 07 00808 g004b

3.3. Model-Building

3.3.1. MBVI

Each MBVI selected for validation, included both spectral and non-spectral predictors (Table 2). The model predictors and intercepts were significant at p = 0.05. Height was the most important non-spectral predictor of AWB, followed by β. Leaf area index was also significant, but exhibited strong multicollinearity with β and was excluded from all of the final models selected for validation. Although EXG, EXGR, and NDI were moderately correlated with rice, alfalfa and cotton, and maize AWB, none of the absolute FVC methods were included in the final models. Spectral predictors from the NIR were the most important for AWB estimation. Rice included HNBs in the 900–1200 nm range. The final model selected was comprised of two untransformed HNBs (963 and 993 nm) and H. The two HNBs only explained an additional 4% of the calibration subset variance after H. More variables in the NIR could have been included, but yielded little additional explanatory power (ΔR2 < 0.01). Alfalfa was the only crop with the optimal MBVI model to include 1st order derivative transformed HNBs. These transformed HNBs were primarily in the visible blue and red portions of the spectrum. It was the only crop that included a relative FVC metric in the final model. Cotton only included three predictor variables: one HNB in the visible green (539 nm), one HNB in the NIR (1134 nm), and H. It explained 88% of the calibration subset variance. Height and β were the primary predictors for maize, as the final model included three NIR HNBs (794, 845, and 865 nm) that only explained an additional 6% of the calibration subset variance.
Table 2. The multiple-band vegetation index with the highest R2 and lowest Bayesian Information Criterion at each step in the regsubset function. H is height (cm), LAI is the leaf area index (m2∙m−2), β is the gap fraction, EXGR is the excess green minus excess red index, and meanG is the mean chromatic greenness. Hyperspectral narrowbands are preceded by “λ”. First derivative transformed bands are identified with ('). Only the first five models are displayed for convenience. X, m, and b indicate the predictor, slope, and intercept of each model. The significance (p) is expressed as (***) for <0.001 and (**) for ≥0.001 and <0.05. Models shaded in gray were selected based on the (1) M/N and (2) ΔR2 criteria (Section 2.5) and plotted in Figure 5.
Table 2. The multiple-band vegetation index with the highest R2 and lowest Bayesian Information Criterion at each step in the regsubset function. H is height (cm), LAI is the leaf area index (m2∙m−2), β is the gap fraction, EXGR is the excess green minus excess red index, and meanG is the mean chromatic greenness. Hyperspectral narrowbands are preceded by “λ”. First derivative transformed bands are identified with ('). Only the first five models are displayed for convenience. X, m, and b indicate the predictor, slope, and intercept of each model. The significance (p) is expressed as (***) for <0.001 and (**) for ≥0.001 and <0.05. Models shaded in gray were selected based on the (1) M/N and (2) ΔR2 criteria (Section 2.5) and plotted in Figure 5.
RiceXmbR2ΔR2pAlfalfaXmbR2ΔR2p
1Height0.057.780.85***1EXGR7.008.340.75***
2Height0.037.640.880.03***2EXGR4.7310.320.780.04***
λ10032.37***λ428'−4257.59***
3Height0.047.580.890.01***3EXGR5.4210.440.810.03***
λ963−42.96**λ438'−14925.85***
λ99345.05**λ478'12192.62***
4Height0.037.980.900.01***4meanG22.144.810.820.01***
λ1165−18.56**λ631'−5177.17**
λ100354.20***λ468'15,276.97***
λ943−34.59***λ438'−16,395.74***
5Height0.047.680.910.00***5meanG15.566.190.830.01***
λ855200.50***λ631'−13,701.27***
λ824−956.69***λ621'26,155.46***
λ814829.14**λ580'−8975.98***
λ783−70.91**λ438'−9906.12***
CottonXmbR2ΔR2pMaizeXmbR2ΔR2p
1Β−6.1613.560.83***1Height0.0111.740.63***
2Height0.049.900.870.04***2Height0.0012.460.680.05***
λ672−17.13***β−1.60***
3Height0.038.810.880.01***3Height0.0011.320.690.01***
λ11345.88***LAI0.09***
λ539−25.70***λ7430.98**
4Height0.028.930.900.01***4Height0.0111.410.710.02***
λ112443.31***LAI0.10***
λ953−40.81***λ1114−12.30**
λ672−21.77***λ96313.39**
5Height0.038.730.900.01***5Height0.0012.560.740.03***
λ11344.73***β−1.70***
λ702−85.48***λ865−422.65***
λ692361.94***λ845560.70***
λ682−298.51***λ794−137.58***

3.3.2. TBVI

The selected TBVIs yielded correlations to AWB that were equal to or higher than the MBVIs on the calibration subset (Table 3). Again, the model predictors and intercepts were significant at p = 0.05. Height explained the most AWB variance of any non-spectral variable, while β (LAI) and relative FVC metrics (meanG and medianG) explained significant, but less variance for maize and alfalfa, respectively. In each case, however, the added performance of the TBVIs was small compared to the added number of predictors used. For the rice model, as before, H was the most important predictor, followed by NIR HNBs (753, 993, 1235, and 1256 nm). The TBVI included two more predictors than the number of predictors included in the MBVIs, but yielded the same explanatory power. The 1st order derivative transformed HNBs were selected for the final alfalfa TBVI. Narrowbands were in the visible blue (438, 458, and 499 nm) and SWIR1 (1508 and 1528 nm). Eight versus four predictor variables in the MBVI explained an additional 7% of calibration subset variance. Cotton included four HNBs in the visible green (539 and 560 nm) and NIR (943 and 963 nm) and explained only an additional 2% of the calibration subset, versus the MBVI selected that consisted of three predictors. The maize model (R2 = 0.73) performed worse than the MBVI, while including three additional variables in the visible blue and NIR.
Table 3. The two-band vegetation index with the highest R2 and lowest Bayesian Information Criterion at each step in the regsubset function. Only the first five models are displayed for convenience. H is height (cm), LAI is the leaf area index (m2∙m−2), β is the gap fraction, EXGR is the excess green minus excess red index, and meanG is the mean chromatic greenness. Hyperspectral narrowbands are preceded by λ. First derivative transformed bands are identified with ('). Only the first five models are displayed for convenience. X, m, and b indicate the predictor, slope, and intercept of each model. The significance (p) is expressed as (***) for <0.001 and (**) for ≥0.001 and <0.05. Models shaded in gray were selected based on the (1) M/N and (2) ΔR2 criteria (Section 2.5) and plotted in Figure 6.
Table 3. The two-band vegetation index with the highest R2 and lowest Bayesian Information Criterion at each step in the regsubset function. Only the first five models are displayed for convenience. H is height (cm), LAI is the leaf area index (m2∙m−2), β is the gap fraction, EXGR is the excess green minus excess red index, and meanG is the mean chromatic greenness. Hyperspectral narrowbands are preceded by λ. First derivative transformed bands are identified with ('). Only the first five models are displayed for convenience. X, m, and b indicate the predictor, slope, and intercept of each model. The significance (p) is expressed as (***) for <0.001 and (**) for ≥0.001 and <0.05. Models shaded in gray were selected based on the (1) M/N and (2) ΔR2 criteria (Section 2.5) and plotted in Figure 6.
RiceXMbR2ΔR2pAlfalfaXmbR2ΔR2p
1H0.057.780.85 ***1λ1588', λ438'−5.998.030.78 ***
2H0.037.840.880.04***2λ1588', λ428'−4.742.890.830.04***
λ1266, λ1165−42.08 *** meanG15.84 ***
3H0.038.250.890.01***3λ1588', λ428'−5.38 0.850.02***
λ993, λ753−2.68 *** meanG49.29 ***
λ1256, λ1235−73.92 ** medianG−39.41 ***
4λ1588', λ428'−3.969.300.860.01***
meanG49.84 ***
medianG−48.77 ***
λ499', λ458'−2.72**
5λ1528', λ438'−16.829.240.890.02***
meanG46.66 ***
medianG−44.81 ***
λ499', λ458'−4.46 ***
λ1508', λ448'14.98 ***
CottonXMbR2ΔR2pMaizeXmbR2ΔR2p
1λ1124, λ550−15.091.070.83***1H0.0111.650.63***
2λ1124, λ539−9.332.960.900.07***2H0.0012.420.680.05***
H0.03***β1.72***
3λ943, λ539−78.35−0.570.900.01***3λ1034, λ42868.1413.990.700.01***
λ963, λ56066.64***λ933, λ42865.48***
H0.03***H0.01***
4λ973, λ550130.150.810.910.01***4λ1003, λ428258.2114.410.720.02***
λ1104, λ57060.57***λ983, λ428255.45***
λ993, λ560180.11***H0.01***
H0.02***LAI0.08***
5λ943, λ539243.811.030.910.01***5λ1003, λ428233.3314.840.730.01***
λ953, λ560235.36***λ983, λ428−251.92***
λ1104, λ570106.05***λ1155, λ42822.15***
λ1064, λ529105.00***H0.01***
H0.02***LAI0.08***

3.4. Model Validation

With the exception of cotton, the MBVIs (Figure 5) performed equally well or better and with lower RMSE than the TBVIs in the validation subset (Figure 6). The cotton and maize distributions appeared nearly bimodal, revealing that more intermediate samples could be introduced in the future to reduce the residuals. The TBVI for cotton explained 91% versus 87% of the validation subset (N = 36) variance. The TBVI for rice explained two percent less variance in the validation subset (N = 37), while the TBVI for alfalfa showed no improvement in the validation subset (N = 37). Although the maize models related well to the calibration subset, they performed only moderately well on the validation subset (N = 38) and worse than the MBVIs.
Figure 5. Scatterplots of observed (validation subset) versus predicted aboveground wet biomass using the multiple-band vegetation indices shaded in gray from Table 2 for rice (A); alfalfa (B); cotton (C); and maize (D). The diagonal line represents a 1:1 relationship. With the exception of alfalfa, which is built on 1st order derivative transformed spectra, equations are built on the untransformed spectra.
Figure 5. Scatterplots of observed (validation subset) versus predicted aboveground wet biomass using the multiple-band vegetation indices shaded in gray from Table 2 for rice (A); alfalfa (B); cotton (C); and maize (D). The diagonal line represents a 1:1 relationship. With the exception of alfalfa, which is built on 1st order derivative transformed spectra, equations are built on the untransformed spectra.
Remotesensing 07 00808 g005
Figure 6. Scatterplots of observed (validation subset) versus predicted aboveground wet biomass using the two-band vegetation indices shaded in gray from Table 3 for rice (A); alfalfa (B); cotton (C); and maize (D). The diagonal line represents a 1:1 relationship. With the exception of alfalfa, which is built on 1st order derivative transformed spectra, equations are built on untransformed spectra.
Figure 6. Scatterplots of observed (validation subset) versus predicted aboveground wet biomass using the two-band vegetation indices shaded in gray from Table 3 for rice (A); alfalfa (B); cotton (C); and maize (D). The diagonal line represents a 1:1 relationship. With the exception of alfalfa, which is built on 1st order derivative transformed spectra, equations are built on untransformed spectra.
Remotesensing 07 00808 g006aRemotesensing 07 00808 g006b

4. Discussion

Combining the explanatory power of non-spectral and spectral predictors improved the accuracy of AWB estimation for each crop versus models built from spectral predictors alone. Hyperspectral narrowbands and non-spectral predictors exhibited multicollinearity, meaning that HNBs were an important predictor of AWB (see [53]), but were less significant in the presence of non-spectral predictors, which explained more and unique variance. The new MBVI approach resulted in an increase of 12% explained variance (R2 = 0.84 versus R2 = 0.72) versus the MBVI developed for rice using the same spectral data presented in [53] (Figure 5A). Lower but significant improvements were seen for alfalfa (Figure 5B) and cotton (Figure 5C): 6% (R2 = 0.86 versus R2 = 0.80) and 7% (R2 = 0.87 versus R2 = 0.80) in explained variance, respectively. The explained variance in maize biomass increased 29% (R2 = 0.59 versus R2 = 0.30) when including non-spectral predictors (Figure 5D), emphasizing the benefit of including non-spectral predictors in the presence of potentially mediocre spectral data (see below). The RMSE in each case was lower in this study than those reported in [53]: 1.51 g∙m−2 (rice), 2.07 g∙m−2 (alfalfa), 1.95 g∙m−2 (cotton), and 1.39 g∙m−2 (maize). The TBVIs for cotton (Figure 6C) and maize (Figure 6D) showed sizable improvements, 14% (R2 = 0.91 versus R2 = 0.77) and 38% (R2 = 0.55 versus R2 = 0.17) in explained variance, while rice (Figure 6A) and alfalfa (Figure 6B) showed only a modest 5% improvement (R2 = 0.82 versus R2 = 0.77, R2 = 0.86 versus R2 = 0.81) in explained variance. The RMSE in each case was lower in this study than those reported in [53] as well: 1.45 g∙m−2 (rice), 2.03 g∙m−2 (alfalfa), 2.05 g∙m−2 (cotton), and 1.44 g∙m−2 (maize). Multiple-band vegetation and two-band vegetation indices were developed to mimic multispectral broadbands by convolving HNBs, but proved ineffective at estimating AWB, revealing that HNBs are essential for developing AWB models that use both spectral and non-spectral predictors.
Other spectral studies, with the exception of maize, yielded correlations similar to or lower than the combination approach taken here, but should be compared with caution, due to different sampling and modeling procedures. Mariotto et al. [61] evaluated spectra against AWB estimated with a spectroradiometer and hyperspectral/broadband spaceborne sensors. The MBVI approach developed from spectrodradiometric data yielded the highest R2. For cotton (N = 150), the R2 was lower than for this study and incorporated HNBs from the NIR and SWIR1. Alfalfa (N = 9) and rice (N = 9) yielded R2s of 1.0 and 0.99 with five and three HNBs, respectively, but these should be scrutinized, due to clear over-fitting. Gnyp et al. and Atzbergeret et al. [62,63] used similar statistical methods to the ones described here and partial least squares regression to estimate AWB and leaf dry weight for rice at major phenological stages, respectively. Similar to our findings, the most highly correlated HNBs were in the NIR and the MBVI approach performed better than the TBVI approach. However, across all stages, the optimal model in the Gynp study (R2 = 0.72) used six 1st order derivative transformed HNBs to explain 12% lower variance than this study, while the optimal model in Nguyen and Lee study (R2 = 0.79) used 12 untransformed HNBs to explain 5% lower variance than the combined model that used only three untransformed predictors presented here. Moran et al. [64] used ground-based broadband reflectance collected over alfalfa in order to evaluate the sensitivity of TBVIs to changes in water status and incoming solar radiation. The most accurate and least sensitive TBVI to changes in water and light was a simple ratio of visible red: NIR. The index used morning hour reflectance and yielded an R2 = 0.86, which is comparable to the results here using three 1st order derivative transformed HNBs and one non-spectral predictor. The strength of the correlations in [64] could be influenced, however, by the small size of their dataset and the presence of leverage points in their correlation plots. Hyperspectral narrowbands were used to estimate AWB for cotton using 1st order derivative transformed single and multiple-band regression in [65]. The models with the highest R2, as in this study, were primarily in the visible and NIR, though SWIR1 bands showed high correlations with AWB in the Bai study as well. The highest ranked model used only one predictor at 502 nm, yielding R2 = 0.89—only two percent lower than the model selected here using four HNBs and one non-spectral predictor. Maize clearly underperformed compared to other studies. In [32], for example, similar statistical methods were used to develop a model employing visible (410 and 496 nm), red-edge (654 nm), and NIR (954 nm) HNBs, yielding an R2 = 0.78, while [61] used 27 data points to generate an MBVI including four HNBs in the NIR and SWIR (R2 = 0.96). The spectral relations over which HNBs were selected to estimate AWB differ between studies, revealing the limitation of empirical models for areas outside the region of interest. In a future study, the empirical models for AWB developed here will be compared with a more process-based, leaf-optic-canopy-reflectance model (PROSAIL: [66]) in areas outside the Central Valley of California to evaluate the robustness of the approach.
The relative insignificance of spectral predictors compared to non-spectral predictors underlines an important new avenue for airborne/spaceborne remote sensing research. Many of the non-spectral predictors cannot be scaled up to remote sensing resolution without performing intensive and extensive field surveys. One alternative would be to derive these metrics with remote sensing data. For example, H, overall the most important predictor of AWB, could be derived over large areas directly from LiDAR, which involves the measurement of light reflected from a surface illuminated by a laser [67].This information could be combined with airborne/spaceborne reflectance to yield more accurate measurements of AWB. Mariotto [61] showed that spaceborne hyperspectral sensors outperform spaceborne multispectral broadband sensors for AWB estimation. The combination of spectra collected from the upcoming HyspIRI mission with LiDAR, therefore, could be a powerful new tool for reducing uncertainty in areal biomass estimates.
Height consistently increased as AWB increases throughout the growing season, highlighting that inconsistency in non-spectral model inputs across phenological phases impacted their significance during model-building. The gap fraction and LAI, unlike H, increased and then decreased as AWB increased for maize and cotton. In addition, H, β, and LAI exhibited correlations with each other for cotton (not shown), which is why β and LAI were not included in the final models for cotton. The importance of β and LAI in the final maize models may be an indication of the quality of the spectral data collected (see below). With the exception of alfalfa, crops were irrigated on a deficit schedule, meaning water was gradually reduced and eventually eliminated beginning in the bud/grain-filling phase. Without water, plants lose their structure, as leaves dry out and turgor pressure reduces. This is particularly true for cotton and maize. For rice, the structure is maintained, because it is densely planted and flooded throughout the season. Alfalfa is on a modified schedule, as water is eliminated at flowering and harvested shortly after, so it never achieves senescence. The absolute RGB methods for FVC were highly inconsistent and appear to be poor tools for AWB prediction, and this could be due to any number of factors, such as poor image thresholding in the presence of soils and shadows. The relative measures tended to do better for alfalfa, but this could have been influenced by the method in which alfalfa data were sampled.
Overall, NIR was the most statistically significant and highly correlated spectral region for AWB estimation, particularly around 794, 845, 865, 943, 963, and 993 nm, followed by the visible at 438, 468, 539, 560, and 631 nm. The NIR and visible portions of the spectrum may explain internal biochemical processes and light geometry affecting plant and cell development unrelated to H. The SWIR1 region was much less important than the NIR and visible, and this could be due to a lack of crop water stress under irrigation. The SWIR2 region exhibited a low S/N ratio, so many of the HNBs in this region were not included in the analysis, potentially obscuring its significance. As in [53], by crop, NIR was the most statistically significant and highly correlated spectral region for rice and maize, while NIR and visible were the most significant spectral regions for alfalfa and cotton. The spectral response to maize AWB, as in [53] was poor. Due to instrumentation and time constraints, many of the maize spectra were collected at a height that was not sufficiently above the canopy. As a consequence, some spectra may have captured individual leaves, instead of the entire canopy. In the follow-up study where broadband and hyperspectral spaceborne sensors will be compared, if poor ground-based maize spectra are the culprit behind the lower performance of the AWB vegetation indices for maize, the spaceborne sensors should yield significantly higher correlations.
The multiple-band vegetation indices incorporating untransformed spectra clearly performed better than the 1st order derivative transformed spectra or the TBVI approach when non-spectral predictors were included. The performance of the MBVI approach is well documented (see [68]). In previous studies, TBVI typically involved only one HNB ratio. Here, we evaluated many two-band ratios together, but found in general that the approach performed worse and used more HNBs than the MBVI approach. Interestingly, a higher correlation between AWB and the TBVI was seen on the cotton validation subset compared to the MBVI, but only marginally. In general, however, the advantage of the TBVIs was seen primarily when only two HNBs were combined with non-spectral predictors to predict AWB, i.e., adding a second two-band ratio added little or no additional explanatory power. The 1st order derivative transformation tends to perform better than the untransformed spectra when transmission data are not available [69]. In our case, since the non-spectral predictors were the most significant predictors and explained the most AWB variance, the advantage of the 1st order derivative transformation could be less important. In addition, the threshold technique used to meet the p < N requirement for MBVI development was biased, in that highly correlated predictors that may have been highly correlated with each other were entered into the regsubsets function. In the future, an alternative data reduction technique that does not require p < N, such as partial least squares regression [70] or artificial neural network models [71], will be compared along with the techniques presented here.

5. Conclusions

This study combined the explanatory power of ground-based spectral and non-spectral predictors to demonstrate a non-destructive alternative for areal AWB estimation in order to facilitate the calibration/validation of remote sensing biomass models with field data. Non-destructive non-spectral and spectral measurements were collected over nearly 1200 quadrats per visits lasting 1–1.5 months with minimal impact to the fields under study. The subset where destructive measurements of AWB were taken in conjunction with non-destructive measurements was used to develop empirical models via (1) stepwise regression with HNBs and (2) stepwise regression with two-band HNB ratios. Based on the findings, it is recommended that field campaigns take measurements of H and HNBs collected from the visible and NIR (400–1000 nm) for effective AWB estimation on a per crop basis, while β and HNBs collected from the SWIR1 (1000–1700 nm) are useful if the research budget permits. The stepwise regression approach with untransformed HNBs was the most effective in this study; however, other empirical and process-based approaches should be compared in the future. The empirical models developed here will be used to create AWB transects over the Central Valley of California to evaluate the performance of several space-borne sensors. Although transferability is typically a concern with empirical models, it is expected that the methods described here can be used to develop AWB transects in other regions of the world to answer a number of research questions where issues of scale are important.

Acknowledgments

This work was funded primarily by the United States Geological Survey (USGS) Climate and Land Use Mission area’s Geographic Analysis and Monitoring and Land Remote Sensing programs in concert with the Mendenhall Research Fellowship Program. Undergraduate and graduate student field support was granted through a cooperative initiative between the USGS and The National Association of Geoscience Teachers. We would like to our extend special thanks to Tony Chang, Jeff Peters, and Bobbijean Freeman who collected the field data and the agricultural extentionists (Mark Keeley, Chris Greer, Cayle Little, Bob Hutmacher, and Shannon Mueller) from the University of California at Davis (UCD) who put us in contact with growers or farm managers interested in our work and permitted us to collect data on UCD managed research farms (West Side and Shafter Research and Extension Centers). We would also like to thank the growers and farm managers (Chuck Mathews, Don Bransford, Jim Casey, Blake Harlan, Mark Boyd, Steve Mellow, Danny Kirshe nmann, John Diener, Marty Roads, Joel Ackerknecht, Scott Schmidt, and Kurt Boeger). Lastly, we would like to thank Deborah Soltesz, Miguel Velasco, Larry Gaffney, and Lois Hales who administered day-to-day logistics while we were in the field.

Author Contributions

Michael Marshall and Prasad Thenkabail conceived and designed the experiments, while major remaining tasks (experiment execution, data analysis, and paper writing) were performed primarily by the former.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Houghton, R.A. Aboveground forest biomass and the global carbon balance. Glob. Chang. Biol. 2005, 11, 945–958. [Google Scholar] [CrossRef]
  2. IPCC. Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change; Cambridge University Press: Cambridge, UK, 2013; p. 1535. [Google Scholar]
  3. Valentini, R.; Arneth, A.; Bombelli, A.; Castaldi, S.; Cazzolla Gatti, R.; Chevallier, F.; Ciais, P.; Grieco, E.; Hartmann, J.; Henry, M.; et al. A full greenhouse gases budget of Africa: Synthesis, uncertainties, and vulnerabilities. Biogeosciences 2014, 11, 381–407. [Google Scholar] [CrossRef] [Green Version]
  4. Schoeneberger, M.; Bentrup, G.; Gooijer, H.; de Soolanayakanahally, R.; Sauer, T.; Brandle, J.; Zhou, X.; Current, D. Branching out: Agroforestry as a climate change mitigation and adaptation tool for agriculture. J. Soil Water Conserv. 2012, 67, 128A–136A. [Google Scholar] [CrossRef]
  5. Sellers, P.J. Canopy reflectance, photosynthesis and transpiration. Int. J. Remote Sens. 1985, 6, 1335–1372. [Google Scholar] [CrossRef]
  6. Challinor, A.J.; Ewert, F.; Arnold, S.; Simelton, E.; Fraser, E. Crops and climate change: Progress, trends, and challenges in simulating impacts and informing adaptation. J. Exp. Bot. 2009, 60, 2775–2789. [Google Scholar] [CrossRef] [PubMed]
  7. Hay, R.K.M. Harvest index: A review of its use in plant breeding and crop physiology. Ann. Appl. Biol. 1995, 126, 197–216. [Google Scholar] [CrossRef]
  8. Prince, S.D.; Haskett, J.; Steininger, M.; Strand, H.; Wright, R. Net primary production of U.S. midwest croplands from agricultural harvest yield data. Ecol. Appl. 2001, 11, 1194–1205. [Google Scholar] [CrossRef]
  9. Doraiswamy, P.C.; Moulin, S.; Cook, P.W.; Stern, A. Crop yield assessment from remote sensing. Photogramm. Eng. Remote Sens. 2003, 69, 665–674. [Google Scholar] [CrossRef]
  10. Lu, D. The potential and challenge of remote sensing-based biomass estimation. Int. J. Remote Sens. 2006, 27, 1297–1328. [Google Scholar] [CrossRef]
  11. Tucker, C.J. A critical review of remote sensing and other methods for non-destructive estimation of standing crop biomass. Grass Forage Sci. 1980, 35, 177–182. [Google Scholar] [CrossRef]
  12. Phinn, S.; Roelfsema, C.; Dekker, A.; Brando, V.; Anstee, J. Mapping seagrass species, cover and biomass in shallow waters: An assessment of satellite multi-spectral and airborne hyper-spectral imaging systems in Moreton Bay (Australia). Remote Sens. Environ. 2008, 112, 3413–3425. [Google Scholar] [CrossRef]
  13. Catchpole, W.R.; Wheeler, C.J. Estimating plant biomass: A review of techniques. Aust. J. Ecol. 1992, 17, 121–131. [Google Scholar] [CrossRef]
  14. Ganguli, A.C.; Vermeire, L.T.; Mitchell, R.B.; Wallace, M.C. Comparison of four nondestructive techniques for estimating standing crop in shortgrass plains. Agron. J. 2000, 92, 1211–1215. [Google Scholar] [CrossRef]
  15. Martin, R.C.; Astatkie, T.; Cooper, J.M.; Fredeen, A.H. A Comparison of methods used to determine biomass on naturalized swards. J. Agron. Crop Sci. 2005, 191, 152–160. [Google Scholar] [CrossRef]
  16. Flombaum, P.; Sala, O.E. A non-destructive and rapid method to estimate biomass and aboveground net primary production in arid environments. J. Arid Environ. 2007, 69, 352–358. [Google Scholar] [CrossRef]
  17. Hutchings, N.J. Spatial heterogeneity and other sources of variance in sward height as measured by the sonic and HFRO sward sticks. Grass Forage Sci. 1991, 46, 277–282. [Google Scholar] [CrossRef]
  18. Paruelo, J.M.; Lauenroth, W.K.; Roset, P.A. Estimating aboveground plant biomass using a photographic technique. J. Range Manag. 2000, 53, 190–193. [Google Scholar] [CrossRef]
  19. Serrano, J.M.; Peça, J.O.; da Silva, J.M.; Shahidian, S. Calibration of a capacitance probe for measurement and mapping of dry matter yield in Mediterranean pastures. Precis. Agric. 2011, 12, 860–875. [Google Scholar] [CrossRef]
  20. Gourley, C.; McGowan, A. Assessing differences in pasture mass with an automated rising plate meter and a direct harvesting technique. Aust. J. Exp. Agric. 1991, 31, 337–339. [Google Scholar] [CrossRef]
  21. Ehlert, D.; Hammen, V.; Adamek, R. On-line sensor pendulum-meter for determination of plant mass. Precis. Agric. 2003, 4, 139–148. [Google Scholar] [CrossRef]
  22. Ollinger, S.V. Sources of variability in canopy reflectance and the convergent properties of plants. New Phytol. 2011, 189, 375–394. [Google Scholar] [CrossRef] [PubMed]
  23. Pinter, P.J.; Hatfield, J.L.; Schepers, J.S.; Barnes, E.M.; Moran, M.S.; Daughtry, C.S.T.; Upchurch, D.R. Remote sensing for crop management. Photogramm. Eng. Remote Sens. 2003, 69, 647–664. [Google Scholar] [CrossRef]
  24. Horler, D.N.H.; Dockray, M.; Barber, J. The red edge of plant leaf reflectance. Int. J. Remote Sens. 1983, 4, 273–288. [Google Scholar] [CrossRef]
  25. Perry, E.M.; Roberts, D.A. Sensitivity of narrow-band and broad-band indices for assessing nitrogen availability and water stress in an annual crop. Agron. J. 2008, 100, 1211–1219. [Google Scholar] [CrossRef]
  26. Ustin, S.L.; Roberts, D.A.; Gamon, J.A.; Asner, G.P.; Green, R.O. Using imaging spectroscopy to study ecosystem processes and properties. BioScience 2004, 54, 523–533. [Google Scholar] [CrossRef]
  27. Goetz, A.F.H. Three decades of hyperspectral remote sensing of the earth: A personal view. Remote Sens. Environ. 2009, 113, S5–S16. [Google Scholar] [CrossRef]
  28. Carter, G.A. Primary and secondary effects of water content of the spectral reflectance of leaves. Am. J. Bot. 1991, 78, 916–924. [Google Scholar] [CrossRef]
  29. Asner, G.P. Biophysical and biochemical sources of variability in canopy reflectance. Remote Sens. Environ. 1998, 64, 234–253. [Google Scholar] [CrossRef]
  30. HyspIRI Mission Study. Available online: http://hyspiri.jpl.nasa.gov/ (accessed on 3 December 2014).
  31. Thenkabail, P.S.; Ward, A.D.; Lyon, J.G. Landsat-5 thematic mapper models of soybean and corn crop characteristics. Int. J. Remote Sens. 1994, 15, 49–61. [Google Scholar] [CrossRef]
  32. Thenkabail, P.S.; Ward, A.D.; Lyon, J.G.; Merry, C.J. Thematic mapper vegetation indices for determining soybean and corn growth parameters. Photogramm. Eng. Remote Sens. 1994, 60, 437–442. [Google Scholar]
  33. Thenkabail, P.S.; Smith, R.B.; de Pauw, E. Hyperspectral vegetation indices and their Relationships with agricultural crop characteristics. Remote Sens. Environ. 2000, 71, 158–182. [Google Scholar] [CrossRef]
  34. Thenkabail, P.S.; Enclona, E.A.; Ashton, M.S.; Legg, C.; de Dieu, M.J. Hyperion, IKONOS, ALI, and ETM+ sensors in the study of African rainforests. Remote Sens. Environ. 2004, 90, 23–43. [Google Scholar] [CrossRef]
  35. CDFA. California Agricultural Statistics Review 2012–2013; California Department of Food and Agriculture: Sacramento, CA, USA, 2013; p. 131. [Google Scholar]
  36. USDA. 2007 Census of Agriculture; United States Department of Agriculture: Washington, DC, USA, 2009; p. 739. [Google Scholar]
  37. Faunt, C.C. Groundwater Availability of the Central Valley Aquifer; U.S. Geological Survey Professional Paper No. 1766; United State Geological Survey: Sacramento, CA, USA, 2009; p. 225. [Google Scholar]
  38. California Department of Water Resources. Available online: http://www.water.ca.gov/ (accessed on 3 December 2014).
  39. California Water Science Center. Available online: http://ca.water.usgs.gov/ (accessed on 3 December 2014).
  40. USDA CropScape. Available online: http://nassgeodata.gmu.edu/CropScape/ (accessed on 3 December 2014).
  41. McCoy, R.M. Field Methods in Remote Sensing; The Guilford Press: New York, NY, USA, 2005. [Google Scholar]
  42. Liu, J.; Pattey, E. Retrieval of leaf area index from top-of-canopy digital photography over agricultural crops. Agric. For. Meteorol. 2010, 150, 1485–1490. [Google Scholar] [CrossRef]
  43. Pérez, A.J.; López, F.; Benlloch, J.V.; Christensen, S. Colour and shape analysis techniques for weed detection in cereal fields. Comput. Electron. Agric. 2000, 25, 197–212. [Google Scholar] [CrossRef]
  44. Woebbecke, D.M.; Meyer, G.E.; von Bargen, K.; Mortensen, D.A. Color indices for weed identification under various soil, residue, and lighting conditions. Trans. ASAE 1995, 38, 259–269. [Google Scholar] [CrossRef]
  45. Meyer, G.E.; Neto, J.C. Verification of color vegetation indices for automated crop imaging applications. Comput. Electron. Agric. 2008, 63, 282–293. [Google Scholar] [CrossRef]
  46. Otsu, N. A threshold selection method from gray-level histograms. Automatica 1975, 11, 23–27. [Google Scholar]
  47. Meyer, G.E.; Hindman, T.W.; Laksmi, K. Machine vision detection parameters for plant species identification. Proc. SPIE 1999, 3543. [Google Scholar] [CrossRef]
  48. Sonnentag, O.; Hufkens, K.; Teshera-Sterne, C.; Young, A.M.; Friedl, M.; Braswell, B.H.; Richardson, A.D. Digital repeat photography for phenological research in forest ecosystems. Agric. For. Meteorol. 2012, 152, 159–177. [Google Scholar] [CrossRef]
  49. Decagon Devices, Inc. Available online: http://www.decagon.com/ (accessed 3 December 2014).
  50. Wilhelm, W.W.; Ruwe, K.; Schlemmer, M.R. Comparison of three leaf area index meters in a corn canopy. Crop Sci. 2000, 40, 1179–1183. [Google Scholar] [CrossRef]
  51. Decagon. In AccuPAR PAR/LAI Ceptometer Model LP-80; Version 10; Decagon Devices, Inc.: Pullman, WA, USA, 2010.
  52. ASD Inc. Available online: http://www.asdi.com/ (accessed 3 December 2014).
  53. Marshall, M.; Thenkabail, P. Biomass modeling of four leading world crops using hyperspectral narrowbands in support of HyspIRI mission. Photogramm. Eng. Remote Sens. 2014, 80, 757–772. [Google Scholar] [CrossRef]
  54. MacArthur, A.; MacLellan, C.J.; Malthus, T. The fields of view and directional response functions of two field spectroradiometers. IEEE Trans. Geosci. Remote Sens. 2012, 50, 3892–3907. [Google Scholar] [CrossRef]
  55. The Comprehensive R Archive Network. Available online: http://cran.r-project.org/ (accessed 3 December 2014).
  56. Serrano, L.; Peñuelas, J.; Ustin, S.L. Remote sensing of nitrogen and lignin in Mediterranean vegetation from AVIRIS data: Decomposing biochemical from structural signals. Remote Sens. Environ. 2002, 81, 355–364. [Google Scholar] [CrossRef]
  57. Thorp, K.R.; Tian, L.; Yao, H.; Tang, L. Narrow-Band and derivative-based vegetation indices for hyperspectral data. Trans. ASAE 2004, 47, 291–299. [Google Scholar] [CrossRef]
  58. Hair, J.F., Jr.; Anderson, R.E.; Tatham, R.L.; Black, W.C. Multivariate Data Analysis, 5th ed.; Prentice Hall: Upper Saddle River, NJ, USA, 1998. [Google Scholar]
  59. Tucker, C.J. Red and photographic infrared linear combinations for monitoring vegetation. Remote Sens. Environ. 1979, 8, 127–150. [Google Scholar] [CrossRef]
  60. Skye Instruments. Available online: http://www.skyeinstruments.com/ (accessed on 3 December 2014).
  61. Mariotto, I.; Thenkabail, P.S.; Huete, H.; Slonecker, T.; Platonov, A. Hyperspectral versus multispectral crop- biophysical modeling and type discrimination for the HyspIRI mission. Remote Sens. Environ. 2013, 139, 291–305. [Google Scholar] [CrossRef]
  62. Gnyp, M.L.; Miao, Y.; Yuan, F.; Ustin, S.L.; Yu, K.; Yao, Y.; Bareth, G. Hyperspectral canopy sensing of paddy rice aboveground biomass at different growth stages. Field Crops Res. 2014, 155, 42–55. [Google Scholar] [CrossRef]
  63. Atzberger, C.; Guérif, M.; Baret, F.; Werner, W. Comparative analysis of three chemometric techniques for the spectroradiometric assessment of canopy chlorophyll content in winter wheat. Comput. Electron. Agric. 2010, 73, 165–173. [Google Scholar] [CrossRef]
  64. Moran, M.S.; Pinter, P.J., Jr.; Clothier, B.E.; Allen, S.G. Effect of water stress on the canopy architecture and spectral indices of irrigated Alfalfa. Remote Sens. Environ. 1989, 29, 251–261. [Google Scholar] [CrossRef]
  65. Bai, J.; LI, S.; Wang, K.; SUI, X.; Chen, B.; Wang, F. Estimating aboveground fresh biomass of different cotton canopy types with homogeneity models based on hyper spectrum parameters. Agric. Sci. China 2007, 6, 437–445. [Google Scholar] [CrossRef]
  66. Jacquemoud, S.; Verhoef, W.; Baret, F.; Bacour, C.; Zarco-Tejada, P.J.; Asner, G.P.; Francois, C.; Ustin, S.L. PROSPECT + SAIL models: A review of use for vegetation characterization. Remote Sens. Environ. 2009, 113, S56–S66. [Google Scholar] [CrossRef]
  67. Lim, K.; Treitz, P.; Wulder, M.; St-Onge, B.; Flood, M. LiDAR remote sensing of forest structure. Prog. Phys. Geogr. 2003, 27, 88–106. [Google Scholar] [CrossRef]
  68. Thenkabail, P.S.; Mariotto, I.; Gumma, M.K.; Middleton, E.M.; Landis, D.R.; Huemmrich, K.F. Selection of Hyperspectral Narrowbands (HNBs) and composition of Hyperspectral two band Vegetation Indices (HVIs) for biophysical characterization and discrimination of crop types using field reflectance and hyperion/EO-1 data. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2013, 6, 427–439. [Google Scholar] [CrossRef]
  69. Tsai, F.; Philpot, W. Derivative analysis of hyperspectral data. Remote Sens. 1998, 66, 41–51. [Google Scholar] [CrossRef]
  70. Huang, Z.; Turner, B.J.; Dury, S.J.; Wallis, I.R.; Foley, W.J. Estimating foliage nitrogen concentration from HYMAP data using continuum removal analysis. Remote Sens. Environ. 2004, 93, 18–29. [Google Scholar] [CrossRef]
  71. Combal, B.; Baret, F.; Weiss, M.; Trubuil, A.; Macé, D.; Pragnère, A.; Myeni, R.; Wang, L. Retrieval of canopy biophysical variables from bidirectional reflectance: Using prior information to solve the ill-posed inverse problem. Remote Sens. Environ. 2003, 84, 1–15. [Google Scholar] [CrossRef]

Share and Cite

MDPI and ACS Style

Marshall, M.; Thenkabail, P. Developing in situ Non-Destructive Estimates of Crop Biomass to Address Issues of Scale in Remote Sensing. Remote Sens. 2015, 7, 808-835. https://0-doi-org.brum.beds.ac.uk/10.3390/rs70100808

AMA Style

Marshall M, Thenkabail P. Developing in situ Non-Destructive Estimates of Crop Biomass to Address Issues of Scale in Remote Sensing. Remote Sensing. 2015; 7(1):808-835. https://0-doi-org.brum.beds.ac.uk/10.3390/rs70100808

Chicago/Turabian Style

Marshall, Michael, and Prasad Thenkabail. 2015. "Developing in situ Non-Destructive Estimates of Crop Biomass to Address Issues of Scale in Remote Sensing" Remote Sensing 7, no. 1: 808-835. https://0-doi-org.brum.beds.ac.uk/10.3390/rs70100808

Article Metrics

Back to TopTop