Tree Species Classiﬁcation and Health Status Assessment for a Mixed Broadleaf-Conifer Forest with UAS Multispectral Imaging

: Automatic discrimination of tree species and identiﬁcation of physiological stress imposed on forest trees by biotic factors from unmanned aerial systems (UAS) o ﬀ ers substantial advantages in forest management practices. In this study, we aimed to develop a novel workﬂow for facilitating tree species classiﬁcation and the detection of healthy, unhealthy, and dead trees caused by bark beetle infestation using ultra-high resolution 5-band UAS bi-temporal aerial imagery in the Czech Republic. The study is divided into two steps. We initially classiﬁed the tree type, either as broadleaf or conifer, and we then classiﬁed trees according to the tree type and health status, and subgroups were created to further classify trees (detailed classiﬁcation). Photogrammetric processed datasets achieved by the use of structure-from-motion (SfM) imaging technique, where resulting digital terrain models (DTMs), digital surface models (DSMs), and orthophotos with a resolution of 0.05 m were utilized as input for canopy spectral analysis, as well as texture analysis (TA). For the spectral analysis, nine vegetation indices (VIs) were applied to evaluate the amount of vegetation cover change of canopy surface between the two seasons, spring and summer of 2019. Moreover, 13 TA variables, including Mean, Variance, Entropy, Contrast, Heterogeneity, Homogeneity, Angular Second Moment, Correlation, Gray-level Di ﬀ erence Vector (GLDV) Angular Second Moment, GLDV Entropy, GLDV Mean, GLDV Contrast, and Inverse Di ﬀ erence, were estimated for the extraction of canopy surface texture. Further, we used the support vector machine (SVM) algorithm to conduct a detailed classiﬁcation of tree species and health status. Our results highlighted the e ﬃ ciency of the proposed method for tree species classiﬁcation with an overall accuracy (OA) of 81.18% (Kappa: 0.70) and health status assessment with an OA of 84.71% (Kappa: 0.66). While SVM proved to be a good classiﬁer, the results also showed that a combination of VI and TA layers increased the OA by 4.24%, providing a new dimension of information derived from UAS platforms. These methods could be used to quickly evaluate large areas that have been impacted by biological disturbance agents for mapping and detection, tree inventory, and evaluating habitat conditions at relatively low costs.


Introduction
Climate-related outbreaks of bark beetle species pose a serious threat to the temperate forests of Europe. Rapidly increasing forest disturbances, such as moths [1], pine sawflies [2], and bark beetles [3], gave rise to substantial economic loss and other values [4]. Expected increases in the earth's temperature will undoubtedly give rise to more frequent storms and increased severity throughout the globe [5], implying increasingly suitable conditions (e.g., increased production of insect-breeding material, They achieved an OA of 90% for classification of two groups (healthy, dead) while the classification of field data to more detailed classes (healthy, infested, dead) decreased the OA to 76%. Even though their results are satisfying, showing the great potential of hyperspectral imagery, yet the acquisition cost of hyperspectral sensors is extremely high. In another study, Klouček et al. [28] also used UAS-based multispectral imagery for the detection of Ips typographus L. in spruce trees in Krkonoše National Park in the Czech Republic. Their study produced an OA of 78%-96% according to the greenness index.
In addition to the benefits of versatility, wide availability, and ergonomics that UASs offer, they could also be employed for the evaluation of plant biological conditions, for example, leaf chlorophyll and nitrogen (N) content, and understanding plant forest health status. Nitrogen is a basic element in enzymes and chlorophyll, thus shortage of N in leaf canopies results in inadequate photosynthetic rates. Unhealthy tree canopies reflect more in the red band than healthy tree canopies. Consequently, the spectra of the problematic tree canopies will significantly differ from the normal canopy reflectance spectra. It has been shown in several studies that there is a high correlation between the chlorophyll content and the red-edge band [29,30]. Multiple vegetation indices (VIs) have been proposed for estimating the amount of N-chlorophyll content in forest tree canopies [31]. Merzlyak and Gitelson [32] proposed the plant senescence reflectance index (PSRI) to determine the stage of leaf senescence, as it is susceptible to carotenoid retention. PSRI values increase during leaf senescence due to the increase in the ratio of carotenoids to chlorophyll. In Northern Europe, a few studies conducted during the growing season showed that PSRI in autumn had similar performance in terms of sensitivity to vegetation dynamics with NDVI [33,34]. Horler et al. [35] was the first who demonstrated the significance of the red-edge band for plant stress detection. Gitelson and Merzlyak [36] exhibited that in the red-edge region of the electromagnetic spectrum, chlorophyll has lower absorption rates and a lesser saturation effect, thus the reflectance remains sensitive to chlorophyll absorption. Wu et al. [37] proposed the replacement of NIR and red bands in the transformed chlorophyll absorption reflectance index/optimized soil adjusted vegetation index (TCARI/OSAVI) [38] and modified chlorophyll absorption in reflectance index (MCARI)/OSAVI [39] by bands respectively at 705 nm and 750 nm. Another study conducted by Clarke et al. [40] revealed the robustness of the normalized difference red-edge (NDRE) in estimating N-chlorophyll content. For estimating chlorophyll content Gitelson et al. [41] proposed the red-edge simple ratio index (RESR). One year later Gitelson et al. [42] suggested the chlorophyll index red-edge (CI). More recently, Perry et al. [43] proposed the modified canopy chlorophyll content index (M3CI) for estimating N content in Pyrus communis L.
From the literature survey, most existing studies on tree species mapping and forest health assessment from UAS-based multispectral images have mainly focused on rural and urban forests rather than managed forests [12,25]. Most are solely based on RGB sensors [14,26] and limited to just spectral analysis [28]. Studies on plant nurseries, with relatively low tree density and low levels of structural complexity [24] and studies conducted on low vegetation [26] are insightful, but forest conditions may pose unique challenges for detection. With respect to the aforementioned studies and their limitations, we propose a more diverse and detailed methodological framework, using a combination of spectral and textural information and support vector machines (SVM) techniques, followed by several statistical techniques for tree species classification and mapping of unhealthy (infested) and dead trees attacked by bark beetles at the tree-level.
SVM is a non-parametric pattern classifier algorithm suitable for both regression and pixel-level classification techniques based on statistical learning theory [44]. The prerequisite for SVM to achieve better results is the appropriate determination of the parameters that play key roles in achieving higher accuracy and better performance [45]. The specified grid search using V-fold cross-validation [46] is the most commonly used method to identify suitable parameters, i.e., epsilon (ε) and capacity (C) with fixed gamma that would produce high-accuracy results. The current literature generally recognizes the SVM for its ability to perform better with limited training samples [47], since it utilizes only the subset of the training samples that define the location of the SVM optimum hyperplane, which is necessary to distinguish between classes [48]. SVM is preferred in the case of several classes and with imbalanced training datasets, thus it has been widely used in several tree species classification studies [49][50][51][52][53][54].
Given the fact that biotic agents have a continuous negative impact on forest resources, affecting several species and forests globally, this study provides new information about the applicability of novel remote sensing methodological approaches, through complex statistical analysis for monitoring-mapping tree health status. Aiming to support precision forest management and decision-making for small-scale areas. In this context, the main objectives of this study were to test how bi-temporal aerial multispectral (5-band) imagery data can be utilized for: (1) testing the sensitivity of specific spectral indices for classifying tree species, and (2) health status assessment by assessing and discriminating crown symptoms.

Study Area
The forest site used in this study is located near the village of Oplany within the territory of the Training School Forest Enterprise in Kostelec, southeast of Prague, Czech Republic ( Figure 1). The study area extends geographically from 49 •  We established an experimental area of 180 × 160 m with an altitude of~430 m a.s.l. and a mean annual temperature of 7.5 • C and mean annual precipitation of 600 mm. The study area is dominated by even-aged (originating from a clear-cut silvicultural treatment) coniferous tree species with an admixture of several broadleaf species, with a small percentage of understory (~20%). Norway spruce (Picea abies L.) and Scots pine (Pinus sylvestris L.) were the dominant conifer species, comprising 70% of all trees, accompanied mainly by Downy oak (Quercus pubescens Willd.), Hornbeam (Carpinus betulus L.), and Silver birch (Betula pendula Roth.), which comprised the remaining 30% of trees on the site. Norway spruce was the dominant species (75% of all conifers), and the three broadleaved species were equally represented. We eliminated all understory trees with a diameter at breast height (1.3 m; DBH) smaller than 15 cm from the field data. Field observations determined that there was evidence of increasing bark beetle outbreaks at alarming rates at the research site.
Remote Sens. 2020, 12, x 4 of 22 Given the fact that biotic agents have a continuous negative impact on forest resources, affecting several species and forests globally, this study provides new information about the applicability of novel remote sensing methodological approaches, through complex statistical analysis for monitoring-mapping tree health status. Aiming to support precision forest management and decision-making for small-scale areas. In this context, the main objectives of this study were to test how bi-temporal aerial multispectral (5-band) imagery data can be utilized for: (1) testing the sensitivity of specific spectral indices for classifying tree species, and (2) health status assessment by assessing and discriminating crown symptoms.

Study Area
The forest site used in this study is located near the village of Oplany within the territory of the Training School Forest Enterprise in Kostelec, southeast of Prague, Czech Republic ( Figure 1). The study area extends geographically from 49°55'6.25" N; 14°50'34.84" E to 49°54'56.34" N, 14°50'54.26" E. We established an experimental area of 180 × 160 m with an altitude of ~430 m a.s.l. and a mean annual temperature of 7.5 °C and mean annual precipitation of 600 mm. The study area is dominated by even-aged (originating from a clear-cut silvicultural treatment) coniferous tree species with an admixture of several broadleaf species, with a small percentage of understory (~20%). Norway spruce (Picea abies L.) and Scots pine (Pinus sylvestris L.) were the dominant conifer species, comprising 70% of all trees, accompanied mainly by Downy oak (Quercus pubescens Willd.), Hornbeam (Carpinus betulus L.), and Silver birch (Betula pendula Roth.), which comprised the remaining 30% of trees on the site. Norway spruce was the dominant species (75% of all conifers), and the three broadleaved species were equally represented. We eliminated all understory trees with a diameter at breast height (1.3 m; DBH) smaller than 15 cm from the field data. Field observations determined that there was evidence of increasing bark beetle outbreaks at alarming rates at the research site.

Reference Data
Field inventory and visual identification of tree species and health status assessments were conducted in August 2019. We assessed discoloration as a typical crown symptom, indicating the level of infestation which was caused by bark beetle infestation. Discoloration refers to any color anomalies occurring in live leaves [55]. Field-based species recognition and health status assessments (including bark visual analyses) were conducted by experienced forest experts in the study area by marking all the trees. Trimble M3 total station (Trimble Inc. Sunnyvale, CA, USA, 1978) was used to measure the tree positions. Initially, for evaluating the relation between reference data and predictor variables, all trees were classified according to tree type and health status into four classes: (i) healthy broadleaves = B, (ii) healthy conifers = C, (iii) infested conifers = I, and (iv) dead conifers = D (general classification). Afterward, depending on the health status and tree type, subgroups were created and the best independent variables were used for predicting the tree species and health status (detailed classification), using the SVM algorithm (

UAS Survey Planning and Data Collection
In total, two flights, conducted on 17 April and 24 July of 2019, were performed. Data from both flights were acquired between 11:00 a.m. and 14:00 p.m. local time under favorable weather conditions for photographing, were conducted using a rotary-wing DJI S900 (Dá-Jiāng Innovations Science and Technology Co. Ltd., Shenzhen, China) UAS equipped with a MicaSense RedEdge-M multispectral camera (MicaSense, Seattle, WA, USA); sensitive to the RGB, red-edge, and NIR spectral regions ( Table 2). Since MicaSense RedEdge-M features global shutters it is less affected by normal vibration, gimbal was not a requirement in our study. Although we did not apply any interior calibration for the camera, the MicaSense calibrated reflectance panel (CRP) was used to convert raw pixel values into reflectance, using an image of CRP taken before and after the flight [56]. This was a necessary step during the processing phase in Metashape professional edition 64-bit V. 1.5.3.8469 (Agisoft LCC, St. Petersburg, Russia). The above mentioned procedure allows image pixel values to accurately describe the material composing the object of interest by compensating lighting and atmospheric conditions. Without the calibration process, data captured over different dates cannot be accurately compared for change detection. Furthermore, the copter was guided semi-automatically, with the route planned for a height of 110 m at a speed of 2 m·s −1 and capture rate of one image per second in all bands simultaneously, recording the position of each image using a high-performance GNSS receiver integrated into the DJI A2 flight control system of the DJI S900. According to the capture rate of the multispectral camera (1 image·s −1 ), flight planning resulted in 85% frontal overlapping and 70% side overlapping. DJI GO Version V. 4.1.9 was the program used for controlling the DJI S900. After approximately eight minutes of flight (based on the pre-defined flight parameters), manual control was regained and the copter was safely landed. For both flights, a single serpentine pattern was used with a nadir-pointing camera.

Photogrammetric Processing and Model Extraction
Multispectral aerial images from both periods were photogrammetrically processed using SfM in Metashape. To ensure optimal photo alignment results during the photogrammetric processing, the automatic image quality estimation feature in Metashape was considered to identify poorly focused images with a quality value of less than 0.5 units. The total number of images included in the alignment process for the first dataset (spring 2019) was 153 out of 153, while for the second dataset (summer 2019) 144 of 144 images taken were used. Feature extraction and matching were performed using the SIFT algorithm [18] to detect image feature points and match them between images. A bundle block adjustment was applied on the matched features between the images, to identify the XYZ location and camera orientation of each image feature resulted in a sparse 3D point cloud.
For photo alignment and dense point cloud generation, we chose ultra-high quality, resulting in little noise during the densification process, while for the depth filtering, the mild option was used. The ground points were then classified and mesh models were created from all point classes. The classification was automatic and it was based on three parameters: (a) maximum distance (m), (b) maximum angle (deg.) and (c) cell size. Texture was applied using the orthophotos for the mapping mode option. Digital terrain models (DTMs) were generated from the ground points and digital surface models (DSMs) were produced from the points of all classes for both datasets, with a resolution of 0.05 m. Additionally, orthophotos were exported for both datasets (Figures 2 and 3). Before any process, we generated normalized canopy height models (nCHMs) by subtracting the DSM from the DTM in ArcMap 10.3.1 by ESRI © .

Photogrammetric Processing and Model Extraction
Multispectral aerial images from both periods were photogrammetrically processed using SfM in Metashape. To ensure optimal photo alignment results during the photogrammetric processing, the automatic image quality estimation feature in Metashape was considered to identify poorly focused images with a quality value of less than 0.5 units. The total number of images included in the alignment process for the first dataset (spring 2019) was 153 out of 153, while for the second dataset (summer 2019) 144 of 144 images taken were used. Feature extraction and matching were performed using the SIFT algorithm [18] to detect image feature points and match them between images. A bundle block adjustment was applied on the matched features between the images, to identify the XYZ location and camera orientation of each image feature resulted in a sparse 3D point cloud.
For photo alignment and dense point cloud generation, we chose ultra-high quality, resulting in little noise during the densification process, while for the depth filtering, the mild option was used. The ground points were then classified and mesh models were created from all point classes. The classification was automatic and it was based on three parameters: (a) maximum distance (m), (b) maximum angle (deg.) and (c) cell size. Texture was applied using the orthophotos for the mapping mode option. Digital terrain models (DTMs) were generated from the ground points and digital surface models (DSMs) were produced from the points of all classes for both datasets, with a resolution of 0.05 m. Additionally, orthophotos were exported for both datasets (Figures 2 and 3). Before any process, we generated normalized canopy height models (nCHMs) by subtracting the DSM from the DTM in ArcMap 10.3.1 by ESRI © .

Spectral and Texture Analysis
Once the positions of all the trees in the study area were identified, a new point shapefile with the tree locations was created. To develop a deeper understanding of the spatial and temporal changes in the pattern ( Figure 4) caused by seasonal changes (spring and summer 2019), several VIs including a variety of wavelength bands (RGB, red-edge, and NIR) were applied (Table 3) to evaluate the amount of vegetation cover change between the two seasons. The PCI geomatica, version 2016.04.01 (PCI Geomatics, Markham, ON, Canada) was used to map the TA within the zonal areas from topographic information (nCHM). For each season, 13 nCHM-based TA layers were mapped, including Mean, Variance, Entropy, Contrast, Heterogeneity, Homogeneity, Angular Second Moment, Correlation, Gray-level Difference Vector (GLDV) Angular Second Moment, GLDV Entropy, GLDV Mean, GLDV Contrast, and Inverse Difference, with a kernel window size of 25 × 25 (32-bit). The main advantage of TA application on nCHMs is related to the fact that seasonal vegetation changes, as well as tree health status, can ultimately affect the structure of tree crown. For this reason, we assumed that the application of TA can lead to an optimum combined classification of forest tree species and health status assessment. A total of 26 TA layers were derived from nCHM for TA. Ultimately, to avoid overlapped tree crown areas and any uncertainties in the extracted values, a buffer zone with 2-m radius was used by means of zonal statistics in ArcMap ( Figure 5). An overview of the methodology can be seen in Figure 6.

Spectral and Texture Analysis
Once the positions of all the trees in the study area were identified, a new point shapefile with the tree locations was created. To develop a deeper understanding of the spatial and temporal changes in the pattern (Figure 4) caused by seasonal changes (spring and summer 2019), several VIs including a variety of wavelength bands (RGB, red-edge, and NIR) were applied (Table 3) to evaluate the amount of vegetation cover change between the two seasons.

Spectral and Texture Analysis
Once the positions of all the trees in the study area were identified, a new point shapefile with the tree locations was created. To develop a deeper understanding of the spatial and temporal changes in the pattern ( Figure 4) caused by seasonal changes (spring and summer 2019), several VIs including a variety of wavelength bands (RGB, red-edge, and NIR) were applied (Table 3) to evaluate the amount of vegetation cover change between the two seasons. The PCI geomatica, version 2016.04.01 (PCI Geomatics, Markham, ON, Canada) was used to map the TA within the zonal areas from topographic information (nCHM). For each season, 13 nCHM-based TA layers were mapped, including Mean, Variance, Entropy, Contrast, Heterogeneity, Homogeneity, Angular Second Moment, Correlation, Gray-level Difference Vector (GLDV) Angular Second Moment, GLDV Entropy, GLDV Mean, GLDV Contrast, and Inverse Difference, with a kernel window size of 25 × 25 (32-bit). The main advantage of TA application on nCHMs is related to the fact that seasonal vegetation changes, as well as tree health status, can ultimately affect the structure of tree crown. For this reason, we assumed that the application of TA can lead to an optimum combined classification of forest tree species and health status assessment. A total of 26 TA layers were derived from nCHM for TA. Ultimately, to avoid overlapped tree crown areas and any uncertainties in the extracted values, a buffer zone with 2-m radius was used by means of zonal statistics in ArcMap ( Figure 5). An overview of the methodology can be seen in Figure 6.  The PCI geomatica, version 2016.04.01 (PCI Geomatics, Markham, ON, Canada) was used to map the TA within the zonal areas from topographic information (nCHM). For each season, 13 nCHM-based TA layers were mapped, including Mean, Variance, Entropy, Contrast, Heterogeneity, Homogeneity, Angular Second Moment, Correlation, Gray-level Difference Vector (GLDV) Angular Second Moment, GLDV Entropy, GLDV Mean, GLDV Contrast, and Inverse Difference, with a kernel window size of 25 × 25 (32-bit). The main advantage of TA application on nCHMs is related to the fact that seasonal vegetation changes, as well as tree health status, can ultimately affect the structure of tree crown. For this reason, we assumed that the application of TA can lead to an optimum combined classification of forest tree species and health status assessment. A total of 26 TA layers were derived from nCHM for TA. Ultimately, to avoid overlapped tree crown areas and any uncertainties in the extracted values, a buffer zone with 2-m radius was used by means of zonal statistics in ArcMap ( Figure 5). An overview of the methodology can be seen in Figure 6.

Spectral Analysis
The results of the multiple comparisons for both seasons showed that there was a significant difference in the amount of VIs between the major classes. As we assumed due to unbalanced reference data, Levene`s test revealed that there was heterogeneity of variances between the classes for the majority of VI classification cases, as described in Table S1.
It is obvious from the post-hoc analysis in Figure 7 that NDVI, SAVI, and MCARI2 performed the best in classifying the reference data from spring 2019, whereas SAVI, MCARI2, and PSRI had the best performance for the summer period ( Figure 8). The above-mentioned VIs were able to successfully distinguish all the major classes from each other. Since the mentioned VIs are based on NIR, red, and red-edge wavelengths, they can detect the changes occurring in the amount of chlorophyll absorption and reflectance of the cellular structure of the leaves [67][68][69]. When the tree is dehydrated or afflicted with a disease, spongy layers within the leaf surface deteriorate, thus, a tree absorbs more NIR compared to a healthy tree [29][30][31].
The post-hoc results highlighted an issue for the differentiation between major classes in spring between healthy broadleaves (B) and dead conifers (D), which likely occurred because of the nature of broadleaves being without leaves during the spring. This issue resulted in confusion in the identification of major classes in the case of RESR, M3CI, CI, PSRI, and NDRE VI layers.
According to the coefficient of determination Eta-squared (η 2 ), MCARI2 and CI respectively had the highest correlation with field data major classes with Eta-squared (η 2 ) of 0.73 in spring 2019 and 0.61 for summer 2019, as described in Table S2. These results confirm the importance of the red-edge band in tree type and health status assessment. Post-hoc results for the summer exhibited better performance of PSRI compared to NDVI in the classification of reference data in different major classes. The results of Eta-squared (η 2 ) were used to select the best independent variables, including VI and TA, for further analysis with the SVM.

Classification Properties for the Distinguishment of Tree Species and Health Classes
In addition to the spectral and texture analysis that was applied for the general classification, the SVM algorithm was adopted for our unbalanced sample data in STATISTICA for conducting detailed classification (Table 1). For improving the accuracy of the output in our study, the SVM classifier type 1 was used. The training and validation samples were randomly selected as 70% for training and 30% for validation of the modeling in each category. For mapping not linearly separable classes, SVM has different kernel functions, among which the radial basis function (RBF) has shown superiority as to other functions [61,62]. RBF kernel was examined in a fixed number of gamma that was calculated according to 1/number of independent variables in SVM [46]. The adjustment of the gamma value influences the shape of the separating hyperplane and its value depends on the data interval and distribution and differs from one dataset to another [63,64]. For selecting the best parameters, V-fold cross-validation (V = 10) was used for minimizing the error function and improving the classification results. The general idea of using V-fold cross-validation is to divide the overall sample into several V-folds (randomly drawn disjoint sub-samples). The same type of SVM analysis is then successively applied to the observations belonging to the V-1 folds, and the results of the analyses are applied to sample V (the sample or fold that was not used to fit the SVM model) to compute the error, usually defined as the sum-of-squared. The results for the V replications are averaged to yield a single measure model error of the stability of the respective model, i.e., the validity of the model for predicting unseen data.

Statistical Analysis
We used SPSS Statistics V. 24 (IBM Corporation, Armonk, NY, U.S.A), STATISTICA V.12 (Dell Inc., Round Rock, TX, USA) and MATLAB ® R2017b professional edition (MathWorks, MA, USA) for the statistical analysis and visualization of the collected data. Because our data were not normally distributed, Levene's test was used to assess the equality of variances for all the tested statistical variables between the major classes; healthy broadleaves = B, healthy conifers = C, infested conifers = I, dead conifers = D for each VI. The following descriptive statistics were calculated for each VI and TA within each tree zonal area: Min, Max, Range, Mean, Standard deviation, Sum, Median, Variance, Mid-range, Skewness, and Kurtosis.
To determine the strength of the association between the categorical (independent) variable and numerical (dependent) variables, the Eta (η) coefficient was used (Equation (1)). Eta-squared (η 2 ) was calculated to measure the coefficient of determination. Multiple comparisons were then applied to compare the means between the variables; the first one was Tuckey's honest significance difference (HSD) test (Equation (2)), assuming equal variance, and the second was Dunnett's T3 test. The latter does not assume any equality of variance for analyzing the differences between the groups (Equation (3)). A significance level (a) of 0.05 was used to determine whether or not the above tests were statistically significant.
where SS effect is the sum-of-squares between the groups and SS total is the total sum-of-squares in the crosstab.
where M i − M j is the difference between the pair of means. To calculate this, M i should be larger than M j , MS w is the mean square within, and n h is the number in the group or treatment.
where, t Dunnett refers to the critical value in the Dunnett critical value table, 2MS S/A is mean square within groups, and n is the sample size. Kappa index, overall accuracy, commission error and user's accuracy used to evaluate the accuracy of SVM predicted results, given by Equations (4)-(7).
where p 0 is the observed agreement and p e is the expected agreement by chance.
where OA stands for overall accuracy, CP stands for correctly predicted values, and N is the total number of test data.
where CE stands for commission error, IP stands for incorrect predicted values and NP is the number of trees predicted by SVM (correct or not correct) for each class, UA stands for user's accuracy. However, according to the study by Shao et al. [65], the use of Kappa is associated with a series of theoretical flaws-uncertainties for accuracy assessment. Thus, for supporting the effectiveness of our methodology the F-score [66], omission error (OE), and producer's accuracy (PA) were calculated using Equations (8)- (10).
where OE stands for omission error, IP stands for incorrect predicted values and NO is the number of observed trees for each class, PA is the producer's accuracy.

Spectral Analysis
The results of the multiple comparisons for both seasons showed that there was a significant difference in the amount of VIs between the major classes. As we assumed due to unbalanced reference data, Levene's test revealed that there was heterogeneity of variances between the classes for the majority of VI classification cases, as described in Table S1.
It is obvious from the post-hoc analysis in Figure 7 that NDVI, SAVI, and MCARI2 performed the best in classifying the reference data from spring 2019, whereas SAVI, MCARI2, and PSRI had the best performance for the summer period ( Figure 8). The above-mentioned VIs were able to successfully distinguish all the major classes from each other. Since the mentioned VIs are based on NIR, red, and red-edge wavelengths, they can detect the changes occurring in the amount of chlorophyll absorption and reflectance of the cellular structure of the leaves [67][68][69]. When the tree is dehydrated or afflicted with a disease, spongy layers within the leaf surface deteriorate, thus, a tree absorbs more NIR compared to a healthy tree [29][30][31].   showing the strength of association between the four major categorical classes (healthy broadleaves = B, healthy conifers = C, infested conifers = I, dead conifers = D) and the most significant descriptive statistics for each VI. Small letters (a,b,c,d) were used to denote significant differences between groups (p-value < 0.05). The post-hoc results highlighted an issue for the differentiation between major classes in spring between healthy broadleaves (B) and dead conifers (D), which likely occurred because of the nature of broadleaves being without leaves during the spring. This issue resulted in confusion in the identification of major classes in the case of RESR, M3CI, CI, PSRI, and NDRE VI layers.
According to the coefficient of determination Eta-squared (η 2 ), MCARI2 and CI respectively had the highest correlation with field data major classes with Eta-squared (η 2 ) of 0.73 in spring 2019 and 0.61 for summer 2019, as described in Table S2. These results confirm the importance of the red-edge band in tree type and health status assessment. Post-hoc results for the summer exhibited better performance of PSRI compared to NDVI in the classification of reference data in different major classes. The results of Eta-squared (η 2 ) were used to select the best independent variables, including VI and TA, for further analysis with the SVM.

Texture Analysis
Similar to the spectral analysis, multiple comparisons result for TA revealed that there was a significant difference in the amount of nCHM-based texture values between the major classes. Due to the unbalanced reference data, Levene's test demonstrated that there was heterogeneity of variances between the classes for the majority of texture values (Table S1). Eta-squared (η 2 ) exhibited that in most cases there was a weak or moderate relationship between predicted and observed values (e.g., Dissimilarity = 0.14, GLDV Angular Second Moment = 0.22, GLDV Entropy = 0.22, GLDV Mean = 0.14, Inverse Difference = 0.15, and Standard Deviation = 0.17), as described in Table S3. nCHM-based Mean was the only variable that showed a stronger relationship with the reference data Remote Sens. 2020, 12, 3722 13 of 21 compared to the rest of variables (Figures 9 and 10), with an Eta-squared (η 2 ) of 0.51 in spring 2019 and 0.28 for summer 2019, as described in Table S3. significant difference in the amount of nCHM-based texture values between the major classes. Due to the unbalanced reference data, Levene`s test demonstrated that there was heterogeneity of variances between the classes for the majority of texture values (Table S1). Eta-squared (η 2 ) exhibited that in most cases there was a weak or moderate relationship between predicted and observed values (e.g., Dissimilarity = 0.14, GLDV Angular Second Moment = 0.22, GLDV Entropy = 0.22, GLDV Mean = 0.14, Inverse Difference = 0.15, and Standard Deviation = 0.17), as described in Table S3. nCHMbased Mean was the only variable that showed a stronger relationship with the reference data compared to the rest of variables (Figures 9,10), with an Eta-squared (η 2 ) of 0.51 in spring 2019 and 0.28 for summer 2019, as described in Table S3.  Similar to the spectral analysis, multiple comparisons result for TA revealed that there was a significant difference in the amount of nCHM-based texture values between the major classes. Due to the unbalanced reference data, Levene`s test demonstrated that there was heterogeneity of variances between the classes for the majority of texture values (Table S1). Eta-squared (η 2 ) exhibited that in most cases there was a weak or moderate relationship between predicted and observed values (e.g., Dissimilarity = 0.14, GLDV Angular Second Moment = 0.22, GLDV Entropy = 0.22, GLDV Mean = 0.14, Inverse Difference = 0.15, and Standard Deviation = 0.17), as described in Table S3. nCHMbased Mean was the only variable that showed a stronger relationship with the reference data compared to the rest of variables (Figures 9,10), with an Eta-squared (η 2 ) of 0.51 in spring 2019 and 0.28 for summer 2019, as described in Table S3.

Detailed Tree Species Classification and Health Status Assessment by SVM
Based on Eta-squared (η 2 ) results, all the VI layers which had strong relation with the reference data were used for modeling process using SVM. Those  (Table 4). Modeling tree species using the above layers successfully distinguished the reference data into detailed classes, producing an F-score of 95.65% for broadleaves, 78.48% for Norway spruce and 53.75% for Scots pine, resulting in OA of 78.82% (Table 4). This result confirms similar UAS-based studies on tree species multispectral classification accuracies reported elsewhere [22][23][24]70]. To improve the model accuracy, fusion of the TA and VI layers was used for modeling the reference data. The criteria for the selection of TA layers for SVM were exclusively based on Eta-squared (η 2 ) results. For TA layers we chose layers with moderate and strong relationship with the reference data, such as Mean (Mean) spring/summer; STD (Range) summer; GLDV Angular Second Moment (Skewness) summer; GLDV Entropy (Range) summer; GLDV Angular Second Moment (Coefficient of Variation) summer. The combination of VI and TA layers successfully improved accuracy of modeling tree species using SVM by 2.36% (OA) ( Table 5). SVM was able to successfully separate broadleaves from the rest of the classes with an accuracy of 100%, resulting in both CE and OE of 0%. Similarly, SVM exhibited the ability to distinguish the Norway spruce from broadleaves and Scots pine trees with PA of 91.89%, resulting in an OE of 8.11%. Consequently, the OA for distinguishing reference data for tree species from SVM was 81.18%. Our method reliably distinguished between broadleaved and coniferous trees ( Table 5). We estimated a Kappa of 0.70, which suggested a strong (substantial) relationship between the reference data and the independent variables. In addition, the combination of selected layers improved the amount of F-score values of 100% for broadleaves, 80.95% for Norway spruce, and 61.90% for Scots pine, respectively. Considering only tree health status, combination of TA and VI layers could not improve the accuracy of health status classification. However, using solely VI layers the SVM algorithm distinguished the classes with a very high OA of 84.71% (Table 6). The result is similar with the result of Näsi et al. [12]. Moreover, the results from CE in all cases validate the significance of our proposed methodology for the classification of the reference data into reliable classes. More specifically, the PA for all tree health classes was found to be ≥80%. Similar to the tree species classification, the Kappa for the assessment of health status was 0.66, which also suggested a substantial relationship between the reference data and the independent variables. In addition, F-score values for broadleaves, Norway spruce, and Scots pine were 72.73%, 88.64%, 81.69%, respectively. Classification of tree species/health using the selected VI layers successfully distinguished the reference data into detailed classes, producing an OA of 65.17% (Table 7). The SVM analysis showed significant performance for dividing the reference data into more complex classes that contained tree species and health status simultaneously. Although SVM analysis showed significant performance for dividing the reference data to a more complex classes, it could not distinguish the different species among the dead trees (Table 7). Our results showed that SVM characteristics such as (a) a high generalization ability and acquisition of high modeling accuracies, (b) convexity of the cost function, (c) powerful in the treatment of high-dimensional data (e.g., hyperspectral data), and (d) limited effort required for architecture design and training step compared to other machine learning algorithms make SVM one of the most powerful non-parametric algorithm that can be significantly helpful for modeling complex sets of data derived from complex environments such as forests [61][62][63][64].
As evident in Table 8, combination of TA and VI layers showed an increase in OA by 4.24%. NSI and SPD classes were successfully distinguished with very high PA, 92% and 100%, respectively, which resulted in an OA of 69.41%. The proposed methodology for classifying SPH yielded poorer, yet comparable, performance when compared to the NSI and SPD. It is evident in Table 8 that a lack of sufficient reference data [71] in some of the classes (e.g., NSD) caused poor performance in the classification process using the SVM classifier. When combining tree species and health status, as expected, the Kappa was poorer yet comparable (0.59), which represented a moderate relationship between reference data and independent variables. As a result of data combination between texture and VI, we were able to improve significantly the distinguishment of coniferous species between healthy trees. The amount of F-score increased by 12.5% for healthy Norway spruce and 6.8% for the healthy Scots pine. Once more this result confirms that non-parametric algorithms such as SVM can be used as a powerful tool for modeling complex environments when a combination of unlimited thematic, spectral, and topographic layers are required [72]. Moreover, based on the result from Tables 7 and 8, the distinguishment of tree species among the dead trees (coniferous) was not possible because of the similarity of reflectance of the dead leaves.

Conclusions
Our primary interest was to investigate the capability of low-cost UAS-based bi-temporal multispectral photogrammetry using the MicaSense RedEdge-M sensor to detect and distinguish tree species and evaluate their health status in mixed forest conditions. Due to the ability of UAS to acquire data in ultra-high spatial and temporal resolutions, our method not only enabled mapping of the tree species and its current health state, but it also aids in monitoring the spread of disease, designating the possibility of either early-stage detection or new disease outbreaks to adjacent healthy trees. This would support precision forest management and local decision-making.
Based on a combination of a bi-temporal integrated spectral and texture information as well as supervised machine learning techniques, we aimed to extract the reflectance values within a buffer zone of 2-m radius from each treetop and not to map the canopy extent of each tree. Using the same buffer zone, we eliminated the impact of crown shape and vague pixels from the overlapped crown areas. The results indicated that the proposed methodology delivered satisfactory outputs and achieved sufficiently accurate and reliable tree species classification with an OA of 81.18%, and a health status assessment with an OA of 84.70%. Chlorophyll index red-edge (CI) had the highest correlation with field data major classes with Eta-squared (η 2 ) of 0.61 compared to the rest of vegetation indices for summer 2019. That confirms the importance of the red-edge band for detail classification of both tree species type and health status during the vegetation period, due to the sensitivity of the red-edge band to small variations in the chlorophyll content and consistent across most species and health status classes.
This study confirmed the high efficiency of SVM with limited and unbalanced training samples and it allows users to combine an unlimited amount of different types of data (e.g., categorical and continuous input layers). Consequently, the combination of VI and TA layers for combined tree species and health status increased the OA by 4.24%. Conclusively, this approach can provide a new dimension of information derived from UAS platforms, which can lead to a better understanding of health status and tree species composition with reliable accuracy. Moreover, the results revealed that textural parameters of healthy tree crowns can be of significant importance in tree species classification. Distinguishment of species between dead trees could not be achieved due to the similar reflectance. This fact might be associated with the amount of error found in tree species classification between conifers, possibly due to the inability of the SVM algorithm to assign the dead conifer trees to one of the two coniferous classes. That is why higher accuracy has been found in broadleaves compared to the coniferous classes.
Limitations in the sampling from certain classes by multispectral sensors account for any of the uncertainties associated with our proposed approach. Identification of physiological stress in earlier stages or detection of smaller spectral differences could be possible with higher reliability and accuracy using hyperspectral sensors in wavelengths that cannot be covered by multispectral sensors. We also strongly believe that the use of thermal sensors could improve the classification of different tree species and health status. Nevertheless, UAS paired with multispectral sensors can significantly reduce the operational costs of monitoring forested areas on small to medium scales.
Author Contributions: A.A. conceived the idea for this manuscript; D.P. and A.A. designed and wrote the manuscript; A.A. processed the remote sensing and terrestrial data. D.P. and A.A. designed and performed the statistical analyses; D.P. and A.A. finalized the manuscript. All authors have read and agreed to the published version of the manuscript.