Mapping fractional cropland covers in Brazil through integrating LSMA and SDI techniques applied to MODIS imagery

MODIS time-series imagery is promising for generating regional and global land cover products. For Brazil, however, accurate fractional cropland covers (FCC) information is difficult to obtain due to frequent cloud coverage and the mixing-pixel problem. To address these problems, this study developed an innovative approach to mapping the FCC of the Mato Grosso State, Brazil through integrating Linear Spectral Mixture Analysis (LSMA) and Seasonal Dynamic Index (SDI) models. With MOD13Q1 time-series EVI imagery, a SDI was developed to represent the phenology of croplands. Furthermore, fractional land covers (e.g., vegetation, soil, and low albedo components) were derived with the LSMA algorithms. A stepwise regression model was established to estimate the FCC at the regional scale . Finally, ground truth cropland cover information was extracted from Landsat TM imagery using a hybrid method. Results indicated that the combination of multiple feature variables produced better results when compared with individual variables. Through cross-validation and comparative analysis, the coefficient of determination (R 2 ) between the reference and estimated FCCs reached 0.84 with a Root Mean Square Error (RMSE) of 0.13. This indicates that the proposed method effectively improved the accuracy of fractional cropland mapping. When compared to the traditional per-pixel “hard” classification, the sub-pixel level maps illustrated detailed cropland spatial distribution patterns.


Introduction 
Croplands play an essential role in the process of land use and cover changes at both regional and global scales [1,2] .During the past decades, the geographic areas of croplands have increased sharply thanks to population growth and technological improvements.
Information of cropland area and spatial distribution is essential for land use analysis, crop yield estimation, soil and water conservation, farmland protection, and agricultural planning [3][4][5] .For cropland detection and mapping at small scales, medium-and high-spatial resolution satellite remote sensing imagery (e.g., Landsat TM (Thematic Mapper), ETM+ (Enhanced Thematic Mapper P lus), and SPOT HRV (High Resolution Visible) have been successfully employed.However, due to their low temporal resolutions, it is difficult to timely obtain large scale cloud-free images during the crop growing season.Due to this, researchers have resorted to coarse-resolution satellite imagery with high temporal resolutions (such as AVHRR (Advanced Very High Resolution Radiometer), MODIS (Moderate-resolution Imaging Spectro-radiometer), etc.) for region-scale cropland detection [7] .MODIS data, in particular, have great potential in regional and global land cover surveys thanks to its high-temporal and moderate-spatial resolutions [6,[8][9][10] .Continuous time-series MODIS data have been widely applied for regional scale agriculture landscape surveys and subsequent applications [10][11][12][13] .A summary of cropland mapping methods based on MODIS imagery is listed in Table 1.
These methods can be generally divided into tw o categories: pixel level "hard" classification and subpixel level "soft" classification.
It has been noticed that previous methods primarily relied on pixel level "hard" classifications.These methods are able to identify cropland qualitatively; however , they have the disadvantage of not being able to extract cropland fractions accurately and quantitatively because the mixed pixel problem is unavoidable in coarse spatial resolution imagery (e.g., MODIS), resulting in poor area estimation and inaccurate s patial patterns [26] .Therefore, sub-pixel level "soft" classification and fractional cropland information extraction has attracted more attention [22][23][24][25][26] .
In the Brazilian Amazon, many studies have illustrated the value of using MODIS time series imagery for mapping cropland change and expansion at the regional scale [11,16,19,25,[27][28][29] .Although these studies have confirmed the importance of using MODIS data for cropland mapping, how to extract fractional cropland cover (FCC) remains a great challenge due to mixed spectral properties and other problems such as coarse resolution, cloud contamination, pure pixels selection and etc.Although linear spectral mixture analysis technologies present an effective way to decompose the spectral reflectance of a pixel into different components, which has proven valuable in medium spatial resolution images such as Landsat or MODIS NDVI (Normalized Difference Vegetation Index) time series data [22,30,31] , the application of this approach is limited in this area.These difficulties can be summarized as follows: 1) In the Brazilian rainforest areas, the crop growth period occurs in the rainy season.It is difficult to obtain time series data without cloud cover during the growing season, resulting in the difficulty of composing a continuous time series enhanced vegetation index (EVI).
2) With the conventional pixel-based ''hard'' classification methods, each pixel is assigned into one category.The mixed pixels problem leads to high uncertainty for cropland covers with small parcel sizes, resulting in information loss of spatial patterns.
3) Unfixed crop sowing time and vertical intensification (e.g., double cropping, triple cropping) have increased the difficulty of automatic cropland extraction.
As the crop growth period occurs in the rainy season, it is difficult to obtain time series EVI data without cloud contamination.
Therefore, the LSMA (Linear Spectral Mixture Analysis) technique is limited due to the difficulty of deriving time series EVI pure pixels.To address this issue, we developed the seasonal dynamic index (SDI) model to construct a relationship between the FCC and crop phenomenon fluctuations of EVI.However, due to the inconsistency of crop sowing and harvest time, the SDI model has some uncertainty in different crop sowing time without considering fallow land.
Therefore, from this perspective, the LSMA components in the dry season are effective supplementary information for the FCC estimation.For 500 m×500 m coarse spatial resolution imagery, the traditional "hard" classification results in the loss of cropland information (e.g., underestimation).Therefore, considering the heterogeneous environmental conditions, agricultural phenology cycles, and mixed pixels problems, this paper proposed a new approach to estimate FCC by integrating multiple feature variables including components developed using the LSMA and the newly developed SDI model.

Spectral similarity
Judges the similarity of the 2 curves by the WCD (Dynamic Time Warping) and others algorithm.
Traditional classification methods and accuracy controllable.
Clear physical meaning and being able to estimate fractional distribution.
Hard to find a proper endmember in larger scale.[22-25]   2 Study area and materials

Study area
The study area is located in the state of Mato Grosso, Brazil (see Figure 1) (longitude 50°-65°W and latitude 10°-20°S).This area is the main soybean producing region in Brazil [19,32] .Depending on the region and the onset of the rainy season, the sowing calendar for soybeans goes from mid-September to late December [19,32] .The conversion rate of wetlands, pasture, and forests to cropland has increased sharply during past few decades [29,33]   , because of the intensive agriculture.According to IBGE 2015 (The Brazilian Institute of Geography and Statistics), the soybean planting area has increased by 5.59 million hm 2 from 1995 (2.34 million hm 2 ) to 2013 (7.93 million hm 2 ).Meanwhile, intensification practices such as double cropping have been widely adopted in this state [20] .

Data collection and preprocessing
MODIS/Terra Vegetation Indices 16-Day L3 Global 250 m SIN Grid (MOD13Q1) for the period 2011-2012 was acquired for this study (h12v10).The MOD13Q1 product provides two vegetation Index (VI) layers including the Normalized Difference Vegetation Index (NDVI) and the Enhanced Vegetation Index (EVI), the two quality assess (QA) layers and surface reflectance bands 1 (Red), 2 (NIR), 3 (Blue), and 7 (MIR), as well as four observation layers.The data composite criteria are low clouds, low view angle and the highest EVI value from all acquisitions between the 16 d periods.All the images were obtained from the United States Geological Survey (USGS) website for free [34] .Because NDVI is easily saturated in high biomass regions and EVI is more sensitive to dense vegetation conditions [19] .This study utilized the enhanced vegetation index (EVI), which is less affected by atmosphere and soil and more suitable for cropland extraction [35 ] .But, the EVI uses the blue band to removal residual atmosphere contamination caused by smoke and thin cloud.Lacking a 250 m blue band, the EVI algorithm uses the 500 m blue band to correct for residual atmospheric effects.The MODIS EVI is defined in Equation (1) [19,35] : where, R, NIR, and B is the red, near infrared, and blue bands, respectively; G = 2.5 (the gains factor); and L = 1, C 1 = 6, and C 2 = 7.5 are the adjusting parameters used to minimize aerosol effects.
The entire MODIS EVI time-series, comprised of images every eight days from 2011 to 2012 were filtered using the Savitzky-Golay (S-G) [19,36,37] filter to remove noise and artifacts caused by thin clouds.

Shuttle Radar Topographic Mission (SRTM)
The SRTM DEM (digital elevation model) data with a spatial resolution of 90 m were downloaded from the website of USGS (United States Geological Survey) [38] .The DEM data will be further employed for calculating the SDI index.

Landsat image
The Landsat TM imagery (path/row: 227/68 and 228/69, acquired in 2011/02/24) was downloaded from USGS Global Visualization Viewer [34] .Data preprocessing and classification proceeded with ENVI software.Landsat TM classification data were used as the reference data source for regression and validation cropland.The cropland imagery with 30 m spatial resolution was aggregated into a new image with a cell size of 250 m for generating fractional cropland data to match the cell size of MODIS EVI data using the mean algorithm.
All datasets where projected to the Lambert Azimuthal Equal Area Projection in order to standardize the cropland area estimates.

Methodology
The framework of the proposed integrative method for FCC estimation from MODIS and Landsat TM imagery is illustrated in

Key identification stage and seasonal dynamic index
The seasonal dynamic index (SDI) assumes that the variation value of EVI is positively related to the proportion of cropland area in a pixel.For Mato Grosso, six of the largest crop type classes (soy-corn, soy-cotton, soy-millet, soy-soy, cotton, and pasture) account for 91.5% of the reported agricultural land area in Mato Grosso.The EVI profiles in the cropland area vary drastically with regularly changing characteristics in a year due to the sowing, growing, and harvest stages in a year (single cropping and double cropping).In contrast, the EVI profiles of forest have almost no changes throughout the year at all and the EVI profiles of Grassland changes are relatively small.In summary, the EVI values fluctuate significantly for croplands, and are relatively stable for grassland and forest areas.
The key identification stage selection is fundamental for constructing a seasonal dynamic index, which is critical for cropland mapping in the Brazilian Amazon region.As almost all images are contaminated by clouds in the rainy season, a feasible means is to use slices of discrete time series data instead of entire continuous time series data for a year for crop mapping [26] .As the EVI profiles in cropland area vary regularly at different stages (e.g., the sowing, growing, and harvest stage), three key identification stages: the sowing (Stage 1, DOY:225-289), growing (Stage 2, DOY:305-001), and harvest (Stage 3, DOY:017-081) seasons are the key identification stages for cropland mapping [26] .The seasonal dynamic index (SDI) model was proposed by Zhu et al..The model can be elaborated by Equations ( 2)-( 8) [26] : Pas mask is a pasture mask.

Linear spectral mixture analysis
Linear spectral mixture analysis (LSMA) assumes that any spectrum measured by the sensor is a linear combination of the spectra of pure uniform targets called endmembers [31,39,40] .It presents a method to decompose the spectral reflectance of a pixel into different proportions [22,31,[41][42][43][44] .Within a pixel, the spectral proportion of endmembers represents a fraction of the area covered by distinct features on the ground, as shown in Equation (9).
where, ρ is the sensor acquired pixel reflectance; A i is the i-th component (distinct feature or endmember) cover area; ρ i is the i-th component reflectance; and ε is the error of residual.
Among the pixel decomposition algorithms, the V-I-S (vegetation, impervious, and soil) conceptua l model is considered to be a classic method.It provides a guideline for decomposing low-resolution images of urban landscapes and linking these components to spectral signatures [45,46] .Inspired by the V-I-S model, a V-S-L (vegetation, soil, and low albedo objects) model was explored and developed in this paper.In the V-S-L model, a constrained least-squares algorithm was applied to decompose six reflective bands into three fractional images and one error distribution image.The endmembers (vegetation, soil, and other low albedo) were selected from the scatterplot in spectral space [47] .
where, ρ is pixel reflectance; C v , C s , and C 0 are components of vegetation, soil, and other low albedo objects, respectively; ρ v is the endmember of vegetation; ρ s is the endmember of soil; ρ 0 is the endmember of low albedo objects , and ε is the error of residuals.
In this paper, the V-S-L was used for unmixing MODIS surface reflectance in dry season.

Cropland classification from Landsat TM imagery
The TM derived FCC was a standard dataset for regression and validation.In this paper, an unsupervised classification algorithm (i.e., ISODATA) was used.The Landsat TM image was first classified into 50 clusters.An analyst assigned each cluster into cropland or non-cropland through visual interpretation on the TM nature color composites images.Therefore, the initial classified image was a binary thematic map with 1 representing cropland area and 0 representing non-cropland area.The binary image was then aggregated to FCC with a spatial resolution of 250 m×250 m by the means algorithm.

Fractional cropland mapping and accuracy validation
The SDI and LSMA components were employed as multiple independent variables.The TM derived FCC was used as a dependent variable.
One thousand samples were randomly selected between the dependent variable and the linked independent variables.Half of them were employed for developing regression models and the rest for validation.
The coefficient of determination, R 2 , was applied to measure the percentage of variances explained by the regression model.The F test was employed to examine whether the regression model was significant or not, and the t test was utilized to examine whether the constant and beta values were significant or not.The highest R 2 value and the significant F and t tests were selected for a further established regression equation.Finally, the regression equation was applied to estimate the FCC in the whole area.In this study, regression models were significant based on the F test at the 95% confidence level.

FCC distribution mapping
Table 2 presents the all regression models with SDI and the LSMA components.It shows that the two single feature linear regression models had similar RMSE and R 2 , where the RMSE was 0.15 and 0.14, and the R 2 was 0.78 and 0.76, respectively.From the regression models listed here, it seems that the SDI based linear model was slightly better than the LSMA based, but the result was not significant.However, the multiple features regression model provided a higher R 2 (0.89) and lower RMSE (0.13) than any of the single variable models.The model accuracy was improved greatly by integrating multiple feature variables.Figure 3 presents the comparative analysis between the reference and estimated data at the two sites.In these two areas, the TM derived FCC data were used for comparison with the estimated FCC by MODIS data.Figures 3(1)-(3) show the comparison of the estimated FCC with the TM derived FCC at validation Site 1. Figures 3(4)-( 6) are the comparison of the estimated FCC with the TM derived FCC at validation Site 2. From Figure 3, the estimation results showed a good performance of the proposed integrative model.The FCC estimated by the integrated method were highly similar to the TM derived FCC on the whole.

Accuracy assessment
Figure 4 shows the residual errors distribution of estimated FCC by SDI and LSMA. Figure 4a is the residual errors of the estimated FCC at Site 1. Figure 4b is the residual errors of the estimated FCC at Site 2. From the residual errors maps, errors were randomly distributed and concentrated on the zero line.The residual plots mainly fluctuated from -0.2 to 0.2.That implies that the model can externalize the relationship between the explanatory variables and the dependent variable.The cropland distribution information can be explained by the estimated FCC rationally.2D scatter plots are a conventional method for testing the effectiveness of the estimate algorithm.Figure 5 shows the scatterplots between the reference and estimate with different methods.In Figure 5, the X axis is the estimated FCC, and the Y axis is the reference data.For each site, 100 samples were selected randomly.
Figures 5a-5c show the scatterplots comparison among the SDI, LSMA, and the hybrid method for Site 1. Figures 5d-5f are the scatterplot comparisons among the different models for Site 2. From Figure 5, we found that the estimated FCC products were in good agreement with the reference data overall by different methods.The integration of SDI and LSMA, however, had a much higher coefficient of determination (R 2 =0.84).This implies that more than 80% of the variances could be explained by the estimation model.Larger estimated errors occurred with an FCC less than 0.2 or more than 0.8.For an FCC lower than 0.2, the SDI model slightly overestimated the FCC as a number of points were located on the X-axis.Similarly, the SDI model slightly underestimated the FCC values when they were close to 1.0.This may be a defect of the model in the estimated FCC on a coarse resolution image at the sub-pixel level.When a mixed pixel was less than 20% of the cropland, the dominant property was non-cropland.
This was easily underestimated as both the elemental spectrum and the phenology had no significant difference with noncropland.Similarly, a pixel component of more than 80%-90% was cropland that was easily overestimated as it was indistinguishable from pure pixels (100%) in both the spectrum and the phenology.

FCC mapping from 2002-2012
More tests were employed to explore the transferability and feasibility of this proposed approach.Time series MODIS images from 2002 to 2012 at two-year intervals were collected.We applied this approach to other years of MODIS data in the same study area for estimating the FCC and monitoring cropland dynamic changes.Figure 6 shows the gradient of FCC maps between 2002 and 2012.
To better illustrate the spatial heterogeneity of croplands on the coarse spatial resolution image, we divided the fraction image into four grades using the threshold segmentation technique based on the scatterplots analysis.In particular, a pixel was assigned to non-cropland with an FCC lower than 0.2, to low fraction cropland with an FCC between 0.2 and 0.4, to moderate fraction cropland with an FCC between 0.4 and 0.6, and to high fraction cropland with an FCC larger than 0.6.On the FCC maps, croplands were largely distributed in the southern (along the BR-163 road) and eastern (Parecis plateau) regions of the study area, where many large farms are distributed.In the northern and western areas, croplands are relatively rare and their patches were small with an FCC less than 50%.Most of these areas are newly reclaimed regions.It was also observed that the main agricultural areas had higher cropland fractions than those in new colonization areas.By the area statistics, cropland areas changed significantly between 2002 and 2012.Specifically, the geographical area of high fraction croplands increased from 0.2 million hm 2 in 2002 to 1.18 million hm 2 in 2012.Similarly, the low fraction cropland area rose from 2.32 million hm 2 to 3.94 million hm 2 , and the moderate cropland area expanded from 1.09 million hm 2 to 1.83 million hm 2 .Consequently, the no cropland area decreased from 2.88 million hm 2 to 2.54 million hm 2 .This reflected that Mato Grosso had been undergoing rapid cropland expansion during the process of Brazil's agricultural intensification.

Discussion
The traditional per-pixel "hard" classification method has not considered the mixture pixel problem.On a coarse resolution remote sensing image though, mixing pixels is a common phenomenon.Arbitrarily classified land covers into cropland or non-cropland are inappropriate and result in small cropland parcel information loss.Therefore, it is not suitable for accurate cropland area statistics.The potential errors associated with the classification of mixed pixels have been widely recognized as a problem that affects the accuracy in image classification [48][49][50] .Therefore, cropland mapping at the sub-pixel level with "soft" classification is a promising method in the future [19,22,24,26,51] .More cropland classification research has focused on fractional mapping using MODIS datasets [22][23][24][25][26] .
Estimated errors mainly come from the following aspects.(1) Unfixed crops calendar.The sowing calendar for crops goes from mid-September to late December, depending on agricultural zoning for different soils, regions, and the onset of the rainy season [19,32] .At the same time, vertical intensification such as double cropping has been widely adopted in Mato Grosso, which is an additional challenge for accurate cropland area mapping [20] .(2) Poor image quality.Even though the SDI model used the time spare resampling algorithm to composite multi-temporal data to one period of data, the inconsistency of crop time information increased the estimation errors of the model [26] .(3) Geometric errors.The difference of spatial resolution between MODIS and Landsat imagery is more than the size of two Landsat pixels (60 m).The mis-registration between Landsat and MODIS can reach a minimum of 50 m (NADIR) [52,53] .In the middle of uniform regions, the cropland proportion tend to be maintained in the same level and similar to the estimated proportion by the regressions, while the borders from one side tend to be lower proportion and the other side of the possible mis-registered map side higher proportion.
Validity and sensitivity of the model, the number of the negative values of residuals is small when FCC is in the range of 0 to 0.2.The number of the positive values of residuals is small when FCC is in the range of 0.8 to 1.0.These residuals can be explained by the imaging mechanism.As we all know, mixed pixels are a common phenomenon on coarse resolution MODIS imagery.
Furthermore, the spectroscopy and phenology characteristic of a pixel was determined by the dominant land cover.Thus, the overestimation and underestimation of the model by MODIS data was inevitable.A pixel component more than 80% was cropland that was easily overestimated as its spectral property was indistinguishable from the pure pixels.Similarly, if a pixel was less than 20% of the cropland, it was easily underestimated because both the elemental spectrum and the phenology characteristics were dominated by non-cropland.In the subpixel, 20% is a threshold, as less or bigger than the threshold cannot be significantly manifested in a mixed pixel.So far, there are no better ways to improve the model's sensitivity as the method relies heavily on the time series, a suitable task for MODIS type data.
Moreover, with the constrained LSMA model, it is assumed that the summation of the fraction of each composition (e.g., vegetation, soil, and low-albedo objects) equals to one (sum-to-one), and each of them is nonnegative (non-negativity).Therefore, although the constrained LSMA model has been widely applied due to its physical soundness, the model overestimated the fractions when the FCC was close to zero, and overestimated the fractions when FCC approached one.That is, negative residuals existed with low FCC values, and positive residuals were dominant with high FCC values.Although this problem may obviously affect the estimation accuracy in low and high FCC areas, it is very difficult to address.Chang et al. [54] proposed an unconstrained LSMA model to address this issue, and pointed out that the unconstrained LSMA was better for object identification, and the constrained LSMA was ideal for quantifying land cover fractions.Recently, Wang et al. [55] integrated LSMA and classification and regression tree (CART) to address this issue.The developed method, however, is ad hoc, and is difficult to apply to such a large study area.Therefore, although with the bias of estimation with low and high FCC, the constrained LSMA was still employed in this study, and the estimation results were still acceptable with an RMSE of 0.13.Moreover, the multivariate regression analysis did help to reduce the bias.
However, the combination method just used a linear regression model.The algorithm structure needs be further improved and enhanced.With the development of artificial intelligence, a large number of intelligent nonlinear algorithms have been proposed, i.e., ANN (Artificial Neural Network), SVM(Support Vector Machine), Deep Learning, and Deep convolutional neural network [11,[19][20][21][56][57][58][59][60] . These nnlinear intelligent algorithms have a better performance in fitting multiple sets of variables.Furthermore, it needs to explore and mine more feature variables expression.This new feature can be better expressed in decimals (the range of value 0.8+ and 0.2-) sensitively.
Nevertheless, this research has shown that the integrative use of feature variables with SDI and LSMA can successfully estimate fractional cropland from MODIS imagery.

Conclusions
This research developed a seasonal dynamic index (SDI) using phenological information, as well as a linear spectral mixture analysis technique, and further applied a multivariate regression analysis to estimate fractional cropland coverage (FCC).Because the SDI and LSMA have a different theoretical basis and hypothesis, these two methods can complement each other for generating better FCC products.
Experiments also show integrating LSMA and SDI technology effectively improved the accuracy of FCC mapping.Comparing the estimation results with the reference data, the coefficient of determination (R 2 ) rose from 0.76 to 0.84 with a much lower RMSE value.Therefore, we can come to the conclusion that a combination of SDI and LSMA provided a better estimation performance than any of the individual feature models.This method may serve as a better alternative for regional FCC mapping for Brazil due to the weather condition.
of China (Granted No. 2017YFB0504201) and the Natural Science Foundation of China (Grant No. 61473286 and 41201460). [References]

2 Figure 1
Figure 1 Study area and the two validation sites

Figure 2 .
The entire process consists of three major steps: (1) generating the SDI indices to represent key growth stages of croplands; (2) calculating the fractions of vegetation, soil, and low albedo (V-S-L) components from the MODIS dry season image through employing LSMA; and (3) estimating FCCs by constructing regression equations of FCCs and SDI, vegetation fraction, soil fraction, and low albedo fraction, and validating estimation accuracy.

Figure 2
Figure 2 Framework of FCC (Fractional Cropland Covers) estimation from MODIS and Landsat data

Figure 3 Figure 4
Figure 3 Local comparative analysis between the estimate of FCC and TM derived reference data

Figure 5
Figure 5 Comparison of the FCC from MODIS estimated and TM derived among different methods

Table 1 Summary of major algorithms for cropland extraction by MODIS data
1 , SD 2 ) × Mask 2 correspond to the seasonal dynamic index at different stages; EVI d , EVI g , and EVI h are cloud-free EVI composites from the dry to wet season transition, the growth, and the harvest season, respectively; EVI 225 , EVI 241 , …, EVI 353 , are the multi-temporal MODIS EVI products, and the number subscript is the acquired day of the year (DOY); Slp mask is the topographic factor mask where the slope is derived from STRM data; Slp mask is a slope mask;

Table 2 Regression models developed from the combination of different feature variables
Fa is fractional cropland, Cs is the soil component of LSMA, R 2 represents the coefficient of determination for the evaluation of the regression model performance.