Fisher–Shannon and detrended fluctuation analysis of MODIS normalized difference vegetation index (NDVI) time series of fire-affected and fire-unaffected pixels

ABSTRACT MODIS-NDVI data from 2002 to 2014 were analysed to evaluate the effect of fire on vegetation in a test site located in Daxing'anling region (Inner Mongolia and Heilongjiang Province). Fire-affected and fire-unaffected areas were processed using two statistical approaches: detrended fluctuation analysis (DFA) and Fisher–Shannon (FS) method. The DFA allows the detection of scaling behaviour in nonstationary signals, whereas the FS method permits to identify the organization/order structure in complex signals. Our findings show that the results obtained by jointly using the two methods are consistent, enabling the characterization and discrimination between the fire-affected and fire-unaffected areas. In particular, among the investigated indices, the mean value of Fisher information measure (FIM) represents the most significant in discriminating between burned and unburned sites; its mean value for burned sites is about 2.5 that is significantly larger than that obtained for unburned sites (∼1.3). FIM is also characterized by the larger effectiveness in discriminating the two classes of sites on the base of its receiver operating characteristic based performance. The scaling exponents estimated by means of the DFA of the fire-affected pixels are averagely higher than those of the fire-unaffected ones, which, furthermore, are characterized by lower organization and higher disorder degree. Both of the two methods would contribute to identify the impact of fires on vegetation.


Introduction
A sustainable and continuous monitoring in terms of type and status of vegetation cover and ongoing changes is a key issue for the management of ecosystem functionality and for implementing effective strategies to mitigate disturbance phenomena as drought, wind, storm, fire, etc. on environment and ecosystems (Yapp et al. 2010;Lawley et al. 2016). Fire is one of the most critical factors of disturbance in worldwide ecosystems (Carvalheiro et al. 2010;Cruz & Alexander 2013;Eastaugh & Hasenauer 2014). The effects of fires on soil and vegetation, depend on many factors, among them are fire frequency, fire severity and plant resistance (French et al. 2008;Keeley 2009). The characterization of vegetation post-fire behaviour is a fundamental issue to model and evaluate the fire resilience, which is the ability of vegetation to recover after fire (Lentile et al. 2006).
During the last decades, the high number and repeated occurrence of severe wildfires recorded in different regions of the world (UNCCD 1994;FAO 2010) have strongly heightened the negative effects of fire on vegetation. The recent increasing occurrence of fire events is due to many reasons, among them are the abandonment of rural areas, variations in land use practice and climate change. The increasing fire frequency can induce significant variations in vegetation fire resilience and this highlights the need of evaluating the vegetation reaction to fires (D ıaz-Delgado et al. 2003). This is of primary importance since changes in the status or types of vegetation play an active role in ecological processes (such as productivity level, creation of altered patches, modification in vegetation structure and shifts in vegetation cover composition), as well as in land surface processes (such as surface energy, water balance and carbon cycle) (Huemmrich et al. 1999).
The assessment of fire impacts on ecological resources requires investigations performed at different temporal and spatial scales, from local up to regional level. In such a context, satellite technologies can be profitably used for investigating the dynamics of vegetation re-growth after fire disturbance at different temporal and spatial scales (Lasaponara 2006a(Lasaponara , 2006b(Lasaponara , 2013. Fire-induced dynamical processes are very difficult to study since they affect the complex soil-surface-atmosphere system, due to the existence of feedback mechanisms involving human activity, ecological patterns and different subsystems of climate (Telesca & Lasaponara 2010). The use of satellite time series along with statistical analysis techniques can be helpful in understanding the functional characteristics of vegetation dynamics and enable the reporting of ongoing trends at a detailed level (Telesca et al. 2007).
The remote sensing of vegetation has been traditionally carried out by using vegetation indices, which are quantitative measures, based on vegetation spectral properties that attempt to measure biomass or vegetative vigour.
The vegetation indices operate by contrasting intense chlorophyll pigment absorption in the red band of the electromagnetic spectrum against the high reflectance of leaf mesophyll in the near infrared band. The simplest form of vegetation index is simply a ratio between two digital values from these two spectral bands. The most widely used index is the well-known normalized difference vegetation index (NDVI) (Rouse et al. 1973(Rouse et al. , 1974Jiang & Huete 2010). NDVI is used to identify vegetation health by measuring the difference of the infrared reflectance with that of the visible red band. NDVI provides information about the spatial and temporal distribution of vegetation communities and vegetation biomass. High values of the vegetation index identify pixels covered by substantial proportions of healthy vegetation (Tucker 1980). NDVI is indicative of plant photosynthetic activity and has been found related to the green leaf area index and the fraction of photosynthetically active radiation absorbed by vegetation. Therefore, variations in NDVI values are indicative of variations in vegetation composition and dynamics. This paper deals with the characterization of vegetation post-fire dynamics carried out in the Daxing'anling region (Heilongjiang Province and Inner Mongolia) in the north-eastern China. This region is selected due to the high frequency of large fires. The investigations are performed by using a remote sensing satellite data-set, which provides a valuable basic data source to monitor change and understand ongoing degradation process. For the purposes of our analyses, the temporal series from 2002 to 2014 of NDVI satellite MODIS data were processed using two statistical approaches: the detrended fluctuation analysis (DFA) and the Fisher-Shannon (FS) method. Both the methods are efficiently used to investigate the time dynamics of complex and nonstationary signals as those analysed in the present paper.
The objectives of this study are: (1) verify the reliability of DFA and FS methods in the detection of the fire-affected pixels and (2) provide a quantitative characterization of the vegetation recovery capability.

Study area and data-set
This study was carried out in the Daxing'anling region (Figure 1), which is located in Heilongjiang Province and Inner , in the north-eastern China. This region is selected because the forest fires have occurred very often and also have severely affected large areas. The total area of Daxing'anling region is about 250,000 km 2 . The forest coverage of Daxing'anling region is about 62% (Liu 2014). There is about 81,300 km 2 of the forest that cover the Inner Mongolia (Wang 2014). The main land covers of this area ( Figure 2) are: coniferous forest (Larix sp., Pinus Sylvestris, Picea sp.), broadleaved deciduous forest (Quercus sp., Betula sp.) and mixed forest (Liu 2014;Wang 2014). This study area experiences a cool continental monsoon climate, which is very cold in winter but warm in summer. The weather is relatively dry and the precipitation amount is low in spring and autumn (65 mmand 80 mm) (Liu 2014;Wang 2014). Thus, fires happen more and more frequently in these seasons as well as in summer due to the high temperature. Several big fires occurred in Daxing'anling region. The fire events are listed in Table 1. In this study, the fires occurred on 26 June 2010 were taken into consideration for analysing the difference between burned area and unburned area. The fires are located approximately at 123.48 E and 51.3 N, which is near the border of Inner Mongolia and Heilongjiang Province.
Data from the MODIS sensor on-board the Aqua satellite were used in this study. In order to reduce the contamination effects that were generated by residual clouds, atmospheric perturbation, variable illumination and viewing geometry in daily NDVI maps, the 16-day composite NDVI product MYD13Q1 that has a 250-m spatial resolution was taken into consideration in this study. Such data can be downloaded from the Land Processes Distributed Active Archive Center (LP DDAC: https://lpdaac.usgs.gov/).
The MODIS data were geometrically rectified and re-projected to UTM WGS84. Then, the image was spatially resized on the study area. According to the analysis of Miller and Thode (2007) and Lanorte et al. (2013), the criterions (formulas (1)-(3)) were used to discriminate the burned area from the unburned area. And then, the burned area and unburned area can be extracted after this procedure by using the region of interest tool of ENVI software. NDVI pre > NDVI post ; (1) NBR ¼ r 2 À r 7 ð Þ = r 2 þ r 7 ð Þ ; (2) dNBR ¼ 1000 Â ðNBR prefire À NBR postfire Þ 400 ; (3) where NDVI pre is the NDVI of pre-fire, while NDVI post is the NDVI of post-fire. NBR is the normalized burn ratio, considered the best index to identify a burned area which is calculated normalizing the difference between the reflectance of the two burn sensitive bands, infrared (NIR) and shortwave infrared (SWIR). At this aim, we used the band 2 (NIR) and Band 7 (SWIR) of MODIS-Aqua daily surface reflectance images at 500-m spatial resolution (MYD02HKM -Level 1B Calibrated Geolocated Radiances). The dNBR is the difference of NBR between pre-fire and post-fire.

Methods
NDVI time series of burned and unburned areas are both characterized by very complex signal variability, and a simple visual inspection or simple standard statistics could not allow to identify the dynamics governing the time series of pixels covering the fire-affected and fire-unaffected areas. The fire lasted until 3 July 2010 and consumed a large amount of fire protection resources. Therefore, robust nonlinear time series tools are necessary to disclose inner features that characterize both types of pixels. In our study, the statistical analysis of the pixel time series is performed by using the FS method and the DFA. The first method allows to identify organizational/order properties in time series (Martin et al. 2001), while the second one permits the robust detection of scaling behaviour in time series in order to reveal persistent/anti-persistent features (Kantelhardt et al. 2002).

The Fisher-Shannon method
The Fisher information measure (FIM) and the Shannon entropy, well-known statistical quantities in the theory of information, are efficiently employed to characterize the complex dynamics in nonstationary time series. The FIM provides information on the organizational structure of a signal quantifying its degree of order (Martin et al. 2001), while the Shannon entropy is used to quantify the degree of uncertainty or disorder of a time series (Shannon 1948). Various geophysical and environmental complex and nonstationary signals, like volcano-related continuous seismic signals , earthquake-related electromagnetic signals ), magnetotellurics , atmospheric particulate matter time series (Telesca et al. 2008) and wind series , were analysed by using the FS method to get an insight into the dynamics underlying their inner time variability: Let f(x) be the probability density of a signal x, then its FIM I is given by and its Shannon entropy H X is defined as (Shannon 1948) Instead of H X , it is generally used the Shannon power entropy as N X, Since both the FIM and the Shannon entropy power are calculated on the basis of the probability density function f(x), it is important that this is estimated with good accuracy, especially when observational time series are investigated. The estimation of the f(x) can be carried out by means of the kernel density estimator technique (Devroye 1987;Janicki & Weron 1994) as shown in the formula below:f wheref M ðxÞ represents the estimation of f(x), with b the bandwidth, M is the number of data and K(u) is the kernel function, which is a continuous non-negative and symmetric function satisfying the two following conditions: Following Troudi et al.'s (2008) and Raykar and Duraiswami's (2006) algorithms, the Gaussian kernel with zero mean and unit variance is used.
3.2. The detrended fluctuation analysis The DFA (Peng et al. 1995) is a well-known method to detect scaling behaviour in nonstationary signals and its effectiveness and easiness of calculation makes DFA one of the most used methods to reveal long-range correlation structures in observational time series in diverse research fields (Zheng et al. 2008;Telesca et al. 2012Telesca et al. , 2016.
The DFA comprises the following steps: (1) the time series x(i), where i = 1, . . . ,N, and N is the length of the series integrated where hxi is the mean value of the series; (2) the integrated series y(k) is divided into boxes of equal length n; (3) for each n-size box, y(k) is fitted by least square line y(n,k) that is the estimated trend in that box, and detrended by subtracting the local trend y(n,k) in each box; (4) for a given n-size box, the root-mean-square fluctuation, F(n), is calculated: (5) the above procedure is repeated for all the available box sizes n to provide a relationship between F(n) and the box size n, which for long-range power-law correlated signals is a power-law (6) the scaling exponent a quantifies the strength of the long-range power-law correlations of the series x(i): if a = 0.5, the series is uncorrelated; if a > 0.5, the series is persistent, where persistence means that a large (small) value (compared to the average) is more likely to be followed by a large (small) value; if a < 0.5, the correlations of the series are anti-persistent, which indicates that a large (small) value (compared to the average) is more likely to be followed by a small (large) value (Mandelbrot 1982;Schertzer & Lovejoy 1987;Marshak et al. 1994;Calif & Schmitt 2014).

Analysis of burned area
The NBR of pre-and post-fire are shown in Figure 3(a,b). It points out that the NBR of some areas in post-fire, which is presented in black colour in Figure 3(b), is smaller than that in pre-fire. This is due to the influence of the forest fire occurred on 28 June 2010. The dNBR between pre-and post-fire is shown in Figure 3(c). It presents clearly that some areas have large variation in the measurement of dNBR, which is depicted in bright white. Thus, it is feasible to identify the fire-affected pixels based on the measurement of dNBR. We extracted 3200 pixels of the fire-unaffected pixels based on the visual inspection for analysing the variation between pre-and post-fire. Figure 3(d) shows the scatter plot of unaffected pixels in NBR of pre-fire versus dNBR. It can be seen that the dNBR Ã 1000 values of most of the fire-unaffected pixels are around from ¡50 to 100 and all of these values are under 400. Thus, the threshold of dNBR Ã 1000 is set reasonably.

The time series NDVI trends
The investigations were carried out in both fire-affected areas (the areas that are depicted in red colour in Figure 1) and fire-unaffected area (the area that is depicted in green colour in Figure 1), which are located in Daxing'anling region. Using the time series NDVI-MODIS data (MYD13Q1), 360 pixels were extracted from the fire-affected and fire-unaffected areas respectively, in the 2002-2014 period. The trends of the time series NDVI data of three pixels belonging to the 360 fire-affected pixels and fire-unaffected pixels are presented in Figure 4(a) and Figure 5(a), respectively. And the NDVI mean values of the fire-affected pixels and unaffected pixels were calculated and shown in Figure 6(a,c). It can be seen that the NDVI of the fire-affected pixels has a sharp decrease when fire occurs. The details of these six pixels are shown in Table 2. All these pixels correspond to the mixed forest vegetation type which consists of mixtures of broadleaf and needle-leaf trees. The NBR of the fire-affected pixels shows greater variations after fire than that of the fire-unaffected pixels according to the column of dNBR. The pixel time variations undergo seasonal periodicities, which have to be removed before applying the FS method and the DFA, in order to avoid any cycle-induced effect on the results. The seasonality was removed by computing for each 16-day composite the normalized departure defined as follows, NDVI d = (NDVI ¡m NDVI )/s NDVI . The composite mean m NDVI is calculated for each calendar day, e.g. 1 January, by averaging over all years in the record. Analogously the standard deviation s NDVI is computed. This filtering procedure for satellite data was already applied by Telesca and Lasaponara (2005, 2006a, 2006b. Figures 4 and 5 show the time variation of the normalized departure. We can see that the seasonal fluctuations (that are indicated by the quasiperiodic behaviour shown by the pixel time series) are not evident anymore, and this indicates that the filtering procedure was correctly done.
The NDVI d of three pixels of the fire-affected and fire-unaffected are shown in Figures 4(b) and 5(b). And the mean values of NDVI d of the fire-affected pixels and unaffected pixels are shown in Figure 6(b,d).
It can be concluded from Figures 3-5 that the NDVI changes periodically every year. It reaches the peak in June or July and goes down to the minimum in January or December. However, the NDVI declines sharply when fires occur, and the vegetation needs some time such as some months or even some years to recover to the level of pre-fire. In this study, we consider both fire-affected and fire unaffected pixels in order to characterize the post-fire dynamics in vegetation through variations in NDVI values. The comparison of the temporal behaviour of the fire-affected and fire unaffected pixels enables us to assess that the variations observed in NDVI values are actually indicative of variations in vegetation conditions and post-fire dynamics rather than other factors such as weather parameters.

Statistics analysis of pixels time series data
To each NDVI d time series, the FS method and the DFA were applied. Figure 7(a) shows the FS information plane, in which each pixel is identified by a point in the plane, whose x-axis is the Shannon entropy power and the y-axis is the FIM. Figure 7(b,c) show the histograms of the FIM and N X values for the burned and unburned pixels. The FIM and Shannon entropy values can be used to classify the fire unaffected (black circles) by the fire-affected (red circles) pixel. The FS analysis suggests that the order/organization structure of the time series of the pixels affected by fire is higher than that of the pixel not affected by fire. Figure 8 shows the histograms of the scaling exponents obtained by applying the DFA on the time series of the same pixel as in Figure 7. The mean ( §standard deviation) of the FIM, N X and a for the fire-affected and unaffected pixels along with the t-value of Welch's test are summarized in Table 3. In particular, the t-value indicates that mean of the parameter's values for the burned pixels is significantly different from the mean of the same parameter's values for the unburned pixels. The values of the DFA scaling exponents and the FS values for the fire-affected and fire-unaffected pixels suggest that it could be possible to reliably discriminate  between the fire-affected and fire-unaffected areas. However, there is a certain overlap between the clouds of points representing both types of pixels that is very probably due to the physical conditions of the investigated pixels. In other words, the pixels corresponding to the overlap could represent sites partially burned, and, thus, their classification is 'on the border' between the fire-affected and fire-unaffected. The parameters derived from the FS method (FIM and N X ) and the DFA (scaling exponent) are useful to describe the dynamics of the pixel time series of the fire-affected and unaffected areas. The consistence between the results obtained by applying the FS and those obtained by applying the DFA to the same satellite pixel time series enables to reliably characterize and discriminate the fire-affected and unaffected areas, on the basis of the DFA scaling exponents and the values of the FIM and Shannon entropy.
The average higher values of the DFA scaling exponent for the fire-affected pixel time series indicate that the temporal fluctuations of the fire-affected vegetation are characterized by a more intense persistent character, which governs the vegetation dynamics by means of a sort of positive feedback mechanisms that aim at destabilizing its status and forcing its recovery after the occurrence of the fire. The ability in recovery of vegetation after fire (resilience) was also evidenced in Telesca and Lasaponara (2005), who investigated the correlation properties of the pixel time series of the 1998-2003 decadal NDVI satellite SPOT VEGETATION data for the fire-affected and fire-unaffected areas. In their study, a significantly higher DFA scaling exponent was estimated for the fire-affected sites, suggesting the role played by fires in forcing a more unstable dynamics for vegetation re-growth processes. Similar vegetational patterns are identified in the time series of post-fire NDVI SPOT VEGETATION pixels (Telesca & Lasaponara 2006b), which are also characterized by lower FIM and higher Shannon entropy power (Lanorte et al. 2014).
A quantitative check of the discrimination performance of the FS and DFA parameters can be carried out by means of the receiver operating characteristic (ROC) analysis that has long been utilized to evaluate the performance of binary classifiers (Kharin & Zwiers 2003). The ROC curve represents the functional relationship between the true positive rate, also called sensitivity (ratio between the number of positives that are correctly classified and the total number of positives) against the false positive rate (ratio between the number of negatives that are incorrectly classified  Figure 8. Histograms of the scaling exponents calculated by using the DFA for the fire-affected and fire-unaffected pixels covering the investigated area. The scaling exponents of the time series of the fire-affected pixels are averagely larger than those of the fire-unaffected ones. and the total number of negatives). Indicating as specificity (SP) the probability of classifying a negative occurrence as negative, the false positive rate is equal to 1-SP. A classifier is considered good or not on the basis of its position in the ROC space. The point (0,1) represents the perfect classifier. One classifier is better than another, if its representative point in ROC space is to the north-west of that of the other. If a classifier is represented by a point on the diagonal line y = x, then the classification is random. A classifier is better than another if its ROC curve is located above that of the other, or if the area under its ROC curve (AUC) is larger. Therefore, the performance of the classifier can be evaluated also from the numerical value of the AUC. From Figure 9, it is clearly visible that FIM is the best classifier of the fire-affected and unaffected pixels, since its ROC curve is above those of the N x and DFA. Its AUC is the largest (0.95) while that of the DFA is the smallest (0.768).

Discussion
Extracting trends from vegetation time series can be challenging due to noise, seasonal variations and other factors that may cover the inter-annual subtle signals. Therefore, in order to assess the reliability of the identified subtle trends, the possibility of using different approaches useful to describe small ongoing vegetational variations is desirable. Herein, in order to assess the complex dynamics of fire effects on vegetation cover, the DFA along with the FS information plane (defined by the FIM and the Shannon entropy power (N X )) were used.
These indices measure different levels of temporal correlation of the investigated time series, thus providing an estimation of fire recovery capabilities in different ways.
In particular, DFA allows the detection of scaling behaviour in nonstationary signals. Shannon's index measures information linked, in this case, to the diverse behaviour after fire events in the system vegetation -fire. Fisher's a is a parameter of the log series which is often used as a diversity index and therefore in the current investigation it pertains to the different post-fire recovery capabilities. In the FS plane, the pixel series seem to aggregate into two different clusters corresponding to two different fire recovery behaviours/capabilities after fire.
In the investigated time series results obtained from both DFA and FS well fit each other, thus confirming the robustness and reliability of the adopted statistical analyses. The main finding of this study is that the joint use of the two robust statistical approaches enables a reliable identification and description of pattern of observed NDVI time series data, allowing us (1) to identify the impact of the fire phenomenon on a given pixel: burned, unburned, partially burned; (2) to assess forecasting schemes (for instance, by extrapolating the identified pattern for predicting future values).
The different effects of fire, even on the same type of vegetation cover, is due to the: (1) diverse fire severity degrees that may affect the samples (pixels) under investigation and (2) the fact that each pixel provides a signal that is an average of its whole spatial cover (for example, only a small part may be actually involved in the fire event). Nevertheless, although the coarse resolution of satellite images used in this study may limit the range of spatial details that can be detected, they have clearly proved to be reliable also for the identification and monitoring of sub-pixels changes. In these cases, if required, further detailed analyses could be conducted using higher resolution data to cope with the specific needs of heterogeneous covers or sub-pixel details. DFA and FS plane jointly used with satellite time series can be reliably used as tools for mapping, quantifying and monitoring fireinduced changes in vegetation behavioural patterns and physical vegetation variations at different space and time scales.

Conclusions
In our analysis, it is evidenced that satellite time series investigated by means of robust statistical tools, as are DFA and FS, could be supportive for a sustainable monitoring and an operative management of vegetation at different temporal scales. In fact, it was evident that the use of DFA and FS enabled us: (1) to identify the physical variation in vegetation that, in the current case under investigation, is the occurrence of fire event and therefore the discrimination of the fire-affected from fireunaffected areas and (2) the different vegetational behaviour and recovery capabilities exhibited by the vegetation over the investigated time window.
As a whole, we would like to highlight that the use of space technologies can optimize the time consuming and expensive field analysis, enabling the monitoring of large areas, helping the decision-makers in the management of burned area and in designing strategies for vegetation recovery.
Furthermore, the proposed integrated approach can be fruitfully applied to similar investigations in several different application fields, such as, for instance, those addressed to extract subtle land degradation-related signals that even if small are, however, significant for mapping changes and provide an early warning of ongoing degradation processes.