2.1.

Study Area

This study covered two different locations: Clinch Mountain in Russell County, Virginia, and the Om Sanctuary, in Asheville, North Carolina. Both sites have associated ground reference data derived from carbon inventories completed by foresters on the ground. These reference data were measured at specific plot locations within the forest at each site. The carbon inventory process details and associated ground reference data will be described in Sec. 2.2.

The Clinch Mountain study area covers $\sim 1646\mathrm{ha}$ and is rugged terrain composed of steep ridges and narrow valleys. Temperate deciduous species dominate the forests here, comprising mostly eastern USA broadleaf species. Reference data existed for 804 stems at this study site. The vegetation species included were predominantly deciduous types, comprising mostly the following species: maple species (Acer negundo, A. rubrum, and A. saccharum), yellow poplar (Liriodendron tulipifera), and oak species (Quercus alba, Q. coccinea, Q. muehlenbergii, Q. prinus, Q. rubra, Q. stellata, and Q. velutina) comprising 23.63%, 22.76%, and 18.66% (65.05% total) of the vegetation present, respectively. There were only three (0.37%) examples of eastern red cedar (Juniperus virginiana), the lone coniferous species, at this study site. This area served as an ideal location to highlight the use of ALS systems for gathering forestry information rapidly in difficult-to-traverse terrain. A map of this location can be seen in Fig. 1.

Fig. 1

The Clinch Mountain study site located near Lebanon, Virginia. The green points are the locations of the measured inventory plots. The majority of data used in this study came from this conservation site.

The Om Sanctuary was far smaller than the Clinch Mountain study area, covering only 17 ha. Similarly, reference data for a far smaller number of trees were available as well (149 trees). A map of the Om Sanctuary can be seen in Fig. 2. This location comprised a more diverse species mixture: yellow poplar (L. tulipifera), sourwood (Oxydendrum arboreum), chestnut oak (Q. prinus), and red maple (A. rubrum) comprising 26.17%, 15.44%, 11.41%, and 11.41% (64.43% total) of the vegetation present, respectively. A notable difference between the species population at the Om Sanctuary when compared with the Clinch Mountain site is the presence of conifers. Two species (Pinus strobus and P. rigidia) were present, totaling 11.41% of the vegetation. This is much larger than the 0.37% coniferous contribution at the other study site. This area served as a supplement to continue development of the sUAS models after data corruption (malfunctioning sensor) reduced the sample size of imaged plots in the Clinch Mountain location.

Fig. 2

The Om Sanctuary study site located just outside of Asheville, North Carolina. This was a small site, only covering 17 ha. The green points are the locations of the measured inventory plots.

This location also proved as challenging as the Clinch Mountain site, since it was located on a steep slope, where the sUAS staging site, to maintain legal line-of-sight (LOS), was separated from the forest by a river. As with the previous site, these forests were also of the temperate deciduous type, being mostly composed of eastern USA broadleaf trees.

2.2.

Reference Data Collection

Reference data were collected in the form of a carbon inventory by ground surveyors at both locations. This was completed in October 2016 at Clinch Mountain, and March 2017 at the Om Sanctuary. At the Clinch Mountain, the audit measured 151, 0.04 ha, fixed-radius plots (in slope-feet), spaced on a regular grid of 305 m across the mountain’s face, for a total inventory of 1949 individual stems. The plots at the Om Sanctuary were not in a regular grid but did utilize the same fixed-radius circular plot size, each covering an area of 0.04 ha. The total inventory measured at the Om Sanctuary was 804 trees.

For each location, diameter at the breast height (DBH) of vegetation $>5^{\prime \prime }$ ($\sim 12.7\mathrm{cm}$) and total height were measured, along with models produced of volume, above ground biomass (AGB), total carbon, and above ground live carbon (AGLC). Plot centers were located using a handheld global positioning system (GPS), with a positional accuracy of $\pm 2\mathrm{m}$. All measurements were collected using standard forest inventory operating procedures by forestry professionals. DBH was measured according to the United States Forest Service Forest Inventory and Analysis protocol, with a tape measure accuracy of $\pm 0.254\mathrm{cm}$ (0.1 in.). Tree height was measured using laser hypsometers and clinometers to derive the tree height from known distances and angles from the base and top of each tree. Tree height accuracy was $\pm 1\mathrm{m}$. These measurement accuracies are important as they are inputs into the allometric equations used to generate the volume, AGB, total carbon, and AGLC values reported as reference data. All aspects for the forest carbon inventory was re-sampled and audited by an independent third party for accuracy and was determined to be within the stated limits. Furthermore, all calculations were conducted to the 2014 California Environmental Protection Agency Air Resources Board Compliance Offset Protocol, US Forest Projects. The inventory resulted in a standard error of 4.98 tonnes/ac giving a 90% confidence interval of $+/-6.580\%$.

2.3.

Data Collection

The manned ALS data collection was completed by a commercial vendor on July 20, 2017, over the Clinch Mountain study site only. A Vexcel color (RGB) camera system and a Leica ALS80 small-footprint discrete return lidar, mounted on a fixed wing aircraft, was flown at a constant altitude $\sim 2400\mathrm{m}$ above mean sea level, and 1478 to 1814 m above the ground level. This ALS has a scanning rate of up to 200 Hz and a pulse rate of 1 MHz, with intensity values for the first three returns of each pulse and up to five returns/pulse being recorded. Point densities achieved using this sensor’s configuration ranged from 30.44 to $99.16\mathrm{pts}/{\mathrm{m}}^2$, with an average point density of $71.19\mathrm{pts}/{\mathrm{m}}^2$. The inconsistency between point densities was attributed to the inconsistent overlap of flight paths at the actual plot locations (plots at the center of the forest exhibited higher densities).

The sUAS data collections were completed on August 11 and 13, 2017, at the Clinch Mountain site and on August 12, 2017, at the Om Sanctuary. All flights were completed between 9h00 and 16h00 local time during (diffuse) overcast or clear sky conditions; there was not a noted performance degradation from the overcast conditions. The platform used for the sUAS was a Dà-Jiāng Innovations (DJI) Matrice M-600 Pro, with a custom payload rig capable of lifting a modular sensor platform; a lidar sensor and an RGB camera were flown for this study. A Velodyne VLP-16 small footprint discrete return lidar was utilized on the sUAS in conjunction with an Applanix APX-15 UAV as the GPS/IMU. The APX-15 employees real-time kinematic (RTK) correction to the GPS positioning data achieve position error bounds of 2 to 5 cm. This RTK correction was applied to the IMU data as well, achieving errors of 0.03 deg for roll and pitch, and 0.18 deg for yaw. System integration was performed by the UAS Research Lab in the Chester F. Carlson Center for Imaging Science at Rochester Institute of Technology. A depiction of the platform, as flown, can be seen in Fig. 3.

Fig. 3

The MX-1 platform developed by RIT. The Velodyne VLP-16 lidar is visible on the forward end of the peripheral mount at the bottom of the airframe.

The flight altitudes above the plots varied from 60 to 90 m above the ground level, with a variability due to a lack of precision in the DEM used by the DJI flight software in the terrain following mode; in-flight variations were $<3\mathrm{m}$. The VLP-16 lidar has a pulse rate of 300 KHz and a scan rate of up to 20 Hz. The lidar system has a maximum range of 100 m, the maximum swath width, assuming $-10\text{-}\mathrm{deg}$ to 10-deg sight lines like that of the ALS80, is 35 m. The VLP-16 is capable of 360-deg scanning, but it is not recommended to use scans at oblique angles for foliage mapping, as large scan angles may introduce error and occlusions.^29,41,42 The VLP-16 is a two return system (first and last return for every pulse), along with associated intensities. With this system, point densities achieved were 497.04 to $2393.93\mathrm{pts}/{\mathrm{m}}^2$, with a mean point density of $1504.2\mathrm{pts}/{\mathrm{m}}^2$. This variance in point density was attributed to the flight altitude, changes in flight overlap from plot-to-plot, and degradation in GPS timing signal, which resulted in returns being removed from the final point cloud. A summary of both lidar system specifications, as flown for both sites, can be found in Table 1.

Table 1

Lidar sensor specifications as flown for both sites.

2.4.

Structure-from-Motion Point Cloud Generation

Contemporaneous imagery for the SfM portion of the study was collected by both platforms. The manned aircraft deployed a four-channel, visible and near infrared, Vexcel UltraCam Eagle with a capture scheme of 60% stereoforward overlap between exposures and 30% side overlap between flight lines. The delivered imagery measured $24267\times 21274\text{pixels}$, with a mean GSD of 9 cm at the ground level. The sUAS system deployed an Allied Vision Mako G-419 RGB camera with a capture scheme of 80% forward overlap between exposures and 80% side overlap between flight lines. The delivered imagery measured $2048\times 2048$, with a mean GSD of 2.5 cm at the ground level. The sUAS utilized a checkerboard fight pattern and image sample rates of 1 Hz to generate a high overlap of stereoimagery for the generation of high-density lidar point clouds. The manned aircraft used a more traditional flight strip pattern of collection, optimized for the lidar sensor. A simple depiction of these two flight paths over a circular target can be seen in Fig. 4.

Fig. 4

(a) Checkerboard flight pattern versus (b) traditional flight strips over a circular forestry plot outline.

Agisoft PhotoScan Professional v1.2.6 was used to implement the SfM algorithms for these two imagery data sets in order to generate 3-D point clouds from stereo imagery.⁴³ In depth details of the processing steps of PhotoScan can be found in the literature,^44–46 but briefly: PhotoScan initially detects key point feature matches between the images using a process based on the scale invariant feature transform;^44,45 3-D positions of the key points and camera positions are then estimated through an iterative bundle adjustment; and this information then is used to generate a dense point cloud from the imagery.³³ In our scenes, ground control points (GCP) were not available for georeferencing the generated point cloud coordinates, thus the direct georeferencing was completed using the GPS information from the APX-15 at the time of each image capture. This is the same method used by Wallace et al. and Turner et al.^33,45 An advantage of using this method for georeferencing was that the SfM closely matched the accuracy of the lidar point clouds, since both data sets drew their georeferencing processing from the same source. Throughout the processing workflow of PhotoScan, high accuracy settings were used for the alignment of photos, alignment of camera positions, and the generation of the dense point clouds. The final point clouds were exported to the American Society for Photogrammetry and Remote Sensing (ASPRS) “.las” file format for ingestion of lidar-specific software analysis.

2.5.

Point Cloud Preprocessing

Each data set had to be processed through a ground classification algorithm for the generation of a DEM. These preprocessing steps were completed using the software suite LAStools (lasground_new function) with a step size of 6 m for the ground classification.⁴⁷ Point clouds were subset to $60\text{-}\times 60\text{-}\mathrm{m}$ grid cells, centered on the field plots. The large plot buffers avoid edge artifacts on the DEM produced, as the ground point classification routine typically contained false-positives along the edges of the point cloud. Circular buffers, matching the fixed radius of the reference plots of 37.5 ft. (11.43 m), were used for the final reduction of the point clouds (QGIS Version 2.18.0).⁴⁸ The clipping of the point clouds with the shape files were completed using the FUSION software suite (polyclipdata).⁴⁹

The height (DEM) normalization of point clouds was completed using LAStools (lasground_new/lasheight). All DEMs and products utilizing the DEMs, e.g., CHMs, were generated on the same 0.5 m raster grid and utilized the same parameters in the software workflow. The SfM for both data sets failed to generate enough ground returns to generate any sort of reliable terrain map, so a hybrid approach using the DEM generated from the corresponding lidar data sets was required for normalizing the SfM point clouds; this will be discussed later as a distinct disadvantage of the more cost-effective SfM approach. This is valid as both the lidar and SfM imagery were captured simultaneously and were georeferenced using the same GPS/IMU data.

2.6.

Functional Carbon Model

A primary metric of concern for sustainable forest management is the measure of carbon stored in the forest and any changes to that sequestered carbon over time. It therefore is imperative to model these data, as it cannot be measured directly from the remote sensing data. It has been shown that the carbon and AGB relationship can be approximated as a linear function, whereby carbon is approximated as 50% of the biomass value.^50–53 This method is common as biomass and carbon can only be directly measured through destructive methods.⁵² Zhao et al. introduced a linear, scale- and shape-invariant, temporally stable model for relating the CHD to biomass density. It is mathematically derived from the allometric relationships of canopy height and DBH, used by ground surveyors to generate biomass estimates from field data.^35,40 Our carbon model is derived from this same framework, since it is a linear model.

The model references a data product referred to the authors as the CHD, denoted as $p(h)$. This can be thought of as the percentage of ground, within the plot being measured, that is obscured by the top canopy at a given height. Effectively, the CHD is a histogram of the CHM height values measured at discrete height bins. Since the model is a function of a function itself, it is referred to as a functional model, i.e., the predictor $p(h)$ is a function or curve and not a group of predictors.^40,54 The relationship between the CHD and biomass, or carbon density, can be described by Eq. (1), with the full derivation from allometric assumptions summarized by Zhao et al.³⁵

A challenge with validating biomass and carbon models from other studies is that each data set covers a different distribution of biomass and carbon and potentially different forest types. Zolkos et al. addressed this issue in their review of 70 published AGB studies and found that expressing the residual standard error (RSE) of a model relative to the mean AGB or carbon value provides a satisfactory means of comparing studies between different regions with different forest CHM (or height) distributions.⁵⁶ The RSE is defined in Eq. (2).

In their review of these studies, it was noted that models assessed on temperate coniferous forest exhibited the lowest error rates [mean RSE(%) = 28.7%], when compared to temperate-deciduous forest [mean RSE(%) = 31%].⁵⁶ In the study by Means et al.,²⁴ it furthermore was found that heteroscedasticity existed in their models of basal area and volume, which also affect carbon and biomass models; they therefore modeled the natural logarithm of these dependent variables.²⁶

In addition to the RSE(%) metric, the RMSE and $r^2$ metrics are reported for each model and are defined in Eqs. (4) and (5), respectively, where $n$, $x_{ip}$, and $x_{\text{if}}$ are defined as they were in Eq. (2).

All CHMs generated in this study used a 0.5-m raster grid for sampling and assigned the maximum height return as the value for each cell. Each modality produced dense enough point clouds that no interpolation was required for empty cells. A depiction of the DEM generation, CHM generation, and CHD generation from a point cloud product can be seen in Fig. 5.

Fig. 5

Depiction of point cloud processing steps. (a) Raw point cloud with elevation values referencing mean sea level, (b) DEM generated from ground-classified points, (c) raw point cloud plotted above DEM for visual reference, (d) normalized point cloud after DEM subtraction with normalized DEM for visual reference, (e) CHM generated from normalized point cloud, and (f) CHD generated from the CHM.

2.7.

Plot Disturbance Detection

An ability to accurately detect forest harvesting or disturbance would facilitate compliance with sustainable forest management principles. Detection of individual harvested trees has been shown in boreal forests with a success rate of 73.4%, utilizing simple differencing of two CHMs taken 2 years apart.⁵⁷ Individual tree detection in closed-canopy, heavily overlapping deciduous forests, such as the sites in this study, is difficult and generally inaccurate, which led to this analysis at the plot level and via single date assessment. This is effectively a binary classification problem, i.e., the forest is disturbed or not, based on ideally simple structural metrics. These structural metrics were derived from the normalized point cloud and from point density values at discrete heights throughout the canopy. The point density values are different from the CHD in that it considers all point values to derive structural information, rather than just the height values of the top canopy. The density values were calculated as the percentage of point returns that fell within each voxel, with each voxel being defined as a $0.5\mathrm{m}\times 0.5\mathrm{m}\times 0.5\mathrm{m}$ cubed volume. Density values are then aggregated in the $x$ and $y$ directions, to generate a vertical structure profile for the plot as a whole, i.e., the percentage of returns within a height bin. These vertical profiles are analogous to the return type seen in a waveform lidar system and have been used to classify forest change temporally.³⁵ Considering this, a few metrics were chosen as predication variables:

For the binary classification operation, a linear support vector machine was trained using the fticsvm function in MATLAB (ver. R2016b).⁵⁵ To validate the classifier, two methods were chosen for comparison: holdout cross validation, where 50% of the observations were used for training and 50% used for validation, and leave-one-out cross validation, where all but one observation was used for training and the remaining observations used for testing.⁵⁹ One thousand trials were conducted, and the mean value of the producer’s accuracy, user’s accuracy, and the Kappa statistic were reported. These two methods of validation were chosen as there were only seven harvested plots in the 49-plot data set; thus, fold numbers were limited to seven for leave-one-out cross validation. It was determined that a comparison of the two validation methods would be the most informative of the performance of the classifier in differentiating between disturbed and undisturbed plots, given the small sample size of harvested plots. Although the sample size arguably is of concern for such a complex classifier, the cross-validation approach used here should at least provide a robust indication of performance.

2.8.

Data Usage

Not all data products were used for both the carbon models and harvest detection, due to a hardware failure and occlusion issues that will be discussed in further details in Sec. 3. Table 2 provides a concise description of data source-analysis pairs, as a reference for readers.

Table 2

Remote sensing data usage for each effort.

3.1.

Canopy Height Models Generation

It was imperative that the generated CHMs are accurate and precise measurements of the forest structure, since the functional carbon model is derived directly from the CHM. However, coordinate alignment of the point clouds was one challenge encountered with generating these models from different remote sensing platforms. We overcame this issue by processing point clouds in local coordinate systems and converting associated shape files used for clipping to match the local coordinates of each data source. This allowed for more accurate clipping of the point clouds for CHM generation, versus converting point clouds to a standard coordinate system, such as UTM 17N. Final alignments showed less than 1 m in alignment disagreement between the two data sources, i.e., manned versus sUAS platforms. This issue, however, does highlight the advantage of including GCPs in the target scene whenever possible, as there would be reference points to correct for these types of errors. A depiction of example CHMs for each data source can be seen in Fig. 6.

The SfM performance of the manned aircraft collection also was deemed unsatisfactory. Occlusion effects, due to the wide flight lines used to collect the lidar data, precluded many plots from generating full coverage of the study site, thus rendering it unusable for any type of modeling of the forest structure or carbon metrics. An example plot depicting the results of the occlusion affects is shown in Fig. 6(b). The sUAS SfM performance was visibly and quantitatively superior, when compared with the manned SfM, as seen in the CHM generated in Fig. 6(d). The sUAS SfM–generated point clouds covered the study area with a high density (1648.14 to $8826.46\mathrm{pts}/{\mathrm{m}}^2$, mean density of $3971.71\mathrm{pts}/{\mathrm{m}}^2$), and had a mean CHM difference of 1.96 m, when compared to the CHM generated by the sUAS lidar across all measured plots. This improvement of performance over the manned SfM data was attributed to the checkboard flight path used to generate imagery from multiple view angles, thereby avoiding occlusion effects. However, neither SfM data set produced adequate terrain measurements and required normalization by the DEMs generated from the lidar point clouds, for their associated sensor platform.

Fig. 6

Sample CHMs derived from the different sensor modalities of plot 113; this plot is one of the disturbed plots with high variation in the canopy structure. (a) CHM derived from lidar data set generated by the manned aircraft system, (b) CHM derived from SfM data set generated by the manned aircraft system, (c) CHM derived from lidar data set generated by the sUAS, and (d) CHM derived from SfM data set generated by the sUAS.

The majority of the CHM error, when comparing the sUAS lidar and sUAS SfM, was attributed to CHM misalignment and the boundary “sharpness” with which the SfM methods were capable of generating 3-D data. This can be seen in Figs. 6(c) and 6(d); plot 113 contains large magnitude, high frequency changes in canopy height, and the lidar-generated CHM has a blurring of these edges when compared to the SfM-generated CHM. Additionally, with the exception of the manned SfM, the maximum height measured within each plot for each data set was within 1 m of the reference data. Recall that this is within the error bounds of the reference measurements, as described in Sec. 2.2. This indicates that the peaks of the tallest trees were captured by the lidar and the imagery, and that the lidar-generated DEMs accurately depicted the topography below the canopy. This metric was used since individual crown delineation was not in the scope of this study, nor were the individual crown locations cataloged during the ground reference collection.

3.2.

Carbon Models

The manned lidar was the only data set able to cover enough sample plots to generate a truly independent model. The sUAS covered two different study sites, and the differences in species were enough to reduce performance to unusable levels [$r^2=0.12,\mathrm{RMSE}=21.72\mathrm{Mg}\mathrm{C}/\mathrm{ha}$, $\mathrm{RSE}(\%)=49.73\%$]. It was noted by the original authors of this modeling method that the $k(h)$ function in Eq. (1) is affected by tree allometry and crown geometries, thereby requiring that models be retrained (calibrated) for each new forest type.³⁵ For example, the increase in conifer species at the Om Sanctuary Site was enough to reduce model performance. When segregating model application to each study site, separately, performance improved distinctly [Clinch Mountain Only: $r^2=0.70,\mathrm{RMSE}=10.02\mathrm{Mg}\mathrm{C}/\mathrm{ha}$, $\mathrm{RSE}(\%)=19.35\%$; Om Sanctuary Only: $r^2=0.82,\mathrm{RMSE}=3.01\mathrm{Mg}\mathrm{C}/\mathrm{ha}$, $\mathrm{RSE}(\%)=11.33\%$]. However, sample sizes were small (9 to 13 samples, depending on sensor) for each site due to loss of data from a hardware subsystem failure during flight (lidar and timing packet corruption occurred).

The manned airborne system’s lidar data performed well when modeling the logarithm of the carbon data [$r^2=0.64,\mathrm{RMSE}=17.14\mathrm{Mg}\mathrm{C}/\mathrm{ha}$, $\mathrm{RSE}(\%)=27.83\%$], falling just short of the performance of the model when applied to a coniferous forest by the original authors ($r^2=0.72,\mathrm{RMSE}=10.75\mathrm{Mg}\mathrm{C}/\mathrm{ha}$).⁴⁰ These results corroborated results from other efforts, with an average of $r^2=0.78$ for the studies reviewed by Brovkina et al.,⁶⁰ noting these focused on more uniform coniferous forest types. This is a positive result, as it is common for biomass and carbon models to exhibit worse performance in deciduous forests, achieving lower $R^2$ values and larger RMSE values than for coniferous sites.⁶¹ Additionally, the RSE(%) attained is below the average RSE(%) reported by Zolkos et al.⁵⁶ for models applied to this forest type. In this case, the sUAS produced better results; however, the smaller sample sizes arguably did result in reduced rigor in terms of sUAS performance.

An alternative method to further compare the manned versus unmanned performance, therefore, was developed. This method applied the sUAS-generated CHD functions to the manned aircraft-trained model at the four sites, which overlapped between the two data sources. The predictions for these four plots closely agreed between the two platforms, with $r^2=0.99,\mathrm{RMSE}=1.77$ Mg C/ha. This demonstrates that both systems are capable of producing similar predictions when imaging the same study sites, and that the sUAS, from a data quality and utility perspective, is a viable alternative to employing a manned system.

The sUAS SfM data produced similar results to the sUAS lidar and underwent a similar testing procedure in terms of applying the SfM CHD to the trained, manned flight’s carbon model. When investigating the four overlapping plots, similar performance again was observed when compared to those for the lidar comparison ($r^2=0.98,\mathrm{RMSE}=2.13\mathrm{Mg}\mathrm{C}/\mathrm{ha}$). This indicated that a hybrid approach, i.e., lidar-derived DEMs and SfM-generated CHDs, is a viable option for this type of modeling. Table 3 provides a summary of carbon model performance on the various data sets utilized in this study.

These results suggest an opportunity for significant cost savings, since procuring and operating an RGB camera is significantly more affordable than lidar sensor packages. Considering that the underlying terrain typically does not change these areas, an approach could be to develop a base DEM from a lidar collection, then to use SfM at later dates for canopy structure measurements, even for temporal (change detection) analyses (Table 4).

Table 3

Carbon model results. All models were fit to the ln(Cf), where Cf is the field-measured carbon data by the ground crew during the audit.

3.3.

Detecting Plot Disturbance

The sUAS lidar did not cover enough examples of disturbed forest to be used for training and testing of a classifier; thus, the manned lidar system was used exclusively for this purpose. However, we would recommend for future studies that sUAS data be used for a similar analysis in future studies. Voxel-based structural data have been used to discriminate between varying stand types over time;⁴⁰ we chose to explore the ability of such data to discriminate between potential harvested forest plots at a single point in time.

We explored the use of binary classification to label plots as disturbed or undisturbed. Three structural metrics were chosen for use as feature descriptors: height of maximum energy, $h_{\max }$, an indicator of dominant canopy height, HOME, and the gap percentage of each observation. The mean vertical profile of each class, disturbed and undisturbed, along with the mean value for HOME and $h_{\max }$ for each class is shown in Fig. 7. Note the clear differences in the vertical profiles and the inverted relationship between HOME and $h_{\max }$ in the disturbed plots, when compared to the disturbed plots.

Fig. 7

(a) Mean vertical profiles for undisturbed observations and (b) disturbed observations. HOME and $h_{\max }$ for each are depicted by the blue line and black dashed line, respectively.

A linear SVM was utilized as the binary classifier for this study. An example of one of the many generated decision boundaries in this study can be seen in Fig. 8. The SVM will generate a different decision boundary for each set of inputs, requiring multiple repeated trials to gain insight into the average performance of the classifier. In this study, the data set is clearly linearly separable, with many of the training sets yielding similarly oriented hyperplanes. These properties allowed for the development of a high performing classifier.

Fig. 8

(a) A depiction of the feature descriptors and decision plane in a 3-D space. (b–d) 2-D depictions of feature descriptors and possible decision boundaries, for the HOME versus gap (b), maximum height versus HOME (c), and maximum height versus gap (d) variable sets.

As discussed before, the unbalanced nature of our data set, with seven samples for the disturbed class and 42 samples for the undisturbed class, necessitated two validation approaches:

Table 4

Results of the disturbance detector method. The detector performed best when applied to the aggressively harvested plots but performed poorly when applied to all disturbed plots.

We expanded the disturbance or nonsustainability detection to also include plots with moderate harvest conditions, 10 years prior to the collection date. These moderate disturbances proved to be difficult for the classifier to identify and were not clearly separable via the structural features we had selected. This can be seen in an example of the hyperplane generated by the detector while operating on the entire data set, depicted in Fig. 9, as well as noting the similarities in mean feature values between the undisturbed and moderately disturbed plots in Table 5. Any of the moderately disturbed observations are not separable utilizing these structural metrics, as they cluster within much of the undisturbed data points. This is due to the lack of aggressive harvesting, which in turn reduced the canopy variation, while also not introducing any large gaps in the canopy. It is worth noting that, in general, the gap feature did not provide as much insight into disturbance classification as HOME and $h_{\max }$. This is likely due to the mixed nature of the forest: when larger vegetation is removed, the understory remains and prevents the lidar first return from reaching the ground and registering as a gap location.

Fig. 9

Example decision plane generated by the linear SVM when applied to full extent of disturbed observations, including both severe and moderate harvest situations. The large number of support vectors is an indication of a poorly conditioned data set with high rates of confusion. The gray plane is the decision boundary generated by the linear SVM, and the points circled in green are the support vectors.

Table 5

Mean feature values for each class of observation.

We concluded, based on these results, that (i) given the set of disturbance indicators, the disturbance detector/classifier was capable of reliably detecting heavily disturbed plots and could provide utility in flagging locations for inspection and (ii) when the moderately harvested plots were included in the algorithm, the detector performance dropped significantly. This implies that the simple structural features derived, while they have a history of providing insight to other predictive models, were not adequate for all harvest instances and severities. Exploration of more intricate structural features or nonlinear classification methods may improve results. However, in the context of sustainable forest management, the detector was successful in detecting those plots with the most severe loss, as can be seen by the decision boundary shown in Fig. 9. Such a severe disturbance arguably is representative of top-cutting or even aggressive timber removal practices, both of which could indicate a reduction in forest sustainability.