Satellite observations of the terrestrial biosphere cover a period of time sufficiently extended to allow the calculation of a reliable climatology. The latter is particularly relevant for studies of land surface response to climate variability. Observations from space of the land surface are hampered by clouds at shorter wavelength and affected by water in the atmosphere in the microwave range. Both polar orbiting and geostationary satellites have a revisit frequency high enough to allow for some redundancy relative to the processes being observed, so that time series where a fraction of observations are removed and the resulting gaps filled are still very useful to monitor land surface processes.

This paper reviews studies done by the authors since the late 1980s on the use of time series analysis techniques to extract concise information from extended time series of multispectral satellite data capturing large areas. Two basic methods have been used: the Fast Fourier Transform, especially in the earlier studies, and Harmonic Analysis in more recent work. Since the first studies, work has been relying on the global radiometric data collected by AVHRR and later by MODIS, FY-2C, AMSR-E, SSM/I and SPOT 4 and SPOT 5, with investigations performed in different continents. It should be noted that available satellite data records are approaching 40 years, since e.g. microwave radiometry data are available since 1978(Nimbus-7/SMMR)and multi-spectral VNIR radiometry has an even longer history(i.e. NOAA-2/VHRR since 1972). These climate data records pave the way for a truly climatological approach to the observation of land surface signals, with the additional strength of the spatial detail provided by imaging sensor systems. Methods for the analysis of time series of image-data appear perfectly suited to extract information from current data assets and build new knowledge on the response of the land surface to weather and climate forcing.

The applications supported by published results are:

(1)The identification and mapping of zones characterized by a similar response of terrestrial vegetation to environmental forcing(Azzali & Menenti, 1999, 2000; González-Loyarte & Menenti, 2008; González-Loyarte et al., 2008);

(2)The determination and characterization of terrestrial vegetation response to climate variability over any period of time covered by available time series of satellite data(Roerink et al., 2000, 2003);

(3)Early warning of anomalies detected during the growth of terrestrial vegetation by using indicators of photosynthetic activity such as Normalized Difference Vegetation Index(NDVI)and Fraction of Absorbed Photosynthetically Active Radiation(fAPAR)(Jia et al., 2012; Lu et al., 2015);

(4)Determination of the time-lag in the response of vegetation growth to precipitation(van Hoek, 2016; van Hoek et al., 2016);

(5)Determination of the phenological response to water availability in relation with water management in a watershed(Jia et al., 2011);

(6)Monitoring of the spectral thermal admittance of the land surface(Ghafarian et al., 2012; Ghafarian, 2015);

(7)Estimation and mapping of air temperature(Alfieri et al., 2013);

(8)Removal of clouds and shadows and reconstruction of multi-spectral images at high spatial resolution(Alfieri et al., 2016);

(9)Monitoring of the thermal load to assess the risk of forest fires(Maffei et al., 2013);

(10)Monitoring of actual evapotranspiration(Jia et al., 2009);

(11)Retrieval of the fractional abundance of water-saturated soil(Shang et al., 2015).

Most of the work mentioned in this study has been reviewed in more detail by Menenti et al.(2010, 2016). Here we aim at providing a concise summary of new knowledge and applications, without entering in the details of the methods developed and applied over time, but conveying a sense of the progress achieved when taking all our past studies together.

The Fast Fourier Transform(FFT)and the Harmonic ANalysis of Time Series algorithm(HANTS)have been developed and applied to time series of satellite observations, e.g. NDVI and Land Surface Temperature(LST), to study vegetation phenology and land surface climate(Menenti et al., 1993, 2010; Verhoef, 1996; Roerink et al., 2000, 2003; Julien et al., 2006; Jia et al., 2011; Moody & Johnson, 2001; Alfieri et al., 2013; Ghafarian et al., 2012). Different from the FFT using all observations regardless of quality, HANTS identifies and removes outliers in data samples.

A global study of the accuracy of HANTS in the reconstruction of NDVI time series has been completed(Zhou et al., 2015). The overall reconstruction error was divided into gaps related error and fitting method related error, then the gaps and fitting-method related errors were quantified independently. The method was applied by Zhou et al.(2015)to map globally the spatial pattern of HANTS performance in the reconstruction of NDVI time series. We have also applied other methods(Ghafarian et al., 2012; Ghafarian, 2015), such as the Multichannel Singular Spectrum Analysis(M-SSA, Broomhead et al., 1987), but we will not describe this method here.

It should be noted that the same method(HANTS, M-SSA or possibly other methods for time-series analysis)can be applied in many different ways to extract fundamentally different information from the signal captured by a time series of observations. Many land surface processes yield periodic signals, resulting from the combination of periodic or quasi-periodic forcing at multiple frequencies. Fourier-based methods to analyze time series of observations of land surface properties offer multiple benefits:(1)the decomposition of a complex signal into its periodic components provides useful insights into both weather and climate forcing and the response of the land surface(Menenti et al., 1993); (2)a model of the observed signal which can be used to identify and remove observations contaminated by clouds and shadows and to reconstruct the time series after removal of outliers(Menenti et al., 1993). The main difference in the nature of the information extracted is related to the length of the time series being analyzed. When a data record of climatological relevance is analyzed, the dynamics of the observed target(pixel)is captured and decomposed into periodic components. For all components with period shorter than the total length of the data record, the mean Fourier coefficients are obtained. If, for example, HANTS is applied to time series of NDVI, the Fourier coefficients provide a metric of phenology and can be used(see e.g. Azzali & Menenti, 2000)to map soil -vegetation -climate complexes, i.e. areas having similar phenology in response to soil conditions and climate. At the other hand of the temporal spectrum, the analysis of short time series, in a near-real time mode, reveals deviations from the long term climatology and can be used to construct indicators for drought early warning(see e.g. Jia et al., 2012; Lu et al., 2015).

In all applications listed above we applied the HANTS algorithm, developed and described in detail by Menenti et al.(1993), Verhoef(1996), Verhoef et al.(1996)and Roerink et al.(2000). The first application of Fourier series to analyze time series of satellite data was described by Menenti(1984)to obtain a constant local solar time in a time series of METEOSAT LST images stretching from Saudi Arabia to Morocco.

The HANTS algorithm is based on the concept of discrete Fourier transform(Menenti et al., 1993; Verhoef, 1996; Verhoef et al., 1996; Roerink et al., 2000)to model time series of satellite data.

A temporal sequence of N observations y_i, i=1 to N can be described by a Fourier series as

$ {{\rm{y}}_i} = {a_0} + \sum\limits_{j - 1}^M {{a_j}} {\rm{cos}}\left({{\omega _j}{t_i} - {\varphi _j}} \right) $

(1)

where ω_j is the frequency of j^th harmonic term in the Fourier series, t_i is the time at which the i^th sample was taken, M is the number of frequencies of the Fourier series(M≤N), α_j and φ_j are the amplitude and the phase of the j^th harmonic term respectively. Because the zero frequency has no phase, the amplitude related to the zero frequency, a₀, is equal to the average of all N observations of y. The harmonic frequencies are a base frequency(i.e. ω₁=2 $ {\rm{\pi }}$/N)and all integer(i.e. i=1 to N)multiples of the base frequency:

$ {\omega _i} =(2\pi /N)\times i,\;\;\;{\rm{where}}\;\quad \quad i = 1,2,\; \ldots,N $

(2)

In the HANTS algorithm after selecting the number of frequencies(M)and the frequencies(ω_j), the unknown parameters of the Fourier series are the amplitudes a_j and the phase φ_j values, which are determined by fitting the time series of observations for each pixel, i.e. generating maps of a_j and φ_j.

From a signal processing point of view, HANTS performs two tasks:(1)identifying and removing outliers and cloudy observations, and(2)gap-filling of the remaining valid observations by temporal interpolation. HANTS handles the Fourier analysis as a curve fitting problem using a weighted least squares solution. Each observation in the time series is assigned a weight of one or zero for good and bad(outliers)data, respectively. In order to find outliers and cloud contaminated observations, HANTS performs the curve fitting iteratively. In the first step, the least square curve fitting is performed using all data in the series. In the second step, observations are compared to the curve determined in the 1st iteration. Observations that deviate more than a pre-defined threshold are removed by assigning a zero weight to them. The remaining data are used to compute the least square curve fitting again, and the outliers are identified and removed again using the same threshold as in the first step. This iteration procedure is repeated until either all the remaining observations are within the pre-defined threshold or the number of remaining data becomes less than the number of parameters by which the curve is described. The accuracy of both the estimated Fourier coefficients and of the reconstructed signal depend on the appropriate setting of several parameters. A detailed explanation on HANTS parameter setting can be found in Roerink et al.(2000).

The applications listed in the introduction can be grouped into fewer and broader categories to clarify the bio-geophysical relevance of the information extracted from the observed signals.

Under this category we group the study by Azzali & Menenti(1996, 1999, 2000), who mapped soil -vegetation -climate complexes in Southern Africa(Fig. 1)using AVHRR/NDVI data over the period 1981—1991 and by González-Loyarte & Menenti(2008), González-Loyarte et al.(2008)who mapped plant communities of shrub lands and open forests species in Western Argentina, also using AVHRR/NDVI data(1982—1991). Both studies showed that the zones mapped as having similar temporal signa-tures correlated with isolines of the Budyko radiational index of dryness(Budyko, 1958), thus confirming the response of the phenological signatures captured by the Fourier spectra to climate forcing.

The value of decomposing a complex temporal signal into its periodic components was documented by Ghafarian et al.(2012)and Ghafarian(2015)who used three years of hourly observations of the Land Surface Temperature(LST)retrieved from Fengyun-2C(FY-2C)/Visible and Infrared Spin Scan Radiometer(S-VISSR)data. A three year time series(2008—2010)of gap-free daily and hourly LST and actual evaporation and soil heat flux derived from geostationary data collected by the FY-2C satellite was reconstructed for a large area including the Qinghai -Tibetan Plateau and the surrounding river basins(Faivre, 2014; Ghafarian, 2015). Two time series were used in this study: soil heat flux(Faivre, 2014)at the land -atmosphere interface(forcing)and land surface temperature(response). The Fourier spectra of both series were calculated, then the land surface thermal admittance(Menenti, 1984)was calculated for all significant periodic components. The thermal admittance is the ratio of the amplitude of a periodic component of the heat flux to the amplitude of LST at the same frequency. If the observed signal includes multiple periodic components, the thermal admittance can be calculated for all components to obtain the land surface spectral thermal admittance. Since periodic components of increasing period penetrate to increasing depths into the soil, the spectral thermal admittance carries information on the vertical structure of the soil. Ghafarian et al.(2012)and Ghafarian(2015)used the amplitudes of the periodic components of LST(t)to map areas having similar thermal behavior. They also inverted the spectral thermal admittance to estimate the thermal properties and thickness of soil layers.

In the case of complex terrains, local spatial variability of climate is not negligible and depends on several factors, among them topography, soil humidity and land use. This leads to the need for detailed spatial patterns of air temperature T_a. The sparse distribution of stations, as well as the frequently insufficient temporal coverage, severely constrains the study of climate-related processes in complex landscapes. Alfieri et al.(2013)developed a new procedure to map time series of air temperature(T_a)at fine spatial resolution using time series analysis of satellite-derived LST observations, specifically all the observations captured daily in the period 2001—2010 by MODIS onboard both the Terra and Aqua satellites. The method assumes that air temperature is known at a single(reference)location, such as in gridded climate data with grid size of the order of 35 km×35 km. Alfieri et al.(2013)modelled the LST spatial and temporal pattern within each grid cell by the pixel-wise ratios r(x, y, t)of the LST at any location to the LST at a reference location. The analysis of these patterns over a decade demonstrated that their intra-annual variability is not negligible, with significant seasonality, even if it is stable throughout the years. The estimates of T_a were evaluated against observations at five weather stations within the 35 km×35 km grid and for the same period 2010—2011. The method has also been evaluated against T_a observations for earlier periods of time(1984—1988), although available data were rather sparse in space and time, slightly larger deviations were obtained.

(1)Europe and North Africa. Work done towards mapping the “soil -vegetation -climate complexes”(Fig. 1)has documented a close relationship between vegetation type, including phenology, and aridity, which was measured by using the Budyko ratio(B). The next question is whether a similar relationship exists between temporal patterns or, in other words, whether the response of vegetation phenology, measured by time series of NDVI(or fAPAR)observations, can be related to climate forcing, measured by time series of Budyko ratio.This relationship is confirmed by the strong correlation of both mean NDVI and the amplitude of the one year component with the value of B(Fig. 2)documented by Roerink et al.(2000, 2003). These results(Fig. 2)were applied(see also Menenti et al., 2010)to estimate the sensitivity of foliar phenology to water availability.

Given the sensitivity of foliar phenology documented above, a step further led to develop metrics of anomalies in photosynthetic activity and LST towards the early detection of drought. Two new drought indicators based on satellite observations of vegetation index and LST, i.e. the Normalized Temperature Anomaly Index(NTAI)and the Normalized Vegetation Anomaly Index(NVAI)were developed by Jia et al.(2012)and applied to monitor drought events in different regions in China and India. The sensitivity of the indicators and of their anomalies to drought conditions and severity was also evaluated against drought assessments by operational drought monitoring services. NTAI and NVAI time series data were extracted for the drought-prone regions in China and India to evaluate the performance of these indicators in detecting and monitoring drought events. The data used by Jia et al.(2012)consisted of 30 MODIS tiles covering most of China and India for every 16-day intervals from 2001 to 2010.

The sensitivity of NTAI and NVAI to drought was documented(Jia et al., 2012)by e.g. the analysis of a drought event in 2007 in Inner-Mongolia(close to Xilinhot). The year started with a slightly negative anomaly in precipitation, above average precipitation in March and persistent below average precipitation from April through July. The NDVI was slightly below average(Fig. 3), while both NTAI and NVAI showed large anomalies, with a fast response to anomalies in precipitation(Fig. 3).

(2)Italy -Campania. The rapid response of thermal anomaly to drought suggested to explore monitoring of thermal anomalies as a predictor and an indicator of forest fire hazard(Maffei et al., 2013). Time series of MODIS daily retrievals of LST during 2000—2006 were used to characterize LST anomalies and thermal load. As a reference, data on about 4400 actual forest fires occurred in the Campania region during the same period were applied. Pixel-wise Fourier coefficients estimated with HANTS, i.e.yearly mean LST value, amplitude and phase of each harmonic were evaluated as predictors of fire occurrence and fire size. The daily LST time series were also used to calculate time series of anomalies and cumulated anomalies. The daily anomaly was defined as the difference between the current LST value at a given pixel and the multi-annual mean LST value at the same pixel given by the Fourier series determined over the entire 2000—2006 period.

The Cumulated Thermal Anomaly(CTA)was a better predictor of fire size(Maffei et al., 2013). The CTA is a measure of heat "accumulated" in a certain area, providing more direct information on the prolonged exposure of vegetation to stress conditions. This is reflected in the prediction of expected mean fire size over two orders of magnitude(Fig. 4), leading to more meaningful fire hazard maps as compared to the ones based on the Thermal Anomaly(TA). CTA is also a good predictor of the proportion of large fires across a larger size range(Maffei et al., 2013).

The results reviewed above document the sensitivity of vegetation to water availability, but raise multiple questions related to how and with which delay does terrestrial vegetation respond to water availability. The time lag between anomalies in precipitation and vegetation activity plays a critical role in early drought detection as crop yield losses are caused by precipitation shortages, at least in rainfed agriculture. van Hoek(2016)and van Hoek et al.(2016)explored a new approach to estimate the time lag between a forcing(precipitation)and a response(NDVI)signal by applying cross-spectral analysis in the frequency domain. The cross-spectral analysis was applied pixel-wise, with the power-spectrum applied as a metric of the cause -effect relationship, to determine the phase-lag for all significant and coherent periodic components. The study by van Hoek(2016)and van Hoek et al.(2016)focused on a dry and a wet year and only the phase lag of the annual component of the forcing and response signal was determined and interpreted, since this was the dominant component. In the wet year with above-average precipitation, the phase lag was rather similar for all land cover types. In the dry year, the phase lag increased by 7.0 days on average, but with specific response signals for the different land cover types(Fig. 5). van Hoek(2016)and van Hoek et al.(2016)interpreted the phase lag against the soil water-holding capacity and observed a slightly higher phase lag in the dry year for soils with a higher water-holding capacity.

To assess water availability and monitor its evolution, the evaporation and transpiration within heterogeneous terrestrial landscapes must be known. Jia et al.(2009)used daily MODIS data and a surface energy balance algorithm to estimate actual ET in the delta of the Yellow River, China. Because of clouds, this approach gave pixel-wise estimates of actual ET but irregularly spaced in time. Moreover, some estimates were based on MODIS observations not identified as cloudy, thus rejected, to yield a rather noisy ET time series(Fig. 6). HANTS was applied to identify and remove outliers, then fill the gaps. The gap-filled time series of actual ET for a "reed -swamp" pixel with unlimited water supply was in good agreement with reference ET. This study also demonstrated that it was not necessary to identify outliers in the time series of input variables(NDVI, LST and albedo)prior to the calculation of ET.

To observe, understand and model the terrestrial water cycle, knowledge of soil water storage is critical. Observations by microwave radiometers are applied to retrieve soil water content, snow water equivalent and water extent, and atmospheric effects are much smaller than at shorter wave-lengths. Shang et al.(2015)applied observations of the Polarization Difference Brightness Temperature(PDBT)at 37 GHz to retrieve the water extent even under dense vegetation. Due to the revisit frequency and swath width of microwave radiometers on polar orbiting satellites, time series of passive microwave data will have gaps of 3—5 days. To generate a consistent, gap-free daily time series, filtering or smoothing algorithms need to be applied. In addition, the effects of rain drops on microwave ground surface observations at 37 GHz are not negligible. Geo-location and resampling introduce random errors in the conversion from swath Brightness Temperature(BT)to the gridded data. Radiometric noise, such as the reference warm-up, also has an additional and unpredictable impact on observed BT. Erroneous observations need therefore to be identified and removed. In this case, we applied the HANTS algorithm described earlier(Menenti et al., 1993; Roerink et al., 2000; Verhoef, 1996)to remove noisy samples and to fill gaps(Fig. 7). The result was a much cleaner signal that was successfully applied to retrieve the fractional abundance of water-saturated soil.

The decomposition of a complex signal into simpler periodic components makes it possible to understand the response of vegetation to different processes determining water availability. Jia et al.(2011)studied the phenological response of vegetation in the arid lower reach of the Heihe River basin in China, where water availability in spring and summer is determined by:(1)snowmelt in the upper reaches, (2)rainfall in the upper and middle reaches and(3)water diversions for irrigation in the middle and lower reaches.

Jia et al.(2011)evaluated whether the South-North(S-N)propagation of the fluctuations in stream flow through the drainage channels and the shallow aquifer system had any influence on the timing of peak green-up, by calculating the mean phase value in the East-West direction, i.e. over all columns of the 12-month phase image for each row, and plotted such mean value in the S-N direction, i.e. along the rows of the same image(Fig. 8).

If we take the entire drainage system into account, we will have a combination of rapid propagation downstream along the drainage channels and slower outwards propagation from each drainage channel into the shallow aquifer. This mechanism is likely to explain the absence of time lag in peak green-up during the wetter years and the fact that a time lag does appear only in the northern, drier portion of the Ejin Oasis and only in drier years.

The review presented here illustrates the results obtained over about 30 years by applying time series analysis techniques to extract concise information from extended time series of large area multispectral satellite data.

The materials summarized above show that observations from space of surface reflectance in the visible, near-infrared and shortwave infrared and of emittance in both the thermal infrared and the microwave spectral range are affected by gaps in the data records and by atmospheric effects, namely clouds in the thermal infrared and precipitation at 37 GHz. This notwithstanding, such observations are still useful, as demonstrated by the applications described in this review. The combination of gaps and outliers in the data record requires additional processing to generate consistent time series of at-surface observations and two aspects need to be mentioned. First, the data record must be somewhat redundant with respect to the process to be observed, i.e. the sampling interval must be shorter than the inherent temporal scale of variability of observed variables and processes. If this is the case, methods like HANTS can be applied to model the observed time series with fewer and irregularly spaced valid observations of the land surface.

The results summarized and reviewed here open a two-fold perspective for further research:

(1)The abundance of space-borne observations of the global land surface which span 30 years or more, has a climatological relevance far superior to what was achievable with the data available 25 years ago when this research began. Therefore, this invaluable source of information deserves the most urgent attention, by a deeper and more robust understanding of the observed response of the terrestrial biosphere to climate variability and, specifically, to inter-annual variability;

(2)Both methods and data have become ripe enough to develop further an efficient and cost-effective early warning system for drought events, based directly on observations of vegetation conditions from space.