Performance Evaluation of Long NDVI Timeseries from AVHRR, MODIS and Landsat Sensors over Landslide-Prone Locations in Qinghai-Tibetan Plateau

: The existence of several NDVI products in Qinghai-Tibetan Plateau (QTP) makes it challenging to identify the ideal sensor for vegetation monitoring as an important factor for landslide detection studies. A pixel-based analysis of the NDVI time series was carried out to compare the performances of ﬁve NDVI products, including ETM+, OLI, MODIS Series, and AVHRR sensors in QTP. Harmonic analysis of time series and wavelet threshold denoising were used for reconstruction and denoising of the ﬁve NDVI datasets. Each sensor performance was assessed based on the behavioral similarity between the original and denoised NDVI time series, considering the preservation of the original shape and time series values by computing correlation coefﬁcient (CC), mean absolute error (MAE), root mean square error (RMSE), and signal to noise ratio (SNR). Results indicated that the OLI slightly outperformed the other sensors in all performance metrics, especially in mosaic natural vegetation, grassland, and cropland, providing 0.973, 0.015, 0.022, and 27.220 in CC, MAE, RMSE, and SNR, respectively. AVHRR showed similar results to OLI, with the best results in the predominant type of land covers (needle-leaved, evergreen, closed to open). The MODIS series performs lower across all vegetation classes than the other sensors, which might be related to the higher number of artifacts observed in the original data. In addition to the satellite sensor comparison, the proposed analysis demonstrated the effectiveness and reliability of the implemented methodology for reconstructing and denoising different NDVI time series, indicating its suitability for long-term trend analysis of different natural land cover classes, vegetation monitoring, and change detection.


Introduction
The most extensive plateau in the world, the Qinghai-Tibetan Plateau (QTP), is covered with a broad range of vegetation from forests to shrublands [1], which are spatially distributed [2] and play an essential role in controlling the amount of transpiration and try inconsistencies (systematic and inherent noise), and ground conditions [46][47][48][49][50][51]. These lead to interruption, loss (missing), or contamination of the NDVI time series due to noise, leading to the inconsistency of NDVI products in QTP. Extensive attempts have been made to reduce noise, gaps, and determine NDVI patterns in QTP on various sensors (MODIS, AVHRR, Landsat, etc.), mainly using spectral frequency methods such as Fourier transform (FT), harmonic analysis of time series (HANTS), and wavelet threshold denoising (WTD) [45,52,53]. However, these studies were limited to specific sensors or NDVI product performance in specific vegetation types (alpine grassland as the predominant vegetation type).
Although individual sensors performance for spatiotemporal monitoring of vegetation dynamics in QTP has been extensively studied, NDVI time series analysis combined with a comparison of multiple satellite sensor modeling performances in different land covers in QTP has received less attention. Therefore, a comprehensive performance comparison between different satellites is essential to identify and illuminate the best satellite product for long-time NDVI time series analysis in the region, especially for landslide-susceptible locations, as an essential basis for the environment and natural disasters (landslide) studies. To fill the existing gap, this study aims explicitly at the pixel-based comparison of the robustness of long NDVI time series derived from five satellite sensors (AVHRR, MOD, MYD, EMT+, and OLI) in the different land covers of the landslide-prone locations in the QTP region during 2000-2020, differing from previous studies only focused on a specific part of the region and on a specific sensor. Identifying the best sensor performance offers an excellent tool for temporal vegetation monitoring and provides an essential and unique baseline for future landslide studies, especially landslide modeling (spatial and temporal) in the region, which requires accurate information from dynamic conditioning factors such as NDVI. Moreover, considering that the NDVI datasets are inherently noisy and contain missing values, we use the two state-of-the-art algorithms, including HANTS for the reconstruction (gap-filling) and the wavelet threshold denoising (WTD) for denoising (smoothing, reducing the noise interference, and data inconsistency) and modeling purposes. The performance assessment was based on the capability of the observed and denoised NDVI time series to preserve the structure of the original shape and time series values. The specific objectives of this study are: i.
Evaluation of the performances of the five NDVI datasets in different land covers as well as over the predominant land cover types in QTP; ii. Identification of the most suitable sensor for NDVI time series analysis across QTP.   NDVI significantly related to increasing temperatures, precipitation effects less pronounced. [44] 2 NPP of Alpine Grassland; effect of topography and human activity NOAA-AVHRR (10-day composite) MODIS (10-day composite)  NPP (Net Primary Productivity) showed greater decreasing trend in high-elevation regions with slope 15-30 • , aspect having little influence and roads having higher effect than residential areas. [54] 3 Grassland vegetation cover and climatic factors GIMMS NOAA-AVHRR (15-day composite) (1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999) NDVI increase in growing season from both advanced growing season and accelerated vegetation activity. Apart from temperate steppe, NDVI exhibited no significant increase in autumn, corresponded to an increase in temperature in spring, and in summer, it was related to temperature and sensitive to precipitation in the spring. There were significant lagged correlations between precipitation and NDVI for alpine grasslands (alpine meadow, alpine steppe). [55] 4 Vegetation greening and elevation SPOT-VGT (10-day maximum value composite) (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) NDVI increased more at lower elevations, but was relatively stable or even decreased at high elevations. Vegetation greening rate is in contrast to the pattern of elevation-dependent warming (EDW) with more significant temperature increases at higher elevations. Decreasing precipitation does not reverse overall increasing trend in NDVI, but it may be a limiting factor. [56] NDVI exhibits spatial variation in soil respiration better than EVI and MSAVI. At the peak growing season of alpine grasslands, soil respiration was generally much higher in the SE Tibetan Plateau and gradually decreased toward NW part.

Study Area
The Qinghai-Tibetan Plateau (QTP), with a total area of about 2,572,400 km 2 , is located in southwestern China (74 • -104 • E, 25 • -40 • N, see Figure 1) [62,63]. It is the highest (average altitude > 4000 m ASL) and largest plateau in the world with unique topography. Due to its structure complexity, the QTP is characterized by low temperatures ranging from −15 • C to 10 • C and low precipitation, with cold and arid climatic conditions in winter controlled by Siberian high and Mongolian high, and semi-arid climate in summer controlled by the South Asian monsoon system [33,64,65]. The vegetation cover in QTP mainly includes steppe, shrub, desert, meadow, forest, barren areas (bare soil), and water bodies (ice and glacier) [38,66]. The major ecoregions in QTP are composed of Karakoram-West Tibetan Plateau alpine steppe, Central Tibetan Plateau alpine steppe, Tibetan Plateau alpine shrublands and meadows, Southeast Tibet shrublands and meadows, North Tibetan Plateau-Kunlun Mountains alpine desert, and the Yarlung Tsangpo arid steppe (Olson et al. 2001). The alpine grasslands (alpine steppe and meadow) cover about 59.15% (approximately 1,521,500 km 2 ) of the total QTP [53,67]. The Himalayan regions around the southern boundary of QTP are considered the landslide hotspots and cause hundreds to thousands of casualties each year [9,11].
This study focuses on the 618 landslide-prone locations (see Figure 1) in the Himalayan region and its surrounding areas. In southern Himalayas (landslide locations), there is a shift of vegetation cover from rainforest westward to needle leaf forest, shrub, meadow, and snow (glacier) (Peng et al., 2012). A total of fourteen major vegetation classes were identified for the 618 landslide locations. The landslide locations were obtained from the NASA Global Landslide Catalog (GLC), which contains information about landslides locations and characteristics from 2007-2016 [11]. The GLC chiefly gathers landslide information from numerous datacenters, including the International Consortium on Landslides; International Landslide Centre, University of Durham; EM-DAT International Disaster Database; International Federation of Red Cross and Red Crescent Societies field reports; Reliefweb; humanitarian disaster information run by the United Nations Office for the Coordination of Humanitarian Affairs; and other online regional and national newspaper articles and media sources. The GLC has the best resolution at 0.2 • [68,69], and it has been used in different landslide studies at the global and regional scale, especially the QTP region [11,[69][70][71][72].

Data Used
In this study, five NDVI source data (AVHRR, MOD, MYD, ETM+, and OLI) covering the QTP during 2000-2020 from Google Earth Engine (GEE) were collected to compare their performances. Overall, about 10,000 satellite imageries were analyzed. All products were corrected, cloud masking was considered in calculations, and cloud pixels were removed. The NDVI is used to estimate vegetation dynamics, and it is the normalized ratio of the near-infrared (NIR) to red (RED) reflectance ρNIR ρRED : where ρNIR is the surface reflectance measured at the near-infrared channel and ρNIR refers to the surface reflectance measured in the red channel of satellite imagery. In general, higher NDVI values indicate healthier vegetation, whereas negative values (close to 0) indicate bare soil or snow [73][74][75]. A pixel-based approach was applied to compare the five NDVI datasets performance in the study area to extract NDVI values for the 618 landslide sites. The long NDVI data products and their specifications are tabulated in Table 2 and described below. The NOAA Climate Data Record (CDR) of the Advanced Very High-Resolution Radiometer (AVHRR V5) sensor has a 1-day temporal resolution and 0.05 • spatial resolution [76,77]. It includes the daily NDVI values obtained from NOAA-AVHRR surface reflectance data [75,78]. The AVHRR sensor is highly sensitive to water vapor in the atmosphere due to the broader bandwidths that affect pixel quality [42,44,79]. Therefore, the Quality Assurance (QA) products were used to select the least cloud-contaminated pixels in the dataset. QA layer provides useful information about each pixel overall cloud contamination and validates the value of the individual channel for each pixel [80]. In this study, based on the QA band information, low-quality pixels were removed, and high-quality pixels were considered for analysis.

MODIS-NDVI (MOD13Q1-MYD13Q1)
The Moderate Resolution Imaging Spectroradiometer (MODIS) vegetation index products are designed to provide a consistent spatiotemporal monitor of vegetation dynamics and have been in production continuously since 2000 [48,81,82]. This study used the MOD13Q1 and MYD13Q1 NDVI products (version 006) with 16-day and 250 m resolutions. Both products contain QA flags that hold information about pixel usefulness, which excludes bad pixels with cloud and shadows [42,83]. The usefulness index was applied so that "good" quality pixels (based on QA band) were selected, and the rest are excluded from both sensors.

Landsat Series (7 and 8)
Different Landsat series (4-8) have been providing information since 1972 to present at different spatial resolutions [84]. Landsat series have near-polar orbits, repeating every 16 days. Two latest sensors (ETM+, OLI) are used to derive long NDVI datasets at 30 m spatial resolution. We used USGS Landsat 7 and 8 Surface Reflectance Climate Data Records (LSR-CDR) in this study. LSR-CDR is processed and orthorectified [75,85] based on the National Aeronautics and Space Administration (NASA)-funded Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) program. Additionally, LSR-CDR data have been atmospherically corrected using LEDAPS and include a cloud, shadow, water, and snow mask produced using CFMASK [48,51] as well as a per-pixel saturation mask [86]. The QA layer obtained from cloud masks [48,51,87] was used to select the highest quality pixels, lowest cloud, and shadow contaminations. Thus, the 16-day NDVI dataset was produced from these corrected reflectance data using Equation (1).

Land Cover Dataset
European Space Agency (ESA) Climate Change Initiative Land Cover (CCI-LC) provides long time series of global data designed explicitly for global landcover characterization [88] and vegetation analysis [89][90][91]. ESA CCI-LC maps were published as the first series of annual 300 m (spatial resolution) global land cover maps for 24 years, from 1992 to 2015, integrating several ground observation products based on ESA GlobCover products with 71.1% accuracy [92,93]. This study used the ESA CCI-LC data (http://maps.elie.ucl.ac.be/CCI/ viewer/download.php; accessed on 18 July 2020) from 2000-2015 (38 land cover classes), and the majority of land cover types for each pixel was considered as the landcover type of the given pixel. Landcover classes were extracted for the 618 landslide locations, and those classes with low data content (Shrubland LCC 120 and evergreen shrubland LCC 121 with one observation) were excluded from the analysis. Hence, a total of twelve landcover classes were identified and retained for further analysis.

Methods
The five long-term NDVI datasets in Table 2 were preprocessed using the harmonic analysis of time series (HANTS) algorithm, specifically to fill and eliminate the gaps. They were then denoised using the wavelet threshold denoising (WTD) algorithm. Figure 2 presents the methodology adopted for the preprocessing and denoising of the five NDVI datasets.

Figure 2.
Schematic of adopted methodology for preprocessing (gap filling using harmonic analysis of time series) and denoising (using wavelet threshold denoising) of the five NDVI time series used in this study.

Harmonic Analysis of Time Series (HANTS) Algorithm
Noise, interruption, loss (missing), or contamination within the original NDVI datasets usually occur due to cloud, atmosphere, and climatic conditions that limit NDVI applications [47][48][49]51]. Thus, it is necessary to derive an effective method (gap filling and de-noising) to extract meaningful information from the original NDVI dataset. In this study, the HANTS algorithm [51] was implemented in the five NDVI datasets to identify and eliminate any sudden drops in observation points caused by different noise sources. Since the NDVI dataset has a generally irregular time interval (due to missing values), the HANT algorithm was explicitly used to fill the dataset's gaps and produce a regular time interval. The HANTS algorithm transforms the time series into the frequency domain. The new time series is then reconstructed using sine and cosine basic functions, allowing the analysis of the signal content (frequency) and generating smooth variable values [45,51,94]. The general equation of HANTS to reconstruct the original data is described as: where NDV I o (t i ) is the observed NDVI value in time t i , and NDV I H (t i ) is the reconstructed NDVI value, with the modelling error of ε(t i ); n f represents the number of frequency components f k , and the A k and θ k are the amplitude and phase of frequency component; α 0 is the average value of time series. The HANTS modeling can be obtained by solving the above equation using the linear least square (LLS) method. Invalid NDVI values (ice and snow flags) were removed to run the HANTS algorithm, and then the HANTS was performed separately on each NDVI dataset. The HANTS algorithm parametrization directly affects its performance; for example, the number of frequency components n f and degree of over-determinedness (DOD) are critical in the analysis [45,51]. The n f varies from 2-5 (three, six, nine, twelve months) in different studies [51,[95][96][97][98]. In this study, the first-four frequency components related to the phenological characteristics were used separately for each year. DOD presents the regression overfitting, which can be calculated as [45]: where MB is the minimum number of observations. After the estimation of HANTS annual coefficients (α 0 , A k , θ k for twelve months), the missing NDVI values for each year were estimated.

Gaussian-Based Wavelet Threshold Denoising (WTD)
After performing the HANTS to derive regular time intervals of NDVI datasets, the noise in them should be removed for the time series analysis. Numerous studies have utilized WTD for denoising and smoothing of NDVI time series [99][100][101][102]. The WTD-based method is also used in this study for handling noise in NDVI datasets.
Discrete wavelet transformation (DWT) decomposes time series into a set of sub-series using basic functions known as "mother wavelet", ψ(t) from dilation (scale), and translation (position) parameters related to low and high range frequency domains [94,103,104]. Given a proper mother wavelet ψ(t) and the decomposition level, NDVI data based on WT can be defined as [105,106]: in which a 0 , b 0 are constants, k is time position factor (translation factor), and j indicates the decomposition level. ψ(u) is the complex conjugate. W NDV I (j, k) is discrete wavelet coefficient of NDVI under position k and level j, including both the approximation coefficient (A j ) and the detail coefficient (D j ) at level j. The orthogonal DWT is used commonly by assigning a 0 = 2 and b 0 = 1 [107] which consists of log 2 n levels given that the analyzed series has the length of n [106,108]. Completed description of DWT can be found in [105,107,[109][110][111].
The DWT which meets the regularity condition can be adopted for NDVI time series reconstruction [106]: Wavelet threshold denoising (WTD) is based on DWT results of NDVI time series and adjusts the wavelet detail coefficients: in which TR represents the thresholding rule, T j is threshold in decomposition level j, and W NDV I (j, k) is the denoised value of W NDV I (j, k). By substituting W NDV I (j, k) into Equation (6), the denoised NDVI time series can be reconstructed. Any residuals of wavelet denoising can be further used to extract information from signal content (noise and outliers) [106,112,113]. The WTD efficiency is highly sensitive to four critical factors, including wavelet basis function, decomposition level choices, threshold estimation, and coefficient thresholding, as they influence the model accuracy [111,114]. The decomposition level is a key factor in wavelet decomposition. Several studies have used different decomposition levels for DWT [115][116][117][118]. Furthermore, Sang [105] proposed an extensive practical procedure to select the suitable decomposition level. The choice of decomposition level in this study was performed based on an empirical examination due to the complex input dataset (different NDVI time series) and the implemented algorithm (twenty decomposition levels were examined in the study). The WTD on detail coefficients was implemented in two steps: i.
Since abrupt changes caused by outliers influence the NDVI time series, it is essential to identify and eliminate these detail coefficients. The Gaussian-based thresholding method was used in this study. Hence, based on the experimental examination, the confidence interval of −σ~σ + µ (68.27%, µ is mean value and σ is standard deviation) was considered as an outlier, and those detail coefficients outside the confidence interval were removed (i.e., set as 0) [119,120]. ii. Typical wavelet thresholding method "RIGRSURE" was performed to estimate threshold T j [121] (in Equation (7) (Figure 3). The mathematical equation of the RIGRSURE threshold [122][123][124] can be defined as: where S j is the standard deviation of the wavelet coefficient at scale j; ω k represents the kth coefficient wavelet square (coefficient at minimal risk) at scale j.

Overall and Localized Model Performance
Performance assessment is essential to evaluate the reliability and effectiveness of the modeled NDVI time series with respect to the original dataset [125,126]. Because the accurate and timely ground reference data are inaccessible for validation, the systematic evaluation of remote sensing data from ground truth data is often challenging [127]. Hence, the performance evaluation and the reliability of the five NDVI datasets can be assessed by the absolute difference between the denoised NDVI time series and the original (observed) time series [51,125,126,128]. Likewise, this difference will point out which original NDVI series are less noisy considering the proposed denoising method. We used the following four performance criteria: correlation coefficient (CC), mean absolute error (MAE), root mean square error (RMSE), and signal to noise ratio (SNR). Let us consider NDV I o i and NDV I D i as the observed and denoised NDVI time series, and NDV I o i and NDV I D i their mean values, with the number N of observations; the CC, MAE, RMSE, and SNR can be derived as follows: In this study, we computed the SNR using the following equation [129]: is the variance of NDV I D i . CC is the rate of collinearity between denoised and observed NDVI time series. Higher values of CC mean higher similarity between the two-time series. RMSE indicates the rate of discrepancy between the modeled and the original NDVI time series: lower values of RMSE refer to a lower rate of dissimilarity between the two time series [130]. SNR is a well-known signal processing metric that evaluates the denoised (modeled) NDVI time series quality. Higher values of SNR signify a higher quality of denoised time series [131].
The overall model performance was provided by box-plot analysis for each sensor, considering the above metrics computed in each location (618 total). The localized performance was similarly computed for the different land cover types by clustering the pixels with similar land covers.

Results
The gap values for each NDVI dataset were estimated and filled using the HANTS algorithm in Equation (2), with an example shown in Figure 4 (panels A-J) for point 18 having "tree cover, broadleaved and evergreen" as landcover type. After that, the DWT algorithm with different mother wavelets and different decomposition levels was examined for denoising the NDVI time series (examples reported in Figure 4), and each NDVI dataset performance was evaluated using the criteria in Equations (9) and (10).

NDVI Time Series Denoising and Evaluation
To select the suitable mother wavelet, we examined different mother wavelet families, including the Daubechies (db), Symlets (sym), Coiflets (coif), BiorSplines (bior), and ReverseBior (rbio) families [132]. The AVHRR dataset was selected as the sample dataset here due to its higher temporal resolution than other datasets. The statistical analysis (boxplots) of the results obtained using different mother wavelets and different decomposition levels are shown in Figures 5 and 6, respectively. They depicted that "db9" exhibited better results (RMSE = 0.096 and CC = 0.741) compared to others. In addition, the boxplot analysis reported in Figure 6 indicated that level 1 shows comparable performance to level 2, with higher levels having worse RMSE and correlation. Thus, decomposition level 2 and "db9" were selected as the most suitable mother wavelet and the best decomposition level for WTD of NDVI datasets. Daubechies family belongs to the orthogonal wavelets with the maximum number of vanishing moments, referring to mother wavelet roots and smoothness [133,134]. From Figure 5, it can be noticed that the near-symmetric "Coiflet 5" showed better performance among the Coiflet's mother wavelets. The Coiflet wavelet has a larger number of vanishing moments in both scaling and wavelet function [135]. "Sym7" performed better among Symlet mother wavelets, which are of the modified Daubechies family with more symmetry [107]. "bior 3.5" with symmetric structure had a higher correlation among BiorSplines [136]. "rbio 2.8" showed a higher correlation than the Rbio family. Rbio is derived from biorthogonal wavelet coupled and guided by biorthogonal spline [137]. Note that more detailed information on the wavelet family refers to wavelet toolbox documentation [132].  panel (B)) values between observed and denoised NDVI time series by using different mother wavelets. The AVHRR sensor was selected as the sample sensor due to its larger dataset (2000-2020). The wavelet "db9" exhibited better RMSE (0.096) and cross-correlation (CC = 0.741) compared to others in each box. In each box, the central mark is the median, and the bottom and top box edges indicate the 25th (Q1) and 75th (Q3) percentiles, respectively; the whiskers extend to the min/max data points.
After selecting the mother wavelet "db9" and the decomposition level 2, the Gaussianbased WTD method was implemented in the detail coefficients to identify and eliminate the high-frequency detail coefficients (out of confidence interval, 68.27%, see Section 2.3.2) caused by outliers. Therefore, each NDVI dataset of the five sensors was denoised and smoothed based on regular wavelet thresholding, as in the example previously shown in Figure 4.  Figure 7.
In Figure 7A,B, the NDVI time series from MYD for the sample point located at 30.0315 • N, 101.0144 • E (China, large-mudslide) is reported. The landslide event occurred on 13 July 2011 at the maximum annual NDVI values. Visually, it did not cause a significant change in the annual phenological trend, and slight differences between the original and denoised NDVI were observed pre-and post-disaster. In Figure 7C was observed post-disaster, a disturbance that the smoothed denoised NDVI tried to follow. Overall, landslides may or may not cause disturbances in NDVI values. Therefore, while considering specific studies mainly based on the use of the NDVI anomalies for landslide detection, the true anomalies or the noisy values could be distinguished by considering a statistical analysis of the differences between the original and the denoised NDVI trends. After denoising the five NDVI time series, an overall and localized (i.e., considering each vegetation class) sensor performance was carried out according to the four criteria (CC, MAE, RMSE, and SNR), and the results are presented in Figures 8 and 9, and Table 3. Each sensor was evaluated considering the behavioral similarity between the original (observed) and the denoised NDVI time series. Specifically, the ability to preserve the original shape (CC) and the time series values (MAE, RMSE, SNR) was quantified. This analysis will also determine the sensor providing the original NDVI series with less noise on the basis of the proposed denoising approach. In Table 3 Table 3). The first-row represents the CC, the second-row shows the MAE, the third-row highlights RMSE, and the fourth-row points to the SNR values for AVHRR, EMT+, OLI, MOD, and MYD sensors, respectively.

Results of AVHRR Dataset
The AVHRR-NDVI product inherently has a noisy data structure as a major drawback in vegetation monitoring applications [42], because it was not initially designed for this purpose [79,138]. Although the QA band eliminates most of the cloud and shadow effects in the data, there are still consecutive local fluctuations with different amplitudes. Besides, the presence of annual and inter-annual patterns revealed that the HANTS could model the overall trend, and then WTD modeled these local patterns as well as local variations (Figure 4A,B). The overall performance metrics from the median quartile (CC = 0.972, MAE = 0.015, RMSE = 0.021, and SNR = 26.236) were higher than MODIS and very comparable to the Landsat 8 ones (Figure 9, and Table 3). The sensor performance ranges were 0.151, 0.013, 0.036, and 7.636 for correlation coefficient, MAE, RMSE, and SNR, respectively. Additionally, the localized performance of AVHRR showed the highest performance of the sensor (CC = 0.974, MAE = 0.014, RMSE = 0.021, and SNR = 27.452) in first-dominant land cover in QTP (needle-leaved, evergreen, closed to open LC% = 22.97%) (Figure 9, Table 3).

Results of Landsat Dataset (ETM+, OLI)
ETM+ analysis and statistical performance ( Figure 4C,D) indicated that both the HANTS and WTD algorithms underestimate the maximum annual NDVI values, which may have reduced sensor performance compared to OLI. The overall sensor performance was acceptable (CC = 0.960, MAE = 0.022, RMSE = 0.030, and SNR = 23.562) as it showed better performance in most land cover classes than the MODIS series ( Figure 8, Table 3). The maximum performance of the sensor (highest correlation coefficient and SNR, with minimum MAE, RMSE) occurred with values 0.967, 24.081, 0.019, and 0.016, respectively, in cropland and grassland classes (Figure 9, Table 3). The sensor performance ranges were 0.374, 0.018, 0.076, and 5.935 for CC, MAE, RMSE, and SNR, respectively.
OLI showed that the annual sensor behavior was relatively smooth, as the HANTS and WTD algorithms modeled the NDVI time series with a minimum degree of fluctuation (smooth local maximum) ( Figure 4E,F). The results also referred to the minimum modeling deviation from the original shape of the data, which leads to higher overall (CC = 0.973, MAE = 0.015, RMSE = 0.022, and SNR = 27.220) and localized performance in different classes (CC > 0.960, MAE ≤ 0.021, RMSE ≤ 0.028, and SNR > 25.821) (Figures 8 and 9, and Table 3). The sensor performance ranges were 0.208, 0.025, 0.0017, 0.083, and 8.290 for CC, MAE, MSE, RMSE, and SNR, respectively. Finally, OLI showed the highest localized performance in the second, third, and fourth dominant land cover classes in QTP (mosaic natural vegetation LC% = 15.37%, grassland LC% = 13.26%, and cropland LC% = 8.90%) ( Table 3).  Table 3) and the grassland and herbaceous cover for the MYD (MAE = 0.034, RMSE = 0.042, CC = 0.944, and SNR = 20.292) (Figure 9, Table 3).
Overall, OLI relatively outperformed the other sensors based on the four metrics: i.e., the original OLI NDVI is less noisy with respect to the original NDVI from the other sensors. It showed the highest localized performance in the second, third, and fourth dominant land cover classes, specifically

Discussion
Landslide events may occur at a different time of the vegetation phenology and can be detected through an NDVI anomaly and change detection techniques between preand post-landside disasters [139,140]. Previous studies have shown that the vegetation disturbance due to the landslide events has a significant impact on NDVI values in the landslide locations [139,141,142]. The presence of the noise in the NDVI time series might be misinterpreted as the anomaly caused by landslide events in NDVI. This anomaly should be statistically clarified by considering the difference between the original and the denoised NDVI time series to identify NDVI anomalies related to landslide events. Furthermore, smooth annual cyclic behaviors of NDVI time series are often influenced by environmental factors that cause defects in its intrinsic behavior, which can be addressed through mathematical and statistical modeling. Among the various methods proposed for NDVI time series analysis, the HANTS (sine and cosine basis functions) is a robust method for reconstructing the NDVI time series and estimating missing values [51,[143][144][145]. Besides, WTD has been proven to be a compelling method for time series analysis of such non-stationary data with different basic functions [111]. Due to the importance of NDVI in different applications and the existence of different datasets for monitoring NDVI and specific attributes of the study area (QTP) in terms of natural hazards (floods, landslides, etc.), reconstruction, modeling, and comparison of the five NDVI datasets were successfully implemented in different land cover classes in this study.
Although the AVHRR has a higher temporal resolution, its inherent characteristics, such as design, significant water vapor absorption due to broader bandwidth than the other sensors, and low spatial resolution, have led to exhibit noisy behavior in the dataset, especially in the QTP region [42,44,61,73,78,146,147]. To overcome this, the use of cloud filters and shadow masks in the QA band can eliminate large numbers of poor-quality pixels and improve the annual cycle behavior of the NDVI time series [148]. However, the QAband-filtered data still contained noise in sequential fluctuations, identified and modeled in a satisfactory way by using HANTS and WTD ( Figure 4A,B). The sensor showed the best performance in the first dominant vegetation class (needle-leaved, evergreen, closed to open; LC = 22.97%) in the QTP. MODIS (MOD, MYD) is commonly used as reference data for long NDVI time series analysis [42], and both MOD and MYD use a similar NDVI compositing algorithm [83,149]. The main difference between the two sensors is their equator crossing time (10:35 am for MOD, 1:35 pm for MYD) [148]. Their localized performance and the range of variability (MOD:CC = 0.148, MAE = 0.034, RMSE = 0.048, SNR = 7.006, and MYD:CC = 0.154, MAE = 0.038, RMSE = 0.052, SNR = 7.139) were almost the same, with slightly better performance of the MOD sensor in most classes, consistent with previous findings [148]. The maximum difference between these two sensors was observed in tree cover, broadleaved, evergreen, and the least difference in mosaic cropland. However, the results showed a relatively higher MAE, RMSE, with lower CC, and SNR values for MODIS compared to the other sensors, which might be related to the higher number of artifacts (abrupt shifts) observed in the original data, that HANTS and WTD handle quite well.
The OLI-NDVI time series were characterized by gaps and cloud contamination (about 35% of imageries) [150,151], which ultimately performed better than the MODIS (MOD, MYD) and AVHRR datasets. Although the cloud masking algorithm (CFMASK) [48,51] created gaps in the dataset, the HANTS and WTD algorithms were able to preserve the original behavior and model the NDVI time series well (see Figures 8 and 9). The overall and localized performance of OLI and AVHRR were higher than the other sensors in all vegetation classes, which indicates the consistency of the NDVI time series, especially in the OLI dataset. As mentioned above, the better performance of the OLI sensor indicates that the original NDVI time series obtained from Landsat 8 can show the interannual and cyclic behavior of the NDVI dataset with a lower amplitude of the noisy fluctuations compared to the other sensors.
Additionally, the overall performance variability in OLI showed a significant agreement with ETM+ with 0.013, 0.007, 0.008, and 3.658 differences in CC, MAE, RMSE, and SNR values, which is consistent with previous studies and shows that two sensors can be used as complementary data [152]. However, as the newest satellite in the Landsat series, OLI performs better than ETM+ due to the main differences in characteristics such as time span, spectral bandwidth (narrower), radiometric resolution, calibration, signal-to-noise ratio [35,47,48,[152][153][154].
In summary, the result shows that the original NDVI datasets are inherently noisy with missing values and outliers. To reduce the inconsistencies of these native NDVI datasets, reconstruction (gap-filling), smoothing, and modeling of the NDVI time series are required in the vegetation studies. The integration of HANTS for reconstruction and WTD for denoising and modeling NDVI datasets is critical for monitoring long-term vegetation dynamics, especially in mountainous zones and landslide-prone areas with sparse vegetation covers. Although different NDVI time series have different spatiotemporal resolutions, in this study, our primary focus was on reconstructing (gap-filling) and denoising the NDVI time series based on the temporal resolution. The pixel-based performance behavior between the denoised and the original NDVI time series was evaluated for each sensor separately, but no intercomparison was done between different NDVI values among different sensors. Hence, the spectral and spatial resolution difference between different sensors has not been considered as a constraint for this study. In this context, the better performance in the temporal modeling of the OLI sensor and its better spatial resolution indicate that this sensor is the best option for vegetation monitoring in the QTP area.
The main limitation of the implemented algorithms is the high sensitivity to the types of landcovers. For example, bare soils with very low NDVI content cannot reflect the behavior of harmonic-cyclic annual and inter-annual NDVI datasets. Therefore, HANTS and WTD may not be able to reconstruct and model time series of such land cover successfully.

Conclusions
Vegetation as a major dynamic conditioning factor plays an important role in landslide detection by investigating the temporal variation in vegetation using NDVI time series. The existence of several NDVI products often makes it challenging to identify the ideal sensor for vegetation monitoring related to landslide studies. Identifying the best-performing sensor is essential for vegetation monitoring, and provides a unique baseline for future landslide studies, especially landslide modeling (spatial and temporal) in the region, which requires accurate information from various landslide conditioning factors (LCFs), including the NDVI as a prominent dynamic factor.
In such cases, application of advanced noise removal techniques, and comparisons between different NDVI datasets, are essential to ascertain the best sensor performance for spatiotemporal NDVI analysis. Hence, a pixel-based analysis of NDVI time series at landslide-prone locations, detecting and modelling the residual effect of noise on atmospherically and environmentally corrected data, was designed in this study to compare the performance of five popular (ETM+, OLI, MOD, MYD, and AVHRR sensors) long-term NDVI datasets. We used the HANTS and WTD algorithms to reconstruct and denoise the five NDVI datasets over 20 years at 618 landslide-prone locations in QTP. The performance assessment was based on the behavioral similarity between the original (observed) and denoised NDVI time series, considering the preservation of the original shape (CC) and time series values (MAE, RMSE, SNR).
The analysis showed that OLI outperformed the other sensors: that is, the native OLI NDVI series are less noisy considering the proposed denoising method. Although statistical and behavioral similarities were found between the Landsat series, MODIS series, and AVHRR-Landsat series, OLI showed the highest localized results in the second, third, and fourth dominant land cover classes in QTP.
OLI showed a better performance in the temporal modeling (gap-filling and denoising), and hence, it can be considered as a better option for long-term NDVI analysis in the QTP region. However, further investigation is recommended to assess the sensor capability for modeling the NDVI behavior in different vegetation classes in other areas. In addition, since the AVHRR NDVI dataset has shown results comparable to OLI and it has a higher temporal resolution, downscaling the AVHRR dataset into finer spatial resolutions can be a viable alternative for NDVI analysis in QTP. By fusing the AVHRR and OLI data, more reliable results of landslide studies in QTP could be obtained in the future. Overall, this study showed that the native NDVI dataset is inherently noisy, with missing values and outliers. The implementation of HANTS for gap-filling and WTD for denoising is fundamental for handling reliable and complete long-term NDVI time series, especially in landslide prone areas. Due to its effectiveness for reconstruction and denoising different NDVI time series, HANTS and WTD can be used in more recent NDVI time series such as Sentinel series, or even Landsat 9 in future, with an increase in the time span, larger number of samples, and a larger number of land cover classes for trend analysis, vegetation monitoring, and change detection at global and regional scales. Data Availability Statement: Data generated of long NDVI time series from AVHRR, MODIS, and Landsat sensors are mainly provided in the form of tables and figures in the content of manuscript documents. Further NDVI time series analysis data and codes are available upon request from the corresponding author. Additionally, NOAA-AVHRR-NDVI-V5 is available at (http://ltdr.nascom.nasa.gov; accessed on 10 July 2020 from GEE). MYD13Q1 and MOD13Q1 are available at (https://lpdaac.usgs.gov/products/myd13q1v006/, and https://lpdaac.usgs.gov/ products/mod13q1v006/; accessed on 10 July 2020 from GEE). The Landsat Series are available at (http://landsat.usgs.gov/CDR_LSR, and http://landsat.usgs.gov/CDR_LSR (accessed on 10 July 2020 from GEE). The landslide locations were obtained from the Global Landslide Catalog (GLC) (accessed on 1 March 2020).

Conflicts of Interest:
The authors declare no conflict of interest.