# Is Time Series Smoothing Function Necessary for Crop Mapping? - Evidence from Spectral Angle Mapper After Empirical Analysis

Ailian Chen, Hu Zhao, Zhiyuan Pei

► To cite this version:

Ailian Chen, Hu Zhao, Zhiyuan Pei. Is Time Series Smoothing Function Necessary for Crop Mapping? - Evidence from Spectral Angle Mapper After Empirical Analysis. 9th International Conference on Computer and Computing Technologies in Agriculture (CCTA), Sep 2015, Beijing, China. pp.335-347, 10.1007/978-3-319-48357-3_33. hal-01557823

HAL Id: hal-01557823

https://inria.hal.science/hal-01557823

Submitted on 6 Jul 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Is time series smoothing function necessary for crop mapping? —Evidence from spectral angle mapper after empirical analysis

Ailian Chen ¹,*, Hu Zhao ¹,ᵇ, Zhiyuan Pei ¹,ᶜ

¹Chinese Academy of Agricultural Engineering, Beijing, 100125, PR China;
*cal-0601@163.com, ᵇzhhoo@126.com, ᶜPeizhiyuan@hotmail.com

**Abstract.** Time series smoothing functions have been frequently applied to fit multi-temporal vegetation index for better extraction of plant seasonal/growing parameters. Questions are raised that whether the smoothing is necessary for crop mapping. Four time series smoothing functions, namely, HANTS, Savitzky-Golay (S-G), double logistics and asymmetric Gaussian, were used to smooth 23 MODIS 16-days composite NDVI images in one year. The effectiveness were compared through visual check, correlation coefficient R, root mean square error (RMSE), and local signal noise ratio (SNR). The best smoothing time series NDVI images, along with the original time series images, were then used to map corn and soybeans by spectral angle mapper (SAM) method and their mapping accuracies were compared. Comparison of smoothing results showed that S-G fitted data got the strongest correlation coefficient R, the lowest RMSE and lower local SNR. Comparison of mapping results further showed that time smoothing function does not improve the classification accuracy obviously with the same training sample and same temporal bands. The whole analysis indicates that it is the mapping method that matters more than time series smoothing function for classification precision.

**Keywords:** MODIS, NDVI, crop phenology, TIMESAT, SAM

## 1 Introduction

High temporal resolution time series (TS) vegetation indices (VIs) have been well used in the research of global environment with respect to vegetation phenology. Large scale environmental parameters, such as global land use/cover change, could be derived with the knowledge of vegetation phenology. However, the high temporal resolution time series VIs often have coarse spatial resolution, which often make the time series VIs noisy and fluctuating due to the complex atmospheric condition within one coarse pixel (such as undetected cloud) and varying sun-sensor-surface geometries (Duggin, 1985). Therefore, a number of noise-reduction functions, including maximum value composite (Huete et al. 2002), Fourier transform (Jentsch and Subba Rao 2014), wavelet transform (Sakamoto et al. 2010), Gaussian transform (Hird and McDermid 2009), etc., were proposed or utilized to solve this problem.

*Corresponding author: Chen Ailian, PhD in Ecology
As for crop mapping, some studies (Sakamoto et al. 2010, Brown et al. 2013, Zhang et al. 2014) have used TS fitting function before mapping, while many others have not (Chang et al. 2007, Wardlow and Egbert 2008, Gusso et al. 2014, de Souza et al. 2015). Some of those which have used time series smoothing functions often conducted more than one smoothing function to the time series images. For example, researchers used MODIS Level 3 VI products always conducted smoothing function upon the already composited VIs by maximum value composite (MVC) (Huete et al. 2002, Daniel de Castro Victoria, et al, 2012).

A scanty few studies have been conducted regarding the benefit of noise reduction for vegetation phenology (Hird and McDermid 2009), or regarding empirical information of noise of the smoothed images themselves, such as Atzberger and Eilers (2011), who have compared the effectiveness of four fitting methods. Fewer studies have investigated the necessity of fitting TS VIs for crop classification. Questions are raised that whether smoothing functions are necessary for crop type mapping. This goals of this paper are: (1) to study the effectiveness of four different time series fitting methods, and then with the best fitting results, (2) to compare the classification accuracy.

# 2 Materials and Methodology

## 2.1 Materials

MODIS 16 days composite NDVI of the year 2013 were used for time series analysis. The time series started at 2013001 as band 1, and ends at 2013353 as band 23, with an interval of 16 days. The value in band 1 represents the max NDVI from the date Jan. 1st to the date Jan. 16th, 2013, and so on. The MODIS NDVI 16 days composite NDVI named "MOD13Q1-Level 3 16-Day Vegetation Indices-250m", is one of the MODIS 16-Day Tiled Land Products (collection 5). The data was downloaded from the U.S. National Aeronautics and Space Administration (NASA). The "MOD13Q1-Level 3 16-Day Vegetation Indices-250m" included both NDVI and enhanced VI, EVI. The NDVI was extracted for research as it was a normalized value. MODIS data are distributed from NASA on a near real-time basis at no cost. It is available at two satellite platforms, Terra and Aqua, which were launched in December 1999 and May 2002, respectively. The spatial resolutions of the MODIS-250 m (MOD13Q1, MYD13Q1) on the sinusoidal projections are 231.7 m. The downloaded tiles were mosaicked, re-projected to Albers projection and resampled to 250m.

Crop data layer (CDL) was used for training and accuracy check. CDL is a specific crop type classification data with resampled resolution of 30m, and is produced and distributed in the name of NASS (national agricultural statistics service), USDA. The data is available through the web of USDA.

One U.S. state, the Ohio, was taken as case study to classify corn and soybeans from others, as in large area mapping, the large area was often divided into small regions. Corn and soybeans are the dominant crop in Ohio, where the yield of corn takes up 4% of the whole state total, the yield of soybeans takes up 7%. Other parts of the
state are mainly forest and shrub lands. The spatial distribution of corn and soybeans
from CDL 2013 was shown in Fig. 1.

<figure><img src="image_2.png" /><figcaption>Fig.1.Corn and soybeans distribution in 30m resolution crop data layer (CDL) data in 2013</figcaption></figure>

2.2 Time series smoothing function and result check

Four popular (Hird and McDermid 2009, Zhang and Ren 2014) smoothing functions
were conducted, including Harmonic analysis of TS (HANTS) (Zhou et al. 2015), S-
G filtering function short for adaptive Savitzky-Golay filtering (Chen et al. 2004),
asymmetric Gaussians (referred to as Gauss herein after) (Jönsson and Eklundh 200)
and double logistic functions (referred to as Logistics herein after) (Beck et al. 2006).
The basic idea of HANTS is to calculate a Fourier series to model a TS of pixel-wise
data, while meanwhile recognizing abnormal data value relative to the TS of Fourier
model. HANTS deletes abnormal data value and substitutes it with the value given by
the model. It is developed with the interactive data language (IDL). The other three
smoothing algorithms are embedded in the free software TIMESAT (Jönsson and
Eklundh 2004).

Visual and statistical check were conducted for the comparison of the four smoothing
functions. Visual check was done first by the screenshot of the original images and all
the smoothed images, and then by compared the time series NDVI profiles of two
pixels, one of which was the corn pixel, and the other was soybeans.

Statistical check included the comparison of correlation coefficient R of one single
band image (Band 14) and 20 sample points, local SNR and RMSE in Band 14, as
Band 14 represents the NDVI value in late June and early July, when the NDVI corn
and soybeans were near the max value(Wardlow and Egbert 2008) .
Correlation coefficient R (whole image and sample point) was calculated with Eq.1

$$R = \frac{\sum_{i=1}^{N} (x_i - \bar{x}) (y_i - \bar{y})}{\sqrt{\sum_{i=1}^{N} (x_i - \bar{x})^2 (y_i - \bar{y})^2}} \quad (1)$$

Where, $x_i, y_i$ represent the pixel value of the whole non-smoothed image and the smoothed image, and $\bar{x}$ and $\bar{y}$ are mean value of the whole non-smoothed image and the smoothed image in Band 14, respectively. For sample points, $x_i, y_i$ represent the pixel value of samples from the non-smoothed image and the smoothed image in Band 14, and $\bar{x}$ and $\bar{y}$ are the mean value of samples from the non-smoothed image and the smoothed image in Band 14 respectively.

As global SNR requires a homogeneous underlying surface, which is always not easy to find, the local SNR was used, which calculated local SNR using the mode of several grids. The local SNR is also used for the determination of HJ-images (Zhu et al. 2010). Sixteen fishnet grid with edge of 10km was used to calculate SNR, and the mode, mean, min and max values were given. Root mean square error (RMSE), local signal noise ratio SNR (Zhu et al. 2010) of 16 fishnet grid of the whole time series images were calculated through the following equations (Eq.2 to Eq.6

RMSE was calculated with Eq.2

$$RMSE = \sqrt{\frac{\sum_{i=1}^{N} (x_i - y_i)^2}{N}} \quad (2),$$

Local SNR is calculated by Eq.3.

$$SNR = \bar{X}/\sigma \quad (3),$$

where, $\bar{X}$ is the mean value of the grid, and $\sigma$ is the standard error, which are calculated by Eq.4 and Eq.5.

$$\bar{X} = \frac{\sum_{i=1}^{N} x_i}{N} \quad (4),$$

$$\sigma = \sqrt{\frac{\sum_{i=1}^{N} (x_i - \bar{x})^2}{N-1}} \quad (5),$$

where, $x_i$ represents the pixel value of $i$ in the analyzing grid, and $N$ is the total pixel numbers in the grid.

For the calculated R, RMSE, and SNR, the criteria for effectiveness is higher R, lower RMSE and lower SNR.

## 2.3 Crop mapping and result check

Spectral Angle Mapping (SAM) in the commercial software ENVI™ was used for classification. SAM “uses an n dimension angle to identify pixels to reference spectra” (Kruse et al. 1993). The spectral angle is reversed by Eq.6.

$$\cos \theta = \frac{XY}{|X||Y|} = \frac{\sum_{i=1}^{N} x_i y_i}{\sqrt{\sum_{i=1}^{N} x_i^2} \sqrt{\sum_{i=1}^{N} y_i^2}} \quad \theta \in (0, \pi) \quad (6),$$
where, n is the total band number. It is 23 in this study, xi and yi are the data value (usually reflectance) of each time series images, and θ is the spectral angle which the mapper would be based on. The algorithm determines the spectral likeness between two spectrum profiles by figuring out the angle between the spectrum profiles in a space with dimensionality equaling to the number of bands. The SAM method is relatively less sensitive to illumination and albedo effects than other methods such as maximum likelihood.

End-member spectra (training data spectra) used in SAM mapping procedure was collected with the aid of CDL. The CDL was resampled to a resolution of 250*250 m, the same as the MODIS NDVI data, with each 250*250 m pixel assigned a value of percent of corn or soybeans. The pure corn and soybeans were extracted for random sampling. Ten percent of Random samples (points) of corn and soybeans were extracted from the pure corn or soybeans layer, respectively. One third of these random samples were used for training of SAM, which was done by performing a MOD function, going by like MOD(ID,3) equaling to 0, and the left were used for accuracy check. The tolerance angle for SAM was set as 0.5 radian.

For the accuracy check of mapping result, three groups were compared. Overall, producer and user accuracy were all derived for comparison. The first group classified all the 23 original and S-G fitted time series images, the second group classified the 12 images from Band 7-18, and the third group classified 7 images including band 9, and Band12-17. Percent and pixel accuracy of corn and soybeans were given, as well as the overall accuracy. The mapping results with best overall accuracy was compared by ratio to the acreage estimated by CDL and acreage from national statistics released by NASS, USDA.

# 3 Results and Discussion

## 3.1 Time series fitting results

### 3.1.1 Visual check

The screen shots of one sub region of the original image and four smoothed images in Band 14 were shown in Fig. 2. The result of HANTS was different from others to the most extent, which has brought more noise to the original image.
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Screenshot of one random pixel in Band 14 before and after fitting</figcaption></figure>

Time series NDVI profiles of one corn pixel and one soybeans pixel from the original images and four smoothed images were shown in Fig.3 and Fig.4, which indicated that HANTS tended to change the max value position, and obscured the start and end of season date. S-G has well kept the original curve.
<figure><img src="image_4.png" /><figcaption>Fig. 3. NDVI Time series profiles of corn pixel before and after smoothing</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. NDVI Time series profiles of soybeans pixel before and after smoothing</figcaption></figure>

### 3.1.2 Correlation coefficient, RMSE and SNR

Correlation coefficient R of the whole image and the 20 sample points in Band 14 between the smoothed and the non-smoothed images were given in Table 1.
<table><caption>**Table 1.** Correlation coefficient between each smoothed and the non-smoothed images</caption><thead><tr><th></th><th>Original whole image</th><th>Original 20 sample points</th></tr></thead><tbody><tr><td>S-G</td><td>0.996</td><td>0.846</td></tr><tr><td>Logistics</td><td>0.995</td><td>0.751</td></tr><tr><td>Gauss</td><td>0.995</td><td>0.708</td></tr><tr><td>HANTS</td><td>0.995</td><td>0.826</td></tr></tbody></table>

The correlation coefficient R showed that the S-G smoothing functions got the best results both for the whole image and for sample points among the four smoothing functions conducted in this study.

Correlation analysis between smoothed data and the original data for the 16 fishnet grids also showed that S-G smoothing function presented higher median and mean correlation coefficient than Logistics, Gauss and HANTS as shown in Table 2:

<table><caption>**Table 2.** Correlation coefficient of each time series smoothing function of 16 fishnet grids</caption><thead><tr><th></th><th>S-G</th><th>Logistics</th><th>Gauss</th><th>HANTS</th></tr></thead><tbody><tr><td>Median</td><td>0.83</td><td>0.81</td><td>0.81</td><td>0.77</td></tr><tr><td>Mean</td><td>0.82</td><td>0.81</td><td>0.81</td><td>0.74</td></tr><tr><td>Min</td><td>0.61</td><td>0.62</td><td>0.60</td><td>0.38</td></tr><tr><td>Max</td><td>0.92</td><td>0.94</td><td>0.96</td><td>0.88</td></tr></tbody></table>

S-G smoothing function presented the lowest RMSE than Logistics, Gauss and HANTS as shown in Table 3:

<table><caption>**Table 3.** RMSE of each time series smoothing function of 16 fishnet grid</caption><thead><tr><th></th><th>S-G</th><th>Logistics</th><th>Gauss</th><th>Hants</th></tr></thead><tbody><tr><td>Median</td><td>962</td><td>1076</td><td>1025</td><td>1308</td></tr><tr><td>Mean</td><td>981</td><td>1046</td><td>1005</td><td>1303</td></tr><tr><td>Min</td><td>703</td><td>770</td><td>734</td><td>859</td></tr><tr><td>Max</td><td>1234</td><td>1271</td><td>1222</td><td>1583</td></tr></tbody></table>

Local SNR analysis showed that the original time series had the highest SNR either in mode, median or in mean value. S-G smoothing function presented lower SNR than HANTS smoothing function, but higher SNR than Logistics and Gauss, as shown in Table 4:

<table><caption>**Table 4.** Local SNR of the original and fitted time series images</caption><thead><tr><th>Statics</th><th>Mode</th><th>Median</th><th>Mean</th></tr></thead><tbody><tr><td>Original</td><td>19.586</td><td>16.306</td><td>14.858</td></tr><tr><td>S-G</td><td>15.383</td><td>12.496</td><td>11.686</td></tr></tbody></table>
<table><tr><td>Logistics</td><td>12.062</td><td>11.196</td><td>9.935</td></tr><tr><td>Gauss</td><td>12.232</td><td>11.716</td><td>10.225</td></tr><tr><td>HANTS</td><td>17.521</td><td>14.042</td><td>13.104</td></tr></table>

## 3.2 Classification results

### 3.2.1 Visual check

Corn and soybeans distribution mapped from time series Band 7 to Band 18 of the original images and S-G smoothed images were shown in Fig. 5 and 6. The coarse resolution of MODIS Vis has made the mapping results less specific compared to the CDL distribution of corn and soybeans as shown in Fig. 1. However, the trend of distribution were generally right where there was corn or soybeans. The red ellipse circle on the map showed slightly different details between the results from the original images and S-G smoothed images.

<figure><img src="image_6.png" /><figcaption>Fig. 5. Spectral angle mapper result from the original time series images</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 6. Spectral angle mapper result from the S-G algorithm smoothed time series images</figcaption></figure>

### 3.2.2 Mapping Accuracy

Overall accuracy, producer and user accuracy of classification were compared in three groups, the first of which classified all 23 original and S-G fitted images, the second classified the 12 images from Band 7-18, and the third group classified 7 images including band 9, and Band12-17. Percent and pixel accuracy of corn and soybeans were given, as well as the overall accuracy, as shown in Table 5:

<table><caption>Table 5. Overall, Producer, and User accuracy of classification of three group</caption><thead><tr><th>Group</th><th>Time series</th><th>Crop type</th><th>Prod. Acc.</th><th>User. Acc.</th><th>Overall Acc.</th></tr></thead><tbody><tr><td rowspan="4">The whole 23 bands</td><td rowspan="2">Original</td><td>Corn</td><td>60.970</td><td>69.000</td><td rowspan="2">64.486%</td></tr><tr><td>Soy</td><td>68.530</td><td>60.480</td></tr><tr><td rowspan="2">S-G smoothing</td><td>Corn</td><td>50.630</td><td>66.220</td><td rowspan="2">59.667%</td></tr><tr><td>Soy</td><td>70.110</td><td>55.150</td></tr><tr><td rowspan="4">Band 7 to Band 18</td><td rowspan="2">Original</td><td>Corn</td><td>65.250</td><td>68.910</td><td rowspan="2">65.613%</td></tr><tr><td>Soy</td><td>66.040</td><td>62.270</td></tr><tr><td rowspan="2">S-G smoothing</td><td>Corn</td><td>65.250</td><td>68.910</td><td rowspan="2">60.594%</td></tr><tr><td>Soy</td><td>66.040</td><td>62.270</td></tr><tr><td rowspan="2">Band9,Band 12 to band 17</td><td rowspan="2">Original</td><td>Corn</td><td>65.250</td><td>68.910</td><td rowspan="2">58.138%</td></tr><tr><td>Soy</td><td>66.040</td><td>62.270</td></tr></tbody></table>
<table><tr><td rowspan="2">S-G smoothing</td><td>Corn</td><td>48.210</td><td>64.130</td><td>58.902%</td></tr><tr><td>Soy</td><td>69.400</td><td>54.180</td><td></td></tr></table>

The ratio of mapped acreage from MODIS time series NDVI images to the acreage of CDL and to the statistics reported by ASD released June 28, 2013, by NASS, USDA were given in Table 6.

<table><caption>**Table 6.** Mapped MODIS result acreage compared to pixel counted from CDL and statistical acreage released June 28, 2013, by NASS, USDA</caption><thead><tr><td>Ratio</td><td>Corn</td><td>Soy</td></tr></thead><tbody><tr><td>Original time series To CDL</td><td>95.273%</td><td>74.977%</td></tr><tr><td>S-G smoothing time series To CDL</td><td>96.463%</td><td>73.690%</td></tr><tr><td>PIF to Statistic</td><td>89.477%</td><td>72.415%</td></tr><tr><td>S-G To Statistic</td><td>90.595%</td><td>71.173%</td></tr></tbody></table>

The mapped acreage of corn was very near to the acreage of both CDL and state statistics released in June 28, 2013, by NASS, USDA.

# 4. Discussion

## 4.1 Time series fitting effectiveness

Visual checking showed that HANTS altered the original data value massively, and thus has shown higher SNR than S-G, Gauss and Logistics. Visual checking by time series profiles showed that S-G, Gauss, and Logistics were similar in the growing season, but differed from each other in other part of the year. Generally, S-G smoothing function got the best results, especially with NDVI peak value neither obviously lagging nor postponing.

Statistical check showed that S-G smoothing algorithms had the highest correlation coefficient either for the whole image or for sample points. S-G smoothing algorithms also presented the lowest RMSE and relatively lower SNR. With these results, S-G was thought to be the best and stably effective smoothing function for crop mapping.

## 4.2 Crop mapping accuracy

With the same methods, the same input data bands, and the same training data, the overall accuracy of the original time series was generally better than the S-G smoothed time series, but the accuracy was not obviously high, no higher than 5%. The ratio of mapped acreage to CDL pixel counted acreage was near 100% for corn, and near 80% for soybeans, and so was the ratio of mapped acreage to national statistical acreage. Besides, the detail in the distribution map also led to different conclusion compared to the accuracy report. Therefore, it was concluded here that the
smoothing functions, the mapping method, the training data and the accuracy check method all had effect on the comparison. It is the goal the matters for the choosing of mapping method.

## 4.3 Limitations

Both time series fitting method and the selection of input bands have effect on the mapping result. This paper only investigated the result of spectral angle mapper mapping method. Other method may got different results as well.

# 5. CONCLUSION

This paper compared the effectiveness of four time series smoothing algorithms through a series of statistical indicators, and used the most effective smoothed result to conduct crop mapping. S-G smoothing was found to be the most effective fitting method. However, the most effective smoothed results did not provide the most accurate mapping results with the same method, the same training data and the same checking data. This findings imply that both mapping methods and the goals or purposes the user require matter more for a better result of crop mapping. It was not necessary to conduct a time series smoothing before a better mapping scheme was founded.

## Acknowledgment

We express our deep gratitude to the free data and software providers: NASA for MODIS, USDA for CDL, GDSC (Geospatial Data Service Centre) for HANTS package, and Lars Eklundh for the software TIMESAT. The whole work was funded by the Civil Space Project in the 12th Five-Year Plan of China (2011-2015)

## References

1. C. Atzberger, and P. H. Eilers, 2011. Evaluating the effectiveness of smoothing algorithms in the absence of ground reference measurements. International Journal of Remote Sensing 32:3689-3709.

2. P. S. A. Beck, C. Atzberger, K. A. Høgda, B. Johansen, and A. K. Skidmore. 2006. Improved monitoring of vegetation dynamics at very high latitudes: A new method using MODIS NDVI, Remote Sensing of Environment 100:321-334.

3. J. C. Brown, J. H. Kastens, A. C. Coutinho, D. d. C. Victoria, and C. R. Bishop, 2013. Classifying multiyear agricultural land use data from Mato Grosso using time-series MODIS vegetation index data, Remote Sensing of Environment 130:39-50.
4. J. Chang, M. C. Hansen, K. Pittman, M. Carroll, and C. DiMiceli, 2007. Corn and soybean mapping in the United States using MODIS time-series data sets, Agronomy Journal 99:1654-1664.

5. J. Chen, P. Jönsson, M. Tamura, Z. Gu, B. Matsushita, and L. Eklundh, 2004. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter, Remote Sensing of Environment 91:332-344.

6. C. H. W. de Souza, E. Mercante, J. A. Johann, R. A. C. Lamparelli, and M. A. Uribe-Opazo, 2015. Mapping and discrimination of soya bean and corn crops using spectro-temporal profiles of vegetation indices, International Journal of Remote Sensing 36:1809-1824.

7. M. Duggin, 1985. Review Article: Factors limiting the discrimination and quantification of terrestrial features using remotely sensed radiance. International Journal of Remote Sensing, 6(1), 3-27.

8. A. Gusso, D. Arvor, J. Ricardo Ducati, M. R. Veronez, and L. G. Da Silveira, 2014. Assessing the MODIS crop detection algorithm for soybean crop area mapping and expansion in the Mato Grosso State, Brazil, The Scientific World Journal 2014.

9. J. N. Hird, and G. J. McDermid, 2009. Noise reduction of NDVI time series: An empirical comparison of selected techniques, Remote Sensing of Environment 113:248-258.

10. A. Huete, K. Didan, T. Miura, E. P. Rodriguez, X. Gao, and L. G. Ferreira, 2002. Overview of the radiometric and biophysical performance of the MODIS vegetation indices, Remote Sensing of Environment 83:195-213.

11. P. K. Jönsson, and L. Eklundh, 2002. Seasonality extraction by function fitting to time series of satellite sensor data, IEEE Transactions on Geoscience and Remote Sensing 40, 1824-1832.

12. P. Jönsson, and L. Eklundh, 2004. TIMESAT—a program for analyzing time-series of satellite sensor data, Computers & Geosciences 30:833-845.

13. C. Jentsch, and S. Subba Rao, 2014. A test for second order stationarity of a multivariate time series, Journal of Econometrics.

14. F. A. Kruse, A. B. Lefkoff, J. W. Boardman, K. B. Heidebrecht, A. T. Shapiro, P. J. Barloon, and A. F. H. Goetz, 1993. The spectral image processing system (SIPS)—interactive visualization and analysis of imaging spectrometer data, Remote Sensing of Environment 44:145-163.

15. T. Sakamoto, B. D. Wardlow, A. A. Gitelson, S. B. Verma, A. E. Suyker, and T. J. Arkebauer, 2010. A Two-Step Filtering approach for detecting maize and soybean phenology with time-series MODIS data, Remote Sensing of Environment 114:2146-2159.

16. B. D. Wardlow, and S. L. Egbert, 2008. Large-area crop mapping using time-series MODIS 250 m NDVI data: An assessment for the US Central Great Plains, Remote Sensing of Environment 112:1096-1116.

17. J. Zhang, L. Feng, and F. Yao, 2014. Improved maize cultivated area estimation over a large scale combining MODIS-EVI time series data and crop phenological information, ISPRS Journal of Photogrammetry and Remote Sensing 94:102-113.

18. J. Zhou, L. Jia, and M. Menenti, 2015. Reconstruction of global MODIS NDVI time series: Performance of Harmonic Analysis of Time Series (HANTS), Remote Sensing of Environment 163:217-228.

19. H. Zhang, and Z.Y. Ren, 2014. Comparison and Application Analysis of Several NDVI Time-Series Reconstruction Methods, Scientia Agricultura Sinica 47(15):2998-3008. (In Chinese)

20. B. Zhu, X.H. Wang, L.L.Tang, and C.R. Li, 2010. Review on Methods for SNR Estimation of Optical Remote Sensing Imagery, Remote Sensing Technology and Application 25(1):303-309. (In Chinese)