Remote Sensing Monitoring of Vegetation Dynamic Changes after Fire in the Greater Hinggan Mountain Area: The Algorithm and Application for Eliminating Phenological Impacts

: Fires are frequent in boreal forests affecting forest areas. The detection of forest disturbances and the monitoring of forest restoration are critical for forest management. Vegetation phenology information in remote sensing images may interfere with the monitoring of vegetation restoration, but little research has been done on this issue. Remote sensing and the geographic information system (GIS) have emerged as important tools in providing valuable information about vegetation phenology. Based on the MODIS and Landsat time-series images acquired from 2000 to 2018, this study uses the spatio-temporal data fusion method to construct reflectance images of vegetation with a relatively consistent growth period to study the vegetation restoration after the Greater Hinggan Mountain forest fire in the year 1987. The influence of phenology on vegetation monitoring was analyzed through three aspects: band characteristics, normalized difference vegetation index (NDVI) and disturbance index (DI) values. The comparison of the band characteristics shows that in the blue band and the red band, the average reflectance values of the study area after eliminating phenological influence is lower than that without eliminating the phenological influence in each year. In the infrared band, the average reflectance value after eliminating the influence of phenology is greater than the value with phenological influence in almost every year. In the second shortwave infrared band, the average reflectance value without phenological influence is lower than that with phenological influence in almost every year. The analysis results of NDVI and DI values in the study area of each year show that the NDVI and DI curves vary considerably without eliminating the phenological influence, and there is no obvious trend. After eliminating the phenological influence, the changing trend of the NDVI and DI values in each year is more stable and shows that the forest in the region was impacted by other factors in some years and also the recovery trend. The results show that the spatio-temporal data fusion approach used in this study can eliminate vegetation phenology effectively and the elimination of the phenology impact provides more reliable information about changes in vegetation regions affected by the forest fires. The results will be useful as a reference for future monitoring and management of forest resources.


Introduction
Forests play an irreplaceable role in maintaining the ecological balance of the terrestrial biosphere due to their wide coverage, complex distribution, and species diversity [1,2], multifunction, and multi-value characteristics [3].
Fires are one of the serious disturbances globally and are particularly prevalent in boreal forests [4]. Forest fires promote dynamic changes in ecosystem structure and function, have positive and negative impacts on ecosystems, and have a profound impact on human life and regional developments [5][6][7][8]. On the one hand, it poses a severe health hazard to people living in the surroundings [9], impacts forest ecosystems, and burns an average of about 350 million ha of forest land per year, which is one of the primary causes for the decline in global forest stocks. The burning of forests severely causes local economic losses. Forest fires also have long term environmental and climate impacts. A certain frequency and intensity of fire can maintain the balance of forest ecology and play an essential role in preserving biodiversity. For example, fires can help regulate fuel accumulation, restore vegetation by removing fungi and microorganisms, control diseases and insects, and gain more energy by exposure to solar radiation, mineral soil, and nutrient release [10]. With climate change and global warming, the frequency of forest fires is increasing and receives increasing attention as an integral part of global environmental change studies [11,12].
Global fires emit about 2.1 PgC (1 PgC = 1015 g of carbon) of carbon flux per year, which is equivalent to 50%-200% of annual terrestrial carbon sinks. Among them, 35% are related to forests [13]. The post-fire forest regeneration process is extremely important. Carbon is emitted from forest fires that are injected into the atmosphere, post-fire vegetation regeneration and carbon sequestration of woody vegetation may help to reduce carbon emission in the atmosphere [14]. Forest disturbance and restoration can, therefore, affect the energy flow and biogeochemical cycles and thus is considered as the primary mechanism for the carbon transfer between the surface and the atmosphere, playing an important role in regional and global carbon cycles [15][16][17]. The detection of forest disturbances and the monitoring of post-fire forest restoration are essential for both ecological research and forest management. Understanding the dynamics of forest regeneration after a fire can help to assess forest resilience and adequately guide forest management after the disturbance. Therefore, information about the spatial patterns and temporal trends of the forest helps in restoration after a fire.
Considering the small spatial coverage, limited sample points, low site accessibility, and high labor costs, site sampling is not suitable for monitoring large-scale vegetation dynamics after a fire, while satellite remote sensing provides an economical and effective tool for monitoring large-scale forest changes [18,19]. Fires can cause profound changes in ecosystems, where vegetation is consumed, leaf chlorophyll is destroyed, the soil is further exposed, and carbonization and moisture changes in vegetation roots cause a large number of spectral changes, which can be detected through satellite data [20,21]. Optical remote sensing data, such as the widely used Landsat images, have been proven to be very suitable for forest interference detection and forest change monitoring because they have the necessary spatial resolution (the resolution of 30 m, consistent with the scale of most local vegetation changes [22]) and spectral coverage (visible, near-infrared, short-wave infrared, and thermal infrared bands) to capture most forest disturbance events caused by natural or artificial management [23]. At the same time, time-series remote sensing data (40-year observations) (e.g., Landsat) provides excellent potential for trajectory monitoring of forest dynamics after fire [24,25], as long-term monitoring of forest recovery after a fire is often required [26,27].
Numerous studies have been carried out on different aspects of post-fire forest restoration [23,[28][29][30][31]. Many scholars have also studied the phenology of vegetation. The experimental results of Frison et al. [32] show that radar data is more accurate for phenological estimation than optical data. Flavio et al. [33] have tested the phenology-based vegetation mapping method and proved it effective. Some studies have calculated the phenological characteristics of mangroves to derive environmental driving factors that affect their growth [34,35]. Vegetation phenology can provide information about the vegetation dynamics and response after forest fires [36]. However, there is less research on analyzing the influence of phenological factors caused by remote sensing data on the monitoring of dynamic vegetation changes. Due to the significant differences in vegetation phenology between different growth stages, in order to avoid the "pseudo-variation" of the timeseries vegetation index caused by the interannual vegetation phenological changes, some studies have chosen to use the images acquired at the time near the vegetation growth peak to monitor the post-fire forest recovery [37]. Remote sensing and geographic information tools have emerged as important tools to study vegetation phenology using long time-series of vegetation indexes to monitor the post-fire forest recovery [37]. However, due to current technical limitations, it is challenging to obtain remote sensing data with high spatial resolution and high temporal resolution simultaneously. The coarse resolution (e.g., MODIS, 250 m/500 m/1000 m) will obscure the details of the features and affect the observation results. The long revisit period (16 days) of satellites (e.g., Landsat), frequent cloud pollution, and other atmospheric conditions limit their application in long time-series detection of surface objects without phenological interference. Therefore, long-term observations by only using images located near the vegetation growth peak in cloudy areas may result in a gap in the study years.
Taking the forest restoration in the Greater Hinggan Mountain area after the "5.6 fire" in 1987 as an example, this study aims at demonstrating the effect of the spatiotemporal fusion algorithm in eliminating the phenological impact when monitoring vegetation restoration using remote sensing images. Here, we used the Landsat and MODIS time-series images to study the vegetation during 2000-2018 based on the spatiotemporal fusion algorithm to eliminate the influence of phenological factors in Landsat images. We compared the band characteristics, NDVI and DI indices prior to and after the elimination of the phenology effect, and further explore the impact of phenology on forest dynamic monitoring. The results of this study prove that the spatiotemporal fusion algorithm can effectively eliminate phenological factors in remote sensing images. The elimination of the phenological effects can provide more reliable information on vegetation restoration. Thus, the present study provides a scientific reference for post-fire forest reconstruction and ecological restoration.

Study Area
The Greater Hinggan Mountain area is located in Heilongjiang Province, in the northern part of Inner Mongolia Autonomous Region which is the watershed of the Mongolian Plateau and the Songliao Plain bounded by latitude 50°10′N to 53°33′N, longitude 121°12′E to 127°00′E ( Figure 1). The area is more than 1200 km long and 200-300 km wide with an average altitude of 1200-1300 m above mean sea level. The Greater Hinggan Mountain area is a typical cold temperate continental monsoon climate with warm summers and cold winters. The annual average temperature of the area is −2.8 °C; the lowest temperature is −52.3 °C. The precipitation, which peaks in summer, is 420 mm annually and is unevenly distributed throughout the year, i.e., more than 60% occurs between June and August [38]. The Greater Hinggan Mountain is the largest modern state-owned forest area with a total area of 8.46 × 10 4 square kilometers and forest coverage of 6.46 × 10 4 square kilometers. Therefore, the forest coverage rate is about 76.4%, and the total storage capacity is about 5.01 × 10 8 m 3 , accounting for 7.8% of the national total [39]. The Greater Hinggan Mountain covers large forest resources, serving as an important stateowned forest area with a vital timber production area in China. At the same time, this forest area has experienced one of the most severe forest fires in China. On 6 May 1987, a severe forest fire occurred in the northern part of the Greater Hinggan Mountain. The burned area was 1.133 × 10 7 km 2 and the area of over-fired forest land was 1.114 × 10 7 km 2 , of which the affected area was 8.17 × 10 5 km 2 . The fire seriously affected the social, economic, and ecological benefits of the forest area, causing unprecedented heavy losses to the country. Since the catastrophic forest fires in the Greater Hinggan Mountains happened in 1987, this place has been one of the areas for research on fire prevention and post-fire forest management [38,40,41].
The burned area of the "5.6 Fire" was extracted in a previous study [42]. The entire burned forest area spanned two Landsat scenes (Path 121/122, Row 23), but it is difficult to acquire the two scenes simultaneously in each year. Considering that around 90% of the burned forest area is within the scene of path 122 row 23, we extracted a sample area ( Figure 1) from Landsat path 122 row 23 as the study area for the recovery monitoring [39].

Data Used and Preprocessing
A total of 16 Landsat surface reflectance data from Path 122, Row 23 with and below 10% cloud cover during the vegetation growth period from 2000 to 2018 was considered. The data was downloaded from the United States Geological Survey (USGS, https://earthexplorer.usgs.gov/); details of the data are given in Table 1. We have tried our best to find the images with the least cloud volume during the vegetation growth period each year. Although the cloud volume of the image in 2001 is 10%, the study area is only part of an image, and most of the clouds are located outside the study area. Therefore, the image from 2001 was still used in the study. The Landsat surface reflectance data was corrected at the sub-pixel-level by topographic and atmospheric correction [43,44]. The FMASK algorithm was used to detect cloud cover and cloud shadow and to generate a mask [45,46]. The MODIS 16d synthetic vegetation index product MOD13Q1 for the periods 2000-2018 was downloaded from the National Aeronautics and Space Administration (NASA) and pre-processed. The zenith BRDF-adjusted reflectance product MCD43A4V006 with a spatial resolution of 500m was obtained from NASA, which is daily reflectance data for spatial and temporal fusion with Landsat data to generate the surface reflectance on the target date. The above MODIS data was converted from Sinusoidal projection to UTM projection with WGS84-51N coordinates. After all the detailed processing, the river, road and building areas in each image were masked based on the 10 m global resolution land cover data [47], supplemented by visual interpretation, and the boundary of the study area was extracted.

Vegetation Phenological Information
Two processes were used to extract phenological information: the smooth reconstruction of the temporal vegetation index and the extraction of phenological parameters. Previous scholars have done a lot of research on the smooth reconstruction of NDVI time-series data. Methods including least squares (i.e., Savitsky-Golay filtering, asymmetric Gaussian function fitting, logistic blending function fitting), and Fourier fitting, Fourier correction algorithms, harmonic analysis method and wavelet analysis based on spectrum analysis technology have been considered. There has also been a large number of studies on the extraction of phenological parameters [48][49][50]. Based on the comparison of all earlier methods, we used an adaptive Savitzky-Golay filter to reconstruct the MOD13Q1 NDVI time-series and used a dynamic threshold method to extract the vegetation phenological index of each year. For the processing of data, we used TIMESAT software [51].
The vegetation index has a large degree of uncertainty showing the highest value of the NDVI peak, and the determination of the beginning and end dates of the growth period is relatively easy; the determination of the mid-point of the growth period is more reliable, which is often located in the peak season of vegetation growth. Vegetation at the midpoint of the growth period in the study area in each year has a relatively consistent growth situation [52]. Comparing the remote sensing indices at this time of each year can effectively eliminate the phenological influence. Therefore, the midpoint of the vegetation growth period in the study area in each year was selected as the date of the image to be synthesized.

Synthesis of Target Image Based on STARFM Fusion Algorithm
Traditional image fusion methods, such as intensity-hue-saturation (IHS) transformation [53], principal component substitution (PCS) [54], and wavelet decomposition [55], focus on combining the spectral properties of low-resolution data with the high spatial resolution of panchromatic images to generate high-resolution multispectral images. These methods are useful for exploiting the different spectral and spatial characteristics of multi-sensor data, but they cannot enhance both the spatial resolution and temporal coverage. In the present study, we have quantitatively captured the changes in radiation measurements (surface reflections) associated with the phenology and study their effects on the monitoring of vegetation restoration associated with fire. The STARFM (spatial and temporal adaptive reflection fusion model) algorithm is used to predict the reflectance of the target date. The algorithm considers the influence of spatial distance on predicted pixels and also considers the spectral difference and temporal difference between pixels. The homogeneous pixels in MODIS data show the relationship with the corresponding Landsat pixels as where ( , ) represents the spatial position of the homogeneous pixel; represents the image acquisition time; ( , , ) represents the reflectivity of Landsat pixels; ( , , ) represents the reflectivity of MODIS pixel; indicates the difference in reflectance between different data. At , the following is the relationship between the MODIS reflectance and the Landsat reflectance of the same pixel: When the ground cover type and the system error between the two types of data remain unchanged, = , the above equation can be expressed as The MODIS pixels are mostly non-homogeneous pixels, and the solar bidirectional reflection changes, and the surface coverage type changes with time, which makes the above ideal conditions challenging to meet. Therefore, the key point of the method is to find similar pixels from neighboring pixels of the target pixel and replace the homogeneous pixels with similar pixels. We used a windowbased threshold method to search for similar pixels from the window. If the pixel in the moving window satisfies the following relationship, the pixel is considered to be a similar pixel of the target pixel.
where represents the size of the moving window; / , / represents the position of the predicted pixel; represents the standard deviation of Landsat surface reflectance; represents the classification number of ground objects in the moving window. Thus, the reflectance value of the predicted pixel can be expressed by the following equation: where represents the number of similar pixels in the window; represents the contribution weighting coefficient of the neighboring pixels to the target pixel. The weighting coefficient can be calculated using three factors: spectral distance, temporal distance, and spatial distance between adjacent pixels and central pixels. Spectral distance is the spectral difference of pixels between simultaneous MODIS and Landsat data in the same location. The MODIS pixel reflectance can be considered as a mixture of multiple Landsat pixel reflectance in the same region. The smaller the spectral distance, the more similar the Landsat pixel and target pixel, and the larger the weight coefficient assigned. Time distance is the difference between MODIS pixel values at different times, which represents the change of the surface coverage status in this period. The smaller the time distance value, the smaller the change of land cover, the larger the contribution of the pixel to the central pixel value, and the larger the weight coefficient assigned. Spatial distance is the distance between the neighboring pixels and the target pixel. The smaller the spatial distance, the larger the weighting coefficient, which is calculated as follows: where is the spectral distance; is the temporal distance; is the spatial distance; ( , ) is the spatial position of similar pixels; is the weight adjustment coefficient, which is a constant. Normalized weight coefficients ( ), and the total weight coefficient ( ) are given as After selection of similar pixels, we filtered to remove the poor-quality pixels. If the spectral and time distance of similar pixels are smaller compared to the target pixel in the center of the moving window, the pixel provides better spectral information and time information compared to the target pixel. Otherwise, the pixel is an unqualified similar pixel. When the uncertainty factors and of Landsat and MODIS surface reflectivity are introduced in the similar pixel screening, the qualified similar pixels must satisfy the following inequality relations: where represents the uncertainty factor between the MODIS and the Landsat reflectance value, and represents the uncertainty factor of the MODIS reflectance at different phases. When all observed pixel reflectance values are independent of each other, and are expressed as After extracting the phenological index of the vegetation, the corresponding date of the midgrowth period of each year can be calculated as the date of the image to be synthesized, and the Landsat and MCD43A4V006 data are fused to construct reflectance data of vegetation with a relatively consistent growth period in different years.

Vegetation Indices
We used NDVI and DI to characterize post-fire vegetation restoration status.

NDVI
As one of the most famous vegetation indices, NDVI shows a good correlation with vegetation regeneration and the photosynthetic effective radiation ratio absorbed by plant canopy, leaf area, and biomass, so it is widely used to study vegetation response to wildfire disturbance [23,26,[56][57][58][59]. NDVI is calculated using Equation (15): where and are the reflectance of the near-infrared band and red wavelengths, respectively.

DI
The calculation of DI is based on the Tasseled Cap transformation [60,61], which is a spectral transformation that converts the original high covariant data into three uncorrelated indices known as brightness (B), greenness (G), and wetness (W). The calculation of DI is based on the observation that disturbed forests usually have higher brightness values and lower green and humidity values compared to undisturbed forest areas [61]. The linear combination of three Tasseled Cap transformation indices include brightness, greenness, and wetness. At the same time, the spectral normalization step is conducted, and the intra-image statistics are used to normalize the radiation variations as where , , and represent the average Tasseled Cap transformation brightness, greenness, and wetness of the "forest in a particular scene"; , , and are the corresponding standard deviations, so , , and represent normalized brightness, greenness, and wetness, respectively. After normalization, the three components are linearly combined to obtain DI as follows: The disturbed forest area usually has a high positive value and low negative values of and , thus showing a high DI value; in contrast, the undisturbed forest area shows a low DI value.

Yearly Composite Image
The date corresponding to the midpoint of the vegetation growth period in each year was obtained from the vegetation index. The reflectance data of vegetation with relatively consistent growth periods in different years were constructed by integrating Landsat and MODIS data. The image acquisition date of each year, the mid-point of the vegetation growth period, and the number of days between them are given in Table 2 (represented by the number of days in a year corresponding to the date). Since the MCD43A4V006 data has data gaps on some dates, in order to have the fusion image as complete as possible, the MCD43A4V006 image nearest to the original date and with the least data gaps was found close to the original date. The adjusted data date was marked in brackets of the original date.  (180)  37  2015  242  196  46  2016  260  213  47  2017  264  202  62  2018  123  202  79 The midpoint of vegetation growth in the study area in each year was almost always around the 200th day of the year ( Table 2)

The Characteristics of Reflectance Prior to and after Eliminating Phenological Influence
The average values of all bands in the study area prior to (original image) and after (fusion image) the elimination of the phenological influence in each year are given in Table 3 and are shown in Figure 2. In the blue and the red bands (Figure 2a,c), the average reflectance value of the study area after eliminating phenological influence is lower than that of the area with phenological influence in   Figure 2d shows that the average reflectance value of the study area after eliminating phenological influence in the study area was greater than that without eliminating phenological influence in the near-infrared band for almost every year except the years 2006 and 2011, indicating that the vegetation has a stronger ability to reflect near-infrared signals during the mid-life phase. At the same time, after eliminating the influence of phenology, the inter-annual change curve of reflectance in the near-infrared band becomes more gradual. The inter-annual reflectance variance of the near-infrared band with phenological influence is 0.0026, while reduced to 0.0002 when the phenological influence is eliminated. Figure 2e,f show that the interannual change of the reflectance in the two short wave infrared bands is more gradual after eliminating the influence of phenology. There is a water absorption band near both short-wave infrared bands, and the one near the second short-wave infrared band (1.9 μm) has stronger water absorption than the one near the first short-wave infrared band (1.4 μm). Figure  2f shows that in the second short-wave infrared band, the average reflectance of the study area after eliminating the phenological influence is almost lower than that prior to eliminating the phenological influence in every year.

NDVI Characteristics Prior to and after Phenological Influence Elimination
The mean values of NDVI in the study area prior to and after the phenological influence elimination in each year are given in Table 4 and shown in Figure 3. The changing trend of NDVI values of each year was more stable after the elimination of phenological influences. The NDVI without eliminating phenological effects shows a significant decreasing trend compared with the previous year for the years in 2002, 2009, and 2016, while the NDVI curve without phenological effects shows a slight increasing trend in these years. At the same time, it can be seen from Table 4 that the difference of NDVI prior to and after eliminating phenological influence was up to 0.4 or more, indicating that the impact of phenology on vegetation monitoring cannot be ignored.

Characteristics of DI Changes Prior to and after Phenological Influence Elimination
The mean values of DI in the study area prior to and after the phenological influence elimination in each year were calculated and are shown in Table 5 and Figure 4. Combined with Table 5 and Figure 4, it can be seen that the changing trend of DI values in each year is more gentle after eliminating the influence of phenology. The DI variations ( Figure 4) without eliminating phenological influence show a significant upward trend in the years 2002, 2016, and 2018 compared with the previous year, while the DI variations after eliminating phenological effects are relatively flat in these years, with no obvious increase or decrease. Table 5 shows the difference in DI prior to and after eliminating phenological influence reached the maximum value of 4.067 in the year 2018, and the corresponding image acquisition date is 79 days prior to the midpoint of the vegetation growth period.  The DI variations after eliminating the impact of phenology in the study areas decreased from

Discussion
The comparison of the band characteristics shows that in the blue band and the red band, the average reflectance values of the study area after eliminating phenological influence was lower than that without eliminating phenological influence in each year, indicating that the vegetation has a stronger absorption ability. In the infrared band, the average reflectance value after eliminating the influence of phenology was greater than the value of the unremoved phenological influence in almost every year. In the second shortwave infrared band, the average reflectance value after eliminating phenological influence was lower than that with phenological influence in almost every year. Since this study used the corresponding date of the midpoint of the vegetation growth period of each year in the study area as the target date of image fusion, the mid-growth period of vegetation is often the period of peak vegetation growth, and most of the acquisition dates of the images cannot be in the peak period of vegetation growth. Therefore, vegetation located at the midpoint of the growing season tends to have better growth compared to the data acquisition date.
At the same time, due to the influence of chlorophyll, plant structure, and water absorption, the vegetation at the midpoint of the growing season has stronger absorption in the blue, red, and shortwave infrared bands and stronger reflection in the near-infrared band compared with the image acquisition date. Meanwhile, in the fused image, the reflectance values of several bands (red band, near-infrared band, shortwave infrared band) tend to be stable, which show a great relationship with the state of vegetation growth, indicating that the method effectively eliminates the disturbance caused by phenology influence in the study of interannual growth and change of vegetation.
When determining the fusion methodology, we considered a number of models. Finally, STARFM, developed by Gao et al. [62] combining Landsat and MODIS data to predict the daily surface reflectance at Landsat spatial resolution and MODIS temporal frequency, was considered. This method was tested in a conifer-dominated region in central British Columbia, Canada, and proved to generate daily surface reflectance with the same spatial resolution as Landsat data. The generated reflectance data is in good agreement with the actual Landsat reflectance data. Zurita-Milla et al. [63] developed another downscaling algorithm based on a linear hybrid model to produce images with medium resolution imaging spectrometer (MERIS) spectral characteristics and similar Landsat time resolution. However, this reduction algorithm requires high resolution land-use data for pixel unmixing and may not be suitable for many applications. The STARFM method does not require any auxiliary data compared to the downscaling algorithm. Zhu et al. [64] developed an enhanced spatial and temporal adaptive reflection fusion model (ESTARFM) based on the STARFM algorithm and tested the simulated and actual satellite data. The results show that ESTARFM improves the accuracy of reflectivity prediction, especially for heterogeneous landscapes. Taking the NIR band as an example, the ESTARFM prediction for a uniform region is slightly better than STARFM (average absolute difference (AAD) 0.0106 vs. 0.0129 reflection units); for complex heterogeneous environments, the prediction accuracy of ESTARFM was further improved compared with STARFM (AAD 0.0135 vs. 0.0194). Although the prediction accuracy of ESTARFM is slightly higher than STARFM, the former requires the input of one image before and after the predicted image, respectively, while the latter allows only one image of the time near the predicted image to be input. Considering that ESTARFM has high requirements for data, it cannot be fully realized in all the years of the study area, and most of the study area is covered by vegetation with a relatively homogeneous ground status, hence the STARFM is found to be more applicable to achieve the objectives of the study. In the future, the ESTARFM methodology can be tried when there is enough data as input.
Although annual 30 m reflectance data of the study area cannot be obtained in this study to verify the accuracy of the fusion image, on the one hand, Gao et al. [62] tested the methodology in central British Columbia, Canada, and found that the daily surface reflectance generated by this method is in good agreement with the actual Landsat data. On the other hand, the 2011 image acquisition date in this study is very close to the target date of image fusion, with only eight days difference. The reflectance value of each band of image acquisition dates is also very close to those of the target date in 2011 (mean absolute difference of each band on two dates in 2011 is 0.001, 0.001, 0.001, 0.004, 0.004, and 0.001, respectively), indicating that the image obtained by the spatio-temporal fusion algorithm has certain reliability.
The analysis results of NDVI and DI values in the study area of each year show that the NDVI and DI curves vary considerably without eliminating the phenological influence, and there is no obvious change in trend. After eliminating the phenological influence, the changing trend of NDVI and DI values in each year is more stable, and on the whole, NDVI shows a slight upward trend, while DI shows a slight downward trend. Therefore, the elimination of phenological influence plays an important role in monitoring vegetation changes. At the same time, after removing the impact of phenology, the NDVI and DI trend curves of the study areas in each year reflect relatively consistent vegetation changes, further illustrating the reliability of the phenology elimination method and the credibility of vegetation monitoring results.
In the quantitative analysis of remote sensing, the relationship between surface property measurements at different spatial resolutions often causes concern [65]. Since vegetation cover can be highly heterogeneous spatially, subpixel variability is likely to introduce uncertainties in the vegetation indices at different resolutions [66]. Several studies have investigated the impact of spatial resolution on NDVI, but with conflicting results. Aman et al. [67] concluded that NDVI derived from the coarse spatial resolution sensor data can be used in lieu of NDVI integrated from fine spatial resolution without introducing significant errors. On the other hand, Price [68] noted that for a region consisting of a mixture of totally vegetated area and non-vegetated area, prominent discrepancies occur between NDVI derived from high-resolution measurements and NDVI derived from low resolution measurements, with the relative difference approaching 30%. This study used Landsat data with 30 m spatial resolution. In future studies, remote sensing data with different resolutions can be used to further explore the impact of eliminating phenological influences on post-fire vegetation restoration monitoring.

Conclusions
Taking the forest restoration in the Greater Hinggan Mountain area after the "5.6 fire" in 1987 as an example,based on the MODIS and Landsat time-series images acquired from 2000 to 2018, this study took the midpoint of the vegetation growth period of each year as the target date and used the STARFM fusion algorithm to construct reflectance images of vegetation with relatively consistent growth periods. The influence of phenology on vegetation monitoring was analyzed using three aspects: band characteristics, NDVI and DI values.
Based on the detailed analysis using remote sensing data, it can be concluded that eliminating phenological influences can more accurately reflect the changes of vegetation within the region, which implies that phenological factors in remote sensing images may affect the observation of vegetation changes. Observation of vegetation changes using remote sensing images of different periods of vegetation growth may cause great errors. The spatio-temporal data fusion method used in this study effectively eliminated the influence of phenological factors during the annual observation of vegetation by establishing vegetation reflection images with relatively consistent growth periods. At the same time, this method is conducive to improving the utilization of remote sensing data because researchers do not need to find remote sensing images with consistent vegetation growth conditions for monitoring but can use images located in different vegetation growth conditions and then transform them to more consistent conditions through the spatiotemporal fusion method, thereby improving the temporal resolution of vegetation monitoring. After eliminating the influence of phenology, the results based on remote sensing indices in the study area showed that although the forest in this region was affected by disturbances in some years, its growth trend is generally better. The conclusion drawn in the present analysis provides a reference for future forest monitoring research and local forest management.