Discriminant Analysis of the Damage Degree Caused by Pine Shoot Beetle to Yunnan Pine Using UAV-Based Hyperspectral Images

Due to the increased frequency and intensity of forest damage caused by diseases and pests, effective methods are needed to accurately monitor the damage degree. Unmanned aerial vehicle (UAV)-based hyperspectral imaging is an effective technique for forest health surveying and monitoring. In this study, a framework is proposed for identifying the severity of damage caused by Tomicus spp. (the pine shoot beetle, PSB) to Yunnan pine (Pinus yunnanensis Franch) using UAV-based hyperspectral images. Four sample plots were set up in Shilin, Yunnan Province, China. A total of 80 trees were investigated, and their hyperspectral data were recorded. The spectral data were subjected to a one-way ANOVA. Two sensitive bands and one sensitive parameter were selected using Pearson correlation analysis and stepwise discriminant analysis to establish a diagnostic model of the damage degree. A discriminant rule was established to identify the degree of damage based on the median value between different degrees of damage. The diagnostic model with R690 and R798 as variables had the highest accuracy (R2 = 0.854, RMSE = 0.427), and the test accuracy of the discriminant rule was 87.50%. The results are important for forest damage caused by the PSB.


Introduction
Pine shoot beetles (PSBs) [1], Tomicus spp. are important forestry pests that wreak havoc on pine trees and pose a threat to ecological stability. The pests are difficult to find because they can remain hidden in trees for long periods. The damage stage of Tomicus spp. can be divided into two stages, which are the shoot-feeding stage and stem-boring stage. From May to November each year, the adults will migrate to the current year's shoots after eclosion, which is called the shoot feeding stage. From November to April of the following year, the adults consume the phloem and xylem of the trunk. Then, they mate and lay eggs in the trunk, which is called the stem-boring stage [2][3][4]. PSBs are among the main pests of Yunnan pine (Pinus yunnanensis Franch). Since the first large-scale outbreak in Yunnan, these pests have damaged Yunnan pine in dozens of cities, killing large numbers of pine forests and affecting the forestry industry and local economies [5].
Currently, visual investigation is a common method used for surveying forest diseases and pests, but this method is expensive in terms of manpower and material resources. It is difficult for researchers to navigate deep into the forest hinterland if they are using visual investigation, and the accuracy is

Data Acquisition
We used a DJI M600 Pro multi-rotor UAV (DJI, Shenzhen, China) equipped with a UHD S185 imaging spectrometer (Cubert GmbH, Ulm, Baden-Württemberg, Germany) to form a UAV-based hyperspectral imaging system. The global position system and inertial measurement unit modules were integrated onto the UAV with a horizontal and vertical position error of approximately 2.0 m and 5.0 m, respectively, with a positioning accuracy of approximately one degree [18,23]. For smallarea image analysis, the relative position is more important than the absolute horizontal and vertical position. Therefore, these error ranges are acceptable in small-area forestry surveys. The UHD S185 spectrometer can obtain accurate hyperspectral cube data over the entire field of view instantaneously with its Snapshot imaging technology. The main parameters of the UHD S185 that we used are listed in Table 2. The data acquisition route was planned by DJI Go Pro (DJI, Shenzhen China) based on Google Earth information.

Data Acquisition
We used a DJI M600 Pro multi-rotor UAV (DJI, Shenzhen, China) equipped with a UHD S185 imaging spectrometer (Cubert GmbH, Ulm, Baden-Württemberg, Germany) to form a UAV-based hyperspectral imaging system. The global position system and inertial measurement unit modules were integrated onto the UAV with a horizontal and vertical position error of approximately 2.0 m and 5.0 m, respectively, with a positioning accuracy of approximately one degree [18,23]. For small-area image analysis, the relative position is more important than the absolute horizontal and vertical position. Therefore, these error ranges are acceptable in small-area forestry surveys. The UHD S185 spectrometer can obtain accurate hyperspectral cube data over the entire field of view instantaneously with its Snapshot imaging technology. The main parameters of the UHD S185 that we used are listed in Table 2. The data acquisition route was planned by DJI Go Pro (DJI, Shenzhen China) based on Google Earth information.
UAV-based data acquisition was carried out in the study areas on 17 November 2019 from 11:00 a.m. to 2:00 p.m. under sunny and windless conditions. Taking into account the actual terrain, the condition of the vegetation in the study area and the study area to be covered, we flew the UAV at an altitude of 100 m above the ground, and the estimated ground sampling distances (GSD) were 0.028 m and 0.56 m for panchromatic and hyperspectral images, respectively. The flight speed was set at 5 m/s. Eight sorties were flown, with the image forward and side overlaps were set to 70% and 60% for each flight, until the entire study area was covered. Before each flight, the UHD S185 scanned the Spectralon ® panel (LabSphere, North Sutton, NH, USA) for giving an all-white calibration image with a reflectance of nearly 100% (Figure 2), and then the lens cap was covered and scanned to give an all-black calibration with 0% reflectivity. The above operations were for the correction of hyperspectral cube data. The acquisition area was approximately 0.23 km 2 . The camera was used to take the photos. Moreover, the position and orientation system (POS) were synchronized to facilitate hyperspectral data correction. When the UAV aerial photography was completed, the aerial photograph was screened to delete the photographs during take-off and landing, and only the photographs of flights within the route plan were retained. Table 2. Main parameters of the UHD S185 (provided by the manufacturer). UAV-based data acquisition was carried out in the study areas on 17 November 2019 from 11:00 a.m. to 2:00 p.m. under sunny and windless conditions. Taking into account the actual terrain, the condition of the vegetation in the study area and the study area to be covered, we flew the UAV at an altitude of 100 m above the ground, and the estimated ground sampling distances (GSD) were 0.028 m and 0.56 m for panchromatic and hyperspectral images, respectively. The flight speed was set at 5 m/s. Eight sorties were flown, with the image forward and side overlaps were set to 70% and 60% for each flight, until the entire study area was covered. Before each flight, the UHD S185 scanned the Spectralon ® panel (LabSphere, North Sutton, NH, USA) for giving an all-white calibration image with a reflectance of nearly 100% (Figure 2), and then the lens cap was covered and scanned to give an all-black calibration with 0% reflectivity. The above operations were for the correction of hyperspectral cube data. The acquisition area was approximately 0.23 km 2 . The camera was used to take the photos. Moreover, the position and orientation system (POS) were synchronized to facilitate hyperspectral data correction. When the UAV aerial photography was completed, the aerial photograph was screened to delete the photographs during take-off and landing, and only the photographs of flights within the route plan were retained. For geometric correction of the hyperspectral data and acquisition of sample trees' regions of interest (ROIs), a Red-Green-Blue (RGB) orthophoto image of the study area were acquired using DJI Phantom 4 RTK (DJI, Shenzhen China). The flight height was 100 m and the speed was 6 m/s. A total of four sorties were flown. Each flight required calibration photographs before take-off and after landing. DJI Go Pro was also used to plan the acquisition range of the RGB orthophoto image, which corresponded to the acquisition range of the hyperspectral data. However, the acquisition range of the RGB orthophoto image was larger than that of hyperspectral data for subsequent geometric correction of hyperspectral data.

Data Preprocessing
Rigorous preprocessing was required for the hyperspectral data in order to derive information from the imagery. The preprocessing of hyperspectral data is shown in Figure 3. The mosaic of RGB images was performed with DJI Terra 2.0 (DJI, Shenzhen, China), which allowed quick stitching of photographs captured by DJI Phantom 4 RTK to obtain RGB orthophoto image. The Cubert-Pilot software (Cubert GmbH, Ulm, Baden-Württemberg, Germany) was used to fuse each hyperspectral image with the corresponding panchromatic image acquired at the same time to obtain the fused For geometric correction of the hyperspectral data and acquisition of sample trees' regions of interest (ROIs), a Red-Green-Blue (RGB) orthophoto image of the study area were acquired using DJI Phantom 4 RTK (DJI, Shenzhen China). The flight height was 100 m and the speed was 6 m/s. A total of four sorties were flown. Each flight required calibration photographs before take-off and after landing. DJI Go Pro was also used to plan the acquisition range of the RGB orthophoto image, which corresponded to the acquisition range of the hyperspectral data. However, the acquisition range of the RGB orthophoto image was larger than that of hyperspectral data for subsequent geometric correction of hyperspectral data.

Data Preprocessing
Rigorous preprocessing was required for the hyperspectral data in order to derive information from the imagery. The preprocessing of hyperspectral data is shown in Figure 3. The mosaic of RGB images was performed with DJI Terra 2.0 (DJI, Shenzhen, China), which allowed quick stitching of photographs captured by DJI Phantom 4 RTK to obtain RGB orthophoto image. The Cubert-Pilot software (Cubert GmbH, Ulm, Baden-Württemberg, Germany) was used to fuse each hyperspectral image with the corresponding panchromatic image acquired at the same time to obtain the fused hyperspectral image. Then, the hyperspectral images were stitched using Agisoft PhotoScan (Agisoft LLC, Saint Petersburg, Russia) with the aid of point cloud data from the panchromatic images to generate hyperspectral orthophoto image. Geometric correction of hyperspectral data was carried out using ENVI's Image Registration Workflow via the RGB orthophoto image of the study area. Radiometric calibration was directly performed for hyperspectral data after geometric correction, and thus, was provided by the manufacturer. The radiometric calibration method described in Equation (1) (the formula is provided by the manufacturer): where ρ i is the standard reflection of band i; G is a collection of pixels in a white area of a mask image; m is the number of pixels in the collection; H i is the spectral radiance of the calibrated image after radiation correction in band i; and F i is the reflectivity of band i after radiation calibration.
Forests 2020, 11, x FOR PEER REVIEW 5 of 22 hyperspectral image. Then, the hyperspectral images were stitched using Agisoft PhotoScan (Agisoft LLC, Saint Petersburg, Russia) with the aid of point cloud data from the panchromatic images to generate hyperspectral orthophoto image. Geometric correction of hyperspectral data was carried out using ENVI's Image Registration Workflow via the RGB orthophoto image of the study area. Radiometric calibration was directly performed for hyperspectral data after geometric correction, and thus, was provided by the manufacturer. The radiometric calibration method described in Equation (1) (the formula is provided by the manufacturer): where ρi is the standard reflection of band i; G is a collection of pixels in a white area of a mask image; m is the number of pixels in the collection; Hi is the spectral radiance of the calibrated image after radiation correction in band i; and Fi is the reflectivity of band i after radiation calibration.

Spectral Reflectance of Individual Yunnan Pine Canopy
Based on the canopy location information (latitude and longitude) and photographs of each sample tree, the location and size of each canopy was determined on the RGB orthophoto image and the canopy range was plotted ( Figure 4) as the ROI. Because the canopy range of a severely damaged tree was difficult to identify, we used the contour of the branches instead of the canopy range. These ROIs were added to the pre-processed hyperspectral image. The average reflectance of the signal ROI was extracted via ENVI (Harris Corporation, Melbourne, FL, USA) as the spectrum of an individual tree for the subsequent analysis.

Spectral Reflectance of Individual Yunnan Pine Canopy
Based on the canopy location information (latitude and longitude) and photographs of each sample tree, the location and size of each canopy was determined on the RGB orthophoto image and the canopy range was plotted ( Figure 4) as the ROI. Because the canopy range of a severely damaged tree was difficult to identify, we used the contour of the branches instead of the canopy range. These ROIs were added to the pre-processed hyperspectral image. The average reflectance of the signal ROI was extracted via ENVI (Harris Corporation, Melbourne, FL, USA) as the spectrum of an individual tree for the subsequent analysis.

First Derivative and Spectral Parameters
The quality of the raw data was poor because the data collection was affected by environmental factors such as soil and light. Therefore, the spectral reflectance curves were smoothed using the Savitzky-Golay (S-G) method [24]. The method eliminates the effects of noise and improves the data quality while preserving as much of the original spectral information as possible. In this study, the original spectral reflectance curve was preprocessed by the 7-point and 3-term S-G method.
The first derivative of the spectral reflectance curve can reflect the position of spectral peaks and troughs very well and plays an extremely important role in the processing of spectral data. Therefore, the first derivative can reflect the slight spectral changes caused by PSB, and thus, it is helpful for monitoring the damage degree. The first derivative [25] was calculated by Equation (2): where λi is the wavelength at i; is the spectral value of wavelength λi; and 2△λ is the difference in the wavelength of λi−1 to λi+1, determined by the spectrum sample interval.
The vegetation index is based on the spectral reflectance characteristics of plants, using the spectral reflectance of two or more wavelengths in various combinations and operations to enhance and extract vegetation information, in order to evaluate the growth condition of vegetation by reflecting information related to vegetation physiology or physics through the vegetation index [26]. Many physiological characteristics of Yunnan pine were altered after PSB infestation. The spectral reflectance information of the canopy was also altered. In this study, 28 parameters were selected with reference to the vegetation indices used in previous studies, combining the canopy spectral curves and first derivative features ( Table 3).
The spectral curves of each sample tree were obtained after data preprocessing. Pearson correlation analysis of spectral data with the damage degree was performed using the R programming language (AT&T Bell Laboratories, Auckland, New Zealand), and stepwise discriminant analysis was carried out using the Statistical Product and Service Solutions (International Business Machines Corporation, Palo Alto, CA, USA).

First Derivative and Spectral Parameters
The quality of the raw data was poor because the data collection was affected by environmental factors such as soil and light. Therefore, the spectral reflectance curves were smoothed using the Savitzky-Golay (S-G) method [24]. The method eliminates the effects of noise and improves the data quality while preserving as much of the original spectral information as possible. In this study, the original spectral reflectance curve was preprocessed by the 7-point and 3-term S-G method.
The first derivative of the spectral reflectance curve can reflect the position of spectral peaks and troughs very well and plays an extremely important role in the processing of spectral data. Therefore, the first derivative can reflect the slight spectral changes caused by PSB, and thus, it is helpful for monitoring the damage degree. The first derivative [25] was calculated by Equation (2): where λ i is the wavelength at i; ρ (λ i ) is the spectral value of wavelength λ i ; and 2 λ is the difference in the wavelength of λ i−1 to λ i+1 , determined by the spectrum sample interval. The vegetation index is based on the spectral reflectance characteristics of plants, using the spectral reflectance of two or more wavelengths in various combinations and operations to enhance and extract vegetation information, in order to evaluate the growth condition of vegetation by reflecting information related to vegetation physiology or physics through the vegetation index [26]. Many physiological characteristics of Yunnan pine were altered after PSB infestation. The spectral reflectance information of the canopy was also altered. In this study, 28 parameters were selected with reference to the vegetation indices used in previous studies, combining the canopy spectral curves and first derivative features (Table 3). Table 3. Spectral parameters. Two types of spectral parameters were used in this study: vegetation indices and hyperspectral parameters. We also list specific calculation formulas and describe specific meanings. Normalized Difference Vegetation Index (NDVI); Red Edge Normalized Difference Vegetation Index (NDVI 705 ); Simple Ratio Index (SR); Photochemical Reflectance Index (PRI).

Variable Categories Variable Definitions and Formulas References
Vegetation indexes Hyperspectral parameters

Reflectance parameters
Rg (Green mountain): Maximum reflectance in the wavelength range of 520-560 nm [30] Rr First derivative parameters Db: The maximum value of the first derivative in the blue edge region (470-520 nm) [30] Dy: The maximum value of the first derivative spectra in the yellow edge region (560-590 nm) [30] Dr: The maximum value of the first derivative spectra in the red-edge region (660-740 nm) [31] Dnir: The maximum value of the first derivative spectra in near-infrared region (760-950 nm) [31] SDb: The sum of the first derivative values in the blue edge region [30] SDy: The sum of the first derivative values in the yellow edge region [30] SDr: The sum of the first derivative values in the red-edge region [30] SDnir: The sum of the first derivative values in the near-infrared region [31] SDr/SDb [31] SDr/SDy [31] SDnir/SDb [31] SDnir/SDr [31] (SDr − SDb)/(SDr + SDb) [31] (SDr − SDy)/(SDr + SDy) [31] Forests 2020, 11, 1258 8 of 22 The spectral curves of each sample tree were obtained after data preprocessing. Pearson correlation analysis of spectral data with the damage degree was performed using the R programming language (AT&T Bell Laboratories, Auckland, New Zealand), and stepwise discriminant analysis was carried out using the Statistical Product and Service Solutions (International Business Machines Corporation, Palo Alto, CA, USA).

Ground Survey and Individual Yunnan Pine Canopy Survey
A field survey of Yunnan pine in the study area was conducted on 18 November 2019. Four sample plots were set up, and 20 trees were selected for each damage degree. There were a total of 80 sample trees, evenly distributed across the study area. For each sample tree, the number of damaged shoots and total shoots was estimated by visual rough estimation, and the latitude and longitude positions of each tree were recorded. To obtain the relationship between the damaged shoot ratio (DSR) of the canopy (Canopy-DSR) and the DSR of ground survey (Ground-DSR), this study used an M200+Z30 ((DJI, Shenzhen, China) camera to collect canopy data of Yunnan pine. A total of 80 samples were collected with clear pictures of the canopy ( Figure 5), which allowed clear visual identification of the damaged shoot. By visual inspection, we obtained the total number of shoots and the number of damaged shoots per canopy. Finally, we calculated the DSR. 80 sample trees, evenly distributed across the study area. For each sample tree, the number of damaged shoots and total shoots was estimated by visual rough estimation, and the latitude and longitude positions of each tree were recorded. To obtain the relationship between the damaged shoot ratio (DSR) of the canopy (Canopy-DSR) and the DSR of ground survey (Ground-DSR), this study used an M200+Z30 ((DJI, Shenzhen, China) camera to collect canopy data of Yunnan pine. A total of 80 samples were collected with clear pictures of the canopy ( Figure 5), which allowed clear visual identification of the damaged shoot. By visual inspection, we obtained the total number of shoots and the number of damaged shoots per canopy. Finally, we calculated the DSR.
The DSR of one individual Yunnan pine can be calculated by Equation (3): where n is the number of shoots affected by PSB and m is the total number of shoots of an individual Yunnan pine. Canopy-DSR represents the DSR of one individual Yunnan pine canopy obtained by visual inspection and Ground-DSR represents the DSR of one individual Yunnan pine obtained by ground survey. Each Canopy-DSR has a corresponding Ground-DSR (Table 4).

Model Building and Accuracy Validation
Seventy percent of the data volume was used for model building and 30% for accuracy validation [32]. A Multivariate Linear Regression (MLR, Equation (4)) model was used to fit the equations [33]: where Y is the dependent variable; X1 are independent variables; β0 is the intercept; βi is the parameters to be estimated by linear regression analysis; and ε is the additive term of the error. In this study, the coefficient of determination (R 2 , Equation (5)), the root-mean-square error (RMSE, Equation (6)), the relative-root-mean-square error (rRMSE, Equation (7)), and the Bias (Equation (8)) were employed as the accuracy verification indicators. The value of R 2 is between 0 The DSR of one individual Yunnan pine can be calculated by Equation (3): where n is the number of shoots affected by PSB and m is the total number of shoots of an individual Yunnan pine. Canopy-DSR represents the DSR of one individual Yunnan pine canopy obtained by visual inspection and Ground-DSR represents the DSR of one individual Yunnan pine obtained by ground survey. Each Canopy-DSR has a corresponding Ground-DSR (Table 4).

Model Building and Accuracy Validation
Seventy percent of the data volume was used for model building and 30% for accuracy validation [32]. A Multivariate Linear Regression (MLR, Equation (4)) model was used to fit the equations [33]: where Y is the dependent variable; X 1 are independent variables; β 0 is the intercept; β i is the parameters to be estimated by linear regression analysis; and ε is the additive term of the error. In this study, the coefficient of determination (R 2 , Equation (5)), the root-mean-square error (RMSE, Equation (6)), the relative-root-mean-square error (rRMSE, Equation (7)), and the Bias (Equation (8)) were employed as the accuracy verification indicators. The value of R 2 is between 0 and 1. If R 2 is close to 1, the reference value of the prediction model is high. The smaller the RMSE value was, the more accurate the prediction model was.
where n is the sample size; y i is the measured value at location i;ŷ i is the predicted value at location i; y is the average value of y; and SST and SSE are the total sum of squares and residual sum of squares, respectively.

The Relationship between Canopy-DSR and Ground-DSR
Pearson correlation analysis was carried out between Canopy-DSR and Ground-DSR [34]. The results showed that Ground-DSR had an extremely significant positive linear correlation with Canopy-DSR, and the Pearson correlation coefficient (PCC, r) was 0.995. Using the proportions of 70% and 30%, 80 sampled trees were divided into a modeling dataset (56 sample trees) and a prediction dataset (24 sample trees). Unitary linear regression was performed to the modeling set with Canopy-DSR as the independent variable and Ground-DSR as the dependent variable; and the equation (y = 0.007 + 1.005x) was obtained (Table 5). We then performed a linear fit between the measured and predicted values, and the results are shown in Figure 6. Table 5. Fitting equation of Ground-DSR (n = 56) and accuracy validation (n = 24). y is the Ground-DSR obtained by the fitting equation, and x is the Canopy-DSR. Damaged shoot ratio (DSR); Root-mean-square error (RMSE); Relative-root-mean-square error (rRMSE). As can be seen from Table 5 and Figure 6, Canopy-DSR can be used to fit Ground-DSR, and the accuracy was high (R 2 > 0.95, RMSE < 0.05 and Bias < 0.001). Therefore, Ground-DSR can be retrieved by Canopy-DSR to determine the damage degree to Yunnan pine. Figure 7a illustrates the spectral reflectance of the Yunnan pine canopy with different damage degrees, showing the spectral changes. The spectral reflectance curves had significant peaks in the wavelength range of 540-560 nm, and the peak decreased with an increase of the damage degree. The magnitudes were in the order of healthy, mild damage, moderate damage, and severe damage. The reflectance curve showed obvious troughs at 660 to 680 nm, and the value of the troughs increased with the severity of the damage. The reflectance of the healthy Yunnan pine canopy was the lowest, and the reflectance of the other three damage degrees showed no significant differences. At 760 to 946 nm, the reflectance of Yunnan pine canopy decreased as the damage degree increased, with clear differences. After 890 nm, the difference between the different damage degrees decreased and the reflectance values gradually converged.

The Relationship between Canopy-DSR and Ground-DSR
Pearson correlation analysis was carried out between Canopy-DSR and Ground-DSR [34]. The results showed that Ground-DSR had an extremely significant positive linear correlation with Canopy-DSR, and the Pearson correlation coefficient (PCC, r) was 0.995. Using the proportions of 70% and 30%, 80 sampled trees were divided into a modeling dataset (56 sample trees) and a prediction dataset (24 sample trees). Unitary linear regression was performed to the modeling set with Canopy-DSR as the independent variable and Ground-DSR as the dependent variable; and the equation (y = 0.007 + 1.005x) was obtained (Table 5). We then performed a linear fit between the measured and predicted values, and the results are shown in Figure 6. Table 5. Fitting equation of Ground-DSR (n = 56) and accuracy validation (n = 24). y is the Ground-DSR obtained by the fitting equation, and x is the Canopy-DSR. Damaged shoot ratio (DSR); Rootmean-square error (RMSE); Relative-root-mean-square error (rRMSE).    As can be seen from Table 5 and Figure 6, Canopy-DSR can be used to fit Ground-DSR, and the accuracy was high (R 2 > 0.95, RMSE < 0.05 and Bias < 0.001). Therefore, Ground-DSR can be retrieved by Canopy-DSR to determine the damage degree to Yunnan pine. Figure 7a illustrates the spectral reflectance of the Yunnan pine canopy with different damage degrees, showing the spectral changes. The spectral reflectance curves had significant peaks in the wavelength range of 540-560 nm, and the peak decreased with an increase of the damage degree. The magnitudes were in the order of healthy, mild damage, moderate damage, and severe damage. The reflectance curve showed obvious troughs at 660 to 680 nm, and the value of the troughs increased with the severity of the damage. The reflectance of the healthy Yunnan pine canopy was the lowest, and the reflectance of the other three damage degrees showed no significant differences. At 760 to 946 nm, the reflectance of Yunnan pine canopy decreased as the damage degree increased, with clear differences. After 890 nm, the difference between the different damage degrees decreased and the reflectance values gradually converged. The first derivative curves of the canopy spectra of Yunnan pine with four damage degrees are shown in Figure 7b. There were significant peaks at 510 to 530 nm, which gradually decreased with the level of damage, but the differences in the values were not significant. There were three troughs and two peaks at 550 to 670 nm, with values at the troughs increasing slightly with increasing damage, but the differences were not significant, and values at the peaks of the four damage degrees were not significantly different. At 700 to 720 nm, there were obvious wave peaks. As the damage degree increased, the peaks decreased significantly, and the position of the peak shifted slightly in the short-wave direction, i.e., there was a "red-edge blue shift."

Spectral Characteristics of the Damaged Yunnan Pine Canopy
The spectral reflectance ratio curves of the healthy/healthy, mild damage/healthy, moderate damage/healthy, and severe damage/healthy ratios can reflect the changes of spectral reflectance of different damage degrees compared to healthy Yunnan pine canopy (Figure 8). From Figure 8, we The first derivative curves of the canopy spectra of Yunnan pine with four damage degrees are shown in Figure 7b. There were significant peaks at 510 to 530 nm, which gradually decreased with the level of damage, but the differences in the values were not significant. There were three troughs and two peaks at 550 to 670 nm, with values at the troughs increasing slightly with increasing damage, but the differences were not significant, and values at the peaks of the four damage degrees were not significantly different. At 700 to 720 nm, there were obvious wave peaks. As the damage degree increased, the peaks decreased significantly, and the position of the peak shifted slightly in the short-wave direction, i.e., there was a "red-edge blue shift." The spectral reflectance ratio curves of the healthy/healthy, mild damage/healthy, moderate damage/healthy, and severe damage/healthy ratios can reflect the changes of spectral reflectance of different damage degrees compared to healthy Yunnan pine canopy (Figure 8). From Figure 8, we can see that the reflectance changes of mild damage, moderate damage, and severe damage in Yunnan pine canopies relative to healthy Yunnan pine canopies were roughly similar. At 510 to 570 nm and 695 to 950 nm, the spectral values of mild damage were lower than those of healthy trees. At 500 to 610 nm and 690 to 950 nm, the spectral values of moderate damage were lower than those of healthy trees. At 495 to 645 nm and 685 to 950 nm, the spectral values of severe damage were lower than those of healthy trees. This indicated the range of spectral reflectance below the healthy Yunnan pine canopy gradually increased and above the healthy Yunnan pine canopy gradually decreased as the damage degree increased. From the overall trend, the spectral reflectance of Yunnan pine canopy of the four damage degrees was in order of healthy, mild damage, moderate damage, and severe damage. Therefore, the damage degree can be determined by the spectral reflectance. of the four damage degrees was in order of healthy, mild damage, moderate damage, and severe damage. Therefore, the damage degree can be determined by the spectral reflectance.

Differential Analysis of Spectra in Yunnan Pine Canopy at Different Damage Degrees
The original spectra, S-G spectra, and first derivatives of the Yunnan pine canopy were analyzed by Pearson correlation with the damage degrees; the bands with the top 10 absolute correlation coefficients (Table 6) were selected for one-way analysis of variance (one-way ANOVA). In the analysis, healthy, mild, moderate, and severe damage were assigned values of 0, 1, 2, and 3, respectively. The results are shown in Figure 9. The bands of the original and S-G spectra were identical, both being 774-810 nm, in the near-infrared band. The band of the original spectrum with the largest correlation coefficient was 798 nm (r = −0.897 **), and those of the S-G spectra were 798 nm and 802 nm, respectively (r = −0.897 **). The bands of the first derivatives with the top 10 correlation coefficients were located around the red-edge region, 698-734 nm, and the band with the largest correlation coefficient was 714 nm (r = −0.922 **). However, the collinearity diagnosis showed that there was severe multicollinearity between these wavelengths, with tolerance approaching 0 and VIF much greater than 10. Therefore, the damage degree estimation model could not be established through these bands. From Figure 9, we could see that the spectral reflectance or first derivative values of the Yunnan pine canopy at different damage degrees were significantly different at these 10 bands (p < 0.05). This means that the damage degree of the Yunnan pine canopy can be differentiated by the values at these wavelengths, and thus, the hazard of an individual Yunnan pine tree can be differentiated.

Differential Analysis of Spectra in Yunnan Pine Canopy at Different Damage Degrees
The original spectra, S-G spectra, and first derivatives of the Yunnan pine canopy were analyzed by Pearson correlation with the damage degrees; the bands with the top 10 absolute correlation coefficients (Table 6) were selected for one-way analysis of variance (one-way ANOVA). In the analysis, healthy, mild, moderate, and severe damage were assigned values of 0, 1, 2, and 3, respectively. The results are shown in Figure 9. The bands of the original and S-G spectra were identical, both being 774-810 nm, in the near-infrared band. The band of the original spectrum with the largest correlation coefficient was 798 nm (r = −0.897 **), and those of the S-G spectra were 798 nm and 802 nm, respectively (r = −0.897 **). The bands of the first derivatives with the top 10 correlation coefficients were located around the red-edge region, 698-734 nm, and the band with the largest correlation coefficient was 714 nm (r = −0.922 **). However, the collinearity diagnosis showed that there was severe multicollinearity between these wavelengths, with tolerance approaching 0 and VIF much greater than 10. Therefore, the damage degree estimation model could not be established through these bands. From Figure 9, we could see that the spectral reflectance or first derivative values of the Yunnan pine canopy at different damage degrees were significantly different at these 10 bands (p < 0.05). This means that the damage degree of the Yunnan pine canopy can be differentiated by the values at these wavelengths, and thus, the hazard of an individual Yunnan pine tree can be differentiated.  A discriminant analysis of the damage degree was carried out via the band with the largest correlation coefficient by using Fisher's discriminant function and leave-one-out cross-validation. The results are shown in Table 7. When using S-G spectral values to distinguish the damage degree, only 798 nm was used for the analysis due to the strong collinearity between 798 nm and 802 nm. The accuracy of all three methods exceeded 70%, with first derivative discrimination being the most effective with an accuracy of 80%.  A discriminant analysis of the damage degree was carried out via the band with the largest correlation coefficient by using Fisher's discriminant function and leave-one-out cross-validation. The results are shown in Table 7. When using S-G spectral values to distinguish the damage degree, only 798 nm was used for the analysis due to the strong collinearity between 798 nm and 802 nm. The accuracy of all three methods exceeded 70%, with first derivative discrimination being the most effective with an accuracy of 80%.

Differential Analysis of Hyperspectral Parameters in Yunnan Pine Canopy at Different Damage Degrees
Pearson correlation analysis was performed for hyperspectral parameters and damage degrees. The hyperspectral parameters with the 10 highest absolute correlation coefficients were selected for one-way ANOVA and the results are shown in Figure 10  . There were significant differences in the four damage degrees of Dr, SDr, SDnir, and SDy (p < 0.05). Therefore, the damage degree of Yunnan pine canopy can be distinguished by these hyperspectral parameters. There were no significant differences in SDb or Db between healthy and mild damage, but these two damage degrees were significantly different from the other two. The values of (Rg − Rr)/(Rg + Rr), SR1, Rg/Rr, and SR2 showed no significant difference between mild damage and moderate damage, but these two damage degrees were significantly different from the other two damage degrees. Although some of the 10 hyperspectral parameters could not provide a good discrimination between all four damage degrees, they were able to distinguish severe damage from other degrees.

Screening of the Damage Degree-Sensitive Variables in Yunnan Pine Canopy
The original spectra (506-602 nm and 690-946 nm, 90 bands), S-G spectra (506-602 nm and 690-946 nm, 90 bands), first derivative values (498-546 nm, 554-618 nm, 626-662 nm, 670-758 nm, 802-830 nm, and 846-946 nm, for a total of 97 bands), and hyperspectral parameters (Dr, SDr, SDnir, SDy, SDb, Db, SR1, Rg/Rr, (Rg − Rr)/(Rg + Rr), SR2, NDVI1, NDVI2, NDVI705, Rg, SDnir/SDb, (SD − SDY)/(SDr + SDy), PRI, Dy, (SDnir − SDb)/(SDnir + SDb), (SDr + SDb), and SDr/SDb, for a total of 21 spectral parameters) of the Yunnan pine canopy, which were significantly correlated with the damage degree, were subjected to stepwise discriminant analysis using Fisher's discriminant function and Wilks's Lambda classification method. In the analysis, variables were introduced or eliminated through corresponding statistical tests. In the final discriminant function, only a small number of variables with high discriminative power were retained, while the variables with the lowest discriminative power were eliminated [35]. The results of the stepwise discriminant analysis are shown in Table 8. The sensitive bands of the original spectra were 690 nm and 798 nm, denoted as R 690 and R 798 , respectively. The sensitive bands of S-G were 690 nm and 802 nm, denoted as R S-G690 and R S-G802 . The sensitive bands of the first derivative were 650 nm and 714 nm, denoted as D 650 and D 714 . The sensitive parameter was Dr.
The leave-one-out cross-validation was applied to the accuracy test; the results are shown in Table 9. The modeling accuracy of all four spectral treatments exceeded 80%, with the highest accuracy of 86.30% for hyperspectral parameters. The first derivative had the highest cross-validation accuracy of 85.00%. According to the modeling accuracy and cross-validation results, the first derivatives were the best discriminated and the S-G spectra were the worst discriminated. The reason for this result may be that the average spectra of the Yunnan pine canopy was used for analysis in this study. The influences of the atmospheric environment, soil, and temperature on the spectral curve were balanced, and the noise level of the curves were small. On this basis, S-G filtering may cause some of the original information of the spectral curve to be lost, which may result in the lowest discriminative accuracy.
(Rg − Rr)/(Rg + Rr): r = −0.699 **, and SR2: r = −0.687 **). There were significant differences in the four damage degrees of Dr, SDr, SDnir, and SDy (p < 0.05). Therefore, the damage degree of Yunnan pine canopy can be distinguished by these hyperspectral parameters. There were no significant differences in SDb or Db between healthy and mild damage, but these two damage degrees were significantly different from the other two. The values of (Rg − Rr)/(Rg + Rr), SR1, Rg/Rr, and SR2 showed no significant difference between mild damage and moderate damage, but these two damage degrees were significantly different from the other two damage degrees. Although some of the 10 hyperspectral parameters could not provide a good discrimination between all four damage degrees, they were able to distinguish severe damage from other degrees.

Screening of the Damage Degree-Sensitive Variables in Yunnan Pine Canopy
The original spectra (506-602 nm and 690-946 nm, 90 bands), S-G spectra (506-602 nm and 690-946 nm, 90 bands), first derivative values (498-546 nm, 554-618 nm, 626-662 nm, 670-758 nm, 802-830 nm, and 846-946 nm, for a total of 97 bands), and hyperspectral parameters (Dr, SDr, SDnir, SDy, SDb, Db, SR1, Rg/Rr, (Rg − Rr)/(Rg + Rr), SR2, NDVI1, NDVI2, NDVI705, Rg, SDnir/SDb, (SD − SDY)/(SDr + SDy), PRI, Dy, (SDnir − SDb)/(SDnir + SDb), (SDr + SDb), and SDr/SDb, for a total of 21 spectral parameters) of the Yunnan pine canopy, which were significantly correlated with the damage degree, were subjected to stepwise discriminant analysis using Fisher's discriminant function and Wilks's Lambda classification method. In the analysis, variables were introduced or eliminated through corresponding statistical tests. In the final discriminant function, only a small number of variables with high discriminative power were retained, while the variables with the lowest discriminative power were eliminated [35]. The results of the stepwise discriminant analysis are shown in Table 8. The sensitive bands of the original spectra were 690 nm and 798 nm, denoted as R690 and R798, respectively. The sensitive bands of S-G were 690 nm and 802 nm, denoted as RS-G690  Table 8. Variable screening results of different spectral processing methods. Screening results were processed through stepwise discriminant analysis. PCCs are the results of the Pearson correlation analysis.

Method
Step Variable PCC The original spectra 1

Modeling the Damage Degree of Yunnan Pine Canopy
Based on the screened original spectra, S-G spectra, first-order derivative sensitive bands and hyperspectral sensitive parameters, MLR models were constructed. Moreover, diagnostic model accuracy was tested via R 2 , RMSE, rRMSE, and Bias.
R 690 and R 798 , R S-G690 and R S-G802 , D 650 and D 714 , and Dr were used to establish MLR models based on the modeling dataset, respectively. The prediction dataset was used for accuracy validation of the models (Table 10). The results showed that the original spectral model fitted best with R 2 , RMSE and Bias of 0.864, 0.416, and −0.068, respectively. The S-G spectral model gave the best prediction with R 2 , RMSE and Bias of 0.827, 0.465, and −0.067, respectively. The measured values were fitted linearly to the predicted values, and the results are shown in Figure 11, with R 2 exceeding 0.8 for all four methods. There were no significant differences among the four models, and all models could better predict the damage degree caused by PSB to Yunnan pine. Table 9. Stepwise discriminant analysis results for different spectral processing methods. Stepwise discriminant analysis was carried out on the results of four spectral data processing methods and the damage degree. The discriminant coefficient of each degree of damage is listed in Fisher's discriminant function. Constant represents the constant term in the discriminant function.

Quantitative Discriminant Rules of the Damage Degree of the Yunnan Pine Canopy
The quantification discriminant rules of the damage degree were established by the diagnostic model. First, the diagnostic model was used to calculate predicted values of the damage degree for the modeling dataset. Then, the average of the predicted values for each damage degree was calculated. The averages of the predicted values of the two damage degrees were used to calculate the median value between these two damage degrees, and this was used as the threshold to establish the discriminant rule. Finally, the damage degree of the prediction dataset was determined by the rule to verify the accuracy of the discriminant rule. The discriminant rules are shown in Table 11, and the accuracy test results are shown in Table 12. The discriminant rule based on the original spectra had the highest test accuracy of 87.50%; however, the first derivative was the lowest at 79.17%. This rule had the best discrimination for healthy Yunnan pine canopy and the worst discrimination for severe damaged Yunnan pine canopy. Nevertheless, healthy trees could be separated from damaged trees according to this rule.

Model
Regression

Discussion
Forest tree mortality is the result of multiple factors, making it difficult to predict and discern the causes of mortality [36]. In recent years, hyperspectral remote sensing has been widely applied in field of agriculture and forestry to rapidly obtain effective information concerning vegetation and crops. In our study, Yunnan pine defoliation and mortality were mainly related to PSB. Thus, this study focused on the spectral characteristics, sensitive variables, and discriminant rules for the Yunnan pine canopy after damage by the PSB. Using UAV-based hyperspectral imaging, we developed diagnostic models and discriminant rules for Yunnan pine following damage by PSB, highlighting the application of UAV tools in forest pest management and monitoring and early warning.

Spectral Characteristics of Damaged Yunnan Pine
In recent decades, a large number of Yunnan pine forests have killed by PSB, which affected the ecological stability and economic outputs of forests [37]. In the study area, PSB has been identified as the most serious forest pest of Yunnan pine. Forest trees affected by pests exhibit specific symptoms, such as yellowing and reddening of leaves, an increased proportion of withering and defoliated trees, a decrease in leaf photosynthetic efficiency, and sudden death [38]. The damage degree was determined by Ground-DSR, but only the vegetation canopy could be observed using satellite and UAV remote sensing. Therefore, we established a relationship between Ground-DSR and Canopy-DSR to determine the damage degree by Canopy-DSR. By analyzing the spectral data, we found that the spectra of damaged Yunnan pine had obvious changes in the green peak (510-560 nm), red valley (640-680 nm), red-edge region (680-740 nm), and near-infrared bands (740-950 nm), consistent with previous studies on the spectral responses of damaged vegetation [6,21,39]. Such spectral changes in leaves and tree canopies indicate that remote sensing could be useful for the monitoring of forest diseases and pests [40].

Sensitive Band and Parameter Screening
The development of hyperspectral technology has brought new vitality to the monitoring of forest tree diseases and pests. The using of hyperspectral data increases the amount of data storage and transmission as well as the rich spectral information, but it also causes some problems. The whole data processing is tedious and time-consuming because of a large amount of data noise and data redundancy [18]. These shortcomings bring obstacles to the popularization and application of hyperspectral technology. Hence, it is very important to streamline data and improve data processing efficiency. New indicators are formed by the combination of different original bands or function transformation, and sensitive bands are selected by statistical analysis or other methods according to the spectral differences of different monitoring targets [41][42][43]. In recent years, several new methods have emerged for the selection of hyperspectral bands or feature extraction [44], such as the screening of sensitive bands through principal component analysis, the successive projection algorithm, or other algorithms [44][45][46][47].
In this study, Pearson correlation analysis and stepwise discriminant analysis were used to screen sensitive bands and spectral parameters in the spectrum of damaged Yunnan pine. All selected bands were located in red and near-infrared regions, and the parameter (Dr) was also related to the red-edge position. These results are consistent with the characteristics of the spectral changes of damaged Yunnan pine. These bands and parameters were used to establish diagnostic models for the damage degree, with all R 2 being greater than 0.8, allowing monitoring of the damage caused by the PSB. Thus, the optimal band and parameter selection method were verified as suitable for the sensitive variable selection of defoliation and yellowing of Yunnan pine needles and can be applied to other pine stands.

Quantitative Discriminant Rule for the Damage Degree
When establishing the model, we assigned different values to the damage degrees, so it is necessary to establish quantitative discriminant rules to determine the prediction results of the model. The averages of the predicted values of the two different damage degrees were used to calculate the median value between these two damage degrees. Then, the median value was used as the threshold to establish the discriminant rule, and the discriminant accuracy of the degrees of damage reached 87.50%. The rule was most effective in identifying healthy Yunnan pine trees and could distinguish healthy Yunnan pine trees from damaged trees. This discriminant rule can also be applied to the monitoring of other diseases and pests. The method of establishing the discriminant rule based on the median value between two different classes as the threshold is simple to implement and it is effective in discriminating the damaged tree.

Monitoring Infestation by Tomicus spp.
With the development of remote sensing technology, more and more studies focus on the use of hyperspectral technology for vegetation health and forest management, especially UAV-based hyperspectral imaging for high-resolution and low-cost monitoring of forest pests [1,[17][18][19]. In addition, hyperspectral sensors are capable of identifying subtle spectral changes. In the study of defoliation of forest vegetation during the Dendrolimus tabulaeformis Tsai et Liu disaster outbreak, the sensitive spectral bands for defoliation rates were located in the green bands and near the red-edge region [18]. Furthermore, in another study of the damage caused by the European spruce bark beetle, it was found that the wavelengths closely related to the different damage classes were in the red and near-infrared bands [19]. In this study, we found that the sensitive bands for the damage degree caused by PSB were located near the red-edge region and in the near-infrared region. It is worth noting from these studies that the red-edge position and the near-infrared band are very sensitive to forest pests. In addition, the red-edge indices are important for the monitoring of forest disease and pest [17]. This has also been confirmed in our study. The spectral characteristic and discriminant analysis of single damaged trees by UAV-based hyperspectral technology will lay the foundation for low-cost, high-resolution monitoring of forest diseases and pests in large areas.
This study aimed to establish a framework to assess the damage degree to Yunnan pine affected by PSB as follows: (1) the selection of an optimal spectral pretreatment method and the establishment of a diagnostic model and quantitative discriminant rule using sensitive bands based on Pearson correlation analysis and stepwise discriminant analysis or (2) the establishment of a diagnostic model and quantitative discriminant rule using sensitive parameters based on the band combination of the spectrum. Although the framework might not be the best, it is superior to visual investigation in identifying the damage degree caused by PSB. In addition, it is both simple and efficient to implement. The highest discriminant accuracy of this framework for determining the damage degree caused by PSB reached 87.50% with only one or two variables. This result can be combined with multi-spectral satellite remote sensing or other effective methods for large-scale monitoring and accurate identification of damage caused by forest diseases and pests in the future.

Outlook
The monitoring framework for the damage degree caused by PSB is composed of some simple methods. We drew the tree canopy ROIs manually to extract the spectra. This process involves human participation, which may lead to a reduction in the identification accuracy, and it is not fully applicable to the identification of forest damage at the canopy level. Therefore, in follow-up research, an algorithm can be used to identify the tree canopy according to the textural features. Furthermore, LiDAR can be used to detect a single tree and obtain the tree canopy, which can be combined with UAV-based hyperspectral data for forest disease and pest monitoring [1]. The three-dimensional structure information of a single tree can be obtained according to the LiDAR data, which is also helpful for the monitoring of forest diseases and pests. Furthermore, LiDAR can also be used to monitor insect behavior [48]. Pearson correlation analysis and stepwise discriminant analysis were adopted in the screening of sensitive bands and parameters; these are simple methods. In the future, different algorithms can be used to acquire sensitive bands and parameters to improve the monitoring accuracy.

Conclusions
This study proposed a low-cost and miniaturized hyperspectral image data analysis system based on UAV remote sensing. Our study shows the potential of UAV-based hyperspectral remote sensing technology for forest disease and pest monitoring. This technology can be used to distinguish the anomalous reflectance characteristics of trees in small areas. Pearson correlation analysis and stepwise discriminant analysis were used to analyze the damage of Yunnan pine caused by PSB in order to obtain the sensitive bands and spectral parameters. First, the relationship between Ground-DSR and Canopy-DSR was verified, and it was found that there was a highly significant positive correlation and a strong linear relationship between the two indexes. We extracted the canopy spectra through the hand-drawn ROIs to analyze the spectral characteristics of damaged Yunnan pine. Then, Pearson correlation analysis and stepwise discriminant analysis were used to select sensitive bands and spectral parameters. Based on these sensitive variables, diagnostic models of the damage degree were established, and the quantitative discriminant rule was established using the median value as a threshold. The results showed that the diagnostic model with R 690 and R 798 as parameters was able to determine the damage degree of Yunnan pine using only two variables, and the test accuracy of the discriminant rule based on this model was the highest: 87.50%. This study established a simple framework to effectively identify the damage caused by PSB. In the future, we expect that the framework can be further improved to realize automatic identification and monitoring of the damage caused by PSB, combining LiDAR data, hyperspectral data, and satellite remote sensing data to achieve large-scale monitoring of PSB damage.

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