Digital elevation models (DEMs) derived from airborne light detection and ranging (LiDAR) technology are becoming the standard in representing terrain surfaces due to its ability to accurately describe relief across large landscapes. LiDAR-derived DEMs have numerous applications in natural resources    . In forestry, DEMs are commonly used in hydrological applications (i.e.,   ) and increasingly in forest operations planning (i.e.,    ). The accuracy of LiDAR-derived DEMs is considerably higher than DEMs derived from alternative sources such as aerial photography or satellite imagery    , which facilitates the creation of high-resolution DEMs (≤1 m) and in turn increases their applications.
The accuracy in altimetry of high-resolution LiDAR-derived DEMs is commonly reported by data providers to be between 15 - 25 cm    . However, elevation errors are typically measured on flat, smooth terrain with no vegetation cover   . LiDAR system parameters such as flight path, scan angle and point density    as well as data processing procedures such as point filtering and interpolation methods     can have an effect on the derived DEM accuracy. Several studies have also demonstrated that vegetation cover, terrain slope, and terrain variability can have a significant effect on elevation errors      . In general, elevation errors increase with vegetation cover and amount of understory as well as with increasing terrain slope and variability    .
Most studies quantifying elevation errors from LiDAR-derived DEMs are based on static modes where errors are calculated using a permanent network of GPS stations that measure in continuous mode. Where not possible, GPS units with real time kinematics (RTK) technique are used to obtain accurate elevation of individual ground control points    . These GPS elevations readings are considered true values and compared with DEM-derived elevations. However, in forests with complex vegetation (dense, multistoried canopies formed by multiple species) and complex terrain conditions (steep, narrow, and short slopes with numerous rock formations), satellite signal strength is weak preventing the use of even high-end GPS units to obtain accurate elevations measures   . Consequently, there is a lack of studies evaluating the accuracy of LiDAR-derived DEMs in areas with complex vegetation and terrain conditions. Such is the case of forests in the Cumberland Plateau of the Appalachian Mountains, in the United States where there is an increasing interest in using LiDAR data for forestry, hydrology, mining, risk assessments, and wildlife applications   .
The objective of this study was to quantify elevation errors of high-resolution, LiDAR-derived DEM in areas with complex vegetation and terrain conditions using an alternative ground truth data collection method. Instead of using GPS units to collect elevations of control points, we used an alternative method that measured multiple relative elevation changes within field terrain plots. Elevation changes were measured between plot center and multiple points located at different distances and azimuths using a star-shaped plot design. This alternative ground truth data collection method is more appropriate to assess the ability of LiDAR-derived DEM to accurately represent terrain surfaces instead of comparing individual point elevations. Several DEM applications require computing elevation changes among different locations (hillshade, slope, flow accumulation, flow direction among others), hence a direct assessment of this relative elevation error provides useful information regarding the error propagation to those applications when using LiDAR data. Lastly, this study overcomes problems of complex vegetation and terrain conditions for the use of GPS with RTK by developing and using this novel DEM sampling method for elevation accuracy assessments.
2.1. Study Area
Research was conducted at The University of Kentucky’s Robinson Forest (Lat. 37.4611, Long. −83.1555), located in the rugged eastern section of the Cumberland Plateau region of southeastern Kentucky in Breathitt, Perry and Knott counties (Figure 1). Due to access limitations, we restricted the study area to the Clemons Fork and Lewis Fork watersheds within Robinson Forest covering an area of almost 1800 ha. Terrain across the study area and the entire Robinson Forest is characterized by a branching drainage pattern, creating narrow ridges with sandstone and siltstone rock formations, curving valleys and benched slopes. The slopes are dissected with many intermittent streams  and are moderately steep ranging from 10 to over 100%, predominately northwest and southeast aspects, and elevation ranging from 252 to 503 meters above sea level. Vegetation is comprised of a diverse contiguous mixed mesophytic forest made
Figure 1. Topography of the study area (1797 ha) within Robinson Forest (4250 ha) located in Breathitt, Knott, and Perry counties in southeastern Kentucky (Lat. 37.4611, Long. −83.1555).
up of approximately 80 tree species with northern red oak (Quercusrubra), white oak (Quercus alba), yellow-poplar (Liriodendron tulipifera), American beech (Fagus grandifolia), eastern hemlock (Tsugacanadensis) and sugar maple (Acer saccharum) as dominant and codominant species, while understory species include eastern redbud (Cerciscanadensis), flowering dogwood (Cornusflorida), spicebush (Lindera benzoin), pawpaw (Asiminatriloba), umbrella magnolia (Magnolia tripetala), and bigleaf magnolia (Magnolia macrophylla)   . Average canopy cover across Robinson Forest is about 93% with small openings scattered throughout. Most areas exceed 97% canopy cover but recently harvested areas have an average cover as low as 63%. After being extensively logged in the 1920’s, most of Robinson Forest is considered second growth forest ranging from 80 - 100 years old, and is protected from commercial logging and mining activities that are typical land uses in the region  .
2.2. LiDAR Datasets
We used two LiDAR datasets covering the study area, collected with the same LiDAR system by the same vendor. One dataset was low density (~1.5 pt∙m−2) collected in the spring of 2013 during leaf-off season for the purpose of acquiring terrain information, as part of a state-wide elevation data acquiring program from the Kentucky Division of Geographic Information. The second dataset was high density (~40 pt∙m−2) acquired in the summer of 2013 during leaf-on season for the purpose of collecting detailed vegetation information by the University of Kentucky’ Department of Forestry. The parameters of the LiDAR system and flight for both datasets are presented in Table 1. The vendor processed both raw LiDAR datasets using the TerraScan software  to classify LiDAR points into ground and non-ground points. A third dataset was also created by combining both low-density and high-density points. For
Table 1. LiDAR data acquisition parameters used for both datasets collected over Robinson Forest.
each of the three LiDAR datasets (low-density, high-density, and combined), the “Create LAS Dataset” tool in ArcGIS (ArcMap, version 10.2) was used to create a LASer (LAS) dataset file. The LAS dataset was then filtered to include ground points only, and the “LAS dataset to Raster” tool in ArcMap was used to create a 1-meter resolution DEM. Four DEMs for each dataset were created considering the average (AVG), inverse distance weighted (IDW), minimum (MIN), and nearest neighbor (NN) interpolation methods. As a result, a total of 12 DEMs covering the study area were created; three LiDAR datasets and four interpolation methods.
2.3. Sampling Design
We expected LiDAR-derived DEM errors to vary with terrain steepness and ruggedness (variability). For the purpose of identifying areas with different levels of terrain steepness and ruggedness in which to locate terrain sample plots, we resampled a 1-m resolution slope raster to a coarser resolution of 36.6 m (120 ft.) using the average value. This 36.6-m resolution was selected to provide a more meaningful scale across the study area and to match the size of field plots used to collect terrain surface information. The coarser resolution slope raster was only used to create three slope classes (low, medium, and high) with relatively similar area across the study site (Table 2, Figure 2(a)). Terrain ruggedness was also calculated at the coarser resolution as the slope variability of the 1-m slope raster, as used by  to determine slope heterogeneity. Slope variability for each cell in the coarser resolution raster was defined as the range (max slope?min slope) of slope values (Figure 2(b)). The study area was then divided into three ruggedness classes (low, medium, and high) with similar surface areas in each class (Table 2). Slope and ruggedness raster layers were then overlaid to identify the nine combined slope/ruggedness classes, in which five field plots were randomly located resulting in a total of 45 field plots (Figure 3).
Figure 2. Coarser 36.6-m resolution raster layers showing percent slope (a) and variability of percent slope, terrain ruggedness (b) across the study area used to create the nine combinations of slope/ruggedness categories in which to locate terrain sample plots.
Figure 3. Location of field plots within the study area. First letter in the abbreviated plot categories indicates level of ruggedness (second letter) and third letter indicates level of slope (fourth letter).
Table 2. Area (ha) under each combination of slope and ruggedness category considered to randomly select field plots to collect surface terrain information.
2.4. Field Plot Data Collection
Prior to the acquisition of LiDAR data, 1.2-meter square plywood boards, painted white were installed at the 45 plot locations determined using a hand-held GPS unit of a 6.1 m accuracy. Boards were installed leveled with their centers placed at the GPS determined location. After LiDAR acquisition, 3D coordinates and intensity values of LiDAR ground points were visually inspected and boards were clearly identified on 32 plots. The exact location of the remaining 13 field plots were found using field triangulation from easily identifiable features in both the LiDAR data and on the field. These features consisted of dominant trees, rock formations, road intersections, and other road features (i.e., cut and fill slope areas, ditch relief culverts, and bends). Azimuth and horizontal distance were measured using an electronic compass (MapStar Compass Module II, Laser Technology Inc.) and a sighted laser range finder (Impulse 200 LR, Laser Technology Inc.) mounted on a tripod.
Once the location of a given plot was identified, eight transects extending 18.3 m (60 ft.) from plot center were established covering 360˚ at 45˚ intervals (four cardinal and four ordinal directions). At six points along each transect, every 3.05 m (10 ft.) horizontal distance from plot center, we measured elevation change (vertical distance) between the plot center and the points (Figure 4). Elevation change for the 48 points (8 transect × 6 points) was measured by mounting the electronic compass and laser range finder on a tripod at the plot center and moving a target, mounted on a pole at the same height above the ground as the range finder, to each point. The location of a given point along a transect was determined by measuring the azimuth using the electronic compass with a precision of 1˚ and the horizontal distance using the sighted laser range finder within 0.06 m.
2.5. Data Analysis
The x- and y-coordinates of the 2,160 points (45 plots × 48 points) were obtained based on the plot center coordinates, horizontal distance from plot center, and azimuth of the corresponding transect. For a given point, the LiDAR-derived DEM elevation change was obtained from the elevation values of the DEM cells containing the location of the point and the plot center. Elevation changes from the field plots were considered true values. Errors in elevation change were calculated as the absolute difference between field plot elevation changes and DEM-derived elevation changes. The total 2160 elevation change errors were arranged by plot number and combination of slope and ruggedness
Figure 4. Diagram of the eight transects and six locations along transects used to collect elevation change information on each field plot.
class. Mean elevation change errors (MECE) were compared by slope and ruggedness class. The calculation of MECE was repeated for the twelve different DEMs to examine the effect of the LiDAR dataset and interpolation method. We tested for significant differences in MECE among slope and ruggedness categories using a 2-way analysis of variance (ANOVA) test in SAS 9.3 (Statistical Analysis Software, SAS Institute, Cary, NC, USA).
3. Results and Discussion
The overall mean error for all DEM-derived elevation changes from all three LiDAR datasets and four interpolation methods was about 73.6 cm (Table 3). Differences in MECE among LiDAR datasets were very small (within 1.5 cm) with the high-density and combined datasets provided the lowest (72.7 cm) and highest (74.1 cm) MECE, respectively. No significant differences were found among LiDAR datasets. When comparing interpolation methods, MECE values were even more similar ranging from 73.4 to 73.8 cm (Figure 5).
As expected, MECE increased with slope and ruggedness level (Figure 6). Values ranged between 42.5 cm and 100.4 cm from low to high slope classes, and between 53.2 cm and 96.0 cm from low to high ruggedness classes. Although MECE values among slope and ruggedness classes were significantly different (non-overlapping 95% confidence intervals), slope seemed to have a larger effect as evidenced by the larger variation. Even more variability in MECE values was observed when considering individual combinations of slope and ruggedness, ranging from 23.2 cm for the low slope/low ruggedness class to 145.5 cm for the high slope/high ruggedness class (Figure 7). Within the high ruggedness class, there is a clear increase in MECE with increasing slope, all of which were signifi-
Table 3. Area (ha) under each combination of slope and ruggedness category considered to randomly select field plots to collect surface terrain information.
Figure 5. Mean elevation error averaged by LiDAR dataset and interpolation method with their respective 95% confidence intervals.
Figure 6. Mean elevation error by slope and ruggedness classes averaged from all LiDAR datasets and interpolation methods with 95% confidence intervals
Figure 7. Mean elevation error by combination of slope and ruggedness classes averaged from all LiDAR datasets and interpolation methods with 95% confidence intervals.
cantly different among classes. Although in the other two ruggedness classes the MECE generally increased with slope, there was a less clear pattern. Within the medium ruggedness class, the high slope class had the highest MECE value which was significantly different than that of the medium and low slope classes. Within the low ruggedness class, the low slope class had the lowest MECE, which was also significantly different than the other two slope classes. A two- way ANOVA that combined MECE from all LiDAR datasets and interpolation methods corroborated these results (Table 4). Differences of the mean MECE values among slope classes (p < 0.0001) and ruggedness classes (p < 0.0001) were all significantly different. We also found a significant interaction between slope and ruggedness (p < 0.0001). In addition, slope classes had the highest F-value of two effects tested, which is further evidence of the greater impact of slope on MECE (Table 4).
The average error in elevation change from the surface terrain represented by the LiDAR-derived DEMs was about 74 cm, which is within values of elevation errors reported in the literature obtained from comparing point elevations using GPS units and LiDAR data (0.1 - 2.7 m)  but higher than errors provided by several studies    . Results showed that both slope and ruggedness had a significant effect on the MECE. Similar relationships have been reported by other studies. For example,  found that elevation error in steeper terrain was about twice as large as those in flatter terrain. Also  found even larger differences, three times as large, between elevation errors in steep and flat slopes. Similar patterns of larger elevation errors have been reported in areas with larger terrain ruggedness   . A likely reason is the misclassification of LiDAR points into ground and non-ground (vegetation) points. Classification algorithms assume that the lowest elevation LiDAR point in a given window (1 × 1 m) represents a ground point and that the slope between adjacent ground points is lower than slope between ground points and adjacent non-grounds point  . In steep and rugged terrain, these assumptions typically result in ground points with higher elevations than the lowest elevation point within the given window being misclassified as low vegetation  . This misclassification of LiDAR points are likely numerous in areas with small rock formations (i.e., outcrops and cliffs) and slopes steep enough that tree crowns are in close proximity to the ground, such as that of our study area. Although, we did not formally test for point misclassification, visual inspection of rugged area showed such cases, which likely resulted in large elevation errors.
Although there were large differences in point densities between the low- and high-density datasets, there were no statistical differences in the MECE. This is
Table 4. Two-way ANOVA using slope and ruggedness.
likely because the high-density data was collected during leaf-on season and only a small proportion of points were able to penetrate the dense canopy and reach the ground level. Figure 8 illustrates the spatial distribution of ground points of the three LiDAR datasets where only a small amount of additional ground points can be observed in the high-density dataset. The similar DEM accuracy level of low- and high-density datasets indicates the importance of the season where data is collected, and thus the importance of vegetation structure on ground elevation error. It seems preferable to acquire LiDAR data during leaf-off season to accurately represent terrain surfaces using LiDAR data. On the contrary, if one is interested in retrieving vegetation attributes, then collecting LiDAR information during leaf-on season seems more appropriate. Several studies have reported similar findings   . Our results also showed no significant differences among interpolation methods with mean elevation error within 1 cm for all three LiDAR datasets. This is likely because of the nominal pulse spacing of both datasets was lower than the 1 m resolution of the LiDAR-derived DEMs. Other studies evaluating different interpolation methods have found significant differences for larger resolution (>5 m) and similar result for fine resolutions   .
Figure 8. Shown is the LiDAR ground point pattern and distribution across a landscape for the high-density (a) low-density (b) and combined (c) LiDAR datasets.
The increasing interest in using LiDAR data for forestry and natural resources applications creates a need to evaluate the accuracy of LiDAR-derived DEM. There are numerous studies quantifying elevation errors and the effect of terrain slope, ruggedness, and vegetation cover. However, there is a lack of studies measuring elevation error in areas with complex terrain and vegetation conditions due to weak GPS satellite signal strength. We used an alternative method to assess DEM accuracy based on measuring elevation changes within 0.1-ha terrain plots. As several DEM applications require the calculation of elevation changes among different locations (slope and flow accumulation among others), there is a need to directly assess relative elevation change errors instead of absolute elevation error from individual control points. Our approach is then more appropriate to assess the ability of LiDAR-derived DEM to accurately represent terrain surfaces. This method is applicable to many areas with complex vegetation and terrain conditions, such as those of the Cumberland Plateau of the Appalachian Mountains as well as other mountainous areas.
Terrain slope and ruggedness had a significant effect on the accuracy of LiDAR-derived DEMs, but slope had a larger effect. MECE ranged from 42.5 cm in areas with flat terrain to 100.4 cm in areas with steep terrain. MECE values ranged between 53.2 cm and 96.0 cm in areas with low to high terrain ruggedness. A larger MECE variation was observed when evaluating the combined slope/ruggedness classes, with values ranging from 23.2 cm in areas with lowest slope and ruggedness classes to 145.5 cm in areas with highest slope and ruggedness classes. Despite large difference in point density among LiDAR dataset (low-, high-density, and combined), errors were similar mainly because of the high-density dataset was collected during leaf-on season and only a small portion of points reached the ground. The small amount of additional ground point in the high-density dataset was not enough to have an effect on the DEM accuracy. As the DEM analyzed had a relatively high resolution of 1 m and above the nominal pulse spacing, interpolation methods also did not have an effect on the MECE. Larger resolution should be included in future research to evaluate their effect, especially in steep and rugged areas under dense vegetation cover. Although LiDAR point misclassification was not formally tested, visual inspection of terrain plots in the rugged terrain class showed clear evidence of such cases, which indicates the need to evaluate alternative classification algorithms to determine the most appropriate for areas with complex terrain and vegetation conditions.
This work was supported by the Department of Forestry at the University of Kentucky and the McIntire-Stennis project KY009026 Accession1001477.