## 3 Method
### Data Process for the Reconstruction of Time Series High Spatial Resolution Data  
A time series of high spatial resolution datasets were generated by combining HJ CCD, GF-1 WFV, Landsat, and MODIS data. The data process includes seven steps: (1) good data selection; (2) atmospheric correction; (3) vegetation index calculation; (4) geometric correction; (5) sensor calibration; (6) spatial and temporal data fusion; and (7) data smoothing. The data acquired during clear sky conditions were selected. Since there are no cloud products for Landsat, HJ, and GF-1 satellites, a time series of the good quality data were identified from these high spatial resolution data using the visual-interpretation method. Moreover, images with snow cover were also visually removed. The 500 m daily MODIS surface reflectance products (MOD09 GA) acquired during clear sky conditions were selected using the data quality of reflectances. Then, the selected Landsat, HJ, and GF-1 satellites data were atmospherically corrected using the Fast Line-of-Sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) Atmospheric Correction Model. MOD09 GA is a level 2G MODIS products that has been atmospherically corrected, so no additional atmospheric correction was needed for those data. These MODIS images were reprojected from the native sinusoidal projection to a UTM-WGS84 reference system, and were resized to the \begin{table}
\begin{tabular}{c c c} \hline \hline
**Sensor** & **Images for Bole** & **Images for Luntai** \\ \hline Huanjing satellite charge coupled device (HJ CCD) & 41 & 55 \\ Gaofen satellite no. 1 wide field-of-view camera (GF-1 WFV) & 0 & 18 \\ Landsat Thematic Mapper (TM)/ Operational Land Imager (OLI) & 6 & 10 \\ Moderate Resolution Imaging Spectroradiometer (MODIS) & 154 & 122 \\ \hline \hline \end{tabular}
\end{table} 
Table 1: The number of the four kinds of satellite data used in Bole and Luntai. \begin{table}
\begin{tabular}{c c c c c c} \hline \hline
**Sensor** & **Band** & **Bandwidth (nm)** & **Resolution (m)** & **Width (km)** & **Return Cycle (Days)** \\ \hline \multirow{4}{*}{HJ CCD} & 1 & 430–520 & & & \\  & 2 & 520–600 & & & \\  & 3 & 630–690 & & & \\  & 4 & 760–900 & & & \\ \hline \multirow{4}{*}{GF-1 WFV} & 1 & 430–520 & & & \\  & 2 & 520–600 & & & \\ \cline{1-1}  & 3 & 630–690 & & & \\ \cline{1-1}  & 4 & 760–900 & & & \\ \hline \hline \end{tabular}
\end{table} 
Table 2: Parameters for HJ satellite constellation CCD sensor and GF-1 satellite WFV sensors. selected study area. The NDVI was then calculated from the time series of Landsat, HJ, GF-1, and MODIS data. The geometric corrections of Landsat, HJ, and GF-1 data were conducted in two steps. First, one set of Landsat NDVI data was georeferenced using a second-order polynomial warping approach based on the selection of 53 ground control points (GCPs) from a 1:10,000 topographic map using the nearest neighbor resampling method, with a position error of within 0.51 Landsat pixels. Next, the warped Landsat NDVI data was used as a reference image in the geometric corrections of other data. The geometric corrections of other data were conducted using the automatic registration module. This module can select GCPs automatically. Then, these automatically-selected GCPs were checked manually to determine whether they were real GCPs. Wrong GCPs were replaced using GCPs selected manually. The 16 m GF-1 data were resized to 30 m with the pixel aggregate resampling method before the geometric corrections. The MODIS data were georeferenced by a second-order polynomial warping approach based on the selection of 43 GCPs on 500 m Landsat images, with a nearest neighbor resampling method in which the position error was within 0.69 500 m Landsat pixels. The 500 m Landsat images were resized from georeferenced Landsat images using the pixel aggregate resampling method.