Find the information such as human life, natural resource,agriculture,forestry, biotechnology, biodiversity, wood and non-wood materials.
Blog List
Thursday, 22 December 2016
Mapping selective logging impacts in Borneo with GPS and airborne lidar
Published Date Forest Ecology and Management 1 April 2016, Vol.365:184–196, doi:10.1016/j.foreco.2016.01.020 Open Access, Creative Commons license
Author
Peter Ellis a,,
Bronson Griscom a,
Wayne Walker b,
Fabio Gonçalves c,
Tina Cormierb,
aThe Nature Conservancy, 4245 Fairfax Drive, Arlington, VA 22203, United States
bWoods Hole Research Center, 149 Woods Hole Rd, Falmouth, MA 02540, United States
cAgrosatélite Geotecnologia Aplicada, 4850 Rod SC-401, Florianópolis, SC 88032-005, Brazil
Received 20 November 2015. Revised 15 January 2016. Accepted 16 January 2016. Available online 3 February 2016.
Highlights
•
Lidar and GPS can be used to monitor logging performance.
•
Lidar maps of logging infrastructure coincide with ground-based GPS maps.
•
Lidar maps of logging disturbance reveal extensive unlogged areas.
•
Fourier methods map biomass differences within 5% of field measurements.
•
GPS maps of logging infrastructure are more affordable for near term monitoring.
Abstract Reduced-impact logging (RIL) is a promising management strategy for biodiversity conservation and carbon sequestration, but incentive mechanisms are hindered by inadequate monitoring methods. We mapped 937 ha of logging infrastructure in a selectively harvested tropical forest to inform a scalable approach to measuring the impacts of discrete management practices (hauling, skidding, and felling). We used a lidar-derived disturbance model to map all skid trails and haul roads within 26 months of the selective harvest of six blocks of dipterocarp forest in five industrial concessions in East Kalimantan, Indonesia. Lidar maps of logging impacts (220 ha) agreed well with ground-based maps (total of 217 ha, RMS error of 6 ha or 3%), but skid trail positions agreed only 59% of the time. Due to rapid forest regeneration, total lidar-derived haul road area was 31% smaller than road area measured in the field; agreement was higher for lidar collections within a year of the harvest. Maps of carbon density generated from Fourier transforms of lidar height profiles estimated skidding and felling biomass losses to within 1–5% of ground-based measurements. Lidar-derived skidding and hauling impact zones covered only 69% of the permitted harvest area; the remaining areas showed no signs of logging disturbance, and available biophysical data did not explain their location. These results emphasize the need for more extensive mapping of logging infrastructure to capture spatial variability in skid trail density and hitherto undetected no-impact zones. While a ground-based GPS is recommended as the most affordable method for wide-scale infrastructure mapping, aerial lidar is an effective tool for remotely quantifying the extent of logging impacts in tropical forests. Keywords
Despite the demonstrated benefits of RIL, opportunities to link ecosystem service incentives with RIL have been limited by the absence of robust verification systems (Putz et al., 2012). Unlike clear-cut harvests and deforestation, selective logging impacts are notoriously difficult to detect and quantify with available satellite imagery (Read, 2003, Frolking et al., 2009, Weishampel et al., 2012 and Réjou-Méchain et al., 2015). To date, performance measurement systems for selective logging have relied on sample-based measurements of infrastructure (roads, log landings, and skid trails) collected from a combination of optical image interpretation, Global Positioning System (GPS) maps, and field inventory data (Griscom et al., 2014 and Pearson et al., 2014). Because these approaches are sample-based (not wall-to-wall), they fail to capture the true spatial distribution of logging infrastructure. And because they rely on two-dimensional overhead imagery, they miss the sub-canopy impacts of selective logging.
The active sensor, three-dimensional assets of airborne lidar (light detection and ranging) provide a solution to this problem. Extensive research has focused on the use of lidar to estimate above-ground biomass or carbon density (Zolkos et al., 2013), but no cost-effective lidar monitoring system currently exists for detecting and measuring selective logging impacts at the logging block scale. In a first step toward such a system, D’Oliveira et al., 2012 and Andersen et al., 2014 devised a novel approach (the Relative Density Model, or RDM) to map haul roads, skid trails, and felling gaps using a dense lidar point cloud.
Here we provide a field-validated test of lidar-derived RDMs to map vegetation damage based on logging infrastructure and apply lidar-based mapping methodologies to quantify the impact of selective logging on the structure and aboveground carbon density of a dipterocarp forest in Borneo. Our results inform logging infrastructure impact parameters for the first applicable methodology to verify carbon emissions reductions from RIL (VCS, 2015), which could contribute 0.16 Gt C yr−1 of climate mitigation (Putz et al., 2008). Our objectives were to: (1) test the applicability of the RDM approach using commercial-grade, small-footprint lidar at medium pulse density (<5 pulses m−1); (2) use large, RDM-mapped skid trail networks to expand the sample size of the latest field-based RIL carbon performance measurement methodology (Griscom et al., 2014) and more accurately map the extent of no-harvest zones; and, (3) estimate differences in forest structure and aboveground carbon density within impact and no-impact zones mapped by the RDM.
2 Materials and methods
2.1 Study site
For five logging concessions in East Kalimantan, Indonesia (B, C, D, G, and H in Griscom et al., 2014), we obtained geospatial coordinate locations of six cutting blocks (314–1750 ha) and 54 nested sub-blocks (42–209 ha, Table 1, Fig. 1) harvested in 2010 and 2011. All concessions were characterized by commercial timber-producing species of Dipterocarpaceae growing on latosol soils at elevations of 201–630 m with average slopes of 41% rise (see methods below). Mean annual rainfall was 225–425 cm yr−1. The most frequently encountered trees ⩾60 cm diameter at breast height (DBH; stem diameter at 1.4 m or above buttresses) were commercially valuable species of Shorea and Dipterocarpus (see Table 1 for timber volume extraction rates). Commercial logging in our study area occurs on 35 year rotations (Sist et al., 2003) and involves the construction of haul roads along ridge lines using Caterpillar D7 bulldozers to skid logs ⩾60 cm DBH (on slopes up to 80%) to log yards that are adjacent to haul roads. Concessions were selected from the nine concessions sampled in Griscom et al. (2014) based on proximity to the closest airport, and represent the same set of landscape conditions defined there.
Table 1. Cutting block spatial statistics for reported sub-blocks. Parenthetical values are the number of sub-blocks in the reported harvest area. Ranges are 95% confidence intervals, calculated from sub-block units (not available for hauling due to low sample size). Skid trail density is calculated as length per skidding/felling zone.
Logging concession
Year cut
Sample area (ha)
Reported harvest area (ha)
Skid trail density (m ha−1)
Hauling extent (%)
Skidding/felling extent (%)
No-impact extent (%)
Average slope (% rise)
Harvest intensity (m3 ha−1)
B (5)
2010
417
374
164 ± 18
2.0
78 ± 8
20 ± 8
57 ± 3
28
B (13)
2011
1308
1061
156 ± 11
2.4
63 ± 10
35 ± 10
49 ± 3
30
C (4)
2011
418
320
189 ± 38
3.1
55 ± 16
42 ± 18
38 ± 13
31
D (4)
2010
378
314
137 ± 37
1.0
33 ± 9
66 ± 9
53 ± 4
32
G (20)
2011
1841
1750
249 ± 12
3.6
77 ± 7
20 ± 8
28 ± 6
43
H (8)
2011
1259
846
158 ± 24
3.7
54 ± 15
43 ± 16
43 ± 12
58
Mean
937±495
777±455
175±32
2.6±0.8
60±13
37±14
41±9
37±9
Fig. 1. Location of cutting blocks sampled in this study.
2.2 Lidar data
A commercial lidar company based in Jakarta, Indonesia (PT Credent Teknologi) was commissioned to acquire 5620 ha of medium-density, discrete-return lidar data across the study site in April 2013, within 26 months of closure of logging operations. Data could not be collected from the southern portion of concession C (1130 ha) due to cloud cover. We ranked cutting blocks on the basis of the “time lag” between harvest and lidar data collection.
Credent collected lidar data at an average pulse density of 4.5 m−2 (above the threshold recommended by Leitold et al., 2015), with 40% swath overlap, scanner angle <20°, and standard flying height of 650 m above ground level. Credent delivered LAS point-cloud files with ground points classified, digital terrain models (DTM) of ground surface (1 m resolution), digital surface models (1 m resolution), percent canopy cover (15 m resolution), canopy height models (1 m resolution), and contemporaneous color orthophotos (∼0.78 m resolution). Slope values (% rise) were calculated from 17 m resolution rasters derived from 1-m DTMs using ArcGIS spatial analyst.
2.3 Lidar disturbance model
To detect skidding and hauling disturbances, we generated 1-m resolution RDMs for the entire lidar sample area, using methods adapted from D’Oliveira et al., 2012. RDM rasters were created from LAS point-clouds using ArcGIS 10.2 model builder and the lasheight, las2las (filter), and lasgrid tools from the LAStools toolbox (http://lastools.org). The RDM model reveals spatial corridors in understory vegetation created by skidding operations. RDM values are calculated by dividing the number of lidar point returns in the RDM height stratum by the total number of returns below the upper RDM limit (Fig. 2 and Fig. 6, D’Oliveira et al., 2012). This analytical method was applicable to the logging landscape in Berau because bulldozers with ∼3.7 m wide blades are used to skid logs from stump to log yards at the margins of haul roads, thereby creating sub-canopy corridors of aboveground biomass disturbance. We used a default RDM lower limit of 1 m and an upper limit of 10 m, identified after testing alternative limit values and assessing the clarity of RDM skidding corridors with our field GPS map of skid trail locations. RDM results were also used to provide a more accurate delineation of disturbance from haul road construction by fully capturing the damage to understory vegetation normally obscured by overhanging tree canopies.
Fig. 2. Illustration of the RDM calculation method adapted from D’Oliveira et al. (2012) with permission of the author. The high RDM value illustrated here indicates a site with a relatively intact understory and therefore low likelihood of skidding impact. RDM values near zero indicate understory damage. An RDM with no data occurs when there are no returns below the lower RDM limit, indicating a very dense over-story that occludes pulse returns from ground level.
2.4 Skid trail and road delineation
Using the default RDM and a DTM for reference, the lead author on-screen digitized the centerlines of all skid trails and haul roads for the six cutting blocks (5620-ha total sample area) at a scale of 1:3500. Haul roads and skid trails were identified by the absence of RDM corridor vegetation (RDM value approaching 0, i.e. almost no lidar returns in the 1–10 m corridor), and visible disturbance of earth in the DTM resulting from bulldozer excavation and subsequent traffic. Haul roads were distinguished from skid trails by the presence of an RDM corridor >15 m wide (Griscom et al., 2014), more pronounced ground disturbance in the DTM, occasional wide areas in the RDM corridors representing log yards, and location on ridgelines. Coincident field GPS data was not consulted during the delineation process.
Although there are many automated linear feature extraction software tools available, (e.g. Clode et al., 2007, White et al., 2010 and Azizi et al., 2014), manual digitization was preferable for this study because: (1) human interpreters are particularly adept at synthesizing a wide range of spatial cues to group features into meaningful line patterns (Quackenbush, 2004); (2) the human interpreter in this study (lead author) was familiar with local landscape conditions and logging practices from extensive field work in the same concessions; and (3) false RDM signals from streams can be more easily ruled out by referring to the DTM.
2.5 Impact zone delineation
We define impact zones as geographic areas of logging disturbance to vegetation; road centerlines define the center of hauling impact zones and skid trail centerlines define the center of skidding/felling impact zones. Because vegetation is completely destroyed during road construction, the hauling impact zone is spatially coincident with the hauling corridor. Due to the diffuse nature of selective logging, the skidding/felling impact zone is more expansive, allowing for felling damage that occurs outside the skidding corridor. We pooled skidding and felling disturbance into a combined impact zone because these disturbances often occupy the same geographic space (e.g. harvest trees are felled onto skid trails). We mapped the extent of all impact zones within logged sub-blocks reported to the Indonesian Ministry of Forestry (reported harvest area).
To delineate hauling impact zones (Fig. 7), we selected all zero-value contiguous RDM pixels that intersected digitized road centerlines. We then grouped the selected pixels to generate vector polygons and re-grouped the individual polygons into a single continuous road polygon for each block by filling high RDM-value anomalous gaps (e.g. created by post-harvest regrowth). To eliminate skidding and felling zones from the final hauling impact zone polygon, we manually pruned segments that intersected digitized skid trails.
To map the skidding/felling zone (Fig. 7), we used the “area accessed” from Griscom et al. (2014) (48 m) to generate a planimetric buffer of digitized skid trails. This buffer distance was determined by adding the maximum measured distance from skid trails centerlines to stumps (N = 124, after removing outliers) to the maximum dipterocarp forest crown radius measured by Ashton et al. (1992). We removed overlapping portions of skidding/felling zones where they intersected hauling impact zones. We categorized the total area within the reported harvest area but outside of the skidding and hauling impact zones as no-impact zones in full recognition that edge effects (e.g., indirect impacts on understory humidity) may extend unknown distances into the adjacent forest.
We generated spatial statistics from skid trail/haul road features and lidar-derived rasters in each of the three zones (Table 1). All haul road, skid trail, and impact zone statistics were calculated based on a horizontal plane. Skid trail density was calculated as the length of skid trail per area skidding/felling zone (m ha−1, Fig. 3, VCS, 2015). To calculate intra-cutting block variability in skid trail density, no-impact zone area, slope, canopy cover, and canopy height, we intersected impact zones with sub-block boundaries in the reported harvest area, using sub-blocks as the unit of replication. Hauling extent was calculated as hauling impact zone area divided by cutting block area and expressed as a percentage. Because haul roads occupy much smaller portions of harvested sub-blocks, and often straddle sub-block boundaries (Fig. 3), we did not evaluate hauling extent at the sub-block scale.
Fig. 3. 1:12:000 map of two entire and portions of 4 adjacent sub-blocks from concession G presenting differences in skidding density, expressed as skid trail length (thin red line) per skidding/felling zone (blue area). Haul road length (thick red line) per hauling impact zone (orange area) is also shown. The skidding densities in the upper left and right sub-blocks shown (green boundaries) are 273 and 224 m ha−1 respectively. No impact zones are shown in white, and were not included in the skidding density calculation.
2.6 Carbon density mapping
To estimate field-based forest biomass of unlogged areas, we collected tree measurements in randomly selected pre-harvest cutting blocks adjacent to sampled logged blocks. At 100 m intervals on two parallel transects separated by 100 m we measured the DBH (stem diameter at 1.4 m) of all ‘in’ trees using a 9.18 m2 ha−1 BAF prism (40 ft2 acre−1), and then tallied the number of ‘in’ trees using a 4.59 m2 ha−1BAF prism (20 ft2 acre−1), following the dual prism sampling method described by Marshall et al. (2004). A total of 27 prism plots were sampled across concessions C, D, and H, and standard errors of prism-based biomass estimates were calculated using Marshall et al. (2004) Eq. (2). We also measured the height of the largest DBH tree per plot with a clinometer. More information on biomass plot field methods can be found in Griscom et al. (2014).
To map carbon density, the discrete-return lidar data was tiled at a resolution of 20 m × 20 m to be consistent in scale with the average size of the variable radius prism plots (i.e., 10 m radius). A waveform was synthesized for each tile as well as for a 10-m radius vertical column centered on each of the 27 prism plots. Waveform synthesis was accomplished by binning the lidar point returns into vertical bins 50 cm wide and counting the number of returns from each bin (Fig. 4). This procedure included removing the influence of ground slope on each lidar point-cloud return by subtracting from its elevation (i.e., z value) the bald-earth elevation estimate derived from the 1-m resolution DTM. A series of filters were also applied to flag and remove tiles with anomalous waveform characteristics (1–4% of the tiles, depending on the concession), including tiles with (1) no ground returns, (2) a point density <1 point/m2, (3) more than 1% of normalized heights >82 m (i.e., 20% greater than the maximum field-observed height), and (4) more than 1% of normalized heights <−3 m.
Fig. 4. Example of a synthesized lidar waveform generated from raw point-cloud data associated with a cylindrical vertical footprint of 10-m radius centered on a single field plot. Each point represents the number of returns (x axis) in each 50-cm height interval (y axis). The gray line provides a loess-smoothed functional representation of the data (for illustration purposes), which suggests a multistoried canopy structure with a limited understory (∼5 m) and mid-story (∼25 m) beneath a more pronounced multi-layer over-story (∼50–70 m). The tallest tree measured from the ground in this plot was 64 m and 305 MgC ha−1 was the field-based estimate of aboveground carbon density.
Following Treuhaft et al., 2010 and Treuhaft et al., 2013 and Gonçalves (2014), we took Fourier transforms of the synthetic lidar waveforms to serve as explanatory variables in carbon density models, exploring the vertical scales at which the canopy is organized explicitly as it relates to aboveground biomass. The normalized Fourier transform γ(f) at spatial frequency f was calculated as (Bracewell, 1986):
equation1
where w(z) is the number of lidar point returns included in the 50-cm bin centered at height z. The peak of the ground return in the synthetic waveform was assumed to be at z = 0, and w(z) was set to zero at altitudes higher than the tallest tree. The waveforms were decomposed into frequencies ranging from 0.3 to 0.01 cyc/m (vertical wavelengths of 3.3–100 m), with a step of 0.01 cyc/m, resulting in a total of 30 complex valued metrics, each described by an amplitude and phase. Regression analysis was used to model the relationship between field-estimated aboveground carbon density (MgC ha−1) and Fourier transforms extracted from the synthesized waveforms. We conducted a data-driven test of potential regression models by fitting all possible explanatory variable subsets with up to 6 metrics per subset and selecting the one that minimized the Bayesian Information Criterion (BIC). A multivariate linear regression including Fourier transforms at vertical wavelengths between 3.3 and 50 m—namely amplitudes at 0.02, 0.03 and 0.13 cyc/m, and phases at 0.11, 0.24 and 0.30 cyc/m—was selected after determining that the relationship between biomass and the structural metrics could be linearized by applying a logarithmic transformation to the response (i.e., aboveground biomass), as suggested by graphical analysis and the Box–Cox procedure (Box and Cox, 1964). This model explained 83% of the variation in aboveground biomass, with a RMS error of 31.4 MgC ha−1 (16.8% of the field estimated mean, Fig. 5a). A repeated random sub-sampling validation test based on 1000 iterations resulted in an average cross-validation RMS error of only 42.1 MgC ha−1 (22.5%), indicating a high predictive ability of the model (Fig. 5b). We note that all statistics were calculated after converting the predictions back to the original scale of measurement, accounting for back-transformation bias. Regression model assumptions (linearity, constant variance, normality, and independence) were evaluated using formal tests and graphical tools, and no significant deviations were identified. The selected model was used to map carbon density across the study site at a spatial resolution of 20 m.
Fig. 5. (a) Relationship between observed (field estimated) and predicted (modeled) aboveground biomass based on data from 28 variable-radius prism plots plus one dummy plot (Z), identified post hoc in the middle of a haul road with zero aboveground carbon density. Error bars represent 95% prediction intervals. (b) Distribution of the RMS error calculated using a Monte Carlo cross-validation approach (1000 iterations). The average cross-validation RMS error was 42.1 MgC ha−1 (22.5%) and the standard deviation was 19.3 MgC ha−1 (10.3%).
2.7 Comparison of lidar-based and GPS infrastructure maps
To evaluate the agreement of impact zones, we calculated regression statistics comparing lidar-based results to impact zones derived from GPS tracks mapped in the field using Garmin 60 CSX GPS receivers as described in Griscom et al. (2014).
To evaluate the positional agreement of skid trails and haul roads, we followed Goodchild and Hunter (1997) to calculate the percentage of lidar-digitized infrastructure lines (Section 2.4) that fall within a specified horizontal distance of the GPS line. We used a GPS line buffer (“specified horizontal distance”) of 10 m to account for the positional accuracy of the GPS receiver (Garmin Ltd., 2006).
To evaluate the agreement of the skidding/felling zones, we used the ArcGIS union tool to delineate an area of commission, where a lidar impact zone is unconfirmed by GPS reference data, and omission, where a lidar impact zone misses impacts observed in GPS reference data. We calculated percent agreement from a confusion matrix (Congalton, 1991). Note that we do not assume lidar or GPS data have inherently higher accuracy, but consider the amount of agreement or disagreement between them as an indicator of the accuracy of both approaches.
To evaluate the agreement of the hauling impact zone maps for each cutting block, we calculated the difference between hauling impact zone area and coincident haul road area from field measurements collected by Griscom et al. (2014). Note that haul road widths included log landings, which we include as part of the road corridor.
All statistical analyses were performed in the R statistical package (R Development Core Team, 2014). We tested for normality and equal variance of the residuals when necessary. Paired comparisons of means were conducted using the Bonferroni–Holm test. All reported ranges reflect 95% confidence intervals unless otherwise noted.
3 Results
3.1 Spatial distribution of logging impacts
The lidar RDM-mapped vegetation corridors 1–10 m above ground, thereby revealing the spatial distribution of all logging impacts, including skidding and felling damage typically occluded from optical sensors by over-hanging vegetation (Figs. 6 and S1). Using RDM as a reference, on-screen digitization resulted in the delineation of 71.4 km of roads and 657.2 km of skid trails in the 5620 ha sample area. Hauling and skidding/felling zones occupied, on average 3% and 60% of the harvest area in the sampled sub-blocks respectively, leaving 37% of the harvest area not directly affected by logging (no-impact zone, Table 1, Fig. 7).
Fig. 6. 1:10,000 1-m resolution map of RDM from concession B’s 2011 cutting block based on lidar data collected 19 months after harvest. RDM values (0–100) are shown in grayscale black to white, with green for no data (no aboveground pulse returns below 10 m). Purple lines delineate sub-block boundaries and the gray line bounds the lidar acquisition. Thick red lines are digitized haul roads, and thin red lines are digitized skid trails. Note the dark sinuous pattern of a natural stream course (blue arrow), vegetation disturbance from post-skidding induced hillside erosion (orange arrow), and patches of white at the edges of skid trails and roads indicating post-harvest regeneration following skidding, felling, and hauling impacts (RDM = 100, green arrow).
Fig. 7. 1:10,000 map of lidar-delineated skidding/felling impact zones (blue), and hauling impact zones (orange) in the concession B 2011 cutting block, with DEM hillshade in background for topographic reference. Note haul roads and major skid trails concentrated on ridgelines, extensive no-impact zones in the central sub-block, and un-harvested sub-block to the north.
For all sub-blocks, we detected no slope difference between skidding/felling zones and no-impact zones. When stratified by concession, cutting blocks B-2011 and D had steeper slopes in no-impact zones than skidding/felling zones, whereas cutting block G showed the opposite pattern (Table 2).
Table 2. Means of sub-block scale slope values calculated from lidar DTM at 17-m resolution, showing t-test results by cutting block. Cutting blocks with significantly different slope values between skidding/felling and no-impact zones are shown in bold.
Cutting block
Skidding/felling zone slope (%)
No-impact zone slope (%)
df
T
p
B-2010
55.0
59.6
4
−0.74
0.482
B-2011
38.4
50.2
12
−5.88
<0.001
C
38.9
33.0
3
0.69
0.539
D
38.9
53.3
3
−4.23
0.024
G
32.3
23.1
19
4.53
<0.001
H
41.2
43.4
7
−0.74
0.482
All cutting blocks
38.2
39.0
53
−0.46
0.649
3.2 Comparison of lidar-based and GPS infrastructure maps
Ninety-eight percent of the variance in GPS-defined impact zone area is explained by the lidar RDM, with a RMS error of 5.8 ha (total GPS impact zones in sample = 217 ha; lidar impact zones in same sample = 220 ha; Fig. 8)
Fig. 8. Relationship between GPS field sample and predicted lidar-based impact zone areas for the 6 sampled cutting blocks (labeled by letter). Dashes delineate the one-to-one line.
Positional agreement was lower than impact zone agreement, with 59% of the RDM skid trail length within 10 m of field reference GPS lines. Overall positional agreement of skidding/felling zones was 70%; the RDM skidding/felling zone was 14% larger than the reference GPS skidding/felling zone due to RDM commission almost twice omission (Fig. 9, Table 3).
Fig. 9. 1:7000 map of GPS sample area from concession G showing GPS skid trails (dotted blue lines), corresponding lidar skid trails (solid black lines), and lidar skid trails from adjacent areas un-sampled in the field (gray). Areas of skidding/felling zone agreement are in green, RDM omission in orange, and RDM commission in red.
Table 3. Results comparing skid trail and haul road lines and zones from the lidar RDM to field GPS lines and impact zones.
Estimated time between harvest and lidar collection, based on communication with concession managers. Ranges reflect different harvest completion dates of nested sub-blocks.
Calculated as field-based haul road area – lidar hauling impact zone area/field-based haul road area.
The lidar hauling impact zone was on average 31% smaller than the field-based measure of hauling road area. The difference between RDM haul road area and field-based haul road area increased as time lag increased (R2 = 0.86, Table 3). When cutting blocks harvested >2 years prior to the April 2013 collection date are removed from the consideration, all positional measures of agreement improve.
3.3 Skid trail density and haul road extent
Using sub-blocks as the unit of replication, a paired comparison of means indicated that only cutting block G’s skid trail density differed from other cutting blocks (p < 0.004, Table 1); its values were highest for both lidar and GPS-based estimates. When cutting block G was excluded, skid trail density remained within 17% of a mean value of 161 m ha−1. Mean lidar skid trail density was higher than the mean from GPS subsamples (175 and 141 m ha−1 respectively, p = 0.047).
Again using sub-blocks as unit of replication, we detected a weak negative relationship between terrain slope and RDM skid trail density (R2 = 0.097, p = 0.02), but this relationship was not evident at the cutting block scale and disappeared entirely if concession G was removed from the analysis (R2 = 0.014, p = 0.50). Skid trail density measured by lidar was not related to harvest intensity (p = 0.74) or committed carbon emissions as measured by Griscom et al. (2014, R2 = 0.312, p = 0.25).
Hauling extent varied considerably among cutting blocks (Table 1), but where the time lag was <2 years, mean hauling extent increased and variability decreased (from 2.6 ± 0.8% to 3.5% ± 0.3%). We detected a weak negative relationship between hauling extent and slope (R2 = 0.55, p = 0.09), and no relationship between haul road length and percent skidding/felling zone or skid trail density (p = 0.51 and 0.70 respectively).
3.4 Forest structure and carbon density
Compared to un-harvested areas in the no-impact zones, lidar-based estimates of canopy height and percent canopy cover in the skidding/felling zones were lower by 4.6 m (15.1%) and 5.3% (5.4%) respectively (p = 0.001 and 0.002). Lidar-based average aboveground live dry carbon density in the no-impact zone was 203 ± 7 MgC ha−1 – within the 95% confidence interval of pre-harvest field plots (187 ± 28 MgC ha−1). Mean carbon density in the skidding/felling zone of the lidar-based map was 178 ± 20 MgC ha−1, which is 13 ± 8% lower. When blocks D and B-2010 are removed to eliminate time lags >2 years, the lidar carbon density difference is 18 ± 8%. We found no difference in carbon density between no-impact zones and unreported zones (Fig. 10 and Fig. 11).
Fig. 10. Distributions of mapped carbon density (MgC ha−1) for each cutting block compared to the carbon density distribution observed in the field.
Fig. 11. 1:10,000 map of carbon density values in concession B (same region shown in Fig. 6). Note the clear signal of the haul road (dark red, black arrow) and skidding/felling impacts (blue circle, see Fig. 7), which reflect a 12% reduction in carbon density.
4 Discussion
4.1 Spatial distribution of logging impacts
In the reported harvest area of the six cutting blocks we studied in East Kalimantan, 3202 ha of the 4665 ha (69%) were directly affected by timber harvesting. Previous studies have insufficiently sampled logging impacts to accurately extrapolate area based emissions factors (Pearson et al., 2014, based on data from six sites across Indonesia, Republic of Congo, and Latin America; and Griscom et al., 2014, in the same study region). In order to use reported harvest areas to extrapolate logging emissions estimates to entire regions using field measurements collected in impact zones (as recommended by these studies and the RIL-C VCS methodology, VCS, 2015), area-based emissions factors need to be adjusted accordingly. Based on a limited field sample of “area accessed” within the same harvest blocks we report on here, Griscom et al. (2014) assumed that impact zones in East Kalimantan occupy 99% of the harvested sub-block and applied an adjustment factor of 1.01. Based on the more extensive mapping of infrastructure in this study, we recommend revising this adjustment factor to 1.45 (1–0.69), resulting in a 31% reduction in the area-based estimate of carbon emissions from selective logging in the region.
We did not detect any relationship between the observed distribution of logging impacts and slope, stream configuration, topographic position, or carbon stocks. Even though slopes were steeper in logged portions of cutting blocks B-2011 and D, the opposite was observed in G (Table 2). This confounds the intuitive assumption that logging crews avoid steep areas because of operational challenges, and confirms the finding of Griscom et al. (2014), who find no relationship between slope and skidding, felling, or hauling emissions. Map interpretation of stream corridors derived from lidar DTMs and flow accumulation models revealed no apparent spatial correlation between no-impact zones and streams or topography. If we accept carbon density as a proxy for standing timber volume, the insignificant difference in carbon density between no impact zones and unreported zones suggests that logging impact distribution cannot be explained by manager selection for areas with higher volumes of commercial timber. More research is needed to understand why concession operators choose to access certain areas and leave others un-logged.
We believe a large opportunity exists to reduce the environmental impacts of logging through informed location of no-impact zones. For example, pre-harvest planning could position impact zones away from sensitive riparian and steep slope areas (Figueiredo, 2007). Also, the spatial distribution of no-impact zones could be optimized to maximize landscape connectivity while minimizing edge effects (Putz and Ruslandi, 2015).
4.2 Comparison of lidar-based and GPS infrastructure maps
The impact zone areas mapped from the RDM lidar model agreed well with GPS impact zones (total predicted impact zone was only 1% larger than the GPS impact zone), but digitized lines did not reliably fall within the 10 m GPS buffer. (Fig. 8 and Fig. 9, Table 3). We believe time lag (time interval between harvest and lidar data acquisition) is the major contributor to this lack of agreement; the strong negative correlation between hauling impact zone accuracy and time lag indicates that regrowth can rapidly obscure lidar-derived disturbance signals in the RDM. This problem is pronounced in dipterocarp forests where regrowth is particularly rapid (Berry et al., 2010 and Banin et al., 2014). The utility of commercial-grade lidar as a tool for delineating logging infrastructure in tropical forests is therefore time-limited. Based on the above results, we recommend collecting lidar data <2 years after harvest completion. It may be feasible to extend the RDM signal by adjusting the RDM height corridor as a function of time since harvest (we used a fixed corridor of 1 to 10 m aboveground), but this would require careful tracking of time since harvest and vegetation regrowth rates in each sub-block. Adjusting corridor heights will also complicate the comparability of data among sub-blocks, so should be done only after additional analyses to calibrate RDM corridor values with field-based metrics and as a function of time lag.
4.3 Skid trail density and haul road extent
Despite apparently large differences in operating conditions (e.g. slopes), skid trail density did not exhibit a clear pattern or relationship to available explanatory variables. The inclusion of concession G suggests a false pattern between slope and skid trail density, because concession G happens to contain the gentlest slopes, and was harvested by sub-contractors with few incentives for skidding efficiency. We believe that the logging sub-contractor’s lack of investment in skid trail design influenced skid trail density more than slope.
Our ability to evaluate the variability of hauling extent with the lidar RDM was limited by inaccuracies associated with time lags discussed above, but we were nonetheless surprised by the low variability among cutting blocks with time lags <2 years (9% of mean at 95% confidence interval). Our field experience indicates that concession managers usually plan haul road networks at the outset, and their driving concern is access to timber, with little regard for environmental conditions (e.g. slope) or impacts (e.g. soil erosion). Therefore, the opportunities to reduce the impacts of haul roads are substantial and warrant attention.
Our assessment of hauling extent follows the RIL-C methodology (Verified Carbon Standard, 2015) to measure haul road widths as “the distance between the nearest tree trunks (>10 cm DBH) perpendicular to any given point along the haul road centerline.” This differs from previous studies (White et al., 2010, Azizi et al., 2014 and Pearson et al., 2014) that delineated haul road corridor widths between from adjacent tree canopies, thereby underestimating the hauling impacts because overhanging vegetation obscures the true disturbance corridor. The more accurate VCS methodology approach requires lidar or field measurements, since 2-dimensional images such as those used in the aforementioned studies cannot perceive understory vegetation structure (see Fig. S1).
4.4 Forest structure and carbon density
Carbon density maps derived from the Fourier-based biomass model accurately describe the natural variation in biomass and detect biomass differences between logged and unlogged areas. Despite the limited calibration sample (28 prism plots), the map-based distribution of biomass for the 6 cutting blocks closely approximates the biomass distribution observed in the field (Fig. 10), and biomass estimation error (17%) is only slightly higher than errors reported in other studies with more robust calibration datasets (10–14%: Mascaro et al., 2011, Meyer et al., 2013 and Réjou-Méchain et al., 2015). Dendritic patterns of reduced biomass shown on the map correspond to roads and skid trails mapped by both GPS and the lidar RDM (Fig. 11), illustrating the sensitivity of the Fourier model to logging-induced biomass losses.
Even more compelling, the Fourier-based biomass model provides an accurate estimate of the magnitude of logging disturbance. When correcting for time lag, the 18 ± 8% lidar-estimated biomass reduction in skidding/felling zones closely approximates the skidding and felling emissions calculated from field measurements by Griscom et al., 2014, who observed that 17 ± 1% of aboveground live biomass is destroyed by skidding and felling. This finding demonstrates the utility of carbon density maps derived from airborne lidar to detect and measure vegetation damage and associated carbon emissions from selectively logged tropical forests. More work is needed to test this lidar-based strategy across a range of tropical ecosystems to operationalize the approach.
4.5 Applications to measurement of RIL performance
The results of this study are relevant to two of the four impact parameters specified by the RIL-C VCS Performance Method Module for North and East Kalimantan (VCS, 2015): SKID and HAUL. SKID is calculated as the product of skidding damage (SKIDdam) per meter of skid trail, and skidding density (SKIDdens), defined as the average length of skid trails in harvest areas (m ha−1). The methodology definition for harvest area is identical to impact zone as defined in this paper, and should not be confused with the total annual cutting block area that includes no-impact zones. Skidding damage measures the performance of skidding operations, while skidding density quantifies the benefits of skid trail planning, as prescribed by many forestry standards (Dykstra and Heinrich, 1996, Forest Stewardship Council, 2014 and Tropical Forest Foundation Indonesia, 2015). Direct field-based measurements of damage are needed to estimate skidding damage, but skidding density can be calculated directly from lidar RDMs and/or GPS tracks once sufficient field data is collected to estimate the mean skid-to-stump buffer distances.
Extensive skidding maps, like those derived in this study, can be used to evaluate the effectiveness of skid trail planning. Our findings emphasize the importance of collecting sufficient GPS or lidar data (as indicated in the VCS Kalimantan module) in order to capture the full extent of no-impact zones and accurately evaluate skidding efficiency. In this study, the high skid trail density in Concession G presents an opportunity for impact reduction through incentivizing owner operation and operator training, capacity building, and supervision. Focusing incentives on measurable results in the form of decreased skid trail density can ensure adequate performance through reduced carbon emissions.
For non-carbon ecosystem services, skidding extent, as mapped by GPS or lidar, can be used to monitor and incentivize RIL improvements for water quality, sediment stabilization, and habitat protection, as well as reduce costs (Broza et al., 2012). Surprisingly, our results suggest minimal influence of environmental conditions on the location and extent of logging impacts. If RIL initiatives can identify priority conservation zones and inform strategic spatial planning of harvest operations, monitored skidding extent (via GPS or lidar) can verify compliance.
In the RIL-C methodology, the HAUL impact parameter is calculated as “the area, in square meters per ha of harvest area, of haul road and log landing corridors.” (VCS, 2015), which is equivalent to this study’s definition of hauling extent divided by 10,000. Our finding that time lag generates considerable error in haul road width confirms the judgement of the RIL-C methodology to derive baseline values from a combination of corridor widths measured in the field and haul road lengths measured from lidar and/or GPS. Our field-based width measurements combined with lidar measurements of haul road length yield a baseline hauling extent of 3.2% (323 m2 ha−1), 23% higher than the hauling extent reported from lidar here, but within the range of error for cutting blocks with time lags <2 years (3.5 ± 0.3%). We believe that future lidar measurements could be used to produce more consistent estimates of haul road width if collected <2 years after harvest completion, but we currently endorse the combined lidar and field-based estimates of the VCS methodology.
4.6 Recommendations for future initiatives to monitor logging infrastructure
Lidar is a powerful tool for mapping and monitoring forest condition and change; much of the information available from lidar is underutilized, particularly for applications to improve forest management. For example, lidar data of the type and quality used in this study can supplement national inventories (Henry et al., 2015), estimate standing timber volumes (Junttila et al., 2015), model water flows (Murphy et al., 2008), assess erosion impacts and siltation (Milodowski et al., 2015), extract forest structure metrics (this study, Palace et al., 2015), quantify carbon stocks (this study, Asner et al., 2012 and Gonçalves, 2014), estimate carbon density change (this study, McRoberts et al., 2015), and map logging infrastructure (this study, D’Oliveira et al., 2012). While agreement with GPS-based infrastructure maps is promising, more research is needed to quantify the uncertainties of lidar based infrastructure mapping, and we recommend that future studies repeat RDM-based infrastructure mapping with separate interpreters to assess variability (Olofsson et al., 2014).
The costs, timing, and logistical challenges of acquisition (in that order) remain substantial barriers to broad-scale use of lidar for infrastructure mapping. While many of the inaccuracies in lidar-based infrastructure mapping reported above can be reduced through minimizing time lags, there are practical limitations to achieving a properly constrained time-window for lidar collection. For this study, we aimed to restrict time lags to <2 years, but were forced to extend it beyond our limit due to logistical and cloud cover challenges that impeded lidar acquisition over our remote study area. For these reasons, we recommend use of ground-based GPS for mapping wall-to-wall infrastructure in similar remote tropical landscapes, because it can be easily integrated into the ground-based monitoring plan needed for other impact parameters (e.g. bucking and felling inefficiencies), is less susceptible to time lag errors (old skid trails are usually easy to detect on the ground), and remains substantially cheaper.
The average rate of infrastructure mapping using ground-based GPS reported in Griscom et al. (2014) was 4 km person−1 day−1, or 250 person-hours to map a 900 ha cutting block with an average of 121 km of infrastructure. Assuming a wage rate of $1.00 h−1 (a generous inflation of the World Bank wage of $0.65 h−1 for Indonesia in 2006, Oostendorp, 2013), GPS-based infrastructure mapping would cost ∼$0.25 ha−1, more than an order-of-magnitude less than the average cost of current-day commercial airborne lidar. Costs would be lower if GPS tracking were integrated into harvest operations. One potential strategy could promote or require that GPS units are attached to bulldozers to map infrastructure during skidding (McDonald et al., 2002).
5 Conclusion
With lidar, we demonstrate the importance of comprehensive infrastructure and carbon density mapping for accurate measurement and monitoring of selective logging impacts. Our analysis reveals that a large proportion of designated harvest blocks (37%) are not directly affected by timber extraction. The distribution of these no-impact zones does not correspond with riparian areas or steep slopes, indicating large opportunities to reduce the impact of logging without reducing harvest volumes. Lidar helped to improve existing performance monitoring in East Kalimantan, Indonesia, and we believe it will likely play a role in logging impact assessments elsewhere in the future. However, ground-based, GPS-supported mapping of logging infrastructure is still the more affordable and efficient near-term option for monitoring selective logging operations.
Acknowledgments
We thank Indonesian logging concession owners, managers, and staff for their cooperation and hospitality. For field work, we are indebted to the Nature Conservancy Indonesia field team and program staff. For logistical support, we thank Delon Martinus, Jeffrey Benz, Lee Hon Chuan and Credent Teknologi. For input into the design and execution of the analysis, we thank Jeffrey Evans, Michael Keller, Rosa Goodman, and Jonathan Leiner. For review of this manuscript, we thank Joe Fargione, Marissa Ahlering, Peter Kareiva, and Jack Putz. Funding for this research was provided by the Australian Department of Agriculture, Fisheries, and Forestry (DAFF) and the Norwegian Agency for Development Cooperation (NORAD). We also thank the National Center for Ecological Analysis and Synthesis for hosting the datasets developed here, as a contribution to the SNAP: Science for Nature and People Forest Sharing or Sparing Working Group.
G.P. Asner, J.K. Clark, J. Mascaro, G.a. Galindo García, K.D. Chadwick, D.a. Navarrete Encinales, G. Paez-Acosta, E. Cabrera Montenegro, T. Kennedy-Bowdoin, Á. Duque, a. Balaji, P. von Hildebrand, L. Maatoug, J.F. Phillips Bernal, a.P. Yepes Quintero, D.E. Knapp, M.C. García Dávila, J. Jacobson, M.F. Ordóñez
High-resolution mapping of forest carbon stocks in the Colombian Amazon
L. Banin, S.L. Lewis, G. Lopez-Gonzalez, T.R. Baker, C.a. Quesada, K.J. Chao, D.F.R.P. Burslem, R. Nilus, K. Abu Salim, H.C. Keeling, S. Tan, S.J. Davies, A. Monteagudo Mendoza, R. Vásquez, J. Lloyd, D.a. Neill, N. Pitman, O.L. Phillips
Tropical forest wood production: a cross-continental comparison
Blaser, J., Sarre, A., Poore, D., Johnson, S., 2011. Status of tropical forest management 2011. ITTO Technical Series, No 38. International Tropical Timber Organization, Yokohama, Japan.
Broza, K.T., Garrastazu, M.C., Braz, E.M., de Mattos, P.P., Rosot, M.A.D., Maldonado, F.D., Sack, L., Tentor, F.R., Fortini, C., 2012. Etapas do planejamento do projeto modeflora em SIG livre. Semin. ATUALIZAÇÃO EM SENSORIAMENTO REMOTO E Sist. INFORMAÇÕES GEOGRÁFICAS Apl. À Eng. Florest. 7.
Estimating forest biomass and identifying low-intensity logging areas using airborne scanning lidar in Antimary State Forest, Acre State, Western Brazilian Amazon
Forest Stewardship Council, 2014. FSC Harmonised Certification Bodies’ Forest Stewardship Standard for the Republic of Indonesia: The Forest Stewardship Standard for the Republic of Indonesia (No. FSC-STD-IDN-01-01-2013 Indonesia).
M. Henry, M. Réjou-Méchain, M.C. Jara, C. Wayson, D. Piotto, J. Westfall, J.M.M. Fuentes, F.A. Guier, H.C. Lombis, E.C. López, R.C. Lara, K.C. Rojas, J. Del Águila Pasquel, Á.D. Montoya, J.F. Vega, A.J. Galo, O.R. López, L.G. Marklund, F. Milla, J. de Jesús Návar Cahidez, E.O. Malavassi, J. Pérez, C.R. Zea, L.R. García, R.R. Pons, C. Sanquetta, C. Scott, M. Zapata-Cuartas, L. Saint-André
An overview of existing and promising technologies for national forest monitoring
F.E. Putz, P.A. Zuidema, T. Synnott, M. Peña-Claros, M.A. Pinard, D. Sheil, Jerome K. Vanclay, P. Sist, S. Gourlet-Fleury, B. Griscom, J. Palmer, R. Zagt
Sustaining conservation values in selectively logged tropical forests: the attained and the attainable
No comments:
Post a Comment