Towards Operational Detection of Forest Ecosystem Changes in Protected Areas

Abstract

This paper discusses the application of the Cross-Correlation Analysis (CCA) technique to multi-spatial resolution Earth Observation (EO) data for detecting and quantifying changes in forest ecosystems in two different protected areas, located in Southern Italy and Southern India. The input data for CCA investigation were elaborated from the forest layer extracted from an existing Land Cover/Land Use (LC/LU) map (time T₁) and a more recent (T₂, with T₂ > T₁) single date image. The latter consist of a High Resolution (HR) Landsat 8 OLI image and a Very High Resolution (VHR) Worldview-2 image, which were analysed separately. For the Italian site, the forest layer (1:5000) was first compared to the HR Landsat 8 OLI image and then to the VHR Worldview-2 image. For the Indian site, the forest layer (1:50,000) was compared to the Landsat 8 OLI image then the changes were interpreted using Worldview-2. The changes detected through CCA, at HR only, were compared against those detected by applying a traditional NDVI image differencing technique of two Landsat scenes at T₁ and T_2. The accuracy assessment, concerning the change maps of the multi-spatial resolution outputs, was based on stratified random sampling. The CCA technique allowed an increase in the value of the overall accuracy: from 52% to 68% for the Italian site and from 63% to 82% for the Indian site. In addition, a significant reduction of the error affecting the stratified changed area estimation for both sites was obtained. For the Italian site, the error reduction became significant at VHR (±2 ha) in respect to HR (±32 ha) even though both techniques had comparable overall accuracy (82%) and stratified changed area estimation. The findings obtained support the conclusions that CCA technique can be a useful tool to detect and quantify changes in forest areas due to both legal and illegal interventions, including relatively inaccessible sites (e.g., tropical forest) with costs remaining rather low. The data obtained through CCA intervention could not only support the commitments undertaken by the European Habitats Directive (92/43/EEC) and the Convention of Biological Diversity (CBD) but also satisfy UN Sustainable Development Goals (SDG).

1. Introduction

1.1. Change Detection Techniques

The present study discusses the application of the Cross Correlation Analysis (CCA) change detection technique, which can detect changes at different spatial scales using a Land Cover/Use (LC/LU) map and a sole recent image for forest monitoring. The approach proposed may reduce both computational and imagery costs as well as require fewer in-field campaigns.

Earth Observation (EO) data and techniques are most promising for monitoring and quantifying forest changes at multiple scales and high frequencies [1,2,3,4,5]. These techniques can provide new products and services for a wide user community including ecologists and decision makers such as those involved in the commitments of Natura 2000 site conservation [6,7,8,9].

Change in forest extent and forest fragmentation has been identified as one of the most important drivers of forest ecosystem services loss and is thus one of the most critical factors to monitor. The selection of the most appropriate change detection technique will depend on the data available about the study area. Actually, when Land Cover/Land use (LC/LU) maps at time T₁ and T₂, (with T₂ > T₁) are not available, the techniques chosen are usually based on direct comparison of calibrated and co-registered image pairs acquired at T₁ and T₂ [10]. In some applications, pixel spectral values in the two images or derived spectral features (e.g., normalized difference vegetation index (NDVI)) have been used in the comparison. In such cases, no semantic information about specific class transitions/modifications, useful for identifying pressures on the area [5,9], can be automatically extracted. In other applications, spectral features are used, involving a higher processing image level. This approach allows detection of changes that may have real-world meaning such as a decrease/increase of the NDVI vegetation index. Whether using the former or the latter technique, in-field campaigns and/or visual inspection of VHR images or aerial orthophoto are to be used to identify specific class transitions/modifications.

Otherwise, when two thematic maps (e.g., LC/LU or habitats maps) independently produced at time T₁ and time T₂ are available, the well-known Post Classification Comparison (PCC) approach [10,11,12,13,14] can be used. The degree of success of this technique depends upon the reliability of the input thematic maps [15,16] since the quality of the output change image is related to the product of the accuracies of the two maps being compared.

Although the two mentioned change detection techniques are very useful, their scale of application depends upon the investigation purposes. For regional decision making, data and maps from high spatial resolution sensors, such as the Landsat series or the European Space Agency’s (ESA) Sentinels optical and radar, appear to be the most appropriate since they can cover greater area and allow more frequent coverage. For more local decision-making, finer scale data and maps such as those provided by airborne/spaceborne Very High Resolution (VHR) sensors (e.g., QuickBird, GeoEye, Worldview-2/3) are being employed. However, since image acquisition must be tasked, operational change detection at local scale can be seriously hampered by lack of their repeated acquisition [6,7,8,9].

1.2. Change Detection with VHR Images

Due to the complexity of class description at fine scales and the huge image information content, change detection based on VHR images is viewed as more difficult to automate compared to coarser spatial resolution data. Thus object based classification techniques are generally suggested [13,17]. Nevertheless, the process still remains computationally expensive when the entire set of thematic classes in the two T₁ and T₂ maps must be processed. Chen et al. [13] suggested applying a stratified change detection approach by considering only one target class of interest (e.g., evergreen forest) at a time in the map. An additional difficulty arises from the high cost of VHR images. This issue may become critical when multi-seasonal VHR acquisitions need to be used in order to produce LC maps having large Overall Accuracy (OA) value and low error rate [18].

1.3. CCA Based Scenario Technique at VHR and HR

The Cross Correlation Analysis (CCA) is a change detection technique developed by the American company Earthsat, Inc. and is fruitful in evaluating the differences between an existing LC/LU map (T₁) and a recent single-date multispectral image (T₂) [19,20,21]. A (CCA) based scenario will be used in this study to evaluate changes at multiple-scale. Our approach includes a pre-existent LC/LU map available at regional/local scales, obtained by visual inspection of orthophotos and validated by in-field campaigns, will be used as reference image at time T₁ in the change detection process. The image T₂ pixels corresponding to the target class in the T₁ map will be analysed by CCA to detect changes at T₂ [19] at both VHR and HR. Places of transitions from the target class can then be identified without a complete classification process at time T₂. Even though the new class at T₂ remains unknown it may be determined by in-field inspection or by visual interpretation of VHR images, when available.

CCA applications to High Resolution (HR) (e.g., Landsat TM) and Medium Resolution (MR) imagery (e.g., MERIS) have been reported [19,20]. More recently, VHR (e.g., WorldView-2) applications for grassland ecosystems changes have been proposed [21]. This study argues that the CCA technique can be attractive for fine scale change detection, since it can reduce change detection costs when: (a) the acquisition of several (multi-seasonal) VHR images at time T₂ (e.g., within a year) is too expensive; and (b) no archival VHR data are available at T₁ for direct image comparison between the T₁ image and a new tasked T₂ image (with T₁ < T₂).

The present study aims at demonstrating: (i) the effectiveness of the CCA techniques for forest change detection at both HR and VHR on two protected areas in different bio-geographical regions, i.e., Southern Italy and Southern India; and (ii) evaluating the CCA accuracy compared to an unsupervised traditional change detection technique based on image differencing of spectral features from available image pair (T₁ and T₂) at HR only. NDVI will be used as suitable spectral feature to detect changes in the images. Estimates of uncertainty in the discrimination of classes (i.e., change and no-change) will be provided based on stratified sampling of ground truth reference samples for validation and the recommendations in [15,22]. The proposed approach could speed up operational forest change detection as well as prove to be cost effective. The findings encourage speculations about the opportunity to implement CCA technique for compliance with both European and United Nation (UN) commitments.

2. Materials and Methods

2.1. Site Description

The CCA analysis and comparison of change detection techniques have been carried out on data obtained from two protected areas, belonging to different biogeographical regions, one located in India and the other in Southern Italy.

2.1.1. Italian Site

The first study site consists of 500 km² located in Apulia, a region of South East Italy, included in the Natura 2000 “Murgia Alta” site (SCI/SPA IT9120007, of the European Union Habitat Directive 92/43/EEC and Bird Directive 147/2009/EC) (Figure 1a). The land consists of a calcareous upland covered by almost 24% semi-natural dry grasslands. The site represents one of the most important areas for the conservation of this ecosystem type in Europe. Highly fragmented evergreen and deciduous forest patches represent what remains of the old natural vegetation coverage otherwise completely destroyed by human activity over time. Agricultural intensification, urbanization, arson and land abandonment represent the main pressure factors on the area biodiversity [23].

For T₁ (2006), an existing LC/LU map (1:5000) based on visual orthophoto interpretation and validated (85% overall accuracy) by in-field campaigns, was available for the site. The authors used this map, originally produced in CORINE and subsequently translated into the FAO-LCCS taxonomy [24], used as image reference (Figure 1b).

To detect changes possibly due to arson, multi-spatial resolution summer images were considered at T₂, according to the procedure described in [19,21]. In particular, a VHR WorldView-2 single date image (spatial resolution 2 m) acquired in July 2012 was considered (Figure 1c). The image was orthorectified in the WGS84/UTM33N projection (EPSG code: 32633), with a RMSE less than one pixel, and calibrated to Top of Atmosphere (TOA) reflectance values. The image was provided at no cost by the European Space Agency (ESA) under the Data Warehouse 2011–2014 policy within the FP7-SPACE BIO_SOS project [25].

For change detection at HR, two images were downloaded from the US Geological Survey [26], i.e., a recent Landsat 8 OLI image [27] acquired in August 2013 (as T₂), and an earlier Landsat 7 ETM+ image dated July 2006 (as T₁). Both images were orthorectified in the WGS84/UTM43N projection, with a RMSE less than one pixel and coregistered. Relative pairwise pixel-based radiometric normalization [28] was applied to the two images in order to reduce disparities in the image acquisition conditions. The Landsat 7 ETM+ image of July 2006 (T₁) was used as reference image for co-registration and radiometric normalization.

2.1.2. Indian Site

The Indian site includes a 540 km² Tiger reserve located in the Western Ghats biodiversity hotspot in Southern India, named “Biligiri Rangaswamy Temple Tiger Reserve” (Figure 2a). The site has a heterogeneous physiography, with hills running in the north–south direction, and elevations ranging from 600 to 1800 m above sea level. The area biophysical conditions, with rainfall in two different seasons allow a distinctive ecosystem to thrive in this area and support a diversity of endemic flora and fauna. The only available map of the region is dated 1997, scale (1:50,000), and includes ten different vegetation types. These range from dry scrub forest to dense wet evergreen forests in the higher elevation areas (Figure 2b). The evergreen forest patches (Figure 2c) are found both in contiguous areas and in dense patches among a mosaic of high elevation grassland areas.

In the present paper, only the wet evergreen forest patches shown in the map were considered as input T₁ layer to CCA. In addition, a Worldview-2 image acquired in March 2013 (summer season) and provided by ESA, within the FP7-SPACE BIO_SOS project [25], was used as VHR image at T₂. The image was orthorectified in the WGS84/UTM43N projection (EPSG code: 32643), with a RMSE less than one pixel, and calibrated to Top of Atmosphere (TOA) reflectance values.

Two additional Landsat images, dated March 1997 (Landsat 5) and March 2016 (Landsat 8 OLI), were considered for direct comparison of NDVI indices [26]. These images were orthorectified in the WGS84/UTM43N projection, and coregistered with a RMSE less than one pixel. Relative pairwise pixel-based radiometric normalization [28] was applied to the two images in order to reduce disparities in the image acquisition conditions. The March 1997 Landsat 5 image was used as reference image.

2.2. Methods

To detect changes in the forest target class of the two protected areas, two experiments were carried out:

• One detected changes at HR comparing the NDVI images from the two Landsat images (Image to Image comparison) acquired at T₁ and T₂, with T₂ > T₁.
• The other detected changes at HR and VHR by Cross Correlation Analysis (Map to Image comparison).

Table 1. Set of experiments and abbreviations used for the Italian site.

Experiment	Input Data at T1	Time T2	Change Method	Change Method Acronym
1	Landsat 7 ETM image: 27 July 2006	Landsat 8 OLS image: 7 August 2013	NDVI direct comparison by image differencing	DIFF_NDVI_HR
2.a	Forest layer from existing Land Cover/Land Use Map dated 2006	Landsat 8 OLS image: 7 August 2013	Cross Correlation Analysis	CCA_HR
2.b		Worldview-2 image: 6 July 2012		CCA_VHR

and

Table 2. Set of experiments and abbreviations used for the Indian site.

Experiment	Input Data at T1	Time T2	Change Method	Change Method Acronym
1	Landsat 5 image: 16 March 1997	Landsat 8 OLS image: 20 March 2016	NDVI direct comparison by image differencing	DIFF_NDVI_HR
2.a	Tropical evergreen forest target class of interest from existing LC/LU map dated 1998	Landsat 8 OLS image: 20 March 2016	Cross Correlation Analysis	CCA_HR
2.b		Worldview-2 image: 14 March 2013		CCA_VHR

synthetize techniques and abbreviations used in the rest of this paper for the two sites, respectively.

Direct image comparison of spectral signatures or derived indexes, such as NDVI, was used when LC/LU maps are not available at T₁ and T₂. After co-registration and relative radiometric rectification of the two images acquired in the same month, of different years, can be compared at pixel level, by direct comparison of the reflectance values of each input pixels (e.g., by image differencing, image regression, change vector analysis, and image rationing) or at features level by image differencing of spectral indexes (e.g., NDVI, and PSRI) [10]. In this paper, the comparison of NDVI indices from Landsat images, acquired at T₁ and T_2, with T₂ > T₁, was used as such NDVI feature introduces some semantic related to the coverage and/or state of vegetation. As a result, the areas where changes in the forest layer occurred were identified and compared with those detected by CCA analysis at T₂.

The CCA change detection technique rests on the evaluation of the differences between an existing LC/LU map (T₁) and a recent single-date multispectral image (T₂) [19,20,21]. In this model of analysis, all pixels of the T₂ image corresponding to a specific thematic layer (target class) in the T₁ map are analysed to determine the expected reference class metrics in the T₂ image (i.e., class average spectral response and standard deviation). Following this, for each pixel in T_2, corresponding to the target class layer target layer at T₁, a statistical measure is computed to evaluate the distance between the pixel spectral signature and the reference target class metrics at T₂. Large values of such measures show evidence of the occurrence of large probable class changes. This information can be used to derive a Z-statistic for each pixel of the T₂ image falling within the target LC/LU class layer. The Z-statistic describes how close the pixel’s response is to the expected spectral response of the target class in T₂. Large Z-statistic values may identify changed pixels while small Z-statistic values may indicate no-change pixels. The Z-statistic values can be computed by Equation (1):

$$Z_{jk}=\sqrt{{\displaystyle \sum }_{i=1}^n{\left(\frac{r_{ijk-}{\mu }_{i{c}_{jk}}}{{\sigma }_{i{c}_{jk}}}\right)}^2}$$

(1)

where ${\mathit{Z}}_{\mathit{j}\mathit{k}}$ is the Z-score for a pixel jk belonging to a given class (stratum); i is the band number in the multispectral image; n is the number of bands; ${\mathit{c}}_{\mathit{j}\mathit{k}}$ is the thematic class (stratum) being considered at T₂, jk is a pixel in the stratum; ${\mathit{r}}_{\mathit{i}\mathit{j}\mathit{k}}$ is the reflectance in band i for pixel jk; ${\mathit{\mu }}_{\mathit{i}{\mathit{c}}_{\mathit{j}\mathit{k}}}$ is the mean reflectance value in band i of all pixels in a given class ${\mathit{c}}_{\mathit{j}\mathit{k}}$; and ${\mathit{\sigma }}_{\mathit{i}{\mathit{c}}_{\mathit{j}\mathit{k}}}$ is the standard deviation of the reflectance value in band i of all pixels in class ${\mathit{c}}_{\mathit{j}\mathit{k}}$.

The selection of a threshold (TH) can thus help to identify the most significant changes [19]. In the paper three values of the threshold (i.e., μ + 1σ, μ + 2σ and μ + 3σ, where μ and σ are the mean and standard deviation of the Z change image) were adopted and the best result considered. Once the location of forest changes has been located, more accurate information about the specific class transitions/modification from the target class (i.e., forest) can only be obtained by either local in-field campaigns or visual inspection of VHR imagery.

2.3. Accuracy Assessment

For the output map validation, sets of reference change and no-change forest polygons were selected through visual inspection of the available WorldView-2 images. Stratified random sampling was applied. When the sampling intensities for each class considered differed, correct calculation of the overall accuracy (OA) would require that the within-class accuracies be weigh according to the proportions of the study area corresponding to the map class [29]. Consequently, OA cannot be calculated as the sum of diagonal counts in the error matrix divided by the total count, as in the case of simple random sampling or systematic sampling design [28]. For this reason, for each experiment, the change error matrix was produced in terms of sample counts. To obtain a more accurate quantification of change OA, the protocol described in [15,22] was adopted. This protocol is based on a more informative presentation of the change error matrix and allows direct computation change accuracy and area estimates.

Given a map with q categories, when the map categories are the rows (i = 1, 2, … q) and the reference categories are the columns (j = 1, 2, … q) of the error matrix, A_tot represents the total area of the map (window), A_m,i is the mapped area (ha) of category i, in the map, and $W_i=\frac{A_{m,i}}{A_{tot}}$ is the proportion of the mapped area as category i, $\text{then }{\widehat{p}}_{ij}$ (i.e., the proportion of area for the population having map class i and reference class j in the change error matrix [22]) can be calculated as:

$${\widehat{p}}_{ij}=W_i\frac{n_{ij}}{n_{i\bullet }}$$

(2)

where $n_{ij}$ represents the sample counts and $n_{i\bullet }$ is the sum of the sample counts for the i-raw computed over the columns of the change error matrix [22].

The unbiased stratified estimator of the area of category j can be obtained as:

$${\widehat{A}}_j=A_{tot}\times {\widehat{p}}_{\bullet j}=A_{tot}{\displaystyle \sum }_i{W}_i\frac{n_{ij}}{n_{i\bullet }}$$

(3)

where ${\widehat{A}}_j$ can be considered as an “error-adjusted” estimator of the change area, because it includes the area of map omission error of category j and leaves out the area of map commission error [29].

The estimated standard error of the estimated area proportion is:

$$S\left({\widehat{p}}_{\bullet j}\right)=\sqrt{{\displaystyle \sum }_{i=1}^q{W}_i^2\frac{\frac{n_{ij}}{n_{i\bullet }}\left(1-\frac{n_{ij}}{n_{i\bullet }}\right)}{n_{i\bullet }-1}}$$

(4)

Therefore, the standard error of the stratified area estimate can be expressed as:

$$S\left({\widehat{A}}_j\right)=A_{tot}\times S\left({\widehat{p}}_{\bullet j}\right)$$

(5)

and an approximate 95% confidence interval for ${\widehat{A}}_j$ is:

$${\widehat{A}}_j\pm 2\times S\left({\widehat{A}}_j\right)$$

(6)

3. Results

Since CCA experiments analysed only the T₂ pixels belonging to a specific target class at T₁, whereas the NDVI technique compares the whole T₁ and T₂ image pixels, for comparison purposes, the shapefile of the specific T₁ target class considered by CCA was also overlaid on the input and output change images from NDVI experiments following the total image analysis and only the changes appearing on the delimited overlaid area have been investigated.

3.1. Italian Site

For the Italian site, all the patches of broadleaved, evergreen and mixed forest areas were merged into a single target forest layer.

Table 3. Change detection matrix for Murgia Alta site, CCA_HR and CCA_VHR results. Producer’s and overall accuracies based on stratified estimation. TH indicates the threshold value applied to the Z-statistic image in CCA experiments. A_m is the mapped changed area.

Change: Transition from Forest to No Forest—at Different TH for an Area of 1518 ha (at T1)
Experiment	Method	TH	Change User’s Acc.%	Change Producer’s Acc.%	No Change User’s Acc.%	No Change Producer’s Acc.%	Overall Acc.%	Am (ha) Mapped Area of Change	Stratified Changed Area Estimate with 95% Conf. Interv. (ha)
1	DIFF_NDVI_HR	DIFF > 0	83.58 ± 1.50	10.87 ± 0.57	49.88 ± 1.09	97.65 ± 0.50	52.18 ± 1.02	371.61	2856.19 ± 111.71
2.a	CCA_HR	CCA > μ + 1σ	72.96 ± 0.07	94.77 ± 0.01	46.49 ± 0.64	11.45 ± 0.05	71.11 ± 0.08	1412.13	1087.13 ± 2.48
		CCA > μ + 2σ	74.92 ± 0.08	75.33 ± 0.03	42.90 ± 0.23	42.35 ± 0.13	65.30 ± 0.09	1062.07	1056.19 ± 2.64
		CCA > μ + 3σ	82.27 ± 0.08	48.01 ± 0.05	41.12 ± 0.13	77.83 ± 0.13	57.50 ± 0.08	604.07	1035.24 ± 2.51
2.b	CCA_VHR	CCA > μ + 1σ	69.44 ± 1.05	92.28 ± 0.15	59.42 ± 5.95	21.78 ± 0.98	68.19 ± 1.18	1364.31	1026.71 ± 36.94
		CCA > μ + 2σ	70.51 ± 1.07	74.63 ± 0.41	52.51 ± 3.74	47.34 ± 1.78	64.47 ± 1.44	1036.44	979.30 ± 45.03
		CCA > μ + 3σ	72.97 ± 1.14	49.25 ± 0.73	45.66 ± 2.27	70.05 ± 2.09	57.12 ± 1.40	654.21	969.44 ± 43.69

presents the results obtained in the different experiments, described in Table 1. A window of the input image including forest patches and the corresponding change output image are shown in Figure 3, Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8 for the different experiments, respectively. More specifically, Figure 3a,b focus on a sub-window of the study area in T₁ and T₂ images, respectively. Figure 3c presents the output change map. Yellow, red and blue circles in the latter Figure evidence some change areas.

The changes detected were interpreted by visually analysing the close-up areas reported in Figure 4. Forest fragmentation by road (Figure 4a–c), tree density modification (Figure 4d–f) and burned forest areas (Figure 4g–i) seem evident in the protected area. The presence of burned areas observed in Figure 4h was confirmed by the State Forestry inventory and related to arson dated 2009.

When CCA was applied to Landsat imagery, CCA_HR, more changes could be detected than in previous experiment. The observable changes appear in Figure 5 and Figure 6.

Figure 7 shows the changed areas (Figure 7c) detected at VHR by CCA_VHR technique with threshold μ + 1σ (Experiment 2.b).

The close-up images of CCA_VHR in Figure 8 confirm the occurrence of the changes identified by CCA_HR (Figure 5) and show more details due to increased image spatial resolution. Table 3 reports the highest overall accuracy value.

3.2. Indian Site

For the Indian site, the CCA analysis focused on the strata of evergreen forest, as extracted from the existing LC/LU map. Quantitative results of the experiments are reported in

Table 4. Change detection matrix for the evergreen tropical forest in Indian site. Results obtained by CCA_HR and CCA_VHR. Producer’s and overall accuracy values are based on stratified estimation. TH indicates the threshold applied to the Z-statistic image in the CCA experiments. A_m is the mapped changed area.

Change: Transition from Evergreen Forest to No Evergreen Forest—at Different TH for an Area of 911 ha (at T1)
Experiment	Method	TH	Change User’s Acc.%	Change Producer’s Acc.%	No Change User’s Acc.%	No Change Producer’s Acc.%	Overall Acc.%	Am (ha) Mapped Area of Change	Stratified Changed Area Estimate with 95% Conf. Interv. (ha)
1	DIFF_NDVI_HR	DIFF > 0	44.29 ± 5.98	9.68 ± 1.88	65.26 ± 2.09	93.31 ± 0.62	63.33 ± 1.98	72.99	334.08 ± 37.31
2.b	CCA_HR	CCA > μ + 1σ	92.03 ± 2.31	34.62 ± 1.78	81.24 ± 1.84	98.95 ± 0.67	82.29 ± 1.67	93.24	247.86 ± 32.01
		CCA > μ + 2σ	98.70 ± 1.30	11.59 ± 0.76	73.54 ± 1.95	99.94 ± 0.16	74.40 ± 1.88	32.49	276.65 ± 36.02
		CCA > μ + 3σ	100.0 ± 0.01	4.66 ± 0.29	68.91 ± 1.98	100.0 ± 0.01	69.37 ± 1.95	14.31	307.36 ± 37.24
2.a	CCA_VHR	CCA > μ + 1σ	69.96 ± 0.26	43.96 ± 0.19	84.56 ± 0.16	94.20 ± 0.14	82.40 ± 0.14	134.53	214.10 ± 2.55
		CCA > μ + 2σ	85.35 ± 0.33	10.95 ± 0.13	72.60 ± 0.17	99.21 ± 0.04	73.09 ± 0.16	34.60	296.76 ± 2.93
		CCA > μ + 3σ	95.93 ± 0.30	2.54 ± 0.05	67.80 ± 0.17	99.95 ± 0.01	68.04 ± 0.17	7.90	298.53 ± 3.01

.

Figure 9 shows input images of experiment 1 (DIFF_NDVI_HR) and corresponding output image.

The DIFF_NDVI_HR experiment underestimated some changes. This is visible in close-ups of Figure 10 compared to the same close up windows shown in subsequent figures related to CCA applications.

In the Indian forest, more changes were detected by CCA_HR. As in the Italian site, the threshold μ + 1σ proved best in detecting most of the changes (Figure 11).

However, as the close-up windows of Figure 10 and Figure 12 may reveal, detecting and interpreting changes of HR (i.e., coarse) imagery, obtained by both techniques, may prove rather difficult, due to the complex features of dense tropical forest.

It seems worth observing that, when using the CCA_VHR technique, the largest OA was obtained with the threshold set at μ + 1σ. Accordingly, stratified changed area estimate with 95% confidence interval was the most accurate (error ±2.55). The inputs and the change images are shown in Figure 13.

The close-ups of changed areas, reported in Figure 14, show clear evidence of forest fragmentation due to the new road (Figure 14c). They also indicate forest degradation, which could be due to either community changes (Figure 14f–i) or other possible disturbances, such as invasion by alien plants. On the one hand, these findings confirm the need for VHR data related to conservation studies, and, on the other hand, they pose the need for important changes in the policies related to both VHR periodic image tasking over protected areas and relative costs.

4. Discussion

4.1. The Italian site

For the Italian site, the analysis of the quantitative results (reported in Table 3) indicate that:

• At comparable scale (HR), when CCA_HR is adopted, the overall accuracy is larger (68.19%) than the accuracy value obtained by DIFF_NDVI_HR (52.18%) with a significant reduction of the error in the stratified changed area estimate (from ±111.71 ha to ±36.9 ha). For example, the change area that is visible in Figure 5c (CCA_HR), close to the red and blue circles, was not detected in Figure 3c (DIFF_NDVI_HR). Moreover, the close-up in Figure 6c (CCA_HR) shows detailed changes due to new roads construction within and in the borders of the forest. These changes, which are clearly visible at VHR in Figure 7b (i.e., the Worldview2 image at T₂), are not detected by DIFF_NDVI_HR (Figure 4c), hence the underestimated changes in this case. In addition, CCA_HR analyses are more cost effective, since only the T₂ pixels that correspond to the specific target class in T₁ must be analysed_. The VHR image allows validation of the changes detected when in-field campaigns may not be feasible (Figure 7b).
• At fine scale, CCA_VHR provides the best overall accuracy (71.11%). This value is very close to the CCA_HR result (68.19%) but the former value is associated with the lowest error in the stratified changed area estimation (±2.48 ha), as already demonstrated in a previous study for grasslands ecosystems [21]. Both forest fragmentation by roads and changes in forest density are confirmed at finer scale by CCA_VHR (Figure 8c). Therefore, the best change results can be achieved by comparing fine scale LC/LU map with one single VHR image. These findings appear to support arguments in favour of the acquisition and storage of VHR images on protected areas at a regular base and, possibly, at a reduced cost for those institutions in charge of the protected areas management [8].

4.2. The Indian Site

For the Indian site, the results reported furnish evidence that:

• The changes identified in the evergreen forest by the DIFF_NDVI_HR are underestimated in comparison to the ones found by CCA_HR at the same scale. For instance, Figure 12c shows changes that are completely missed in Figure 10c. In addition, in the Indian experiment, the overall accuracy of CCA_HR is larger (82.29%) than the one from the DIFF_NDVI_HR (63.33%) with comparable errors in the stratified area estimate (Table 4). At HR, the results appear to confirm the inadequacy of direct NDVI image comparison for change detection in tropical forest (Figure 10i).
• For validation purposes, CCA was carried out also at VHR resolution, even if the scale of the existing LC/LU map is coarse (1:50,000). The output change images in the Figure 14c,f,i clearly validate the changes detected at HR by both CCA_HR. The changes observed may be mainly due to both deforestation and construction of new roads. Although the findings reported appear to be very interesting, more detailed information on the conditions of the changes in the tropical forest would be required. These may involve complementary in-field inspection and the collaboration of different scientific expertise.
• It must be recalled that, in the Indian site, the evergreen forest layer is adjacent to different vegetable classes named “woodland to savanna woodland (tall)” and “tree savanna”. Figure 10d,e evidence both such conditions and the presence of some changes in the vegetable classes which remained un-detected by CCA (Figure 10e). This discrepancy can be ascribed to the a priori class selection made to perform CCA analysis. In order to have inclusive detection of overall changes, additional CCA processing steps should be applied. However, the more inclusive process would be time consuming. On the other hand, DIFF_NDVI_HR can yield more comprehensive detection of possible changes in one single processing step, but the resulting overall accuracy would be lower than the one obtained through CCA. Undoubtedly, the choice of the most appropriate change detection technique will depend on specific user’s requirements as well as comprehensive cost effectiveness.

5. Conclusions

The study carried out offers grounds to suggest that CCA can be effectively applied for detection of changes in forest ecosystems at multiple scales. In fact, when an existing fine scale LC/LU map is available, the comparison of the map with a single recent image (T₂) can provide reliable results at both high and very high spatial resolution at rather low costs. The results obtained in the two study areas under investigation are comparable at both scales in terms of coverage (ha) of changed areas even though with reduced error at fine scale. However, for fine scale studies, the cost of VHR image acquisition is still very high for use by public bodies and decision makers. As clearly evidenced in [6,21], and in [3], agreements between space agency and national authorities should be encouraged to reduce such costs for a widespread EO data application to ecosystems monitoring.

CCA, applied to coarser images (e.g., Landsat 8 OLI) of both sites investigated, appears to provide better accuracy values than the ones obtained by a traditional technique (DIFF_NDVI_HR).

Estimates of uncertainty in the analysis were made on the basis of stratified sampling and recommendations found in [1,2]. However, it seems worth noting that the quality of the input LC/LU map and the semantic of the LC/LU classes can influence the accuracy of the change detection analysis.

The findings reported, prove that CCA technique can yield useful results for operational purposes, such as detection of change in specific classes (e.g., forest, and grassland) at different scale [30]. They also indicate that by using different threshold values, output change maps can be produced with different levels of change details and related accuracy values.

From the above discussion it can be concluded that CCA technique investigated could allow long-term monitoring of natural ecosystems in support to conservation management. This would require quantifying LC/LU change dynamics that can play an essential role in estimating forest accounts. The latter could provide a thorough measure of forest assets, flows of forest-related services and yield information on how these variables can change through time [31]. Forest accounts information are linked to traditional indicators such as the gross domestic products. They can also be extended to include other forest products such as fuel wood and ecosystem services.

The information produced can help design and monitor strategies for implementing the following United Nations (UN) Sustainable Development Goals (SDG) [32]: sustainable forest management (SDG 15), sustainable energy from fuel wood (SDG 7), and combat climate change and its impacts (SDG 13).