Seasonal Pattern of Vegetative Cover from NDVI Time-Series

1. Introduction

Indonesia manages various forested land utilizations, for instance natural forest and plantations. In the past, natural forest had been exploited throughout the country, mainly in the islands of Sumatera and Borneo (Nawir & Rumboko, 2007). Greater criticisms on forest exploitation lead to a moratorium which needs to be monitored frequently. Most of forest concessions were converted into plantations (Kartodihardjo & Supriono, 2000). Popular plantations developed in the country have been rubber and oil palm. Successful management of the plantation as well as natural forest has been under careful examination. Assessment of woody vegetation could be taken using field surveys or remote sensing. Although having possibilities to capture detailed datasets, field survey is lacking in terms of implementation: lengthy data capture and least favorable to remote areas. In these situations, remotely sensed data play a crucial role.

Current practice of the vegetative assessment has been single temporal analysis assisted by remote sensing data; some with extension to multi-temporal analysis. Mono-temporal analysis was shown useful for biomass estimation, using optical (Foody et al., 2001) or Synthetic Aperture Radar (SAR) data (Raimadoya et al., 2005). Nonetheless, changing climate reveals difficulties to determine exact condition of the vegetative land cover. Season is probably the most intriguing factor in this case, especially in long-term land cover changes (Lambin, 1999). The study suggested the importance of seasonality in remote sensing data analysis to account fluctuating forest state under climate change.

In this study, we employed seasonal adjustment of time-series statistical method to understand phenology and detect disturbance on some woody vegetation. Seasonal adjustment methods have been developed particularly to forecast and vastly applied on financial studies. Nonetheless, this has been rarely applied to environmental studies, except in meteorology and hydrology. X12-ARIMA is one of the latest versions of X method that dealing with seasonal adjustment time-series. It decomposes single value of time-series data into three components, namely trend, seasonal and irregular. Seasonal components can be used to test seasonality or impact of season on the vegetation index. Test on global trend of the index and irregular aspect can also be inferred from the method. This research employed NDVI time-series of SPOT VEGETATION. NDVI has been a popular index in vegetation related studies. Foody et al. (2001) found that NDVI has some weaknesses, for instance: (1) relationship between NDVI and biomass was usually asymptotic that limits its ability to represent large amount of biomass accurately, and (2) sensitivity of NDVI to represent vegetation greenness varied upon environmental circumstances. Despite its weaknesses, the advantages of using NDVI to understand vegetative cover phenomena is due to its long-term availability, easy in calculation, and its relationship with some vegetative parameters such as greenness, and biomass (Hess et al., 1996).

2. Remote sensing of tropical forest cover

Remote sensing application on tropical forest was first demonstrated by Tucker et al. (1984). Using Advanced Very High Resolution Radiometer (AVHRR), the study indicated large clearing which was associated with large-scale dwelling. Since then, tropical forest gains particular attentions, employing multispectral, hyperspectral and radar datasets. It appears that Landsat was one of popular multispectral sensor employed for tropical forest assessment and monitoring. Sader et al. (1989) pioneered Landsat applications in tropical forest biomass estimation. Although not highly successful, the research informed that fairly mature successional forest were undetectable by means of Normalized Difference Vegetation Index (NDVI), which has led to many attempts to improve the achievement.

It is well understood that tropical regions have severe cloud cover. Although there have been numerous attempts to reduce or eliminate atmospheric disturbance (for instance a paper by Richter, 1996), the issue persists. Longer wavelength sensor such as Synthetic Aperture Radar has been exploited to minimize atmospheric disturbance (Reddy, 2008). At first, single polarimetric SAR datasets were employed, mainly on VV (vertical transmit and receive) or HH (horizontal transmit and receive). C-band and L-band SAR systems were largely reported, however, higher wavelength (L- or P-band) is preferable (Trisasongko, 2009). Exploiting the advance of processing technology, many SAR sensors have been mounted into spaceborne platform. Currently, users have options to obtain dual polarimetric SAR; some SAR systems allow data capture in fully polarimetric mode. Despite limited sensors capable to acquire in fully polarimetric mode, several applications of fully polarimetric SAR were successfully demonstrated for tropical forest monitoring (for instance Trisasongko et al., 2007; Trisasongko, 2010).

Forests, including the ones in tropical regions, reveal complexity in monitoring. Seasonality has been one of major trends in forest applications of remote sensing, especially in conjunction with climate studies. Studies on seasonal pattern of vegetative cover have been published in attempt to understand impacts of climate change to vegetation growth. Brando et al. (2010), for instance, studied climate variations seasonally and interannually to understand drought in Amazon using time-series data of MODIS Enhanced Vegetation Index. In the same region, Koltunov et al. (2009) presented a study confirming association between logging activity and forest phenology in dry seasons of 2000-2002. Using the time-series of NDVI (Normalized Difference Vegetation Index) of AVHRR or MODIS some scholars showed that vegetation phenology has been affected by climate dynamic particularly precipitation and/or temperature (Wang et al., 2008; Meng et al., 2011).

3. Time-series analysis

Time-series analysis of remotely sensed data, as shown earlier, has gained special attention supported by availability of wide-coverage, high temporal satellite data. Earlier works exploited descriptive analysis, and focusing on visual analysis on graphical representation of remote sensing data of selected sites. Goetz et al. (2006), for instance, presented useful graphical representations of NDVI time-series data to assess post-fire forest regrowth in Canada. Present advance using X12-ARIMA demonstrated that time-series data were indispensable to characterize fire spot in tropical region (Panuju et al., 2010).

Analysis of time-series data has been employed for following aims: (1) to describe characteristics of the data; (2) to explain variance of any variable associated to surrogate variables; (3) to predict a variable beyond the time span; or (4) to control processes (Chatfield, 1984). Underlying assumptions of time-series analysis are variance homogeneity and stationary of the data. Lu et al. (2001) successfully employed decomposition technique to visually separate herbaceous from woody vegetation of temperate region.

Basic formula of X12-ARIMA is ARIMA and the X technique. ARIMA is a standard time series introduced by Box and Jenkins (1970), while the X technique is a seasonal adjustment procedure started from X0 to the recent form the-X12 (Shiskin et al., 1965). Since the background of developing procedure is financial application, most of terminology to process the decomposition represents business and financial related activities. The application of the ARIMA model on environmental study was widely reported by hydrological studies for instance by Hipel and McLeod (1994). The ARIMA model is a modification of ARMA process which is a mixed model of autoregressive and moving average (Wei, 2006):

$${\varphi }_p(B)Z_t={\theta }_q(B)a_t$$ E1

where Ф_p(B)=1-Ф_pB-….- Ф_pB^p, and θ_q(B)=1-θ_qB-….- θ_qB^q. Meanwhile the ARIMA is denoted as follows:

$${\varphi }_p(B){(1-B)}^d{Z}_t={\theta }_0+{\theta }_q(B)a_t$$ E2

where the difference between both equation are on the differencing process represents by 1-B)^dZ_t.

4. Methodology

4.1. Test site

The study was located in Jambi Province, Indonesia (Fig. 1), between 101.15 E – 104.43 E and 0.76 S – 2.77 S. The province is famous for two key national parks: The Berbak National Park, largest swamp forest in South-East Asia which is also a Ramsar Convention on Wetlands site (www.ramsar.org), and The Kerinci-Seblat National Park, a natural habitat of gigantic flowers such as Amorphophallus and Rafflesia. Jambi also hosts Bukit Tigapuluh NP where aboriginal tribes of Talang Mamak and Orang Rimba reside. The province is also an important region for Indonesia which has contributed in oil and gas (BPS, 2010) and oil palm estate (http://aplikasi.deptan.go.id/bdsp/newlok.asp). The important estates include rubber and oil palm, which are managed by locals, government and private companies. Rubber and oil palm estates are found throughout the province.

Figure 1.

Site location

4.2. Remote sensing and field data

In this research, the main remote sensing datasets were time-series NDVI SPOT VEGETATION (SPOT-VGT) and high resolution ALOS AVNIR-2 imageries. SPOT-VGT NDVI data were collected from the Flemish Institute for Technological Research (VITO), Belgium. VITO website (http://free.vgt.vito.be) provided 10-day NDVI composite images, which were then constructed into single database from April 1998 to December 2010. To assist time-series analysis on selected land cover types, a set of ALOS AVNIR-2 imageries provided by Indonesian Ministry of Agriculture were used. The imageries were georeferenced according to Indonesian Base Map specification (WGS 1984, geographic coordinate system). Google Earth (http://earth.google.com) was also utilized to assure locations for samples.

Field surveys were designated to obtain current land use in selected regions. Selection was made considering vastness of the province. First field visit was conducted in October 2010, which covered following regencies (kabupatens): Batanghari, Bungo, Tebo, Muaro Jambi, Sarolangun, Tanjung Jabung Barat and Merangin (Fig. 2). Another field visit was arranged to focus on specific areas of Bungo and Kerinci regencies in July 2011.

Figure 2.

Field survey which successfully collected data in most of regencies of Jambi Province. (a) Oil palm plantation (Tanjung Jabung Barat), and; (b) Rubber plantation (Merangin).

4.3. Data analysis

First step on this research was samples identification for forest, rubber and oil palm area. We selected two regions of interest (ROI) for each land cover to develop time series NDVI of SPOT-VGT. Each land cover was delineated for about 25 pixels of SPOT-VGT or 25 kilometers. All those samples were named Forest1, Forest2, Rubber1, Oilpalm1, Rubber2 and Oilpalm2 (Fig. 3). Forest1 was a wetland forest (Berbak National Park) while forest2 was taken from mountainous region (Kerinci Seblat National Park). Rubber1 was selected from Tebo Regency (101.82 E-102.82 E; 0.88 S-1.91 S) while rubber2 was from Merangin Regency (101.55 E-102.57 E; 1.63 S-2.77 S). Then, the oilpalm1 and oilpalm2 were from Tanjung Jabung Barat (102.64 E-103.66 E; 0.75 S-1.50 S) and Merangin Regency respectively.

NDVI from each 25 samples were then averaged. Since the NDVI of SPOT-VGT was the 10-days synthesis product (S10) and procedure of X12-ARIMA only accommodated monthly data, we computed maximum value of the S10 NDVI. Some researchers including Jia et al. (2002) noted that maximum of NDVI will produce best representation of monthly data and anticipate effect of cloud cover which reducing NDVI value.

Time series analysis was then processed using X12-ARIMA procedure. It was utilized to determine seasonal model of the NDVI series of each sample. The model usually was represented as (p, d, q) (P, D, Q), where p symbolizes autoregressive order, d denotes differencing and q stands for moving average order. According to Wei (2006), seasonal ARIMA model is denoted as follow:

$${\Phi }_P(B^s){\varphi }_p(B){(1-B)}^d{(1-B^s)}^D{Z}_t={\theta }_q(B){\Theta }_Q(B^s)a_t$$ E3

Where Ф is autoregressive function and θ expresses the moving average process. Further discussion on seasonal model refers to Findley et al. (1998). The model was analyzed using ARIMA auto-modeling software developed by Eurostat which freely accessible in 2002. Based on the auto-modeling, the model will be selected from seasonal or non seasonal model with some alternatives order of autoregressive, differencing and moving average order. Since the X12-ARIMA decomposing series data into three components, we used the decomposition process to compare differences among land covers. Even though there are multiplicative and log alternatives, all series were only performed the additive model option. The model was tested based on residual analysis using BIC (Bayesian Information Criterion) or SBC (Schwartz’s Bayesian Criterion). The parameter was calculated using this equation:

$$SBC(M)=N\ln {\widehat{\sigma }}_a^2+M\ln N$$ E4

where M denotes number of parameters and N refers to number of series data. The best model will produce the smallest BIC or SBC. Furthermore, to test randomness of error L-jung Box suggested probability calculation of χ². If χ² is less than 10%, then ARIMA automodel will reject the option. We employed TRAMO (Time series Regression with ARIMA noise, Missing value and Outliers) option to understand history of series. The TRAMO generated possible option to explore any noise such as missing and extreme value due to any disturbances into the series. Indicative option will lead conclusion into the certain time of evidences perturbing NDVI.

Figure 3.

Sampling locations.

5. Results and discussion

5.1. Comparison of forest phenologies

Study on phenology of vegetative cover is important to understand plant and animal behavior. Seasonality of vegetation behavior is apparently a resultant of micro- and macro-climatic aspects as well as other living species activities. Studies on phenology by biologists mostly focus on physiological aspects of plants and phenophases, including leafing, flowering and fruiting (van Schaik et al., 1993). In this study, we discuss a different perspective to investigate phenology. Physiological aspect combined with biotic and abiotic factors as well as sensors condition and disturbance were assumed simultaneously epitomized by vegetation index (NDVI). Decomposition of NDVI series was expected to explain plant characteristics and activities occurring in the tropical forest. Since the spatial resolution of SPOT Vegetation is considered coarse, decomposition process was intended to deliver some details of the series. Lu et al. (2001) showed that decomposition of NDVI time series of AVHRR, particularly trend and seasonal components, was able to separate woody and herbaceous vegetation by using STL (seasonal trend decomposition based on loess). This research extended Lu et al. (2001) results to decompose NDVI series of forested and plantation areas.

Natural forest cover in Jambi is usually managed as a National Park. Under this scheme, forest cover is expected to be undisturbed and preserved to play its ecological function. Nonetheless, threads remain, mostly for plantations or other anthropogenic activities such as shifting cultivation. Some reports note that local inhabitants surrounding national parks usually take advantage of forest for their daily life. This phenomenon should be understood as a function of forest to support people nearby. In Jambi, involvement of people on forest utilization apparently depends on ecosystem itself and access into the forest. Therefore, it was necessary to obtain at least two different types of forest, wetland and montane forests. In these areas, we attempted to explore decomposition of time series NDVI to understand evidences on both ecosystems.

The general parameters of ARIMA model is displayed in Table 1. Apparently, both forests produced similar model with some parameters were a slight better on Kerinci-Seblat National Park (KSNP). The (p, d, q) indicates that both wetland and montane forests fit into ARIMA model. Seasonal model is shown fairly obvious on both forest types, indicated by (P, D, Q) units. This suggests that seasonality was quite apparent on both forests, hence monitoring solely based on single date should be avoided. The results also advise that different seasons should be taken into account for selecting imageries in multitemporal analysis, for instance by means of higher optical or SAR sensors. Due to lacking of available datasets, imageries of different season in multitemporal analysis such as land use detection and change have been utilized. In this case, as this research suggests, detailed assessment are needed to avoid bias in detection. The standard error (SE) of residual and BIC of KSNP series were smaller than from Berbak National Park (BNP). The model is modestly better fitted on KSNP.

Furthermore, the decomposition of NDVI in two Jambi forests is presented in Figure 4. It is shown that original series as well as their decomposition of both forest ecosystems were quite different. The original series of BNP wetland forest which is located nearby the Eastern coastline was more stable rather than samples obtained from the mountainous forest of KSNP. According to van Schaik et al. (1993), the phenology of tropical plants has been affected mostly by biotic and abiotic factors, such as water, light, and nutrients availability.

Parameters	BNP forest	KSNP forest
(p, d, q) (P, D, Q)12	(0, 0, 0) (0, 1, 1)	(0, 0, 0) (0, 1, 1)
SE(res)	0.083	0.077
Q-val	12.431	12.489
BIC	-254.704	-271.326

Table 1.

Parameters of seasonal time series model for forests.

Figure 4.

Original and decomposition series of forest taken from Berbak National Park (a) and Kerinci Seblat National Park (b).

Continuous water availability in wetland area was suspected to be a primary cause on the stability. Original NDVI amplitude of BNP is about 0.5 with minimum value of 0.4 and the largest at about 0.9. Meanwhile, the NDVI amplitude of KSNP is around 0.7 with minimum value 0.2 and the maximum at 0.9. It suggests that disturbances were severe at the montane forest. Second field visit supports this argument. During the visit, several shifting cultivations in the KSNP were observed. Local inhabitants have planted several economic commodities such as cinnamons, duku (Lansium domesticum Corr.) and pinang (Areca catechu L.). It appears that there has been a mutual agreement between the government and local inhabitants to access some of KSNP areas as agroforestry to support local economies. In turn, local inhabitants are asked to assist KSNP protection program.

Trend series of both types of forests were similar. However, it is discernible that seasonal factor of BNP illustrates larger amplitude than KSNP. This might be related to the fact that BNP is located in wetland areas (altitude of 0-20 meters), where tidal dynamics (about 2-2.5 meters, according to Wetlands International, Indonesia Office) play a role at different season. Another interesting finding was that the irregular components of KSNP fluctuated more than BNP. This would be the subject of future investigation.

5.2. Comparison of plantation phenologies

Jambi has been one of important provinces who manage an immense coverage of plantations. Oil palm and rubber have been the popular estate commodities in Jambi, along with cinnamon, tea and Aquilaria (local name: gaharu). However, only oil palm and rubber are discussed in this chapter. Results of seasonal time series modeling are presented in following table.

Parameters	Oilpalm1	Oilpalm2	Rubber1	Rubber2
(p, d, q) (P, D, Q)12	(0, 0, 0) (0, 1, 1)	(0, 1, 1) (0, 1, 1)	(0, 0, 0) (0, 1, 1)	(0, 0, 0) (0, 1, 1)
SE(res)	0.073	0.084	0.101	0.085
Q-val	18.722	19.684	21.851	23.272
BIC	-296.633	-240.155	-202.475	-253.079

Table 2.

Parameters of seasonal time series model for plantations.

It is clear that seasonality was important to both oil palm and rubber at all cases. This also suggests that remote sensing data analysis involving more than two acquisitions should be carefully taken to avoid preconceived notions. At some points, oil palm produced a better certainty than rubber. Some parameters such as standard error of residual and BIC of oil palm were smaller than of the rubber (Table 2).

To save space, we presented only Merangin region for the assessment of oil palm and rubber plantations (Fig. 5). The results show that the original series of NDVI of rubber generated smaller amplitude than the oil palm. NDVI of rubber spans from 0.6 to 0.9 while oil palm slightly varies, from 0.5 to 0.9. This is due to gaps often found in the oil palm blocks due to various reasons, including diseases. The moderate to large oil palm plantations (about tens of thousand hectares) are usually set up factories on-site, which would severely reduce NDVI values. The trend components of both plants were fairly similar. It tends to increase by time depending on age. Seasonal trend shows different pattern during the periods. Irregular components indicate noise and any disturbances on specific times, which cannot be describe further at this time.

Figure 5.

Original and decomposition series of oil palm (a) and rubber (b) plantations.

5.3. Disturbances on forest and plantations

On previous subsections, disturbances could be indicated by time series datasets. To observe more on this disturbance, the TRAMO was employed in identification of any unusual observation of the series. The unusual observation could be due to biotic or abiotic causes, which was not determined at the moment. Instead, our focus was in detection ability of TRAMO. Result of TRAMO detection is shown in Table 3.

Outlier type	Forest1	Forest2	Oilpalm1	Oilpalm2	Rubber1	Rubber2
AO	(Oct-2006)		-	-	-	-
TC	-	(Jan-2002); (Dec-2003)	-	-	-	(Dec-2003)
LS	-	-	-	-	-	-
IO	-	-	-	-	-	-

Table 3.

Detection of disturbance by TRAMO.

Comparison with higher remote sensing data on Forest1 (BNP) indicated false-alarm produced by TRAMO. Using Landsat ETM+ images dated June 3, 2007 (LE71250612007154SGS01) and July 5, 2007 (LE71250612007186PFS00), the forested area has been shown intact with minimum disturbance by human. Detection on KSNP (Forest2) could be considered successful as indicated by second field trip. The survey observed that agroforestry estates (cinnamons, duku and pinang) were at mature stage (various ages). Merangin site which hosts Rubber2 is considered dynamic regency where several estate conversions (changing commodities) were reported.

6. Conclusion

NDVI has been proven useful to study the dynamic of vegetative covers using time-series approach. In this study, seasonality in NDVI time-series data was explored to understand various responses provided by woody vegetation. Using SPOT-VGT data, we showed that seasonality should be taken into account in multitemporal or time-series analyses. We found that NDVI time series of wetland and montane forests could be modeled by seasonal ARIMA, although level of confidence was found higher on montane forest (KSNP). Seasonal models of both types of forest were similar as well. Exploration on seasonal ARIMA model also found that NDVI series of wetland forest (BNP) tend to be more stable, which indicated a stable environment with least deviations on water, light and/or nutrients. Higher NDVI span on KSNP signified variable response from vegetative cover. Field surveys indicated some contributions from agroforestry (cinnamons, duku and pinang) in KSNP which might be contributing factors to more fluctuating NDVI values. NDVI ranges were found quite similar between oil palm and rubber plantations. Slightly bigger range on oil palm might be related to some gaps in plantation; one of them is Crude Palm Oil (CPO) on-site production facilities. This research also discovered some disturbances from time-series data. Although a false-alarm fault was observed, TRAMO was shown useful to identify disruptions in KSNP which was related to agroforestry.

Acknowledgments

We are indebted to VITO who provided SPOT-VGT datasets, Ms. Gusmaini for her assistance in downloading the datasets and Ms. Lili Suryani in field surveys. Several Landsat ETM+ images were downloaded from Earth Explorer provided by the USGS.