ABSTRACT

Land cover changes (LCC) are an important component of Global Change. LCC can be described not only by its occurrence, but also by the land cover replacement, causal agent and change duration or recuperation. Nowadays, remote sensing offers the opportunity to assemble reliable time series, however this fails to make a characterization of LCC since the series represents dynamics due to the combination of several processes occurring simultaneously. In this article we proposed an approach to the study of LCC using wavelet transform (WT) and MODIS vegetation time series. Through this work we have demonstrated the capacity of this tool in order to recognize and characterize four different LLC documented in scientific publications, presenting the results divided in frequency scales as interannual, seasonal and rapid changes. The information decomposed in frequency allows the interpretation of each involved process without the interference of others. The uses of WT in an image time series give us the possibility of joining temporal and spatial dimension in a single raster. Layers generated with WT might be used to pattern recognition in LCC and to improve an image classification.

1. Introduction

All Earth’s ecosystems are in a continuous state of change. Transitions from one state to another are originated by natural and/or anthropogenic causes and occurs at several scales of time and space [1-3]. Vegetation transitions are usually catalogued as vegetation phenology when they are driven by the photosynthetic activity through seasons [4], or Land Cover Change (LCC) when the variation of the vegetation is slight or abrupt [1,3]. LCC can be described not only by its occurrence but also by the land cover replacement (e.g. forest to bare land: deforestation, forest to agriculture: agriculturization, grassland to forest: afforestation); causal agent (e.g. cropping, harvesting, natural fires, floodings); and change duration or recuperation. In every case, LCC monitoring is important to understand the vegetation dynamics [4], and, in consequence, to establish relations between policy decisions, regulatory actions or subsequent land use activities [2].

Remote sensing (RS) offers a unique opportunity to collect valuable information to study and understand LCC [1]. Satellite sensor platforms have shown the ability to cover large geographic regions [2], providing multispectral data, with a revisit time ranged from minutes to months, at different spatial resolutions, and available on line for free. According to these intrinsic characteristics, RS have shown its capacity to detect LCC and assemble reliable time series [5-7]. Despite the achievement of regular data, characterizing LCC using RS data is still complex because of the combination of several processes occurring at the same time: seasonal changes, abrupt changes, climate alterations and acquisition errors [8]. The Normalized Difference Vegetation Index (NDVI), calculated from near infrared and red reflectances, was described firstly by Rouse et al. [9] and is widely the most used vegetation index. This index is generally recognized as a good indicator of vegetation activity [6, 10-12] and has been related with several biophysical variables such as: fraction of absorbed photosynthetically active radiation [13,14], leaf area index [14,15], primary production [16,17], among others. Also, NDVI is especially useful in multi temporal datasets because they permit to describe vegetation phenology [5,12,18,19].

In order to detect the occurrence of LCC and to characterize their dynamics, the methodology applied must indentify variability with one temporal scale, while identifying changes with another scale. This means that the methodology must separate changes in ecosystems into three frequency components: seasonal changes (e.g. succession, seasonal precipitation), gradual changes (e.g. interannual temperatures tendencies, land degradation) and abrupt changes (e.g. fire events, flooding, farming) [7]. Signal digital processing offers a wide variety of procedures to detect changes in remotely sensed time series and in monitor vegetation dynamics. The techniques more commonly used are statistical approaches such as principal component analysis [6,11], curve fitting [4,5,20], frequency transformation techniques such as Fourier [12,21], and most recently, Wavelet transform (WT) which has been successfully used to characterize vegetation dynamics [3,22] and expansion-intensification of agriculture [23].

Wavelet transform (WT) is a powerful mathematical tool that has been used since the beginning of the 1980s, but still does not have the diffusion of its precursor: the Fourier transform [24]. The principal advantage of WT is that it has not an unique basis, it is based on a family of functions [25]. As a result, WT has the capacity to disaggregate a signal into its component frequencies (gradual change, seasonal change, etc.), and afterwards describes how each component evolves without interference from the others [26]. Thus, this kind of decomposition introduces the possibility of studying time and frequency domains simultaneously [3,22,27].

The objective of this article is to investigate the potential of time series analysis based on wavelet transform to characterize different LCC at several frequency scales (interannual changes, seasonal changes and rapid changes). We have developed an algorithm capable of applying this transformation to data obtained from MODIS sensor products, allowing us to evaluate a complementary approach to four case studies. These cases were taken from recent publications and they encompass diverse LCC such as fire events, inundations, farming and urbanization from different locations of the globe.

2. Materials and Methods

2.1. Study Sites

In order to embrace different LCC around the globe, we have selected four papers published in the recent years, where localization, region description and vegetation dynamics studies were well documented and published (Table 1).

The four cases encompass: 1) a seasonal inundation in the Mekong Delta where the Tonle Sap Lake increases its area 3 - 6 times in wet seasons [28]; 2) a fire event occurred in August 2006 in the northwest of Spain where almost 930 km² were burned [29]; 3) a deforestation process n Chaco, Argentina, related to soybean and planted pasture expansion [30]; and 4) a rural to urban conversion in Xi’an region, China [31]. For more details of the studied regions, the mentioned article in each case is recommended (Table 1).

2.2. Satellite Data

In this study, the NDVI time-series have been obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS). MODIS product, MOD13, provides a 16 day NDVI composite, minimizing errors from cloudy days, viewing geometry and atmospheric aerosol presence [19]. Two sources of this product have been used: MODIS Global Subset page (http://daac.ornl. gov/MODIS/modis.shtml) to obtain single pixel ASCII time series; and MODIS WIST (https://wist.echo.nasa.gov/api/) to download complete images in HDF format to obtain time series images (Table 2). In all cases, time-series started on the 18th of February of 2000 and ended on the 26th of June of 2010, spanning 239 measures. We also used two high resolution images from LANDSAT-TM corresponding to the 10th of February of 2001 and the 16th of February of 2009. These LANDSAT-TM images were downloaded from INPE website (http://www.inpe. br/).

2.3. Wavelet Transform Concepts

Due to the complex mathematical description of WT, a brief explanation is presented. For further theory and formulae, these publications are recommended: [3,25-27, 32]. The WT is a mathematical tool which decomposes a signal using a family of functions based in a mother

wavelet [32]:

(1)

where is the family of functions based in dilations and shifts of a mother wavelet; is a coefficient that regulates the dilation and adjusts the location. Once the mother wavelet is chosen the WT calculate wavelets coefficients using the following expression:

(2)

where the wavelet coefficients of a function are represented by. For each scale and dilation a coefficient is calculated. In those cases where the spectral information of is similar to the wavelet the coefficient will have higher values [3].

2.3.1. Discrete Wavelet Transform

Discrete wavelet transform (DWT) is a computing efficient implementation of the WT, where scales are based on powers of two (2ⁿ). DWT splits a signal into details (D) and approximations (A), using a high frequency pass (HP) and a low pass (LP) frequency pass filters, respectively. Afterwards DWT coefficients are obtained applying a down-sampling with a 2:1 relation to D and A (Figure 1).

High-pass transfer function is determined by the mother wavelet function ψ, which has to satisfy the admissibility condition and has unitary energy [27]. The low-pass filter is determined by a so called scaling function φ, which is related to some mother wavelets. The selection of a mother wavelet is a difficult task because none of these wavelets has all its properties optimized [26]. Meyer orthogonal discrete wavelet (Figure 2) has been chosen due to its infinitely regularity, its symmetry, and because its mother wavelet permits an exact reconstruction of the original function [33].

2.3.2. Multi-Resolution Analysis

The idea of a multi-resolution analysis (MRA) is the iterative application of the DWT to decompose a signal into several frequencies scales. If we apply the DWT “m” times, the signal f is decomposed into “m + 1” DWT coefficients. After these coefficients are produced, a reconstruction of each frequency scale is made, using up-sampling factors with a 1:2 relation and two reconstruction filters HP’ and LP’ (Figure 3). These reconstruction filters are those which assemble, together with decomposition filters (HP and LP), a quadrature mirror filters system [33].

2.4. NDVI Time-Series Processing

Two types of time-series processing have been made: the

analysis of two related series per study site and the study of image-time series of a region.

In the first case, two time-series of NDVI from all the locations described previously were processed via DWT. The first signal corresponded to a pixel where a LLC occurs and the second one corresponds to a pixel located nearby of the first one which will provide a control signal. Both signals were decomposed in five levels (m = 4) using the Meyer mother wavelet, where each component meaning is described in Table 3.

In the second case, the images acquired were also attacked by the DWT using five levels and the Meyer mother wavelet. In this topic it is necessary to explain that this study did not perform a bi-dimensional WT. This research decomposed the image time series into five image time series where each one will represent variations at a different frequency scale. As a full representation of image time series will require a video to appreciate the decomposition, we used the time series quadratic sum of some frequency scales to represent the changes. The quadratic sum of a signal can be understood as an indicator of changes, where bigger changes are represented by higher values.

3. Results and Discussion

3.1. Single Point Analysis

3.1.1. Seasonal Inundation Analysis in Cambodia & Vietnam

Contrast between dry land and floodable land was already

noticeable in the original time-series f, remaining almost constant in the first one and showing cycles in the other (Figures 4(a) and (b)). The NDVIs₄ decomposition, which represents the inter-annual changes, showed a mean ± std NDVI value of 0.79 ± 0.021 in the dry land and 0.47 ± 0.078 in the floodable land, implicating a lower photosynthetic activity due to the presence of water in the floodable period. The NDVIr₄ reconstruction, which represents slow seasonal changes, had an almost constant value in the dry land and a sinusoidal-like dynamic in the floodable area. This behavior is associated to seasonal variation of NDVI due to inundation, with higher peaks in the first semester and lower values in the second semester, in concordance with the flooding beginning in May-June described by Sakamoto et al. [28]. Finally, the sum of NDVIr₁, NDVIr₂, and NDVIr₃ represented all rapid changes plus rapid seasonal changes in both signals. It was possible to identify that floodable land has more variability than dryland, because of the occurrence of rapid changes originated by the dryland-wetland transition. These variations changed every year and can be related to rainfall variability.

3.1.2. Fire Event Analysis in Spain

Original time series f showed a sudden variation of the NDVI in the year 2006 (Figures 5(a) and (b)) in concordance with a fire event occurred in Galicia in August [29]. Interannual level NDVIs₄ illustrates the effect of fire in the location showing a decrease between 2006 and 2007, and vegetation regrowth starting in 2007 and finishing in 2009. Observing the mid-term changes in NDVIr₃₊₄ it is possible to find a discontinuity in the seasonal cycle during 2006 and 2007, denoted by the decrease of higher peaks in the second semester. It is also noticeable the temporal location of the fire-event given by the rapid changes, described in the high frequency scale reconstruction (NDVIr₂₊₁). Other rapids changes showed in NDVIr₁₊₂ are provoked by spurious image acquisition and intraweek changes.

3.1.3. Farming Analysis in Argentina

Differences between the control signal and the farming signal appear since 2005 (Figures 6(a) and (b)). Clark et al. [30] declares that the area was characterized by a rapid deforestation related to soybean and planted pasture expansion between 2002 and 2006. The interannual component NDVIs₄ of the agriculturized area showed a decrease of NDVI in 2005, from a mean value of 0.68 to 0.53, meanwhile the control signal remained almost constant. The NDVIr₃₊₄ signal also shows an increment in the peak to peak amplitude due to the replacement of the deciduous forest by crops, meanwhile the maximum amplitude before 2005 was of 0.33, the maximum amplitude after the faming increase to 0.64. In this component it is easy to identify the peaks and the date where they took place, telling us that the forest was replaced by a summer crop. Most rapid changes NDVIr₁₊₂ also increased after 2005 due to the irregular curves of crops alternation.

3.1.4. Urbanization Analysis in China

Both original signals do not showed discrepancy until 2004 (Figures 7(a) and (b)) where the original time-series fluctuated from bimodal curves to a sudden decrease in 2004. The NDVIs₄ reconstruction showed the interannual

change from a mean value of 0.48 to 0.22, as a consequence of the replacement of crops by urbanized structures. The sum NDVIr₃₊₄ reconstructions and the sum NDVIr₁₊₂ also expose a change during from 2004. In the seasonal component the NDVI goes from two crop cultivation per year to a single vegetation peak in August, and the reduction of seasonality is denoted by the decrease of the peak-to-peak magnitude. The rapid changes NDVIr₁₊₂ were also reduced as the urbanization reduced the most rapid changes in vegetation.

3.2. Image Time Series Analysis

In the frame of this research we also developed a scrpit to process time series analysis from raster data (Section 2.4). As a result, we generated images that joint spatial dimension with temporal analysis. In Figure 8 we represented the ratio between each pixel most rapid changes quadratic sum and the quadratic sum of the original time series. As we have seen previously in item 3.1.3, agriculture has an impact in rapid changes, therefore it was expected that the cultivated places would have a higher quadratic sum (i.e. the higher values in the image). A visual comparison between two high resolution images of the site (where the anthropic activity was noticeable and the energy image was achieved via WT) showed a good recognition of changes. These type of images can be used

as input layers in a classification, because they resume both spatial and temporal information.

4. Conclusions

The study of land cover use and land cover changes, as part of global change, demands the integration of data sources, models and methodologies, and technologies which provide analyzed data according the necessities (e.g. monitoring, alarms, variability research). In response to these requirements, we proposed an approach based in WT. Through this work we had demonstrated the capacity of this tool to study and characterize four different LLC and submitted the results separated in frequency scales (interannual, seasonal and rapid changes). This disaggregated information allowed us to interpret each involved process without the interference of others. The use of WT in image time series gave us the possibility to characterize different change levels in spatial dimension. Layers generated with WT might be used to pattern recognition in LCC and to improve image classification.

It must be said that in this article we worked using NDVI time series, but this tool can be applied to any time series (e.g. leave area index time series, surface temperature time series and others). The use of multiple nature time series might not only describe the LCC at any time scale but also determine its causes.

5. Acknowledgements

Alfredo N. Campos has an INTA post-grade scholarship. This work was supported by INTA (AERN 294001) and by a grant from the InterAmerican Institute for Global Change Research (IAI) CRN2031, which is supported by the US National Science Foundation (Grant GEO045- 2325).

REFERENCES