## 3 Method
### Data Process for the Reconstruction of Time Series High Spatial Resolution Data
Systematic biases were corrected in the NDVI time series from MODIS, Landsat, HJ, and GF-1 observations because these sensor systems have differences in orbit parameters, bandwidth, acquisition time, and spectral response function. Several studies had addressed the problem of sensor difference using linear sensor-adjustment method [24,30]. However, because the difference between the reflectance of vegetation and non-vegetation in HJ data is much smaller than that of Landsat satellites, a large error occurs in using linear sensor-adjustment methods. In this case, we calibrated the sensor differences one-by-one according to land cover types. The sensor calibration process includes three steps. First, vegetation and non-vegetation were classified using a threshold method for NDVI data. Pixels with NDVI greater than a dynamic threshold were classified as vegetation, while other pixels were classified as non-vegetation. The dynamic threshold value was calculated using regions of interest (ROIs) selected by the visual-interpretation method. Second, linear sensor-adjustment model of vegetation and non-vegetation land cover types were built using the images acquired at the same date, respectively. Landsat-HJ NDVI image pairs acquired on 22 April 2011 in Bole and 19 August 2013 in Luntai were used to build those models for the calibration of Landsat data, respectively. The GF-HJ NDVI image pair acquired on 7 October 2013 in Luntai was used to build the sensor-adjustment model for the calibration of GF-1 data. Finally, these models were used to calibrate all the Landsat and GF-1 data, due to the number of HJ data is higher than the number of Landsat and GF-1 data. To fill the temporal gaps, the STDFA was used to fuse MODIS with Landsat, HJ, and GF-1 NDVI data using the following formula [36]:  
\[r(x,y,t_{k})=r(x,y,t_{0})+\overline{r}(c,t_{k})-\overline{r}(c,t_{0}) \tag{1}\]  
where \(r(x,y,t_{k})\) and \(r(x,y,t_{0})\) are the fine-resolution reflectances or NDVI of pixels at location (x,y) from time \(t_{0}\) to \(t_{k}\); \(\overline{r}(c,t_{0})\) and \(\overline{r}(c,t_{k})\) are the mean reflectances or NDVI of land cover class \(c\) from time \(t_{0}\) to \(t_{k}\) which can be calculated using the Ordinary Least Squares technique according linear mixing model [42]:  
\[R(x,y,t)=\sum_{c=0}^{M}f_{c}\left(x,y,c\right)\times\overline{r}(c,t)+\xi \left(x,y,t\right) \tag{2}\]  
\[\text{Constrained:}\sum_{c=0}^{M}f_{c}(x,y,c)=1\text{ and }1\geq f_{c}(x,y,c)\geq 0 \text{ for all.}\]  
where \(\xi(x,y,t)\) is the residual error term; \(f_{c}(x,y,c)\) is the fractional cover of class \(c\) in coarse pixel at location (x,y); \(R(x,y,t)\) is the coarse spatial reflectance at time \(t\); and cover types. By inputting two medium-resolution data and time-series MODIS data, time series synthetic medium-resolution data can be generated by STDFA. One medium-resolution data is used for the base image. The other medium-resolution data is used for land cover mapping. To generate synthetic medium-resolution data from time \(t_{1}\) to \(t_{n}\), one actual medium-resolution data with the shortest date interval to this period were used as the base image. The medium-resolution data with the second shortest date interval to this period was selected as another set of medium-resolution data. MODIS data acquired from time \(t_{1}\) to \(t_{n}\) were used as the time series MODIS data. Due to this, data fusion is more suitable to generate time series of NDVI data than reflectances [37, 40], daily 30 m NDVI data were generated by fusion MODIS with Landsat, HJ, and GF-1 NDVI data using STDFA. Finally, the time series of high spatial resolution NDVI data sets at each pixel were smoothed using a Savitzky-Golay filter.