Analysis of Multifractal and Organization/Order Structure in Suomi-NPP VIIRS Normalized Difference Vegetation Index Series of Wildfire Affected and Unaffected Sites by Using the Multifractal Detrended Fluctuation Analysis and the Fisher–Shannon Analysis

The analysis of vegetation dynamics affected by wildfires contributes to the understanding of ecological changes under disturbances. The use of the Normalized Difference Vegetation Index (NDVI) of satellite time series can effectively contribute to this investigation. In this paper, we employed the methods of multifractal detrended fluctuation analysis (MFDFA) and Fisher–Shannon (FS) analysis to investigate the NDVI series acquired from the Visible Infrared Imaging Radiometer Suite (VIIRS) of the Suomi National Polar-Orbiting Partnership (Suomi-NPP). Four study sites that were covered by two different types of vegetation were analyzed, among them two sites were affected by a wildfire (the Camp Fire, 2018). Our findings reveal that the wildfire increases the heterogeneity of the NDVI time series along with their organization structure. Furthermore, the fire-affected and fire-unaffected pixels are quite well separated through the range of the generalized Hurst exponents and the FS information plane. The analysis could provide deeper insights on the temporal dynamics of vegetation that are induced by wildfire.


Introduction
The natural phenomenon of vegetation dynamics can reflect the status, types, and changes of terrestrial biomass. The analysis of long-term monitoring data of vegetation helps to better understand the variations in vegetation communities and ecosystems under the abnormal disturbances, especially the wildfires [1]. As one of the most destructive natural hazards, the frequency and severity of wildfires have significantly increased in recent decades with the global warming and climate changing [2,3], which can greatly perturb the vegetation growth and dynamics [4,5]. Therefore, the investigation of the characteristics of vegetation reactions that are provoked by wildfire is imperative for the assessment of vegetation dynamics and resilience [6,7].

Data
The vegetation index NDVI [14] is formulated, as follows: where ρ NIR refers to the atmospherically corrected reflectance of the near-infrared band, while ρ red represents that of the visible red band. This index has been the most widely used, because it efficiently measures the status of growth, distribution, and change of the surface biomass [31]. The higher the NDVI value, the better the conditions of vegetation. Since the vegetation can be highly reflectant in the near-infrared spectral domain while being strongly absorbant in the visible red domain, the NDVI value significantly decreases when the vegetation biomass is affected by wildfires as the reflectance of near-infrared band decreases and that of visible red band increases [32].
The NDVI time series analyzed in this study are obtained from the Suomi-NPP VIIRS vegetation indices Version 1 data product VNP13A1 [15,33]. The product has been designed to provide global NDVI values since the year 2012 by selecting the optimal available pixel in a frequency of 16 days at the spatial resolution of 500 m. The NDVI were firstly pre-processed to filter out the anomalous pixels that may have been affected by atmosphere interference, cloud contamination, or bidirectional reflectance distribution function (BRDF) effects or a combination of these factors [34]. This pre-processing procedure is based on the combination of the quality information implemented in the VNP13A1 product, which is the pixel reliability metric informing about the quality and usefulness of pixels, the Quality Assessment (QA) metadata objects, and QA science data sets (SDS) that can provide information of VI quality assurance.
In addition, we investigated the departure NDVI (NDVI d ) series instead of the original NDVI series for the purpose of eliminating the effects of seasonal variations in the NDVI time series. The departure data NDVI d are widely used in the satellite data processing and filtering [19,[35][36][37]. The NDVI d is calculated, as follows: where µ NDVI is the average of all samples in all available years that were measured in the same calendar day.

Study Sites
The Camp Fire occurred in California, USA on November 8, 2018 with the initial coordinates of 39.813 • N and 121.435 • W. It was the deadliest and most destructive wildfire in California history, being caused by the Pacific Gas and Electric Co. (PG&E) transmission lines and finally killing 85 people [38]. We analyzed the VIIRS NDVI time series from 2012 to 2019 of four study sites to characterize the time dynamics and the fire-affected status of vegetation; two sites were affected by the Camp Fire (L1 and L3), while the other two were not affected (L2 and L4) and located outside the burned area, as shown in Figure 1. For each study site, we extracted 3 × 3 adjacent pixels (P1 to P9). The central coordinates of the four study sites are: L1 (39. [39,40]). Programme (IGBP) classification scheme generated from the MODIS Land Cover Product MCD12Q1 [39,40]).

Multifractal Detrended Fluctuation Analysis
The Multifractal Detrended Fluctuation Analysis (MFDFA) [22] is a well-known method for analyzing the multifractal features of non-stationary signals, like the vegetational series. Given the series x(i) where i = 1, 2, ..., N, where N is its size, xave its mean, the profile y(i) is calculated by integration: The profile y(i) is then divided into int( ) non-overlapping boxes of equal size m.
Becasue N is usually not a multiple of scale m, a small portion of the series at the end would be left. The same data division is performed from the end of the series, thereby finally generating 2Nm segments, in order not to ignore this portion. For each of the 2Nm segments, the polynomial local trend is calculated by a least square fit, and the following variance is obtained: where yν(i) is the fitting polynomial of p-th degree in the segment ν. The fitting polynomial of degree p eliminates all of the trends in the profile of order up to p, thus up to p − 1 in the original series [22]. Subsequently, the q-th order fluctuation function Fq(m) is computed averaging over all the segments, as follows: Figure 1. Location of the four study sites (L1, L2, L3, and L4) and the Camp Fire (dark red area).

Multifractal Detrended Fluctuation Analysis
The Multifractal Detrended Fluctuation Analysis (MFDFA) [22] is a well-known method for analyzing the multifractal features of non-stationary signals, like the vegetational series. Given the series x(i) where i = 1, 2, . . . , N, where N is its size, x ave its mean, the profile y(i) is calculated by integration: The profile y(i) is then divided into N m = int( N m ) non-overlapping boxes of equal size m. Becasue N is usually not a multiple of scale m, a small portion of the series at the end would be left. The same data division is performed from the end of the series, thereby finally generating 2N m segments, in order not to ignore this portion. For each of the 2N m segments, the polynomial local trend is calculated by a least square fit, and the following variance is obtained: where y ν (i) is the fitting polynomial of p-th degree in the segment ν. The fitting polynomial of degree p eliminates all of the trends in the profile of order up to p, thus up to p − 1 in the original series [22]. Subsequently, the q-th order fluctuation function F q (m) is computed averaging over all the segments, as follows: In the fluctuation function F q (m), the moment order q can assume any real value other than zero, while for q = 0, the fluctuation function is calculated by a logarithmic averaging procedure: The large fluctuations of signal are enhanced by the positive values of q, while the small fluctuations by the negative values. F q (m) increases with m as a power-law if the series is characterized by long-range correlations: Plotting F q (m) on log-log scales and fitting it by a least square method versus m, the slope of the fitting line is the so-called generalized Hurst exponent h q , which represents the scaling properties of the series at a certain moment order q. In particular, for stationary series, h 2 is the Hurst exponent H. If the series is monofractal, the exponent h q are nearly constant with q, which suggests that the scaling behavior of the small and large variations is approximately identical; however, if the series is multifractal, the exponent h q decreases with q, which suggests that the small and large fluctuations scale differently and the series has a more complex structure.

Fisher-Shannon Analysis
The Fisher-Shannon (FS) method permits the analysis of complex time series by jointly using the Fisher Information Measure (FIM) and the Shannon entropy (SE). The FIM and SE describe the characteristics of the probability density function of a series, respectively, at a local and global scale [24], and are commonly used to investigate the complexity of non-stationary time series in terms of order and organization (FIM) and disorder and uncertainty (SE) [41]. The FIM and SE are calculated, as follows: where f (x) is the probability density function of the series x. Because SE can also be negative, the exponential transformation of Shannon entropy is generally applied to obtain the so-called Shannon entropy power N X that is commonly utilized in statistical analysis: According to the isoperimetric inequality FIM · N X ≥ D [42], where D is the dimension of the space (1 in case of time series), the FIM and the N X are interrelated, which suggests that a better description of the dynamics of a time series would be given by using jointly both the measures. The equality stands in case of Gaussian processes. The use of both measures provides a more complete description of the time series, since the FIM focuses on the local properties of the probability density function, while the Shannon entropy on the global properties. It was also shown in Martin et al. (1999) [43] that FIM allowed for the detection of some non-stationary behavior in situations where the Shannon entropy showed a limited dynamics.
As the calculation of FIM and N X depends on the probability density function, attention has to be paid to its good estimation. In this study, we used the kernel-based approach to estimate f (x), which has been shown to have a better performance than the discrete-based approach in calculating the value of FIM and SE for the Gaussian distributed series [44]. The kernel-based approach for estimating the probability density function is based on the kernel density estimator technique [45,46]: where b refers to the bandwidth, M represents the number of data, and K(u) is a continuous non-negative and symmetric kernel function that satisfies the following two conditions: The estimation of f (x) uses an optimized integrated method [43] that is based on Troudi et al.'s [47] and Raykar and Duraiswami's [48] algorithms, where a Gaussian kernel with zero mean and unit variance is adopted:f The isoperimetric inequality enables the application of the Fisher-Shannon (FS) information plane to explore the dynamics of a series [49], in which the coordinate axes are N X and FIM. For scalar signals, the line FIM·N X = 1 divides the FS information plane into two parts, and each signal is represented by a point that lies exclusively in the half-space of FIM·N X > 1.

Results
We analyzed four study sites, among which two were affected by the Camp Fire (L1 and L3), and two were not affected (L2 and L4), as they were located outside the area burned by the wildfire. According to the data pre-processing procedure of outlier filtering described in Section 2.1, very few outliers were omitted for each site, in particular 2.54% (L1), 2.54% (L2), 1.98% (L3), and 1.98% (L4). L1 and L2 were covered by evergreen needleleaf forests, while L3 and L4 by woody savannas. For each study site we extracted the NDVI time series of 3×3 adjacent pixels (P1 to P9). For each pixel, we calculated NDVI d that we analyzed by the MFDFA and the FS method. Figure 2 shows, as an example, the NDVI and the NDVI d of some pixels of each site. Although the NDVI of needleleaf forests (Figure 2a,c) appears to be higher than that of woody savannas (Figure 2e,g), however the departure NDVI d of both vegetation covers do not show significant difference, thus confirming the effectiveness of the seasonal removal by the procedure of Equation (2). Furthermore, even by visual inspection, it is clearly seen that both NDVI and NDVI d drop sharply down at the time of Camp Fire occurrence, towards increase, then, gradually through time to the pre-fire conditions. Entropy 2020, 22, x 6 of 16 where b refers to the bandwidth, M represents the number of data, and K(u) is a continuous non-negative and symmetric kernel function that satisfies the following two conditions: The estimation of f(x) uses an optimized integrated method [43] that is based on Troudi et al.'s [47] and Raykar and Duraiswami's [48] algorithms, where a Gaussian kernel with zero mean and unit variance is adopted: The isoperimetric inequality enables the application of the Fisher-Shannon (FS) information plane to explore the dynamics of a series [49], in which the coordinate axes are NX and FIM. For scalar signals, the line FIM·NX = 1 divides the FS information plane into two parts, and each signal is represented by a point that lies exclusively in the half-space of FIM·NX > 1.

Results
We analyzed four study sites, among which two were affected by the Camp Fire (L1 and L3), and two were not affected (L2 and L4), as they were located outside the area burned by the wildfire. According to the data pre-processing procedure of outlier filtering described in Section 2.1, very few outliers were omitted for each site, in particular 2.54% (L1), 2.54% (L2), 1.98% (L3), and 1.98% (L4). L1 and L2 were covered by evergreen needleleaf forests, while L3 and L4 by woody savannas. For each study site we extracted the NDVI time series of 3×3 adjacent pixels (P1 to P9). For each pixel, we calculated NDVId that we analyzed by the MFDFA and the FS method. Figure 2 shows, as an example, the NDVI and the NDVId of some pixels of each site. Although the NDVI of needleleaf forests (Figure 2a,c) appears to be higher than that of woody savannas (Figure 2e,g), however the departure NDVId of both vegetation covers do not show significant difference, thus confirming the effectiveness of the seasonal removal by the procedure of Equation (2). Furthermore, even by visual inspection, it is clearly seen that both NDVI and NDVId drop sharply down at the time of Camp Fire occurrence, towards increase, then, gradually through time to the pre-fire conditions. In order to select the optimal value of p, we calculated and compared the fluctuation functions obtained by using different degrees p since the MFDFA depends on the degree of the detrending polynomial p. Since the detrending of the time series is performed by subtracting the polynomial from the profile, different orders of DFA differ in their possibility to eliminate trends in the series. In MFDFA of order p, trends of order p are eliminated from the profile. Thus, comparing the fluctuation functions for different orders p allows for one to estimate the type of the polynomial trend in the time series and to select the optimal one, which is the order that is enough to remove those trends in the series. Figure 3 shows, as an example, the fluctuation functions for two values of the moment order q, −5 ( Figure 3a) and 5 (Figure 3b), for the pixel P8 of site L1, for different orders of p from 1 to 5. Since the fluctuation functions tend to overlap for value of p higher than 4 in both cases, the optimal In order to select the optimal value of p, we calculated and compared the fluctuation functions obtained by using different degrees p since the MFDFA depends on the degree of the detrending polynomial p. Since the detrending of the time series is performed by subtracting the polynomial from the profile, different orders of DFA differ in their possibility to eliminate trends in the series. In MFDFA of order p, trends of order p are eliminated from the profile. Thus, comparing the fluctuation functions for different orders p allows for one to estimate the type of the polynomial trend in the time series and to select the optimal one, which is the order that is enough to remove those trends in the series.  (Figure 3b), for the pixel P8 of site L1, for different orders of p from 1 to 5. Since the fluctuation functions tend to overlap for value of p higher than 4 in both cases, the optimal degree of the detrending polynomial is 4; thus, the multifractal analysis was performed by using a 4th degree detrending polynomial. degree of the detrending polynomial is 4; thus, the multifractal analysis was performed by using a 4 th degree detrending polynomial.
(a) (b) The multifractality of a time series is revealed by the power-law behavior of the fluctuation functions for any moment order q. Plotting the fluctuation function in log-log scales, the multifractality of the series is indicated by the linear behavior of the fluctuation functions, whose slope, estimated by using a least square method, is the so-called generalized Hurst exponent hq. Figure 4 shows as an example the fluctuation functions of the NDVId of pixel P1 of site L3 for q ranging from −5 to 5.  Figure 5 shows the generalized Hurst exponents hq that were calculated for all NDVId series. Most of the hq decrease with the increase of the moment order q for all of the pixels; however, for a few pixels (P1, P2, P3, and P8 of site L4), hq fluctuates around the average value. We applied the MFDFA to one hundred shuffles for each NDVId series and calculated their hq in order to recognize the source of multifractality (long-range correlations or type of distribution). Since the shuffling just removes the correlations but preserves the distribution of the series, we can discriminate between the multifractality due to the existence of long-range correlations (the hq of the shuffles fluctuates around the value of 0.5) or to the probability density function of the series (the hq of the shuffles behave in nearly the same manner as those of the original series). Figure 6 shows, as an example, the generalized Hurst exponents hq (black dots) of NDVId series and the 95% confidence band (red dotted lines) of the hq of the shuffles for the pixel P7 of site L1; we can see that the hq of the original series are well beyond the 95% confidence band of the shuffles, whose generalized Hurst exponents, however, slightly decrease with the moment order q; this indicates that the multifractality of the The multifractality of a time series is revealed by the power-law behavior of the fluctuation functions for any moment order q. Plotting the fluctuation function in log-log scales, the multifractality of the series is indicated by the linear behavior of the fluctuation functions, whose slope, estimated by using a least square method, is the so-called generalized Hurst exponent h q . Figure 4 shows as an example the fluctuation functions of the NDVI d of pixel P1 of site L3 for q ranging from −5 to 5. The multifractality of a time series is revealed by the power-law behavior of the fluctuation functions for any moment order q. Plotting the fluctuation function in log-log scales, the multifractality of the series is indicated by the linear behavior of the fluctuation functions, whose slope, estimated by using a least square method, is the so-called generalized Hurst exponent hq. Figure 4 shows as an example the fluctuation functions of the NDVId of pixel P1 of site L3 for q ranging from −5 to 5.  Figure 5 shows the generalized Hurst exponents hq that were calculated for all NDVId series. Most of the hq decrease with the increase of the moment order q for all of the pixels; however, for a few pixels (P1, P2, P3, and P8 of site L4), hq fluctuates around the average value. We applied the MFDFA to one hundred shuffles for each NDVId series and calculated their hq in order to recognize the source of multifractality (long-range correlations or type of distribution). Since the shuffling just removes the correlations but preserves the distribution of the series, we can discriminate between the multifractality due to the existence of long-range correlations (the hq of the shuffles fluctuates around the value of 0.5) or to the probability density function of the series (the hq of the shuffles behave in nearly the same manner as those of the original series). Figure 6 shows, as an example, the generalized Hurst exponents hq (black dots) of NDVId series and the 95% confidence band (red dotted lines) of the hq of the shuffles for the pixel P7 of site L1; we can see that the hq of the original series are well beyond the 95% confidence band of the shuffles, whose generalized Hurst exponents, however, slightly decrease with the moment order q; this indicates that the multifractality of the  Figure 5 shows the generalized Hurst exponents h q that were calculated for all NDVI d series.
Most of the h q decrease with the increase of the moment order q for all of the pixels; however, for a few pixels (P1, P2, P3, and P8 of site L4), h q fluctuates around the average value. We applied the MFDFA to one hundred shuffles for each NDVI d series and calculated their h q in order to recognize the source of multifractality (long-range correlations or type of distribution). Since the shuffling just removes the correlations but preserves the distribution of the series, we can discriminate between the multifractality due to the existence of long-range correlations (the h q of the shuffles fluctuates around the value of 0.5) or to the probability density function of the series (the h q of the shuffles behave in nearly the same manner as those of the original series). Figure 6 shows, as an example, the generalized Hurst exponents h q (black dots) of NDVI d series and the 95% confidence band (red dotted lines) of the h q of the shuffles for the pixel P7 of site L1; we can see that the h q of the original series are well beyond the 95% confidence band of the shuffles, whose generalized Hurst exponents, however, slightly decrease with the moment order q; this indicates that the multifractality of the NDVI d series could be due to both the presence of long-range correlations in the series and the distribution of the values.  We use the range of the generalized Hurst exponent (hq-range) that is the difference between the maximum and the minimum generalized Hurst exponent to quantify the multifractality. Figure 7 shows the hq-range of all the analyzed NDVId series. From a visual inspection, for evergreen needleleaf forests, the hq-range of burned and unburned sites is approximately the same (Figure 7a), while a certain difference seems to appear in the case of savannas (Figure 7b), indicating that the wildfire could have changed the multifractal characteristics of savannah, but not for evergreen needleleaf vegetation. In order to check whether the averages of hq, <hq>, are significantly different, we applied the two sample Student's t-test to the group of hq calculated for L1 and L2 and for L3 and   We use the range of the generalized Hurst exponent (hq-range) that is the difference between the maximum and the minimum generalized Hurst exponent to quantify the multifractality. Figure 7 shows the hq-range of all the analyzed NDVId series. From a visual inspection, for evergreen needleleaf forests, the hq-range of burned and unburned sites is approximately the same (Figure 7a), while a certain difference seems to appear in the case of savannas (Figure 7b), indicating that the wildfire could have changed the multifractal characteristics of savannah, but not for evergreen needleleaf vegetation. In order to check whether the averages of hq, <hq>, are significantly different, we applied the two sample Student's t-test to the group of hq calculated for L1 and L2 and for L3 and We use the range of the generalized Hurst exponent (h q -range) that is the difference between the maximum and the minimum generalized Hurst exponent to quantify the multifractality. Figure 7 shows the h q -range of all the analyzed NDVI d series. From a visual inspection, for evergreen needleleaf forests, the hq-range of burned and unburned sites is approximately the same (Figure 7a), while a certain difference seems to appear in the case of savannas (Figure 7b), indicating that the wildfire could have changed the multifractal characteristics of savannah, but not for evergreen needleleaf vegetation. In order to check whether the averages of h q , <h q >, are significantly different, we applied the two sample Student's t-test to the group of h q calculated for L1 and L2 and for L3 and L4. We found that at 95% confidence, L1 and L2 are not significantly distinguishable (p-value = 0.381106), while L3 and L4 are significantly distinguishable (p-value = 0.010356), which indicates that the fire-affected and fire-unaffected sites are distinguishable for savanna vegetation covers, and not for evergreen needleleaf forests.   Figure 8 shows the generalized Hurst exponents and their 95% confidence band based on one hundred shuffles. Although the FGN is a monofractal process, its simulations show a slight multifractality that is probably due to the short length. However, their hq-range is 0.15 for FGN03, 0.18 for FGN05, and 0.16 for FGN08; the small value of the hq-range confirms the monofractal character of the three time series. Furthermore, the generalized Hurst exponents of FGN03 and FGN08 are well beyond the 95% confidence band that is based on the shuffles, and this strengthens the robustness of the MFDFA in detecting multifractality, even in short time series. Although the pixels extracted for each site are spatially located close to each other, a certain local variability could still exists that might influence the multifractal character of the NDVId; thus, we also analyzed the mean NDVId, <NDVId>, which is the mean of NDVId over the nine pixels for each site, and that might represent the behavior of vegetation in a site in a more global manner. We calculated the generalized Hurst exponents for three simulations of Fractional Gaussian Noise (FGN) with three different values of the Hurst exponent (0.3, 0.5, and 0.8), indicated respectively as FGN03, FGN05, and FGN08, and same length as our time series, in order to check the robustness of the results against the short length of the time series (395 values). Figure 8 shows the generalized Hurst exponents and their 95% confidence band based on one hundred shuffles. Although the FGN is a monofractal process, its simulations show a slight multifractality that is probably due to the short length. However, their h q -range is 0.15 for FGN03, 0.18 for FGN05, and 0.16 for FGN08; the small value of the h q -range confirms the monofractal character of the three time series. Furthermore, the generalized Hurst exponents of FGN03 and FGN08 are well beyond the 95% confidence band that is based on the shuffles, and this strengthens the robustness of the MFDFA in detecting multifractality, even in short time series.
Entropy 2020, 22, x 10 of 16 L4. We found that at 95% confidence, L1 and L2 are not significantly distinguishable (p-value = 0.381106), while L3 and L4 are significantly distinguishable (p-value = 0.010356), which indicates that the fire-affected and fire-unaffected sites are distinguishable for savanna vegetation covers, and not for evergreen needleleaf forests. We calculated the generalized Hurst exponents for three simulations of Fractional Gaussian Noise (FGN) with three different values of the Hurst exponent (0.3, 0.5, and 0.8), indicated respectively as FGN03, FGN05, and FGN08, and same length as our time series, in order to check the robustness of the results against the short length of the time series (395 values). Figure 8 shows the generalized Hurst exponents and their 95% confidence band based on one hundred shuffles. Although the FGN is a monofractal process, its simulations show a slight multifractality that is probably due to the short length. However, their hq-range is 0.15 for FGN03, 0.18 for FGN05, and 0.16 for FGN08; the small value of the hq-range confirms the monofractal character of the three time series. Furthermore, the generalized Hurst exponents of FGN03 and FGN08 are well beyond the 95% confidence band that is based on the shuffles, and this strengthens the robustness of the MFDFA in detecting multifractality, even in short time series. Although the pixels extracted for each site are spatially located close to each other, a certain local variability could still exists that might influence the multifractal character of the NDVId; thus, we also analyzed the mean NDVId, <NDVId>, which is the mean of NDVId over the nine pixels for each site, and that might represent the behavior of vegetation in a site in a more global manner. Although the pixels extracted for each site are spatially located close to each other, a certain local variability could still exists that might influence the multifractal character of the NDVI d ; thus, we also analyzed the mean NDVI d , <NDVI d >, which is the mean of NDVI d over the nine pixels for each site, and that might represent the behavior of vegetation in a site in a more global manner. Figure 9 shows the generalized Hurst exponents h q and the h q -range of the <NDVI d > of the four sites. A general pattern seems to be evidenced from the analysis of <NDVI d >; in fact, for both the evergreen needleleaf forests and woody savannas, the <NDVI d > is characterized by a value of h q -range for burned sites (L1 and L3) larger than that of unburned sites (L2 and L4), although for savanna covers the difference between the burned and unburned sites is much larger.  Figure 9 shows the generalized Hurst exponents hq and the hq-range of the <NDVId> of the four sites. A general pattern seems to be evidenced from the analysis of <NDVId>; in fact, for both the evergreen needleleaf forests and woody savannas, the <NDVId> is characterized by a value of hq-range for burned sites (L1 and L3) larger than that of unburned sites (L2 and L4), although for savanna covers the difference between the burned and unburned sites is much larger.
The Fisher-Shannon method was applied to all of the investigated NDVId series and the results are shown in the Fisher-Shannon information plane ( Figure 10). It can be seen that a clear separation exists between the burned and unburned sites for both vegetation covers. The two sample Student's t-test for the NX and FIM indicates a significant difference between the two groups of the obtained parameters of burned and unburned sites (p-value < 0.05) ( Table 1). The Fisher-Shannon method was applied to all of the investigated NDVI d series and the results are shown in the Fisher-Shannon information plane ( Figure 10). It can be seen that a clear separation exists between the burned and unburned sites for both vegetation covers. The two sample Student's t-test for the N X and FIM indicates a significant difference between the two groups of the obtained parameters of burned and unburned sites (p-value < 0.05) ( Table 1).  Figure 9 shows the generalized Hurst exponents hq and the hq-range of the <NDVId> of the four sites. A general pattern seems to be evidenced from the analysis of <NDVId>; in fact, for both the evergreen needleleaf forests and woody savannas, the <NDVId> is characterized by a value of hq-range for burned sites (L1 and L3) larger than that of unburned sites (L2 and L4), although for savanna covers the difference between the burned and unburned sites is much larger.
The Fisher-Shannon method was applied to all of the investigated NDVId series and the results are shown in the Fisher-Shannon information plane ( Figure 10). It can be seen that a clear separation exists between the burned and unburned sites for both vegetation covers. The two sample Student's t-test for the NX and FIM indicates a significant difference between the two groups of the obtained parameters of burned and unburned sites (p-value < 0.05) ( Table 1).

Discussion
The NDVI is one of most widely used vegetation indices for characterizing the vegetation status. The NDVI can also be used to discriminate between different types of vegetation covers; in our case, the variability of NDVI is smaller for evergreen needleleaf forests and larger for woody savannas. The NDVI is very sensitive to situations where the vegetation is stressed, like when being affected by a wildfire. In this study, we analyzed the vegetation dynamics of two different vegetation covers to investigate the impact of fire. We used two different methods, the MFDFA and the FS method, obtaining a set of parameters (hq-range, FIM and N X ), by which the discrimination between burned and unburned sites was investigated.
The h q -range is used to quantify the degree of multifractality of a series. A larger multifractality indicates that the series is characterized by a larger heterogeneity that means the large fluctuations of the series scale differently from the small fluctuations since the MFDFA highlights the scaling behavior of the time dynamics of a series. The NDVI of unburned sites of savannas type are featured by a relatively small value of h q -range, which indicated that the NDVI is rather monofractal and the scaling behavior of the large and small variations is nearly the same; the vegetation is probably only subjected to fluctuations induced by the seasonal and climatic effects. The NDVI of fire-affected sites of woody savannas covers appear, instead, more heterogeneous, suggesting that small and large fluctuations scale differently; the wildfire impacted on this vegetation that responded to such stress by intermittently fluctuating and thus increasing its multifractality. Such a difference between burned and unburned sites covered by woody savannas becomes much more evident when analyzing the <NDVI d >, which might more globally represent the status of vegetation of a site, filtering out the local variability that could affect single pixels. The <NDVI d > of woody savannas, in fact, allows for very clearly discriminating between burned and unburned sites, with the last being signaled by a lower multifractality degree. Additionally, the <NDVI d > of unburned sites of evergreen needleleaf forests is characterized by a relatively lower multifractalty degree than burned sites, although no significant difference has been found between the single pixels of fire-affected and fire-unaffected sites. Such not so clear multifractal difference between burned and unburned sites that are covered by evergreen needleleaf forests could be due to the fact that fire mainly affects the grass, whose vegetation fine texture makes it more sensitive to the burning than is true of leaves and stems of forests, as was observed in [50].
The FS information plane allows for investigating the complex structure of a time series in terms of organization/order. Generally speaking, the bottom right part of the FS plane is occupied by series characterized by less ordered and less organized series, while the top left part of it by more ordered and more organized series. In the FS information plane, the pixels of the burned sites are well discriminated from those of the unburned sites for both evergreen forests and woody savannas. In particular, pixels of burned sites have smaller N X and larger FIM values, whereas those of the unburned sites are characterized by larger N X and smaller FIM in agreement with previous researches [21]. A larger FIM and a smaller N X quantifies a larger organization and less disorder of the series, where, for larger organization and less disorder, we mean that the distribution of the values of the series is rather peaked, and some values are more frequent than others. Since the vegetation undergoes the disturbance of wildfire, the sudden decrease of NDVI soon after the wildfire and its gradual increase during the post-fire vegetation recovery makes the distribution of the values less uniform than that characterizing the NDVI of the unburned sites; thus, the vegetation does not uniformly fluctuate and the FIM and N X are larger and smaller, respectively.

Conclusions
The Suomi-NPP VIIRS NDVI series of two types of vegetation covers were analyzed using the multifractal detrended fluctuation analysis and the Fisher-Shannon method. Pixel time series of burned and unburned sites were investigated. Several parameters, in particular the range of the generalized Hurst exponent (h q -range), Shannon entropy power (N X ), and Fisher Information Measure (FIM), were utilized to measure the degree of heterogeneity and organization of the vegetation series in conjunction with their separating ability to discriminate fire-affected and fire-unaffected pixels.
In general, the NDVI of the fire-affected sites are more heterogeneous in comparison to that of the fire-unaffected sites, especially for woody savannas. In addition, the vegetation dynamics of fire-affected sites are characterized by a larger organization and smaller disorder, due to the less uniform behavior that characterizes vegetation when it is stressed by a fire and during its post-fire recovery process.
Furthermore, the FS information plane has revealed its potential in being sensitive to wildfire impact on different types of vegetation more than MFDFA; in fact, the MFDFA has allowed for discriminating the vegetation dynamics between the burned and unburned sites more significantly in woody savannas covers than in evergreen needleleaf forests.

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