Mapping Cork Oak Mortality Using Multitemporal High-Resolution Satellite Imagery

Abstract

In the Mediterranean region, a significant decline in the vitality of vegetation has been observed in the last two decades, with a high forest mortality rate for several species. The increase in mortality has been attributed to water stress resulting from an increase in temperature and long periods of drought. To detect and quantify the impact of these events on tree mortality, an efficient and easy-to-use methodology for rapid damage assessment is required. Our study aims to assess the potential of high spatial resolution multispectral images from the Pleiades constellation to detect and map cork oak mortality in a pasture environment with multiple forest species. An approach based on change detection and the use of an unsupervised classifier is proposed to detect mortality at the cork oak level. The change in the values observed for three vegetation indices, NDVI, RGI, and GNDVI, between two epochs is used in an unsupervised classification algorithm to estimate the dead tree class. The classification results are accurate, with precision and recall values higher than 90%. Detailed cork oak mortality mapping is of significant use in comprehending ecosystem change as a result of tree mortality and for the implementation of mitigation mechanisms for the ongoing desertification process.

1. Introduction

Climate warming, which is associated with an increased frequency of drought and heat events, is positively related to forest mortality in several species worldwide [1,2,3,4]. Climate-induced tree mortality rapidly alters the ecosystem structure and function, having an impact on the ecology and biosphere–atmosphere dynamics [5,6]. According to the latest assessment report of the International Panel on Climate Change (IPCC) [7], over the last two decades in the Mediterranean region, the mean annual temperature has been increasing while precipitation has been decreasing, predicting an increased frequency of extreme drought and precipitation events [8,9]. The impact on water resources is causing groundwater levels to decrease, which consequently hinders cork oak growth since it depends on deep-rooting to access the water table [10]. In Portugal, a significant decline in cork oak vitality has been observed in the last two decades, with trees having a high mortality rate [11,12]. The increase in mortality has been attributed to water stress, of which the cork oak is subject due to the increase in temperature and the long periods of drought [13]. Besides land-use disturbances, which are mainly associated with grazing intensity/type and shrub encroachment, tree water stress is an additional cause of decline [14,15]. Cork oak is a special tree with a long life span (more than one hundred years), reaching its mature production of cork past 40 years. Cork oak was established as Portugal’s National Tree in 2011, mostly due to its economic, social, and environmental relevance. Besides the importance of the cork industry, it also promotes local population employment and contributes to the preservation of a wide range of flora and fauna habitats. Legislation protecting cork oak stands forbids the felling of trees, meaning that only dead or diseased trees can be cut down and only with permission from the National Nature Conservation and Forest Authority (ICNF).

The decline of the cork oak is a slow process that can take several years, wherein the first signs of vitality loss are expressed by a gradual loss of foliage [16]. The detection of foliage loss is not always possible and depends on several factors. For example, reduced foliage loss is undetected in the reflectance registered at the satellite level, or the spatial resolution of the sensor does not allow the individualization of the trees [17]. In the case of the cork oak, it should be noted that foliage may easily recover with the re-establishment of the hydrological cycle. Also, sometimes the foliage loss is only partial, with the remaining crown foliage intact. Detecting cork oak foliage decay is a complex process that can take several years. Multitemporal image analysis enables us to overcome some of the aforementioned limitations as long as some of the images are acquired in the tree’s healthy stage, just before the onset of vitality loss, so that the pre-death conditions of the tree can be quantified.

Given this, the goal is to identify, in the images time series, areas where a change in the tree’s vitality condition occurred, and to characterize the spectral signature of dying trees. As the cork oak is an evergreen tree, it performs photosynthesis for a longer period during the year when compared to other broadleaf trees that lose their leaves during winter. This makes the use of vegetation indices (VI) time series suitable for monitoring the phenology of the cork oak and identifying leaf phenological events. The Normalized Difference Vegetation Index (NDVI) has proved to be a good indicator of cork oak vitality, allowing the characterization of cork oak mortality in a time series of Sentinel-2 satellite images [18]. Other authors [17] have reported that the NDVI shows a good relationship with cork production at the stand level, revealing that one-tenth of the area of cork and holm oak stands in mainland Portugal has been significantly declining over the past 34 years. This trend, which was identified using an NDVI time series with images acquired during summer (July and August) by the Landsat and MODIS satellites, resulted not from short-term events, such as droughts or debarking, but reflects a long-term evolution of forest stand health.

In studies where greater detail is required for cork oak mortality analysis, for example at the tree level, a spatial resolution compatible with the size of the tree crown is required, i.e., less than 1 m [16,19]. Although freely available low–medium resolution satellite data have been successfully used to monitor cork oak decline on a national scale, high-resolution images are required to develop an operational system to be used by the national forestry associations to automatically detect unhealthy or dead trees on an annual basis, which implies detailed information at the tree level. This automatic detection largely decreases the efforts in the field by the forest owners for the identification of dead trees and quickens the permission request to ICNF for cutting down cork oak trees. Currently, several observation missions are available that fulfill these requirements of spatial and spectral resolution (Pleiades, QuickBird, WorldView). However, since the images acquired by these space missions are not freely available, their use is limited to a few images per year. High spatial resolution images have been used to detect changes in forests, such as burnt areas due to wildfires [20,21] and the occurrence of small patch-cutting (or isolated dead trees) in dense forests [22]. Clear-cuts in dense forests represent a good indicator of vitality loss, but may not be suitable for open forests, such as cork oak forests. In this case, studies will be necessary to assess the usefulness of high-resolution images to quantify cork oak mortality in areas of “montado” (cork oak in a pasture environment), where mortality is due to different causes such as pests, pathogenesis, or plants water stress. It is important to mention that low-medium resolution satellite data allow a national-scale analysis of cork oak mortality whereas high-resolution satellite data provide detailed information at the tree level.

This study aims to assess the potential of high spatial resolution multispectral images from the Pleiades constellation to detect and map cork oak mortality in a pasture environment with multiple forest species. An approach based on change detection and the use of an unsupervised classifier is proposed for the detection of mortality at the cork oak level. The rationale for the development of this approach was that it should be as simple as possible to be implemented as an operational system to be used by national forestry associations, and, most importantly, to be independent of in-situ data and to rely on few annual commercial high-resolution satellite data—due to its high cost. Three Pleiades images acquired in 2018, 2019, and 2020 in a cork oak forest area in central Portugal were used in this study. As a first step, the images are pan-sharpened, geolocalised, and radiometrically normalized. Then, the images are used to spectrally characterize dead and healthy trees to evaluate their spectral separability. The change in the values observed for three vegetation indices, the NDVI, the RGI (Simple Ratio Red/Green Index), and the GNDVI (Green Normalized Difference Vegetation Index), between 2018 and 2020 is used in an unsupervised classification algorithm to estimate the dead tree class. Good classification results were achieved with the developed approach, with both precision and recall exhibiting values higher than 90%. Detailed cork oak mortality mapping is of significant use in comprehending the ecosystem change as a result of tree mortality and for the implementation of mitigation mechanisms for the ongoing desertification process.

2. Materials and Methods

2.1. Study Area

The study was undertaken in one agro–silvo–pastoral system in central Portugal, namely Herdade da Machoqueira do Grou (HMG website) (Figure 1). Herdade da Machoqueira do Grou is located at the left bank tributary of the Tagus River, approximately 100 km from Lisbon. This region is characterized by a mild Mediterranean climate with Atlantic influences. The bioclimate is considered semi-arid to sub-humid, with a mean annual rainfall of 662.5 mm and a mean annual temperature of 16.3 °C.

Herdade da Machoqueira do Grou is a private agro-silvo-pastoral system with 2423 ha. Cork oak is the dominant species, occupying an area of 1017 ha and an additional 464 ha of a mixed-tree stand with stone pine. In the study area, the cork oak mortality is attributed to the lowering of groundwater levels, which consequently harms cork oak growth since it relies on deep-rooting to access the groundwater table. The entire watershed is on deep Miocene sands and the main soil types include fluvisols, leptosols, and podzols. The landscape is mostly plain, with altitudes ranging from 79 to 173 m, and slopes between 0% and 5%—with some exceptional examples being up to 35%. Cork oak’s average density is 90 trees/ha and the average production of dry weight cork is 1.3 ton/ha.

2.2. Pleiades Satellite Images

Pleiades 1A/1B are very high-resolution twin satellites offering 0.5 m spatial resolution products. The constellation is highly reactive and has been designed for daily revisit anywhere on Earth with stereo and tri-stereo capacity. The available output bands of the Pleiades block are: red (0.590–0.710 µm), green (0.500–0.620 µm), blue (0.430–0.550 µm), NIR (0.740–0.940 µm), and panchromatic (0.470–0.830 µm). The spatial resolution is 0.5 m for panchromatic and 2 m for multispectral bands. A single 100 km² orthoimage in a UTM/WGS84 projection covers the study area. For this study, three cloud-free images acquired on 9 October 2018, 21 October 2019, and 7 October 2020 were used. In early autumn, the humidity reaches its lower values, and, consequently, the trees’ activity is reduced until the winter dormancy period, thereby allowing for better discrimination between health and dead trees. Level-3 products provide orthorectified bottom-of-atmosphere reflectance in a JPEG2000 format. Multispectral images were pan-sharpened with the panchromatic band, resulting in three multispectral images with a 0.5 m spatial resolution [23]. The pan-sharpened transformation was performed within Orpheo Toolbox (© Copyright 2020 CNES) (Figure 2) and co-registered to the 2018 image using the SURF algorithm [24]. All the processed images have the same size: 14,301 columns, 10,869 lines, and 4 bands.

2.3. Radiometric Normalization

The ground reflectance recorded by Earth observation satellites is affected by the atmosphere during the radiation path between the Sun, the surface, and towards the satellite [25]. The correction of the atmospheric effects may be performed with atmospheric models [26] provided that the atmosphere properties at the time of image acquisition are known. When atmosphere properties are unknown, a simple radiometric normalization of the time series may be applied as an alternative [27]. Relative radiometric normalization generally assumes that the relationship between a pair of multispectral images, acquired at two different epochs, is assumed to be a linear function. The critical aspect lies in the identification of the time-invariant pixels to be used for the linear function parameter estimation [28]. In Canty & Nielsen [29], an automatic algorithm is proposed to determine the time-invariant pixels based on the invariance properties of the Multivariate Alteration Detection (MAD) transformation [30]. In our study, we applied an improved version of MAD called Iterative Reweighted MAD (IR-MAD) [31]. This new version of MAD performs better in isolating invariant pixels in images with large seasonal radiometric changes.

Radiometric normalization was performed using the 2018 reflectance values as a reference, being that the 2018–2019 and 2018–2020 image pairs were normalized independently. The green, red, and infrared bands were used. The normalization parameters of the 2020 invariant pixels are listed in

Table 1. Radiometric normalization parameters for the 2020 invariant pixels: interception, slope, and correlation. Values are digital numbers (DN).

Channel	Intercept	Std.err.	Slope	Std.err.	Correlation	RMSE
Red	−118.7	0.093	1.42	0.00035	0.99	9.60
Green	−241.5	0.254	1.57	0.00079	0.98	13.91
NIR	−51.5	0.200	1.08	0.00037	0.99	34.66

. In the 2020 image, there is a scaling effect in the visible bands, with values of 1.42 and 1.57, which is not observed in the infrared band. The reduction of this scale effect, as well as the bias (−118.7 and −241.5), independent of the band, will have a direct effect on the value of the vegetation indices. The radiometric normalization of the 2020 image mitigated the unmodeled atmospheric effects and the radiometric changes arising from the orthorectification and mosaicking process.

The Pleiades satellites, like most high spatial resolution satellites, are limited to four spectral bands—generally blue, green, red, and infrared. These four bands allow the calculation of multiple vegetation indices suitable for the detection of tree mortality. Some authors highlighted that the NDVI decay is an indicator of decreased photosynthesis, and therefore, is extremely useful in discriminating healthy from dying trees [32]. Other greenness-based indices are also mentioned as good indicators of leaf greenness level, being considered as differentiators of healthy and dead trees as well [33]. In this study, the indices NDVI [34], RGI [35], and GNDVI [36] were used to detect tree mortality.

The study area includes other land cover elements such as agricultural land, lakes, paths, and artificial areas that can be masked out from the images. As the objective is to detect trees in loss of vitality or dead between 2018 and 2020, only pixels within forest areas were considered for the analysis. A detailed land cover map, provided by the landowner (Figure 1, delimitation in yellow), was used to subset the forest areas, which include all the tree species present in the field (pine, holm oak, and cork oak). Pixels outside the forest area were masked out in the 2018 image.

2.4. Individual Tree Crowns Identification

Tree detection and crown delineation algorithms based on high resolution (0.5 to 1 m) multispectral images make use of the spectral reflectance properties of sparse forest areas (or “montados”), which are characterized by high canopy intensities in contrast to low surrounding ground intensities in all spectral bands. The approach consists of detecting the outer limits of the canopy and then the local maxima as a representative of each crown’s center. For this purpose, several methods have been proposed that combine image analysis methods with probabilistic methods [37] or morphological analysis [38,39] or by defining separability thresholds between trees and the ground (for a review, see [40]). The vast majority of studies use the visible and near-infrared bands and aim to identify trees in forests of moderate density [41]. With the increase in the spatial resolution of satellite sensors, and with the massive use of UAVs (Unmanned Aerial Vehicles), new algorithms for individual tree identification in open forests and orchards have emerged [18,38]. The identification of individual trees depends on several factors, of which the spatial resolution of the sensor or the distance between trees stands out—in particular when there is great variability in the size and shape of the trees, as is the case with cork oak stands. Cork oak stands are characterized by great diversity in the dimension and shape of the trees and by the irregularity of their distribution, often with the presence of tangled branches, which constrains the identification of individual trees.

The detection of trees in multispectral images is usually preceded by a spectral dimensionality reduction, transforming them into a grey-level image [42] or a vegetation index like the NDVI—or any of the abovementioned vegetation indexes. The vegetation indices have the advantage of implicitly promoting image segmentation into levels of photosynthetic activity related to the structure of the leaves. In the case of cork oak forests, it is expected that healthy trees with a higher leaf density will exhibit higher vegetation index intensities than the surrounding soil. However, in high-density cork oak forests, the reduced distance between trees (or groups of trees) results in the spectral mixing of soil and canopy, compromising the detection of the trees. In this study, the NDVI index was used. The NDVI is a vegetative index related to the leaf of the tree in the sense that it highlights the difference in spectral reflectance between the pigment and the leaf structure. The pigment, chlorophyll, strongly absorbs radiation in the red band while the leaf structure reflects more than 50% of the radiation in the near-infrared band. The use of NDVI seems to be suited to identify the contrast between trees and the surrounding soil, whether it corresponds to areas of little or no vegetation or areas with sparse herbaceous vegetation. However, in cork oak forests, the soil is often occupied by pasture or small shrubs with a significant photosynthetic activity, which decreases the soil/tree contrast.

The NDVI index image is then subjected to an object-based segmentation according to the local thresholding method of Sauvola & Pietikainen [43], in which a local separability threshold is defined for each pixel and calculated with the value of the mean and standard deviation of the vegetation index of the neighboring pixels. The choice of dimension is critical in this method. Given the spatial resolution of the images (0.5 m) and the average crown size (2 m), a 5 × 5 pixel window was adopted. To reduce the number of isolated pixels or small clusters of pixels, a low-pass median filter is applied to the binary image. The final result for Herdade da Machoqueira do Grou is presented in Figure 3b.

The filtered binary image is composed of pixel clusters that correspond to single trees, tree clusters, and noise. Here, noise is interpreted as clusters of pixels that do not correspond to trees. At this step, no method for noise detection/reduction was applied, except for the removal of pixel clusters with an area inferior to 8 m², assuming this value as the minimum area of the crown of an adult cork oak. Next, the binary image was vectorized, resulting in a set of polygons that are expected to correspond to trees or clusters of trees. The polygons resulting from the binarization of the image are divided into two groups according to their area and further processed distinctly. It was then assumed that polygons with an area larger than 20 m², which likely correspond to clusters of trees, and polygons related to individual trees should be dealt with differently. However, in both cases, the same strategy was used to identify the centers of the trees. For each polygon, the Euclidean distance map is calculated using the binarized image. The algorithm measures the distance of each pixel inside the polygons to its border and considers the pixel with the longest distance from the border as the polygon centroid. One can assume that the centroid corresponds to the tree crowns centre. Figure 3c shows the Euclidean distance map for the binary image. In this case, the Euclidean distance inside each polygon has a maximum value of 25 pixels, which is about 0.50 × 25 = 12.5 m.

The center of each tree is determined as the maximum distance inside the corresponding polygon. For polygons with multiple trees, the multiple tree centers are identified the same way, but on the condition that the local maxima are at a greater distance than the minimum distance between trees, i.e., 2 m. In this case, the process is executed recursively for all the local maxima of the distance function, starting with the highest value of the local maximum. Figure 3d shows the final results of the identification of the centers of the trees in Herdade Machoqueira do Grou. The final result has the tree crown centers’ locations.

A tree is assumed to belong to the class “dead tree” whenever there is a significant loss of photosynthetic activity, expressed by an abrupt reflectance reduction in the infrared band, or when a tree has been cut down. In the 2020 image, 897 dead trees were identified out of a total of 7572 trees automatically identified in 2018. The remaining 6675 trees were considered healthy. The identification of dead trees was performed by visual interpretation of the 2020 image in its false-color composite (IR, R, G) complemented with the visualization of Microsoft’s Bing Maps RGB images (dated 2021). The two classes of trees (healthy/dead trees) were used to characterize the spectral signature of each class and to validate the algorithm by calculating accuracy metrics. It is important to emphasize that the knowledge of the trees’ conditions is not incorporated into the image classification process once an unsupervised classifier is applied. This is another advantage of the developed approach, considering that the identification of dead trees is itself a source of error in the supervised classification. Although the vitality of the trees was consistently evaluated by a trained operator, the definition of a percentage of defoliation and discolored branches may be quite challenging and can be a source of error in the tree class assignment.

2.5. Tree Mortality Detection

The reduction in cork oak vitality is related to crown structure and leaf loss. This change can be detected in satellite images acquired at different stages of tree evolution [33,44]. The temporal difference of a VI image, relatively to a reference VI image, was used to characterize the temporal variation of the tree’s vitality [33]. The temporal difference of the indices was calculated as:

ΔVI (t) = VI (t = 0) − VI (t)

(1)

where ΔVI is the difference between a VI image calculated at a time (t) and its corresponding reference VI image (t = 0) so that the variation in the vegetation index is positive if the index value increases in time. In this study, the 2020 VI images are compared to the 2018 VI reference images.

The differences calculated for the three VI considered in this study are then combined into a 3-band image. The 3-band image is used in an unsupervised classification process that aimed to segment the image into classes. These algorithms aggregate the pixels into clusters with similar spectral characteristics according to a given model. In the case of k-means, pixels are aggregated according to their distance in the feature space from the center of the cluster. The feature space in this case is defined as the three-dimensional space of variation of the vegetation indices. The number of clusters is defined a priori. One of the estimated classes is expected to be the “dead trees” class, which also includes cut trees and trees with loss of vitality. The best model, defined by the combination of bands and the number of clusters, is determined by the analysis of the classification accuracy metrics.

2.6. Rating Evaluation Metrics

The performance of the classification model was evaluated by the computation of several accuracy metrics. For this purpose, confusion matrices (also referred to as double-entry contingency tables) were generated, where the rows correspond to the reference data (ground truth) and the columns refer to the predicted classes (retrieved by the classifier), to quantify the agreement between the automatic classification and the reality observed on the ground. In this study, the interest lies in the evaluation of the ability of the proposed algorithm to correctly estimate the “dead trees”. The most appropriate metrics to evaluate the classification accuracy in a binary classification (“dead trees” vs. “healthy trees”) are precision, recall, and F1-score.

Precision, or user’s accuracy, is the ability of the model to return relevant instances, i.e., the proportion of correctly classified plots relative to the total number of plots (correctly or incorrectly) classified in that class (column’s total):

$$\mathrm{precision}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}}$$

where TP are true positives, i.e., trees correctly classified as dead, and FP are trees that are classified as dead but are healthy in the reference. Precision is related to the commission error (=1 − precision).

Recall, or producer’s accuracy, refers to the completeness of the model, i.e., it is the proportion of correctly classified trees compared to the total amount of trees observed on the ground for a certain class. For each class, the number of correctly classified trees is divided by the total number of trees in that class (row’s total):

$$\mathrm{recall}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}$$

where TP are the true positives, i.e., trees correctly classified as dead, and FN are the trees that are classified as healthy but are dead trees. Recall is related to the omission error (=1 − recall).

The F1-score is a weighted average of precision and recall, with equal weight for both values (it also varies between 0 and 1). The expression to calculate the F1-score is:

$$\mathrm{F}1-\mathrm{score}=2\frac{\mathrm{precision}\ast \mathrm{recall}}{\mathrm{precision}+\mathrm{recall}}$$

In most cases, the question of a possible trade-off between precision and recall is to be considered. A generic model will maintain a balance between these two metrics that can be evaluated with the F1-score metric.

3. Results

3.1. Analysis of the Time Series of Vegetation Indices

A comparison of the vegetation index values for healthy trees and trees with vitality loss or those that are dead shows that healthy trees have higher NDVI and GNDVI values and lower RGI values (Figure 4a–c).

One can notice that mortality occurred after 2019 with a significant decrease in the NDVI and GNDVI values and a significant increase in the RGI values. The NDVI index demonstrates a higher ability to discriminate healthy from dead cork oaks trees when compared to the other two VI indices. According to the cumulative distribution function, for an NDVI value below 0.45, the probability of detecting a dead tree is 90% (Figure 4d). Only 10% of the dead trees have NDVI values higher than 0.5, which can be related to a misinterpretation of the trees’ conditions in the field or to the existence of shrub vegetation that has replaced the former tree. The results are consistent for the three dates under analysis, although there is a greater variability for the 2019 indices, thereby resulting in a lower ability to discriminate between dead and healthy trees. This is likely due to the presence of clouds in the 2019 image, in particular the existence of thin, non-perceptible clouds over the study area. The GNDVI index behavior is very similar to the NDVI index, which is justified by the strong relationship between the vegetation vitality and the infrared reflectance. There is, however, a lower ability to discriminate between healthy and dead trees. For a GNDVI value below 0.40, the probability of occurrence of dead cork oak is 82% (Figure 4f). The RGI index has a lower ability to discriminate between dead and healthy trees, with a confusion region between the two classes that equals an error of at least 15% (probability of finding a dead tree with an RGI value below 0.9).

It is expected that radiometric differences in the image remain due to atmospheric variability and lighting conditions (shadows) as well as changes in land cover resulting in reflectance variations.

For each pair of indices, the relationship between the vegetation indices differences between the reference image (2018) and the 2020 image was analyzed (Figure 5). For the ΔNDVI-ΔRGI pair (Figure 5a), a negative relationship was observed in which pixels with an ΔNDVI positive variation have an ΔRGI negative variation, and vice versa. The two classes are logically distributed along the correlation axis, with dead trees concentrated in the upper left quadrant, while healthy trees are mainly concentrated in the lower left quadrant yet spreading to the center of the graph. The highest density is concentrated at the origin of the frame, that is, around zero. The level of separability between the two classes is defined by the confidence ellipse at a 90% level, which is represented in black in Figure 5a. The confidence ellipse defines the region that contains 90% of the data for each class, assuming that the variation of the vegetation indices has a normal distribution with a mean and a standard deviation. The correlation coefficient is identical for dead trees and healthy trees, −0.59 and −0.54, respectively. Curiously, the two classes have a distinct distribution pattern in the feature space. Unhealthy and dead trees have a larger dispersion than healthy trees, as shown by the area of the confidence ellipse, which is twice as large for dead trees. Unhealthy trees start to lose leaves, which leads to an increase in the spectral response from the ground surrounding the trees, resulting in omission errors that can only be minimized by increasing the commission error. From the analysis of the graph, it is clear that the two classes are disjointed in the ΔNDVI versus ΔRGI indices feature space.

Contrarily, for the ΔNDVI-ΔGNDVI pair, a positive relationship is observed between dead and healthy trees, with a correlation coefficient of 0.85 and 0.81, respectively. As in the previous pair, the two classes are distributed along the correlation axis, with healthy trees occupying the positive quadrant and dead trees occupying the negative quadrant. Likewise, a high level of separability is observed between the two classes for a 90% confidence level, which is represented in black in Figure 5b.

3.2. Detection of Tree Mortality

The former three vegetation indices proved to be suitable for the detection of dead trees and were then used as inputs for an unsupervised classification process. To decide on the optimal vegetation indices combination for cork oak mortality detection, three combinations were evaluated: (ΔNDVI, ΔRGI, ΔGNDVI), (ΔNDVI, ΔRGI), and (ΔNDVI, ΔGNDVI). For each combination, tests were carried out with different numbers of clusters: 3, 4, 5, and 6. The tree mortality detection performance was assessed using the following accuracy metrics: precision, recall, and F1-score, assuming that true positives were the correct estimation of dead trees.

For all three datasets, an increase in the number of clusters is positively related to precision, while the opposite is verified for recall (Figure 6). F1-scores do not seem to be related to the number of clusters. An increase in the number of clusters results in the segmentation of the image into more homogeneous regions, which is expected to minimize the area of each class in the feature space. Thus, increasing the number of clusters will allow the reduction of the cluster area corresponding to dead trees. Consequently, a smaller number of trees will be classified as dead trees, decreasing the number of false positives and therefore increasing accuracy. For the same reason, a greater number of false negatives is observed, and therefore there is a decrease in recall.

The behavior of the midpoint (center) of clusters corresponding to the dead tree class was analyzed as a function of the increase in the number of clusters. We verified that, with the increase in the number of clusters, the dead tree cluster centers shift away from the origin of the feature space (Figure 7), leading to an increase in precision, and, consequently, to a decrease in the number of healthy trees misclassified as dead trees. For the same reason, dead trees located near the separability threshold between the two classes are likely to be erroneously classified as healthy trees.

In this study, the main objective is to develop a methodology for dead tree detection from the producer’s perspective, i.e., maximizing the estimated number of dead trees while keeping a low commission error. The appropriate metric to evaluate the accuracy of the estimate against the ground truth is the recall, or, equivalently, the omission error. According to the results in Figure 6, an increase in the recall is achieved by reducing the number of clusters as long as the accuracy is not significantly reduced. The F1-score metric, which combines both the precision and recall, was used to define the optimal number of clusters. For all indices’ combinations, the F1-score was maximal for four clusters, so that was the number of clusters adopted for the k-means classification.

Accordingly, the class “dead tree” is associated with the cluster centered in the upper left quadrant, and, in case there is more than one, with the cluster with a center in the upper left quadrant farthest from the origin of the feature space. The remaining trees are associated with the class “healthy trees”.

The result of the mortality classification using the three indices (ΔNDVI, ΔRGI, ΔGNDVI) and the k-means algorithm with four clusters is presented in Figure 8. Out of a total of 7572 trees identified in the study area, 882 trees were classified as dead. The accuracy metrics were calculated from the confusion matrix data in

Table 2. Confusion matrix for the result of the k-means classifier using the three vegetation indices. Columns represent the classification results and rows represent the reference data.

		Classification
		Dead	Healthy	Recall
Reference	Dead	825	72	0.92
	Healthy	57	6618	0.99
	Precision	0.94	0.99
	F1-score	0.93	0.99	OA = 0.98
	Macro F1-score	0.96

. According to the data in Table 2, dead trees were identified with a precision of 94%, a recall of 92%, and an F1-score of 93%, meaning that only 6% of trees were incorrectly classified as dead trees (commission error) and only 8% of dead trees were omitted from the classification (omission error).

3.3. Spatial Pattern of Mortality

Mortality has an irregular spatial distribution pattern characterized by high concentration areas of dead trees, with 21 dead trees per ha, and scattered dead trees throughout the study area, with a density of about 2.2 dead trees per ha (Figure 8). The spatial distribution of dead trees is related to forest management activities and loss of tree vitality.

The highest density of dead trees is observed to the west of the property in an area of pine trees (Figure 8). Other smaller concentrations of dead trees are also observed on the perimeter of the study area, also corresponding to areas of pine trees. The comparative analysis of areas with high dead trees density in the 2018 and 2020 images revealed that those are areas where logging took place. Figure 9a,b show, respectively, the 2018 and 2020 images in a false-color composite (NIR, R, G) for the westernmost area of pine trees in Figure 8. Comparing the two images, a reduced number of trees are notorious in the 2020 image due to tree felling. Figure 9c presents the results of the k-means classification for the class “dead tree” in light blue overlaid on the 2020 image. The visual analysis of the classification results demonstrates the good performance of the classification approach in detecting land cover changes. It is perceptible that the detection of changes in pixels adjacent to the crown of trees results from differences in the geometry of image acquisition due to—albeit slight—differences in the tree crown’s position. One can also observe the detection of changes resulting from variations in vegetation cover as a consequence of annual rainfall variability.

Dispersed mortality is characterized by trees with a significant change in the spectral response resulting from a reduction in the tree’s photosynthetic activity. Cork oak is predominant in this group. Figure 9d,e show images for 2018 and 2020, respectively, in a false-color composite (NIR, R, G) for a cork oak trees area. In Figure 9f the results of the k-means classification for the class “dead tree” is presented in light blue overlaid on the 2020 image. The cork oak trees have an irregular distribution pattern, with great spatial variability in density, and can vary between isolated trees to clusters of trees. This spatial variability, associated with its small leaf size, may constrain the characterization of the spectral response of healthy (or dead) trees and the subsequent detection of spectral changes over time. Figure 9d shows three trees where a decrease in the spectral response was observed between 2018 and 2020, which were correctly assigned to the class ‘dead trees’ by the classifier assigned (see Figure 9f). The trees, still visible in the 2020 image, exhibit a significant reduction in the NIR band reflectance and an increase in the red band reflectance, and therefore reveal a brownish color. Also, in Figure 9f, a fair agreement between the spectral response changes detected by the classifier in light blue and the identification of dead trees is evident.

4. Discussion

We aimed to evaluate the potential of high spatial resolution multispectral images for the detection of cork oak trees for loss of vitality or trees that were dead. For this purpose, we used the temporal variation of three vegetation indicators (NDVI, RGI, and GNDVI) calculated from two Pleiades images acquired in 2018 and 2020 to classify areas with the presence of dead trees or trees with decreased foliar density, and also areas where trees have been cut down. The results show that the proposed methodology accurately maps the trees’ vitality, enabling the identification of the dead trees’ location and their quantification, and, consequently, making the request for permission from the national forestry entity to cut down dead cork oak trees easier and promoting more efficient forest management. However, in the proposed approach, some issues are essential to achieve high performance for tree mortality detection.

First, the identification of the trees’ location is of extreme relevance. The knowledge of the trees’ positions limits the analysis of the vegetation indices variability to the corresponding canopies. This is critical in multispectral studies where herbaceous vegetation and soil can introduce noise into the process of evaluating the spectral separability between dead and healthy trees [45]. The process of tree identification can become difficult or even impossible in areas of low contrast between the tree crowns and the surrounding vegetation [12,46]. To increase the contrast between the tree crowns and the surrounding herbaceous vegetation, one should use the images acquired at the end of the dry season, before the first rains of the season (in the Mediterranean, between the second half of September and the first half of October). At this time of the year, the shrub vegetation is dry, with reduced photosynthetic activity, and therefore, exhibiting a high red radiation reflectance in contrast to a high red radiation absorption for the healthy trees [18].

In the case of cork oak forests, where the density is generally low and the tree crowns are individualized, the adaptive binarization algorithm proved to be efficient in identifying the canopy of isolated trees [43]. However, in areas of dense forest or where groups of cork oaks canopies occur, the algorithm fails to identify the individual trees. A new algorithm was proposed in which the center of each crown is identified as a local maximum of the Euclidean distance of the binary mask resulting from the adaptive binarization. The segmentation of tree crowns clusters into a set of separate trees enables the detection and analysis of vitality at the tree level.

The second crucial issue is the radiometric calibration of the images. In multi-temporal studies, the meteorological conditions at the time of image acquisition differ and may interfere with the radiance recorded at the sensor level. The MAD algorithm identifies a set of pixels based on which a linear relation is established for each band [29]. The radiometric calibration of the image has a direct effect on the intensity of the spectral bands and an indirect effect on the vegetation indices calculation. Vegetation indices are band ratios or ratios of algebraic band operations in which the effect of global linear calibration is mitigated or even eliminated. The calibration of individual bands has the effect of changing the vegetation indices, as was verified for the 2019 image calibration (see Figure 4). It should be noted that small variations in the vegetation indices can disguise the trees’ vitality trend in the classification process. Therefore, the images’ radiometric calibration and their co-registration are critical in the process of change detection and subsequent mortality classification.

Another significant issue is the selection of the vegetation indices. Although band differences might be used for change detection, it is recognized that vegetation indices are more suitable for vegetation vitality detection/monitoring [33,47]. The four spectral bands in the Pleiades images: blue, green, red, and infrared, which are characteristic of all high spatial resolution Earth observation satellites and also of the UAVs, were used to compute the NDVI, RGI, and GNDVI. These three indices proved to be good indicators of tree mortality (Figure 5). The NDVI and GNDVI indices are the ones that showed a higher ability to discriminate between dead and healthy trees with a 90% probability for the detection of dead trees for a given threshold [33]. For 2020 the separability is around 10%. The relationship between the vegetation indices behavior and the evolution of the evergreen trees’ vitality is well-founded theoretically. The chlorophyll present in the leaves absorbs essentially the radiation in the red band, to carry out photosynthesis, and strongly reflects the radiation in the infrared band (60%) and to a lesser extent in the green band (about 20%) [48]. When the tree loses vitality, the reduced photosynthetic activity causes an increase in the red band reflectance and a decrease in the infrared band reflectance. Therefore, the RGI index increases, and the NDVI and GNDVI decrease, for trees in loss of vitality relatively to healthy trees. During the process of cork oak vitality loss, trees lose their leaves leading to a reduction of the canopy area and consequently to a decrease in the latter two vegetation indices values. It is important to emphasize that cork oak vitality loss is a slow process that can last for at least 5 up to 10 years [17]. Small variations in vegetation indices may indicate the beginning of cork oak vitality loss, which should be confirmed for subsequent years. According to the previous rationale, it is expected that in case of consistent trees’ vitality loss over the years, the RGI values will increase and the NDVI and GNDI values will decrease. There is another important aspect. As tree crown area decreases due to leaf loss, the ground reflectance over time becomes dominant generating misclassifications. The use of the GNDVI index is critical in distinguishing green leaves from shrub vegetation. Also critical is the choice of the images’ date, which should be acquired in the period between the end of the dry season and the beginning of the rainy season to minimize the influence of the surrounding vegetation [49].

Finally, the fourth and last issue is the choice of the classifier. The k-means algorithm is an intuitive unsupervised classifier that is easy to interpret and analyze and is available on most computer platforms. The reason for choosing an unsupervised classification algorithm is that it does not require data for training the classifier [50]. This aspect is important for the operationalization of the classification process and its use by technicians not familiar with supervised classification processes. Supervised classification processes have been widely and successfully used in the classification of tree mortality [51]. However, it requires the acquisition of training data, which is always a time-consuming and expensive process and not always susceptible to spatial and temporal transferability. That is, the model trained for a given location and date may not apply to another location or other dates. It is essential to underline that the proposed methodology does not rely specifically on the k-means algorithm, but on any unsupervised classification algorithm. Other unsupervised classification algorithms such as the maximum expectation [44,52] or Isodata [53] might have been used, though some preliminary evaluations did not recognize any advantages over the k-means algorithm for detecting cork oak mortality. Although no experiments were performed, we strongly assume that the proposed methodology can be applied to very high-resolution images acquired by UAVs [54]. In this case, the issues regarding different acquisition geometries, different illumination conditions, and co-registration are of particular relevance.

The methodology presented here can be automated for the annual monitoring of cork oak forests with high spatial resolution satellite images or UAV images, either as a sole monitoring tool or as a support to other field data-collecting instruments, or even to other techniques using lower resolution images.

5. Conclusions

In this study, we present a methodology for cork oak tree mortality detection using multitemporal high-resolution Pleiades images. The methodology combines unsupervised classification with temporal differences for three vegetation indices, NDVI, RGI, and GNDVI, for the detection of dead or unhealthy trees. The analysis of the vegetation indices showed that healthy and unhealthy trees can be discriminated with a 90% probability using the NDVI vegetation index. Also, the behavior expressed by the NDVI and the GNDVI is similar, while the RGI has an opposite behavior. Regarding the bi-temporal differences in the vegetation indices, the results showed a negative correlation between the ΔNDVI and the ΔRGI, whereas a positive correlation was observed between the ΔNDVI and the ΔGNDVI. The correlation coefficients are −0.54 and 0.85, respectively. The three vegetation indices were considered relevant for the detection of dead trees and therefore used in unsupervised classification. The experiments pointed out that an increase in the number of clusters in the unsupervised classification has the effect of increasing the precision (lower commission errors), at the expense of decreasing the recall (higher omission errors). The F1-score does not appear to be related to the number of clusters. The best classification performance was achieved from the combination of the three bi-temporal vegetation indices using 4 clusters, with 92% and 94% for recall and precision, respectively, and an F1-score of 93%. The proposed methodology proved to be efficient in areas that have experienced a natural cork oak vitality decline, including areas where trees have been cut down. This is of utmost importance, especially for dead cork oak tree detection in mixed-species woodlands, given that it provides forestry owners/managers an operational monitoring tool for the quick detection of cork oak vitality loss/tree mortality and allows for a timely forest management adaptation.