Monitoring Multiple Cropping Index of Henan Province, China Based on MODIS-EVI Time Series Data and Savitzky-Golay Filtering Algorithm

Multiple cropping index (MCI) is a very important indicator in crop production and agricultural intensification, which represents the utilizing degree of agriculture resources at time scale and the effective utilization situation of arable land. The objective of this paper is monitoring multiple cropping index of Henan province of China according to the time series of MODIS (Moderate-Resolution Imaging Spectroradiometer) EVI (Enhanced Vegetation Index) after Savitzky-Golay filter processing from the year 2006 to 2011.The results revealed that this method could provide an effective way to monitor multiple cropping index, and the method of no additional authentication data is independent and reliable. The result was accurate and stable, the slope of linear regression of the multiple cropping index between the statistical results and the remote sensing results was 1.0136 (R=0.779). The precision of sample areas validation was 97.91%. Suggesting that the time series MODIS-EVI which after Savitzky-Golay filtering processed, could provide an effective way to extract spatial information of multiple cropping index for management department of agriculture.


Introduction
In agriculture, multiple cropping is the practice of growing two or more crops in the same space during a single growing season, which is one of effective ways to ensure region grain security. Multiple cropping is found in many agricultural traditions and prevalent in China. Multiple cropping is measured by multiple cropping index. The multiple cropping index (MCI) is an important indicator of land utilization intensity, and the temporal-spatial dynamics can also help to understand the coupling effects of human activities and the ecological environment [Li, Liu, Sun et al. (2018)]. The MCI is computed as: MCI(%)= annual total sown area of the crops×100/annual total cultivated land area. Verburg et al. Wu (2004)] extracted MCI of china using SPOT/VEGETATION NDVI dataset from the year 1999 to 2002. Based on the method of Fan and Wu, Yan et al. [Yan, Cao, Liu et al. (2005)] added some limiting condition and extracted multiple cropping information from 8 km 10 days composite AVHRR/NDVI. Zhu et al. [Zhu, Li, Shen et al. (2008)] using SPOT-VGTNDVI, the MCI of 17 provinces in northern China form the year 1999 to 2004 was extracted. Zuo et al. [Zuo, Wang, Liu et al. (2013)] explored spatial explicit multiple cropping efficiency of China in 2005 by coupling time series remote sensing data with an econometric model-stochastic frontier analysis. Henan province has well-developed agricultural, which located in the middle and lower reaches of Yellow river basin and between latitudes 31°23′~36°22′N and longitudes 110°21′~116°39′ E in central China (Fig. 1). The annual mean temperature fluctuates between 12 and 16°C and the yearly precipitation varies between 500 and 900 mm, decreasing from south to north. Landscape sub-plains, hilly land and mountainous areas three categories, the plains account for about 55.7% of the total area, mountain and hilly land account for about 44.3% of the total area. Henan is the second largest province of cultivated land area in all provinces of China, the total arable lands is about 8.15×10 4 km 2 , covers 6.24% of the total cultivated land in China. Henan is the most important production base of wheat, cotton and oil plants in China. The grain, cotton and oil plants output account for 1/9, 1/6 and 1/7 national output respectively. It is the leading typical two crops a year region of China, the mainly summer crops are winter wheat and the mainly autumn crops is maize. The total grain planting area of 2010 in Henan reach 1.417×10 5 km 2 , among all of the grain crops in the province, wheat ranks first in planting acreage, accounting for 54% of total arable land in the province, and the yield of wheat of Henan accounting for more than 20% of all china wheat yields, occupy the first in China provinces. Maize is the second largest grain crop, which is the most important autumn grain crop. The crop yield of Henan province directly affects the people's standard of living and social stability of Henan and even the whole country. Consequently, monitoring the spatial and temporal changes of MCI of Henan is particularly significance for the whole nation, which can provide accurate and timely decision support information for Henan agricultural decision-making departments and agricultural enterprises. While the research on Henan MCI spatial distribution and space-time changes rarely involved. The important limitation of monitoring MCI using time series remote sensing data was noise disturbances, especially cloud cover, which were commonly observed in the remote sensing datasets [Qiu, Zhong, Tang et al. (2014)]. The original VI time series datasets with cloud cover and abnormally low values could be identified as troughs. To overcome the scarcity, this paper proposes an efficient methodology, Savitzky-Golay filtering, to reconstruct and smooth the time series data to filter residual noise and provide more reliable data for study. , vol.119, no.2, pp. 331-348, 2019 Figure 1: The location and false-color composited imagery of the study area

CMES
The objective of this study is to extract the MCI and its spatial and temporal changes information of Henan province using 250 m time series MODIS-EVI after Savitzky-Golay smoothing method processed from 2006 to 2011. The results will assist government agencies in decision-making for agricultural policy and administer.

Materials
The data has been collected which contains time series of MODIS-EVI data and Land cover data to extract MCI.

Time series of MODIS-EVI data
Moderate Resolution Imaging Spectroradiometer (MODIS) images have been utilized in recent years because they offer a distinctive capability in maintaining both spatial and temporal density for crop mapping from regional to global scales [Biggs, Thenkabail, Gumma et al. (2006); Pan, Li, Zhang et al. (2012);Singh, Dutta, Stein et al. (2012)]. In particular, MODIS time-series datasets, with high temporal and intermediate spatial resolution, offered a unique ability for crop mapping from regional to global scales [Arvor, Jonathan, Meirelles et al. (2011);Biradar and Xiao (2011)]. The MODIS-EVI time series data used in this study was derived from 16-day composite MODIS-EVI reflectance data at 250m spatial resolution from the MOD13Q1(Terra) and MYD13Q1(Aqua) spanning from 2006 to 2011, ordered through NASA EOSDIS. The MOD13Q1 and MYD13Q1 dataset contains 12 layers, including NDVI, EVI, red reflectance, blue reflectance, NIR reflectance, MIR reflectance, VI quality, view zenith, sun zenith, relative azimuth angle, QA, and 16 days composite day of year. The enhanced vegetation index (EVI) was developed to enhance the vegetation signal by reducing influences from the atmosphere and canopy background and to improve sensitivity in high biomass regions [Huete, Didan, Miura et al. (2002)]. Unlike NDVI, EVI also incorporates a soil adjustment factor as well as an atmosphere resistance term using the blue band in its formulation [Sjöström, Ardö, Arneth et al. (2011)]. Several Previous studies have shown that EVI performs better than NDVI in estimating crop area [Houborg, Soegaard and Boegh (2007); Huete, Didan, Miura et al. (2002); Justice, Townshend, Vermote et al. (2002)]. In addition, EVI is less susceptible than NDVI to biases resulting from cloud and haze contamination [Miura, Huete, Yoshioka et al. (2001); Waring, Coops, Fan et al. (2006)]. The MODIS-EVI data were re-projected from Sinusoidal to Albers Equal Area projection using the MODIS Reprojection Tool. The MODIS data were line stretched to 0～255. The MOD13Q1 and MYD13Q1 products have 8 d time interval, combine MOD13Q1 along with MYD13Q1, the product is similar to the 8 d composite time series data.

Land cover data
Vegetation Index is not a characteristic index value unique to crops. To interpret and distinguish natural vegetation and crops in the case of a large amount of natural vegetation in the research area. The study extracts the cultivated land information of Henan province based on HJ-1A/1B remote sensing images. The HJ-1A/1B images were acquired on 2010 and the images were geometric correction and atmospheric correction. HJ-1A/B is a new generation of small Chinese civilian Earth-observing optical RS satellites. The widecoverage multispectral CCD camera has four bands of blue, green, red and shortwave infrared spectral wavelengths (B1:0.43-0.52 μm, B2:0.52-0.60 μm, B3:0.63-0.69 μm, B4:0.76-0.90 μm). The CCD camera has nadir pixel resolution of 30 m, width of view of 360 km and central-pixel matching accuracy of 0.3 pixels [Wang, Wu, Li et al. (2010)]. Assisted by field observation samples, the multi-temporal HJ-1A/B images covering the whole Henan were classified using eCognition Developer software (Trimble, Inc.), which is the original object-based image analysis software. Chubey et al. [Chubey, Franklin and Wulder (2006)] demonstrated that the object-oriented method could improve classification accuracy better than traditional classification methods. Qi et al. [Qi, Yeh, Li et al. (2012)] demonstrated that the object-based method could reduce speckle effects substantially by combining textural information, getting a satisfactory result with an overall accuracy of 86.6%. The land cover classification process in Henan province includes remote sensing data preprocessing, image multi-scale image segmentation (pre-defined parameters: scale, color/shape and smoothness/compactness), classification features and rule establishment, image classification and the post-processing. The land cover classification results are shown in Fig. 2. In this study, land cover is divided into two types: cultivated land and noncultivated land. The overall accuracy of agricultural land is larger than 95% after humancomputer interaction interpretation, calibrated based on field investigate points. Resampling the land-use data spatial resolution to 250 m, that coincides with the resolution of MODIS-EVI.

Methods
In this study, we used the time series of MODIS-EVI data which was reconstructed and smoothed to extract the MCI of Henan.

MODIS-EVI data smoothing
Due to the effect of sensor, cloud and atmospheric conditions, there are serious residual noise in time series data, which seriously affected the monitoring results [Cihlar, Ly, Li et al. (1997)]. Therefore, prior to application, it is essential to reconstruct and smooth the time series data to filter residual noise and provide more reliable data for study. Numerous algorithms have been developed in the literature for smoothing time series data. The smoothing algorithms are Maximum Value Composite [Holben (1986)], Best Index Slope Extraction Algorithm [Lovel and Graetz (2001)], Asymmetric Gaussian Function Fitting Approach [Per and Lars (2002)], Local Maximum Fitting [Sajia, Ipshita, Kazi et al. (2007)], Fourier Transform [Roerink, Menenti and Verhoef (2000)], Harmonic Analysis Algorithm [Immerzeel, Quiroz and Jong (2005)], Savitzky-Golay filtering [Savitzky and Golay (1964)]. Each approach has advantages and defects, and some aspects need to be improved [Li, Zhang, Liu et al. (2009)]. Considering the advantages and disadvantages of each approach, the Savitzky-Golay was applied to the MODIS-EVI time series data to minimize the effects of cloud cover and other sources of noise [Chen, Jönsson, Tamura et al. (2004)].The Savitzky-Golay is simple and intuitive in theory, which is a low-pass filtering in time-domain method that smooth time series data by using local polynomial regression model. It can be thought of as a generalized moving average. The filter coefficients are derived by performing an unweighted linear least square fit using a polynomial of a given degree. Assume we want to smooth a series of data point's fi, i=0, 1… M. This filter replaces each data value fi by a linear combination gi of itself and some number of nearby Where nL is the number of points used "to the left" of a data point i, while nR is the number of points used "to the right". The idea of Savitzky-Golay filtering is to find filter coefficients cn that approximate the underlying function within the moving window not by a constant (whose estimate is the average), but by a polynomial of higher order, typically quadratic or quartic. The original MODIS-EVI profile and curve smoothing by Savitzky-Golay filtering were shown in Fig. 3. The blue curve shows single cropping (one harvest a year), and the red curve shows double cropping (two plantings a year). As we can see, the MODIS-EVI time spectrum curves become smooth, retain actual details information and variation trend more clearly, ensure the emergence time of feature points (minimum and maximum value) of original curve after S-G filtering smoothing, thereby data more favorable to extract multiple cropping index.

Methodology of multiple cropping index extraction
Researches show the curve of MODIS-EVI time series data is the record of dynamics of crop cultivate, directly presents the physical process of planting, seeding, heading, and harvest of crop in one year [Jonsson and Eklundh (2002); Sakamoto, Yokozawa, Toritani et al. (2005)]. The crest of MODIS-EVI time series data curve corresponding to crop heading stage and trough corresponding to after crop harvest stage [Jiang, Wang, Yang et al. (2002)]. The crest of time series MODIS-EVI profile indicates the colony ground biomass reaches the tiptop, so MCI can be deemed equal to the number of crests of time series MODIS-EVI profile. Consequently, to extract MCI of cultivated land, only need to extract the crest number of time series of MODIS-EVI curve in one year. According to the mathematical implication of crest, the crest must correspond to maximum point on the curve function. Using two times differencing method trace all crests of the time series MODIS-EVI curve. Let EVIi is the EVI value of time phase i on the time series MODIS-EVI curve. Firstly, calculation the difference between adjacent phases MODIS-EVI, which is defined sequence d1. Then the sign of element in sequence d1 need to be determined, if the element sign is negative, and then denoted by -1. If the element sign is positive, denoted by 1, sequence d2 is obtained. Calculation the difference between adjacent elements of sequence d2, which is defined sequence d3. If the element value in sequence d3 is -2, and the before and after elements are 0, the position is crest. If the element value in sequence d3 is 2, and the before and after elements are 0, the position is wave trough. They were calculated as the following equation: After tracing the positions of all crests and wave troughs, calculate the wavelength of each wave and determine each wave meet the characteristics of single crop (the interval between two adjacent crests is at least two months, that is 8 phase) or not. If not, dislodge the interference crest. Then, calculation the curve number of time series EVI curve of each cultivated land pixel in a year, MCI is obtained.

Multiple cropping index of regional
The MCI of regional is computed as the following equation: where P is the MCI of regional, Pi is the MCI of each cultivated land pixel (the number of crest of the pixel time-series vegetation index curve). n is the number of cultivated land pixel in the statistics region.

Multiple cropping index distribution in Henan Province
The spatial distribution of MCI in Henan from 2006 to 2011 is given in Fig. 4. As we can see, the spatial distribution of MCI in Henan has obvious differences. The double-cropped possessed the largest percentage in Henan, with 76% of the entire agriculture land area, which distribution on the north China plain (the east of Henan, i.e., Shangqiu, Zhoukou, Zhumadian) and Nanyang basin. The typical crop rotation system is winter wheat-summer maize (Fig. 5), which is one of the dominant crop planting pattern in Henan and other Northern China. The summer maize is usually sowing in early June and harvested in late September. Winter wheat is sowing in October and harvested in June of the following year. Its curve is show in red curve in Fig. 3. While single-cropped took up 24% of Henan's agriculture land and distribution on hilly and mountain. For instance, in Jiyuan and Sanmenxia (Taihang mountain and FuNiu mountain), part of cultivated land only plant one harvest corn or soybean within a year. In Xinyang (South of Huaihe River), part of the cultivated land plant one harvest middle rice (Fig. 6). The rice is usually seedling emergence in late April and harvested in late September. Triple-cropped accounted for very few (less than 0.1%) agriculture land area, generally for three season vegetables or noise in time series of MODIS-EVI data. The spatial distribution of MCI of each Henan City is shown in Fig. 7 Fig. 4 and Fig. 7 showed that the six-year MCI distribution tendencies were accordant and it in accordance with the system of agricultural regionalization in china, which based on the temperature, moisture, topography and socio-economic factors.

Accuracy assessment
To validate the reliability of monitoring results, this study compared the results to statistical results and use sampling area to verify the results accuracy. According to Henan province statistical yearbook, we computed the MCI of each Henan City and total province from 2006 to 2011. The statistical MCI of cultivated land in total Henan province from the year 2006 to 2011 is 176.57, 177.73, 178.92, 179.11, 179.76, and 179.88 respectively, and the relative error is -1.48%, 0.22%, 0.09%, -2.95%, -3.30%, and -3.34% respectively. The average relative error is -1.79%. The linear relationship correlation between the statistical results and the remote sensing results at city level is shown in Fig. 8. It can be seen from the Fig. 8, the remote sensing results agreed reasonably well with the statistical results, with an r-squared value of 0.779. The slope of linear regression is 1.0136. It demonstrates the remote sensing results and statistical results in good goodness of fit on the city scale. The experimental area is a typical two crops a year region, the mainly summer crop is winter wheat and the mainly autumn crop is maize. The land border and crop type of each field block in May and September are recognized using high resolution remote sensing images as shown in Fig. 9 and Fig. 10. Overlay the summer and autumn crop mapping, the MCI of experimental area is 188.67 in 2008, the remote sensing result is 184.73, and relative accuracy is 97.91%.

Discussion
For a long time, the MCI which used to assess and evaluate the efficiency and sustainability of cropping systems was generally computed from the data collected by traditional survey method that are time consuming and non-spatial. It is urgent to develop a new method extracting MCI independent of statistic data. Nowadays, with the development of remote sensing, lots of vegetation indexes were developed to reflect the growing situation of vegetation. In addition, MODIS data have a sufficiently high temporal resolution. It is a good method to extract MCI from MODIS-EVI time series data. The results of the methodology presented here suggested that it is capable of monitoring MCI across large regions. This method of this study provides an independent and reliable way to make it possible that extraction agricultural information timely, rapid, low cost and high-precision. Instead of directly applying the original MODIS-EVI profile, the methodology proposed in this study overcomes the effect of cloud and atmospheric conditions. However, this approach also has some limitations. Such as the affect of mixed pixels and the error caused by the natural vegetation growth in fallow land. This study selects the MODIS-EVI data with 250 m spatial resolution. On the plains and basin which fragmentation degree of agricultural fields is small, the mixed pixels have almost no obvious effect to the study result. However, for the hilly and mountainous area with cultivated land fragmentation, there are some mixed pixels in 250 m spatial scale remote sensing images, which is one of the important factors that affect the accuracy. The southern region of Henan province (i.e., Xinyang) has better hydrothermal conditions, there may have weeds grew in fallow land. These natural vegetation and crops are indistinguishable in remote sensing images, which may introduce error into the estimation result. In addition, the ground validation of remote sensing results is a complicated task. It needs a lot of ground-based measurements to assess the accuracy of remote sensing monitoring result in different regions, then, identify problems and improve remote sensing algorithm.

Conclusions
The increase of MCI plays an important role in increase food production. Based on analysis and comparison of the various data reconstruction algorithm, this study uses the Savitzky-Golay filter technology smooth the MODIS-EVI time series. Results show that the Savitzky-Golay filter can minimize the effects of cloud cover and other sources of noise. A novel procedure for monitoring the multiple cropping index based on MODIS time series data was proposed in this study. The study extracted MCI of Henan province of China using 8 days composite 250 m time series MODIS-EVI after Savitzky-Golay smoothing filter from 2006 to 2011. Compare the results to statistical results, and use sampling area to verify the results accuracy. The slope of linear regression of the MCI between the statistical results and the remote sensing results was 1.0136 (R 2 =0.779). The total precision of sampling areas was 97.91%. Suggesting that the method of this study is accurate and reliable, and the time series MODIS-EVI could provide an effective way to extract the spatial information of MCI for management department of agriculture. Overall, this method has important practical significance and potential applications, it will be further improved the extraction of MCI more accurate, faster, and provide timely information about the spatial distribution of farming systems to the management department of agriculture. Especially in the macro-scale, high accuracy can be achieved to monitor the spatial distribution of cultivated land cropping index.