Estimation of Maize Residue Cover Using Landsat-8 OLI Image Spectral Information and Textural Features

The application of crop residue has become increasingly important for providing a barrier against water and wind erosion and improving soil organic matter content, infiltration, evaporation, temperature, and soil structure. The objectives of this work were to: (i) estimate maize residue cover (MRC) from Landsat-8 OLI images using seven vegetation indices (VIs) and eight textural features; and (ii) compare the VI method, textural feature method, and combination method (integration of textural features and spectral information) for estimating MRC with partial least squares regression (PLSR). The results showed that the normalized difference tillage index (NDTI), simple tillage index (STI), normalized difference index 7 (NDI7), and shortwave red normalized difference index (SRNDI) were significantly correlated with MRC. The MRC model based on NDTI outperformed (R = 0.84 and RMSE = 12.33%) the models based on the other VIs. Band3mean, Band4mean, and Band5mean were highly correlated with MRC. The regression between Band3mean and MRC was stronger (R = 0.71 and RMSE = 15.21%) than those between MRC and the other textural features. The MRC estimation accuracy using the combination method (R = 0.96 and RMSE = 8.11%) was better than that based on only the VI (R = 0.88 and RMSE = 11.34%) or textural feature (R = 0.90 and RMSE = 9.82%) methods. The results suggest that the combination method can be used to estimate MRC on a regional scale. OPEN ACCESS Remote Sens. 2015, 7 14560


Introduction
The application of crop residue to agricultural fields after harvest is important for protecting the soil surface against water and wind erosion [1][2][3][4].This agricultural practice also helps improve the soil structure and increase the organic matter content in soil [5], thereby facilitating the biodegradation of pollutants [6].Thus, management of crop residues is an integral part of many conservation tillage systems.The Conservation Technology Information Center (CTIC) has classified different tillage management strategies (intensive tillage, reduced tillage, and conservation tillage) based on the percentage of crop residue cover [2].Long-term use of conservation tillage practices can lead to increased soil organic carbon, improved soil structure, and increased aggregation compared with intensively tilled soils [5].Crop residue estimates are also a critical parameter for soil carbon estimations and in modeling and monitoring improvements in carbon sequestration that result from changes in land management approaches.Therefore, it is important to estimate crop residue cover (CRC) on agricultural land at a regional scale.However, the traditional survey-based methods for quantifying CRC over large areas are time consuming, unsuitable, and costly [7].
Remote sensing is a useful tool to efficiently quantify CRC on a regional scale and provides spatial information of CRC within crop fields [7][8][9][10][11].Early attempts to estimate CRC using remote sensing techniques began in 1975 [12].Over the past few decades, several methods have been used to estimate CRC in different regions based on multispectral and hyperspectral remote sensing data and spectroradiometric ground measurements of soil type, residue type, and residue percentage [7][8][9][10][11]13,14].McNairn and Protz [15] analyzed the relationship between CRC and the normalized difference index (NDI) under different soil conditions using Landsat-5 Thematic Mapper (TM) data; the results presented a strong linear relationship between NDI and maize residue cover (R 2 = 0.65-0.84).van Deventer et al. [16] studied the relationships between the simple tillage index (STI), normalized difference tillage index (NDTI), and CRC using bands five and seven of Landsat-5 TM data, and good relationships were exhibited between these parameters.Daughtry [8] proposed the cellulose absorption index (CAI) based on cellulose and lignin absorption features near 2100 nm to discriminate crop residue from soil using in situ spectra acquired with the GER 3700 spectroradiometer.Furthermore, Daughtry et al. [9] showed that CRC was linearly related to CAI (R 2 = 0.89).High estimation accuracy (R 2 = 0.91) between the normalized difference senescent vegetation index (NDSVI) and fractional senescent vegetation cover was achieved using a two-component linear unmixing analysis [17].Crop residue cover has been estimated and mapped on agricultural land using hyperspectral and IKONOS data to demonstrate the strong relationship between predicted and measured CRC (index of agreement (D) = 0.92) [10].The shortwave infrared normalized difference residue index (SINDRI) derived from ASTER bands six and seven is strongly related to CRC (R 2 = 0.74) [14].Pacheco and McNairn [7] estimated CRC based on SPOT and Landsat-5 TM data, and good consistency between the predicted and measured CRC was obtained using the spectral unmixing method.Zheng et al. [18,19] applied the minimum NDTI (minNDTI) to estimate CRC based on a time series of Landsat TM and Enhanced TM Plus (ETM+) data (R 2 = 0.89).Bocco et al. [20] estimated CRC using neural network (NN) and crop residue index multiband (CRIM) models from Landsat TM and ETM+ image data, and the results showed that both methodologies were appropriate for estimating CRC from Landsat data (R 2 = 0.95 for the NN model and R 2 = 0.87 for the CRIM model).
The textural feature method can be used to extract textural features through image processing techniques and to analyze texture either quantitatively or qualitatively [21].Image texture analysis involves measuring heterogeneity in the tonal values of pixels within a defined area of an image [22], and this method has been used to estimate parameters of vegetative structure [23,24].Empirical models have been developed to study the relationships between vegetation indices (VIs) and CRC on regional scales.However, the use of textural features for estimating CRC has not been thoroughly investigated.Darvishzadeh et al. [25] reported that using full spectral subsets or the maximum amount of spectral information available will require more computation time but is not likely to increase estimation accuracy.Moreover, several features are sufficient for obtaining the essential information to estimate LAI, nitrogen concentration, and vegetation classification [26][27][28][29].
Maize is the main crop in the Songnen Plain, which is an important food production base in China.The estimation of maize residue cover (MRC) for this crop is important for agricultural production and management in this region.The distribution of MRC throughout an agricultural field can be highly heterogeneous, so it is difficult to effectively obtain MRC at a regional scale.Moreover, previous studies have failed to clearly describe the CRC estimation accuracy using Landsat-8 OLI imagery data.Therefore, the main purpose of this study was to accurately estimate MRC using Landsat-8 OLI imagery data in the Songnen Plain.The specific objectives of this study were to: (1) investigate the potential of Landsat-8 OLI images for MRC estimation using seven VIs and eight textural features; and (2) compare the VI, textural feature, and combination (combination of VIs and textural features) methods for estimating CRC with partial least squares regression (PLSR).The results provide a guideline for estimating MRC based on spectral information and textural features on a regional scale.This method to measure percent MRC is, thus, adequate for policy and program development, land management decision-makers, and carbon modelers.

Study Area
The study area was located between Changchun (43°5′N-45°15′N, 124°18′E-127°2′E), Jilin Province, and Suihua (45°3′N-48°2′N, 124°13′ E-128°30′E) in Heilongjiang Province, Northeast China (Figure 1).This area is part of the Songnen Plain, which is the main distribution area of spring maize.Maize cultivation covers about 80% [30] of the arable land in the Songnen Plain.This area is under a typical semi-humid northern temperate monsoon continental climate.The average maximum temperature is 23.5 °C in the summer, and the average minimum temperature is −13.3 °C in the winter.The average annual precipitation ranges from 522 to 615 mm, and the average frost-free period lasts for 120 to 150 d.Local maize cultivars are planted and harvested in May and October, respectively.

Field Measurements
The MRC at field sites was measured using the line transect method [31] from 31 October to 6 November 2014 (Figure 2).We used a 100-ft (29.85 m) measuring tape, which could be easily divided into 100 parts with 1 ft intervals, shown as red markings, to make the measurements.At each sampling site, the tape was stretched diagonally across rows [32] and the number of markings intersecting crop residue was counted.Percent cover was calculated as the average of the marking numbers.In addition, we used a UniStrong G3 GPS unit (positional accuracy < 5 m) to record the location of each measurement site, acquire photographs, and record information for each site.We measured a total of 36 fields.In this study, each sample plot is an adequate representation of a much greater area (>50,000 m 2 ).The area of each sample plot was ~1781.68 m 2 .The MRC data were randomly divided into two parts using SPSS software (16.0,SPSS, Chicago, IBM, Armonk, NY, USA): a calibration dataset with 24 samples and a validation dataset with 12 samples.The calibration dataset was used to establish the regression equations between MRC and vegetation indices.The validation dataset was used to evaluate the consistency between the predicted and measured MRC.The statistics of each subset for MRC are summarized in Table 1.As shown in Table 1, the MRC for the calibration dataset ranged from 10.00% to 95.00%, with an average of 60.04% and a standard deviation of 27.19%.Similarly, the corresponding values for the validation dataset were 10.90%, 89.00%, 59.60%, and 27.69%, respectively.Table 1 showed that the calibration dataset and the validation dataset differed.The range of the calibration dataset was larger than the range of the validation dataset.The mean of the calibration dataset was higher than the mean of the validation dataset.However, the standard deviation of the calibration dataset was higher than the standard deviation of the validation dataset.

Remote Sensing Data and Pre-Processing
Two Landsat-8 OLI satellite images were obtained from the United States Geological Survey (USGS) website (http://glovis.usgs.gov/).The corresponding lines and rows were 118 and 028, and 118 and 029 on 4 November 2014.The Landsat-8 OLI satellite images were processed (including radiometric calibration, atmospheric correction, and geometric correction) using ENVI 4.7 software.Land surface reflectance was then retrieved using the Fast Line-of-Sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) module.Each Landsat-8 OLI image was co-registered with a geo-referenced Landsat TM image at a geometric accuracy of <0.5 pixels.Since the area of each sample plot was ~1781.68 m 2 , nine pixels of Landsat-8 OLI satellite image data were corrected to match the ground sample plot.

Vegetation Indices
Our study included seven vegetation indices (VIs).Six were derived from previous studies, and the remaining one was developed in this study (Table 2, [15][16][17]33]).The previous studies indicated that the six VIs were correlated with CRC.

Textural Features
Textural features were statistically obtained with a gray level co-occurrence matrix (GLCM).Many such metrics can be derived from this matrix.In this study, we selected eight textural feature metrics, including the mean, variance, homogeneity, contrast, dissimilarity, entropy, second moment, and correlation [22].We calculated the textural features of each Landsat-8 OLI image by analyzing the eight metrics with a window size of 3 × 3 pixels, which matched the scale at which our MRC data were collected from the ground.

Extraction of Maize Cultivation Area
The maximum likelihood method was used to extract the maize cultivation area from the Landsat-8 OLI satellite image data and remove the effect of built-up areas, water, forests, and roads.This technique assumes that the statistics for each class of each band are normally distributed and calculates the probability that a given pixel belongs to a specific class [34].Unless a probability threshold is selected, all pixels are classified.Each pixel is assigned to the class with the highest probability.If the highest probability is smaller than the threshold specified, the pixel will remain unclassified.The ENVI 4.7 software completed the classification according to 280 field survey points with known land cover types.The overall accuracy of the extracted maize cultivation area was 94%, which satisfied the requirements of this study.

Partial Least Squares Regression
The PLSR technique generalizes and combines regression modeling (multiple linear regressions), data structure analysis (principal component analysis, PCA), and correlation analysis between two sets of variables (canonical correlation analysis).The method is a powerful modeling tool that reduces large numbers of measured collinear VIs to a few non-correlated latent variables or factors [35,36].The noise and collinearity in the original spectra are thus significantly eliminated from the condensed components.The optimal number of factors in PLSR analysis is determined by minimizing the prediction residual error sum of the squares statistic.Compared with PLSR, PCA only summarizes the information of independent variable, and does not consider the relationships between the independent variables and dependent variables; therefore, it loses useful information when extracting the principal component.The PLSR technique not only considers the information of the independent variable, but it also considers the relationships between the independent and dependent variables.Therefore, the stability of PLSR is better than that of PCA.
PLSR models are an extension of multiple linear regression models (e.g., multiple or general stepwise regressions).This method is particularly useful when one needs to predict a set of dependent variables from a large set of independent variables.In its simplest form, a linear model specifies the (linear) relationship between a dependent (response) variable y and a set of predictor (x) variables [37]: where b0 is the regression coefficient for the intercept and the bi values are regression coefficients (for variables 1 to p) computed from the data.

Statistical Analysis
The relationships between MRC and VIs/textural features were analyzed.The coefficient of determination (R 2 ) and root mean square error (RMSE) were used to quantify the amount of variation explained by the developed relationships and the accuracy of the relationships.Generally, the performance of a model is determined by comparing the differences in the R 2 and RMSE values.Higher R 2 values correspond to lower RMSE values and a higher precision and accuracy of a model for predicting MRC [38].

Relationship between Maize Residue Cover and Vegetation Indices
We found significant relationships between the VIs and MRC based on linear and nonlinear regressions (Table 3).The NDTI-MRC and MCRC-MRC regressions had the highest and lowest R 2 values (0.84 and 0.50, respectively).The order of the VIs from highest to lowest with respect to the R 2 values of the MRC regression relationships was NDTI, STI, NDI7, SRNDI, NDI5, NDSVI, and MCRC.Of the R 2 values, one was above 0.80, four were above 0.70, and two were below 0.6.The MCRC was fit to a linear equation; NDI5, NDI7, SRNDI, and NDSVI were fit to exponential equations; and STI and NDTI were fit to power equations.These R 2 values were adjusted, as the adjusted R 2 values better demonstrate the relationships between MCR and VIs.To validate the estimation accuracy of the MRC regression equations, we compared the predicted and measured MRC values.The RMSE values ranged from 12.33% to 21.16% (Table 3).These results suggested that NDTI, STI, NDI7, and SRNDI can be used to estimate MRC (Table 3 and Figure 3).

Relationship between Maize Residue Cover and Textural Features
Eight textural feature indicators (i.e., mean, variance, homogeneity, contrast, dissimilarity, entropy, second moment, and correlation) were calculated to examine the regression relationships between each textural feature indicator and MRC.The results showed that 10 textural feature indicators were significantly correlated with MRC (Table 4).All of the textural feature indicators were highly correlated with MRC, with the exception of Band3dissimilarity and Band3homogeneity.Band3mean and Band3dissimilarity had the strongest (R 2 = 0.71) and weakest (R 2 = 0.21) relationships with MRC, respectively.The R 2 values of the 10 textural feature metrics ranged from 0.21 to 0.67, and the corresponding order of the textural feature-MRC regression relationships from lowest to highest was Band3dissimilarity, Band3homogeneity, Band3entropy, Band3second moment, Band3correlation, Band6mean, Band2mean, Band5mean, Band4mean, and Band3mean.Again, these R 2 values were adjusted, as the adjusted R 2 values better demonstrate the relationship between MCR and textural features.To validate the estimation accuracy of the MRC models, the predicted and measured values were compared.The RMSE values ranged from 15.21% to 26.35% (Table 4).The RMSE values indicate the prediction power of the estimation equation for MRC.The lower RMSE values show that the estimation equation for MRC had better prediction accuracy.These statistical criteria (R 2 and RMSE) were comprehensively considered to select the best fitting regression equation for estimating MRC.Therefore, Band3mean-Band5mean were used to estimate MRC (Table 4 and Figure 4).

Estimating MRC via Partial Least Squares Regression (PLSR)
Three methods (VI, textural feature, and combination method) for estimating MRC via PLSR were compared (Table 5).For the VI methods, NDTI, STI, NDI7, and SRNDI were selected to estimate MRC using PLSR because they had high R 2 and low RMSE values, and the R 2 and RMSE values of the regression equation were 0.87 and 11.36%, respectively.The R 2 and RMSE values of the MRC model based on all the VIs were 0.88 and 11.34%, respectively.Thus, the addition of MCRC, NDI5, and NDSVI to the model did not significantly improve the MRC estimation accuracy.For the textural feature method, Band3mean, Band4mean and Band5mean were used to estimate MRC by PLSR because they had high R 2 and low RMSE values, and the R 2 and RMSE values were 0.83 and 12.32%, respectively.The R 2 and RMSE values of the MRC regression equation based on all the textural features were 0.90 and 9.82%, respectively.The MRC estimation accuracy was further improved using PLSR based on all the textural features.For the combination method (the integration of textural features and spectral information), four VIs (NDTI, STI, NDI7, and SRNDI) and three textural feature indicators (Band3mean-Band5mean) were selected to estimate MRC via PLSR, and the R 2 and RMSE values for the model were 0.95 and 8.43%, respectively.The combination method based on all the VIs and textural feature indicators were selected to estimate MRC via PLSR, and the R 2 and RMSE values for the model were 0.96 and 8.11%, respectively.The MRC estimation accuracy was higher using the combination method than the VI and textural feature methods alone (Table 5 and Figure 5).Thus, the combination method can be used to estimate MRC.

MRC Mapping
The spatial distribution of MRC was estimated based on all VIs and texture feature indicators via the combination method (Figure 6).The MRC ranged from 7% to 98% over the entire Songnen Plain.The MRC occupying 14.92%, 19.15%, and 20.47% of the maize cultivation area ranged from 7% to 25%, from 25% to 50%, and from 50% to 75%, respectively.Of the cultivated area, 22.89% had MRC ranging from 75% to 85% and 22.57% had MRC ranging from 85% to 98%.These results were consistent with our field observations.The highest and lowest overall MRC was measured in Suihua and Changchun, respectively, and the measured MRC in Harbin was between these two values.

Discussion
In this study, the sample data was divided into a calibration set and validation set, which is a method that has commonly been used in previous studies.Many scientists have used this method to carry out related studies [39][40][41][42][43][44][45].These studies have indicated that this method can be used to obtain good estimation results.In this paper, the amount of sample data was suitable for the statistical requirements.There were some differences between the two groups, which reflect the stability of the estimation model.Finally, the results of the estimation model show that the results of this study are representative.
All of the investigated VIs were significantly correlated with MRC, particularly NDTI, STI, NDI7, and SRNDI (Table 3).Daughtry et al. [9] indicated that two major water absorption bands exist near 1450 and 1960 nm, and oven-dried and water-saturated residues have the highest and lowest spectrum, respectively.Murray and Williams [39] suggested that a broad absorption feature near 2100 nm is also evident in the reflectance spectra of all dry crop residues, which is probably associated with lignin and cellulose in crop residue.Thus, spectral information at 1450-1960 nm and near 2100 nm is highly correlated with CRC and can be used to estimate CRC [4,10,14,16,18,19].NDTI and STI were calculated from bands six (1560-1660 nm) and seven (2100-2300 nm) of the Landsat-8 OLI image data, and the NDTI-MRC and STI-MRC regressions yielded R 2 values of 0.84 and 0.78, respectively (Table 3).Our results were consistent with those of Zheng et al. [18,19] and Galloza et al. [11], indicating that it is feasible to estimate MRC using NDTI or STI from Landsat-8 OLI images.The NDI7 values were calculated from bands five (845-885 nm) and seven of the Landsat-8 OLI imagery data.Previous results have indicated that band five is highly correlated with CRC because of the sensitivity of near infrared wavelengths to plant structure [46], meaning that crop residue has a higher reflectance than bare soil in this band [15].Therefore, NDI7 was significantly related to CRC in this study.Qi et al. [17] proposed SRNDI based on bands four (630-680 nm) and seven, and our results showed that SRNDI was better correlated with CRC than NDSVI (Table 3).Daughtry et al. [8,47] indicated that the cellulose absorption index (CAI) is a better proxy than the other VIs for estimating crop residues based on cellulose and lignin absorption features near 2100 nm.This may explain why SRNDI performed better than NDSVI in this study.Our results also suggested that NDTI, STI, NDI7, and SRNDI can be used to estimate MRC on a regional scale.
Previous findings have indicated that textural features can be used to estimate LAI and improve the estimation accuracy of aboveground biomass in forests [48][49][50].In the current study, 10 textural feature metrics presented a good relationship with CRC.Band3mean, Band4mean, and Band5mean were significantly correlated with CRC (Table 4), which is consistent with the results of previous studies [22,51].Thus, the 10 textural feature indicators investigated in the current study can be used to estimate MRC, which further confirms the results of previous studies [48,50].
The results showed that NDTI, STI, NDI7 and SRNDI were highly correlated with MRC using PLSR (Table 5).This may be because the PLSR method fully considers the relationships between the covariance of VIs and MRC variables by applying data compression into regression factors.Previous studies have shown that PLSR can be used to improve the estimation accuracy of biophysical and biochemical parameters [11,52].However, our results showed that the addition of MCRC, NDI5, and NDSVI to the PLSR did not significantly improve the CRC estimation accuracy.Li et al. [26] and Mutanga et al. [27] indicated that the main spectral information (sensitive spectral bands and VIs) is sufficient to obtain high estimation accuracy, and our results were in agreement with their studies.The addition of Band2mean, Band3homogeneity, Band3dissimilarity,Band3entropy, Band3second moment, Band3correlation, and Band6mean further impoved the estimation accauracy of MRC based on Band3mean, Band4mean, and Band5mean using PLSR (Table 5).These results suggested that Band2mean, Band3homogeneity, Band3dissimilarity, Band3entropy, Band3second moment, Band3correlation, and Band6mean include useful information for MRC estimations.The results also demonstrated that textural feature indicators can be used to estimate MRC by PLSR.The results showed that the estimation accuracy of MRC was improved by the textural feature method based on all the textural feature indicators compared with the vegetation index methods based on all the VIs (Table 5).This may be because (1) the number of textural feature indicators was more than the number of VIs; (2) the addition of Band2mean, Band3homogeneity, Band3dissimilarity, Band3entropy, Band3second moment, Band3correlation, and Band6mean contributed more useful information to MRC estimations than the addition of MCRC, NDI5, and NDSVI.The MRC estimation accuracy based on the VI or textural feature methods was improved with PLSR.This suggests that PLSR fully consider the covariance of VIs or textural features and MRC.The combination method obtained higher estimation accuracy of MRC based on the four VIs (NDTI, STI, NDI7, and SRNDI) and three textural feature indicators (Band3mean, Band4mean, and Band5mean) or all VIs and textural feature indicators using PLSR.The results show that the estimation accuracy of MRC was slightly improved based on all the VIs and textural feature indicators.In addition, the four VIs (NDTI, STI, NDI7, and SRNDI) and three textural feature indicators (Band3mean, Band4mean, and Band5mean) were enough to obtain a high estimation accuracy of MRC.The results suggest that PLSR provides a better explanation in the same time of MRC, VIs, and textural feature indicators.Thus, PLSR is a useful tool that can be used to estimate several response variables simultaneously, while accounting for collinear variables.The combination method resulted in a stronger MRC model than that based on only spectral information (Table 5).Berberoglu et al. [48] and Kelsey and Neff [51] indicated that the combination of textural features and spectral information can further enrich the effective information.Thus, the combination method can be used to further improve MRC estimation accuracy (Table 5).In addition, this proved that it is feasible to estimate MRC by combining textural features and spectral information.Finally, the spatial distribution of MRC was mapped using the combination method.The crop residue cover mapping was used to classify tillage management strategies based on the spatial distribution of MRC.In soil carbon models, the spatial distribution of MRC, when implemented within a geographic information system, provides important boundary conditions on the dynamics of soil organic matter across landscapes.Our results provide a method for estimating MRC based on spectral information and textural features from remotely sensed signals in a regional scale.In some countries, the crop residue is spread evenly as the crop is harvested with a combine harvester.In this study, the maize residue cover (MRC) was estimated based on spectral information and textural features from Landsat-8 OLI images.Previous results have indicated differences in the spatial distribution of crop residue even when crops were harvested with a combine harvester, [2,10,14,47].Therefore, this method can be used to estimate the uniform distribution of crop residue on a field scale.It is important to note that this study was limited to the two Landsat-8 OLI images obtained between Changchun in Jilin Province and Suihua in Heilongjiang Province, China.Future studies should obtain more Landsat-8 OLI images to verify the applicability of our results to different crops and environments.

Figure 2 .
Figure 2. Maize residue cover sampling design in this study.

Figure 5 .
Figure 5. Relationships between predicted and measured MRC: (a) VI method; (b) textural feature method; (c) combination method; and (d) combination method based on all VIs and textural feature indicators.

Figure 6 .
Figure 6.Spatial distribution of MRC in the Songnen Plain.

Table 1 .
Descriptive statistics of maize residue cover (%) in the study area.

Table 2 .
Summary of selected satellite optical vegetation indices, wavebands, and references for crop residue cover estimation.

Table 3 .
Relationships between maize residue cover and vegetation indices.

Table 4 .
Relationships between maize residue cover and textural feature indicators.Band2 mean represents the mean texture of Band2; Band3 mean , Band3 homogeneity , Band3 dissimilarity , Band3 entropy , Band3 second moment , and Band3 correlation represent the mean, homogeneity, dissimilarity, entropy, second moment, and correlation texture of Band3, respectively; Band4 mean represents the mean texture of Band4; Band5 mean represents the mean texture of Band5; Band6 mean represents the mean texture of Band6; ** and * indicate a significant relationship at the 0.05 and 0.01 probability levels, respectively.

Table 5 .
Comparison of three methods for estimating of maize residue cover using partial least squares regression.